---
title: "Complete Clinical Analysis Walkthrough"
subtitle: "From Kaplan-Meier baseline to validated multivariable model"
author: "John Ehrlinger"
date: last-modified
format:
  html:
    code-fold: false
    toc: true
    toc-depth: 3
    number-sections: true
vignette: >
  %\VignetteIndexEntry{Complete Clinical Analysis Walkthrough}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---

```{r}
#| include: false
# R CMD build installs the package into R_BUILD_TEMPLIB. A build launched
# through callr (rcmdcheck, devtools) passes Quarto's R session a profile
# that resets .libPaths() without it, and every chunk skipped (#276).
templib <- Sys.getenv("R_BUILD_TEMPLIB")
if (dir.exists(templib)) .libPaths(c(templib, .libPaths()))
has_pkg <- requireNamespace("TemporalHazard", quietly = TRUE)
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
# On CI a package that cannot load must fail the build, not skip every
# chunk and pass over nothing (#276). loadNamespace() says why it failed.
if (nzchar(Sys.getenv("CI"))) for (p in c("TemporalHazard", "ggplot2")) loadNamespace(p)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  eval = has_pkg
)
```

The previous vignettes covered each piece of the analytical workflow
in isolation: how to fit a model, how to predict from it, how to
validate it. This vignette runs all of those pieces together on a
single dataset, in the order you'd actually run them in a real
clinical analysis. The point isn't to introduce new functions; it's
to show how the pieces chain into a disciplined sequence that turns
raw covariates into a fitted, validated, defensible hazard model.

The sequence mirrors the one in the original SAS HAZARD system, which
codified what cardiothoracic-surgery biostatisticians had been doing
informally for decades:

1. **Nonparametric baseline**: Kaplan-Meier life table to see what
   the data is saying before any model intervenes.
2. **Shape fitting**: start with a single-distribution Weibull,
   build up to a multiphase decomposition when the KM curve
   demands it.
3. **Variable screening**: univariable association tests and
   calibration plots to decide which covariates enter the model and
   in what functional form.
4. **Multivariable model**: covariates layered onto a hazard shape
   you already trust.
5. **Prediction**: patient-specific risk profiles for clinical
   reporting and decision support.
6. **Validation**: decile-of-risk calibration to verify the model
   is honest across the risk spectrum, not just on average.

This corresponds to the SAS programs `ac.*` → `hz.*` → `lg.*` →
`hm.*` → `hp.*` → `hs.*`. If you're familiar with the SAS workflow,
the section structure should feel familiar; if not, treat this as
the canonical analytical sequence to follow on any new dataset.

We use the **AVC** dataset (310 patients, atrioventricular canal
repair) which has rich covariates and two identifiable hazard
phases. That is fewer phases than the 3-phase CABG example you've seen in
other vignettes, but AVC has the covariate complexity needed to
exercise the screening, multivariable-fit, and validation steps.

## Data preparation

Load the package and the AVC dataset, drop incomplete rows so the
design matrix is rectangular for the multivariable fits to come, and
inspect the resulting column types and ranges. The `na.omit()` step
is conservative (losing rows is a real cost), but for a walkthrough
we want every later fit to use the same patient set so the comparisons
are apples-to-apples.

```{r}
#| label: setup
library(TemporalHazard)
if (has_ggplot2) library(ggplot2)

data(avc)
avc <- na.omit(avc)
str(avc)
```

The AVC dataset contains `r if (exists("avc", inherits = TRUE)) nrow(avc) else 310` patients after removing incomplete
cases.  Key covariates include age at repair (`age`, in months), NYHA
functional class (`status`), presence of malalignment (`mal`), and
interventricular communication (`com_iv`).

## Step 1: Nonparametric baseline

Before fitting any parametric model, establish the empirical survival curve
using the Kaplan-Meier estimator.  This is the benchmark against which all
parametric fits will be compared.

`hzr_kaplan()` wraps `survival::survfit()` and adds logit-transformed exact
confidence limits that respect the `[0, 1]` boundary and are more accurate in
the tails than the default Greenwood intervals.  This matches the SAS
`kaplan.sas` macro output structure.

```{r}
#| label: km-life-table
km <- hzr_kaplan(time = avc$int_dead, status = avc$dead)
head(km)
```

