Dear R-help, I am using R 2.14.1 on Windows 7. I would like to produce a plot like the attached - although simplified to actual vs. Predicted survival time with distinguishing marks for censored and observed points. I have a dataset and have fitted a Cox model to it. In an attempt to visualise how accurate the model is it would be ideal if I could plot the actual survival times against the predicted survival times. I have been looking on the internet to see if there are ways to do this in R. The only post I found (https://stat.ethz.ch/pipermail/r-help/2009-February/189888.html) that seemed directly relevant suggested that I shouldn't be generating survival times at all. Given that, I was concerned about proceeding but I would like to have access to a plot to make a decision on its usefulness. I appreciate that there are predict.coxph and predict.cph options available to me. My first attempt was as follows: # fit Cox model # coxfita = coxph(Surv(tsecond,seccens)~stroke(smess)+rels(smess)+asleep(smess)+eeg1(smess)+eeg2(smess)+ct1(smess)+ct2(smess)+treat(smess),data=smess) # Find censored and observed groups # messcens <- subset(smess,seccens==1) messobs <- subset(smess,seccens==0) # Obtain predicted survival times # explp <- exp(predict(coxfita,type="lp")) explp2 <- mean(ssmess$tsecond,na.rm=TRUE)*explp smess2 <- data.frame(ssmess,explp2) # Find censored and observed groups # smesscens <- subset(smess2,seccens==1) smessobs <- subset(smess2,seccens==0) # Produce plot # plot(smesscens$explp2,messcens$tsecond,pch=4,col="blue",ylab="Actual Survival Time",xlab="Predicted Survival Time",main="Survival Times",xlim=c(0,3500),ylim=c(0,3500)) points(smessobs$explp2,messobs$tsecond,pch=4,col="red") This leads to the attached plot. It doesn't seem correct though as the predicted times all start over 500 days. Any suggestions would be very welcome. Many thanks, Laura -------------- next part -------------- A non-text attachment was scrubbed... Name: Actual vs. Survival LJB.pdf Type: application/pdf Size: 12927 bytes Desc: Actual vs. Survival LJB.pdf URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20120207/0c6e1e57/attachment.pdf>
Actual vs. predicted survival times
4 messages · David Winsemius, Bonnett, Laura, alanm (Alan Mitchell)
On Feb 7, 2012, at 10:32 AM, Bonnett, Laura wrote:
Dear R-help, I am using R 2.14.1 on Windows 7. I would like to produce a plot like the attached - although simplified to actual vs. Predicted survival time with distinguishing marks for censored and observed points. I have a dataset and have fitted a Cox model to it. In an attempt to visualise how accurate the model is it would be ideal if I could plot the actual survival times against the predicted survival times. I have been looking on the internet to see if there are ways to do this in R. The only post I found (https://stat.ethz.ch/pipermail/r-help/2009-February/189888.html ) that seemed directly relevant suggested that I shouldn't be generating survival times at all. Given that, I was concerned about proceeding but I would like to have access to a plot to make a decision on its usefulness. I appreciate that there are predict.coxph and predict.cph options available to me. My first attempt was as follows: # fit Cox model # coxfita = coxph(Surv(tsecond,seccens)~stroke(smess)+rels(smess) +asleep(smess)+eeg1(smess)+eeg2(smess)+ct1(smess)+ct2(smess) +treat(smess),data=smess) # Find censored and observed groups # messcens <- subset(smess,seccens==1) messobs <- subset(smess,seccens==0) # Obtain predicted survival times # explp <- exp(predict(coxfita,type="lp"))
That gives you relative risks, not survival times.
explp2 <- mean(ssmess$tsecond,na.rm=TRUE)*explp
Why are you multiplying times by relative risks? That makes no sense.
David. > smess2 <- data.frame(ssmess,explp2) > > # Find censored and observed groups # > smesscens <- subset(smess2,seccens==1) > smessobs <- subset(smess2,seccens==0) > > # Produce plot # > plot(smesscens$explp2,messcens$tsecond,pch=4,col="blue",ylab="Actual > Survival Time",xlab="Predicted Survival Time",main="Survival > Times",xlim=c(0,3500),ylim=c(0,3500)) > points(smessobs$explp2,messobs$tsecond,pch=4,col="red") > > This leads to the attached plot. It doesn't seem correct though as > the predicted times all start over 500 days. > > Any suggestions would be very welcome. > > Many thanks, > Laura > > <Actual vs. Survival > LJB.pdf>______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD West Hartford, CT
Dear David, Thank you for your input. I was under the impression that it was not possible to obtain survival times using the predict.coxph function. In a post to the R help group I noticed someone had suggested calculating survival times in the manner I have used although obviously this is not appropriate. Are you able to provide any suggestions as to how I might achieve the graph of actual vs. Predicted survival times? Many thanks, Laura -----Original Message----- From: David Winsemius [mailto:dwinsemius at comcast.net] Sent: 07 February 2012 16:59 To: Bonnett, Laura Cc: 'r-help at r-project.org' Subject: Re: [R] Actual vs. predicted survival times
On Feb 7, 2012, at 10:32 AM, Bonnett, Laura wrote:
Dear R-help, I am using R 2.14.1 on Windows 7. I would like to produce a plot like the attached - although simplified to actual vs. Predicted survival time with distinguishing marks for censored and observed points. I have a dataset and have fitted a Cox model to it. In an attempt to visualise how accurate the model is it would be ideal if I could plot the actual survival times against the predicted survival times. I have been looking on the internet to see if there are ways to do this in R. The only post I found (https://stat.ethz.ch/pipermail/r-help/2009-February/189888.html ) that seemed directly relevant suggested that I shouldn't be generating survival times at all. Given that, I was concerned about proceeding but I would like to have access to a plot to make a decision on its usefulness. I appreciate that there are predict.coxph and predict.cph options available to me. My first attempt was as follows: # fit Cox model # coxfita = coxph(Surv(tsecond,seccens)~stroke(smess)+rels(smess) +asleep(smess)+eeg1(smess)+eeg2(smess)+ct1(smess)+ct2(smess) +treat(smess),data=smess) # Find censored and observed groups # messcens <- subset(smess,seccens==1) messobs <- subset(smess,seccens==0) # Obtain predicted survival times # explp <- exp(predict(coxfita,type="lp"))
That gives you relative risks, not survival times.
explp2 <- mean(ssmess$tsecond,na.rm=TRUE)*explp
Why are you multiplying times by relative risks? That makes no sense.
David. > smess2 <- data.frame(ssmess,explp2) > > # Find censored and observed groups # > smesscens <- subset(smess2,seccens==1) > smessobs <- subset(smess2,seccens==0) > > # Produce plot # > plot(smesscens$explp2,messcens$tsecond,pch=4,col="blue",ylab="Actual > Survival Time",xlab="Predicted Survival Time",main="Survival > Times",xlim=c(0,3500),ylim=c(0,3500)) > points(smessobs$explp2,messobs$tsecond,pch=4,col="red") > > This leads to the attached plot. It doesn't seem correct though as > the predicted times all start over 500 days. > > Any suggestions would be very welcome. > > Many thanks, > Laura > > <Actual vs. Survival > LJB.pdf>______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD West Hartford, CT
Have you tried plotting obvserved survival against X*Beta? I believe the usual predictions from a cox model are just monotonic transformations of this. -Alan -----Original Message----- From: Bonnett, Laura [mailto:L.J.Bonnett at liverpool.ac.uk] Sent: Wed 2/8/2012 1:52 AM To: 'David Winsemius' Cc: 'r-help at r-project.org' Subject: Re: [R] Actual vs. predicted survival times Dear David, Thank you for your input. I was under the impression that it was not possible to obtain survival times using the predict.coxph function. In a post to the R help group I noticed someone had suggested calculating survival times in the manner I have used although obviously this is not appropriate. Are you able to provide any suggestions as to how I might achieve the graph of actual vs. Predicted survival times? Many thanks, Laura -----Original Message----- From: David Winsemius [mailto:dwinsemius at comcast.net] Sent: 07 February 2012 16:59 To: Bonnett, Laura Cc: 'r-help at r-project.org' Subject: Re: [R] Actual vs. predicted survival times
On Feb 7, 2012, at 10:32 AM, Bonnett, Laura wrote:
Dear R-help, I am using R 2.14.1 on Windows 7. I would like to produce a plot like the attached - although simplified to actual vs. Predicted survival time with distinguishing marks for censored and observed points. I have a dataset and have fitted a Cox model to it. In an attempt to visualise how accurate the model is it would be ideal if I could plot the actual survival times against the predicted survival times. I have been looking on the internet to see if there are ways to do this in R. The only post I found (https://stat.ethz.ch/pipermail/r-help/2009-February/189888.html ) that seemed directly relevant suggested that I shouldn't be generating survival times at all. Given that, I was concerned about proceeding but I would like to have access to a plot to make a decision on its usefulness. I appreciate that there are predict.coxph and predict.cph options available to me. My first attempt was as follows: # fit Cox model # coxfita = coxph(Surv(tsecond,seccens)~stroke(smess)+rels(smess) +asleep(smess)+eeg1(smess)+eeg2(smess)+ct1(smess)+ct2(smess) +treat(smess),data=smess) # Find censored and observed groups # messcens <- subset(smess,seccens==1) messobs <- subset(smess,seccens==0) # Obtain predicted survival times # explp <- exp(predict(coxfita,type="lp"))
That gives you relative risks, not survival times.
explp2 <- mean(ssmess$tsecond,na.rm=TRUE)*explp
Why are you multiplying times by relative risks? That makes no sense.
David. > smess2 <- data.frame(ssmess,explp2) > > # Find censored and observed groups # > smesscens <- subset(smess2,seccens==1) > smessobs <- subset(smess2,seccens==0) > > # Produce plot # > plot(smesscens$explp2,messcens$tsecond,pch=4,col="blue",ylab="Actual > Survival Time",xlab="Predicted Survival Time",main="Survival > Times",xlim=c(0,3500),ylim=c(0,3500)) > points(smessobs$explp2,messobs$tsecond,pch=4,col="red") > > This leads to the attached plot. It doesn't seem correct though as > the predicted times all start over 500 days. > > Any suggestions would be very welcome. > > Many thanks, > Laura > > <Actual vs. Survival > LJB.pdf>______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD West Hartford, CT