## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 6, fig.height = 4)
set.seed(20260821)
library(metaGLMM)

## ----data---------------------------------------------------------------------
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

## ----basic-random-fit---------------------------------------------------------
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)
coef(gaussian_fit)
confint(gaussian_fit, method = "wald")
stopifnot(is.finite(gaussian_fit$tau), gaussian_fit$tau > 0)

## ----basic-forest-------------------------------------------------------------
gaussian_studies <- as_metafor_data(
  gaussian_fit, labels = gaussian_dat$study
)
head(gaussian_studies)

forest(
  gaussian_fit,
  labels = gaussian_dat$study,
  xlab = "Effect estimate",
  ci_methods = c("Wald", "profile", "SBC")
)

## ----fixed-fit----------------------------------------------------------------
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)

## ----random-fit---------------------------------------------------------------
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)
random_fit$tau2
random_fit$tau
logLik(random_fit)
nobs(random_fit)
stopifnot(is.finite(random_fit$tau), random_fit$tau > 0)

## ----log-parameterization, eval=FALSE-----------------------------------------
# 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-contract---------------------------------------------------------
gaussian_dat$design <- factor(c("A", "A", "B", "B", "A", "B", "A", "B"))
interaction_fit <- metaGLMM(
  estimate ~ 0 + design + moderator:design,
  data = gaussian_dat,
  vi = gaussian_dat$vi,
  ni = gaussian_dat$ni,
  tau2 = 0,
  tau2_var = FALSE,
  family = gaussian(link = "identity"),
  fast = TRUE
)

colnames(model.matrix(interaction_fit))
names(coef(interaction_fit))

