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 ######################### ______________________________________________ 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.
