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

[Shi, Tao]

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