Publication bias and sensitivity in R

Practical · Wednesday 15:22–15:55 · FRB-CESAB training course

By the end of this practical, you will be able to

1. Draw and read a contour-enhanced funnel plot.
2. Test for small-study effects and time-lag bias in a multilevel model.
3. Run leave-one-study-out, influence and study-quality sensitivity analyses.
4. Build a robustness table for your paper.

Time	What	Section
15:22	Load data, fit the main model, funnel plot (code given)	1 · 2
15:27	Classic vs multilevel Egger, PET-PEESE	3
15:37	Sensitivity: leave-one-study-out, validity	4
15:47	Your turn — or the (optional) parts: time-lag, Cook, fail-safe N	6 · 5
15:55	Debrief with the whole group	7

Sections marked (optional) are for those who are ahead.

1 Setup

Open meta_analysis_course.Rproj (the course kit), then td_bias.R (exercises with blanks ___). Same data as in the previous practical.

library(dplyr)
library(metafor)
library(ggplot2)

# In the participant pack: data_file <- "intercropping_biodiversity.csv"
ic <- read.csv(data_file)

grp <- "Pests"        # <- in section 6 you will change this word only

dat <- ic |>
  filter(functional_group == grp, !sd_missing, !zero_mean) |>
  escalc(measure = "ROM", m1i = b_mean_t, sd1i = b_sd_t, n1i = b_n_t,
         m2i = b_mean_c, sd2i = b_sd_c, n2i = b_n_c, data = _) |>
  mutate(sei = sqrt(vi),                          # standard error
         se_n = sqrt(1 / b_n_t + 1 / b_n_c),     # precision based on sample size only
         inv_n = 1 / b_n_t + 1 / b_n_c,
         year_c = year - mean(year))              # centred publication year

pct <- function(x) round(100 * (exp(x) - 1), 1)

# the main model (from the previous practical)
m <- rma.mv(yi, vi, random = ~ 1 | study_id / es_id, data = dat,
            test = "t", dfs = "contain")
pct(c(m$b, m$ci.lb, m$ci.ub))

[1] -26.9 -39.9 -11.1

2 Contour-enhanced funnel plot

Exercise 1. Draw the funnel with significance contours centred on zero (fill in one number, then spend your time on the interpretation).

funnel(dat$yi, dat$vi, yaxis = "sei",
       level = c(90, 95, 99), shade = c("white", "gray75", "gray90"),
       refline = ___, legend = TRUE, xlab = "lnRR")

Solution

funnel(dat$yi, dat$vi, yaxis = "sei",
       level = c(90, 95, 99), shade = c("white", "gray75", "gray90"),
       refline = 0, legend = TRUE, xlab = "lnRR")

Interpret (2 min)

Is the funnel symmetric? On which side, and in which zone (white = non-significant), are the imprecise comparisons?

Answer

The bottom of the funnel (imprecise comparisons) is roughly symmetric around 0, with points in the white (non-significant) zone on both sides: no obvious gap where suppressed non-significant results would be. At the very top, 83 comparisons have SE < 0.02 — suspiciously precise (two studies; possibly SEs recorded as SDs). Tiny variances are always worth checking at extraction.

3 Testing small-study effects and time-lag bias

Exercise 2. First the classic Egger test, which treats every comparison as independent:

regtest(rma(yi, vi, data = dat))

Regression Test for Funnel Plot Asymmetry

Model:     mixed-effects meta-regression model
Predictor: standard error

Test for Funnel Plot Asymmetry: z =  3.5610, p = 0.0004
Limit Estimate (as sei -> 0):   b = -0.4245 (CI: -0.5092, -0.3398)

Now the multilevel version (Nakagawa et al. 2022): add a precision measure as a moderator of the main model.

The lnRR trap

The SE of a lnRR is computed from the means themselves, which creates an artificial correlation between effect and SE. For lnRR (and SMD), use a precision measure based on sample size only: se_n $= \sqrt{1/n_T + 1/n_C}$.

Predict first: will the small-study effect look stronger with sei or with se_n?

pet_se <- rma.mv(yi, vi, mods = ~ sei,  random = ~ 1 | study_id / es_id,
                 data = dat, test = "t", dfs = "contain")
