metaGLMM fits generalized linear mixed-effects
likelihoods to aggregate study summaries. Each row contains a response,
its sampling variance vi, and a study size or exposure
ni. A formula supplies the fixed-effect model matrix, so
interactions, factors, and offsets use ordinary R model syntax.
The initial release provides four built-in families:
| Family | Typical response | Common links |
|---|---|---|
gaussian() |
continuous effect estimate | identity, log |
binomial() |
proportion | logit, probit,
cloglog |
poisson() |
rate | log, identity |
Gamma() |
positive mean | log, inverse |
The package uses the fast numerical path for built-in families when
fast = TRUE. An R-level custom family uses QMC integration
with fast = FALSE; see the custom-family vignette.
The following data are deliberately small so that the example is quick to run and easy to reproduce.
dat <- data.frame(
study = paste0("Study ", 1:6),
estimate = c(-0.80, -0.30, 0.10, 0.70, 1.10, -0.50),
moderator = c(0, 0, 1, 1, 0, 1),
vi = c(0.04, 0.05, 0.03, 0.06, 0.04, 0.05),
ni = c(100, 90, 120, 80, 110, 95)
)
dat
#> study estimate moderator vi ni
#> 1 Study 1 -0.8 0 0.04 100
#> 2 Study 2 -0.3 0 0.05 90
#> 3 Study 3 0.1 1 0.03 120
#> 4 Study 4 0.7 1 0.06 80
#> 5 Study 5 1.1 0 0.04 110
#> 6 Study 6 -0.5 1 0.05 95Fit a random-effects model with an intercept and a moderator. The study estimates intentionally span a broad range, so the between-study component is estimable rather than collapsing to the boundary:
library(metaGLMM)
fit <- metaGLMM(
estimate ~ moderator,
data = dat,
vi = dat$vi,
ni = dat$ni,
tau2 = NA,
family = gaussian(link = "identity"),
tau2_var = TRUE,
fast = TRUE
)
fit
#> Aggregate-data generalized linear mixed-effects meta-analysis
#> Family:gaussian(identity)
#> Formula: estimate ~ moderator
#>
#> Fixed effects:
#> (Intercept) moderator
#> 0.002241 0.093444
#>
#> Heterogeneity: tau^2 = 0.3996 , tau = 0.6322
summary(fit)
#> Aggregate-data generalized linear mixed-effects meta-analysis
#>
#> Call:
#> metaGLMM(formula = estimate ~ moderator, data = dat, vi = dat$vi,
#> ni = dat$ni, tau2 = NA, family = gaussian(link = "identity"),
#> tau2_var = TRUE, fast = TRUE)
#>
#> Family:gaussian(identity)
#> Random-effect structure: row_intercept
#> Integration: closed_form
#>
#> Fixed effects:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.002241 0.384242 0.006 0.995
#> moderator 0.093444 0.544338 0.172 0.864
#>
#> Heterogeneity:
#> tau^2: 0.3996
#> tau: 0.6322
#>
#> Fit diagnostics:
#> Optimizer: L-BFGS-B
#> Convergence code: 0
#> Boundary tau^2: no
#> Positive-definite Hessian: yescoef() and vcov() expose fixed effects
only. The estimated heterogeneity is available separately as
tau2 and tau:
coef(fit)
#> (Intercept) moderator
#> 0.00224059 0.09344402
vcov(fit)
#> (Intercept) moderator
#> (Intercept) 0.1476419 -0.1476450
#> moderator -0.1476450 0.2963039
fit$tau2
#> [1] 0.3996363
fit$tau
#> [1] 0.6321679
confint(fit, method = "wald")
#> lower upper
#> (Intercept) -0.7508598 0.755341
#> moderator -0.9734389 1.160327
stopifnot(is.finite(fit$tau), fit$tau > 0)The object also retains the model frame and model matrix used for fitting. This makes factor contrasts and coefficient names explicit:
formula(fit)
#> estimate ~ moderator
terms(fit)
#> estimate ~ moderator
#> attr(,"variables")
#> list(estimate, moderator)
#> attr(,"factors")
#> moderator
#> estimate 0
#> moderator 1
#> attr(,"term.labels")
#> [1] "moderator"
#> attr(,"order")
#> [1] 1
#> attr(,"intercept")
#> [1] 1
#> attr(,"response")
#> [1] 1
#> attr(,".Environment")
#> <environment: R_GlobalEnv>
#> attr(,"predvars")
#> list(estimate, moderator)
#> attr(,"dataClasses")
#> estimate moderator
#> "numeric" "numeric"
model.frame(fit)
#> estimate moderator
#> 1 -0.8 0
#> 2 -0.3 0
#> 3 0.1 1
#> 4 0.7 1
#> 5 1.1 0
#> 6 -0.5 1
colnames(model.matrix(fit))
#> [1] "(Intercept)" "moderator"For a Gamma likelihood, provide positive aggregate means and their sampling variances. Here the response is on the original positive scale and the log link makes the treatment coefficient multiplicative.
gamma_dat <- data.frame(
study = paste0("Study ", 1:6),
mean = c(9.6, 9.9, 13.8, 16.5, 8.2, 14.6),
sd = c(0.7, 8.3, 4.2, 7.4, 4.3, 2.2),
n = c(20, 25, 23, 18, 75, 20),
treatment = c(1, 1, 1, 1, 1, 0)
)
gamma_dat$vi <- (gamma_dat$sd / gamma_dat$mean)^2
gamma_fit <- metaGLMM(
mean ~ treatment,
data = gamma_dat,
vi = gamma_dat$vi,
ni = gamma_dat$n,
tau2 = NA,
family = Gamma(link = "log"),
tau2_var = TRUE,
start_tau2 = 0.1,
fast = TRUE
)
summary(gamma_fit)
#> Aggregate-data generalized linear mixed-effects meta-analysis
#>
#> Call:
#> metaGLMM(formula = mean ~ treatment, data = gamma_dat, vi = gamma_dat$vi,
#> ni = gamma_dat$n, tau2 = NA, family = Gamma(link = "log"),
#> tau2_var = TRUE, fast = TRUE, start_tau2 = 0.1)
#>
#> Family:Gamma(log)
#> Random-effect structure: row_intercept
#> Integration: ghq
#>
#> Fixed effects:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 2.6816 0.2367 11.330 <2e-16 ***
#> treatment -0.2652 0.2649 -1.001 0.317
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Heterogeneity:
#> tau^2: 0.04776
#> tau: 0.2185
#>
#> Fit diagnostics:
#> Optimizer: L-BFGS-B
#> Convergence code: 0
#> Boundary tau^2: no
#> Positive-definite Hessian: yes
stopifnot(is.finite(gamma_fit$tau), gamma_fit$tau > 0)Use the family-specific vignettes for binary and Poisson summaries, then see the inference and prediction vignettes for intervals, new-study predictions, and study-level displays.