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

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)
```

## Before you fit

Most of the decisions that make or break a MAIHDA analysis are taken before any
model is fitted: which dimensions define the strata, how many cells that implies,
and whether the sample can populate them. This vignette walks through those
design decisions, with small runnable checks you can do on your own data first to
evaluate the dimensions, the number of strata, and the analytic sample.

```{r lib}
library(MAIHDA)
data("maihda_health_data")
```

## Is MAIHDA the right tool?

MAIHDA is for questions of the form *"how much of the variation in an outcome lies
between people's intersectional social positions, and how much of that is more than
the sum of its parts?"* It is well suited when:

- you have several **categorical** social dimensions (gender, race/ethnicity,
  education, class, ...) whose **joint** categories define the strata;
- the outcome is measured at the **individual** level;
- you have enough individuals to populate the cells (see below).

## The central tradeoff: more dimensions means emptier cells

Strata are the **cross-product** of the dimensions, so cell counts fall off fast as
you add dimensions. `make_strata()` builds the strata and returns a `strata_info`
table of counts you can inspect *before* modelling:

```{r strata2}
s2 <- make_strata(maihda_health_data, vars = c("Gender", "Race"))
nrow(s2$strata_info)                       # number of strata
summary(s2$strata_info$n)              # cell-size distribution
```

Add education and the same sample splits into many more, smaller cells:

```{r strata3}
s3 <- make_strata(maihda_health_data, vars = c("Gender", "Race", "Education"))
nrow(s3$strata_info)
summary(s3$strata_info$n)
sum(s3$strata_info$n < 10)             # how many strata have < 10 people
```

Each extra dimension multiplies the number of strata and divides the people among
them. Small cells are not fatal, (partial pooling shrinkage is exactly what
protects MAIHDA against noisy small strata) but they have consequences (next
section). A useful rule: choose the fewest dimensions that answer your
question, and look at the cell-size distribution before committing.

## One call: `maihda_describe()`

The checks above — and the rest of a pre-model "Table 1" — come bundled in
`maihda_describe()`. It builds the strata with the same machinery as
`fit_maihda()` (identical stratum IDs, labels, and counts), then reports the
total and complete-case analytic sample, the observed vs. expected strata
(empty cells enumerated, small cells flagged), each dimension's distribution, a
family-aware outcome summary per stratum (proportions for a binary outcome,
never a Gaussian mean/SD), missingness, and a set of data-quality warnings:

```{r describe}
desc <- maihda_describe(BMI ~ Age + (1 | Gender:Race:Education),
                        data = maihda_health_data, flag_stratum_n = 20)
desc
```

Every element is an export-ready data frame — `desc$strata` is the per-stratum
Table 1, ready for `write.csv()` or `knitr::kable()`:

```{r describe-strata}
head(desc$strata[order(desc$strata$n),
                 c("label", "n", "n_analytic", "outcome_mean", "small")])
```

The stratum-size view makes the sparsity of the strata space visible at a
glance (the dashed line is `flag_stratum_n`):

```{r describe-plot, fig.alt="Histogram of stratum sizes with the small-stratum threshold marked"}
plot(desc, type = "stratum_size")
```

For a contextual cross-classified design, pass `context =` (as in
`fit_maihda()`) and `maihda_describe()` also tabulates the context units and
flags a weakly identified context (few levels) before you fit. A fitted
`maihda_model` or `maihda()` analysis can be described the same way —
`maihda_describe(model)` — to document the exact analytic sample post hoc.

## What sparse cells do: singular fits

When cells get very small the maximum-likelihood (`lme4`) estimate of the
between-stratum variance can collapse to the boundary ( a singular fit) and
report a VPC of (near) zero with no uncertainty. The package records this and
surfaces it in a "Fit diagnostics" note rather than letting it pass silently:

```{r singular}
over <- fit_maihda(
  BMI ~ 1 + (1 | Gender:Race:Education),
  data = maihda_health_data[1:60, ]       # deliberately too few people per stratum
)
over
```

If you see a singular-fit note, do not read the VPC as a clean zero. The solution is
to collapse dimensions or categories (fewer, larger cells). `engine = "brms"` does
not repair this: it never puts the variance at zero, but with cells this small its
estimate depends more on the prior than on the data, and its posterior interval
shows how little the cells say (the
[Bayesian sparse vignette](bayesian_sparse_maihda.html)). The interval is what `brms`
adds: its PCV is a point estimate without one, and with few strata it is understated
(see `?calculate_pcv`).

## Continuous variables and the analytic sample

- **Keep continuous variables out of the strata.** A continuous variable in the
  grouping term gives one stratum per value. `make_strata()` will auto-bin a numeric
  dimension into tertiles (with a `message()`), but a continuous *covariate* belongs
  in the fixed part of the formula, not the strata. 

## What the summaries can and cannot tell you

| Quantity | Answers | Does *not* answer |
|---|---|---|
| **VPC/ICC** | share of variance between strata | the *amount* of between-stratum variation (a share can rise just because the residual fell) |
| **PCV** | additive share of the between-stratum variance | a causal decomposition; a negative PCV is not proof of hidden inequality |
| **Discriminatory accuracy (AUC/MOR)** | how well strata predict the *individual* outcome | how large the *group* differences are (a high VPC can go with modest AUC) |

The PCV of a binary or Poisson count outcome fitted with `lme4` compares
maximum-likelihood variances. With few or small strata these run low, the adjusted
model's most, so that PCV runs high. `maihda(estimation = "REML")` and
`calculate_pcv(estimation = "REML")` compare REML-type variances instead, which
corrects this when strata are few but not small (see `?calculate_pcv`).

## Which engine, which design?

- **`lme4` (default)** -- fast frequentist fits for adequately-sized cells.
- **`brms`** -- Bayesian; posterior intervals for the VPC and the variance shares.
  With few-level dimensions the variance components of
  a crossed-dimensions fit depend on the prior on the random-effect standard
  deviations, and its realised-strata block depends little on it when each stratum
  has enough data (see `vignette("cross_classified", package = "MAIHDA")`). Its
  two-model PCV has no interval and is understated when strata are few.

For extensions beyond the cross-sectional case, see the
[crossed random effects](cross_classified.html) (dimensions/contexts) and
[longitudinal](longitudinal.html) vignettes.

## A suggested learning path

1. [Introduction to MAIHDA](introduction.html) -- the end-to-end workflow.
2. [Interpreting MAIHDA plots and diagnostics](interpreting_plots.html).
3. [Finding interaction patterns](finding_interactions.html).
4. [Reporting MAIHDA results](reporting_results.html) -- tidy output and tables.
5. Specialised designs: [binary outcomes](binary_outcomes.html),
   [group comparison](group_comparison.html),
   survey weights, and [Bayesian / sparse](bayesian_sparse_maihda.html).

## References

- Evans, C. R., Leckie, G., Subramanian, S. V., Bell, A., & Merlo, J. (2024). A
  tutorial for conducting intersectional multilevel analysis of individual
  heterogeneity and discriminatory accuracy (MAIHDA). *SSM - Population Health*, 26,
  101664.