## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
if (!"package:drmTMB" %in% search()) {
  library(drmTMB)
}

## -----------------------------------------------------------------------------
set.seed(101)

K <- 30
mu_true <- 0.40 # pooled effect
tau_true <- 0.30 # between-study SD

# Known sampling variances: larger studies (smaller vi) and smaller studies.
vi <- runif(K, 0.02, 0.10)

# True study effects scatter around mu_true with SD tau_true.
theta_i <- rnorm(K, mean = mu_true, sd = tau_true)

# Observed effect sizes: each true effect seen with its known sampling error.
yi <- rnorm(K, mean = theta_i, sd = sqrt(vi))

dat <- data.frame(study = factor(seq_len(K)), yi = yi, vi = vi)
head(dat)

## -----------------------------------------------------------------------------
fit <- drmTMB(
  bf(yi ~ 1 + meta_V(V = vi), sigma ~ 1),
  family = gaussian(),
  data = dat
)

summary(fit)

## -----------------------------------------------------------------------------
is_converged(fit)                 # optimizer convergence
is_converged(fit, include_hessian = TRUE) # also requires a positive-definite Hessian

## -----------------------------------------------------------------------------
diagnostics <- check_drm(fit)
diagnostics[, c("check", "status", "value", "message")]

## -----------------------------------------------------------------------------
mu_hat <- coef(fit, "mu")[["(Intercept)"]]
mu_hat

confint(fit, parm = "mu:(Intercept)")[, c("parm", "lower", "upper")]

## -----------------------------------------------------------------------------
tau_hat <- sigma(fit)[1]
c(tau = unname(tau_hat), tau_squared = unname(tau_hat^2))

## -----------------------------------------------------------------------------
w <- 1 / dat$vi
v_typical <- ((K - 1) * sum(w)) / (sum(w)^2 - sum(w^2))
I2 <- tau_hat^2 / (tau_hat^2 + v_typical)

c(
  tau_squared = unname(tau_hat^2),
  typical_v = v_typical,
  I2_percent = unname(100 * I2)
)

## -----------------------------------------------------------------------------
if (requireNamespace("metafor", quietly = TRUE)) {
  rma_fit <- metafor::rma(yi = yi, vi = vi, method = "ML", data = dat)

  comparison <- data.frame(
    quantity = c("pooled mu", "tau^2", "I^2 (%)"),
    drmTMB = c(mu_hat, tau_hat^2, 100 * I2),
    metafor = c(as.numeric(rma_fit$beta), rma_fit$tau2, rma_fit$I2)
  )
  print(comparison, row.names = FALSE, digits = 4)
}

## -----------------------------------------------------------------------------
fit_reml <- drmTMB(
  bf(yi ~ 1 + meta_V(V = vi), sigma ~ 1),
  family = gaussian(),
  data = dat,
  REML = TRUE
)

data.frame(
  estimator = c("ML", "REML"),
  pooled_mu = c(coef(fit, "mu")[[1]], coef(fit_reml, "mu")[[1]]),
  tau = c(sigma(fit)[1], sigma(fit_reml)[1]),
  tau_squared = c(sigma(fit)[1]^2, sigma(fit_reml)[1]^2)
)

## -----------------------------------------------------------------------------
set.seed(202)
dat$dose <- scale(runif(K, 1, 10))[, 1] # a study-level moderator
# Give the effect size a genuine dependence on the moderator.
dat$yi <- dat$yi + 0.25 * dat$dose

fit_mr <- drmTMB(
  bf(yi ~ 1 + dose + meta_V(V = vi), sigma ~ 1),
  family = gaussian(),
  data = dat
)

coef(fit_mr, "mu")

## -----------------------------------------------------------------------------
c(
  residual_tau_no_moderator = unname(sigma(fit)[1]),
  residual_tau_with_moderator = unname(sigma(fit_mr)[1])
)

## ----layered-meta-syntax, eval=FALSE------------------------------------------
# # Study-level location-SD regression (LSS): z_study is constant within study.
# drmTMB(
#   bf(
#     yi ~ x + (1 | study) + meta_V(V = V),
#     sigma ~ z,
#     sd(study) ~ z_study
#   ),
#   family = gaussian(), data = dat
# )
# 
# # Nested effect-level location-SD regression (LSSS): effect is nested in study
# # and has repeated rows.
# drmTMB(
#   bf(
#     yi ~ x + (1 | study) + (1 | effect) + meta_V(V = V),
#     sigma ~ z,
#     sd(study) ~ z_study,
#     sd(effect) ~ z_effect
#   ),
#   family = gaussian(), data = dat
# )

## -----------------------------------------------------------------------------
set.seed(303)
n_dense <- 8
dat_dense <- data.frame(
  yi = 0.25 + 0.10 * seq_len(n_dense) + stats::rnorm(n_dense, sd = 0.04),
  x = seq_len(n_dense)
)
V_dense <- 0.012 * outer(
  seq_len(n_dense), seq_len(n_dense),
  function(i, j) 0.55^abs(i - j)
)

# A useful preflight: V is numeric, n by n, symmetric, and PSD.
stopifnot(
  is.numeric(V_dense),
  identical(dim(V_dense), c(nrow(dat_dense), nrow(dat_dense))),
  isTRUE(all.equal(V_dense, t(V_dense))),
  min(eigen(V_dense, symmetric = TRUE, only.values = TRUE)$values) >= 0
)

fit_dense <- drmTMB(
  bf(yi ~ x + meta_V(V = V_dense), sigma ~ 1),
  family = gaussian(),
  data = dat_dense
)
check_drm(fit_dense)