# load data rhdnase <- read.table("http://publicifsv.sund.ku.dk/~pka/epidata/rhdnase.txt", header = TRUE) library(survival) fit_survival <- survfit(Surv(start, stop, status !=0)~trt, data = subset(rhdnase, etype == 1)) # extract information for Nelson-Aalen estimates na_all_1 <- cumsum(fit_survival[1]$n.event[fit_survival[1]$n.risk !=0] / fit_survival[1]$n.risk[fit_survival[1]$n.risk !=0]) na_all_2 <- cumsum(fit_survival[2]$n.event[fit_survival[2]$n.risk !=0] / fit_survival[2]$n.risk[fit_survival[2]$n.risk !=0]) plot(na_all_1 ~ fit_survival[1]$time[fit_survival[1]$n.risk !=0], type = "s", xlab = "Days", ylab = "Cumulative all-cause hazard", cex = 1.3, lwd = 2) lines(na_all_2 ~ fit_survival[2]$time[fit_survival[2]$n.risk !=0], type = "s", lwd = 2, col = 2) legend("topleft", legend = c("trt 0", "trt 1"), col = c(1, 2), lwd=2, bty = "n", cex = 1.3) # 2 fit1 <- coxph(Surv(start, stop, status !=0) ~ factor(trt==0), data = subset(rhdnase, etype == 1)) summary(fit1) fit2 <- coxph(Surv(start, stop, status !=0) ~ factor(trt==0) + factor(enum1, levels = c(5, 1:4)), data = subset(rhdnase, etype == 1)) summary(fit2) # 3 rhdnase$wait <- rhdnase$stop - rhdnase$start fit3 <- coxph(Surv(wait, status !=0) ~ factor(trt==0), data = subset(rhdnase, etype == 1)) summary(fit3) # almost identical fit4 <- coxph(Surv(wait, status !=0) ~factor(trt==0) + factor(enum1, levels = c(5, 1:4)), data = subset(rhdnase, etype == 1)) summary(fit4) # slightly higher HRs # 4 library(reda) fit5 <- mcf(Survr(id, stop, status !=0, origin = startold, check = FALSE) ~ trt, data = subset(rhdnase, etype == 1)) plot(fit5) # 5 fit6 <- coxph(Surv(startold, stop, status !=0) ~ frailty(id) + factor(trt==0), data = subset(rhdnase, etype == 1)) summary(fit6)