[R] coxreg vs coxph: time-dependent treatment

Ehsan Karim wildscop at hotmail.com
Mon Sep 12 04:31:46 CEST 2011


Dear List,

After including cluster() option the coxreg (from eha package)
produces results slightly different than that of coxph (from survival)
in the following time-dependent treatment effect calculation (example
is used just to make the point). Will appreciate any explaination /
comment.

cheers,

Ehsan


############################
require(survival)
require(eha)
data(heart)

# create weights
follow <- heart$stop - heart$start
fit <- glm(transplant ~ age + surgery + year + follow,
  data=heart, family = binomial)
heart$iptw <- ifelse(heart$transplant == 0,
  1 - predict(fit, type = "response"),
  predict(fit, type = "response"))
summary(heart$iptw)

# no weights (basic calculation)
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
transplant1
  0.3782417
> exp(coef(fit4))    # HR from coxph
transplant1
  0.3782417
> exp(coef(fit5))[1] # HR from coxreg having cluster option
transplant1
  0.3645507
> exp(coef(fit6))[1] # HR from coxreg
transplant1
  0.3782417
############################



More information about the R-help mailing list