pet    <- rma.mv(yi, vi, mods = ~ ___, random = ~ 1 | study_id / es_id,
                 data = dat, test = "t", dfs = "contain")
round(c(slope_p_SE = pet_se$pval[2], slope_p_n = pet$pval[2]), 3)
robust(pet, cluster = study_id, clubSandwich = TRUE)   # cluster-robust check

Solution

pet_se <- rma.mv(yi, vi, mods = ~ sei,  random = ~ 1 | study_id / es_id,
                 data = dat, test = "t", dfs = "contain")
pet    <- rma.mv(yi, vi, mods = ~ se_n, random = ~ 1 | study_id / es_id,
                 data = dat, test = "t", dfs = "contain")
round(c(slope_p_SE = pet_se$pval[2], slope_p_n = pet$pval[2]), 3)

slope_p_SE  slope_p_n 
     0.018      0.103 

robust(pet, cluster = study_id, clubSandwich = TRUE)

Multivariate Meta-Analysis Model (k = 382; method: REML)

Variance Components:

            estim    sqrt  nlvls  fixed          factor 
sigma^2.1  0.2078  0.4559     28     no        study_id 
sigma^2.2  0.1893  0.4351    382     no  study_id/es_id 

Test for Residual Heterogeneity:
QE(df = 380) = 100480.3902, p-val < .0001

Number of estimates:   382
Number of clusters:    28
Estimates per cluster: 1-80 (mean: 13.64, median: 8)

Test of Moderators (coefficient 2):¹
F(df1 = 1, df2 = 13.38) = 1.9876, p-val = 0.1814

Model Results:

         estimate      se¹     tval¹     df¹    pval¹    ci.lb¹    ci.ub¹     
intrcpt   -0.5648  0.1699   -3.3246   16.78   0.0041   -0.9236   -0.2060   ** 
se_n       0.6900  0.4894    1.4098   13.38   0.1814   -0.3643    1.7443      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1) results based on cluster-robust inference (var-cov estimator: CR2,
   approx t/F-tests and confidence intervals, df: Satterthwaite approx)

With the SE, the slope looks “significant”; with the sample-size measure (and even more with cluster-robust inference) it is not: no robust evidence of a small-study effect. The SE-based result was largely the lnRR artefact.

PET-PEESE: would the conclusion change if there were bias? The intercept of pet estimates the effect of an infinitely large study (PET); PEESE uses inv_n instead. Rule (Stanley and Doucouliagos 2014): if the PET intercept is significant, report PEESE.

peese <- rma.mv(yi, vi, mods = ~ inv_n, random = ~ 1 | study_id / es_id,
                data = dat, test = "t", dfs = "contain")
rbind(main  = pct(c(m$b, m$ci.lb, m$ci.ub)),
      PET   = pct(c(coef(pet)[1],   pet$ci.lb[1],   pet$ci.ub[1])),
      PEESE = pct(c(coef(peese)[1], peese$ci.lb[1], peese$ci.ub[1])))

      intrcpt            
main    -26.9 -39.9 -11.1
PET     -43.2 -60.8 -17.6
PEESE   -35.6 -50.5 -16.3

The adjusted estimates are, if anything, stronger reductions: the conclusion does not rest on publication bias. With strong heterogeneity, treat them as sensitivity checks, not as “the true effect”.

Exercise 3 (optional). Time-lag bias: do effects change with publication year?

tl <- rma.mv(yi, vi, mods = ~ ___, random = ~ 1 | study_id / es_id,
             data = dat, test = "t", dfs = "contain")
tl

Solution

tl <- rma.mv(yi, vi, mods = ~ year_c, random = ~ 1 | study_id / es_id,
             data = dat, test = "t", dfs = "contain")
tl

Multivariate Meta-Analysis Model (k = 382; method: REML)

Variance Components:

            estim    sqrt  nlvls  fixed          factor 
sigma^2.1  0.2168  0.4656     28     no        study_id 
sigma^2.2  0.1908  0.4368    382     no  study_id/es_id 

Test for Residual Heterogeneity:
QE(df = 380) = 96736.6562, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 26) = 0.0257, p-val = 0.8740

Model Results:

         estimate      se     tval  df    pval    ci.lb    ci.ub     
