Sections marked (optional) are for those who are ahead: skip them if you are short of time.
1 Setup
Open meta_analysis_course.Rproj (the course kit), then td_models.R (it contains every exercise with blanks ___). Same data and same effect sizes as this morning: pest abundance, intercropping vs monoculture. As announced this morning, comparisons with a zero mean or no usable SD are set aside (!sd_missing, !zero_mean) — in a real review, report how many and test the conclusion with imputed SDs.
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 onlysub <- ic |>filter(functional_group == grp, !sd_missing, !zero_mean)dat <-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 = sub)c(comparisons =nrow(dat), studies =n_distinct(dat$study_id))
comparisons studies
382 28
A small helper to turn lnRR into a % change:
pct <-function(x) round(100* (exp(x) -1), 1)
2 Common vs random effects
Exercise 1. Fit a common-effect model (method = "EE") and a random-effects model (default estimator REML). test = "knha" (Knapp–Hartung) widens the CI slightly to account for the uncertainty in τ² — recommended. Compare the estimates and CIs in %.
m_ee <-rma(yi, vi, data = dat, method ="___")m_re <-rma(yi, vi, data = dat, test ="knha")pct(c(m_ee$b, m_ee$ci.lb, m_ee$ci.ub))pct(c(m_re$b, ___, ___))
Solution
m_ee <-rma(yi, vi, data = dat, method ="EE")m_re <-rma(yi, vi, data = dat, test ="knha")pct(c(m_ee$b, m_ee$ci.lb, m_ee$ci.ub))
[1] -17.2 -17.3 -17.0
pct(c(m_re$b, m_re$ci.lb, m_re$ci.ub))
[1] -27.7 -32.7 -22.3
Exercise 2. Which study dominates? Compute the share of the total weight taken by each study in both models.
In the common-effect model, two studies carry almost all the weight (study 525, Zhou et al. 2013, and study 645, Maluleke et al. 2005). In the random-effects model the weights are more equal — but study 645, which has 80 comparisons, now counts for about a quarter of the result because it has many rows. That is the non-independence problem of section 4.
Look at dat[dat$study_id == 525, c("b_mean_t", "b_mean_c", "b_n_t", "vi")]: n = 100 subsamples and the same treatment means reused against several controls give tiny variances. Always inspect your smallest variances — they may reveal pseudo-replication at extraction.
3 Heterogeneity and prediction interval
m_re
Random-Effects Model (k = 382; tau^2 estimator: REML)
tau^2 (estimated amount of total heterogeneity): 0.3386 (SE = 0.0296)
tau (square root of estimated tau^2 value): 0.5819
I^2 (total heterogeneity / total variability): 99.92%
H^2 (total variability / sampling variability): 1227.22
Test for Heterogeneity:
Q(df = 381) = 159242.4776, p-val < .0001
Model Results:
estimate se tval df pval ci.lb ci.ub
-0.3244 0.0364 -8.9153 381 <.0001 -0.3960 -0.2529 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Read the output (2 min)
Find τ², I² and Q in the output above. What does I² = 99.9% mean — and not mean?
Answer
τ² ≈ 0.34 (τ ≈ 0.58 on the lnRR scale): true effects vary a lot between comparisons. I² ≈ 99.9% means that almost all the observed variance is real heterogeneity rather than sampling error — because sampling errors are small here, not necessarily because effects are extremely different. I² is relative; τ and the prediction interval tell you how much effects vary.
Exercise 3. Use predict() to get the mean, its CI and the prediction interval, back-transformed to %.
mean ci.lb ci.ub pi.lb pi.ub
-27.7 -32.7 -22.3 -77.0 127.5
In a typical field, pests decrease with intercropping, and the CI excludes 0. But the prediction interval spans from a strong decrease to a strong increase: in a new field, the effect could go either way. (exp(mean) − 1 is the change in a typical — median — field; the average % change, exp(μ + τ²/2) − 1, is smaller here.)
4 Multilevel model and robust inference
Exercise 4. Fit a three-level model: comparisons (es_id) nested in studies (study_id). Use test = "t", dfs = "contain" so that the degrees of freedom follow the number of studies.
m_ml <-rma.mv(yi, vi, random =~1| ___ / ___, data = dat, test ="t", dfs ="contain")m_mlpct(c(m_ml$b, m_ml$ci.lb, m_ml$ci.ub))
Solution
m_ml <-rma.mv(yi, vi, random =~1| study_id / es_id, data = dat, test ="t", dfs ="contain")m_ml
Multivariate Meta-Analysis Model (k = 382; method: REML)
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.2055 0.4533 28 no study_id
sigma^2.2 0.1909 0.4369 382 no study_id/es_id
Test for Heterogeneity:
Q(df = 381) = 159242.4776, p-val < .0001
Model Results:
estimate se tval df pval ci.lb ci.ub
-0.3138 0.0954 -3.2875 27 0.0028 -0.5096 -0.1179 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pct(c(m_ml$b, m_ml$ci.lb, m_ml$ci.ub))
[1] -26.9 -39.9 -11.1
The two variance components are sigma^2.1 (between studies) and sigma^2.2 (between comparisons within studies). The mean is similar to the random-effects model, but the CI is more than twice as wide: we now have 28 studies, not 382 independent observations.
(optional) How is the total variance split between the two levels? Multilevel I² (Nakagawa and Santos 2012): each variance component divided by the total, where the total also includes a “typical” sampling variance ṽ (a weighted average of the vi):
Code given — just run it:
w <-1/ dat$vik <-length(w)vt <- (k -1) *sum(w) / (sum(w)^2-sum(w^2)) # "typical" sampling varianceround(100* m_ml$sigma2 / (sum(m_ml$sigma2) + vt), 1) # % between studies, % within studies
[1] 51.8 48.1
Exercise 5. Add cluster-robust inference (a safety net if the dependence structure is misspecified).
Multivariate Meta-Analysis Model (k = 382; method: REML)
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.2055 0.4533 28 no study_id
sigma^2.2 0.1909 0.4369 382 no study_id/es_id
Test for Heterogeneity:
Q(df = 381) = 159242.4776, p-val < .0001
Number of estimates: 382
Number of clusters: 28
Estimates per cluster: 1-80 (mean: 13.64, median: 8)
Model Results:
estimate se¹ tval¹ df¹ pval¹ ci.lb¹ ci.ub¹
-0.3138 0.0951 -3.2990 26.09 0.0028 -0.5092 -0.1183 **
---
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-test and confidence interval, df: Satterthwaite approx)
pct(c(m_rob$b, m_rob$ci.lb, m_rob$ci.ub))
[1] -26.9 -39.9 -11.2
Close to the multilevel result: good news, the model-based CI was not over-optimistic.
(optional) Correlated sampling errors: shared control plots
When several treatments are compared with the same control, their sampling errors are correlated. A pragmatic “working model” (called CHE + RVE, Pustejovsky and Tipton (2022)): assume a correlation (e.g. ρ = 0.5) between all rows of a study, build the variance–covariance matrix with vcalc(), and keep robust() on top. RVE needs enough studies (≳ 10–20 clusters; clubSandwich = TRUE adds the small-sample correction). The ρ assumption is strong; with good group identifiers, vcalc() can instead compute the covariance of shared controls exactly (arguments grp1, grp2, w1, w2).
V <-vcalc(vi, cluster = study_id, obs = es_id, rho =0.5, data = dat)m_v <-rma.mv(yi, V, random =~1| study_id / es_id, data = dat, test ="t", dfs ="contain")pct(c(m_v$b, m_v$ci.lb, m_v$ci.ub))
Try ρ = 0.2 and 0.8: if the conclusion does not change, it is robust to this assumption.
5 Show the result: an orchard plot
Exercise 6. The plotting code is given: run it and match each element to the output of predict() (pred, ci.lb/ci.ub, pi.lb/pi.ub). Then answer: what share of the individual comparisons falls outside the prediction interval? Is that what you expected?
pr_ml <-predict(m_ml)set.seed(1)ggplot(dat, aes(x = yi, y =1)) +geom_jitter(aes(size =1/sqrt(vi)), height =0.3, alpha =0.3, colour ="#56B8B9") +annotate("segment", x = pr_ml$pi.lb, xend = pr_ml$pi.ub, y =1, yend =1, linewidth =1) +annotate("segment", x = pr_ml$ci.lb, xend = pr_ml$ci.ub, y =1, yend =1, linewidth =4) +annotate("point", x = pr_ml$pred, y =1, size =6, shape =21, fill ="#BCCF00") +geom_vline(xintercept =0, linetype =2) +scale_size_continuous(range =c(0.5, 5), guide ="none") +scale_x_continuous(sec.axis =sec_axis(~100* (exp(.) -1), name ="% change",breaks =c(-95, -75, -50, 0, 100, 300))) +coord_cartesian(xlim =c(-4, 2)) +labs(x ="lnRR (pest abundance)", y =NULL,subtitle =paste0("k = ", nrow(dat), " comparisons, ", n_distinct(dat$study_id), " studies")) +theme_minimal(base_size =13) +theme(axis.text.y =element_blank(), panel.grid.major.y =element_blank())
The orchaRd package (Nakagawa et al. 2021) draws these plots directly (orchard_plot()), including for moderators; it is installed from GitHub, so we build it by hand here.
Answer
mean(dat$yi < pr_ml$pi.lb | dat$yi > pr_ml$pi.ub)
[1] 0.1151832
About 12% here — more than the naive 5%. The prediction interval is for true effects; observed values add their own sampling error on top, and a few studies are extreme. Worth a look: which studies are they?
6 Your turn: natural enemies
Predict first: will intercropping increase or decrease natural enemies?
Then go back to section 1, replace grp <- "Pests" by grp <- "Natural enemies", and rerun the script from the top. Compare the multilevel estimate and prediction interval with those for pests.
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 = _)m_ne <-rma.mv(yi, vi, random =~1| study_id / es_id, data = ne, test ="t", dfs ="contain") # same as rerunning with grp <- "Natural enemies"pr_ne <-predict(m_ne)c(k =nrow(ne), studies =n_distinct(ne$study_id))
mean ci.lb ci.ub pi.lb pi.ub
8.3 -9.1 29.0 -57.0 172.8
Natural enemies tend to be more abundant with intercropping, but the CI includes zero: the evidence is weaker than for pests. Comparing both groups properly requires a meta-regression (moderator = functional group): that is for the Advanced methods session.
7 Debrief and key messages
Weights follow precision; τ² tames very precise studies, the multilevel model stops studies with many rows from dominating.
Choose the model a priori: random effects, multilevel when studies give several effect sizes.
Report the mean with its CI, τ² and the prediction interval; add robust() as a safety net.
Show all the data: orchard plots.
Next session (15:00): are these 28 studies a representative sample — or is there publication bias?
