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()` 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` powers the causal-inference examples in the final sections.

## Simulated cohort

    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.

    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()`:

    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.

    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.

    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.

    # 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()` 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.

    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` 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.

    # ── 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.

    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.

    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.

    # 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

    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
    )