[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