---
title: "Inference and model comparison (MLE)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Inference and model comparison (MLE)}
  %\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
)
```

Maximum-likelihood fits support classical inference on the **full Bernstein
likelihood** (not Cox partial likelihood).

## Before you model

```{r setup}
library(spsurv)
library(generics)
library(survival)
library(ggplot2)
data(veteran)
```

We compare a null model (intercept-only baseline) with a model including Karnofsky
score — the same question as in the paper: is performance status worth including?

## Fit and summarize

```{r 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
)
```

```{r tidy}
td <- tidy(fit1, conf.int = TRUE, exponentiate = TRUE)
td
```

```{r glance}
glance(fit1)
```

## Visual check — forest plot

```{r 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()
```

## Visual check — model comparison

```{r 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
```

```{r anova}
anova(fit0, fit1)
anova(fit1, fit2)
anova(fit2)
```

Pairwise comparisons report the likelihood-ratio statistic on the simpler model
row. A single fitted model produces a sequential table (`NULL`, then each term
in formula order), matching `survival::survreg` conventions.

### Why this is insightful

Forest plots translate coefficients into decision-relevant effect sizes. The LR
test and AIC answer whether Karnofsky improves fit enough to justify the extra
regression parameter (Bernstein degree held constant).

### What to decide

- HR CI excludes 1 and significant LR → include covariate.
- AIC favors complex model but effect negligible → balance parsimony vs fit.
- Compare models at the **same Bernstein degree** — `df` includes all `gamma`
  parameters.

## Degrees of freedom (visual)

```{r 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
```

## Confidence intervals and variance

```{r confint}
confint(fit1, level = 0.95)
```

```{r vcov}
dim(vcov(fit1))
dim(vcov(fit1, bp.param = TRUE))
```

## `estimates()`, `se()`, and `extractAIC()`

```{r estimates}
head(estimates(fit1))
extractAIC(fit1)
```

## Null model

```{r null}
summary(fit0)
logLik(fit0)
```