R help - Nov 2010 - Is there an equivalent to predict(..., type="linear") of a Proportional hazard model for a Cox model instead?

Hi all,

Is there an equivalent to predict(...,type="linear") of a Proportional
hazard
model for a Cox model instead?

For example, the Figure 13.12 in MASS (p384) is produced by:

(aids.ps <- survreg(Surv(survtime + 0.9, status) ~ state + T.categ + 
pspline(age, df=6), data = Aidsp))

zz <- predict(aids.ps, data.frame(state = factor(rep("NSW", 83),
levels =
levels(Aidsp$state)),
    T.categ = factor(rep("hs", 83), levels = levels(Aidsp$T.categ)),
age =
0:82), se = T, type = "linear")
plot(0:82, exp(zz$fit)/365.25, type = "l", ylim = c(0, 2), xlab =
"age", ylab =
"expected lifetime (years)")
lines(0:82, exp(zz$fit+1.96*zz$se.fit)/365.25, lty = 3, col = 2)
lines(0:82, exp(zz$fit-1.96*zz$se.fit)/365.25, lty = 3, col = 2)
rug(Aidsp$age + runif(length(Aidsp$age), -0.5, 0.5), ticksize = 0.015)

Is it possible to achieve something similar with a Cox model instead?

Is there a more detailed explanation of the "type" option for
predict.coxph than
what's in the help of predict.coxph? e.g. type=c("lp",
"risk", "expected",
"terms")

thanks in advance,

Ben

I manage to achieve similar results with a Cox model as follows but I don't 
really understand why we have to take the inverse of the linear prediction with 
the Cox model and why we do not need to divide by the number of days in the year
anymore?

Am I getting a similar result out of pure luck?

thanks in advance,

Ben

# MASS example with the proportional hazard model
par(mfrow = c(1, 2));
(aids.ps <- survreg(Surv(survtime + 0.9, status) ~ state + T.categ + 
pspline(age, df=6), data = Aidsp))

zz <- predict(aids.ps, data.frame(state = factor(rep("NSW", 83),
levels =
levels(Aidsp$state)),
    T.categ = factor(rep("hs", 83), levels = levels(Aidsp$T.categ)),
age =
0:82), se = T, type = "linear")
plot(0:82, exp(zz$fit)/365.25, type = "l", ylim = c(0, 2), xlab =
"age", ylab =
"expected lifetime (years)")
lines(0:82, exp(zz$fit+1.96*zz$se.fit)/365.25, lty = 3, col = 2)
lines(0:82, exp(zz$fit-1.96*zz$se.fit)/365.25, lty = 3, col = 2)
rug(Aidsp$age + runif(length(Aidsp$age), -0.5, 0.5), ticksize = 0.015)


# same example but with a Cox model instead
(aids.pscp <- coxph(Surv(survtime + 0.9, status) ~ state + T.categ + 
pspline(age, df=6), data = Aidsp))
zzcp <- predict(aids.pscp, data.frame(state = factor(rep("NSW",
83), levels =
levels(Aidsp$state)),
    T.categ = factor(rep("hs", 83), levels = levels(Aidsp$T.categ)),
age =
0:82), se = T, type = "lp")
plot(0:82, 1/exp(zzcp$fit), type = "l", ylim = c(0, 2), xlab =
"age", ylab =
"expected lifetime (years)")
lines(0:82, 1/exp(zzcp$fit+1.96*zzcp$se.fit), lty = 3, col = 2)
lines(0:82, 1/exp(zzcp$fit-1.96*zzcp$se.fit), lty = 3, col = 2)
rug(Aidsp$age + runif(length(Aidsp$age), -0.5, 0.5), ticksize = 0.015)

1. survreg() does NOT fit a proportional hazards model, a mistake
repeated multiple times in your post

2. The coxph function operates on the risk scale: large values of Xbeta
= large death rates = bad
   The survreg operates on the time scale: large values of xbeta longer
liftetime = good.

3. predict(fit, type='risk') = exp(predict(fit, type='linear'))
in a Cox
model returns an estimate of the relative risk for each subject.  That
is, his/her predicted death rate as compared to the others in the
sample.  It has no units of "years" or "days" or anything
else.  The
predicted survival TIME for a subject is something else entirely.

  predict(fit, type='response') in a survreg model does give predicted
survvival times.  

  If you really want to understand the interrelationships of these
things more deeply I think you need some textbook time.  Read the book
by Kalbfleisch and Prentice for accelerated failure time models, or even
better Escobar and Meeker which comes from the industrial reliability
view.  For predicted survival from a Cox model see Chapter 10 of
Therneau and Grambsch.  The answers to your specific questions would be
a document rather than an email.

Terry Therneau