Re: [R] calculating martingale residual on new data using predict.coxph
This feature has been added in survival 2.36-1, which is now on CRAN. (2.36-2 should appear in another day or so) Terry T. -begin included message I was trying to use predict.coxph to calculate martingale residuals on a test data, however, as pointed out before http://tolstoy.newcastle.edu.au/R/e4/help/08/06/13508.html predict(mycox1, newdata, type=expected) is not implemented yet. __ R-help@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.
Re: [R] calculating martingale residual on new data using predict.coxph
Thank you, Terry! - Original Message From: Terry Therneau thern...@mayo.edu To: Shi, Tao shida...@yahoo.com Cc: r-help@r-project.org; dieter.me...@menne-biomed.de; r_ting...@hotmail.com Sent: Mon, November 22, 2010 6:11:15 AM Subject: Re: calculating martingale residual on new data using predict.coxph This feature has been added in survival 2.36-1, which is now on CRAN. (2.36-2 should appear in another day or so) Terry T. -begin included message I was trying to use predict.coxph to calculate martingale residuals on a test data, however, as pointed out before http://tolstoy.newcastle.edu.au/R/e4/help/08/06/13508.html predict(mycox1, newdata, type=expected) is not implemented yet. __ R-help@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.
Re: [R] calculating martingale residual on new data using
Thank you for the advice, Frank! ...Tao - Original Message From: Frank Harrell f.harr...@vanderbilt.edu To: r-help@r-project.org Sent: Sun, November 21, 2010 5:49:36 AM Subject: Re: [R] calculating martingale residual on new data using The tendency is to use residual-like diagnostics on the entire dataset that was available for model development. For test data we typically run predictive accuracy analyses. For example, one of the strongest validations is to show, in a high-resolution calibration plot, that absolute predictions (e.g., probability of survival at 2 years) are accurate. Frank - Frank Harrell Department of Biostatistics, Vanderbilt University -- View this message in context: http://r.789695.n4.nabble.com/calculating-martingale-residual-on-new-data-using-predict-coxph-tp3050712p3052377.html Sent from the R help mailing list archive at Nabble.com. __ R-help@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. __ R-help@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.
Re: [R] calculating martingale residual on new data using predict.coxph
Thank you, David! You've been always helpful! ...Tao - Original Message From: David Winsemius dwinsem...@comcast.net To: Shi, Tao shida...@yahoo.com Cc: r-help@r-project.org; dieter.me...@menne-biomed.de; r_ting...@hotmail.com Sent: Sun, November 21, 2010 5:50:31 AM Subject: Re: [R] calculating martingale residual on new data using predict.coxph On Nov 21, 2010, at 3:42 AM, Shi, Tao wrote: Hi David, Thanks, but I don't quite follow your examples below. I wasn't really sure they did anything useful anyway. The residuals you calculated are still based on the training data from which your cox model was generated. I'm interested in the testing data. The survest function in rms and the survfit function in survival will calculate survival probabilities given a model and newdata, and depending on your definition of residual you could take the difference between the calculation and validation data. That must be what happens (at least at a gross level of description) when Harrell runs his validate function on his cph models in the rms/Design package, although I don't know if something that you would recognize as a martingale residual is an identifiable intermediate. If you are using survfit, it would appear from my reading that you would need to set the individual parameter to TRUE. I'm assuming you planned to calculate these (1- expected) at the event times of the validation cohort (which it appears the default method fixes via the censor argument)? --David Best, ...Tao - Original Message From: David Winsemius dwinsem...@comcast.net To: David Winsemius dwinsem...@comcast.net Cc: Shi, Tao shida...@yahoo.com; r-help@r-project.org; dieter.me...@menne-biomed.de; r_ting...@hotmail.com Sent: Fri, November 19, 2010 10:53:26 AM Subject: Re: [R] calculating martingale residual on new data using predict.coxph On Nov 19, 2010, at 12:50 PM, David Winsemius wrote: On Nov 19, 2010, at 12:32 PM, Shi, Tao wrote: Hi list, I was trying to use predict.coxph to calculate martingale residuals on a test data, however, as pointed out before What about resid(fit) ? It's my reading of Therneau Gramsch [and of help(coxph.object) ] that they consider those martingale residuals. The manner in which I _thought_ this would work was to insert some dummy cases into the original data and then to get residuals by weighting the cases appropriately. That doesn't seem to be as successful as I imagined: test1 - list(time=c(4,3,1,1,2,2,3,3), weights=c(rep(1,7), 0), + status=c(1,1,1,0,1,1,0,1), + x=c(0,2,1,1,1,0,0,1), + sex=c(0,0,0,0,1,1,1,1)) coxph(Surv(time, status) ~ x , test1, weights=weights)$weights Error in fitter(X, Y, strats, offset, init, control, weights = weights, : Invalid weights, must be 0 # OK then make it a small number test1 - list(time=c(4,3,1,1,2,2,3,3), weights=c(rep(1,7), 0.01), + status=c(1,1,1,0,1,1,0,1), + x=c(0,2,1,1,1,0,0,1), + sex=c(0,0,0,0,1,1,1,1)) print(resid( coxph(Surv(time, status) ~ x , test1,weights=weights) ) ,digits=3) 1 2 3 4 5 6 7 8 -0.6410 -0.5889 0.8456 -0.1544 0.4862 0.6931 -0.6410 0.0509 Now take out constructed case and weights test1 - list(time=c(4,3,1,1,2,2,3), + status=c(1,1,1,0,1,1,0), + x=c(0,2,1,1,1,0,0), + sex=c(0,0,0,0,1,1,1)) print(resid( coxph(Surv(time, status) ~ x , test1) ) ,digits=3) 1 2 3 4 5 6 7 -0.632 -0.589 0.846 -0.154 0.486 0.676 -0.632 Expecting approximately the same residuals for first 7 cases but not really getting it. There must be something about weights in coxph that I don't understand, unless a one-hundreth of a case gets up indexed inside the machinery of coxph? Still think that inserting a single constructed case into a real dataset of sufficient size ought to be able to yield some sort of estimate, and only be a minor perturbation, although I must admit I'm having trouble figuring out ... why are we attempting such a maneuver? The notion of residuals around constructed cases makes me statistically suspicious, although I suppose that is just some sort of cumulative excess/deficit death fraction. http://tolstoy.newcastle.edu.au/R/e4/help/08/06/13508.html predict(mycox1, newdata, type=expected) is not implemented yet. Dieter suggested to use 'cph' and 'predict.Design', but from my reading so far, I'm not sure they can do that. Do you other ways to calculate martingale residuals on a new data? Thank you very much! ...Tao --David
Re: [R] calculating martingale residual on new data using predict.coxph
Hi David, Thanks, but I don't quite follow your examples below. The residuals you calculated are still based on the training data from which your cox model was generated. I'm interested in the testing data. Best, ...Tao - Original Message From: David Winsemius dwinsem...@comcast.net To: David Winsemius dwinsem...@comcast.net Cc: Shi, Tao shida...@yahoo.com; r-help@r-project.org; dieter.me...@menne-biomed.de; r_ting...@hotmail.com Sent: Fri, November 19, 2010 10:53:26 AM Subject: Re: [R] calculating martingale residual on new data using predict.coxph On Nov 19, 2010, at 12:50 PM, David Winsemius wrote: On Nov 19, 2010, at 12:32 PM, Shi, Tao wrote: Hi list, I was trying to use predict.coxph to calculate martingale residuals on a test data, however, as pointed out before What about resid(fit) ? It's my reading of Therneau Gramsch [and of help(coxph.object) ] that they consider those martingale residuals. The manner in which I _thought_ this would work was to insert some dummy cases into the original data and then to get residuals by weighting the cases appropriately. That doesn't seem to be as successful as I imagined: test1 - list(time=c(4,3,1,1,2,2,3,3), weights=c(rep(1,7), 0), +status=c(1,1,1,0,1,1,0,1), +x=c(0,2,1,1,1,0,0,1), +sex=c(0,0,0,0,1,1,1,1)) coxph(Surv(time, status) ~ x , test1, weights=weights)$weights Error in fitter(X, Y, strats, offset, init, control, weights = weights, : Invalid weights, must be 0 # OK then make it a small number test1 - list(time=c(4,3,1,1,2,2,3,3), weights=c(rep(1,7), 0.01), +status=c(1,1,1,0,1,1,0,1), +x=c(0,2,1,1,1,0,0,1), +sex=c(0,0,0,0,1,1,1,1)) print(resid( coxph(Surv(time, status) ~ x , test1,weights=weights) ) ,digits=3) 12 3 45 6 78 -0.6410 -0.5889 0.8456 -0.1544 0.4862 0.6931 -0.6410 0.0509 Now take out constructed case and weights test1 - list(time=c(4,3,1,1,2,2,3), +status=c(1,1,1,0,1,1,0), +x=c(0,2,1,1,1,0,0), +sex=c(0,0,0,0,1,1,1)) print(resid( coxph(Surv(time, status) ~ x , test1) ) ,digits=3) 1 2 3 4 5 6 7 -0.632 -0.589 0.846 -0.154 0.486 0.676 -0.632 Expecting approximately the same residuals for first 7 cases but not really getting it. There must be something about weights in coxph that I don't understand, unless a one-hundreth of a case gets up indexed inside the machinery of coxph? Still think that inserting a single constructed case into a real dataset of sufficient size ought to be able to yield some sort of estimate, and only be a minor perturbation, although I must admit I'm having trouble figuring out ... why are we attempting such a maneuver? The notion of residuals around constructed cases makes me statistically suspicious, although I suppose that is just some sort of cumulative excess/deficit death fraction. http://tolstoy.newcastle.edu.au/R/e4/help/08/06/13508.html predict(mycox1, newdata, type=expected) is not implemented yet. Dieter suggested to use 'cph' and 'predict.Design', but from my reading so far, I'm not sure they can do that. Do you other ways to calculate martingale residuals on a new data? Thank you very much! ...Tao --David Winsemius, MD West Hartford, CT __ R-help@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.
Re: [R] calculating martingale residual on new data using
The tendency is to use residual-like diagnostics on the entire dataset that was available for model development. For test data we typically run predictive accuracy analyses. For example, one of the strongest validations is to show, in a high-resolution calibration plot, that absolute predictions (e.g., probability of survival at 2 years) are accurate. Frank - Frank Harrell Department of Biostatistics, Vanderbilt University -- View this message in context: http://r.789695.n4.nabble.com/calculating-martingale-residual-on-new-data-using-predict-coxph-tp3050712p3052377.html Sent from the R help mailing list archive at Nabble.com. __ R-help@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.
Re: [R] calculating martingale residual on new data using predict.coxph
On Nov 21, 2010, at 3:42 AM, Shi, Tao wrote: Hi David, Thanks, but I don't quite follow your examples below. I wasn't really sure they did anything useful anyway. The residuals you calculated are still based on the training data from which your cox model was generated. I'm interested in the testing data. The survest function in rms and the survfit function in survival will calculate survival probabilities given a model and newdata, and depending on your definition of residual you could take the difference between the calculation and validation data. That must be what happens (at least at a gross level of description) when Harrell runs his validate function on his cph models in the rms/Design package, although I don't know if something that you would recognize as a martingale residual is an identifiable intermediate. If you are using survfit, it would appear from my reading that you would need to set the individual parameter to TRUE. I'm assuming you planned to calculate these (1- expected) at the event times of the validation cohort (which it appears the default method fixes via the censor argument)? -- David Best, ...Tao - Original Message From: David Winsemius dwinsem...@comcast.net To: David Winsemius dwinsem...@comcast.net Cc: Shi, Tao shida...@yahoo.com; r-help@r-project.org; dieter.me...@menne-biomed.de; r_ting...@hotmail.com Sent: Fri, November 19, 2010 10:53:26 AM Subject: Re: [R] calculating martingale residual on new data using predict.coxph On Nov 19, 2010, at 12:50 PM, David Winsemius wrote: On Nov 19, 2010, at 12:32 PM, Shi, Tao wrote: Hi list, I was trying to use predict.coxph to calculate martingale residuals on a test data, however, as pointed out before What about resid(fit) ? It's my reading of Therneau Gramsch [and of help(coxph.object) ] that they consider those martingale residuals. The manner in which I _thought_ this would work was to insert some dummy cases into the original data and then to get residuals by weighting the cases appropriately. That doesn't seem to be as successful as I imagined: test1 - list(time=c(4,3,1,1,2,2,3,3), weights=c(rep(1,7), 0), +status=c(1,1,1,0,1,1,0,1), +x=c(0,2,1,1,1,0,0,1), +sex=c(0,0,0,0,1,1,1,1)) coxph(Surv(time, status) ~ x , test1, weights=weights)$weights Error in fitter(X, Y, strats, offset, init, control, weights = weights, : Invalid weights, must be 0 # OK then make it a small number test1 - list(time=c(4,3,1,1,2,2,3,3), weights=c(rep(1,7), 0.01), +status=c(1,1,1,0,1,1,0,1), +x=c(0,2,1,1,1,0,0,1), +sex=c(0,0,0,0,1,1,1,1)) print(resid( coxph(Surv(time, status) ~ x , test1,weights=weights) ) ,digits=3) 12 3 45 6 78 -0.6410 -0.5889 0.8456 -0.1544 0.4862 0.6931 -0.6410 0.0509 Now take out constructed case and weights test1 - list(time=c(4,3,1,1,2,2,3), +status=c(1,1,1,0,1,1,0), +x=c(0,2,1,1,1,0,0), +sex=c(0,0,0,0,1,1,1)) print(resid( coxph(Surv(time, status) ~ x , test1) ) ,digits=3) 1 2 3 4 5 6 7 -0.632 -0.589 0.846 -0.154 0.486 0.676 -0.632 Expecting approximately the same residuals for first 7 cases but not really getting it. There must be something about weights in coxph that I don't understand, unless a one-hundreth of a case gets up indexed inside the machinery of coxph? Still think that inserting a single constructed case into a real dataset of sufficient size ought to be able to yield some sort of estimate, and only be a minor perturbation, although I must admit I'm having trouble figuring out ... why are we attempting such a maneuver? The notion of residuals around constructed cases makes me statistically suspicious, although I suppose that is just some sort of cumulative excess/deficit death fraction. http://tolstoy.newcastle.edu.au/R/e4/help/08/06/13508.html predict(mycox1, newdata, type=expected) is not implemented yet. Dieter suggested to use 'cph' and 'predict.Design', but from my reading so far, I'm not sure they can do that. Do you other ways to calculate martingale residuals on a new data? Thank you very much! ...Tao --David Winsemius, MD West Hartford, CT David Winsemius, MD West Hartford, CT __ R-help@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.
Re: [R] calculating martingale residual on new data using predict.coxph
On Nov 19, 2010, at 12:32 PM, Shi, Tao wrote: Hi list, I was trying to use predict.coxph to calculate martingale residuals on a test data, however, as pointed out before What about resid(fit) ? It's my reading of Therneau Gramsch [and of help(coxph.object) ] that they consider those martingale residuals. http://tolstoy.newcastle.edu.au/R/e4/help/08/06/13508.html predict(mycox1, newdata, type=expected) is not implemented yet. Dieter suggested to use 'cph' and 'predict.Design', but from my reading so far, I'm not sure they can do that. Do you other ways to calculate martingale residuals on a new data? Thank you very much! ...Tao David Winsemius, MD West Hartford, CT __ R-help@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.
Re: [R] calculating martingale residual on new data using predict.coxph
Hi David, Thank you for the quick reply! resid(fit) only gives the residuals on the training data not on test data. ...Tao - Original Message From: David Winsemius dwinsem...@comcast.net To: Shi, Tao shida...@yahoo.com Cc: r-help@r-project.org; dieter.me...@menne-biomed.de; r_ting...@hotmail.com Sent: Fri, November 19, 2010 9:50:06 AM Subject: Re: [R] calculating martingale residual on new data using predict.coxph On Nov 19, 2010, at 12:32 PM, Shi, Tao wrote: Hi list, I was trying to use predict.coxph to calculate martingale residuals on a test data, however, as pointed out before What about resid(fit) ? It's my reading of Therneau Gramsch [and of help(coxph.object) ] that they consider those martingale residuals. http://tolstoy.newcastle.edu.au/R/e4/help/08/06/13508.html predict(mycox1, newdata, type=expected) is not implemented yet. Dieter suggested to use 'cph' and 'predict.Design', but from my reading so far, I'm not sure they can do that. Do you other ways to calculate martingale residuals on a new data? Thank you very much! ...Tao David Winsemius, MD West Hartford, CT __ R-help@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.
Re: [R] calculating martingale residual on new data using predict.coxph
On Nov 19, 2010, at 12:50 PM, David Winsemius wrote: On Nov 19, 2010, at 12:32 PM, Shi, Tao wrote: Hi list, I was trying to use predict.coxph to calculate martingale residuals on a test data, however, as pointed out before What about resid(fit) ? It's my reading of Therneau Gramsch [and of help(coxph.object) ] that they consider those martingale residuals. The manner in which I _thought_ this would work was to insert some dummy cases into the original data and then to get residuals by weighting the cases appropriately. That doesn't seem to be as successful as I imagined: test1 - list(time=c(4,3,1,1,2,2,3,3), weights=c(rep(1,7), 0), + status=c(1,1,1,0,1,1,0,1), + x=c(0,2,1,1,1,0,0,1), + sex=c(0,0,0,0,1,1,1,1)) coxph(Surv(time, status) ~ x , test1, weights=weights)$weights Error in fitter(X, Y, strats, offset, init, control, weights = weights, : Invalid weights, must be 0 # OK then make it a small number test1 - list(time=c(4,3,1,1,2,2,3,3), weights=c(rep(1,7), 0.01), + status=c(1,1,1,0,1,1,0,1), + x=c(0,2,1,1,1,0,0,1), + sex=c(0,0,0,0,1,1,1,1)) print(resid( coxph(Surv(time, status) ~ x , test1,weights=weights) ) ,digits=3) 1 2 3 4 5 6 7 8 -0.6410 -0.5889 0.8456 -0.1544 0.4862 0.6931 -0.6410 0.0509 Now take out constructed case and weights test1 - list(time=c(4,3,1,1,2,2,3), + status=c(1,1,1,0,1,1,0), + x=c(0,2,1,1,1,0,0), + sex=c(0,0,0,0,1,1,1)) print(resid( coxph(Surv(time, status) ~ x , test1) ) ,digits=3) 1 2 3 4 5 6 7 -0.632 -0.589 0.846 -0.154 0.486 0.676 -0.632 Expecting approximately the same residuals for first 7 cases but not really getting it. There must be something about weights in coxph that I don't understand, unless a one-hundreth of a case gets up indexed inside the machinery of coxph? Still think that inserting a single constructed case into a real dataset of sufficient size ought to be able to yield some sort of estimate, and only be a minor perturbation, although I must admit I'm having trouble figuring out ... why are we attempting such a maneuver? The notion of residuals around constructed cases makes me statistically suspicious, although I suppose that is just some sort of cumulative excess/deficit death fraction. http://tolstoy.newcastle.edu.au/R/e4/help/08/06/13508.html predict(mycox1, newdata, type=expected) is not implemented yet. Dieter suggested to use 'cph' and 'predict.Design', but from my reading so far, I'm not sure they can do that. Do you other ways to calculate martingale residuals on a new data? Thank you very much! ...Tao -- David Winsemius, MD West Hartford, CT __ R-help@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.
