On Nov 25, 2010, at 10:08 AM, Ben Rhelp wrote:
Hi David,
Thank you for your reply. See below for more information.
From: David Winsemius
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.
I understand (i think) the difference in model structures between a
Cox (coxph)
and Proportional hazard model (survreg).
I couldn't tell whether this means you decided that those citations
answered your question. If not, then refer to Therneau's or Lumley's
replies in rhelp to a similar question earlier this month.:
https://stat.ethz.ch/pipermail/r-help/2010-November/259796.html
https://stat.ethz.ch/pipermail/r-help/2010-November/259747.html
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 .
Both models are run from the same data, so I am not expecting any
differences in
time scales.
To get similar results, I need to actually run the following
equations:
expected_lifetime_in_years = exp(fit)/365.25 --->
Linear
prediction of the Proportional hazard model
expected_lifetime_in_years = 1/exp(fit)
---> Linear
prediction of the Cox model
where fit come from the linear prediction of each models,
respectively.
Actually, in the code below, I re-run the models and predictions
based on a
yearly sampling time (instead of daily).
Again, to get similar results, I now need to actually run the
following
equations:
expected_lifetime_in_years = exp(fit)
---> Linear
prediction of the Proportional hazard model
expected_lifetime_in_years = 1/exp(fit)
---> Linear
prediction of the Cox model
I think I understand the logic behind the results of the
proportional hazard
model, but not for the prediction of the Cox model.
Cox models are PH models.
Thank you for your help. I hope this is not a too stupid hole in my
logic.
Here is the self contained R code to produce the charts:
library(MASS);
library(survival);
#Same data but parametric fit
make.aidsp <- function(){
cutoff <- 10043 # July 1987 in julian days
btime <- pmin(cutoff, Aids2$death) - pmin(cutoff, Aids2$diag)
atime <- pmax(cutoff, Aids2$death) - pmax(cutoff, Aids2$diag)
survtime <- btime + 0.5*atime
status <- as.numeric(Aids2$status)
data.frame(survtime, status = status - 1, state = Aids2$state,
T.categ = Aids2$T.categ, age = Aids2$age, sex = Aids2$sex)
}
Aidsp <- make.aidsp()
# 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)
# Change the sampling time from daily to yearly
par(mfrow = c(1, 1));
(aids.ps <- survreg(Surv((survtime + 0.9)/365.25, 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), type = "l", ylim = c(0, 2), xlab = "age",
ylab =
"expected lifetime (years)")
(aids.pscp <- coxph(Surv((survtime + 0.9)/365.25, 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")
lines(0:82, 1/exp(zzcp$fit),col=4)
David Winsemius, MD
West Hartford, CT
______________________________________________
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.