The returned data frame is the life table: one row per event time with
`n_risk`, `n_event`, `survival`, logit-transformed 95% CLs
(`cl_lower` / `cl_upper`), and a KM-based cumulative hazard
(`-log(survival)`; use `hzr_nelson()` if you need the Nelson-Aalen
estimator instead).  The plot uses the survival column for the step
curve and the two confidence limits for the band:

```{r}
#| label: fig-km-baseline
#| fig-cap: "Kaplan-Meier survival estimate for AVC patients (n = 305)"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
km_df <- data.frame(
  time     = km$time,
  survival = km$survival * 100,
  lower    = km$cl_lower * 100,
  upper    = km$cl_upper * 100
)

ggplot(km_df, aes(time, survival)) +
  geom_step(colour = "#D55E00", linewidth = 0.8) +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              stat = "identity", alpha = 0.15, fill = "#D55E00") +
  labs(x = "Months after AVC repair", y = "Survival (%)") +
  coord_cartesian(ylim = c(0, 100)) +
  theme_minimal()
```

The KM curve shows a sharp early drop (operative mortality) followed by a
roughly constant attrition rate.  This two-phase pattern suggests an early
CDF phase plus a constant phase, with no obvious late rising hazard.

## Step 2: Shape fitting (simple to complex)

### 2a. Single-phase Weibull

The discipline here is to start with the simplest parametric model
and earn any added complexity. A single-distribution Weibull
intercept-only fit gives us a smooth two-parameter curve to compare
against the KM baseline. If the Weibull tracks the KM step function
reasonably well, we may be done; if it can't, the gap tells us
exactly where a richer model needs to go.

```{r}
#| label: weibull-fit
fit_weib <- hazard(
  survival::Surv(int_dead, dead) ~ 1,
  data  = avc,
  dist  = "weibull",
  theta = c(mu = 0.05, nu = 0.5),
  fit   = TRUE
)
summary(fit_weib)
```

The Weibull forces a monotone hazard shape.  Let's overlay it on the
Kaplan-Meier to see where it fits well and where it doesn't.

```{r}
#| label: fig-weibull-overlay
#| fig-cap: "Single Weibull vs. Kaplan-Meier"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
t_grid <- seq(0.01, max(avc$int_dead) * 0.9, length.out = 200)
surv_weib <- predict(fit_weib,
                     newdata = data.frame(time = t_grid),
                     type = "survival") * 100

ggplot() +
  geom_step(data = km_df, aes(time, survival),
            colour = "#D55E00", linewidth = 0.6) +
  geom_line(data = data.frame(time = t_grid, survival = surv_weib),
            aes(time, survival), colour = "#0072B2", linewidth = 0.8) +
  labs(x = "Months after AVC repair", y = "Survival (%)") +
  coord_cartesian(ylim = c(0, 100)) +
  theme_minimal()
```

The single Weibull typically misses the sharp early drop; it compromises
between the early and late time frames.

### 2b. Two-phase model (early CDF + constant)

The KM curve drops steeply in the first few months (operative
mortality) and then settles into a roughly linear decline. The
Weibull tried to capture both with one shape and ended up
compromising. Splitting the hazard into two phases lets each
mechanism have its own parameterization: a `"cdf"` early phase that
saturates as operative risk resolves, plus a `"constant"` phase that
carries the steady background attrition. Two phases is what AVC
actually needs; CABG with longer follow-up would need a third
late-rising phase to handle graft deterioration, but the AVC
follow-up window doesn't extend far enough to identify one.

```{r}
#| label: multiphase-fit
fit_mp <- hazard(
  survival::Surv(int_dead, dead) ~ 1,
  data   = avc,
  dist   = "multiphase",
  phases = list(
    early    = hzr_phase("cdf", t_half = 0.5, nu = 1, m = 1,
                          fixed = "shapes"),
    constant = hzr_phase("constant")
  ),
  fit     = TRUE,
  control = list(n_starts = 5, maxit = 1000)
)
summary(fit_mp)
```

Note the use of `fixed = "shapes"`. We fix the temporal shape parameters
and only estimate the scale (log_mu) for each phase.  This matches the
standard HAZARD workflow: shapes are set from clinical knowledge or
preliminary exploration, then scales and covariates are estimated.

