## ----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)
veteran$celltype <- factor(veteran$celltype)

## ----eda-summary--------------------------------------------------------------
eda <- data.frame(
  n = nrow(veteran),
  events = sum(veteran$status),
  censored = sum(1 - veteran$status),
  pct_censored = round(100 * mean(1 - veteran$status), 1)
)
eda

## ----eda-km, fig.cap = "Kaplan-Meier survival for the veteran data (overall)."----
km_all <- survfit(Surv(time, status) ~ 1, data = veteran)
km_df <- data.frame(
  time = km_all$time,
  surv = km_all$surv
)
ggplot(km_df, aes(x = time, y = surv)) +
  geom_step(linewidth = 0.6) +
  labs(x = "Time (days)", y = "Survival probability", title = "Overall KM") +
  theme_bw()

## ----fit----------------------------------------------------------------------
fit <- bpph(
  Surv(time, status) ~ karno + celltype,
  degree = 5,
  data = veteran,
  approach = "mle",
  init = 0
)
summary(fit)

## ----martingale-data----------------------------------------------------------
mart <- residuals(fit, type = "martingale")
csnell <- residuals(fit, type = "cox-snell")
fitted_ch <- csnell

## ----martingale-plot, fig.cap = "Martingale residuals vs fitted cumulative hazard (Cox-Snell)."----
mart_df <- data.frame(
  fitted = fitted_ch,
  martingale = mart,
  status = factor(veteran$status),
  celltype = veteran$celltype
)
ggplot(mart_df, aes(x = fitted, y = martingale, color = status)) +
  geom_point(alpha = 0.6, size = 1.8) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey40") +
  geom_smooth(aes(group = 1), method = "loess", se = FALSE, color = "black", linewidth = 0.5) +
  coord_cartesian(ylim = c(-4, 2)) +
  labs(
    x = "Fitted cumulative hazard (Cox-Snell)",
    y = "Martingale residual",
    color = "Event"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

## ----martingale-facet, fig.cap = "Martingale residuals by cell type."---------
ggplot(mart_df, aes(x = fitted, y = martingale)) +
  geom_point(alpha = 0.6, size = 1.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey40") +
  geom_smooth(method = "loess", se = FALSE, color = "black", linewidth = 0.5) +
  facet_wrap(~ celltype, scales = "free_x") +
  coord_cartesian(ylim = c(-4, 2)) +
  labs(x = "Fitted cumulative hazard", y = "Martingale residual") +
  theme_bw()

## ----coxsnell-calibration, fig.cap = "Cox-Snell calibration: KM of exp(-CS) vs Exponential(1) reference."----
km_cs <- survfit(Surv(csnell, veteran$status) ~ 1)
cs_df <- data.frame(
  cs = km_cs$time,
  cumhaz = -log(pmax(km_cs$surv, .Machine$double.eps))
)
ggplot(cs_df, aes(x = cs, y = cumhaz)) +
  geom_step(linewidth = 0.6) +
  geom_abline(slope = 1, intercept = 0, color = "steelblue", linetype = "dashed") +
  labs(
    x = "Cox-Snell residual",
    y = "Nelson-Aalen cumulative hazard of exp(-CS)",
    title = "Cox-Snell calibration"
  ) +
  theme_bw()

## ----deviance-plot, fig.cap = "Deviance residuals (largest magnitudes highlighted)."----
dev <- residuals(fit, type = "deviance")
dev_df <- data.frame(
  index = seq_along(dev),
  deviance = dev,
  highlight = abs(dev) >= quantile(abs(dev), 0.95)
)
ggplot(dev_df, aes(x = index, y = deviance, color = highlight)) +
  geom_point(size = 1.8, alpha = 0.7) +
  scale_color_manual(values = c("FALSE" = "grey40", "TRUE" = "firebrick"), guide = "none") +
  labs(x = "Observation index", y = "Deviance residual") +
  theme_bw()

## ----cox-compare, fig.cap = "spsurv vs Cox martingale residuals."-------------
cox_fit <- coxph(Surv(time, status) ~ karno + celltype, data = veteran)
mart_cox <- residuals(cox_fit, type = "martingale")
ggplot(data.frame(bp = mart, cox = mart_cox), aes(x = bp, y = cox)) +
  geom_point(alpha = 0.5, size = 1.8) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey40") +
  labs(x = "spsurv martingale", y = "Cox martingale") +
  theme_bw()

## ----vcov-stability-----------------------------------------------------------
v <- vcov(fit, bp.param = TRUE)
stab <- data.frame(
  gamma_information_stable = attr(v, "gamma_information_stable"),
  kappa_gamma = signif(attr(v, "gamma_information_kappa"), 4)
)
stab