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