---
title: "Bayesian MAIHDA for sparse intersections"
author: "Hamid Bulut"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Bayesian MAIHDA for sparse intersections}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.2,
  warning = FALSE,
  message = FALSE
)
library(MAIHDA)
library(ggplot2)
# brms needs Stan and a compiler and takes minutes, so it cannot run on CRAN's or
# pkgdown's builders. Its results are precomputed in data-raw/precompute_sparse_vignette.R
# and read from a small cache here; the fast lme4 fits run live.
pc <- readRDS("sparse_precomputed.rds")
pct <- function(x) sprintf("%.0f%%", 100 * x)
```

## Sparse intersectional cells

Intersectional MAIHDA decomposes the between-stratum variation into an additive
part (the constituent dimensions' main effects) and an interaction part (the
model-implied departure from additivity). That conclusion is defensible only
when each intersectional stratum is populated well enough to estimate its
effect. In practice it often is not. Crossing several dimensions multiplies the strata while the
sample size stays fixed, so cell counts fall away quickly. Estimating an interaction variance from
cells like these is hard.

## A dataset with a known interaction

`maihda_sparse_data` is simulated for exactly that: 4 dimensions form 36 intersectional strata,
and the interaction accounts for **`r pct(pc$truth$gaussian$interaction_share)` of the
variance between the 36 stratum effects** -- on a Gaussian outcome `y` and on the latent
scale of a binary outcome `event` (`?maihda_sparse_data` gives the split). Individuals are
allocated to strata unevenly, reproducing the skewed cell sizes of real survey data.

```{r data}
data(maihda_sparse_data)
d <- maihda_sparse_data

cells <- table(interaction(d$gender, d$ethnicity, d$education, d$age_group, drop = TRUE))
summary(as.numeric(cells))                    # cell-size distribution
sum(cells < 5)                                # how many strata have < 5 people
attr(d, "truth")$gaussian$interaction_share   # the true interaction share: 0.41
```

A third of the strata hold fewer than five individuals. This is not an unusual case, as it is often the situation
once more than two or three dimensions are crossed, and it is
precisely where maximum likelihood becomes unreliable.

## What lme4 reports

We fit the crossed-dimensions decomposition with `engine = "lme4"` on both outcomes.

```{r lme4}
m_g <- maihda(y ~ 1 + (1 | gender:ethnicity:education:age_group),
              data = d, decomposition = "crossed-dimensions", engine = "lme4")
m_b <- maihda(event ~ 1 + (1 | gender:ethnicity:education:age_group),
              data = d, decomposition = "crossed-dimensions",
              engine = "lme4", family = "binomial")

c(gaussian = m_g$decomposition$interaction_share,
  binary   = m_b$decomposition$interaction_share)
