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

## ----setup--------------------------------------------------------------------
library(spsurv)
library(KMsurv)
library(survival)
library(ggplot2)
library(generics)
data(larynx)
larynx$stage <- factor(larynx$stage)

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

## ----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")

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

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

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

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

## ----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")

## ----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()

## ----print--------------------------------------------------------------------
print(fit, what = "summary")