Software 7.5.1 (The partial likelihood by hand). Example 7.2.1: the log partial likelihood written as a sum over events of \(\beta x_i - \log \sum _{\text {risk set}}e^{\beta x_l}\), maximised with optimize; then coxph:

> library(survival)
> tm <- c(2, 3, 5, 6, 8, 9, 11, 14)           # eight subjects, no ties
> ev <- c(1, 1, 0, 1, 1, 0, 1, 1)
> x  <- c(1, 0, 1, 1, 0, 1, 0, 0)             # a binary covariate
> logPL <- function(b) {
+   s <- 0
+   for (i in which(ev == 1)) {               # one term per event
+     risk <- tm >= tm[i]                     # the risk set at that time
+     s <- s + b * x[i] - log(sum(exp(b * x[risk])))
+   }
+   s
+ }
> o <- optimize(logPL, c(-5, 5), maximum = TRUE)
> round(c(beta = o$maximum, logPL = o$objective), 4)
   beta   logPL
 0.4442 -7.6200
> round(coef(coxph(Surv(tm, ev) ~ x)), 4)    # coxph agrees
     x
0.4442

Reading the output. The maximum is at \(\hat \beta = 0.4442\), and coxph gives the same value. The baseline hazard never appears: only the covariates of those at risk at each event time matter.

Python: lifelines.CoxPHFitter().fit(df, ’time’, ’event’).

Software 7.5.2 (Sex as a predictor). Example 7.3.1: the coefficient table, the hazard ratio with its interval, the three tests of \(\beta = 0\), and the log-rank statistic for comparison:

> fit <- coxph(Surv(time, status) ~ sex, data = lung)
> round(summary(fit)$coefficients, 4)
      coef exp(coef) se(coef)       z Pr(>|z|)
sex -0.531     0.588   0.1672 -3.1764   0.0015
> round(summary(fit)$conf.int, 4)             # HR for women vs men, and CI
    exp(coef) exp(-coef) lower .95 upper .95
sex     0.588     1.7007    0.4237     0.816
> round(c(wald = summary(fit)$waldtest[1], lr = summary(fit)$logtest[1],
+         score = summary(fit)$sctest[1]), 3)
 wald.test    lr.test score.test
    10.090     10.634     10.325
> round(survdiff(Surv(time, status) ~ sex, data = lung)$chisq, 3)  # log-rank
[1] 10.327

Reading the output. \(\hat \beta = -0.531\), HR \(0.588\) (\(0.424\) to \(0.816\)), \(p = 0.0015\) by the Wald test. exp(-coef) \(= 1.70\) is the hazard ratio the other way round: men’s hazard was \(1.7\) times women’s. The score test, \(10.325\), is the log-rank test up to the treatment of ties.

Python: cph.print_summary() shows exp(coef) and its interval.

Software 7.5.3 (Several covariates). Example 7.3.2: hazard ratios, intervals and \(p\)-values for age, sex and ECOG score, and the hazard ratio for ten years of age:

> lg <- na.omit(lung[, c("time", "status", "age", "sex", "ph.ecog")])
> m <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lg)
> round(cbind(HR = exp(coef(m)), exp(confint(m)),
+             p = summary(m)$coefficients[, 5]), 4)
            HR  2.5 % 97.5 %      p
age     1.0111 0.9929 1.0297 0.2324
sex     0.5754 0.4142 0.7994 0.0010
ph.ecog 1.5900 1.2727 1.9864 0.0000
> round(exp(10 * coef(m)["age"]), 4)          # HR for ten years of age
  age
1.117

Reading the output. After adjustment, each point of ECOG score multiplies the hazard by \(1.59\), and being female by \(0.58\); ten years of age multiply it by \(1.12\), not significantly. These are the values of Figure 7.

Python: the same model with CoxPHFitter and a formula ’age + sex + ph_ecog’.