[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