bio-data-visualization-forest-funnel-plots

What this skill does

## Version Compatibility

Reference examples tested with: metafor 4.4+, forestplot 3.1+, ggforestplot 0.1+ (subgroup forests), ggforest from survminer 0.4.9+, MendelianRandomization 0.10+.

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('<pkg>')` then `?function_name`

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

# Forest and Funnel Plots

**"Summarize effects across studies / subgroups"** -> Render each effect estimate (HR, OR, RR, β) as a square (size = inverse variance / weight), horizontal bar (95% CI), and label, with an optional summary diamond at the bottom from a meta-analysis pool (fixed-effect or random-effects). The funnel plot diagnoses publication bias by plotting effect size vs precision; asymmetry indicates missing small-study-with-null-result publications (Egger 1997).

- R: `metafor::forest`, `metafor::funnel`, `forestplot::forestplot`, `survminer::ggforest` (Cox HR forests), `MendelianRandomization::mr_forest`

## The Single Most Important Modern Insight -- Heterogeneity Is the First Question

A pooled effect estimate is meaningless if the underlying studies are heterogeneous. The Higgins I² statistic (Higgins-Thompson 2002 *Stat Med* 21:1539) quantifies between-study variance; the conventional 25/50/75% interpretation tiers come from the Cochrane Handbook §10.10.2 (Higgins et al editors), NOT the original Higgins-Thompson paper which cautioned against rigid cutoffs. A meta-analysis with I² > 75% and a pooled effect must explain the heterogeneity (subgroup analysis, meta-regression) — pooling without explanation is statistically defensible but biologically unhelpful.

A forest plot's bottom must report: pooled estimate + 95% CI + I² + τ² (between-study variance) + Q-test p-value. Without these, the plot is a list of effects, not a meta-analysis.

## Decision Tree by Analysis Type

| Analysis | Tool | Pooling model | Forest method |
|----------|------|---------------|---------------|
| Single-trial subgroup HRs | survminer::ggforest | None (subgroup display) | Coxph object |
| Meta-analysis of binary outcomes | metafor::rma -> forest() | DerSimonian-Laird or REML random-effects | Standard forest |
| Meta-analysis of continuous outcomes | metafor::rma(yi, vi) | REML random-effects | Standard forest |
| Mendelian randomization | MendelianRandomization::mr_forest | Multiple MR methods | MR-specific forest |
| Subgroup forest with interaction p | metafor::rma + addpoly + interaction model | Subgroup REML | Nested forest |
| Network meta-analysis | netmeta::forest.netmeta | Bayesian or frequentist | Network forest |
| Cumulative meta-analysis (over time) | metafor::cumul + forest | – | Cumulative forest |

## Fixed-Effect vs Random-Effects Meta-Analysis

| Model | Assumption | When appropriate | Pooled estimate weight |
|-------|------------|-------------------|------------------------|
| Fixed-effect (Mantel-Haenszel, IVS) | All studies estimate the same true effect | Single mechanism, homogeneous design | 1 / within-study variance |
| Random-effects (DerSimonian-Laird, REML) | Studies estimate distinct true effects from a common distribution | Heterogeneous designs / populations | 1 / (within + between variance) |

**Use random-effects by default.** Fixed-effect assumes all studies estimate the *same* parameter, which is almost never true across multi-center trials with different populations. REML is the modern default (Viechtbauer 2005); DerSimonian-Laird is the older default still commonly seen.

## metafor::rma + forest -- The Reference Implementation

**Goal:** Pool study-level effect estimates with random-effects meta-analysis; render a forest plot with study weights, individual effect+CI, and pooled summary diamond.

**Approach:** Compute per-study yi (effect) and vi (sampling variance); fit REML random-effects model; pass to forest() with prediction interval if heterogeneity is non-trivial.

```r
library(metafor)

# Input: per-study effect (yi) and variance (vi)
# For OR: yi = log(OR), vi = SE(log(OR))^2
# For HR: yi = log(HR), vi = SE(log(HR))^2
res <- rma(yi = log_or, vi = log_or_se^2,
           data = studies, slab = paste(author, year),
           method = 'REML')

