Gaussian aggregate outcomes

Weighted study estimates

A Gaussian analysis uses one aggregate estimate and its sampling variance per row. The sample size ni is retained in the fitted object and can be used by downstream code, while the Gaussian weighted likelihood is driven by vi.

gaussian_dat <- data.frame(
  study = paste0("Study ", 1:8),
  estimate = c(-0.80, -0.30, 0.10, 0.70, 1.10, -0.50, 0.40, 1.30),
  moderator = c(0, 1, 0, 1, 0, 1, 0, 1),
  vi = c(0.04, 0.05, 0.03, 0.06, 0.04, 0.05, 0.03, 0.05),
  ni = c(100, 90, 120, 80, 110, 95, 130, 85)
)
gaussian_dat
#>     study estimate moderator   vi  ni
#> 1 Study 1     -0.8         0 0.04 100
#> 2 Study 2     -0.3         1 0.05  90
#> 3 Study 3      0.1         0 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  95
#> 7 Study 7      0.4         0 0.03 130
#> 8 Study 8      1.3         1 0.05  85

Basic random-effects analysis and forest plot

For a single-group meta-analysis, fit an intercept-only model. Each input row is already one study estimate, so as_metafor_data() can pass the study estimates and sampling variances directly to forest().

gaussian_fit <- metaGLMM(
  estimate ~ 1,
  data = gaussian_dat,
  vi = gaussian_dat$vi,
  ni = gaussian_dat$ni,
  tau2 = NA,
  family = gaussian(link = "identity"),
  tau2_var = TRUE,
  fast = TRUE
)

summary(gaussian_fit)
#> Aggregate-data generalized linear mixed-effects meta-analysis
#> 
#> Call:
#> metaGLMM(formula = estimate ~ 1, data = gaussian_dat, vi = gaussian_dat$vi, 
#>     ni = gaussian_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.2484     0.2502   0.993    0.321
#> 
#> Heterogeneity:
#>   tau^2: 0.4572 
#>   tau:   0.6762 
#> 
#> Plug-in prediction interval for a future underlying true effect:
#>     fit  lower upper
#>  0.2484 -1.165 1.662
#> 
#> Fit diagnostics:
#>   Optimizer: L-BFGS-B 
#>   Convergence code: 0 
#>   Boundary tau^2: no 
#>   Positive-definite Hessian: yes
coef(gaussian_fit)
#> (Intercept) 
#>    0.248429
confint(gaussian_fit, method = "wald")
#>                  lower     upper
#> (Intercept) -0.2419502 0.7388083
stopifnot(is.finite(gaussian_fit$tau), gaussian_fit$tau > 0)

The study rows show the supplied Gaussian estimates and their sampling intervals. The lower rows compare the three supported fixed-effect intervals and the plug-in prediction interval.

gaussian_studies <- as_metafor_data(
  gaussian_fit, labels = gaussian_dat$study
)
head(gaussian_studies)
#>     yi   vi slab
#> 1 -0.8 0.04    1
#> 2 -0.3 0.05    2
#> 3  0.1 0.03    3
#> 4  0.7 0.06    4
#> 5  1.1 0.04    5
#> 6 -0.5 0.05    6

forest(
  gaussian_fit,
  labels = gaussian_dat$study,
  xlab = "Effect estimate",
  ci_methods = c("Wald", "profile", "SBC")
)

Fixed-effect reference

Setting tau2_var = FALSE with tau2 = 0 requests the exact no-random-effect model. Its fixed-effect estimates can be compared with the weighted least squares reference in base R.

fixed_fit <- metaGLMM(
  estimate ~ moderator,
  data = gaussian_dat,
  vi = gaussian_dat$vi,
  ni = gaussian_dat$ni,
  tau2 = 0,
  tau2_var = FALSE,
  family = gaussian(link = "identity"),
  fast = TRUE
)

X <- model.matrix(~ moderator, data = gaussian_dat)
W <- sweep(X, 1, gaussian_dat$vi, "/")
beta_wls <- solve(crossprod(X, W),
                  crossprod(X, gaussian_dat$estimate / gaussian_dat$vi))
names(beta_wls) <- colnames(X)

cbind(metaGLMM = coef(fixed_fit), weighted_least_squares = beta_wls)
#>               metaGLMM           
#> (Intercept) 0.20714286 0.20714286
#> moderator   0.07546584 0.07546584

This no-random-effect fit is an analytical reference only. The next section uses the same study summaries while estimating a non-zero between-study component.

Estimating heterogeneity

Set tau2 = NA and tau2_var = TRUE to estimate the non-negative between-study variance. tau is its square root.

random_fit <- metaGLMM(
  estimate ~ moderator,
  data = gaussian_dat,
  vi = gaussian_dat$vi,
  ni = gaussian_dat$ni,
  tau2 = NA,
  tau2_var = TRUE,
  family = gaussian(link = "identity"),
  fast = TRUE
)

summary(random_fit)
#> Aggregate-data generalized linear mixed-effects meta-analysis
#> 
#> Call:
#> metaGLMM(formula = estimate ~ moderator, data = gaussian_dat, 
#>     vi = gaussian_dat$vi, ni = gaussian_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.20051    0.34997   0.573    0.567
#> moderator    0.09754    0.49933   0.195    0.845
#> 
#> Heterogeneity:
#>   tau^2: 0.455 
#>   tau:   0.6745 
#> 
#> Fit diagnostics:
#>   Optimizer: L-BFGS-B 
#>   Convergence code: 0 
#>   Boundary tau^2: no 
#>   Positive-definite Hessian: yes
random_fit$tau2
#> [1] 0.4549558
random_fit$tau
#> [1] 0.6745041
logLik(random_fit)
#> 'log Lik.' -1.231958 (df=3)
nobs(random_fit)
#> [1] 8
stopifnot(is.finite(random_fit$tau), random_fit$tau > 0)

The alternative tau2_param = "log_tau2" uses an unconstrained optimization parameterization while reporting tau2 on its original scale. Both parameterizations expose the same fixed-effect API.

random_fit_log <- metaGLMM(
  estimate ~ moderator,
  data = gaussian_dat,
  vi = gaussian_dat$vi,
  ni = gaussian_dat$ni,
  tau2 = NA,
  tau2_var = TRUE,
  tau2_param = "log_tau2",
  family = gaussian(link = "identity"),
  fast = TRUE
)
summary(random_fit_log)

Formula and contrasts

The model matrix is the source of truth for coefficient names. This remains true for no-intercept models, factors, and interactions.