# The variance of each dimension and of the interaction
round(rbind(
  gaussian = c(m_g$decomposition$per_dim, interaction = m_g$decomposition$interaction_var),
  binary   = c(m_b$decomposition$per_dim, interaction = m_b$decomposition$interaction_var)
), 3)
```

In both fits maximum likelihood has placed a dimension's variance at zero
(gender's in the Gaussian fit; gender's and ethnicity's in the binary one), and the
interaction variance is positive in both. The interaction share is understated all
the same, at `r pct(pc$lme4$gaussian$share)` for the Gaussian outcome and
`r pct(pc$lme4$binary$share)` for the binary one, against a true
`r pct(pc$truth$gaussian$interaction_share)`.

Only part of that gap is the sparse cells. The share is the interaction variance over
the sum of the four dimension variances and its own, and a dimension's variance,
estimated from two or three levels, describes the population its levels come from:
together the four are larger than the spread of the main effects among the 36 strata.
Fitted to the same stratum effects with 1,000 individuals in every stratum, lme4 puts
the share at `r pct(pc$lme4_large$gaussian$share)` for the Gaussian outcome and
`r pct(pc$lme4_large$binary$share)` for the binary one. The two outcomes' effects differ
only in scale: the gap between these two is that of REML (Gaussian) and maximum
likelihood (binary). The Gaussian `r pct(pc$lme4$gaussian$share)` is close to
`r pct(pc$lme4_large$gaussian$share)`; on the binary outcome the sparse cells make a
substantial interaction look absent.
Neither estimate comes with an interval in this output. Asked for one
(`bootstrap = TRUE`), `maihda()` simulates outcomes from the fitted model and refits it;
the interval for the interaction share then runs from 0% to about 90% on the Gaussian
outcome and over the whole range on the binary one. The sparse cells barely constrain
the split. The draws are simulated from the fit itself, so none of them has a gender
effect.

## What brms adds: posterior uncertainty

Refitting with `engine = "brms"` and a weakly-informative prior on the random-effect
standard deviations changes the output in one important respect: the posterior yields a
credible interval in place of a single number. No standard deviation is reported at
zero, since each is a posterior median, which is above zero whatever the prior; the
`normal(0, 0.5)` prior has its mode at zero and pulls the standard deviations toward it.
The call is shown but not evaluated here (`brms` requires Stan and several minutes; see
the note on precomputation below).

```{r brms-call, eval = FALSE}
set.seed(1)
m_g_brms <- maihda(
  y ~ 1 + (1 | gender:ethnicity:education:age_group),
  data = d, decomposition = "crossed-dimensions", engine = "brms",
  prior = brms::set_prior("normal(0, 0.5)", class = "sd"),
  chains = 4, iter = 2000, warmup = 1000, cores = 4,
  control = list(adapt_delta = 0.97)
)
m_g_brms$decomposition$interaction_share
m_g_brms$decomposition$interaction_share_ci
```

The cached posterior summaries:

```{r brms-summary}
bg <- pc$brms$gaussian; bb <- pc$brms$binary
data.frame(
  outcome      = c("Gaussian", "Binary"),
  true_share   = c(pc$truth$gaussian$interaction_share, pc$truth$binary_latent$interaction_share),
  brms_share   = c(bg$share, bb$share),
  ci_low       = c(bg$ci[1], bb$ci[1]),
  ci_high      = c(bg$ci[2], bb$ci[2]),
  max_rhat     = c(bg$diag$max_rhat, bb$diag$max_rhat),
  divergences  = c(bg$diag$divergences, bb$diag$divergences)
)
```

The posterior median is low as well: `brms` reports the same ratio of variance
parameters, which is not the `r pct(pc$truth$gaussian$interaction_share)` of the drawn
strata, and with two or three levels its value depends on the prior on the dimensions'
standard deviations. The credible interval, however, runs from near zero to past the true value, reflecting
posterior uncertainty about the interaction share. This is the substantive difference.
Maximum likelihood's default output is a single point; `brms` reports an
interval that includes the true share in this simulation and makes the uncertainty
explicit. The diagnostics (`max_rhat` close to 1, no divergent transitions) suggest
the posterior is well-sampled.

These summaries also depend on the prior on the random-effect standard deviations,
`normal(0, 0.5)` above and `normal(0, 1)` for the binary fit: each of the four
dimensions has two or three levels, and the variance of such a dimension follows that
prior (see the caveats in `vignette("cross_classified", package = "MAIHDA")`).
Refitted under brms's default prior instead, the interaction share is about 3% for
both outcomes, and the intervals end between 25% and 31% (the upper limit moves with
the sampler's seed), short of the true `r pct(pc$truth$gaussian$interaction_share)`.

## Comparison

```{r fig-share, fig.height = 3.8}
fig <- data.frame(
  outcome = factor(rep(c("Gaussian", "Binary"), each = 2), levels = c("Gaussian", "Binary")),
  method  = rep(c("lme4", "brms (Bayesian)"), 2),
  share   = c(pc$lme4$gaussian$share, pc$brms$gaussian$share,
              pc$lme4$binary$share,   pc$brms$binary$share),
  lo      = c(NA, pc$brms$gaussian$ci[1], NA, pc$brms$binary$ci[1]),
  hi      = c(NA, pc$brms$gaussian$ci[2], NA, pc$brms$binary$ci[2])
)
truth <- pc$truth$gaussian$interaction_share
large <- data.frame(
  outcome = factor(c("Gaussian", "Binary"), levels = c("Gaussian", "Binary")),
  method  = "lme4",
  share   = c(pc$lme4_large$gaussian$share, pc$lme4_large$binary$share)
)

ggplot(fig, aes(method, share, colour = method)) +
  geom_hline(yintercept = truth, linetype = "dashed") +
  geom_pointrange(aes(ymin = lo, ymax = hi), na.rm = TRUE, linewidth = 0.8) +
  geom_point(size = 2.5) +
  geom_point(data = large, shape = 1, size = 3.5) +
  facet_wrap(~ outcome) +
  scale_y_continuous("Interaction share of between-stratum variance",
                     labels = function(x) paste0(round(100 * x), "%")) +
  labs(x = NULL, colour = NULL,
       title = "A singular ML point estimate vs. a brms posterior interval",
       subtitle = paste0("Dashed line: true interaction share (", pct(truth),
                         "); open circle: lme4 with 1,000 per stratum")) +
  theme_minimal() + theme(legend.position = "none")
```

Both point estimates fall below the true value and below lme4's share with ample data
(the open circles). They differ in what is reported alongside the point. The
`brms` credible interval spans the true `r pct(truth)`, whereas `lme4` reports a single value unless
it is bootstrapped. Reading the `lme4` point at face value (an interaction of roughly
3--23%) is the error the singular fit invites.

## Why the interval, not the point, is what changes

The per-stratum effects (the stratum random intercepts) are similar under the two
engines as both apply partial pooling, shrinking thinly-estimated strata toward the
additive prediction. What differs is the treatment of uncertainty. Maximum likelihood
commits to a single variance decomposition and gives no interval for it unless one is
bootstrapped from that decomposition. `brms` propagates the full posterior, so the interaction share is
accompanied by an interval whose width reflects how little the sparse cells constrain it. Summarising
the posterior by its point discards exactly the information that matters: the degree of
uncertainty. When intersectional cells are sparse,
a regularised Bayesian fit with an explicit interval is the more defensible
summary, read with the priors it rests on.

The two-model PCV is the exception. `maihda(engine = "brms")` with the default
`decomposition = "two-model"` fits the null and the adjusted model separately, so
there is no joint posterior to take a PCV interval from. The PCV it reports is a point
estimate, one minus the ratio of the two models' posterior-median between-stratum
variances, and with few strata it is understated -- the reverse of the
maximum-likelihood PCV, which is overstated (see `?calculate_pcv`).