Skip to contents

Choosing PH, PO, and AFT models

Source: vignettes/model-families.Rmd
model-families.Rmd

spsurv implements three standard survival regression frameworks under one likelihood and syntax:

Family Function Covariate effect Typical interpretation
PH bpph() Multiplicative on hazard Hazard ratio
PO bppo() Multiplicative on cumulative odds Odds ratio
AFT bpaft() Additive on log-time scale Time ratio

Before you model 

library(spsurv)
library(survival)
library(ggplot2)
data(veteran)
veteran$celltype <- factor(veteran$celltype)
km_ct <- survfit(Surv(time, status) ~ celltype, data = veteran)
km_long <- data.frame(
  time = km_ct$time,
  surv = km_ct$surv,
  celltype = rep(levels(veteran$celltype), km_ct$strata)
)
ggplot(km_long, aes(x = time, y = surv, color = celltype)) +
  geom_step(linewidth = 0.6) +
  labs(x = "Time (days)", y = "Survival probability", color = "Cell type") +
  theme_bw() +
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 45, hjust = 1))
Kaplan-Meier by cell type.

Kaplan-Meier by cell type. 

ggplot(veteran, aes(x = factor(status), y = karno, fill = factor(status))) +
  geom_violin(alpha = 0.6) +
  geom_boxplot(width = 0.15, outlier.alpha = 0.4) +
  scale_x_discrete(labels = c("Censored", "Event")) +
  labs(x = NULL, y = "Karnofsky score", fill = NULL) +
  theme_bw() +
  theme(legend.position = "none")
Karnofsky score by event status.

Karnofsky score by event status.

The paper examples sometimes restrict to prior == 0; we use the full veteran dataset here for simplicity.

Why this is insightful

KM by cell type reveals crossing and separation patterns. Karnofsky differs between events and censoring — a strong candidate covariate.

What to decide

  • Crossing KM curves → PO or AFT may fit better than PH.
  • Karnofsky separation → include in all three family fits for fair comparison.

Fit and summarize 

fml <- Surv(time, status) ~ karno + celltype
deg <- 5L
fit_ph <- bpph(fml, degree = deg, data = veteran, approach = "mle", init = 0)
fit_po <- bppo(fml, degree = deg, data = veteran, approach = "mle", init = 0)
fit_aft <- bpaft(fml, degree = 3, data = veteran, approach = "mle", init = 0)

AFT uses degree = 3 here so the Bernstein gamma block stays numerically stable on veteran; PH and PO use degree = 5. The overlay below plots point estimates of survival; when gamma_information_stable is FALSE, bands are still computed but survfit() warns that they may be unreliable.

summary(fit_ph)$coefficients[, c("coef", "exp(coef)", "Pr(>|z|)")]
#>                          coef exp(coef)     Pr(>|z|)
#> karno             -0.03061037 0.9698534 1.951228e-09
#> celltypesmallcell  0.73804877 2.0918498 3.294217e-03
#> celltypeadeno      1.15357958 3.1695182 8.086335e-05
#> celltypelarge      0.33021792 1.3912713 2.317345e-01

Visual check — survival overlay

Same covariate profiles under PH, PO, and AFT:

newdata <- data.frame(karno = c(30, 70), celltype = "squamous")
plot_times <- seq(0, max(veteran$time), length.out = 121)

pr_ph <- predict(fit_ph, newdata = newdata, times = plot_times)
pr_po <- predict(fit_po, newdata = newdata, times = plot_times)
pr_aft <- predict(fit_aft, newdata = newdata, times = plot_times)

pr_ph$family <- "PH"
pr_po$family <- "PO"
pr_aft$family <- "AFT"
pr_all <- rbind(pr_ph, pr_po, pr_aft)
pr_all$label <- paste0(
  pr_all$family, " (Karnofsky ",
  newdata$karno[match(as.character(pr_all$id), as.character(seq_len(nrow(newdata))))],
  ")"
)

ggplot(pr_all, aes(x = time, y = surv, color = label, linetype = family)) +
  geom_line(linewidth = 0.5) +
  labs(x = "Time (days)", y = "Survival probability", color = NULL, linetype = "Family") +
  theme_bw() +
  theme(legend.position = "bottom")
Predicted survival: PH vs PO vs AFT (Karnofsky 30 and 70, squamous).

Predicted survival: PH vs PO vs AFT (Karnofsky 30 and 70, squamous). 

aic_tbl <- data.frame(
  model = c("PH", "PO", "AFT"),
  AIC = c(AIC(fit_ph), AIC(fit_po), AIC(fit_aft))
)
aic_tbl
#>   model      AIC
#> 1    PH 1450.026
#> 2    PO 1438.685
#> 3   AFT 1464.954

Why this is insightful

The same covariates on different regression scales produce different survival curves. AIC summarizes fit but does not replace clinical judgment when curves diverge meaningfully.

What to decide

  • PH vs PO curves similar → PH may suffice for reporting hazard ratios.
  • Large divergence or known crossing in KM → prefer PO or check diagnostics.
  • AFT when time ratios are the estimand of interest.
  • Lower AIC ≠ automatic winner if curves disagree clinically.

When to use which family

  • PH — constant hazard ratios over time.
  • PO — proportionality on cumulative odds; useful when hazard ratios are not stable.
  • AFT — direct interpretation on the time scale.

All three share the same Bernstein baseline machinery and Stan backend.