Meta Analysis

Meta-Analysis

Systematic review and meta-analysis for clinical and preclinical evidence. Covers protocol registration, search strategy, screening, PRISMA 2020 flow diagrams, risk of bias, pooling with fixed and random effects models, small-study effects, sensitivity diagnostics, network meta-analysis, and GRADE/CINeMA certainty rating. Uses metafor, netmeta, PRISMA2020, synthesisr, and robvis.

When to Use This Skill

Activate when the user requests:

• A systematic review or meta-analysis protocol
• Search strategy construction for PubMed, Embase, or Cochrane CENTRAL
• PRISMA 2020 flow diagram generation
• Deduplication of records across databases
• Title/abstract or full-text screening workflows
• Inter-rater agreement between screeners
• Risk of bias assessment (RoB 2, ROBINS-I, ROBINS-E, QUADAS-2, Newcastle-Ottawa)
• Risk of bias visualization (traffic light or summary plots)
• Data extraction templates for RCTs or observational studies
• Effect size computation (OR, RR, RD, SMD, MD, ROM, HR)
• Pooling with equal-effects or random-effects models
• Heterogeneity quantification (tau^2, I^2, Q) and prediction intervals
• Subgroup analysis and meta-regression
• Hazard ratio reconstruction from published Kaplan-Meier curves
• Forest and funnel plots
• Publication bias or small-study effect assessment
• Leave-one-out, influence, or cumulative sensitivity analysis
• Network meta-analysis, transitivity, inconsistency, or treatment ranking
• GRADE or CINeMA certainty rating

Inputs

Environment

# Core, all actively maintained (versions verified 2026-08)
install.packages(c("metafor", "meta", "PRISMA2020", "robvis", "synthesisr", "irr"))
#   metafor     5.0-1    effect sizes and models
#   meta        8.5-0    alternative interface, GRADE-friendly output
#   PRISMA2020  1.1.4    flow diagrams
#   robvis      0.3.1    risk of bias plots
#   synthesisr  0.4.1    bibliographic import and deduplication
#   irr         0.85     Cohen's and Fleiss' kappa

# Network meta-analysis and bias sensitivity
install.packages(c("netmeta", "metasens", "gemtc", "multinma"))
#   netmeta     3.6-1    frequentist NMA. NOTE: pairwise() is in `meta`
#   metasens    1.5-3    Copas selection model, limit meta-analysis
#   gemtc       1.1-1    Bayesian NMA, JAGS backend
#   multinma    0.9.1    Bayesian NMA, Stan; aggregate + individual patient data

# Programmatic search
install.packages(c("rentrez", "easyPubMed"))
#   rentrez     1.2.4
#   easyPubMed  3.1.6

Packages to avoid, and why:

  revtools    0.4.1, last released 2019-12. Superseded by synthesisr, which is
              by the same authors and actively maintained.
  metagear    0.7, last released 2021-02.
  metaviz     0.3.1, last released 2020-04. metafor and meta cover the plots.
  esc         0.5.1, last released 2019-12. Use metafor::escalc().
  pcnetmeta   2.8, last released 2022-08.
  nmadb       1.2.0, last released 2019-12.

Not an R package at all:
  CINeMA      web application for certainty in network meta-analysis.
              Do not search CRAN for it.

Not on CRAN, install from source if needed:
  ASySD       deduplication, higher sensitivity than reference-manager dedup
  litsearchr  search term discovery from a naive search
  dmetar      companion package to the Doing Meta-Analysis book

Protocol and Registration

Register before screening starts. A protocol written after seeing results is not a protocol.

Where to register:
  PROSPERO        health-related reviews with a health outcome. Free.
                  Registration before data extraction begins.
  OSF Registries  anything PROSPERO will not take (preclinical, methods,
                  animal studies, scoping reviews).
  INPLASY         alternative when PROSPERO turnaround is too slow.

What must be pre-specified, because changing it later is a protocol deviation
that has to be reported:
  - PICO elements and the review question
  - Eligibility criteria, including study designs and language limits
  - Databases and the planned search date range
  - Primary and secondary outcomes, defined precisely
  - Effect measure (OR, RR, HR, MD, SMD)
  - Synthesis model (fixed vs random effects) and heterogeneity handling
  - Subgroup and sensitivity analyses, declared in advance
  - Risk of bias tool

