---
title: "Forest plots for regression models"
format: html
vignette: >
  %\VignetteIndexEntry{Forest plots for regression models}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
knitr:
  opts_chunk:
    collapse: true
    comment: "#>"
    fig.width: 8
    fig.height: 4.5
    out.width: "100%"
---

```{r setup}
#| include: false
library(forrest)
has_broom            <- requireNamespace("broom",            quietly = TRUE)
has_marginaleffects  <- requireNamespace("marginaleffects",  quietly = TRUE)
has_matchit          <- requireNamespace("MatchIt",          quietly = TRUE)
```

Forest plots are not limited to meta-analyses. Any table of point estimates
with confidence intervals can be displayed as a forest plot — results from
regression models, dose-response analyses, multi-outcome studies, and more.

This vignette shows how to go from a fitted R model to a forest plot in a few
lines of code.
[`broom::tidy()`](https://broom.tidymodels.org/reference/tidy.lm.html) extracts
parameters; `broom` is not a hard dependency and may be replaced by any function
that returns a data frame with the right columns.
[`marginaleffects`](https://marginaleffects.com) powers the causal-inference
examples in the final sections.

---

## Simulated cohort

```{r sim-data}
set.seed(42)
n <- 800

cohort <- data.frame(
  age      = round(runif(n, 30, 75)),
  female   = rbinom(n, 1, 0.52),
  bmi      = rnorm(n, 26, 4.5),
  smoker   = rbinom(n, 1, 0.20),
  activity = rbinom(n, 1, 0.45),   # 1 = physically active
  educ     = sample(c("Low", "Medium", "High"), n,
                    replace = TRUE, prob = c(0.25, 0.45, 0.30))
)

# Systolic blood pressure (mmHg): continuous outcome with true signal
cohort$sbp <- 110 +
  0.40 * cohort$age    -
  4.50 * cohort$female +
  0.50 * cohort$bmi    -
  2.00 * cohort$smoker +
  2.50 * cohort$activity +
  rnorm(n, sd = 8)

# Hypertension (SBP ≥ 140): binary outcome
cohort$htn <- as.integer(cohort$sbp >= 140)

# Education as ordered factor (for dose-response example)
cohort$educ <- factor(cohort$educ,
                      levels = c("Low", "Medium", "High"), ordered = FALSE)
cohort$educ <- relevel(cohort$educ, ref = "Low")
```

---

## Multiple predictors from one model

The simplest workflow: fit a model, call `broom::tidy(conf.int = TRUE)`, and
pass the result to `forrest()`. Column names `estimate`, `conf.low`, and
`conf.high` map directly to the `estimate`, `lower`, and `upper` arguments.

```{r single-model}
#| eval: !expr has_broom
#| fig-height: 5
#| fig-width: 9
fit_lm <- lm(sbp ~ age + female + bmi + smoker + activity, data = cohort)

coefs <- broom::tidy(fit_lm, conf.int = TRUE)
coefs <- coefs[coefs$term != "(Intercept)", ]
coefs$term <- c(
  "Age (per 1 year)",
  "Female sex",
  "BMI (per 1 kg/m\u00b2)",
  "Current smoker",
  "Physically active"
)

# Formatted coefficient + CI text for right-hand column
coefs$coef_ci <- sprintf(
  "%.2f (%.2f, %.2f)",
  coefs$estimate, coefs$conf.low, coefs$conf.high
)

forrest(
  coefs,
  estimate = "estimate",
  lower    = "conf.low",
  upper    = "conf.high",
  label    = "term",
  header   = "Predictor",
  cols     = c("Coef (95% CI)" = "coef_ci"),
  widths   = c(3, 4, 2.5),
  xlab     = "Regression coefficient (95% CI)",
  stripe   = TRUE
)
```

For **logistic regression**, set `exponentiate = TRUE` in `broom::tidy()` to
get odds ratios directly, then set `log_scale = TRUE` and `ref_line = 1` in
`forrest()`:

```{r logistic}
#| eval: !expr has_broom
#| fig-height: 5
fit_glm <- glm(htn ~ age + female + bmi + smoker + activity,
               data = cohort, family = binomial)

or_coefs <- broom::tidy(fit_glm, conf.int = TRUE, exponentiate = TRUE)
or_coefs <- or_coefs[or_coefs$term != "(Intercept)", ]
or_coefs$term <- c(
  "Age (per 1 year)",
  "Female sex",
  "BMI (per 1 kg/m\u00b2)",
  "Current smoker",
  "Physically active"
)

or_coefs$or_text <- sprintf(
  "%.2f (%.2f, %.2f)",
  or_coefs$estimate, or_coefs$conf.low, or_coefs$conf.high
)

forrest(
  or_coefs,
  estimate   = "estimate",
  lower      = "conf.low",
  upper      = "conf.high",
  label      = "term",
  log_scale  = TRUE,
  ref_line   = 1,
  header     = "Predictor",
  cols       = c("OR (95% CI)" = "or_text"),
  widths     = c(3, 4, 2),
  xlab       = "Odds ratio (95% CI)"
)
```

---

## Same predictor across progressive adjustment models

A common reporting pattern is to show how an association changes as
covariates are progressively added.

```{r adjustment-models}
#| eval: !expr has_broom
#| fig-height: 4
m1 <- glm(htn ~ activity,
          data = cohort, family = binomial)
m2 <- glm(htn ~ activity + age + female,
          data = cohort, family = binomial)
m3 <- glm(htn ~ activity + age + female + bmi + smoker,
          data = cohort, family = binomial)

extract_activity <- function(mod, label) {
  row         <- broom::tidy(mod, conf.int = TRUE, exponentiate = TRUE)
  row         <- row[row$term == "activity", ]
  row$model   <- label
  row
}

adj_dat <- rbind(
  extract_activity(m1, "Model 1: unadjusted"),
  extract_activity(m2, "Model 2: + age, sex"),
  extract_activity(m3, "Model 3: + BMI, smoking")
)

adj_dat$or_text <- sprintf(
  "%.2f (%.2f, %.2f)",
  adj_dat$estimate, adj_dat$conf.low, adj_dat$conf.high
)

forrest(
  adj_dat,
  estimate   = "estimate",
  lower      = "conf.low",
  upper      = "conf.high",
  label      = "model",
  log_scale  = TRUE,
  ref_line   = 1,
  header     = "Adjustment set",
  cols       = c("OR (95% CI)" = "or_text"),
  widths     = c(4, 3.5, 2),
  xlab       = "OR for physical activity on hypertension (95% CI)"
)
```

---

## Same predictor across multiple outcomes

To compare one exposure against several outcomes, loop over outcomes, extract
the exposure row from each model, and stack into a single data frame. Use
`group` to colour by outcome domain.

```{r multi-outcome}
#| eval: !expr has_broom
#| fig-height: 6
outcomes <- list(
  list(var = "sbp",      label = "Systolic BP (mmHg)",    domain = "Continuous"),
  list(var = "htn",      label = "Hypertension",          domain = "Binary"),
  list(var = "bmi",      label = "BMI (kg/m\u00b2)",      domain = "Continuous"),
  list(var = "smoker",   label = "Current smoker",        domain = "Binary")
)

rows <- lapply(outcomes, function(o) {
  if (o$domain == "Continuous") {
    fit <- lm(as.formula(paste(o$var, "~ activity + age + female")),
              data = cohort)
    r   <- broom::tidy(fit, conf.int = TRUE)
  } else {
    fit <- glm(as.formula(paste(o$var, "~ activity + age + female")),
               data = cohort, family = binomial)
    r   <- broom::tidy(fit, conf.int = TRUE, exponentiate = TRUE)
  }
  r          <- r[r$term == "activity", ]
  r$outcome  <- o$label
  r$domain   <- o$domain
  r
})

multi_out <- do.call(rbind, rows)
multi_out$est_text <- sprintf(
  "%.2f (%.2f, %.2f)",
  multi_out$estimate,
  multi_out$conf.low,
  multi_out$conf.high
)

forrest(
  multi_out,
  estimate   = "estimate",
  lower      = "conf.low",
  upper      = "conf.high",
  label      = "outcome",
  group      = "domain",
  ref_line   = 0,
  legend_pos = "topright",
  header     = "Outcome",
  cols       = c("Estimate (95% CI)" = "est_text"),
  widths     = c(3.5, 3.5, 2.5),
  xlab       = "Effect of physical activity (coef or OR, 95% CI)"
)
```

> **Note:** because continuous and binary outcomes live on different scales,
> a single reference line may not be ideal for all rows. Consider splitting
> into two separate panels, or use `ref_line = NULL`.

---

## Dose-response (exposure categories)

For categorical exposures (quartiles, tertiles, and so on), the reference category
has `estimate = NA` because there is no contrast to estimate. `forrest()`
renders this row without a CI or point but keeps the label, making the
reference position visually clear. Supply a text column to print `"1.00 (ref)"`
in the right-hand panel. Set `ref_label = TRUE` to have `forrest()` append
`" (Ref.)"` to the label automatically.

```{r dose-response}
#| eval: !expr has_broom
#| fig-height: 5
# Quartiles of age
cohort$age_q <- cut(
  cohort$age,
  breaks         = quantile(cohort$age, 0:4 / 4),
  include.lowest = TRUE,
  labels         = c("Q1 (\u226445 y)", "Q2 (45\u201356 y)",
                     "Q3 (56\u201366 y)", "Q4 (>66 y)")
)
cohort$age_q <- relevel(cohort$age_q, ref = "Q1 (\u226445 y)")

fit_q   <- glm(htn ~ age_q + female + bmi + smoker + activity,
               data = cohort, family = binomial)
q_coefs <- broom::tidy(fit_q, conf.int = TRUE, exponentiate = TRUE)
q_coefs <- q_coefs[grep("^age_q", q_coefs$term), ]

q_plot <- rbind(
  data.frame(
    label     = "Q1 (\u226445 y, ref)",
    estimate  = NA_real_,
    conf.low  = NA_real_,
    conf.high = NA_real_,
    or_ci     = "1.00 (ref)",
    is_sum    = FALSE
  ),
  data.frame(
    label     = gsub("^age_q", "", q_coefs$term),
    estimate  = q_coefs$estimate,
    conf.low  = q_coefs$conf.low,
    conf.high = q_coefs$conf.high,
    or_ci     = sprintf("%.2f (%.2f, %.2f)",
                        q_coefs$estimate,
                        q_coefs$conf.low,
                        q_coefs$conf.high),
    is_sum    = FALSE
  )
)

forrest(
  q_plot,
  estimate   = "estimate",
  lower      = "conf.low",
  upper      = "conf.high",
  label      = "label",
  is_summary = "is_sum",
  log_scale  = TRUE,
  ref_line   = 1,
  header     = "Age quartile",
  cols       = c("OR (95% CI)" = "or_ci"),
  widths     = c(3, 4, 2),
  xlab       = "OR for hypertension (95% CI)"
)
```

---

## Marginal effects: g-computation

Parametric g-computation estimates the *average treatment effect* (ATE) by
predicting outcomes under two counterfactual worlds — everyone treated vs. no
one treated — and averaging the difference.
[`marginaleffects::avg_comparisons()`](https://marginaleffects.com/chapters/gcomputation.html)
implements this in one line.

The example below estimates the marginal effect of physical activity on SBP,
overall and stratified by sex and education level.

```{r gcomputation}
#| eval: !expr (has_marginaleffects && has_broom)
#| fig-height: 6
#| fig-width: 10
library(marginaleffects)

# Outcome model with treatment × covariate interactions so that
# g-computation captures heterogeneous treatment effects
fit_sbp <- lm(
  sbp ~ activity * (age + female + bmi + smoker),
  data = cohort
)

# Helper: extract one avg_comparisons result into a plot row
make_row <- function(me, label, is_sum = FALSE) {
  data.frame(
    label    = label,
    estimate = me$estimate,
    lower    = me$conf.low,
    upper    = me$conf.high,
    is_sum   = is_sum
  )
}

# Overall ATE
ate_sbp <- avg_comparisons(fit_sbp, variables = "activity")

# ATE by sex: subset newdata to each stratum
sex_rows <- lapply(
  list(`  Female` = 0, `  Male` = 1),
  function(f) avg_comparisons(fit_sbp, variables = "activity",
                               newdata = subset(cohort, female == f))
)

# ATE by education level
educ_rows <- setNames(
  lapply(c("Low", "Medium", "High"), function(e)
    avg_comparisons(fit_sbp, variables = "activity",
                    newdata = subset(cohort, educ == e))),
  c("  Low", "  Medium", "  High")
)

# Stack into a single data frame with a stratum column for section headers.
# No manual header or spacer rows needed.
gcomp_dat <- rbind(
  make_row(ate_sbp, "Overall (ATE)", is_sum = TRUE) |>
    transform(stratum = "Overall"),
  do.call(rbind, lapply(
    list(list(label = "Female", f = 0), list(label = "Male", f = 1)),
    function(s) {
      r <- make_row(avg_comparisons(fit_sbp, variables = "activity",
                                    newdata = subset(cohort, female == s$f)),
                   label = s$label)
      r$stratum <- "Sex"
      r
    }
  )),
  do.call(rbind, lapply(
    list(list(label = "Low", e = "Low"),
         list(label = "Medium", e = "Medium"),
         list(label = "High",   e = "High")),
    function(s) {
      r <- make_row(avg_comparisons(fit_sbp, variables = "activity",
                                    newdata = subset(cohort, educ == s$e)),
                   label = s$label)
      r$stratum <- "Education level"
      r
    }
  ))
)

gcomp_dat$est_text <- sprintf(
  "%.2f (%.2f, %.2f)",
  gcomp_dat$estimate, gcomp_dat$lower, gcomp_dat$upper
)

forrest(
  gcomp_dat,
  estimate   = "estimate",
  lower      = "lower",
  upper      = "upper",
  label      = "label",
  section    = "stratum",
  is_summary = "is_sum",
  cols       = c("Marginal effect (95% CI)" = "est_text"),
  widths     = c(3.5, 4, 3),
  ref_line   = 0,
  header     = "Stratum",
  xlab       = "Marginal effect on SBP (mmHg): physical activity (g-computation)"
)
```

---

## Marginal effects: inverse probability weighting (IPW)

IPW re-weights observations by the inverse of their propensity score so that
the weighted sample resembles a randomised trial with respect to observed
confounders.
The example below uses
[`marginaleffects`](https://marginaleffects.com/bonus/ipw.html) to implement
the Hajek estimator and compares it to g-computation for five strata (overall
+ sex + education). Using `group` and `dodge`, both sets of estimates are
displayed side-by-side in one figure.

```{r ipw}
#| eval: !expr (has_marginaleffects && has_broom)
#| fig-height: 5
#| fig-width: 10
# ── Step 1: propensity score model ────────────────────────────────────────
ps_mod       <- glm(activity ~ age + female + bmi + smoker + educ,
                    data = cohort, family = binomial)
cohort$ps    <- predict(ps_mod, type = "response")
cohort$ipw   <- ifelse(cohort$activity == 1,
                       1 / cohort$ps,
                       1 / (1 - cohort$ps))

# Trim extreme weights at the 1st / 99th percentile for numerical stability
cohort$ipw_t <- pmin(pmax(cohort$ipw,
                           quantile(cohort$ipw, 0.01)),
                     quantile(cohort$ipw, 0.99))

# ── Step 2: IPW-weighted outcome model ────────────────────────────────────
fit_sbp_ipw <- lm(
  sbp ~ activity * (age + female + bmi + smoker),
  data = cohort, weights = ipw_t
)

# ── Step 3: build comparison rows — no header rows, just labelled strata ──
# Each stratum appears twice: once for each method.  The dodge layout
# places the two symbols side-by-side so CIs can be read off the same axis.
strata <- list(
  list(label = "Overall",         nd = cohort,                       is_sum = TRUE),
  list(label = "  Female",        nd = subset(cohort, female == 0),  is_sum = FALSE),
  list(label = "  Male",          nd = subset(cohort, female == 1),  is_sum = FALSE),
  list(label = "  Low educ.",     nd = subset(cohort, educ == "Low"), is_sum = FALSE),
  list(label = "  High educ.",    nd = subset(cohort, educ == "High"),is_sum = FALSE)
)

rows <- do.call(rbind, lapply(strata, function(s) {
  # G-computation: standard avg_comparisons on the interaction model
  gc  <- avg_comparisons(fit_sbp,     variables = "activity",
                          newdata = s$nd)
  # IPW (Hajek): weighted outcome model + Hajek normalisation via `wts`
  ipw <- avg_comparisons(fit_sbp_ipw, variables = "activity",
                          wts = "ipw_t", newdata = s$nd)
  rbind(
    data.frame(label = s$label, method = "G-computation",
               estimate = gc$estimate,  lower = gc$conf.low,
               upper    = gc$conf.high, is_sum = s$is_sum),
    data.frame(label = s$label, method = "IPW (Hajek)",
               estimate = ipw$estimate, lower = ipw$conf.low,
               upper    = ipw$conf.high, is_sum = s$is_sum)
  )
}))

# Separate text columns per method for cols_by_group display
rows$text_gc  <- ifelse(rows$method == "G-computation",
                         sprintf("%.2f (%.2f, %.2f)",
                                 rows$estimate, rows$lower, rows$upper), "")
rows$text_ipw <- ifelse(rows$method == "IPW (Hajek)",
                         sprintf("%.2f (%.2f, %.2f)",
                                 rows$estimate, rows$lower, rows$upper), "")

forrest(
  rows,
  estimate      = "estimate",
  lower         = "lower",
  upper         = "upper",
  label         = "label",
  group         = "method",
  is_summary    = "is_sum",
  dodge         = TRUE,
  cols_by_group = TRUE,
  cols          = c("G-comp (95% CI)" = "text_gc",
                    "IPW (95% CI)"    = "text_ipw"),
  widths        = c(3, 3.5, 2.5, 2.5),
  ref_line      = 0,
  legend_pos    = "topright",
  header        = "Subgroup",
  xlab          = "Marginal effect on SBP (mmHg): physical activity (95% CI)"
)
```

---

## Time-varying exposures and longitudinal studies

Any analysis that produces estimates at multiple time points — distributed lag
models, life-course analyses, repeated measures, cumulative exposure windows —
yields a data frame where time is just another dimension. The same
`section` / `dodge` / `group` building blocks that handle predictor groups and
subgroup analyses apply here without modification.

### Distributed lag models

A distributed lag model (DLM) estimates the association between an exposure and
an outcome separately at each lag. The result is one row per lag per exposure.
Use `section = "exposure"` to group lags under each exposure, and add an
overall (cumulative) estimate as an `is_summary = TRUE` diamond at the bottom
of each section.

```{r dlm}
#| fig-height: 9
set.seed(2025)

# Simulate lag-specific estimates for three environmental exposures
# Each exposure has a distinct lag-response shape
make_lags <- function(pattern) {
  se  <- runif(length(pattern), 0.03, 0.06)
  data.frame(
    estimate = pattern + rnorm(length(pattern), sd = 0.02),
    lower    = pattern - 1.96 * se,
    upper    = pattern + 1.96 * se
  )
}

pm25  <- make_lags(c(0.04, 0.11, 0.10, 0.06, 0.03,  0.01,  0.00))
noise <- make_lags(c(0.09, 0.05,  0.02, 0.01, 0.00, -0.01,  0.00))
green <- make_lags(c(-0.02, -0.04, -0.07, -0.09, -0.10, -0.09, -0.08))

lag_labels <- paste("Lag", 0:6, "(weeks)")

dlm_dat <- rbind(
  data.frame(exposure = "PM2.5",        lag = lag_labels, pm25,  is_sum = FALSE),
  data.frame(exposure = "PM2.5",        lag = "Cumulative", estimate = 0.35,
             lower = 0.18, upper = 0.52, is_sum = TRUE),
  data.frame(exposure = "Road noise",   lag = lag_labels, noise, is_sum = FALSE),
  data.frame(exposure = "Road noise",   lag = "Cumulative", estimate = 0.16,
             lower = 0.04, upper = 0.28, is_sum = TRUE),
  data.frame(exposure = "Green space",  lag = lag_labels, green, is_sum = FALSE),
  data.frame(exposure = "Green space",  lag = "Cumulative", estimate = -0.49,
             lower = -0.68, upper = -0.30, is_sum = TRUE)
)

dlm_dat$est_text <- ifelse(
  dlm_dat$is_sum,
  sprintf("%.2f (%.2f, %.2f)", dlm_dat$estimate, dlm_dat$lower, dlm_dat$upper),
  sprintf("%.2f (%.2f, %.2f)", dlm_dat$estimate, dlm_dat$lower, dlm_dat$upper)
)

forrest(
  dlm_dat,
  estimate   = "estimate",
  lower      = "lower",
  upper      = "upper",
  label      = "lag",
  section    = "exposure",
  is_summary = "is_sum",
  ref_line   = 0,
  header     = "Lag",
  cols       = c("Coef (95% CI)" = "est_text"),
  widths     = c(3.5, 4, 2.8),
  xlab       = "Regression coefficient per IQR increase (95% CI)"
)
```

### Life-course analysis: same exposure across developmental periods

A life-course design estimates the same exposure–outcome association at multiple
developmental stages. Use `section = "life_stage"` to organise the stages as
named groups, with one row per exposure within each stage.

```{r life-course}
#| fig-height: 7
set.seed(99)

make_row <- function(exposure, life_stage, est, se) {
  data.frame(
    life_stage = life_stage,
    exposure   = exposure,
    estimate   = est  + rnorm(1, sd = 0.01),
    lower      = est  - 1.96 * se,
    upper      = est  + 1.96 * se
  )
}

lc_dat <- rbind(
  make_row("Noise exposure",         "Prenatal",        0.08, 0.03),
  make_row("Green space access",     "Prenatal",       -0.05, 0.04),
  make_row("Traffic-related air",    "Prenatal",        0.06, 0.03),
  make_row("Noise exposure",         "Early childhood", 0.12, 0.04),
  make_row("Green space access",     "Early childhood",-0.09, 0.05),
  make_row("Traffic-related air",    "Early childhood", 0.10, 0.04),
  make_row("Noise exposure",         "Mid-childhood",   0.15, 0.05),
  make_row("Green space access",     "Mid-childhood",  -0.13, 0.05),
  make_row("Traffic-related air",    "Mid-childhood",   0.14, 0.05),
  make_row("Noise exposure",         "Adolescence",     0.09, 0.04),
  make_row("Green space access",     "Adolescence",    -0.07, 0.05),
  make_row("Traffic-related air",    "Adolescence",     0.08, 0.04)
)

forrest(
  lc_dat,
  estimate = "estimate",
  lower    = "lower",
  upper    = "upper",
  label    = "exposure",
  section  = "life_stage",
  group    = "exposure",
  ref_line = 0,
  header   = "Exposure",
  xlab     = "Coefficient per IQR increase in SBP (95% CI)"
)
```

### Life-course analysis: comparing stages side by side (dodge)

When the question is *how does the association at one developmental stage
compare to another for the same exposure*, flip the structure: `label` is the
exposure, `group` is the life stage, and `dodge = TRUE` places all stages
side by side within each exposure row.

```{r life-course-dodge}
#| fig-height: 5
#| fig-width: 10
# Same data — one row per exposure × life stage
forrest(
  lc_dat,
  estimate = "estimate",
  lower    = "lower",
  upper    = "upper",
  label    = "exposure",
  group    = "life_stage",
  dodge    = TRUE,
  ref_line = 0,
  header   = "Exposure",
  xlab     = "Coefficient per IQR increase in SBP (95% CI)"
)
```

---

## Saving plots

```{r save}
#| eval: false
save_forrest(
  file   = "regression_forest.pdf",
  plot   = function() {
    forrest(
      coefs,
      estimate = "estimate",
      lower    = "conf.low",
      upper    = "conf.high",
      label    = "term",
      stripe   = TRUE,
      xlab     = "Regression coefficient (95% CI)"
    )
  },
  width  = 9,
  height = 5
)
```