```{r}
#| label: fig-multiphase-overlay
#| fig-cap: "Two-phase parametric model vs. Kaplan-Meier"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
surv_mp <- predict(fit_mp,
                   newdata = data.frame(time = t_grid),
                   type = "survival") * 100

ggplot() +
  geom_step(data = km_df, aes(time, survival),
            colour = "#D55E00", linewidth = 0.6) +
  geom_line(data = data.frame(time = t_grid, survival = surv_mp),
            aes(time, survival), colour = "#0072B2", linewidth = 0.8) +
  labs(x = "Months after AVC repair", y = "Survival (%)") +
  coord_cartesian(ylim = c(0, 100)) +
  theme_minimal()
```

### 2c. Decomposed hazard

The overlay plot shows the *total* fit tracks KM well, but it
doesn't show whether each phase is doing the job we asked of it.
Decomposing the cumulative hazard into per-phase contributions is
the diagnostic that answers that question: the early phase should
account for most of the steep early drop, the constant phase should
carry the rest, and neither phase should be doing implausible work
in the wrong time window.

```{r}
#| label: fig-decomposed
#| fig-cap: "Phase decomposition of cumulative hazard"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
decomp <- predict(fit_mp,
                  newdata = data.frame(time = t_grid),
                  type = "cumulative_hazard",
                  decompose = TRUE)

decomp_df <- data.frame(
  time  = t_grid,
  total = decomp[, "total"],
  early = decomp[, "early"],
  const = decomp[, "constant"]
)

ggplot(decomp_df, aes(x = time)) +
  geom_line(aes(y = total), linewidth = 0.9, colour = "black") +
  geom_line(aes(y = early), linewidth = 0.7, colour = "#E69F00",
            linetype = "dashed") +
  geom_line(aes(y = const), linewidth = 0.7, colour = "#56B4E9",
            linetype = "dashed") +
  labs(x = "Months after AVC repair", y = "Cumulative hazard") +
  theme_minimal()
```

The early phase contributes most of its risk in the first few months, then
levels off.  The constant phase accumulates linearly over time.

## Step 3: Variable screening

Before entering covariates into the hazard model, screen them for
association with the outcome.  This step corresponds to the `lg.*` programs
in the SAS workflow.

### 3a. Univariable logistic screening

The cheapest screening tool we have is a univariable logistic
regression of the event indicator on each covariate. It throws away
the time-to-event structure but it answers a binary question
quickly: does this covariate have *any* association with mortality
at all? Covariates whose univariable p-value is huge can be
deprioritized, but not discarded outright, because adjusting for
other covariates can reveal an association the univariable test
misses. Covariates with small
univariable p-values deserve closer functional-form inspection
before they enter the formula.

```{r}
#| label: screening
covariates <- c("age", "status", "mal", "com_iv", "orifice",
                "inc_surg", "opmos")

screen_results <- data.frame(
  variable = covariates,
  coef     = numeric(length(covariates)),
  p_value  = numeric(length(covariates))
)

for (i in seq_along(covariates)) {
  v <- covariates[i]
  fml <- as.formula(paste("dead ~", v))
  lg <- glm(fml, data = avc, family = binomial)
  s <- summary(lg)$coefficients
  if (nrow(s) >= 2) {
    screen_results$coef[i] <- s[2, 1]
    screen_results$p_value[i] <- s[2, 4]
  }
}

screen_results <- screen_results[order(screen_results$p_value), ]
screen_results$significant <- ifelse(screen_results$p_value < 0.05,
                                      "*", "")
screen_results
```

### 3b. Functional form assessment

For continuous covariates that appear significant, examine whether the
relationship with outcome is monotone on the logit scale.  `hzr_calibrate()`
implements the SAS `logit.sas` / `logitgr.sas` macros: bin the continuous
covariate into quantile groups, compute the observed event rate per bin,
and transform it to the logit scale.

```{r}
#| label: age-calibrate
cal_age <- hzr_calibrate(x = avc$age, event = avc$dead,
                          groups = 10, link = "logit")
cal_age
```

