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) > 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.
