---
title: "Gaussian aggregate outcomes"
author: "metaGLMM authors"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Gaussian aggregate outcomes}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 6, fig.height = 4)
set.seed(20260821)
library(metaGLMM)
```

## 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`.

```{r 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-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()`.

```{r 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)
```

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.

```{r 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-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.

```{r 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)
```

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.

```{r 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)
```

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.

```{r 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 and contrasts

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

```{r 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))
```