Plot the logit-transformed event rate against the bin centre (`mean`
column).  A straight line means the covariate can enter the model
untransformed; monotone-but-curved shapes suggest a log or square-root
transform; U-shapes call for a quadratic term.

```{r}
#| label: fig-calibrate-age
#| fig-cap: "Decile calibration: age vs. logit(P(death))"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
ggplot(cal_age, aes(mean, link_value)) +
  geom_point(size = 3, colour = "#0072B2") +
  geom_line(colour = "#0072B2", alpha = 0.4) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE,
              colour = "#D55E00", linetype = "dashed") +
  labs(x = "Age at repair (months, decile midpoint)",
       y = "logit(P(death))") +
  theme_minimal()
```

The blue points are the observed decile logits; the dashed orange line is a
linear reference.  Deviations from linearity flag where a transform is
warranted.

## Step 4: Multivariable model

### 4a. Manual specification

We've already established the two-phase shape and screened the
covariates; the multivariable model is the synthesis. Drop the
covariates that passed screening onto the right-hand side of the
formula and refit. The shape parameters stay fixed (we've earned
that decision from §2); the optimizer estimates the two phase
scales and one coefficient per covariate, jointly. This is the
working model you'd report from this analysis.

```{r}
#| label: multivariable-fit
fit_mv <- hazard(
  survival::Surv(int_dead, dead) ~ age + status + mal + com_iv,
  data   = avc,
  dist   = "multiphase",
  phases = list(
    early    = hzr_phase("cdf", t_half = 0.5, nu = 1, m = 1,
                          fixed = "shapes"),
    constant = hzr_phase("constant")
  ),
  fit     = TRUE,
  control = list(n_starts = 5, maxit = 1000)
)
summary(fit_mv)
```

The coefficient table shows phase-specific covariate effects.  A positive
coefficient means higher hazard (worse prognosis).  Interpretation:
covariates scale the phase-specific hazard multiplicatively via
exp(x * beta).

### 4b. Automated stepwise selection

The manual approach works when screening has already narrowed the
candidate pool, but with a larger pool it helps to let the model
choose.  `hzr_stepwise()` runs forward, backward, or two-way
selection against an existing `hazard` fit, scoring candidates with
the score statistic, Wald p-values or delta-AIC.  Defaults match SAS `PROC HAZARD`
(`SLENTRY = 0.30`, `SLSTAY = 0.20`).

We start from a baseline with no covariates, offer the same screened
candidate pool, and let the algorithm decide.  For multiphase models
the scope is a named list of one-sided formulas keyed by phase, so
each phase can advertise its own candidate set; stepping operates
per `(variable, phase)` pair, letting a covariate enter one phase
and not another.  Single-distribution models accept either a flat
one-sided formula or a character vector of names.

`hzr_stepwise()` defaults to `criterion = "score"`, which reproduces
SAS's `SELECTION` Q statistic by testing each candidate at the current
estimates without a per-candidate refit.  Here we pass `criterion = "wald"`
deliberately, so the `$steps` table below reports the refit-based Wald
p-values.

```{r}
#| label: stepwise-forward
#| eval: false
base_mp <- hazard(
  survival::Surv(int_dead, dead) ~ 1,
  data   = avc,
  dist   = "multiphase",
  phases = list(
    early    = hzr_phase("cdf", t_half = 0.5, nu = 1, m = 1,
                          fixed = "shapes"),
    constant = hzr_phase("constant")
  ),
  fit     = TRUE,
  control = list(n_starts = 3, maxit = 500)
)

fit_step <- hzr_stepwise(
  base_mp,
  scope     = list(
    early    = ~ age + status + mal + com_iv,
    constant = ~ age + status + mal + com_iv
  ),
  data      = avc,
  direction = "both",
  criterion = "wald",
  trace     = FALSE,
  control   = list(n_starts = 2, maxit = 500)
)
```

```{r}
#| label: stepwise-forward-load
#| include: false
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE)
# Refitting the chunk above takes long enough to push R CMD check past
# CRAN's time limit, so data-raw/vignette-fit-step.R stores its result and
# test-vignette-fit-step.R refits and compares. Without the file, refit.
stored <- system.file("extdata", "vignette-fit-step.rds",
                      package = "TemporalHazard")
if (nzchar(stored)) {
  fit_step <- readRDS(stored)
} else {
  eval(parse(text = knitr::knit_code$get("stepwise-forward")))
}
```

