[R] coxreg vs coxph: time-dependent treatment
Ehsan Karim
wildscop at hotmail.com
Mon Sep 12 04:43:15 CEST 2011
Sorry: there was an error in the weight calculation, fixed version is
the following, but still the final estimates differ as explained in
the original email:
#########################
require(survival)
require(eha)
data(heart)
head(heart)
follow <- heart$stop - heart$start
fit <- glm(transplant ~ age + surgery + year + follow,
data=heart, family = binomial)
heart$wt <- ifelse(heart$transplant == 0,
(1 - predict(fit, type = "response")),
(predict(fit, type = "response")))
heart$iptw <- unlist(tapply(1/heart$wt, heart$id, cumprod))
summary(heart$iptw)
# no weights
fit0 <- coxph(Surv(start,stop,event)~transplant, data=heart)
fit0 # fit with coxph without case-weights
fit1 <- coxreg(Surv(start,stop,event)~transplant, data=heart)
fit1 # fit with coxreg from eha without case-weights
# coxph
fit2 <- coxph(Surv(start,stop,event)~transplant + cluster(id),
data=heart, weights = iptw, robust = T)
fit2 # fit with coxph having robust and cluster option
fit3 <- coxph(Surv(start,stop,event)~transplant + cluster(id),
data=heart, weights = iptw)
fit3 # fit with coxph having cluster option
fit4 <- coxph(Surv(start,stop,event)~transplant,
data=heart, weights = iptw)
fit4 # fit with coxph
# coxreg
fit5 <- coxreg(Surv(start,stop,event)~transplant + cluster(id),
data=heart, weights = iptw)
fit5 # fit with coxreg from eha having cluster option
fit6 <- coxreg(Surv(start,stop,event)~transplant,
data=heart, weights = iptw)
fit6 # fit with coxreg from eha
exp(coef(fit3)) # HR from coxph having cluster option
exp(coef(fit4)) # HR from coxph
exp(coef(fit5))[1] # HR from coxreg having cluster option
exp(coef(fit6))[1] # HR from coxreg
#########################
> exp(coef(fit3)) # HR from coxph having cluster option
transplant1
17.94681
> exp(coef(fit4)) # HR from coxph
transplant1
17.94681
> exp(coef(fit5))[1] # HR from coxreg having cluster option
transplant1
20.06519
> exp(coef(fit6))[1] # HR from coxreg
transplant1
17.94681
#########################
More information about the R-help
mailing list