Report the PROSPERO ID in the manuscript. Reviewers check it, and deviations between the registration and the paper are a common reason for rejection.

Search Strategy

Translating PICO into a query

PICO -> query blocks, combined with AND across blocks and OR within blocks:

  P (Population)    disease terms, controlled vocabulary + free text
  I (Intervention)  drug/exposure names, including synonyms and brand names
  C (Comparator)    usually omitted from the search; it over-restricts
  O (Outcome)       usually omitted; outcomes are poorly indexed and
                    including them loses relevant records

Standard practice is to search P AND I only, and apply C and O at screening.
Searching all four is the most common cause of a search that misses studies.

Controlled vocabulary and free text

Every block needs both. MeSH alone misses recent records that have not been indexed yet; free text alone misses records where the concept is only in the MeSH.

PubMed
  "Carcinoma, Non-Small-Cell Lung"[Mesh] OR nsclc[tiab] OR
  "non small cell lung"[tiab]

  [Mesh]        controlled vocabulary, auto-explodes to narrower terms
  [Mesh:NoExp]  suppress explosion
  [tiab]        title/abstract
  [tw]          text word, broader than tiab
  [pt]          publication type

Embase (Ovid syntax)
  exp Carcinoma, Non-Small Cell Lung/ OR nsclc.ti,ab,kw.

  exp   explode the Emtree term
  /     Emtree term marker
  .ti,ab,kw.  title, abstract, keyword

Cochrane CENTRAL
  [mh "Carcinoma, Non-Small-Cell Lung"] OR nsclc:ti,ab,kw

Truncation differs by platform. PubMed uses * and requires at least four characters before it; Ovid uses $ or *. PubMed does not support left truncation.

Validated study design filters

Do not write your own RCT filter. Use a validated one and cite it.

Cochrane Highly Sensitive Search Strategy (CHSSS), sensitivity-maximizing
version for PubMed. Published in the Cochrane Handbook, chapter 4.

  (randomized controlled trial[pt] OR controlled clinical trial[pt] OR
   randomized[tiab] OR placebo[tiab] OR clinical trials as topic[mesh:noexp] OR
   randomly[tiab] OR trial[ti]) NOT (animals[mh] NOT humans[mh])

Note the animal exclusion is written as NOT (animals NOT humans), not
NOT animals[mh]. The latter drops human studies that also used animal models.

For observational designs there is no equivalent gold standard. Filters exist but lose sensitivity, so most reviews search without a design filter and exclude at screening.

Programmatic PubMed search

library(rentrez)   # v1.2.4

query <- paste(
  '("Carcinoma, Non-Small-Cell Lung"[Mesh] OR nsclc[tiab])',
  'AND ("Immunotherapy"[Mesh] OR pembrolizumab[tiab] OR nivolumab[tiab])',
  'AND ("2015/01/01"[PDAT] : "2026/12/31"[PDAT])'
)

# use_history keeps results server-side; required above ~10k records
res <- entrez_search(db = "pubmed", term = query, use_history = TRUE, retmax = 0)
res$count

# Fetch in batches. NCBI throttles to 3 requests/second without an API key,
# 10/second with one. Set it via ENTREZ_KEY or set_entrez_key().
recs <- character()
for (start in seq(0, res$count - 1, by = 200)) {
  recs <- c(recs, entrez_fetch(
    db = "pubmed", web_history = res$web_history,
    rettype = "medline", retmode = "text",
    retstart = start, retmax = 200
  ))
  Sys.sleep(0.34)
}
writeLines(recs, "pubmed_records.nbib")

The search string, the database, the platform, the date run, and the number of hits must all be recorded per database. PRISMA 2020 requires the full strategy for at least one database in the manuscript or supplement.

Deduplication

Records overlap heavily across databases. PubMed and Embase alone typically overlap 40-60%.

library(synthesisr)   # v0.4.1

files <- c("pubmed.nbib", "embase.ris", "central.ris")
refs <- read_refs(files, tag_naming = "best_guess", return_df = TRUE)

# Exact match on DOI first: fast and safe
refs <- refs[!duplicated(refs$doi) | is.na(refs$doi), ]

