Code
# load libraries
library(MASS)
library(kableExtra)
library(metafor)
library(glmmTMB)Fitting meta-analysis with the glmmTMB R package
This webpage is supporting material to the paper “Meta-analysis with the glmmTMB R package” (DOI: XXX).
Below, we show how to fit a meta-analysis using the equalto covariance structure in glmmTMB R package (Brooks et al., 2017), which enables to specify the known sampling error variance-covariance (VCV) matrix.
Contact
For any questions or to raise an error, please contact Coralie Williams: coralie.williams@unsw.edu.au
Last updated: 04 April 2026
1 Introduction
To demonstrate how to fit a meta-analysis using the glmmTMB R package (Brooks et al., 2017), we will fit a multilevel meta-analytic model assuming multiple effect sizes per study. This model synthesises effect sizes across studies by accounting for variation at both the study and observation (effect size) levels, while also modelling correlated sampling errors within studies.
Given \(j\) effect sizes within \(i\) studies, we specify the model as:
\[ y_{ij} = b_0 + u_i + m_{ij} + e_{ij}, \]
where \(b_0\) is the overall mean effect, \(u_i \sim N(0, \sigma^2_{u}I)\) represents the among-study random effect, \(m_{ij} \sim N(0, \sigma^2_{m}I)\) captures the within-study random effect (at the observation-level), and \(e_{ij} \sim N(0, \text{VCV})\) are the sampling errors at the observation-level.
\(\text{VCV}\) represents a block-diagonal variance-covariance matrix including known (or calculated) sampling variances \(\sigma^2_{e}=v_i\) on the diagonal and accounts for correlations among sampling errors from the same study.
Load librairies
Check if the equalto covariance structure is available (it should be in position 15 of this output) :
Code
glmmTMB:::.valid_covstruct diag us cs ar1 ou exp gau mat toep rr
0 1 2 3 4 5 6 7 8 9
homdiag propto hetar1 homcs homtoep equalto
10 11 12 13 14 15 Simulate some data
Now, let’s simulate some data following the structure of the above multilevel meta-analysis.
Code
# set up indices
k.studies <- 10 # number of studies
k.per.study <- 3 # number of effect sizes per study
study <- rep(seq_len(k.studies), times=k.per.study) # study id
k <- length(study) # total number of effect sizes
id <- seq_len(k) # id of between study effect sizes
es.id <- unlist(lapply(k.per.study, seq_len)) # id of within study effect sizes
# fixed effect coeff true value (i.e. the overall mean)
b0 <- 0
# simulate sampling errors variances
set.seed(123)
vi <- rbeta(k, 2, 20)
# get VCV matrix assuming within-study correlation of rho=0.5
V <- matrix(0, nrow = k, ncol = k)
for (i in 2:k) {
for (j in 1:i) {
if (study[i] == study[j]) {
V[i,j] <- 0.5 * sqrt(vi[i] * vi[j])
}
}
}
V[upper.tri(V)] <- t(V)[upper.tri(V)] # fill in upper diagonal
diag(V) <- vi # fill in diagonal
# random effects
sigma2.u <- 0.2
sigma2.m <- 0.3
u <- rnorm(k.studies, 0, sqrt(sigma2.u))[study]
m <- rnorm(k, 0, sqrt(sigma2.m))
# simulate multivariate normal sampling errors with specific VCV
set.seed(123)
e <- mvrnorm(n = 1, mu = rep(0, length(vi)), Sigma = V)
#e <- rnorm(k, 0, sqrt(vi))[study] #without within-study correlation (diag VCV)
# compute y
y <- b0 + u + m + e
# combine into dataframe
dat <- data.frame(y = y,
vi = vi,
study = study,
id = id,
es.id = es.id,
obs = as.factor(id), #unique ID for propto and equalto (really important it is a factor variable)
g = 1 #group ID for propto and equalto (a constant here)
)Fit glmmTMB model
Before we fit the model we need to obtain the within-study error variance-covariance matrix (which accounts for within-study correlation) assuming a constant correlation of \(\rho=0.5\). The correlation between effect sizes from different studies is assumed to be zero.
We use the metafor vcalc function to obtain a block diagonal VCV matrix for the sampling errors:
Code
VCV <- metafor::vcalc(vi, cluster=study, obs=obs, data=dat, rho=0.5)Now we fit the glmmTMB model with equalto covariance structure.
We specify:
the study-level random effect \(u_i\) as
(1|study).the sampling errors \(e_{ij}\) are specified using the covariance structure
equalto(0 + obs|g,VCV)allowing to specify the VCV matrix of the sampling errors.the within-study random effect \(m_{ij}\) that captures the variability at the observation level is modelled as the residual error in this model.
Note
Important note: the equalto multivariate random effect variable (here it is called obs) must be a factor and the corresponding VCV matrix must have rows and columns matching the levels of the obs variable.
Code
# Fit glmmTMB model with equalto ----------------------------------
fit.equalto <- glmmTMB(y ~ 1 + (1|study) + equalto(0 + obs|g,VCV),
data=dat,
REML=T)
# > summary(fit.equalto)
#
# Family: gaussian ( identity )
# Formula: y ~ 1 + (1 | study) + equalto(0 + obs | g, VCV)
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 61.3 65.5 -27.6 55.3 27
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study (Intercept) 0.03440 0.1855
# g obs1 0.05545 0.2355
# obs2 0.07455 0.2730 0.00
# obs3 0.33583 0.5795 0.00 0.00
# obs4 0.09652 0.3107 0.00 0.00 0.00
# obs5 0.10144 0.3185 0.00 0.00 0.00 0.00
# obs6 0.13395 0.3660 0.00 0.00 0.00 0.00 0.00
# obs7 0.02728 0.1652 0.00 0.00 0.00 0.00 0.00 0.00
# obs8 0.04926 0.2220 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# obs9 0.06157 0.2481 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# obs10 0.12738 0.3569 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual 0.25008 0.5001
# Number of obs: 30, groups: study, 10; g, 1
#
# Dispersion estimate for gaussian family (sigma^2): 0.25
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 0.005619 0.131220 0.043 0.966Fit metafor model
Let’s compare with the metafor package rma.mv function that fits a multilevel model with random effects. Here the model is specified with two random effects: (1|study) corresponding to the among-study (\(u_i\)) and (1|id) corresponding to the within-study (\(m_{ij}\)) random effects. The sampling errors \(e_{ij}\) VCV matrix is specified using the V argument.
Code
# Fit metafor model --------------------------------------------
fit.metafor <- rma.mv(yi = y,
V = VCV,
random = list(~1 | study, ~1 | id),
control=list(REMLf=FALSE), #set this to obtain the exact same likelihood as glmmTMB
data = dat) #uses REML by default
# > summary(fit.metafor)
#
# Multivariate Meta-Analysis Model (k = 30; method: REML)
#
# logLik Deviance AIC BIC AICc
# -27.6455 55.2911 61.2911 65.3930 62.2511
#
# Variance Components:
#
# estim sqrt nlvls fixed factor
# sigma^2.1 0.0344 0.1855 10 no study
# sigma^2.2 0.2501 0.5001 30 no id
#
# Test for Heterogeneity:
# Q(df = 29) = 207.0443, p-val < .0001
#
# Model Results:
#
# estimate se zval pval ci.lb ci.ub
# 0.0056 0.1310 0.0429 0.9658 -0.2511 0.2624
#
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Then, we extract the model results for each model and compare outputs:
Code
### Compare output
equalto <- data.frame(model = "glmmTMB",
est = coef(summary(fit.equalto))$cond[,"Estimate"],
se = coef(summary(fit.equalto))$cond[,"Std. Error"],
zval = coef(summary(fit.equalto))$cond[,"z value"],
pval = coef(summary(fit.equalto))$cond[,"Pr(>|z|)"],
sigma.u = VarCorr(fit.equalto)$cond$study[,"(Intercept)"], ##among study variance estimate
sigma.m = sigma(fit.equalto)^2) ##within study variance estimate
metafor <- data.frame(model = "metafor",
est = coef(fit.metafor)[1],
se = se(fit.metafor)[1],
zval = fit.metafor$zval,
pval = fit.metafor$pval,
sigma.u = fit.metafor$sigma2[1], ## among study variance
sigma.m = fit.metafor$sigma2[2]) ## within study variance
output <- rbind(equalto, metafor)
kable(output)| model | est | se | zval | pval | sigma.u | sigma.m | |
|---|---|---|---|---|---|---|---|
| 1 | glmmTMB | 0.0056193 | 0.1312196 | 0.0428234 | 0.9658423 | 0.0344011 | 0.2500832 |
| intrcpt | metafor | 0.0056193 | 0.1310007 | 0.0428949 | 0.9657853 | 0.0344011 | 0.2500831 |
The overall mean estimate and random effect variance estimates are identical across the two packages ✅
Extra analyses
When the number of studies is small (<10) it is preferable to use the \(t\)-distribution and adjusted degrees of freedom to avoid Type I error rate of coefficient tests and low coverage rates of coefficient confidence intervals (read more about the recommendations here).
In metafor you can do this directly in the rma.mv function when fitting a multilevel/multivariate model with the arguments test="t" and dfs="contain" (actually if dfs="contain" is set this automatically sets test="t").
Code
# Fit metafor model --------------------------------------------
fit.metafor.t <- rma.mv(yi = y,
V = VCV,
random = list(~1 | study, ~1 | id),
test = "t",
dfs = "contain",
data = dat)
summary(fit.metafor.t)
Multivariate Meta-Analysis Model (k = 30; method: REML)
logLik Deviance AIC BIC AICc
-25.9449 51.8899 57.8899 61.9918 58.8499
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0344 0.1855 10 no study
sigma^2.2 0.2501 0.5001 30 no id
Test for Heterogeneity:
Q(df = 29) = 207.0443, p-val < .0001
Model Results:
estimate se tval df pval ci.lb ci.ub
0.0056 0.1310 0.0429 9 0.9667 -0.2907 0.3020
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1With glmmTMB you can derive the \(t\)-test statistic, adjusted degrees of freedom and corresponding P value the following way:
Code
# summary stats of fixed effects
fef <- as.data.frame(summary(fit.equalto)$coefficients$cond)
b0.est <- fef$Estimate
b0.se <- fef$`Std. Error`
# if the coefficient of interest varies at the level of corresponding to a random effect (i.e. here it is `study`) then the degree of freedom if the number of unique levels minus the number of estimated coefficients in the model
df.contain <- length(unique(dat$study)) - nrow(fef)
# otherwise it is the total nuber of estimates in the analysis minus the total number of model coefficients
#df.adj <- nrow(dat) - nrow(fef)
# get t-test statistic and p-value
t <- b0.est / b0.se
p <- 2 * pt(abs(t), df = df.contain, lower.tail = FALSE) # two-sided p value
# 95% CI
crit <- qt(0.975, df = df.contain)
ci <- c(b0.est - crit * b0.se, b0.est + crit * b0.se)
# glmmTMB t-test output
fef.t <- data.frame(Estimate = b0.est,
SE = b0.se,
tval = t,
df = df.contain,
pval = p,
ci.lb = ci[1],
ci.ub = ci[2])
fef.t Estimate SE tval df pval ci.lb ci.ub
1 0.005619274 0.1312196 0.04282345 9 0.9667773 -0.29122 0.30245852 Illustrative examples
To illustrate how to fit meta-analyses with glmmTMB, we re-analyse three published datasets from distinct research domains:
- Medicine — vaccine effectiveness.
- Ecology & evolution — plant diversity effects.
- Social sciences — school-based learning programs.
Across these case studies, we show model specification, output, and interpretation using glmmTMB.
2.1 Bivariate meta-analysis
Meta-analysis is widely used in medicine to combine evidence across clinical trials and observational studies. As many outcomes are events in medicine, common effect sizes measures used are odds ratios, risk ratios, risk differences, and their log transforms.
We will use the Colditz et al. (1994) dataset that is directly avaiable in metafor via the metadat R package. This published meta-analysis assessed the effectiveness of the Bacillus Calmette-Guerin (BCG) vaccine against tuberculosis across 13 studies. We will fit a bivariate meta-analysis that treats the two vaccination groups as separate outcomes, and estimates, for each group, the overall mean effect and among-study variance, as well as the correlation of study-level effects between groups (i.e., how strongly they co-vary).
Code
dat <- dat.colditz1994Reshape to long format so that each arm has its own row (“v” = vaccinated, “c” = unvaccinated/control), and setting control as the reference group.
Code
dat_long <- to.long(measure="OR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
levels(dat_long$group) <- c("c", "v")
dat_long$group <- relevel(dat_long$group, ref="c")We compute the logit-transformed event proportion (PLO) as the effect-size measure. And we also need to create a constant grouping dummy variable g for glmmTMB model specification (this needs to be a constant value for all rows).
Code
dat_long <- escalc(measure="PLO", xi=out1, mi=out2, data=dat_long)
dat_long$obs <- as.factor(1:nrow(dat_long))
dat_long$g <- 1 #dummy variable
head(dat_long)
trial author year tpos tneg cpos cneg ablat alloc study group out1
1 1 Aronson 1948 4 119 11 128 44 random 1 v 4
2 1 Aronson 1948 4 119 11 128 44 random 1 c 11
3 2 Ferguson & Simes 1949 6 300 29 274 55 random 2 v 6
4 2 Ferguson & Simes 1949 6 300 29 274 55 random 2 c 29
5 3 Rosenthal et al 1960 3 228 11 209 42 random 3 v 3
6 3 Rosenthal et al 1960 3 228 11 209 42 random 3 c 11
out2 yi vi obs g
1 119 -3.3928 0.2584 1 1
2 128 -2.4541 0.0987 2 1
3 300 -3.9120 0.1700 3 1
4 274 -2.2458 0.0381 4 1
5 228 -4.3307 0.3377 5 1
6 209 -2.9444 0.0957 6 1 Derive the variance-covariance matrix (VCV) with in diagonal the sampling error variances and zero on the off-diagonals (assume studies are independent):
Code
VCV <- diag(dat_long$vi)
VCV[1:5,1:5] #show first five measures of VCV matrix [,1] [,2] [,3] [,4] [,5]
[1,] 0.2584034 0.00000000 0.00 0.00000000 0.0000000
[2,] 0.0000000 0.09872159 0.00 0.00000000 0.0000000
[3,] 0.0000000 0.00000000 0.17 0.00000000 0.0000000
[4,] 0.0000000 0.00000000 0.00 0.03813239 0.0000000
[5,] 0.0000000 0.00000000 0.00 0.00000000 0.3377193To fit a bivariate meta-analysis with glmmTMB using equalto supply the sampling-error variance–covariance matrix (here it is called VCV), and fix the dispersion at zero (dispformula = ~ 0) since we don’t want to model any within-study variance effect.
Code
bv_fit <- glmmTMB(yi ~ group + (0 + group|study) + equalto(0 + obs|g,VCV),
dispformula = ~ 0,
REML = TRUE,
data=dat_long)
# > summary(bv_fit)
# Family: gaussian ( identity )
# Formula: yi ~ group + (0 + group | study) + equalto(0 + obs | g, VCV)
# Dispersion: ~0
# Data: dat_long
#
# AIC BIC logLik -2*log(L) df.resid
# 78.1 84.4 -34.1 68.1 21
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study groupc 2.617300 1.61781
# groupv 1.548603 1.24443 0.95
# g obs1 0.258403 0.50833
# obs2 0.098722 0.31420 0.00
# obs3 0.170000 0.41231 0.00 0.00
# obs4 0.038132 0.19528 0.00 0.00 0.00
# obs5 0.337719 0.58114 0.00 0.00 0.00 0.00
# obs6 0.095694 0.30934 0.00 0.00 0.00 0.00 0.00
# obs7 0.016203 0.12729 0.00 0.00 0.00 0.00 0.00 0.00
# obs8 0.004112 0.06412 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# obs9 0.030502 0.17465 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# obs10 0.021450 0.14646 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Number of obs: 26, groups: study, 13; g, 1
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -4.0960 0.4530 -9.042 < 2e-16 ***
# groupv -0.7414 0.1892 -3.918 8.94e-05 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#
#
emmeans::emmeans(bv_fit, ~ group) |>
summary(infer=TRUE, type="response") group emmean SE df asymp.LCL asymp.UCL z.ratio p.value
c -4.10 0.453 Inf -4.98 -3.21 -9.042 <0.0001
v -4.84 0.353 Inf -5.53 -4.14 -13.691 <0.0001
Confidence level used: 0.95 Code
# bivariate multilevel => compare with brms
bv_ml_fit <- glmmTMB(yi ~ group + (0 + group|study) + equalto(0 + obs|g,VCV),
dispformula = ~ (0 + group|study),
REML = TRUE,
data=dat_long)To obtain the variance component estimates and corresponding correlation \(\rho\) (i.e. rho) for each group we use the VarCorr function:
Code
vc <- VarCorr(bv_fit)
var_con <- vc$cond$study["groupc", "groupc"]
var_exp <- vc$cond$study["groupv", "groupv"]
rho.tmb <- vc$cond$study["groupc", "groupv"] / sqrt(var_con * var_exp)Then we can extract the model results as follows:
Code
bv_out <- data.frame(est.b0 = as.numeric(unlist(fixef(bv_fit))[[1]]),
se.b0 = as.numeric(sqrt(vcov(bv_fit)[[1]]))[1],
est.b1 = as.numeric(unlist(fixef(bv_fit))[[2]]),
se.b1 = as.numeric(sqrt(vcov(bv_fit)[[1]]))[4],
tau2_1 = var_con, #study-level heterogeneity for group=control
tau2_2 = var_exp, #study-level heterogeneity for group=exp
rho = rho.tmb)Warning in sqrt(vcov(bv_fit)[[1]]): NaNs produced
Warning in sqrt(vcov(bv_fit)[[1]]): NaNs producedCode
bv_out est.b0 se.b0 est.b1 se.b1 tau2_1 tau2_2 rho
1 -4.095985 0.4529815 -0.7413877 0.1892434 2.6173 1.548603 0.9450462This is equivalent to the following model in metafor:
Code
bv_fit_rma <- rma.mv(yi, vi, mods = ~ group,
random = ~ group | study,
struct = "UN",
data=dat_long)We can extract the model results as follows:
Code
bv_out_rma <- with(bv_fit_rma, data.frame(est.b0 = as.numeric(b[1]),
se.b0 = se[1],
est.b1 = as.numeric(b[2]),
se.b1 = se[2],
tau2_1 = tau2[1],
tau2_2 = tau2[2],
rho = rho))
bv_out_rma est.b0 se.b0 est.b1 se.b1 tau2_1 tau2_2 rho
1 -4.095985 0.4529369 -0.7413877 0.1880277 2.6173 1.548603 0.94504622.2 Phylogenetic meta-analysis
In meta-analysis of ecological and evolutionary studies, effect sizes are often drawn from multiple species. Because species share evolutionary history, these effect size estimates are not independent. To account for this phylogenetic relatedness, we can include a species-level random effect structured by a phylogenetic correlation (which is derived from a phylogenetic tree and assumed model of evolution). We can fit such a random effect in glmmTMB with the propto covariance structure.
We will use a published meta-analysis dataset of plant diversity effects on leaf traits (Felix et al., 2023) to showcase how to use glmmTMB to fit a phylogenetic meta-analysis.
Code
library(readr)
ef <- read_csv("felix_data/effect_sizes.csv")
options(na.action = "na.pass")Data set-up:
Code
# load libraries
suppressPackageStartupMessages({
library(rotl)
library(ape)
library(dplyr)
})
#remove empty rows in Nfixing
ef <- ef[which(ef$Nfixing!=is.na(ef$Nfixing)),]
# construct phylogentic tree matrix for use as random factor
species <- unique(ef$species) # list of unique species in meta-analysis
species <- as.character(species) # change to character object
taxa <- tnrs_match_names(species)
tree <- rotl::tol_induced_subtree(taxa$ott_id)
tree$tip.label <- strip_ott_ids(tree$tip.label, remove_underscores = TRUE) # change ids to the names from dataset
# calculate correlations between all species = cor
tree2 <- compute.brlen(tree)
phylo.cor <- vcv(tree2, cor = T)
# add phylo random factor
ef$phylo <- factor(ef$species, levels = rownames(phylo.cor))
# set up phylo.cor row and column names
phylo.cor <- phylo.cor[rownames(phylo.cor) %in% levels(ef$phylo),
colnames(phylo.cor) %in% levels(ef$phylo)]
# add g variable
ef$g <- 1
# derive effect size level variable => make sure it is a factor
ef$id <- as.factor(1:nrow(ef))
# set up sampling errors variance-covariance matrix
# assuming within-study a constant correlation of 0.5
V <- metafor::vcalc(vi, cluster=ACC, obs=id, data=ef, rho=0.5)
# make sure row matrix and id levels match
rownames(V) <- colnames(V) <- levels(ef$id)Let’s view the five rows of the species correlation matrix phylo.cor and the sampling error variance-covariance matrix V (which assumes no within-study or among effect size correlations):
Code
phylo.cor[1:5,1:5] Rhus chinensis Anacardium excelsum
Rhus chinensis 1.000000 0.990099
Anacardium excelsum 0.990099 1.000000
Choerospondias axillaris 0.980198 0.980198
Cedrela odorata 0.960396 0.960396
Swietenia macrophylla 0.960396 0.960396
Choerospondias axillaris Cedrela odorata
Rhus chinensis 0.980198 0.960396
Anacardium excelsum 0.980198 0.960396
Choerospondias axillaris 1.000000 0.960396
Cedrela odorata 0.960396 1.000000
Swietenia macrophylla 0.960396 0.990099
Swietenia macrophylla
Rhus chinensis 0.960396
Anacardium excelsum 0.960396
Choerospondias axillaris 0.960396
Cedrela odorata 0.990099
Swietenia macrophylla 1.000000Code
V[1:5,1:5] 1 2 3 4 5
1 0.13633355 0.06777343 0.06837032 0.06765727 0.06870997
2 0.06777343 0.13476472 0.06797580 0.06726687 0.06831349
3 0.06837032 0.06797580 0.13714894 0.06785929 0.06891514
4 0.06765727 0.06726687 0.06785929 0.13430315 0.06819640
5 0.06870997 0.06831349 0.06891514 0.06819640 0.13851498We fit a multilevel model using experimental site (), study ID (), individual effect ID () and plant species () were included as random factors to account for non-independence among effect sizes.
This multilevel phylogenetic meta-analysis uses two covariance structures in glmmTMB:
propto: to model shared evolutionary history across species assuming a proportional phylogenetic correlation matrixequalto: to incorporate known sampling-errors variance-covariance
Code
fit_phylo <- glmmTMB(yi ~ Nfixing-1 +
(1 | experiment) + (1 | ACC) +
(1 | species) + propto(0 + phylo|g, phylo.cor) +
equalto(0 + id|g, V),
data = ef,
REML=TRUE)
# Family: gaussian ( identity )
# Formula: yi ~ Nfixing - 1 + (1 | experiment) + (1 | ACC) + (1 | species) +
# propto(0 + phylo | g, phylo.cor) + equalto(0 + id | g, V)
# Data: ef
#
# AIC BIC logLik -2*log(L) df.resid
# 2856.4 2890.5 -1421.2 2842.4 948
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# experiment (Intercept) 2.044e-01 0.4521456
# ACC (Intercept) 6.535e-08 0.0002556
# species (Intercept) 1.063e-01 0.3260582
# g phyloRhus chinensis 2.958e-01 0.5438523
# phyloAnacardium excelsum 2.958e-01 0.5438523 0.99
# phyloChoerospondias axillaris 2.958e-01 0.5438523 0.98 0.98
# phyloCedrela odorata 2.958e-01 0.5438523 0.96 0.96 0.96
# phyloSwietenia macrophylla 2.958e-01 0.5438523 0.96 0.96 0.96 0.99
# phyloSapindus mukorossi 2.958e-01 0.5438523 0.89 0.89 0.89 0.89 0.89
# phyloKoelreuteria bipinnata 2.958e-01 0.5438523 0.89 0.89 0.89 0.89 0.89 0.99
# phyloAcer saccharum 2.958e-01 0.5438523 0.89 0.89 0.89 0.89 0.89 0.94 0.94
# phyloAcer pseudoplatanus 2.958e-01 0.5438523 0.89 0.89 0.89 0.89 0.89 0.94 0.94 0.99
# phyloAcer rubrum 2.958e-01 0.5438523 0.89 0.89 0.89 0.89 0.89 0.94 0.94 0.98 0.98
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# g.1 id1 1.363e-01 0.3692337
# id2 1.348e-01 0.3671031 0.50
# id3 1.371e-01 0.3703363 0.50 0.50
# id4 1.343e-01 0.3664739 0.50 0.50 0.50
# id5 1.385e-01 0.3721760 0.50 0.50 0.50 0.50
# id6 1.526e-01 0.3905846 0.50 0.50 0.50 0.50 0.50
# id7 5.509e-01 0.7422410 0.00 0.00 0.00 0.00 0.00 0.00
# id8 9.201e-01 0.9591954 0.00 0.00 0.00 0.00 0.00 0.00 0.50
# id9 1.176e+00 1.0846430 0.00 0.00 0.00 0.00 0.00 0.00 0.50 0.50
# id10 1.908e+00 1.3812190 0.00 0.00 0.00 0.00 0.00 0.00 0.50 0.50 0.50
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual 4.256e-01 0.6523670
# Number of obs: 955, groups: experiment, 39; ACC, 39; species, 102; g, 1
#
# Dispersion estimate for gaussian family (sigma^2): 0.426
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# NfixingN-fixing neighbours 0.10552 0.32404 0.326 0.745
# NfixingNo N-fixing neighbours 0.08583 0.31441 0.273 0.785We can extract and view the fixed effect estimates, showing the two groups for Nfixing (whether neighboring plants have nitrogen fixing or not):
Code
summary(fit_phylo)$coefficients$cond Estimate Std. Error z value Pr(>|z|)
NfixingN-fixing neighbours 0.10552475 0.3240432 0.3256502 0.7446890
NfixingNo N-fixing neighbours 0.08583235 0.3144103 0.2729948 0.7848572
attr(,"ddf")
[1] "asymptotic"The variance component estimates can be obtained:
Code
# get random-effects varcor parameters
re_est <- as.data.frame(confint(fit_phylo, parm="theta_"))
# summarise random-effects variance estimates
data.frame(sigma2.experiment = round(re_est$Estimate[1]^2,4),
sigma2.ACC = round(re_est$Estimate[2]^2,4),
sigma2.species = round(re_est$Estimate[3]^2,4),
sigma2.phylo = round(re_est$Estimate[4]^2,4),
sigma2.id=round(sigma(fit_phylo)^2,4)) sigma2.experiment sigma2.ACC sigma2.species sigma2.phylo sigma2.id
1 0.2044 0 0.1063 0.2958 0.4256This model is equivalent to the following model fitted with metafor:
Code
fit_phylo_rma <- rma.mv(yi, vi,
mods = ~ Nfixing-1,
random = list( ~ 1 | experiment, ~ 1 | ACC,
~1 | species, ~1 | phylo,
~1 | id),
R = list(phylo = phylo.cor),
data = ef,
method="REML") #default is REML
summary(fit_phylo_rma)
Multivariate Meta-Analysis Model (k = 955; method: REML)
logLik Deviance AIC BIC AICc
-1367.5406 2735.0813 2749.0813 2783.0986 2749.1998
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0588 0.2424 39 no experiment no
sigma^2.2 0.0000 0.0001 39 no ACC no
sigma^2.3 0.1651 0.4064 102 no species no
sigma^2.4 0.4570 0.6760 102 no phylo yes
sigma^2.5 0.0258 0.1608 955 no id no
Test for Residual Heterogeneity:
QE(df = 953) = 1433.1312, p-val < .0001
Test of Moderators (coefficients 1:2):
QM(df = 2) = 0.7668, p-val = 0.6815
Model Results:
estimate se zval pval ci.lb
NfixingN-fixing neighbours 0.2715 0.3561 0.7626 0.4457 -0.4263
NfixingNo N-fixing neighbours 0.2996 0.3498 0.8563 0.3918 -0.3861
ci.ub
NfixingN-fixing neighbours 0.9694
NfixingNo N-fixing neighbours 0.9853
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We can store time and compare with the model fitted with metafor:
Code
## Fit glmmTMB model and store runtime
time_glmmTMB <- system.time(
fit_phylo <- glmmTMB(
yi ~ Nfixing - 1 +
(1 | experiment) + (1 | ACC) +
(1 | species) +
propto(0 + phylo | g, phylo.cor) +
equalto(0 + id | g, V),
data = ef,
REML = TRUE
)
)[["elapsed"]]
## Fit rma.mv model and store runtime
time_rma <- system.time(
fit_phylo_rma <- rma.mv(
yi, vi,
mods = ~ Nfixing - 1,
random = list(
~ 1 | experiment,
~ 1 | ACC,
~ 1 | species,
~ 1 | phylo,
~ 1 | id
),
R = list(phylo = phylo.cor),
data = ef,
method = "REML"
)
)[["elapsed"]]
## Put runtimes into a tidy object
runtime_df <- data.frame(
model = c("glmmTMB_phylo", "rma.mv_phylo"),
runtime_sec = c(time_glmmTMB, time_rma)
)
runtime_df model runtime_sec
1 glmmTMB_phylo 33.11
2 rma.mv_phylo 467.67Some notes about speeding up model fitting with metafor: https://www.metafor-project.org/doku.php/tips:speeding_up_model_fitting
2.3 Location-scale meta-analysis
In a standard random-effects meta-analysis, the between-study variance \(\tau^2\) is assumed constant across studies. This is often unrealistic. Location–scale models relax this assumption by modelling both the location (the mean effect) and the scale (the amount of heterogeneity) as functions of study-level moderators. In practice, you can include different predictors in the mean and variance parts to test whether heterogeneity varies with study features (e.g., outcome type, design) or differs across subgroups. This yields a more flexible meta-analytic model that can explain, not just accommodate, between-study variability.
We will use a social-science meta-analysis of 48 studies on school-based writing-to-learn interventions. Each study contrasts an intervention group (greater emphasis on writing tasks) with a control group. The effect size measures used in this dataset is the standardised mean differences, where positive values indicate better academic performance under the intervention.
Code
dat <- dat.bangertdrowns2004
dat$id <- as.factor(dat$id)
dat$g <- rep(1, nrow(dat))Define V - the sampling error variance-covariance matrix assuming a within-study correlation of 0.6:
Code
V <- diag(dat$vi)
row.names(V) <- colnames(V) <- levels(dat$id)We fit a simple location-scale meta-analysis modelling the study sample size variable (ni) in the location and scale parts. We use dispformula in glmmTMB to specify the scale part of the model.
Code
fit_ls <- glmmTMB(yi ~ ni + equalto(0 + id|g, V),
dispformula = ~ ni,
data = dat,
REML = TRUE)
# > summary(fit_ls)
# Family: gaussian ( identity )
# Formula: yi ~ ni + equalto(0 + id | g, V)
# Dispersion: ~ni
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 54.6 62.0 -23.3 46.6 44
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# g id1 0.070 0.26458
# id2 0.126 0.35496 0.00
# id3 0.042 0.20494 0.00 0.00
# id4 0.019 0.13784 0.00 0.00 0.00
# id5 0.022 0.14832 0.00 0.00 0.00 0.00
# id6 0.009 0.09487 0.00 0.00 0.00 0.00 0.00
# id7 0.106 0.32558 0.00 0.00 0.00 0.00 0.00 0.00
# id8 0.007 0.08367 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# id9 0.040 0.20000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# id10 0.052 0.22804 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual NA NA
# Number of obs: 48, groups: g, 1
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 0.3016814 0.0667402 4.520 6.18e-06 ***
# ni -0.0005527 0.0001987 -2.781 0.00542 **
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Dispersion model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -0.960439 0.331888 -2.894 0.00381 **
# ni -0.004587 0.002540 -1.806 0.07095 .
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#
# View the estimates of the location-part of the model:
Code
summary(fit_ls)$coefficients$cond Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.3016813623 0.0667402453 4.520232 6.177203e-06
ni -0.0005526633 0.0001987119 -2.781230 5.415341e-03
attr(,"ddf")
[1] "asymptotic"View the estimates of the scale-part of the model:
Code
summary(fit_ls)$coefficients$disp Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.960439271 0.331888336 -2.893863 0.003805345
ni -0.004587144 0.002540198 -1.805821 0.070946268
attr(,"ddf")
[1] "asymptotic"We can fit the same model with metafor with the rma function and using the scale argument:
Code
fit_ls_rma <- rma(yi, vi,
mods = ~ ni,
scale = ~ ni,
control=list(REMLf=FALSE), #this provides the same log-likelihood as glmmTMB
data = dat)By default metafor uses the variance for the scale, whereas glmmTMB uses the standard deviation. Hence, to compare both outputs we need to times by x2 the estimates from the glmmTMB model (given the scale model uses a log link).
3 Useful notes and links
Changing optimization method in glmmTMB
In some cases, the default optimization method may not converge or can yield warning messages. In such cases, you can try alternative optimization methods by modifying the control = glmmTMBControl() argument. By default, glmmTMB uses the nlminb non-linear optimizer for parameter estimation.
Further information on optimization methods and troubleshooting convergence and warning messages can be found in the glmmTMB documentation:
Covariance structures in glmmTMB
The glmmTMB package offers a variety of covariance structures that can be used to model different types of dependencies in the data. A list of available covariance structures and their applications can be found in the package documentation:
4 References
Bangert-Drowns, R. L., Hurley, M. M., & Wilkinson, B. (2004). The effects of school-based writing-to-learn interventions on academic achievement: A meta-analysis. Review of Educational Research, 74, 29-58.
Brooks, M. E., Kristensen, K., van Benthem, K. J., Magnusson, A., Berg, C. W., Nielsen, A., Skaug, H. J., Mächler, M., & Bolker, B. M. (2017). glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R Journal, 9(2), 378–400. https://doi.org/10.32614/RJ-2017-066
Colditz, G. A., Brewer, T. F., Berkey, C. S., Wilson, M. E., Burdick, E., Fineberg, H. V., & Mosteller, F. (1994). Efficacy of BCG vaccine in the prevention of tuberculosis: Meta-analysis of the published literature. Journal of the American Medical Association, 271(9), 698–702. https://doi.org/10.1001/jama.1994.03510330076038
Felix, J. A., Stevenson, P. C., & Koricheva, J. (2023). Plant neighbourhood diversity effects on leaf traits: A meta-analysis. Functional Ecology, 37(12), 3150–3163. https://doi.org/10.1111/1365-2435.14441
Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1–48. https://doi.org/10.18637/jss.v036.i03
Viechtbauer,W.,&López-López, J.A. (2022). Location-scale models for meta-analysis. Research Synthesis Methods, 13(6), 697–715. https://doi.org/10.1002/jrsm.1562
White T, Noble D, Senior A, Hamilton W, Viechtbauer W (2022). metadat: Meta-Analysis Datasets. R package version 1.2-0, https://CRAN.R-project.org/package=metadat.
5 Session info
Code
library(sessioninfo)
library(details)
suppressWarnings(
details::details(session_info(), summary = 'Current session info', open = FALSE)
)Current session info
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.5.2 (2025-10-31 ucrt)
os Windows 11 x64 (build 26100)
system x86_64, mingw32
ui RTerm
language (EN)
collate English_Australia.utf8
ctype English_Australia.utf8
tz Australia/Sydney
date 2026-04-04
pandoc 3.6.3 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
quarto NA @ C:\\Users\\z5394590\\AppData\\Local\\Programs\\Quarto\\bin\\quarto.exe
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
ape * 5.8-1 2024-12-16 [1] CRAN (R 4.5.2)
bit 4.6.0 2025-03-06 [1] CRAN (R 4.5.2)
bit64 4.6.0-1 2025-01-16 [1] CRAN (R 4.5.2)
boot 1.3-32 2025-08-29 [2] CRAN (R 4.5.2)
cli 3.6.5 2025-04-23 [1] CRAN (R 4.5.2)
clipr 0.8.0 2022-02-22 [1] CRAN (R 4.5.2)
coda 0.19-4.1 2024-01-31 [1] CRAN (R 4.5.2)
codetools 0.2-20 2024-03-31 [2] CRAN (R 4.5.2)
crayon 1.5.3 2024-06-20 [1] CRAN (R 4.5.2)
curl 7.0.0 2025-08-19 [1] CRAN (R 4.5.2)
desc 1.4.3 2023-12-10 [1] CRAN (R 4.5.2)
details * 0.4.0 2025-02-09 [1] CRAN (R 4.5.2)
digest 0.6.39 2025-11-19 [1] CRAN (R 4.5.2)
dplyr * 1.2.0 2026-02-03 [1] CRAN (R 4.5.2)
emmeans 2.0.2 2026-03-05 [1] CRAN (R 4.5.2)
estimability 1.5.1 2024-05-12 [1] CRAN (R 4.5.2)
evaluate 1.0.5 2025-08-27 [1] CRAN (R 4.5.2)
farver 2.1.2 2024-05-13 [1] CRAN (R 4.5.2)
fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.5.2)
generics 0.1.4 2025-05-09 [1] CRAN (R 4.5.2)
glmmTMB * 1.1.14 2026-01-15 [1] CRAN (R 4.5.2)
glue 1.8.0 2024-09-30 [1] CRAN (R 4.5.2)
hms 1.1.4 2025-10-17 [1] CRAN (R 4.5.2)
htmltools 0.5.9 2025-12-04 [1] CRAN (R 4.5.2)
htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.5.2)
httr 1.4.8 2026-02-13 [1] CRAN (R 4.5.2)
jsonlite 2.0.0 2025-03-27 [1] CRAN (R 4.5.2)
kableExtra * 1.4.0 2024-01-24 [1] CRAN (R 4.5.2)
knitr 1.51 2025-12-20 [1] CRAN (R 4.5.2)
lattice 0.22-7 2025-04-02 [2] CRAN (R 4.5.2)
lifecycle 1.0.5 2026-01-08 [1] CRAN (R 4.5.2)
lme4 1.1-38 2025-12-02 [1] CRAN (R 4.5.2)
magrittr 2.0.4 2025-09-12 [1] CRAN (R 4.5.2)
MASS * 7.3-65 2025-02-28 [2] CRAN (R 4.5.2)
mathjaxr 2.0-0 2025-12-01 [1] CRAN (R 4.5.2)
Matrix * 1.7-4 2025-08-28 [1] CRAN (R 4.5.2)
metadat * 1.4-0 2025-02-04 [1] CRAN (R 4.5.2)
metafor * 4.8-0 2025-01-28 [1] CRAN (R 4.5.2)
mgcv 1.9-4 2025-11-07 [1] CRAN (R 4.5.2)
minqa 1.2.8 2024-08-17 [1] CRAN (R 4.5.2)
multcomp 1.4-29 2025-10-20 [1] CRAN (R 4.5.2)
mvtnorm 1.3-3 2025-01-10 [1] CRAN (R 4.5.2)
nlme 3.1-168 2025-03-31 [2] CRAN (R 4.5.2)
nloptr 2.2.1 2025-03-17 [1] CRAN (R 4.5.2)
numDeriv * 2016.8-1.1 2019-06-06 [1] CRAN (R 4.5.2)
otel 0.2.0 2025-08-29 [1] CRAN (R 4.5.2)
pillar 1.11.1 2025-09-17 [1] CRAN (R 4.5.2)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.5.2)
png 0.1-8 2022-11-29 [1] CRAN (R 4.5.2)
prettyunits 1.2.0 2023-09-24 [1] CRAN (R 4.5.2)
progress 1.2.3 2023-12-06 [1] CRAN (R 4.5.2)
R6 2.6.1 2025-02-15 [1] CRAN (R 4.5.2)
rbibutils 2.4.1 2026-01-21 [1] CRAN (R 4.5.2)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.5.2)
Rcpp 1.1.1 2026-01-10 [1] CRAN (R 4.5.2)
Rdpack 2.6.6 2026-02-08 [1] CRAN (R 4.5.2)
readr * 2.2.0 2026-02-19 [1] CRAN (R 4.5.2)
reformulas 0.4.4 2026-02-02 [1] CRAN (R 4.5.2)
rentrez 1.2.4 2025-06-11 [1] CRAN (R 4.5.2)
rlang 1.1.7 2026-01-09 [1] CRAN (R 4.5.2)
rmarkdown 2.30 2025-09-28 [1] CRAN (R 4.5.2)
rncl 0.8.9 2026-01-21 [1] CRAN (R 4.5.2)
rotl * 3.1.1 2026-01-15 [1] CRAN (R 4.5.2)
rstudioapi 0.18.0 2026-01-16 [1] CRAN (R 4.5.2)
sandwich 3.1-1 2024-09-15 [1] CRAN (R 4.5.2)
scales 1.4.0 2025-04-24 [1] CRAN (R 4.5.2)
sessioninfo * 1.2.3 2025-02-05 [1] CRAN (R 4.5.2)
stringi 1.8.7 2025-03-27 [1] CRAN (R 4.5.2)
stringr 1.6.0 2025-11-04 [1] CRAN (R 4.5.2)
survival 3.8-3 2024-12-17 [2] CRAN (R 4.5.2)
svglite 2.2.2 2025-10-21 [1] CRAN (R 4.5.2)
systemfonts 1.3.1 2025-10-01 [1] CRAN (R 4.5.2)
textshaping 1.0.4 2025-10-10 [1] CRAN (R 4.5.2)
TH.data 1.1-5 2025-11-17 [1] CRAN (R 4.5.2)
tibble 3.3.1 2026-01-11 [1] CRAN (R 4.5.2)
tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.5.2)
TMB 1.9.19 2025-12-15 [1] CRAN (R 4.5.2)
tzdb 0.5.0 2025-03-15 [1] CRAN (R 4.5.2)
vctrs 0.7.1 2026-01-23 [1] CRAN (R 4.5.2)
viridisLite 0.4.3 2026-02-04 [1] CRAN (R 4.5.2)
vroom 1.7.0 2026-01-27 [1] CRAN (R 4.5.2)
withr 3.0.2 2024-10-28 [1] CRAN (R 4.5.2)
xfun 0.55 2025-12-16 [1] CRAN (R 4.5.2)
XML 3.99-0.20 2025-11-08 [1] CRAN (R 4.5.2)
xml2 1.5.1 2025-12-01 [1] CRAN (R 4.5.2)
xtable 1.8-4 2019-04-21 [1] CRAN (R 4.5.2)
yaml 2.3.12 2025-12-10 [1] CRAN (R 4.5.2)
zoo 1.8-15 2025-12-15 [1] CRAN (R 4.5.2)
[1] C:/Users/z5394590/AppData/Local/R/win-library/4.5
[2] C:/Program Files/R/R-4.5.2/library
* ── Packages attached to the search path.
──────────────────────────────────────────────────────────────────────────────