## ----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(survival)
library(ggplot2)
data(veteran)

## ----eda-event-times, fig.cap = "Distribution of observed event times."-------
ggplot(subset(veteran, status == 1), aes(x = time)) +
  geom_histogram(bins = 25, fill = "steelblue", color = "white", alpha = 0.8) +
  labs(x = "Time (days)", y = "Count", title = "Event times only") +
  theme_bw()

## ----fit-degrees--------------------------------------------------------------
fit_low <- bpph(
  Surv(time, status) ~ karno,
  degree = 3,
  data = veteran,
  approach = "mle",
  init = 0
)
fit_high <- bpph(
  Surv(time, status) ~ karno,
  degree = 8,
  data = veteran,
  approach = "mle",
  init = 0
)
c(low = length(fit_low$bp.param), high = length(fit_high$bp.param))

## ----aic-compare--------------------------------------------------------------
AIC(fit_low, fit_high)

## ----baseline-overlay, fig.cap = "Cumulative baseline hazard: degree 3 vs degree 8."----
times <- seq(0, max(veteran$time), length.out = 100)
tau <- max(veteran$time)

baseline_curve <- function(fit, times, tau) {
  bb <- bp.basis(times, degree = length(fit$bp.param), tau = tau)
  data.frame(
    time = times,
    cumhaz = as.vector(bb$G %*% fit$bp.param),
    degree = paste0("m = ", length(fit$bp.param))
  )
}

bl_df <- rbind(
  baseline_curve(fit_low, times, tau),
  baseline_curve(fit_high, times, tau)
)

ggplot(bl_df, aes(x = time, y = cumhaz, color = degree)) +
  geom_line(linewidth = 0.7) +
  labs(x = "Time (days)", y = "Cumulative baseline hazard", color = NULL) +
  theme_bw()

## ----plot-baseline, fig.cap = "Estimated baseline hazard components (degree 8)."----
plot(fit_high, graph = "baseline", cumulative = FALSE)

## ----plot-basis, fig.cap = "Bernstein basis functions on [0, tau]."-----------
plot(fit_high, graph = "basis", cumulative = FALSE)

## ----kappa-table--------------------------------------------------------------
kappa_tbl <- data.frame(
  degree = c(3, 8),
  gamma_information_stable = c(
    attr(vcov(fit_low, bp.param = TRUE), "gamma_information_stable"),
    attr(vcov(fit_high, bp.param = TRUE), "gamma_information_stable")
  ),
  kappa_gamma = c(
    signif(attr(vcov(fit_low, bp.param = TRUE), "gamma_information_kappa"), 4),
    signif(attr(vcov(fit_high, bp.param = TRUE), "gamma_information_kappa"), 4)
  )
)
kappa_tbl

## ----bp-basis-----------------------------------------------------------------
times <- seq(0, max(veteran$time), length.out = 50)
bb <- bp.basis(times, degree = 5, tau = max(veteran$time))
dim(bb$g)

## ----pw-basis-----------------------------------------------------------------
pw <- pw.basis(degree = 5)
dim(pw)