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

David Winsemius dwinsemius at comcast.net
Thu Nov 25 15:14:20 CET 2010


On Nov 25, 2010, at 7:27 AM, Ben Rhelp wrote:

> 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

Different parameterization. You can find expanded answer(s) in the  
archives and in the documentation of survreg.distributions.


> and why we do not need to divide by the number of days in the year
> anymore?

Here I'm guessing (since you don't offer enough evidence to confirm)  
that the difference is in the time scales used in your Aidsp$survtime  
versus some other example to which you are comparing .

 > data(Aidsp)
Warning message:
In data(Aidsp) : data set 'Aidsp' not found

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


David Winsemius, MD
West Hartford, CT



More information about the R-help mailing list