---
title: "Prediction & Visualization"
subtitle: "Survival curves, phase decomposition, and risk profiles"
format:
  html:
    code-fold: false
    toc: true
    toc-depth: 3
    number-sections: true
vignette: >
  %\VignetteIndexEntry{Prediction & Visualization}
  %\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) &&
  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  = "#>",
  eval     = has_pkg
)
```

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

A fitted hazard model is only useful insofar as you can extract
answers from it. This vignette covers the predict-and-visualize half
of the workflow: pulling survival probabilities, cumulative hazards,
and patient-level risk scores out of a fitted model, attaching
delta-method uncertainty to those predictions, and plotting them in
the formats that come up in actual clinical analyses.

The structure mirrors how a real analysis unfolds. We start with the
Kaplan-Meier baseline (the nonparametric reference every parametric
fit gets judged against), enumerate the prediction types the model
exposes, plot the parametric survival curve against KM as a
goodness-of-fit check, layer on confidence bands, then move to the
multiphase decomposition and patient-specific risk profiles where the
covariate machinery pays off. The last section reuses the workflow
across multiple endpoints on the same cohort. If you haven't seen the
fitting workflow yet, `vignette("fitting-hazard-models")` is the
prerequisite.


## Kaplan-Meier baseline

The Kaplan-Meier curve is what the data itself says about survival
before any parametric model intervenes. It makes no assumption about
the shape of the hazard, handles censoring honestly, and gives you a
nonparametric step function that's the gold-standard reference for any
follow-up fit. Every parametric prediction that follows gets compared
back to this curve. If a model's smooth prediction can't track the KM
step function, the model is missing something the data has been
telling you all along.

```{r}
#| label: km-baseline
data(avc)
avc <- na.omit(avc)
km <- survfit(Surv(int_dead, dead) ~ 1, data = avc)
```

```{r}
#| label: fig-km
#| fig-cap: "Kaplan-Meier survival estimate: death after AVC repair"
#| fig-width: 7
#| fig-height: 4
km_df <- data.frame(time = km$time, survival = km$surv * 100)

ggplot(km_df, aes(time, survival)) +
  geom_step(linewidth = 0.6) +
  scale_y_continuous(limits = c(0, 100)) +
  labs(x = "Months after repair", y = "Freedom from death (%)") +
  theme_minimal()
```


## Prediction types

The `predict()` method exposes four quantities through its `type=`
argument, each useful for a different downstream question:

- `"linear_predictor"`: $x^\top \beta$, the covariate-driven
  log-hazard shift. Use it to rank patients on relative risk without
  pinning down absolute survival probabilities.
- `"hazard"`: $\exp(x^\top \beta)$, the multiplicative hazard ratio
  relative to the baseline. Use it when you want hazard ratios for
  reporting or for further calculation.
- `"cumulative_hazard"`: $H(t \mid x)$, the integrated hazard up to
  time $t$ for each row. Use it when you need the raw integrated
  intensity (decomposing it by phase, for example, or computing
  derived quantities like the expected number of events).
- `"survival"`: $S(t \mid x) = \exp(-H(t \mid x))$, the probability
  of surviving past time $t$. This is what clinicians and patients
  actually want to see; it's also what we plot for diagnostics.

The two predictions we extract below are `"survival"` and
`"cumulative_hazard"`: survival for the clinical communication, and
the cumulative hazard so we can sanity-check the relationship
$S = \exp(-H)$ holds row-by-row. We fit a multivariable Weibull on the
AVC death endpoint first.

```{r}
#| label: fit-mv
fit <- hazard(
  Surv(int_dead, dead) ~ age + status + mal + com_iv,
  data  = avc,
  dist  = "weibull",
  theta = c(mu = 0.20, nu = 1.0, rep(0, 4)),
  fit   = TRUE,
  control = list(maxit = 500)
)
```

Pick a representative patient (here a median-aged, mid-status,
no-malalignment, no-interventricular-communication profile) and evaluate the
prediction types over a fine time grid. The result is a data frame
with one row per time point, ready for plotting or downstream
analysis.

```{r}
#| label: pred-types
t_grid <- seq(0.01, max(avc$int_dead) * 0.95, length.out = 200)
profile <- data.frame(
  time   = t_grid,
  age    = median(avc$age),
  status = 2,
  mal    = 0,
  com_iv = 0
)

surv   <- predict(fit, newdata = profile, type = "survival")
cumhaz <- predict(fit, newdata = profile, type = "cumulative_hazard")

profile$survival          <- surv
profile$cumulative_hazard <- cumhaz