# Then fuzzy match on title for records lacking a DOI
dups <- find_duplicates(
  refs$title,
  method = "string_osa",
  to_lower = TRUE, rm_punctuation = TRUE,
  threshold = 7          # max edit distance; raise to catch more, at the
)                        # cost of false merges
deduped <- extract_unique_references(refs, matches = dups)

nrow(refs) - nrow(deduped)   # duplicates removed, needed for PRISMA

Deduplication is not a solved problem. Two failure modes, opposite directions:

  Under-merging   the same trial published as a conference abstract and a
                  full paper has different titles and no shared DOI. These
                  are duplicate STUDIES, not duplicate RECORDS, and must be
                  linked at data extraction, not here.

  Over-merging    companion papers reporting different outcomes of one trial
                  have near-identical titles. Merging them loses an outcome.

Always eyeball the merged pairs before accepting. A threshold that removes
"too many" duplicates is worse than one that removes too few, because the
lost records are invisible downstream.

Reference-manager deduplication (EndNote, Zotero) has lower sensitivity than dedicated tools. If dedup accuracy matters, ASySD reports sensitivity 0.95-0.99 with specificity above 0.99.

Screening

Two reviewers, independently

Title/abstract screening
  Two reviewers screen all records independently against the eligibility
  criteria. Liberal inclusion: if either reviewer says maybe, it advances.
  Reconcile disagreements by discussion, with a third reviewer to break ties.

Full-text screening
  Same two-reviewer process. This is the stage where every exclusion needs a
  recorded REASON, because PRISMA 2020 requires reporting them with counts.

Pilot first
  Both reviewers screen the same 50-100 records, compare, and refine the
  criteria before screening the rest. Most criteria ambiguity surfaces here.

Measuring agreement

library(irr)   # v0.85

# screening: data.frame with one column per reviewer, one row per record
kappa2(screening[, c("reviewer_1", "reviewer_2")])       # two reviewers
kappam.fleiss(screening[, c("r1", "r2", "r3")])          # three or more

Interpreting kappa for screening:
  < 0.40   poor. The criteria are ambiguous. Stop and rewrite them.
  0.40-0.60  moderate. Usually fixable with a calibration round.
  0.60-0.80  substantial. Acceptable.
  > 0.80   almost perfect.

Kappa is deflated when inclusion is rare, which it always is at title/abstract
(typical inclusion 2-5%). A low kappa with high raw agreement is the expected
pattern, not necessarily a problem. Report both.

Screening automation

Active-learning tools rank records by predicted relevance so screening can stop early. They are decision aids, not replacements for a second reviewer.

ASReview    open source, active learning. Screens in relevance order and
            plateaus once the recall curve flattens.
Rayyan      web based, free tier, supports blinded two-reviewer workflow.

If you use one, report it: the tool, the model, the stopping rule, and
whether a human screened every record or only until the stopping rule fired.
An unreported stopping rule is a reproducibility gap.

PRISMA 2020 Flow Diagram

PRISMA 2020 replaced the 2009 statement and changed the diagram structure. The current version separates records identified from databases and registers, and distinguishes reports from studies.

library(PRISMA2020)   # v1.1.4

# The package expects a specific set of row names. Start from the template
# shipped with the package rather than building the data frame by hand.
template <- system.file("extdata", "PRISMA.csv", package = "PRISMA2020")
counts <- read.csv(template)

# Key fields, all as integers:
#   database_results, register_results       identification
#   duplicates, excluded_automatic, excluded_other
#   records_screened, records_excluded       screening
#   dbr_sought_reports, dbr_notretrieved_reports
#   dbr_assessed, dbr_excluded               eligibility, with reasons
#   new_studies, new_reports                 included

data <- PRISMA_data(counts)

plot <- PRISMA_flowdiagram(
  data,
  interactive = FALSE,
  previous = FALSE,        # TRUE only for a review update with prior studies
  other = TRUE,            # records found outside database searching
  detail_databases = TRUE, # break identification down per database
  side_boxes = TRUE
)

PRISMA_save(plot, filename = "prisma_flow.pdf", filetype = "PDF", overwrite = TRUE)

The arithmetic must reconcile, and reviewers check it:

  records_screened = (database_results + register_results)
                     - duplicates - excluded_automatic - excluded_other

  dbr_assessed = dbr_sought_reports - dbr_notretrieved_reports

  new_studies <= dbr_assessed - dbr_excluded

