One model’s coefficients

ggforest() draws the hazard ratio and 95% confidence interval for every term of a single Cox model, with the reference level, N, and p-value laid out beside each row and a global model summary in the caption:

ggforest(cox, data = colon2)

Reading a forest plot. The vertical reference line at HR = 1 marks “no effect”: a box to its left (HR < 1) means that level lowers the hazard relative to the reference, to its right (HR > 1) means it raises the hazard. The box is the point estimate and the whiskers its 95% confidence interval; a whisker that crosses 1 is not statistically significant. The box size grows with the precision of the estimate. (The same reference line appears in every forest below.)

ggforest() has arguments to tailor the layout – columns chooses which statistic columns appear (any subset of c("N", "hr", "ci", "p")), and palette colours the points and intervals (by variable, or a single accent) instead of the default black.

Is the effect the same in every subgroup?

A single model gives one averaged treatment effect. What clinicians publish is the treatment effect within each subgroup, next to a test of whether the effect actually differs across a subgroup’s levels. ggforest_subgroup() draws exactly that:

ggforest_subgroup(
  cox, data = colon2,
  treatment = "rx",
  subgroups = c(Sex = "sex", "Age group" = "agegrp", Differentiation = "differ")
)

Read it in two parts, because it combines two different models:

• each row is the treatment hazard ratio estimated within that subgroup level (a Cox model refit on that subset), so you can see where the benefit is larger or smaller;
• each header carries one interaction p-value for the whole subgroup variable – from a single model with a treatment × subgroup term (a likelihood-ratio test), not one p per level. A small interaction p is evidence of genuine effect modification; a large one says the apparent differences between levels are within noise.

This two-model structure is why a subgroup forest cannot be read off a single ggforest(). Two details worth stating: the per-level hazard ratios are the treatment-only effect within each subset (not adjusted for the other covariates), and the No. (%) column counts the subgroup size – so “Overall” here is 619, the full cohort, which is larger than the 594 complete cases the multivariable ggforest() above was fit on.

One covariate across several models

When you add covariates to a model, the estimate for a key variable moves. The classic reporting device – unadjusted versus adjusted, or naive versus time-varying – puts that one covariate’s hazard ratio on one row per model. ggforest_models() does this:

cox_crude <- coxph(Surv(time, status) ~ rx, data = colon2)

ggforest_models(
  list("Unadjusted" = cox_crude, "Adjusted" = cox),
  term = "rx"
)

Here the treatment hazard ratio moves little after adjustment (0.69 unadjusted → 0.66 adjusted), consistent with limited confounding by the measured covariates – though part of the small shift is a change of sample, since the adjusted model drops the few cases with missing covariates. The hazard ratio, confidence interval and p-value are read from each model with a robust variance when one is present, and a model missing the term is dropped with a warning.

Can you trust the hazard ratio? The proportional-hazards check

Every hazard ratio above assumes proportional hazards – that a covariate’s effect is a constant multiplier over all follow-up. When that fails, the single HR is a time-average that can hide a treatment whose benefit grows, shrinks, or reverses. ggcoxnph() gathers the standard non-proportional-hazards diagnostics into one panel. We use the veteran lung-cancer trial and a performance-status split (Karnofsky ≤ 60 vs > 60), a textbook example where the effect is strongly time-varying:

veteran$perf <- factor(ifelse(veteran$karno <= 60, "Poor", "Good"),
                       levels = c("Poor", "Good"))
cox_perf <- coxph(Surv(time, status) ~ perf, data = veteran)

ggcoxnph(cox_perf, variable = "perf", palette = "jco")

• A. Complementary log-log survival by group: parallel curves support proportional hazards; here they clearly converge, warning against it.
• B. Scaled Schoenfeld residuals with a smooth of β(t) and the Grambsch-Therneau p-value: a flat line (large p) is consistent with proportional hazards; here β(t) drifts and p = 0.00032, so PH is rejected.
• C. Hazard-ratio trend exp(β(t)): the shape of the effect over time – here the protective hazard ratio sits far below 1 early and drifts back toward it, so the benefit of good performance status is largest early and attenuates over follow-up (read it as shape, not magnitude).
• D. RMST difference versus horizon: the absolute benefit as the follow-up window lengthens – a proportional-hazards-free view that previews the restricted mean survival time from the Absolute effect measures vignette.

When the assumption fails

That is exactly this case: panel B rejects proportional hazards (p = 0.00032) and panel A’s curves converge, so a single hazard ratio for performance status is not a faithful summary – report an absolute effect measure instead. The Absolute effect measures vignette covers restricted mean survival time, milestone survival, and landmark analysis for this situation.

Where to next

• Get started – presets, themes, and the composable curve.
• Absolute effect measures – RMST, milestone, and landmark analysis for when proportional hazards fails.
• Competing risks – cumulative incidence when subjects can fail in more than one way.