head(profile[, c("time", "survival", "cumulative_hazard")])
```


## Parametric survival with KM overlay

The fundamental diagnostic for any parametric survival fit is whether
its predictions track the Kaplan-Meier step function. Plot the
median-profile survival curve from the previous chunk on top of the
KM estimate from the raw cohort. Where the parametric curve hugs the
steps, the model is faithful to the data; where it drifts, it's
imposing a shape the data doesn't support. This isn't a hypothesis
test; it's an eyeball check, and it's the single most informative
thing you can do with a fitted model before trusting its
predictions.

```{r}
#| label: fig-surv-overlay
#| fig-cap: "Weibull parametric survival vs. Kaplan-Meier (AVC death)"
#| fig-width: 7
#| fig-height: 4.5
ggplot() +
  geom_step(data = km_df, aes(time, survival, colour = "Kaplan-Meier"),
            linewidth = 0.5) +
  geom_line(data = profile,
            aes(time, survival * 100, colour = "Parametric (Weibull)"),
            linewidth = 1) +
  scale_colour_manual(
    values = c("Parametric (Weibull)" = "#0072B2",
               "Kaplan-Meier"         = "#D55E00")
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  labs(x = "Months after repair", y = "Freedom from death (%)",
       colour = NULL) +
  theme_minimal() +
  theme(legend.position = "bottom")
```

## Confidence limits on predictions

A point estimate without an uncertainty band tells you what the model
thinks, not how confident it is. For patient-level survival
predictions, a survival probability of 0.85
with a tight band around it is a different number than
a survival probability of 0.85 with a 0.4–0.97 band around it, even
though both round to "85%". Delta-method confidence limits let you
plot the second case honestly instead of pretending it's the first.

`predict.hazard()` accepts `se.fit = TRUE` (default
`FALSE`) and `level = 0.95` to return delta-method standard errors
and confidence limits alongside the point estimate. The return value
changes shape: a plain numeric vector with `se.fit = FALSE`, a data
frame with columns `fit`, `se.fit`, `lower`, and `upper` with
`se.fit = TRUE`.

```{r}
#| label: pred-with-se
# Build a clean newdata frame (the earlier chunk appended result
# columns to `profile`, which would confuse predict()'s column count).
profile_ci <- data.frame(
  time   = t_grid,
  age    = median(avc$age),
  status = 2,
  mal    = 0,
  com_iv = 0
)
surv_ci <- predict(fit, newdata = profile_ci,
                   type = "survival", se.fit = TRUE, level = 0.95)
head(surv_ci)
```

Confidence limits use SAS-matched transformations chosen so the
interval respects the natural range of each quantity: log-scale for
`hazard` and `cumulative_hazard` (keeps the lower bound positive),
and `log(-log(S))` for `survival` (keeps the interval inside
$[0, 1]$). The linear predictor uses symmetric natural-scale CLs.
A symmetric CI on
survival probability would frequently dip below 0 or rise above 1,
giving nonsense values that the transformed form avoids by construction.

```{r}
#| label: fig-surv-ci
#| fig-cap: "Parametric survival with 95% delta-method confidence band"
#| fig-width: 7
#| fig-height: 4.5
ci_df <- data.frame(
  time     = profile_ci$time,
  survival = surv_ci$fit * 100,
  lower    = surv_ci$lower * 100,
  upper    = surv_ci$upper * 100
)

ggplot() +
  geom_step(data = km_df, aes(time, survival, colour = "Kaplan-Meier"),
            linewidth = 0.5) +
  geom_ribbon(data = ci_df,
              aes(time, ymin = lower, ymax = upper),
              fill = "#0072B2", alpha = 0.15) +
  geom_line(data = ci_df,
            aes(time, survival, colour = "Parametric + 95% CI"),
            linewidth = 1) +
  scale_colour_manual(
    values = c("Parametric + 95% CI" = "#0072B2",
               "Kaplan-Meier"        = "#D55E00")
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  labs(x = "Months after repair", y = "Freedom from death (%)",
       colour = NULL) +
  theme_minimal() +
  theme(legend.position = "bottom")
```

Weibull and multiphase use a
closed-form Jacobian for the delta-method calculation. Exponential,
log-logistic, and log-normal fall back to `numDeriv::jacobian()` on a
per-call cumulative-hazard closure. The results are identical; the
fallback is only slightly slower because the Jacobian is computed
numerically per prediction point. The user-facing API is the same
regardless of distribution.

## Decomposed multiphase hazard

The multiphase model's whole point is that the total hazard is a sum
of phase contributions. `predict()` with `decompose = TRUE` and
`type = "cumulative_hazard"` exposes those contributions row by row,
returning a data frame with one column per phase plus a `total`
column. Numerically differentiating each column gives you the
instantaneous hazard rate for each phase, which is the diagnostic
that tells you whether the multiphase model is doing what you asked
of it. The early phase should dominate near $t = 0$ and fall off, the
constant phase should be a flat floor, and the late phase should sit
near zero early and rise after a lag.

```{r}
#| label: fit-mp
data(cabgkul)

fit_mp <- hazard(
  Surv(int_dead, dead) ~ 1,
  data   = cabgkul,
  dist   = "multiphase",
  phases = list(
    early    = hzr_phase("cdf", t_half = 0.2, nu = 1, m = 1,
                          fixed = "shapes"),
    constant = hzr_phase("constant"),
    late     = hzr_phase("g3",  tau = 1, gamma = 3, alpha = 1, eta = 1,
                          fixed = "shapes")
  ),
  fit     = TRUE,
  control = list(n_starts = 5, maxit = 1000)
)
```

```{r}
#| label: fig-decomposed-hazard
#| fig-cap: "Additive phase decomposition: total hazard rate (solid) = early + constant + late (dashed)"
#| fig-width: 7
#| fig-height: 4.5
t_mp <- seq(0.01, max(cabgkul$int_dead) * 0.95, length.out = 200)
nd   <- data.frame(time = t_mp)

decomp <- predict(fit_mp, newdata = nd, type = "cumulative_hazard",
                  decompose = TRUE)

# Numerical differentiation: h(t) ≈ ΔH(t) / Δt
num_hazard <- function(cumhaz, time) {
  dt <- diff(time)
  dH <- diff(cumhaz)
  c(dH[1] / dt[1], dH / dt)
}

h_long <- rbind(
  data.frame(time = t_mp, hazard = num_hazard(decomp$early, t_mp),
             Phase = "Early"),
  data.frame(time = t_mp, hazard = num_hazard(decomp$constant, t_mp),
             Phase = "Constant"),
  data.frame(time = t_mp, hazard = num_hazard(decomp$late, t_mp),
             Phase = "Late"),
  data.frame(time = t_mp, hazard = num_hazard(decomp$total, t_mp),
             Phase = "Total")
)
h_long$Phase <- factor(h_long$Phase,
                       levels = c("Total", "Early", "Constant", "Late"))

ggplot(h_long, aes(time, hazard, colour = Phase, linetype = Phase)) +
  geom_line(aes(linewidth = Phase)) +
  scale_colour_manual(values = c(Total = "#222222", Early = "#E69F00",
                                 Constant = "#56B4E9", Late = "#CC79A7")) +
  scale_linetype_manual(values = c(Total = "solid", Early = "dashed",
                                   Constant = "dashed", Late = "dashed")) +
  scale_linewidth_manual(values = c(Total = 1.3, Early = 0.7,
                                    Constant = 0.7, Late = 0.7)) +
  labs(x = "Months after CABG", y = "Hazard rate",
       colour = "Phase", linetype = "Phase", linewidth = "Phase") +
  theme_minimal() +
  theme(legend.position = "bottom")
```

Each phase is doing its job. The early (orange) phase captures the
steep post-operative risk that peaks within the first months. The
constant (blue) phase represents the ongoing background mortality
that persists once patients are past the operative window. The late
(pink) phase captures the gradually increasing risk of late
attrition (graft failure, comorbidity progression, the aging cohort).
The early and late shapes were set by hand: `fixed = "shapes"` holds
`t_half`, `nu` and `m` for the early phase and `tau`, `gamma`,
`alpha` and `eta` for the late phase at the values in the
`hzr_phase()` calls. The optimizer estimated only the three scales,
one per phase, which set how much each phase contributes. If a phase looked wrong here (a constant phase dominating
where we expected an early peak, or a late phase that never rose)
that would be a signal to revisit either the starting values, the
shape parameters, or the data itself.


## Multiphase survival with KM overlay

The decomposed-hazard plot above tells us the model has the right
*shape* per phase. To check that the *total* fit also tracks the
data, we collapse back to the overall survival curve and overlay it
on the Kaplan-Meier estimate for the same CABGKUL cohort. This is the
same diagnostic as the single-Weibull case earlier in this vignette,
but on the CABG data the multiphase model was fit to, and against a
model that has the structural flexibility to follow the cohort's
actual mortality pattern.

```{r}
#| label: fig-mp-surv
#| fig-cap: "Multiphase parametric survival vs. Kaplan-Meier (CABGKUL death)"
#| fig-width: 7
#| fig-height: 4.5
km_cabg    <- survfit(Surv(int_dead, dead) ~ 1, data = cabgkul)
km_cabg_df <- data.frame(time = km_cabg$time, survival = km_cabg$surv * 100)
surv_mp    <- predict(fit_mp, newdata = nd, type = "survival") * 100

ggplot() +
  geom_step(data = km_cabg_df,
            aes(time, survival, colour = "Kaplan-Meier"),
            linewidth = 0.5) +
  geom_line(data = data.frame(time = t_mp, survival = surv_mp),
            aes(time, survival, colour = "Multiphase (3-phase)"),
            linewidth = 1) +
  scale_colour_manual(
    values = c("Multiphase (3-phase)" = "#0072B2",
               "Kaplan-Meier"         = "#D55E00")
  ) +
  scale_y_continuous(limits = c(0, 100)) +
  labs(x = "Months after CABG", y = "Freedom from death (%)",
       colour = NULL) +
  theme_minimal() +
  theme(legend.position = "bottom")
```


## Patient-specific risk profiles

This is the payoff of having covariates in the model. A
population-level curve is fine for cohort-level reporting, but the
question a clinician usually wants answered is: *what does the model
predict for this particular patient, given their risk factors?*
Holding the fitted model fixed and varying the covariate profile
generates patient-specific survival curves you can show side by side.
Below we score three plausible profiles drawn from the cohort's own
quantiles: a low-risk patient (young, mild status, no malalignment,
no interventricular communication), a median patient, and a high-risk
patient (older, severe status, with both malalignment and
interventricular communication).

```{r}
#| label: fig-risk-profiles
#| fig-cap: "Predicted survival by risk profile"
#| fig-width: 7
#| fig-height: 4.5
profiles <- list(
  "Low risk"  = data.frame(age = quantile(avc$age, 0.25),
                            status = 1, mal = 0, com_iv = 0),
  "Median"    = data.frame(age = median(avc$age),
                            status = 2, mal = 0, com_iv = 0),
  "High risk" = data.frame(age = quantile(avc$age, 0.90),
                            status = 4, mal = 1, com_iv = 1)
)

curves <- do.call(rbind, lapply(names(profiles), function(nm) {
  nd <- profiles[[nm]][rep(1, length(t_grid)), ]
  nd$time <- t_grid
  data.frame(time = t_grid,
             survival = predict(fit, newdata = nd, type = "survival") * 100,
             Profile = nm)
}))
curves$Profile <- factor(curves$Profile,
                         levels = c("Low risk", "Median", "High risk"))

ggplot(curves, aes(time, survival, colour = Profile)) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = c("Low risk" = "#009E73",
                                 "Median"   = "#0072B2",
                                 "High risk" = "#D55E00")) +
  scale_y_continuous(limits = c(0, 100)) +
  labs(x = "Months after AVC repair", y = "Freedom from death (%)",
       colour = NULL) +
  theme_minimal() +
  theme(legend.position = "bottom")
```

The gap between the three curves is the model's *prognostic
discrimination*: a wider spread means the covariates are doing real
work, separating patients with different actual risk levels. A
collapsed plot (all three curves on top of each other) means the
covariate effects are too weak to distinguish patients meaningfully,
even if they're statistically significant.


## Multi-endpoint visualization: valves

When a cohort has multiple clinical endpoints (death, valve
endocarditis, reoperation), putting them on the same survival axis
gives an immediate comparative picture of which event the cohort is
most at risk for, and over what time horizon. The `valves` dataset
has separate event-time pairs for each endpoint, so any of them can
go on shared axes as a Kaplan-Meier curve without needing a
parametric fit at all. Below we plot two: death and prosthetic valve
endocarditis (PVE). (You'd parametrize them in the next step, one
endpoint per `hazard()` call, as we did in
`vignette("fitting-hazard-models")`.)

```{r}
#| label: fig-valves-endpoints
#| fig-cap: "Freedom from death and PVE after valve replacement"
#| fig-width: 7
#| fig-height: 4.5
data(valves)
valves <- na.omit(valves)

km_death <- survfit(Surv(int_dead, dead) ~ 1, data = valves)
km_pve   <- survfit(Surv(int_pve, pve) ~ 1, data = valves)

ep_df <- rbind(
  data.frame(time = km_death$time, survival = km_death$surv * 100,
             Endpoint = "Death"),
  data.frame(time = km_pve$time, survival = km_pve$surv * 100,
             Endpoint = "PVE")
)

ggplot(ep_df, aes(time, survival, colour = Endpoint)) +
  geom_step(linewidth = 0.7) +
  scale_y_continuous(limits = c(0, 100)) +
  scale_colour_manual(values = c("Death" = "#D55E00", "PVE" = "#0072B2")) +
  labs(x = "Months after valve replacement",
       y = "Freedom from event (%)", colour = NULL) +
  theme_minimal() +
  theme(legend.position = "bottom")
```