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
