Code
# load libraries
library(MASS)
library(kableExtra)
library(metafor)
library(glmmTMB)
library(tidyverse)Fitting meta-analysis with the glmmTMB R package
Supplementary material for Research Synthesis Methods
Coralie Williams 
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: 16 July 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 libraries
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 ✅
Derive adjusted \(t\)-distribution p-values
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 |
|---|---|---|---|---|---|---|
| 0.0056193 | 0.1312196 | 0.0428234 | 9 | 0.9667773 | -0.29122 | 0.3024585 |
2 Illustrative examples
To illustrate the range of meta-analytic models that can be fitted with glmmTMB, we re-analyse four published datasets from different research domains:
- Bivariate meta-analysis — applied to a medical dataset on vaccine effectiveness.
- Network meta-analysis — applied to a psychology dataset comparing multiple interventions.
- Phylogenetic meta-analysis — applied to an ecology and evolution dataset on the effects of plant diversity.
- Location–scale meta-analysis — applied to a social science dataset evaluating school-based learning programs.
For each case study, we demonstrate model specification, output, and interpretation using glmmTMB. We also fit equivalent models using the metafor package (when possible) and compare the resulting estimates across packages.
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)
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 | out2 | yi | vi | obs | g |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Aronson | 1948 | 4 | 119 | 11 | 128 | 44 | random | 1 | v | 4 | 119 | -3.392829 | 0.2584034 | 1 | 1 |
| 1 | Aronson | 1948 | 4 | 119 | 11 | 128 | 44 | random | 1 | c | 11 | 128 | -2.454135 | 0.0987216 | 2 | 1 |
| 2 | Ferguson & Simes | 1949 | 6 | 300 | 29 | 274 | 55 | random | 2 | v | 6 | 300 | -3.912023 | 0.1700000 | 3 | 1 |
| 2 | Ferguson & Simes | 1949 | 6 | 300 | 29 | 274 | 55 | random | 2 | c | 29 | 274 | -2.245832 | 0.0381324 | 4 | 1 |
| 3 | Rosenthal et al | 1960 | 3 | 228 | 11 | 209 | 42 | random | 3 | v | 3 | 228 | -4.330733 | 0.3377193 | 5 | 1 |
| 3 | Rosenthal et al | 1960 | 3 | 228 | 11 | 209 | 42 | random | 3 | c | 11 | 209 | -2.944439 | 0.0956938 | 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
#
# # We can also get the estimated marginal means
# emmeans::emmeans(bv_fit, ~ group) |>
# summary(infer=TRUE, type="response")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(coef(summary(bv_fit))$cond[,"Estimate"][1]),
se.b0 = as.numeric(coef(summary(bv_fit))$cond[,"Std. Error"][1]),
est.b1 = as.numeric(coef(summary(bv_fit))$cond[,"Estimate"][2]),
se.b1 = as.numeric(coef(summary(bv_fit))$cond[,"Std. Error"][2]),
tau2_1 = var_con, #study-level heterogeneity for group=control
tau2_2 = var_exp, #study-level heterogeneity for group=exp
rho = rho.tmb)
bv_out| est.b0 | se.b0 | est.b1 | se.b1 | tau2_1 | tau2_2 | rho |
|---|---|---|---|---|---|---|
| -4.095985 | 0.4529815 | -0.7413877 | 0.1892434 | 2.6173 | 1.548603 | 0.9450462 |
This 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)
# summary(bv_fit_rma)
#
# Multivariate Meta-Analysis Model (k = 26; method: REML)
#
# logLik Deviance AIC BIC AICc
# -31.5017 63.0033 73.0033 78.8936 76.3367
#
# Variance Components:
#
# outer factor: study (nlvls = 13)
# inner factor: group (nlvls = 2)
#
# estim sqrt k.lvl fixed level
# tau^2.1 2.6173 1.6178 13 no c
# tau^2.2 1.5486 1.2444 13 no v
#
# rho.c rho.v c v
# c 1 - 13
# v 0.9450 1 no -
#
# Test for Residual Heterogeneity:
# QE(df = 24) = 5270.3863, p-val < .0001
#
# Test of Moderators (coefficient 2):
# QM(df = 1) = 15.5470, p-val < .0001
#
# Model Results:
#
# estimate se zval pval ci.lb ci.ub
# intrcpt -4.0960 0.4529 -9.0432 <.0001 -4.9837 -3.2082 ***
# groupv -0.7414 0.1880 -3.9430 <.0001 -1.1099 -0.3729 ***
#
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1We 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 |
|---|---|---|---|---|---|---|
| -4.095985 | 0.4529369 | -0.7413877 | 0.1880277 | 2.6173 | 1.548603 | 0.9450462 |
We can also fit the GLMM specification of the same bivariate model by taking the joint distribution of the probability of an event in the treatment and control groups of each study. This model is sometimes called the “Van Houwelingen bivariate” model
Code
dat_long <- dat_long |>
mutate(vacc = as.integer(group == "v"),
cont = as.integer(group == "c"))Now fit the model with glmmTMB function:
Code
bv_joint_fit <- glmmTMB(cbind(out1,out2) ~ factor(vacc) + (cont+vacc-1|study),
family = binomial(link = "logit"),
data=dat_long)
# > summary(bv_joint_fit)
#
# Family: binomial ( logit )
# Formula: cbind(out1, out2) ~ factor(vacc) + (cont + vacc - 1 | study)
# Data: dat_long
#
# AIC BIC logLik -2*log(L) df.resid
# 257.4 263.7 -123.7 247.4 21
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study cont 2.438 1.561
# vacc 1.436 1.198 0.95
# Number of obs: 26, groups: study, 13
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -4.1234 0.4374 -9.427 < 2e-16 ***
# factor(vacc)1 -0.7417 0.1834 -4.043 5.27e-05 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Code
vc <- VarCorr(bv_joint_fit)
var_cont <- vc$cond$study["cont", "vacc"]
var_vacc <- vc$cond$study["vacc", "vacc"]
rho.tmb <- vc$cond$study["cont", "vacc"] / sqrt(var_cont * var_vacc)
bv_joint_out <- data.frame(est.b0 = as.numeric(coef(summary(bv_joint_fit))$cond[,"Estimate"][1]),
se.b0 = as.numeric(coef(summary(bv_joint_fit))$cond[,"Std. Error"][1]),
est.b1 = as.numeric(coef(summary(bv_joint_fit))$cond[,"Estimate"][2]),
se.b1 = as.numeric(coef(summary(bv_joint_fit))$cond[,"Std. Error"][2]),
tau2_1 = var_cont, #study-level heterogeneity for group=control
tau2_2 = var_vacc, #study-level heterogeneity for group=vaccine
rho = rho.tmb)
bv_joint_out| est.b0 | se.b0 | est.b1 | se.b1 | tau2_1 | tau2_2 | rho |
|---|---|---|---|---|---|---|
| -4.123398 | 0.4374109 | -0.7417176 | 0.1834433 | 1.768658 | 1.43569 | 1.10992 |
2.2 Network meta-analysis
When several treatments are compared across a set of studies, a network meta-analysis can be used. It allows treatments to be compared even when they were not directly compared in the same study, provided they are connected through other treatment comparisons. This allows all treatments in the network to be compared within one model.
We use the Hasselblad (1998) dataset from the metadat R package. It contains 24 studies comparing four smoking-cessation conditions: no contact, self-help, individual counselling and group counselling. For each study arm, the dataset records the number of individuals who remained abstinent for 6 to 12 months and the total number of individuals. Most studies compared two conditions, while two studies compared three.
Code
dat <- metadat::dat.hasselblad1998
head(dat)| id | study | authors | year | trt | xi | ni |
|---|---|---|---|---|---|---|
| 1 | 1 | Reid et al. | 1974 | no_contact | 75 | 731 |
| 2 | 1 | Reid et al. | 1974 | ind_counseling | 363 | 714 |
| 3 | 2 | Cottraux et al. | 1983 | no_contact | 9 | 140 |
| 4 | 2 | Cottraux et al. | 1983 | ind_counseling | 23 | 140 |
| 5 | 2 | Cottraux et al. | 1983 | grp_counseling | 10 | 138 |
| 6 | 3 | Slama et al. | 1990 | no_contact | 2 | 106 |
Before fitting the model, we visualise the network of direct treatment comparisons. Each node represents a treatment, and each edge represents a direct comparison between two treatments. Wider edges indicate that more studies contributed to that comparison.
Code
library(igraph, warn.conflicts=FALSE)
pairs <- data.frame(do.call(rbind, sapply(split(dat$trt, dat$study), function(x) t(combn(x,2)))),
stringsAsFactors=FALSE)
lvls <- c("no_contact", "self_help", "ind_counseling", "grp_counseling")
pairs$X1 <- factor(pairs$X1, levels=lvls)
pairs$X2 <- factor(pairs$X2, levels=lvls)
tab <- table(pairs[,1], pairs[,2])
g <- graph_from_adjacency_matrix(tab, mode = "plus", weighted=TRUE, diag=FALSE)
vertex_attr(g, "name") <- c("No Contact",
"Self-Help",
"Individual\nCounseling",
"Group\nCounseling")
plot(g,
edge.curved = FALSE,
edge.width = E(g)$weight,
edge.color = "gray70",
layout = layout_on_grid,
vertex.size = 45,
vertex.color = "lightgray",
vertex.frame.color = "gray40",
vertex.label.color = "black",
vertex.label.font = 2,
vertex.label.family = "sans",
vertex.label.cex = 0.9) 
We then format the treatment variable and use no contact as the reference condition.
Code
# formatting for output
dat$trt <- factor(dat$trt, levels=c("no_contact", "self_help", "ind_counseling", "grp_counseling"))
levels(dat$trt) <- paste0(" ", levels(dat$trt))We fit a one-stage arm-based network meta-analysis using a binomial GLMM with a logit link. The model uses the numbers of abstinent (xi) and non-abstinent individuals (ni - xi) in each study arm. The treatment effects compare each counselling condition with no contact, while the study-level random effects account for variation between studies.
This is an arm-based (AB) likelihood: the raw arm-level counts are modelled directly via an exact binomial distribution, rather than first reducing each study to estimated log-odds ratios pooled under a normal approximation as in the traditional contrast-based (CB) likelihood (Salanti, 2012; White et al., 2019).
Model 1: independent arm-level random effects
The first model uses treatment-specific between-study variances with no correlation between arms of the same study (using the diag covariance structure).
Code
nma_fit <- glmmTMB(cbind(xi, ni-xi) ~ trt + diag(0 + trt | study),
family = binomial(link = "logit"),
data=dat)
# > summary(nma_fit)
#
#
# Family: binomial ( logit )
# Formula: cbind(xi, ni - xi) ~ trt + diag(0 + trt | study)
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 402.6 417.9 -193.3 386.6 42
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev.
# study trt no_contact 0.2161 0.4648
# trt self_help 0.6245 0.7902
# trt ind_counseling 0.5108 0.7147
# trt grp_counseling 0.3930 0.6269
# Number of obs: 50, groups: study, 24
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -2.5651 0.1363 -18.813 < 2e-16 ***
# trt self_help 0.6118 0.3683 1.661 0.096652 .
# trt ind_counseling 0.8857 0.2203 4.020 5.82e-05 ***
# trt grp_counseling 1.1390 0.3211 3.548 0.000389 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Model 2: compound-symmetry random effects
The second model adds a single common correlation \(\rho\) between arms of the same study (using the cs compound-symmetry covariance structure), capturing the shared study baseline.
Code
nma_cs_fit <- glmmTMB(cbind(xi, ni-xi) ~ trt + cs(0 + trt | study),
family = binomial(link = "logit"),
data=dat)
# > summary(nma_cs_fit)
#
#
# Family: binomial ( logit )
# Formula: cbind(xi, ni - xi) ~ trt + cs(0 + trt | study)
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 401.3 418.5 -191.7 383.3 41
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study trt no_contact 0.1858 0.4311 0.49 (cs)
# trt self_help 0.4847 0.6962
# trt ind_counseling 0.5443 0.7377
# trt grp_counseling 0.5051 0.7107
# Number of obs: 50, groups: study, 24
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -2.4902 0.1273 -19.567 < 2e-16 ***
# trt self_help 0.5268 0.2912 1.809 0.070452 .
# trt ind_counseling 0.7796 0.1855 4.204 2.63e-05 ***
# trt grp_counseling 1.0707 0.3163 3.385 0.000711 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Code
anova(nma_fit, nma_cs_fit)| Df | AIC | BIC | logLik | deviance | Chisq | Chi Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| nma_fit | 8 | 402.5668 | 417.8629 | -193.2834 | 386.5668 | NA | NA | NA |
| nma_cs_fit | 9 | 401.3110 | 418.5193 | -191.6555 | 383.3110 | 3.255704 | 1 | 0.0711754 |
The likelihood ratio test weakly favours the compound-symmetry model (χ²₁ = 3.26, p = 0.071; ΔAIC = 1.3), and we retain it as the primary analysis since arms within a study plausibly share a common baseline. Note that the slightly larger estimates under the diagonal model reflect within-study correlation being misattributed to treatment contrasts when ignored.
Code
est <- exp(coef(summary(nma_cs_fit))$cond[,1])
ci_prof <- exp(confint(nma_cs_fit, parm="beta_", method = "profile"))
results_nma <- data.frame(
Estimate = est,
Lower_CI = ci_prof[, 1],
Upper_CI = ci_prof[, 2]
)
results_nma| Estimate | Lower_CI | Upper_CI | |
|---|---|---|---|
| (Intercept) | 0.082893 | 0.0619236 | 0.1055801 |
| trt self_help | 1.693519 | 0.8972999 | 3.3346474 |
| trt ind_counseling | 2.180618 | 1.5205723 | 3.2687439 |
| trt grp_counseling | 2.917375 | 1.4199303 | 6.1399210 |
The above confint() call computes profile likelihood confidence intervals, which are based on the shape of the likelihood rather than the Wald (normal) approximation and are more reliable in moderate samples.
Exponentiating the profile likelihood intervals gives odds ratios versus no contact of 1.69 (95% CI 0.90–3.33) for self-help, 2.18 (1.52–3.27) for individual counselling, and 2.92 (1.42–6.14) for group counselling. Both counselling formats show evidence of benefit, while the interval for self-help includes 1, so its evidence is inconclusive.
2.3 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 |
|---|---|---|---|---|
| 0.2044 | 0 | 0.1063 | 0.2958 | 0.4256 |
This 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
-1368.8279 2737.6557 2751.6557 2785.6731 2751.7743
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0615 0.2479 39 no experiment no
sigma^2.2 0.0000 0.0001 39 no ACC no
sigma^2.3 0.2125 0.4610 102 no species no
sigma^2.4 0.1784 0.4224 102 no phylo yes
sigma^2.5 0.0266 0.1631 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.8209, p-val = 0.6633
Model Results:
estimate se zval pval ci.lb
NfixingN-fixing neighbours 0.1759 0.2459 0.7154 0.4744 -0.3060
NfixingNo N-fixing neighbours 0.2091 0.2361 0.8856 0.3758 -0.2536
ci.ub
NfixingN-fixing neighbours 0.6579
NfixingNo N-fixing neighbours 0.6718
---
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 63.72
2 rma.mv_phylo 1076.51Some notes about speeding up model fitting with metafor: https://www.metafor-project.org/doku.php/tips:speeding_up_model_fitting
2.4 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
Robust variance estimation
vcovHC() is a function in glmmTMB that computes robust (sandwich) variance-covariance matrices for the fixed effects estimates. This can be particularly useful when there is concern about model misspecification or when the data exhibit heteroscedasticity or clustering that is not fully accounted for by the random effects structure.
More information is provided here: Post-model-fitting procedures with glmmTMB
Important note: currently glmmTMB doesn’t support small sample adjustments (though it may be added in future developments).
Changing optimization method
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:
Other covariance structures
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:
Orchard plots with glmmTMB
Orchard plots are extensions of forest plots that display model estimates, confidence and prediction intervals, and the distribution of individual effect sizes (Nakagawa, 2021). To demonstrate how to produce orchard plots with the orchaRd package, we use the phylogenetic meta-analysis example and plot the estimated effects for the two plant groups: N-fixing and non-N-fixing species.
Code
## As of July 2026: you need to use the following branch
#detach("package:orchaRd", unload = TRUE)
#remove.packages("orchaRd")
#devtools::install_github("daniel1noble/orchaRd", ref = "glmmtmb-rma-adapter", force = TRUE)
library(orchaRd)
tmb.res <- glmmTMB_to_rma(fit_phylo,
yi = "yi",
vi = "vi",
data = ef)
tmb.Nfixing <- orchaRd::mod_results(
tmb.res,
mod = "Nfixing",
group = "experiment"
)
orchaRd::orchard_plot(
tmb.Nfixing,
xlab = "lnRR",
angle = 45,
g = FALSE
)
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
Hasselblad, V. (1998). Meta-analysis of multitreatment studies. Medical Decision Making, 18(1), 37–43. https://doi.org/10.1177/0272989X9801800110
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. , and 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
Salanti, G. (2012). Indirect and mixed-treatment comparison, network, or multiple-treatments meta-analysis: Many names, many benefits, many concerns for the next generation evidence synthesis tool. Research Synthesis Methods, 3(2), 80–97.
Nakagawa, S., Lagisz, M., O’Dea, R. E., Pottier, P., Rutkowska, J., Senior, A. M., Yang, Y., & Noble, D. W. A. (2023). orchaRd 2.0: An R package for visualising meta-analyses with orchard plots. Methods in Ecology and Evolution, 14(8). https://doi.org/10.1111/2041-210X.14152
White, I. R., Turner, R. M., Karahalios, A., & Salanti, G. (2019). A comparison of arm-based and contrast-based models for network meta-analysis. Statistics in Medicine, 38(27), 5197–5213.
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 26200)
system x86_64, mingw32
ui RTerm
language (EN)
collate English_Australia.utf8
ctype English_Australia.utf8
tz Australia/Sydney
date 2026-07-16
pandoc 3.8.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)
beeswarm 0.4.0 2021-06-01 [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.1 2026-05-25 [1] CRAN (R 4.5.3)
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.3 2026-04-09 [1] CRAN (R 4.5.3)
estimability 2.0.0 2026-06-26 [1] CRAN (R 4.5.3)
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)
forcats * 1.0.1 2025-09-25 [1] CRAN (R 4.5.2)
generics 0.1.4 2025-05-09 [1] CRAN (R 4.5.2)
ggbeeswarm 0.7.3 2025-11-29 [1] CRAN (R 4.5.2)
ggplot2 * 4.0.2 2026-02-03 [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)
gtable 0.3.6 2024-10-25 [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)
igraph * 2.2.1 2025-10-27 [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)
labeling 0.4.3 2023-08-29 [1] CRAN (R 4.5.2)
latex2exp 0.9.8 2026-01-09 [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)
lubridate * 1.9.4 2024-12-08 [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)
orchaRd * 2.1.3 2026-07-15 [1] Github (daniel1noble/orchaRd@73c1bf7)
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)
purrr * 1.2.1 2026-01-09 [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.31 2026-03-26 [1] CRAN (R 4.5.3)
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.19.0 2026-06-11 [1] CRAN (R 4.5.3)
S7 0.2.1 2025-11-14 [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)
tidyr * 1.3.2 2025-12-19 [1] CRAN (R 4.5.2)
tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.5.2)
tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.5.2)
timechange 0.3.0 2024-01-18 [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)
vipor 0.4.7 2023-12-18 [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.3 2026-06-19 [1] CRAN (R 4.5.3)
xfun 0.57 2026-03-20 [1] CRAN (R 4.5.3)
XML 3.99-0.20 2025-11-08 [1] CRAN (R 4.5.2)
xml2 1.5.2 2026-01-17 [1] CRAN (R 4.5.3)
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.
──────────────────────────────────────────────────────────────────────────────