---
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
has_pkg <- requireNamespace("TemporalHazard", quietly = TRUE)
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
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 — fewer phases than the 3-phase CABG example you've seen in
other vignettes, but with 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 --- 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).  Plotting just requires the survival column:

```{r}
#| label: fig-km-baseline
#| fig-cap: "Kaplan-Meier survival estimate for AVC patients (n = 310)"
#| 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 --- 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: an `"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 will not
suddenly become significant in the multivariable hazard model
either — those can be deprioritized. 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 red 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
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.

```{r}
#| label: stepwise-forward
#| eval: !expr requireNamespace("TemporalHazard", quietly = TRUE)
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)
)
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:

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

When the screening and stepwise agree on the same covariate set the
log-likelihoods match exactly; when stepwise selects a leaner model
AIC drops.  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 (red)"
#| 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) indicates
adequate calibration --- the model's predictions are consistent with
observed event rates across the risk spectrum.

### Conservation of events check

The conservation-of-events principle states that a well-fitting model
should predict the same total number of events as actually observed.
`hzr_gof()` tracks this cumulatively over time --- the residual
(expected minus observed) should stay near zero.

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

```{r}
#| label: fig-gof-survival
#| fig-cap: "Parametric (blue) vs. Kaplan-Meier (orange) survival"
#| 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()
```

A residual that hovers near zero across the full time range indicates the
model conserves events well.  Persistent positive residual means the model
overpredicts risk; persistent negative means it underpredicts.

## 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 prediction | Approximate | Exact |
| Smooth extrapolation beyond data | No | Yes |
| Conservation of events check | No | Yes (`hzr_gof()`) |
| 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, LOESS 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.