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)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