```{r}
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE)
fit_step
```

The `$steps` table records every accepted action with its Wald
statistic and p-value:

```{r}
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE)
fit_step$steps[, c("step_num", "action", "variable", "phase",
                   "p_value", "aic")]
```

The final model is the fit at the top of the object, reachable
through the usual hazard accessors.  AIC counts only the estimated
parameters.  Both models hold the early phase's three shape parameters
fixed (`fixed = "shapes"`), and `fit$fit$fixed_mask` flags them, so we
leave them out of the count, as the `Final model` line of `fit_step`
does:

```{r}
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE)
logLik_manual <- fit_mv$fit$objective
logLik_step   <- fit_step$fit$objective
n_free <- function(fit) sum(!fit$fit$fixed_mask)
c(manual = logLik_manual, stepwise = logLik_step,
  free_manual = n_free(fit_mv), free_step = n_free(fit_step),
  aic_manual = 2 * n_free(fit_mv) - 2 * logLik_manual,
  aic_step   = 2 * n_free(fit_step) - 2 * logLik_step)
```

The two models need not agree.  Here stepwise enters `status` in both
phases but `age`, `mal` and `com_iv` in the early phase only, so it
carries fewer free parameters than the manual fit (`free_step` against
`free_manual` above).  Its log-likelihood is lower and its AIC is lower
too (`aic_step` against `aic_manual`), because the dropped constant-phase
coefficients cost more in parameters than they bought in fit.  Even on
an identical covariate set the two log-likelihoods can differ, since
the fits use different `n_starts` and the multiphase likelihood can
have more than one optimum.  For an AIC-driven run use `criterion = "aic"`; for a
forward-only sweep `direction = "forward"`.  See `?hzr_stepwise` for
the full argument surface including `force_in`, `force_out`, and the
`max_move` oscillation guard.

## Step 5: Prediction

### 5a. Baseline survival curve

Generate the survival curve for a reference patient (median age, NYHA
class 1, no malalignment, no interventricular communication).

`predict(..., se.fit = TRUE)` adds a delta-method standard error and 95%
confidence band around every prediction (log(-log) scale for survival so
the band is guaranteed to lie in `[0, 1]`).  This is the parametric
analogue of the KM confidence limits from Step 1, a patient-specific
uncertainty estimate the nonparametric curve cannot provide.

```{r}
#| label: fig-prediction
#| fig-cap: "Reference patient survival with 95% delta-method CI and KM overlay"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
ref_patient <- data.frame(
  time   = t_grid,
  age    = median(avc$age),
  status = 1,
  mal    = 0,
  com_iv = 0
)

surv_ref <- predict(fit_mv, newdata = ref_patient,
                    type = "survival", se.fit = TRUE)

ref_df <- data.frame(
  time     = t_grid,
  survival = surv_ref$fit * 100,
  lower    = surv_ref$lower * 100,
  upper    = surv_ref$upper * 100
)

ggplot() +
  geom_step(data = km_df, aes(time, survival),
            colour = "#D55E00", linewidth = 0.6, alpha = 0.7) +
  geom_ribbon(data = ref_df, aes(time, ymin = lower, ymax = upper),
              fill = "#0072B2", alpha = 0.15) +
  geom_line(data = ref_df, aes(time, survival),
            colour = "#0072B2", linewidth = 0.8) +
  labs(x = "Months after AVC repair", y = "Survival (%)") +
  coord_cartesian(ylim = c(0, 100)) +
  theme_minimal()
```

### 5b. Sensitivity analysis (risk factor comparison)

Coefficient tables tell you about effect *direction and magnitude*
on the log-hazard scale. They don't directly tell a clinician what
to expect at the bedside. Sensitivity analysis closes that gap by
scoring two clinically meaningful profiles (a low-risk patient
drawn from the favorable end of each covariate, and a high-risk
patient drawn from the unfavorable end) through the fitted model
and plotting the resulting survival curves with delta-method
confidence bands. The vertical and horizontal gaps between the two
curves are the model's clinical-impact statement, in the units
(survival probability and time) that the audience actually
recognizes.

