---
title: "Choosing PH, PO, and AFT models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Choosing PH, PO, and AFT models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  dpi = 96,
  message = FALSE,
  warning = FALSE
)
```

`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

```{r setup}
library(spsurv)
library(survival)
library(ggplot2)
data(veteran)
veteran$celltype <- factor(veteran$celltype)
```

```{r eda-km, fig.cap = "Kaplan-Meier by cell type."}
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))
```

```{r eda-karno, fig.cap = "Karnofsky score by event status."}
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")
```

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

```{r formula}
fml <- Surv(time, status) ~ karno + celltype
deg <- 5L
```

```{r fit-ph}
fit_ph <- bpph(fml, degree = deg, data = veteran, approach = "mle", init = 0)
```

```{r fit-po}
fit_po <- bppo(fml, degree = deg, data = veteran, approach = "mle", init = 0)
```

```{r fit-aft}
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.

```{r summaries}
summary(fit_ph)$coefficients[, c("coef", "exp(coef)", "Pr(>|z|)")]
```

## Visual check — survival overlay

Same covariate profiles under PH, PO, and AFT:

```{r survival-overlay, fig.cap = "Predicted survival: PH vs PO vs AFT (Karnofsky 30 and 70, squamous)."}
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) +  # Bernstein survival is continuous; do not use geom_step
  labs(x = "Time (days)", y = "Survival probability", color = NULL, linetype = "Family") +
  theme_bw() +
  theme(legend.position = "bottom")
```

```{r aic-compare}
aic_tbl <- data.frame(
  model = c("PH", "PO", "AFT"),
  AIC = c(AIC(fit_ph), AIC(fit_po), AIC(fit_aft))
)
aic_tbl
```

### 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.