Studies vs reports: one study can produce several reports. The included box
reports both counts, and they are usually different. Conflating them is the
most common PRISMA diagram error.

previous = TRUE is only for review updates. Leaving it at the default on a
new review produces boxes for prior studies that do not exist.

Data Extraction

Extract in duplicate

Two extractors, independently, into the same template, then reconcile. Single extraction has a documented error rate high enough to change pooled estimates.

What to capture

Study level
  citation, PROSPERO/trial registration ID, country, funding source,
  conflicts of interest, design, follow-up duration

Population
  n randomized, n analysed, age, sex, disease stage, line of therapy,
  key prognostic factors

Intervention and comparator
  agent, dose, schedule, duration, co-interventions

Outcomes, per outcome and per timepoint
  definition as reported, timepoint, n analysed
  binary        events and total, per arm
  continuous    mean, SD, n, per arm
  time-to-event HR with CI, or the numbers needed to derive one

Always record what was NOT reported. "Not reported" and "zero" are different
and are handled differently downstream.

Deriving what is missing

Trials frequently report the wrong summary statistic. Convert rather than dropping the study, and record every conversion.

Median and IQR -> mean and SD        Wan et al. 2014, Luo et al. 2018
SE -> SD                            SD = SE * sqrt(n)
95% CI -> SD                        SD = sqrt(n) * (upper - lower) / 3.92
                                    3.92 = 2 * 1.96; use the t quantile if n < 60
p value -> SE                       back-calculate from the test statistic

Every derived value is an assumption. Flag them and test them in a
sensitivity analysis that excludes derived data.

Risk of Bias

Choosing the tool

What design are you assessing?

  Randomized trial
    -> RoB 2. Five domains, signalling questions, per-outcome not per-study.
       Assess each outcome separately; a trial can be low risk for mortality
       and high risk for a subjective outcome.

  Non-randomized study of an INTERVENTION
    -> ROBINS-I (2016). Seven domains, judged against a target trial.
       ROBINS-I V2 was posted 2025-11-20 but is still a DRAFT and is
       "subject to change". Use the 2016 version for work you intend to
       publish, and state which version you used.

  Non-randomized study of an EXPOSURE
    -> ROBINS-E. Same logic as ROBINS-I, adapted for exposures.

  Diagnostic accuracy study
    -> QUADAS-2.

  Prognostic factor study
    -> QUIPS.

  Missing evidence in the synthesis itself
    -> ROB ME.

  Newcastle-Ottawa Scale
    -> Widely used for cohort and case-control studies and often demanded by
       journals, but it produces a numeric score, and Cochrane recommends
       against collapsing risk of bias into a score. If a journal requires
       NOS, report it alongside ROBINS-I rather than instead of it.

RoB 2 and ROBINS-I are domain-judgement tools, not checklists. Each domain resolves to Low / Some concerns (or Moderate) / High / Critical via an algorithm from the signalling questions. Do not average domains: the overall judgement is driven by the worst domain.

Visualization

library(robvis)   # v0.3.1

# One row per study. Columns: Study, D1..Dk, Overall, and optionally Weight.
# Judgement strings must match the tool's expected levels exactly.
rob_traffic_light(data = rob_data, tool = "ROB2", psize = 10)

rob_summary(data = rob_data, tool = "ROB2", overall = TRUE,
            weighted = TRUE)   # weight bars by study weight from the model

# Supported: "ROB2", "ROB2-Cluster", "ROBINS-I", "ROBINS-E",
#            "QUADAS-2", "QUIPS", "Generic"
rob_tools()   # confirms the list for the installed version

Use tool = "Generic" for anything else, including the Newcastle-Ottawa Scale, which robvis does not model directly.

The robvis CRAN package is maintained (0.3.1, June 2026), but riskofbias.info
states they are no longer able to support the robvis web app or the Excel tool
implementations. Use the R package rather than the hosted app.

Risk of bias feeds the synthesis. Studies at high risk are not silently dropped: they are either excluded in a pre-specified sensitivity analysis, or retained with the sensitivity analysis reported alongside.

Effect Measures

Choose the measure before extraction, because it determines what has to be extracted.

