Software 11.5.1 (Cumulative incidence). In R’s survival package, a factor event whose first level is censoring makes survfit compute the Aalen–Johansen estimates. The mgus2 data, with the time to progression or to death (in months), the estimated probabilities at \(10\) and \(20\) years of being event-free, progressed or dead, and the incorrect one minus Kaplan–Meier for progression:

> library(survival); options(show.signif.stars = FALSE)
> m <- mgus2
> m$etime <- with(m, ifelse(pstat == 1, ptime, futime))      # months
> m$event <- factor(with(m, ifelse(pstat == 1, 1, 2 * death)), 0:2,
+                   labels = c("censor", "pcm", "death"))
> table(m$event)

censor    pcm  death
   409    115    860
> aj <- survfit(Surv(etime, event) ~ 1, data = m)           # Aalen-Johansen
> s <- summary(aj, times = c(120, 240))
> round(cbind(months = s$time, s$pstate), 4)
     months
[1,]    120 0.4045 0.0637 0.5318
[2,]    240 0.1762 0.0998 0.7240
> km <- survfit(Surv(etime, event == "pcm") ~ 1, data = m)  # wrong: 1 - KM
> round(1 - summary(km, times = c(120, 240))$surv, 4)
[1] 0.0952 0.2096

Reading the output. At \(120\) months, \(40\%\) are alive without progression, \(6.4\%\) have progressed and \(53\%\) have died without progression; the rows add to one. The naive estimate of progression, \(9.5\%\), is too high.

Python: lifelines.AalenJohansenFitter, with event_of_interest set to the cause.

Software 11.5.2 (Cause-specific and Fine–Gray models). A cause-specific Cox model for progression, the Fine–Gray model fitted with finegray (which builds the weighted data set) and a weighted coxph, and the cause-specific model for death:

> cs <- coxph(Surv(etime, event == "pcm") ~ age + sex + mspike, data = m)
> fg <- finegray(Surv(etime, event) ~ ., data = m, etype = "pcm")
> sd <- coxph(Surv(fgstart, fgstop, fgstatus) ~ age + sex + mspike,
+             data = fg, weights = fgwt)
> round(rbind(cause_specific = exp(coef(cs)),
+             subdistribution = exp(coef(sd))), 4)
                   age   sexM mspike
cause_specific  1.0165 0.9950 2.4208
subdistribution 0.9833 0.8081 2.4312
> cd <- coxph(Surv(etime, event == "death") ~ age + sex + mspike, data = m)
> round(exp(coef(cd)), 4)                       # cause-specific, for death
   age   sexM mspike
1.0674 1.4754 0.9425

Reading the output. The effect of the spike is the same in both models; the age and sex effects differ, because age and sex act strongly on death, the competing risk.

Python: lifelines fits cause-specific models with CoxPHFitter; the Fine–Gray model is not in lifelines.

Software 11.5.3 (Frailty and robust standard errors). The kidney data with a marginal model and robust standard errors, and with a shared gamma frailty for each patient, with the estimated frailty variance and its test:

> k1 <- coxph(Surv(time, status) ~ age + sex + cluster(id), data = kidney)
> k2 <- coxph(Surv(time, status) ~ age + sex +
+               frailty(id, distribution = "gamma"), data = kidney,
+             control = coxph.control(iter.max = 100))
> round(rbind(marginal = exp(coef(k1)), frailty = exp(coef(k2))), 3)
           age   sex
marginal 1.002 0.436
frailty  1.005 0.204
> round(summary(k1)$coef["sex", c("robust se", "Pr(>|z|)")], 4)
robust se  Pr(>|z|)
   0.4835    0.0863
> round(summary(k2)$coefficients[1:2, c("se(coef)", "p")], 4)
    se(coef)      p
age   0.0119 0.6587
sex   0.4606 0.0006
> round(c(theta = k2$history[[1]]$theta,              # frailty variance
+         p_frailty = summary(k2)$coefficients[3, "p"]), 4)
    theta p_frailty
   0.4078    0.0404

Reading the output. The frailty model’s sex effect is stronger (\(0.20\) against \(0.44\)) and more significant; \(\hat \theta = 0.41\) with \(p = 0.04\) shows the patients’ two infection times are dependent.

Python: CoxPHFitter with cluster_col gives robust errors; gamma frailty is not in lifelines.

Software 11.5.4 (Recurrent events). The cgd trial in counting-process form: the Andersen–Gill and PWP models with robust errors, the number of infections in each arm, and a model for the first infection only:

> ag <- coxph(Surv(tstart, tstop, status) ~ treat + cluster(id), data = cgd)
> pwp <- coxph(Surv(tstart, tstop, status) ~ treat + strata(enum) +
+                cluster(id), data = cgd)
> cols <- c("exp(coef)", "se(coef)", "robust se", "Pr(>|z|)")
> round(rbind(AG = summary(ag)$coef[, cols],
+             PWP = summary(pwp)$coef[, cols]), 4)
    exp(coef) se(coef) robust se Pr(>|z|)
AG     0.3344   0.2610    0.3119   0.0004
PWP    0.4231   0.2802    0.2919   0.0032
> with(cgd, tapply(status, treat, sum))         # infections in each arm
placebo  rIFN-g
     56      20
> first <- coxph(Surv(tstart, tstop, status) ~ treat, data = cgd,
+                subset = enum == 1)            # first infection only
> round(summary(first)$coef[, c("exp(coef)", "Pr(>|z|)")], 4)
exp(coef)  Pr(>|z|)
   0.3349    0.0011

Reading the output. All three analyses put the reduction in infections at \(58\%\) to \(67\%\); the robust standard error is about \(20\%\) larger than the naive one in the Andersen–Gill model.

Python: CoxTimeVaryingFitter with start–stop data in lifelines.