9.4 Checking the assumption in software
Software 9.4.1 (Tests of proportional hazards). cox.zph for the veteran model: a test for each covariate (with \(3\) degrees of freedom for the four-level cell type) and a global test:
> library(survival); options(show.signif.stars = FALSE)
> vet <- veteran; vet$trt <- factor(vet$trt, labels = c("standard", "test"))
> m1 <- coxph(Surv(time, status) ~ trt + celltype + karno, data = vet)
> z <- cox.zph(m1, transform = "km") # tests of proportional hazards
> round(z$table, 4)
chisq df p
trt 0.2063 1 0.6497
celltype 14.5317 3 0.0023
karno 12.8133 1 0.0003
GLOBAL 22.8597 5 0.0004
Reading the output. Treatment satisfies the assumption (\(p = 0.65\)); cell type (\(p = 0.002\)) and the Karnofsky score (\(p = 0.0003\)) do not, and the global test rejects proportional hazards for the model (\(p = 0.0004\)). The treatment comparison, the trial’s purpose, is still valid if the model is fixed for the other covariates, for example by stratifying by cell type and letting the Karnofsky effect vary with time.
Python: cph.check_assumptions(df) in lifelines.
Software 9.4.2 (A Schoenfeld residual by hand). The Schoenfeld residual for the Karnofsky score at the first death time shared by no other death, computed as the dying patient’s score minus the risk-weighted mean score of those at risk, against the value from residuals with type set to schoenfeld; then the slope of the scaled residuals on log time:
> r <- residuals(m1, type = "schoenfeld")[, "karno"] > dt <- sort(vet$time[vet$status == 1]) # death times, in order > k <- which(!(dt %in% dt[duplicated(dt)]))[1] # the first untied one > tj <- dt[k]; i <- which(vet$time == tj & vet$status == 1) > risk <- vet$time >= tj; lp <- predict(m1, type = "lp") > w <- exp(lp[risk]) / sum(exp(lp[risk])) # weights in the risk set > c(time = tj, by_hand = vet$karno[i] - sum(w * vet$karno[risk]), + coxph = unname(r[k])) time by_hand coxph 2.000000 -4.624778 -4.624778 > sr <- residuals(m1, type = "scaledsch")[, "karno"] # beta(t), roughly > round(coef(lm(sr ~ log(dt))), 5) # a trend in beta(t)? (Intercept) log(dt) -0.07203 0.01004 > round(c(beta_hat = unname(coef(m1)["karno"])), 5) beta_hat -0.03127
Reading the output. At day \(2\) the patient who died had a Karnofsky score \(4.62\) points below the weighted average of those at risk, the value R reports. The scaled residuals rise by \(0.010\) per unit of log time: \(\beta (t)\) increases, from strongly negative towards zero, while the fitted constant is \(-0.031\).
Python: cph.compute_residuals(df, ’schoenfeld’).
Software 9.4.3 (Time-dependent coefficients). The follow-up split at \(100\) days with survSplit, a separate Karnofsky coefficient in each period (shown as hazard ratios for \(10\) points, with \(p\)-values); then a coefficient varying as \(\beta + \gamma \log t\) with tt():
> vs <- survSplit(Surv(time, status) ~ ., data = vet, cut = 100,
+ episode = "period") # before and after 100 days
> vs$karno_early <- vs$karno * (vs$period == 1)
> vs$karno_late <- vs$karno * (vs$period == 2)
> m2 <- coxph(Surv(tstart, time, status) ~ trt + celltype + karno_early +
+ karno_late, data = vs)
> round(exp(10 * coef(m2)[c("karno_early", "karno_late")]), 4) # 10 points
karno_early karno_late
0.6508 1.0427
> round(summary(m2)$coefficients[c("karno_early", "karno_late"), 5], 6)
karno_early karno_late
0.000000 0.689911
> m3 <- coxph(Surv(time, status) ~ trt + celltype + karno + tt(karno),
+ data = vet, tt = function(x, t, ...) x * log(t))
> round(summary(m3)$coefficients[c("karno", "tt(karno)"), c(1, 5)], 5)
coef Pr(>|z|)
karno -0.08087 0.00000
tt(karno) 0.01350 0.00175
Reading the output. Before \(100\) days, \(10\) points lower the hazard by a third (HR \(0.65\)); after, there is no effect (HR \(1.04\), \(p = 0.69\)). In the smooth model, \(\hat \gamma = 0.0135\), \(p = 0.002\), so the effect shrinks steadily over time.
Python: lifelines.CoxTimeVaryingFitter with data in start–stop form.
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.