[R] calculating martingale residual on new data using "predict.coxph"

[Shi, Tao]

Thank you, David!  You've been always helpful!

...Tao




----- Original Message ----
> From: David Winsemius <dwinsemius at comcast.net>
> To: "Shi, Tao" <shidaxia at yahoo.com>
> Cc: r-help at r-project.org; dieter.menne at menne-biomed.de; r_tingley at 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 <dwinsemius at comcast.net>
> >>  To: David Winsemius <dwinsemius at comcast.net>
> >>  Cc: "Shi, Tao" <shidaxia at yahoo.com>; r-help at r-project.org;
> >> dieter.menne at menne-biomed.de; r_tingley at 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
> 
>