---
title: "Getting started with spsurv"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with spsurv}
  %\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** fits semi-parametric survival regression for right-censored data using
Bernstein-polynomial (BP) baselines. Three model families share one interface:

- **PH** — proportional hazards (`bpph`, `model = "ph"`)
- **PO** — proportional odds (`bppo`, `model = "po"`)
- **AFT** — accelerated failure time (`bpaft`, `model = "aft"`)

## Before you model

```{r setup}
library(spsurv)
library(KMsurv)
library(survival)
library(ggplot2)
library(generics)
data(larynx)
larynx$stage <- factor(larynx$stage)
```

```{r eda-censoring}
censor_tbl <- do.call(rbind, lapply(split(larynx, larynx$stage), function(d) {
  data.frame(
    stage = as.character(d$stage[1]),
    n = nrow(d),
    events = sum(d$delta),
    censored = sum(1 - d$delta),
    stringsAsFactors = FALSE
  )
}))
censor_tbl$pct_censored <- round(100 * censor_tbl$censored / censor_tbl$n, 1)
censor_tbl
```

```{r eda-km, fig.cap = "Kaplan-Meier survival by larynx cancer stage."}
km_stage <- survfit(Surv(time, delta) ~ stage, data = larynx)
km_long <- data.frame(
  time = km_stage$time,
  surv = km_stage$surv,
  stage = rep(levels(larynx$stage), km_stage$strata)
)
ggplot(km_long, aes(x = time, y = surv, color = stage)) +
  geom_step(linewidth = 0.6) +
  labs(x = "Time (years)", y = "Survival probability", color = "Stage") +
  theme_bw() +
  theme(legend.position = "bottom")
```

### Why this is insightful

Stage separation in the KM plot motivates including `stage` as a covariate.
Censoring percentages set expectations for precision — heavy censoring or tiny
stage groups mean wider confidence intervals.

### What to decide

- Clear stage ordering → PH with `stage` is a plausible starting point.
- Heavy censoring or small groups → interpret CIs cautiously; consider simpler
  models or Bayesian analysis.
- If KM curves cross → see the *Choosing PH, PO, and AFT models* vignette before
  fitting Cox PH blindly.

## Fit and summarize

```{r fit}
fit <- bpph(
  Surv(time, delta) ~ age + stage,
  degree = 5,
  data = larynx,
  approach = "mle",
  init = 0
)
summary(fit)
```

The generic `spbp()` function is equivalent when you pass `model` explicitly:

```{r spbp}
fit2 <- spbp(
  Surv(time, delta) ~ age + stage,
  degree = 5,
  data = larynx,
  model = "ph",
  approach = "mle",
  init = 0
)
```

## Specifying the Bernstein degree

```{r bernstein-spec}
fit3 <- bpph(
  Surv(time, delta) ~ age + stage,
  data = larynx,
  approach = "mle",
  dist = bernstein(5),
  init = 0
)
length(fit3$bp.param)
```

## The `spbp` object

```{r object-parts}
names(fit)[names(fit) %in% c("coefficients", "bp.param", "n", "nevent")]
fit$call$model
fit$call$approach
```

| Component | Meaning |
|-----------|---------|
| `coefficients` | Regression estimates (on original covariate scale) |
| `bp.param` | Bernstein baseline coefficients (`gamma`) |
| `degree` | Bernstein polynomial degree used in the fit |
| `bp.param` length | Same as `degree` |
| `call$model` | `"ph"`, `"po"`, or `"aft"` |
| `call$approach` | `"mle"` or `"bayes"` |

## Visual check — KM (step) vs Bernstein survival (smooth)

Kaplan--Meier is nonparametric, so it is drawn with **steps**. Bernstein
polynomial survival is a continuous function of time, so model curves use
**lines** (`ggplot2::geom_line` or `plot(survfit(fit))`).

```{r km-bp-overlay, fig.cap = "Kaplan-Meier by stage (steps) vs Bernstein PH at median age (smooth dashed)."}
newdata <- data.frame(
  age = median(larynx$age),
  stage = factor(levels(larynx$stage), levels = levels(larynx$stage))
)
plot_times <- seq(0, max(larynx$time), length.out = 121)
pr <- predict(fit, newdata = newdata, times = plot_times)
pr$stage <- newdata$stage[match(as.character(pr$id), as.character(seq_len(nrow(newdata))))]

ggplot() +
  geom_step(
    data = km_long,
    aes(x = time, y = surv, color = stage),
    linewidth = 0.5
  ) +
  geom_line(
    data = pr,
    aes(x = time, y = surv, color = stage),
    linetype = "dashed",
    linewidth = 0.7
  ) +
  labs(x = "Time (years)", y = "Survival probability", color = "Stage") +
  theme_bw() +
  theme(legend.position = "bottom")
```

## Visual check — forest plot

```{r forest, fig.cap = "Exponentiated coefficients (hazard ratios) with 95% CIs."}
td <- tidy(fit, conf.int = TRUE, exponentiate = TRUE)
td$term <- factor(td$term, levels = rev(td$term))
ggplot(td, aes(x = estimate, y = term, xmin = conf.low, xmax = conf.high)) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey50") +
  geom_pointrange(linewidth = 0.4) +
  labs(x = "Hazard ratio", y = NULL) +
  theme_bw()
```

### Why this is insightful

Forest plots translate coefficients into decision-relevant effect sizes. Stage
effects relative to the reference level show how much each stage shifts hazard
after adjusting for age.

### What to decide

- CI excludes 1 → evidence of effect at conventional levels.
- Large age effect → clinically meaningful covariate; keep in model.
- For more covariate-specific ribbons → see the *Survival
  prediction and ggplot* vignette.

## Printing

```{r print}
print(fit, what = "summary")
```

## Next steps

- *Choosing PH, PO, and AFT models* — compare regression families
- *Survival prediction and ggplot* — smooth survival curves
- *Bayesian analysis with Stan* — priors and posterior summaries