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

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

## ----tidy---------------------------------------------------------------------
td <- tidy(fit1, conf.int = TRUE, exponentiate = TRUE)
td

## ----glance-------------------------------------------------------------------
glance(fit1)

## ----forest, fig.cap = "Karnofsky hazard ratio with 95% CI."------------------
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()

## ----criterion-compare, fig.cap = "AIC and log-likelihood: null vs Karnofsky model."----
cmp <- data.frame(
  model = c("Null", "Karnofsky"),
  AIC = c(AIC(fit0), AIC(fit1)),
  logLik = c(as.numeric(logLik(fit0)), as.numeric(logLik(fit1))),
  df = c(attr(logLik(fit0), "df"), attr(logLik(fit1), "df"))
)
cmp

## ----anova--------------------------------------------------------------------
anova(fit0, fit1)
anova(fit1, fit2)
anova(fit2)

## ----df-table-----------------------------------------------------------------
df_tbl <- data.frame(
  Quantity = c(
    "logLik() / AIC df",
    "summary() logtest df",
    "anova() nested df diff"
  ),
  fit0 = c(
    attr(logLik(fit0), "df"),
    "\u2014",
    "\u2014"
  ),
  fit1 = c(
    attr(logLik(fit1), "df"),
    summary(fit1)$logtest["df"],
    attr(logLik(fit1), "df") - attr(logLik(fit0), "df")
  ),
  Counts = c(
    "Regression + length(bp.param)",
    "Regression terms only",
    "Total parameter difference"
  )
)
df_tbl

## ----confint------------------------------------------------------------------
confint(fit1, level = 0.95)

## ----vcov---------------------------------------------------------------------
dim(vcov(fit1))
dim(vcov(fit1, bp.param = TRUE))

## ----estimates----------------------------------------------------------------
head(estimates(fit1))
extractAIC(fit1)

## ----null---------------------------------------------------------------------
summary(fit0)
logLik(fit0)