# I^2 and tau^2 in the summary
summary(res)
# I^2 (residual heterogeneity)
# tau^2 (estimated amount of (residual) heterogeneity)
# Q-test for heterogeneity

forest(res,
       atransf = exp,                           # display OR on natural scale
       at = log(c(0.25, 0.5, 1, 2, 4)),         # ticks at meaningful OR values
       refline = 0,                              # log(1) for OR/HR/RR
       xlab = 'Odds Ratio (95% CI)',
       header = c('Study', 'OR [95% CI]'),
       mlab = bquote(paste('RE Model (Q = ', .(round(res$QE, 2)),
                            ', df = ', .(res$k - 1),
                            ', p = ', .(format.pval(res$QEp, digits = 2)),
                            '; ', I^2, ' = ', .(round(res$I2, 1)), '%)')),
       addpred = TRUE)                          # prediction interval per Higgins 2009
```

`addpred = TRUE` adds a 95% prediction interval — where a new study's effect is expected to fall (Higgins-Thompson-Spiegelhalter 2009 *JRSS-A*). This is the most honest summary when I² > 30%.

## ggforest for Cox Subgroup Forests

```r
library(survminer)
fit <- coxph(Surv(time, status) ~ treatment + age + sex + stage, data = df)
ggforest(fit,
         data = df,
         main = 'Subgroup HRs',
         cpositions = c(0.02, 0.22, 0.4),
         fontsize = 0.7,
         refLabel = 'Reference',
         noDigits = 2)
```

ggforest produces a publication-ready subgroup forest from a coxph object. For pre-specified subgroup analyses (treatment × subgroup interaction), test interaction explicitly and annotate the p-value.

## Small-k Regime -- When Meta-Analysis Asymptotics Break

For k < 5 studies, the REML-based 95% CI from `metafor::rma()` is severely anti-conservative — it relies on chi-square asymptotics that fail with few studies. Use **Hartung-Knapp-Sidik-Jonkman (HKSJ)** adjustment (`test = 'knha'`):

```r
res_hksj <- rma(yi = log_or, vi = log_or_se^2, data = studies,
                method = 'REML', test = 'knha')           # HKSJ for k<5
```

HKSJ uses a t-distribution with k-1 degrees of freedom and adjusts SE via the Q-statistic — well-calibrated even at k=3. For k < 3 a meta-analysis is not advisable; report individual study effects in a forest plot without a pooled summary.

Also: I² is uninterpretable below k = 5 (Borenstein 2017 *Res Synth Methods* 8:5); the point estimate has wide CI dominated by k itself, not heterogeneity. **Do not report I² for k < 5.**

## Funnel Plot and Egger Test

**Goal:** Diagnose publication bias by visual asymmetry of effect size vs precision.

**Approach:** Plot effect (x) vs SE (inverted y); under no bias, points form a symmetric inverted funnel with the pooled estimate at the apex. Asymmetry suggests missing small-N null-result studies. Egger's regression test (Egger 1997 *BMJ* 315:629) formalizes the asymmetry.

```r
funnel(res,
       xlab = 'log(OR)',
       refline = res$b)

# Egger's test
regtest(res, model = 'lm', predictor = 'sei')
# significant p indicates asymmetry; suggests publication bias

# Trim-and-fill (Duval-Tweedie 2000) -- adjusts for asymmetry
res_tf <- trimfill(res)
forest(res_tf)
funnel(res_tf)
```

**Contour-enhanced funnel plot** (Peters 2008 *J Clin Epidemiol* 61:991) overlays significance contours (p < 0.10, < 0.05, < 0.01); asymmetry concentrated in "non-significant" regions indicates publication bias more specifically than generic asymmetry.

```r
funnel(res, level = c(90, 95, 99), shade = c('white', 'gray55', 'gray75'),
       refline = 0, legend = TRUE)
```

## Per-Method Failure Modes

### Pooling under high heterogeneity without explanation

**Trigger:** Random-effects meta-analysis pooled with I² >