[R] calculating martingale residual on new data using "predict.coxph"
Shi, Tao
shidaxia at yahoo.com
Tue Nov 23 02:36:54 CET 2010
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
>
>
More information about the R-help
mailing list