```{r}
#| label: fig-sensitivity
#| fig-cap: "Survival by risk profile: low-risk (blue) vs. high-risk (orange)"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
low_risk <- data.frame(
  time   = t_grid,
  age    = quantile(avc$age, 0.25),
  status = 1,
  mal    = 0,
  com_iv = 0
)

high_risk <- data.frame(
  time   = t_grid,
  age    = quantile(avc$age, 0.90),
  status = 4,
  mal    = 1,
  com_iv = 1
)

surv_lo <- predict(fit_mv, newdata = low_risk,
                   type = "survival", se.fit = TRUE)
surv_hi <- predict(fit_mv, newdata = high_risk,
                   type = "survival", se.fit = TRUE)

sens_df <- data.frame(
  time     = rep(t_grid, 2),
  survival = c(surv_lo$fit, surv_hi$fit) * 100,
  lower    = c(surv_lo$lower, surv_hi$lower) * 100,
  upper    = c(surv_lo$upper, surv_hi$upper) * 100,
  profile  = rep(c("Low risk", "High risk"), each = length(t_grid))
)

ggplot(sens_df, aes(time, survival, colour = profile, fill = profile)) +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              alpha = 0.15, colour = NA) +
  geom_line(linewidth = 0.8) +
  scale_colour_manual(values = c("Low risk" = "#0072B2",
                                  "High risk" = "#D55E00")) +
  scale_fill_manual(values = c("Low risk" = "#0072B2",
                                "High risk" = "#D55E00")) +
  labs(x = "Months after AVC repair", y = "Survival (%)",
       colour = NULL, fill = NULL) +
  coord_cartesian(ylim = c(0, 100)) +
  theme_minimal() +
  theme(legend.position = "bottom")
```

Non-overlapping confidence bands between the two profiles indicate the
risk-factor combination is well-identified; overlap around the bands'
centre-of-mass flags where the model can't distinguish the profiles with
the available sample size.

## Step 6: Validation (decile-of-risk calibration)

Partition patients into deciles by predicted risk and compare observed
vs. expected event rates.  Good calibration means the two track each
other.  This implements the workflow of the SAS `deciles.hazard.sas`
macro.

```{r}
#| label: fig-deciles
#| fig-cap: "Observed vs. expected event rates by risk decile"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
cal <- hzr_deciles(fit_mv, time = max(avc$int_dead))
print(cal)
```

```{r}
#| label: fig-decile-plot
#| fig-cap: "Decile calibration: observed (bars) vs. expected (points) event rates"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
ggplot(cal, aes(x = group)) +
  geom_col(aes(y = observed_rate), fill = "#56B4E9", alpha = 0.7) +
  geom_point(aes(y = expected_rate), colour = "#D55E00", size = 3) +
  geom_line(aes(y = expected_rate), colour = "#D55E00", linewidth = 0.5) +
  scale_x_continuous(breaks = seq_len(nrow(cal))) +
  labs(x = "Risk decile (1 = lowest)", y = "Event rate",
       caption = paste("Overall chi-sq =",
                        round(attr(cal, "overall")$chi_sq, 2),
                        ", p =",
                        format.pval(attr(cal, "overall")$p_value,
                                     digits = 3))) +
  theme_minimal()
```

A non-significant overall chi-square statistic (p > 0.05) means the test
found no evidence of poor calibration. It does not prove the model is
calibrated: with ten groups and few events in each, the test has little
power. Read the plot as well. Calibration is adequate when the observed and
expected rates track each other across the deciles.

### Conservation of events check

The conservation-of-events principle says a model fit by maximum
likelihood predicts as many events as were observed: add up every
patient's cumulative hazard at the end of their follow-up and you get
the event count back.  `hzr_gof()` is the R version of the SAS
`hazplot` display, a running tally of observed and expected events
over time.  Each patient's expected count is their own cumulative
hazard, from their own covariates, at the time they leave follow-up.

```{r}
#| label: gof-summary
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE)
gof <- hzr_gof(fit_mv)
print(gof)
```

Expected 68.0 against 68 observed, an E/O of 1, so the fit conserves
events.  `fit_mv` is a multiphase fit with Conservation of Events
applied (the default), which holds the total to the event count
exactly.  A Weibull fit gets there at its maximum, to within the
optimizer's tolerance.

