Work in pairs. Only two blanks (___) to fill: the rest of the script is given, and the work is to read and discuss the outputs. The questions are in the script (# Q:). | 16:55 | Debrief (until 17:00) | 5 |
1 Setup
Open meta_analysis_course.Rproj (the course kit), then td_advanced.R (exercises with blanks ___).
library(dplyr)library(metafor)library(ggplot2)# In the participant pack: data_file <- "intercropping_biodiversity.csv"ic <-read.csv(data_file)dat <- ic |>filter(functional_group %in%c("Pests", "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 = _)pct <-function(x) round(100* (exp(x) -1), 1)fit <-function(d, mods =~1) # the multilevel model of this afternoonrma.mv(yi, vi, mods = mods, random =~1| study_id / es_id,data = d, test ="t", dfs ="contain")table(dat$functional_group)
Natural enemies Pests
195 382
2 Meta-regression
Predict first: will intercropping have the same effect on pests and on their natural enemies?
Exercise 1. Add functional_group as a moderator.
m0 <-fit(dat) # no moderatorm1 <-fit(dat, mods =~ ___) # with moderatorm1 # find: Test of Moderators (F, p)robust(m1, cluster = study_id, clubSandwich =TRUE) # few studies: check with CR2
Multivariate Meta-Analysis Model (k = 577; method: REML)
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.1211 0.3480 34 no study_id
sigma^2.2 0.1720 0.4147 577 no study_id/es_id
Test for Residual Heterogeneity:
QE(df = 575) = 160501.6025, p-val < .0001
Number of estimates: 577
Number of clusters: 34
Estimates per cluster: 1-80 (mean: 16.97, median: 11.5)
Test of Moderators (coefficient 2):¹
F(df1 = 1, df2 = 10.04) = 20.3527, p-val = 0.0011
Model Results:
estimate se¹ tval¹ df¹ pval¹ ci.lb¹
intrcpt 0.1269 0.0893 1.4218 25.19 0.1674 -0.0569
functional_groupPests -0.4348 0.0964 -4.5114 10.04 0.0011 -0.6495
ci.ub¹
intrcpt 0.3108
functional_groupPests -0.2202 **
---
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)
Test of Moderators: F is large, p < 0.001 → the effect differs between groups.
With cluster-robust (CR2) tests the denominator df drop to about 10 (only a few studies carry the contrast), but p stays ≈ 0.001: the difference survives.
Exercise 2 (code given — run and read). One estimate per group (intercept removed with 0 +) and the share of heterogeneity explained. How large is the pseudo-R²? At which level does the moderator explain something?
Pests decrease, natural enemies tend to increase — consistent with the enemies hypothesis. But the moderator explains only about 10% of the heterogeneity overall, and almost all of it within studies (~15%) rather than between studies (~2%). Pseudo-R² is rough (it can even be negative) — both models must be fitted on the same rows, as here.
3 Within-study check and study quality
Exercise 3. Split the moderator into a within-study part (does a study that measured both groups see a difference?) and a between-study part (do studies of pests differ from studies of enemies?). Only a few studies measured both groups.
dat <- dat |>mutate(pest =as.numeric(functional_group =="Pests")) |>group_by(study_id) |>mutate(pest_m =mean(pest), # share of pest rows in the study (between)pest_w = pest - ___) |># deviation from the study mean (within)ungroup()sum(tapply(dat$pest, dat$study_id, function(x) length(unique(x))) ==2) # studies with bothm_wb <-fit(dat, mods =~ pest_w + pest_m)pct(coef(m_wb)[-1])m_wb
Solution
dat <- dat |>mutate(pest =as.numeric(functional_group =="Pests")) |>group_by(study_id) |>mutate(pest_m =mean(pest), pest_w = pest - pest_m) |>ungroup()sum(tapply(dat$pest, dat$study_id, function(x) length(unique(x))) ==2)
Within studies (same field, year, crop): pests about −36% relative to natural enemies, clearly different from 0.
Between studies: about −28%, and uncertain (p ≈ 0.09).
Same direction and similar size: the difference is not just an artefact of which studies measured which group. The within-study estimate is the one that controls for study context — it is the stronger evidence.
(optional)Exercise 4. For natural enemies only, is study quality (validity_biodiv) a moderator? How many comparisons could not be appraised?
Solution
ne <-filter(dat, functional_group =="Natural enemies")sum(is.na(ne$validity_biodiv))
Comparisons meeting all validity criteria show a clear increase of natural enemies; the others show none. But the “not valid” level rests on only 6 studies (93 rows) and the “valid” one on 12 studies: a between-study association with few clusters. Poor designs may mask a real benefit — or valid studies differ in other ways (crop, region, taxa). And 1 comparison in 5 could not be appraised at all.
4 Missing SDs
Back to pests, this time keeping the comparisons without a usable SD (but with non-zero means).
pz <- ic |>filter(functional_group =="Pests", !zero_mean)c(total =nrow(pz), sd_missing =sum(pz$sd_missing))
total sd_missing
440 58
Exercise 5. Check the CVs first (code given — run and read). The coefficient of variation (SD / mean) of each study with SDs. Which studies look implausible?
Two studies (722 and 729) have a mean CV above 3 — study 729 with n = 1,200 per group. In the database, these SDs were rebuilt as SE × √n with n = the number of subsamples, not of independent plots (pseudo-replication): the “SD” is inflated and so is n. Go back to the paper, and recompute from it — never just patch the number. Here we set these two studies aside, both from the CV pool and from the analysis.
Exercise 6. Impute and compare (code given — run it and read the output).
suspect <-c(722, 729) # studies flagged in Exercise 5# typical CV = mean CV weighted by n, per group, then squared# (Nakagawa et al. 2023; what metafor does with vtype = "AV")weighted.mean(cvs$cv_t, cvs$b_n_t)^2# without the screen: absurd
[1] 28.09928
ok <-filter(cvs, !study_id %in% suspect)cv2t <-weighted.mean(ok$cv_t, ok$b_n_t)^2cv2c <-weighted.mean(ok$cv_c, ok$b_n_c)^2c(cv2t, cv2c) # with the screen: plausible
Does the conclusion change? What would you write in the methods and results sections?
Answer
Nearly identical estimates (about −28% vs −30%): the conclusion is robust to the missing SDs. But look where the imputed rows come from: 45 of 53 from a single study (525). The first step is to go back to that paper. Without the CV screen, the squared typical CV is about 30 instead of about 0.2, and imputation would have given these rows almost no weight. Report the number of comparisons with imputed variances, the method (Nakagawa et al. 2023), the CV screening (and the studies set aside), and both results.
5 Debrief
Meta-regression tests whether a moderator explains heterogeneity: report the omnibus test, coefficients and residual heterogeneity (pseudo-R²).
Between-study moderator effects are associations; within-study contrasts are stronger.
Missing SDs: check CVs by study, recompute from the paper when you can, impute transparently, compare with complete cases.
Tomorrow (group projects): the biodiversity–yield pairs (biodiversity_yield_pairs.csv) are a good playground for a multivariate question.