Binary outcome
  OR   odds ratio. Symmetric, works with case-control, but is misread as a
       risk ratio whenever the event is common (> ~10%).
  RR   risk ratio. Interpretable, preferred for cohort and trial data.
       Not estimable from case-control designs.
  RD   risk difference. Absolute scale, so it transports poorly across
       populations with different baseline risk. Usually more heterogeneous.
  PETO one-step OR. Only for rare events with balanced arms. Biased when
       arms are unbalanced or effects are large.

Continuous outcome
  MD   mean difference. Use when every study measured the SAME instrument
       on the same scale.
  SMD  standardised mean difference (Hedges' g in metafor). Use when studies
       used DIFFERENT instruments for the same construct.
  ROM  ratio of means. Use for ratio-scale outcomes where a proportional
       change is more natural than an absolute one.

Time-to-event
  HR   hazard ratio, pooled on the log scale.

Pick one and keep it. Switching measure after seeing the pooled result is a form of analytic flexibility that inflates false positives.

Computing Effect Sizes

escalc() computes the effect size yi and its sampling variance vi, which is what every model consumes.

library(metafor)   # v5.0-1
data(dat.bcg, package = "metadat")

# Binary: 2x2 table per study
dat <- escalc(measure = "RR", data = dat.bcg,
              ai = tpos, bi = tneg,     # events / non-events, treatment
              ci = cpos, di = cneg)     # events / non-events, control

# Continuous
# escalc(measure = "SMD", m1i =, sd1i =, n1i =, m2i =, sd2i =, n2i =, data = )

# Time-to-event: supply log(HR) and its standard error directly
# dat$yi <- log(dat$hr);  dat$vi <- ((log(dat$hr_upper) - log(dat$hr_lower)) / 3.92)^2

metafor 5.0 changed two escalc() defaults. Code written for 4.x runs without
error and returns DIFFERENT numbers.

  correct = TRUE is now the default for "ROM", "ROMC", "CVR" and "CVRC".
  The second-order Taylor bias correction is applied unless you pass
  correct = FALSE. Any ROM or CVR meta-analysis run on 4.x will not reproduce
  on 5.x at the default.

  The default `add` value changed to 0 for eight measures where bias
  corrections are now applied.

  pi.type was renamed predtype. The old name still works but is deprecated.

If you are reproducing a published analysis, pin the metafor version and say
which one you used.

Zero cells

# add = 1/2, to = "only0" is the default: add 0.5 only to studies with a zero cell
dat <- escalc(measure = "OR", ai = ai, bi = bi, ci = ci, di = di, data = dat)

# Double-zero studies contribute nothing and are dropped by default
# drop00 = TRUE removes them explicitly

The 0.5 continuity correction is a convenience, not a solution. It biases the
estimate toward the null and the bias grows as arms become unbalanced.

For rare events, prefer a method that does not need it:
  - Peto OR, when events are rare AND arms are roughly balanced
  - a beta-binomial or exact model via rma.glmm()

Never "fix" zero cells by deleting the studies. That is informative deletion.

Fitting the Model

The terminology trap

rma(yi, vi, data = dat, method = "EE")   # equal-effects
rma(yi, vi, data = dat, method = "FE")   # fixed-effects

"EE" and "FE" produce IDENTICAL numbers and mean different things.

  EE (equal-effects)  assumes one single true effect underlies every study.
                      Differences between studies are sampling error only.
                      This is what most people mean when they write
                      "fixed-effect meta-analysis".

  FE (fixed-effects)  makes no such assumption. Inference is conditional on
                      the studies actually included: it estimates the average
                      effect IN THIS SET, and does not generalize beyond it.

Older metafor used "FE" for what is now "EE". Papers saying "fixed effect"
almost always mean EE. State which model you fitted and what you claim from it.

Random-effects, the default

res <- rma(yi, vi, data = dat,
           method = "REML",   # tau^2 estimator; the default
           test = "knha")     # Knapp-Hartung; NOT the default, must be asked for
summary(res)

Two choices carry most of the weight.

tau^2 estimator: REML
  metafor's default and its author's recommendation, because REML gives
  approximately unbiased estimates of heterogeneity. DerSimonian-Laird is
  the historical default in older software and underestimates tau^2, which
  makes confidence intervals too narrow. Use REML unless reproducing an
  older analysis, and then say so.
  Available: "REML", "ML", "DL", "PM", "EB", "SJ", "HS", "HSk", "HE", "GENQ".

Knapp-Hartung: test = "knha"
  Default is test = "z", which uses a normal distribution and produces
  intervals that are too narrow when the number of studies is small.
  test = "knha" uses a t-distribution with k - p degrees of freedom.
  metafor's author calls it "highly recommended".

  Honest caveat: simulation work (IntHout 2014, Jackson 2017) shows coverage
  is slightly BELOW nominal when heterogeneity is low (I^2 < 30%) and study
  sizes are very uneven. It still beats DerSimonian-Laird across most of the
  parameter space. Report that you used it.

Heterogeneity

res                  # prints Q, its p-value, tau^2, I^2, H^2
confint(res)         # confidence interval for tau^2 and I^2 — report it
predict(res, digits = 3)   # pooled estimate with a PREDICTION interval

What each quantity actually tells you:

  Q       a test of whether heterogeneity exceeds sampling error. Badly
          underpowered with few studies, and trivially significant with many.
          A non-significant Q does NOT establish homogeneity.

  tau^2   the variance of true effects, on the analysis scale. The only one
          of these that is a magnitude rather than a proportion.

  I^2     the PERCENTAGE OF VARIABILITY due to heterogeneity rather than
          chance. It is NOT the amount of heterogeneity. I^2 rises as studies
          get larger even when tau^2 is unchanged, because sampling error
          shrinks. Two meta-analyses with identical tau^2 can have I^2 of 25%
          and 90%.

  The 25/50/75% "low/moderate/high" thresholds are explicitly described in the
  Cochrane Handbook as rough and context-dependent. Do not treat them as rules.

Prediction interval
  The confidence interval describes the AVERAGE effect. The prediction
  interval describes where the effect of a NEW study would fall. With
  substantial tau^2 the prediction interval routinely crosses the null while
  the confidence interval does not. Report both, or the review overstates
  what is known.

Subgroup Analysis

# Subgroups are a moderator, not separate meta-analyses
res_sub <- rma(yi, vi, mods = ~ factor(alloc), data = dat, test = "knha")
res_sub          # QM = omnibus test of the moderator; this is the test that matters

# Pooled estimate within each level, with a shared tau^2
predict(res_sub, newmods = rbind(c(0,0), c(1,0), c(0,1)))

The mistake that shows up in most published subgroup analyses:

  Running a separate meta-analysis in each subgroup and comparing whether one
  is significant and the other is not. That is not a comparison. A subgroup
  can be significant purely because it has more studies.

  The correct question is whether the SUBGROUP DIFFERENCE is non-zero, which
  is the QM test above.

Subgroup analyses are observational even in a review of randomized trials.
Studies were not randomized to subgroups, so a subgroup difference is a
hypothesis, not an effect. Pre-specify them, keep them few, and report how
many you ran.

Meta-Regression

res_mr <- rma(yi, vi, mods = ~ ablat + year, data = dat, test = "knha")
res_mr
# R^2 in the output = proportion of tau^2 explained by the moderators

regplot(res_mr, mod = "ablat", xlab = "Absolute latitude", las = 1)

Power rule of thumb: at least 10 studies per moderator, and that is a floor,
not a target. Meta-regression on 8 studies with 2 moderators is curve-fitting.

Aggregation bias: a study-level covariate is not a patient-level covariate.
A relationship between mean age and effect size across studies does not imply
the same relationship across patients. This is ecological inference, and it
is the single most over-claimed result in meta-regression.

Hazard Ratios from Published Curves

When a trial reports a Kaplan-Meier curve but no hazard ratio, the HR can be reconstructed.

Preferred order:

  1. HR and CI reported directly            use them
  2. Reconstruct from reported statistics   Parmar/Tierney methods, using
                                            O-E and variance, or the log-rank
                                            p-value with events per arm
  3. Digitize the KM curve                  Guyot algorithm reconstructs
                                            individual patient data from the
                                            curve plus numbers at risk

Digitizing is a last resort. It requires the numbers-at-risk table to be
printed; without it the reconstruction is unreliable. The R implementation
(IPDfromKM) was last released in 2020, so validate its output against any
reported median survival before pooling.

Whatever you use, record the method per study and run a sensitivity analysis
excluding reconstructed estimates.

Forest Plots

forest(res,
       slab = paste(dat$author, dat$year, sep = ", "),
       atransf = exp,               # display on the ratio scale, model fitted on log
       at = log(c(0.05, 0.25, 1, 4)),
       xlab = "Risk Ratio (log scale)",
       header = "Author(s) and Year",
       mlab = "")
addpoly(res, row = -1, atransf = exp, mlab = "RE Model (REML, KNHA)")

Fit ratio measures on the log scale and transform only for display. Pooling raw ratios is wrong: the sampling distribution is skewed and the variance formula assumes the log scale.

Publication Bias and Small-Study Effects

"Publication bias" is one explanation for funnel plot asymmetry. It is not the only one, and the tests cannot distinguish between them.

Why a funnel plot can be asymmetric, all producing the same picture:

  Publication bias      small negative studies never published
  True heterogeneity    small studies done in higher-risk populations with
                        genuinely larger effects
  Poorer methods        small studies less likely to be blinded or allocation
                        concealed, inflating their effects
  Outcome reporting     the significant outcome reported, the others not
  Chance                with 10 studies, asymmetry happens

The Cochrane Handbook's term for this family is SMALL-STUDY EFFECTS, which is
the honest label. Use it in the manuscript rather than asserting publication
bias you cannot demonstrate.

Funnel plot and Egger's test

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

regtest(res, model = "lm")   # Egger's regression test
ranktest(res)                # Begg's rank correlation, low power

Do not run Egger's test with fewer than 10 studies. The Cochrane Handbook is
explicit, and running it anyway is one of the most common errors in published
reviews: with k < 10 the test cannot distinguish asymmetry from chance.

A significant Egger's test does not establish publication bias. It establishes
asymmetry. The five causes above all produce it.

Agreement between Begg's test, Egger's test and trim-and-fill is empirically
only weak to moderate, so a "negative" result from one is not reassurance.

For ratio measures use the log scale, and prefer a funnel plot against the
standard error rather than the sample size.

Trim-and-fill, and why to be careful with it

tf <- trimfill(res)
tf              # imputed studies and the "adjusted" estimate
funnel(tf)

Cochrane Handbook, on trim-and-fill:

  it is "built on the strong assumption that there should be a symmetric
  funnel plot", it "does not take into account reasons for funnel plot
  asymmetry other than publication bias", and "'corrected' intervention
  effect estimates from this method should be interpreted with great caution".

Practical reading: use trim-and-fill as a SENSITIVITY ANALYSIS, to ask whether
the conclusion survives a pessimistic scenario. Do not report the filled
estimate as the result, and do not report the number of imputed studies as if
it were a count of suppressed trials. It is not.

Selection models, the principled alternative

# Model the publication process explicitly rather than assuming symmetry
sel <- selmodel(res, type = "step", steps = c(0.025, 0.5))
sel

library(metasens)   # v1.5-3
limitmeta(m)        # limit meta-analysis: estimate as study size -> infinity
copas(m)            # Copas selection model, sensitivity across selection strengths

Selection models state their assumption about how publication depends on the p-value, which makes the assumption arguable. Trim-and-fill hides its assumption inside a symmetry requirement. Prefer the former when you have enough studies to fit it.

Sensitivity and Influence Diagnostics

leave1out(res)                 # refit dropping each study in turn
inf <- influence(res); plot(inf)   # Cook's distance, DFFITS, hat, tau^2 delete-1
baujat(res)                    # heterogeneity contribution vs influence on the pooled effect
cumul(res, order = year)       # cumulative meta-analysis, in time order

What each answers:

  leave1out   does any single study drive the result? Report the range of
              pooled estimates, not just the full-data one.
  influence   which studies are outliers AND influential. Being an outlier
              is not enough; a small outlier changes nothing.
  baujat      separates "inflates heterogeneity" from "moves the estimate".
              The top-right quadrant is what to investigate.
  cumul       shows whether the effect stabilized or is still moving. An
              estimate that shifts with each new trial is not settled.

Pre-specify sensitivity analyses. Running leave-one-out and reporting only the
version that reaches significance is p-hacking with extra steps.

Network Meta-Analysis

NMA compares three or more treatments by combining direct and indirect evidence. It answers questions no single trial asked.

Transitivity, the assumption everything rests on

Indirect evidence for A vs C comes from A vs B and B vs C trials.

That is valid ONLY if the A-vs-B and B-vs-C trials are similar enough that a
patient in one could plausibly have been enrolled in the other. Formally:
effect modifiers must be distributed similarly across comparisons.

Check it BEFORE fitting anything, and check it with clinical data, not
statistics: tabulate age, disease severity, line of therapy, year, and dose
by comparison. If early trials used a lower dose of B than later ones, the
network is not transitive and no amount of modelling repairs it.

Statistical inconsistency tests have low power. A non-significant
inconsistency test is not evidence of transitivity.

Fitting

library(netmeta)   # v3.6-1
library(meta)      # pairwise() lives HERE, not in netmeta

# Long arm-based data -> contrast-based pairwise comparisons,
# with correlated multi-arm trials handled correctly
p <- pairwise(treat = treatment, event = events, n = total,
              studlab = study, data = dat, sm = "OR")

net <- netmeta(p, common = FALSE, random = TRUE,
               reference.group = "placebo",
               small.values = "desirable",   # do lower values mean benefit?
               prediction = TRUE)
summary(net)

Multi-arm trials are the trap. A three-arm trial contributes three pairwise
comparisons that are CORRELATED, because they share an arm. Feeding those in
as if independent double-counts patients and understates the standard error.

pairwise() constructs the correct correlation structure. Building the
comparison table by hand in a spreadsheet does not.

Network geometry

netgraph(net, plastic = FALSE, thickness = "number.of.studies",
         number.of.studies = TRUE)

Look at it before interpreting anything. A comparison supported by one small trial can sit between two well-studied nodes and carry the entire indirect estimate.

Inconsistency

netsplit(net)        # direct vs indirect per comparison, with a p-value
decomp.design(net)   # design-by-treatment interaction, the global test
netheat(net)         # which designs contribute the inconsistency

netsplit is local: it asks whether direct and indirect disagree for each
comparison. decomp.design is global.

Both are underpowered. Report them, and treat a disagreement as a reason to
re-examine transitivity rather than as a number to correct.

If direct and indirect genuinely conflict, the network estimate is a weighted
average of two incompatible things. Do not report it as if the conflict were
resolved.

Ranking, and its limits

netrank(net, small.values = "desirable")   # P-scores, frequentist analogue of SUCRA

P-scores and SUCRA are the most over-interpreted output in the whole field.

  A treatment can rank first on a single small trial. The ranking metric does
  not encode uncertainty in a way readers perceive.

  Rankings are unstable: adding one study can reorder the top three.

  "Ranked first" is not "better than second". Report the effect estimates with
  confidence intervals alongside any ranking, and never report a league table
  of ranks without them.

Report certainty too. A first-ranked treatment supported by very-low-certainty
evidence should be described that way.

Bayesian NMA

library(gemtc)     # v1.1-1, JAGS backend
# mtc.network() -> mtc.model() -> mtc.run(), then gelman.diag() for convergence

library(multinma)  # v0.9.1, Stan backend; supports individual patient data

Use Bayesian NMA when you need full posterior distributions, want to incorporate prior information on heterogeneity, or need to combine aggregate and individual patient data. Check convergence (R-hat, trace plots) before reading any estimate.

Certainty of Evidence

Pairwise meta-analysis -> GRADE
  Start high for RCTs, low for observational studies, then rate down for:
    risk of bias, inconsistency, indirectness, imprecision, publication bias
  and up (observational only) for:
    large effect, dose-response, plausible confounding working against the effect

  Produce a Summary of Findings table. GRADEpro is the standard tool.

Network meta-analysis -> CINeMA
  GRADE does not extend cleanly to networks, because each estimate mixes
  direct and indirect evidence in different proportions. CINeMA
  (Confidence In Network Meta-Analysis) handles this across six domains:
  within-study bias, reporting bias, indirectness, imprecision, heterogeneity,
  and incoherence.

  CINeMA is a WEB APPLICATION, not an R package. Do not look for it on CRAN.
  netmeta can export the data it needs.

Certainty is assessed per outcome and per comparison, not once for the review. A review with high certainty for the primary outcome and very low certainty for harms must say so.

Output Specification

Validation Checks

Search
  Every known-relevant "seed" study is retrieved by the final search.
  Build a seed set of 5-10 papers you already know qualify, and confirm the
  search f

…(truncated)