```{r}
#| label: fig-gof-survival
#| fig-cap: "Parametric survival at the covariate means (blue) vs. Kaplan-Meier (orange)"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
ggplot(gof, aes(x = time)) +
  geom_step(aes(y = km_surv * 100), colour = "#D55E00", linewidth = 0.6) +
  geom_line(aes(y = par_surv * 100), colour = "#0072B2", linewidth = 0.8) +
  labs(x = "Months after AVC repair", y = "Survival (%)") +
  coord_cartesian(ylim = c(0, 100)) +
  theme_minimal()
```

```{r}
#| label: fig-gof-events
#| fig-cap: "Cumulative observed (orange) vs. expected (blue) events"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
ggplot(gof, aes(x = time)) +
  geom_line(aes(y = cum_observed), colour = "#D55E00", linewidth = 0.8) +
  geom_line(aes(y = cum_expected), colour = "#0072B2",
            linewidth = 0.8, linetype = "dashed") +
  labs(x = "Months after AVC repair", y = "Cumulative events") +
  theme_minimal()
```

```{r}
#| label: fig-gof-residual
#| fig-cap: "Conservation-of-events residual (expected minus observed) over time"
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
ggplot(gof, aes(x = time, y = residual)) +
  geom_line(linewidth = 0.7, colour = "grey30") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey60") +
  labs(x = "Months after AVC repair",
       y = "Residual (expected - observed)") +
  theme_minimal()
```

The first plot is the one place the mean patient appears.  It sets the
survival curve for a patient with average covariates against the
Kaplan-Meier estimate for the whole cohort, and that curve runs above
it from the first weeks on.  The mean patient fares better than the
cohort, because the mean patient's cumulative hazard is not the average
of the patients' cumulative hazards.  That is why `hzr_gof()` builds
its expected count patient by patient rather than from this curve.

The other two plots tally events as patients leave follow-up.  A
patient adds their expected count when they leave, by death or
censoring, and that count is the hazard they built up over their whole
follow-up.  Most deaths come in the first weeks, to patients who had
little time to build up hazard: 47 deaths in the first half month
against about 4.8 expected.  The residual bottoms out near -47.5 at
about two months.  The patients followed longest, most of them
censored, carry the rest of the expected count, and it arrives only as
they leave follow-up years later, so the residual climbs back and ends
at zero.  Read the shape of the curve as a record of when deaths and
follow-up happen in this cohort, not as a month-by-month calibration
check.  The end point is the conservation-of-events identity, and the
decile table above is the calibration check.

## Why not Cox regression?

For this kind of analysis the temporal parametric approach buys you
several things that Cox proportional hazards does not:

| Feature | Cox (`coxph`) | TemporalHazard |
|:---|:---:|:---:|
| Proportional hazards required | Yes | No |
| Baseline hazard specified | No | Yes (parametric) |
| Multiple hazard phases | No | Yes (additive) |
| Patient-specific survival curve | Step function (Breslow baseline) | Smooth, parametric |
| Smooth extrapolation beyond data | No | Yes |
| Conservation of events | Yes (martingale residuals sum to 0) | Yes (per-subject cumulative hazards sum to the event count) |
| Interval censoring | Not standard | Supported |

Many clinical outcomes show **non-proportional
hazards**: the risk factors that dominate early mortality (e.g., surgical
complexity) differ from those driving late attrition (e.g., age, comorbidity).
The multiphase model captures this through phase-specific
covariate effects.

## Summary

The complete analytical workflow follows a disciplined sequence:

1. **Kaplan-Meier baseline**: establish the empirical survival pattern
2. **Shape fitting**: match the temporal shape, simple to complex
3. **Variable screening**: logistic screening, decile-logit functional form
4. **Multivariable model**: enter covariates with fixed shapes
5. **Prediction**: reference curves, sensitivity analysis
6. **Validation**: decile calibration, conservation-of-events check

This sequence ensures that the temporal shape is established before
covariates are introduced, and that the final model is validated against
the observed data.  See `vignette("getting-started")` for the minimal
API workflow and `vignette("mf-mathematical-foundations")` for the
mathematical details.