intrcpt   -0.3131  0.0983  -3.1839  26  0.0037  -0.5153  -0.1110  ** 
year_c    -0.0029  0.0181  -0.1602  26  0.8740  -0.0402   0.0344     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Slope ≈ 0 and far from significant: no evidence of a decline effect — but with 28 studies the test has little power.

4 Sensitivity analyses

Exercise 4. Leave-one-study-out. The loop is given: complete it so that each refit drops one whole study. (≈ 1 min to run.)

ids <- unique(dat$study_id)
loo <- sapply(ids, function(i) {
  f <- rma.mv(yi, vi, random = ~ 1 | study_id / es_id,
              data = dat[dat$study_id ___ i, ], test = "t", dfs = "contain")
  c(est = f$b[1], ci.lb = f$ci.lb, ci.ub = f$ci.ub)
})
pct(range(loo["est", ]))     # range of the pooled effect
any(loo["ci.ub", ] > 0)      # does any CI cross zero?

Solution

ids <- unique(dat$study_id)
loo <- sapply(ids, function(i) {
  f <- rma.mv(yi, vi, random = ~ 1 | study_id / es_id,
              data = dat[dat$study_id != i, ], test = "t", dfs = "contain")
  c(est = f$b[1], ci.lb = f$ci.lb, ci.ub = f$ci.ub)
})
pct(range(loo["est", ]))

[1] -29.3 -22.9

any(loo["ci.ub", ] > 0)

[1] FALSE

We leave out studies, not rows: rows of a study are not independent.

Exercise 5 (optional). Influence (Cook’s distance, by study) (Viechtbauer and Cheung 2010):

cd <- cooks.distance(m, cluster = dat$___)
sort(round(cd, 2), decreasing = TRUE)[1:5]

Solution

cd <- cooks.distance(m, cluster = dat$study_id)
sort(round(cd, 2), decreasing = TRUE)[1:5]

1323  742  522  763  645 
0.32 0.14 0.11 0.11 0.07 

No study stands out from the others (the largest distance is small and close to the next ones): no single study drives the result.

Exercise 6. Study quality. How many comparisons could not be appraised (NA)? Then refit the model with only the comparisons that meet all validity criteria (validity_biodiv == 1, coded by the authors of the database), and test validity as a moderator (comparing two CIs by eye is not a test).

m_val <- rma.mv(yi, vi, random = ~ 1 | study_id / es_id,
                data = filter(dat, validity_biodiv == ___), test = "t", dfs = "contain")
pct(c(m_val$b, m_val$ci.lb, m_val$ci.ub))

Solution

m_val <- rma.mv(yi, vi, random = ~ 1 | study_id / es_id,
                data = filter(dat, validity_biodiv == 1), test = "t", dfs = "contain")
pct(c(m_val$b, m_val$ci.lb, m_val$ci.ub))

[1] -27.4 -39.1 -13.4

m_val$s.nlevels[1]                  # number of studies left

[1] 13

sum(is.na(dat$validity_biodiv))     # comparisons that could not be appraised

[1] 130

m_vmod <- rma.mv(yi, vi, mods = ~ factor(validity_biodiv), random = ~ 1 | study_id / es_id,
                 data = dat, test = "t", dfs = "contain")
m_vmod$pval[2]                      # does validity change the effect?

[1] 0.7367947

Same estimate, validity does not moderate the effect for pests. Note that a third of the comparisons could not be appraised: report it.

5 (optional) Fail-safe N: why not

fsn(yi, vi, data = dat)

Fail-safe N Calculation Using the Rosenthal Approach

Observed Significance Level: <.0001
Target Significance Level:   0.05

Fail-safe N: 4098124

Discuss

Millions of “missing” studies would be needed to cancel the effect. Does that make you confident? What does this number assume about the missing studies, about heterogeneity and about dependence?

Answer

It assumes missing studies have exactly zero effect, ignores heterogeneity and the dependence between comparisons, and says nothing about the size of any bias. It is almost always “reassuring”, hence uninformative: do not report it (Nakagawa et al. 2022).

(optional) Trim-and-fill, as a sensitivity analysis only

Trim-and-fill cannot handle dependent effect sizes; apply it to one aggregated effect per study (here assuming ρ = 0.5 within studies, and a common effect within each study — a strong simplification):

agg <- aggregate(dat, cluster = study_id, rho = 0.5)
tf  <- trimfill(rma(yi, vi, data = agg))
tf

