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
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