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

[David Winsemius]

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