Estimated number of missing studies on the right side: 0 (SE = 2.6471)

Random-Effects Model (k = 28; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 0.7345 (SE = 0.2107)
tau (square root of estimated tau^2 value):      0.8570
I^2 (total heterogeneity / total variability):   99.98%
H^2 (total variability / sampling variability):  5517.89

Test for Heterogeneity:
Q(df = 27) = 20289.3510, p-val < .0001

Model Results:

estimate      se     zval    pval    ci.lb    ci.ub     
 -0.4634  0.1665  -2.7826  0.0054  -0.7898  -0.1370  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here no study is imputed. Report it, if at all, as a sensitivity check — never as “the corrected effect”.

6 Your turn: natural enemies

Go back to section 1, set grp <- "Natural enemies", and rerun sections 1–3 only (skip the slow leave-one-out). Is there a small-study effect with se_n? Would bias explain the (small, non-significant) increase? Then run Exercise 6: does validity matter here?

Solution

ne <- ic |>
  filter(functional_group == "Natural enemies", !sd_missing, !zero_mean) |>
  escalc(measure = "ROM", m1i = b_mean_t, sd1i = b_sd_t, n1i = b_n_t,
         m2i = b_mean_c, sd2i = b_sd_c, n2i = b_n_c, data = _) |>
  mutate(se_n = sqrt(1 / b_n_t + 1 / b_n_c))
f_ne <- function(mods) rma.mv(yi, vi, mods = mods, random = ~ 1 | study_id / es_id,
                              data = ne, test = "t", dfs = "contain")
m_ne   <- f_ne(~ 1)
pet_ne <- f_ne(~ se_n)
val_ne <- f_ne(~ factor(validity_biodiv))
rbind(main = pct(c(m_ne$b, m_ne$ci.lb, m_ne$ci.ub)),
      PET  = pct(c(coef(pet_ne)[1], pet_ne$ci.lb[1], pet_ne$ci.ub[1])))

     intrcpt           
main     8.3  -9.1 29.0
PET    -10.4 -35.0 23.3

c(slope_p = pet_ne$pval[2], validity_p = val_ne$pval[2])

   slope_p validity_p 
0.13048793 0.01407015 

coef(val_ne)

                 intrcpt factor(validity_biodiv)1 
             -0.07274506               0.30890264 

The increase of natural enemies is small and not significant; the small-study slope is not significant with se_n, and the PET estimate is close to zero: the evidence is inconclusive, whatever the bias. But validity matters: comparisons meeting all validity criteria show a clearly larger increase of natural enemies (about +27%) than the others (about −7%) — poorly designed comparisons may be masking a real benefit. That is a moderator — the topic of the next session.

7 Debrief: your robustness table

Fill it in for pests — this is what reviewers expect:

Analysis	Estimate [95% CI]	Conclusion changes?
Main multilevel model
Only valid comparisons
Leave-one-study-out (range)
Small-study effect (se_n) / PEESE
Time-lag (slope)

1. Bias is best prevented by the search, not corrected by statistics.
2. Asymmetry ≠ publication bias: read the direction; use contour-enhanced funnels.
3. With several effect sizes per study, use multilevel tests — for lnRR with a sample-size-based precision; drop fail-safe N; trim-and-fill only as sensitivity.
4. Report all sensitivity analyses in a robustness table.

Next session (16:15): explaining heterogeneity — meta-regression. (The multilevel Egger test you just ran was already a meta-regression!)

Going further in the online notebook

• Publication bias and sensitivity
• Data extraction and critical appraisal
• Reading a meta-analysis critically

7.1 References

Nakagawa, S., M. Lagisz, M. D. Jennions, J. Koricheva, D. W. A. Noble, T. H. Parker, A. Sánchez-Tójar, Y. Yang, and R. E. O’Dea. 2022. Methods for testing publication bias in ecological and evolutionary meta-analyses. Methods in Ecology and Evolution 13:4–21.

Stanley, T. D., and H. Doucouliagos. 2014. Meta-regression approximations to reduce publication selection bias. Research Synthesis Methods 5:60–78.

Viechtbauer, W., and M. W.-L. Cheung. 2010. Outlier and influence diagnostics for meta-analysis. Research Synthesis Methods 1:112–125.