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

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

## Why binary outcomes?

Many of the original MAIHDA applications target a **binary** health outcome --
obese vs. not, hypertensive vs. not, screened vs. not. Merlo's framing of
*discriminatory accuracy* is, at heart, a question about a binary classifier: how
well do the intersectional strata alone separate people who do and do not have the
outcome?

For a binary outcome the model is a multilevel **logistic** regression:

\[
\text{logit}\,\Pr(y_{ij} = 1) = \beta_0 + u_j, \qquad u_j \sim N(0, \sigma_u^2),
\]

where \(j\) indexes the intersectional stratum. MAIHDA reads off two
complementary quantities from this model:

- the **VPC/ICC**, the share of the (latent) variation that lies *between* strata;
- the **discriminatory accuracy** (e.g. the AUC / C-statistic), how well stratum
  membership predicts the individual outcome.

A high between-stratum VPC can still go with only *moderate* discriminatory
accuracy at the individual level -- that contrast is the whole point of the "DA"
in MAIHDA, and this vignette shows how to get both numbers.

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

# A two-level outcome: obese (Yes) vs. not (No)
table(maihda_health_data$Obese)
```

## Fitting a logistic MAIHDA model

`fit_maihda()` detects a two-level outcome automatically. If you leave `family` at
its default, the function switches to `family = "binomial"`, recodes the response
to 0/1 for `glmer()`, and warns you so the change is never silent:

```{r autodetect}
model_null <- fit_maihda(
  Obese ~ 1 + (1 | Gender:Race:Education),
  data = maihda_health_data
)
```

The warning is a feature, not a problem. To be explicit (and to silence the
warning), pass `family = "binomial"` yourself -- this is the recommended form for
scripts:

```{r explicit, eval=FALSE}
model_null <- fit_maihda(
  Obese ~ 1 + (1 | Gender:Race:Education),
  data   = maihda_health_data,
  family = "binomial"
)
```

> If you actually want a linear probability model on a 0/1 outcome, ask for it
> explicitly with `family = "gaussian"`; the auto-switch only fires for the
> default family.

## The VPC is on the latent scale

`summary()` reports the VPC/ICC for the logistic model. There is no observed-scale
residual variance for a Bernoulli outcome, so MAIHDA uses the standard
**latent-scale** approximation: the level-1 (within-stratum) variance is fixed at
\(\pi^2/3 \approx 3.29\) for the logit link (and the corresponding constant for a
probit link). The VPC is therefore the between-stratum share of variation on that
underlying latent scale, not on the probability scale.

```{r summary-null}
summary_null <- summary(model_null)
print(summary_null)
```

Read from this null model, the VPC is the **total** between-stratum share of
latent variation in the odds of obesity. As in the Gaussian case, the stratum
random effects capture the *combined* additive + interaction differences across
strata; they isolate the *pure* intersectional (interaction) component only once
the additive main effects of the strata variables are in the model.

> For a binomial model `summary()` reports the discriminatory accuracy (AUC /
> Median Odds Ratio) automatically, so the printed summary above already carries that
> block -- the [Discriminatory accuracy](#discriminatory-accuracy-auc-and-median-odds-ratio)
> section below explains how to read it and how to obtain the pieces on their own.

For an interpretable **probability-scale** complement to the latent-scale VPC, the
package provides `maihda_vpc_response()`, which estimates the response-scale VPC by
simulation (Goldstein, Browne & Rasbash 2002):

```{r vpc-response}
maihda_vpc_response(model_null, seed = 1)
```

Report it *alongside* -- not instead of -- the latent-scale VPC: it depends on the
overall prevalence, and for adjusted models it is evaluated at the mean covariate
profile, so it is best read from the null model.

## Adjusted model and PCV

A "Model 2" adds individual-level covariates to ask how much of the
between-stratum variation they account for. Here we adjust for age, fit on the
**same analytic sample and strata** as the null model, and compare with
`calculate_pcv()`. PCV compares variances across models, so both must use the same
complete-case sample:

```{r adjusted}
health_complete <- maihda_health_data[complete.cases(
  maihda_health_data[, c("Obese", "Age", "Gender", "Race", "Education")]
), ]

model_null2 <- fit_maihda(
  Obese ~ 1 + (1 | Gender:Race:Education),
  data = health_complete, family = "binomial"
)

# Model 2: adjust for an individual-level covariate (age)
model_adj <- fit_maihda(
  Obese ~ Age + (1 | Gender:Race:Education),
  data = health_complete, family = "binomial"
)

pcv <- calculate_pcv(model_null2, model_adj)
print(pcv)
```

The same caveats as in the continuous case apply, with one extra wrinkle: for
non-Gaussian models the latent residual variance is fixed by the link, so a change
in the between-stratum variance is partly a **rescaling** of the latent scale, not
only "variance explained". Interpret the PCV as a model-dependent, descriptive
change, not a causal decomposition.

`glmer()` estimates these variances by maximum likelihood. With few or small strata
they run low, and lower the more fixed effects a model has, so the PCV runs high --
most of all when the second model adds the strata's own main effects.
`estimation = "REML"` compares REML-type variances instead, the counterpart of what
the default does for a continuous outcome. It needs strata large enough for their
own estimates to be precise (see `?calculate_pcv`):

```{r adjusted-reml}
calculate_pcv(model_null2, model_adj, estimation = "REML")
```

> **A note on adjusting for the strata's own categories.** If instead you add the
> *categorical* main effects that define the strata (e.g.
> `Obese ~ Gender + Race + Education + ...`) to recover the additive-vs-interaction
> split shown in the [introduction](introduction.html), `glmer()` may warn that the
> model "failed to converge" or "is nearly unidentifiable", most often when the model
> also has a covariate far from zero or in large units, such as `Age` in years. lme4
> takes its checks in the units of the covariates and with a tight tolerance, so they
> can fire at the maximum itself. `fit_maihda()` then measures the fit with the fixed
> effects centred and scaled, counts it as converged within a tenth of a standard
> error of the maximum, and `print()` says what lme4 warned of and how far the fit is
> from the maximum. Centring and scaling the covariate, together with
> `control = lme4::glmerControl(optimizer = "bobyqa")`, usually avoids the warning.

## What "interaction" means on the logit scale

Once those main effects *are* in the model, the stratum random effect is the pure
intersectional component, and `maihda_interactions()` flags the strata whose
component is credibly non-zero:

```{r interactions}
model_int <- fit_maihda(
  Obese ~ Gender + Race + Education + (1 | Gender:Race:Education),
  data = health_complete, family = "binomial"
)

maihda_interactions(model_int)
```

Read the numbers on the scale they were computed on. In the Gaussian case the
stratum BLUP \(u_j\) is the part of the stratum's mean outcome attributable to
interaction, in the outcome's own units. Here it is a deviation in **log-odds** --
so what is flagged is *multiplicative*, odds-scale interaction. Evans et al. (2024,
sec. 2.5.1) say it plainly: in a logistic MAIHDA you can "no longer directly
interpret \(u_j\) as the change in mean outcomes (i.e., shift in probabilities)
attributable to interaction effects". The printed output says so too, which is why
the binomial fit above carries a line the Gaussian one does not -- and why it prints
a `prob_diff` column the Gaussian one does not: the same departure expressed in
probability points, so the log-odds figure is never the only number on screen. That
column is a display convenience; `scale = "response"` below is how to get it, with
its interval, as data.

There are in fact **three** quantities in play here, and this analysis reports two of
them. Keeping them apart is the whole of the difficulty.

**1. The multiplicative interaction, \(u_j\)** -- whether a dimension multiplies the
odds by the same factor at every level of the others. That is what was just printed,
and what `flagged` is about.

**2. The same departure in probability points, \(\pi^{B}_{j}\)** -- the units to
report if your reader thinks in risks. More on it in a moment.

**3. The additive (risk-difference) interaction** -- whether a dimension adds the
same number of *percentage points* of risk at every level of the others. Neither 1
nor 2 reports this, and it is generally non-zero even when \(u_j\) is exactly zero.
Take a model with a \(-2\) baseline and \(+0.7\) for each of two dimensions and no
interaction whatsoever:

```{r scale-dependence}
p <- function(a, b) plogis(-2 + 0.7 * a + 0.7 * b)

c(`risk added by B, at A = 0` = p(0, 1) - p(0, 0),
  `risk added by B, at A = 1` = p(1, 1) - p(1, 0))
```

The odds ratio for B is identical at both levels of A -- that is what \(u_j = 0\)
means -- yet B adds nine and a half points of risk in one group and fourteen in the
other. The excess is a property of the logistic curve, which is steeper in the middle
than in the tails, not a finding about the world. So **"no strata flagged" supports
"no credible multiplicative interaction", not "no interaction"**: quantity 3 is
generally non-zero regardless, and it is often the one a policy audience cares about.

*Generally*, not always -- worth knowing, because the exception is a perfectly
ordinary design rather than a curiosity. The risk-difference interaction is positive
below the curve's midpoint and negative above it, so it must pass through zero, and
it does so where the two comparisons straddle the midpoint symmetrically:

```{r scale-dependence-zero}
rd_int <- function(b0, a, b) {
  (plogis(b0 + a + b) - plogis(b0 + a)) - (plogis(b0 + b) - plogis(b0))
}

vapply(c(`-2.0` = -2, `-0.7` = -0.7, `0.0` = 0), rd_int, numeric(1), a = 0.7, b = 0.7)
```

At an intercept of \(-(a + b)/2 = -0.7\) the risks add exactly. That is one
configuration out of many, and nothing in the flags tells you which one you are in --
which is the point: they are not evidence about quantity 3 in either direction.

Quantity 2, then. Evans et al. (2024) define it as the gap between a stratum's total
predicted probability and the probability implied by the additive main effects
alone,

\[
\pi^{B}_{j} = \pi_j - \pi^{A}_{j}, \qquad
\pi_j = \mathrm{logit}^{-1}(x_j'\beta + u_j), \qquad
\pi^{A}_{j} = \mathrm{logit}^{-1}(x_j'\beta),
\]

and rank-plot \(\pi^{B}_{j}\) where the linear case plots \(u_j\). Ask for it with
`scale = "response"`:

```{r interactions-response}
maihda_interactions(model_int, scale = "response")
```

Same strata, same evidence, readable units: a flagged stratum's interaction is now a
difference in **percentage points** of predicted probability rather than a log-odds
departure.

That correspondence is exact, not approximate. Write \(g\) for the map from a
stratum's BLUP to its \(\pi^{B}_{j}\); \(g\) is strictly increasing and
\(g(0) = 0\), so the estimate, both interval endpoints and zero carry across
together. `flagged`, `direction` and the p-values are identical under either scale,
and the response-scale interval is the exact image of the link-scale one -- no
simulation, and no delta-method approximation. It inherits the same conditionality
(fixed effects and variance components held at their point estimates). Two things do
change: `se` is dropped, because the resulting interval is deliberately asymmetric
about the estimate and no single standard error would reproduce it; and the ranking
can differ, since the same log-odds departure buys more probability near
\(\pi = 0.5\) than out in the tail.

A `rope` is then read in the same units, which is often where the argument becomes
concrete -- "are more than two percentage points of this stratum's risk attributable
to its interaction?" is a question a reader can actually answer:

```{r interactions-rope}
table(maihda_interactions(model_int, scale = "response", rope = 0.02)$decision)
```

**If it is quantity 3 you want**, neither scale will give it to you -- `scale =
"response"` re-expresses the multiplicative interaction, it does not replace it with
the risk-difference one. Fit a linear probability model instead: `family = "gaussian"`
on the 0/1 outcome, the option `fit_maihda()` points at when it auto-detects a binary
outcome. There the BLUP *is* the risk-difference interaction by construction, at the
usual costs (predictions outside \([0, 1]\), heteroskedastic residuals).

## Discriminatory accuracy (AUC and Median Odds Ratio)

The VPC summarises *variation*; discriminatory accuracy summarises *prediction*.
For a binomial model this is reported **automatically**: `summary()` of a binomial
`maihda_model` carries a `discriminatory_accuracy` slot, and `maihda(...,
family = "binomial")` surfaces it on its summaries and headline `print()`. The
explicit `maihda_discriminatory_accuracy()` below is the same quantity on its own --
it bundles the two individual-level summaries for a logistic MAIHDA model: the
**AUC / C-statistic** (how well the predicted probabilities separate cases from
non-cases) and the **Median Odds Ratio (MOR)** (the between-stratum heterogeneity
expressed on the odds-ratio scale). The strata-only model's AUC is the
discriminatory accuracy of the intersectional strata *themselves* -- Merlo's central
quantity. Comparing it with the adjusted model shows whether individual information
beyond stratum membership sharpens classification:

```{r da}
da_null <- maihda_discriminatory_accuracy(model_null2)
da_adj  <- maihda_discriminatory_accuracy(model_adj)

da_null
da_adj
```

You can also call the pieces directly: `maihda_auc(prob, y)` on any vector of
predicted probabilities and 0/1 outcomes (it equals the Mann-Whitney U statistic,
so it needs no extra package), and `maihda_mor(model)` for the Median Odds Ratio.

```{r auc-direct}
prob_null <- predict_maihda(model_null2, type = "individual", scale = "response")
y_obs     <- as.numeric(lme4::getME(model_null2$model, "y"))

maihda_auc(prob_null, y_obs)
```

An AUC of 0.5 is chance. Here both models sit around 0.6: even a non-trivial
between-stratum VPC translates into only **modest** accuracy at the *individual*
level -- exactly the cautionary message that motivates discriminatory-accuracy
reporting in MAIHDA. Note the adjusted model barely moves the AUC: the
categorical covariates that *define* the strata are already captured by stratum
membership, so only the continuous covariate (`Age`) adds genuinely new
individual-level information.

> **The MOR needs the logit link.** The Median Odds Ratio is an odds-ratio-scale
> quantity, so it is defined only for `binomial(link = "logit")` (the default). For a
> `binomial(link = "probit")` fit, `maihda_mor()` errors and
> `maihda_discriminatory_accuracy()` reports the AUC with `mor = NA` -- the AUC,
> being rank-based, is link-agnostic.

## Plots adapt to the binomial family

The standard plots recognise the binomial family and switch to the probability
scale and to deviance-based diagnostics automatically:

```{r plot-predicted}
# Predicted probabilities per stratum with intervals
plot(model_adj, type = "predicted")
```

```{r plot-vpc}
# Latent-scale variance partition
plot(model_adj, type = "vpc")
```

```{r plot-prediction-deviation, warning = FALSE}
# For binomial fits the dashboard highlights the largest absolute
# deviance residuals rather than raw deviations from the mean.
plot(model_adj, type = "prediction_deviation")
```

See the [plot interpretation vignette](interpreting_plots.html) for how to read
each of these.

## Count outcomes work the same way

A Poisson (count) outcome follows the identical pattern -- pass
`family = "poisson"` to `fit_maihda()`. The VPC then uses the log-link
latent-scale residual variance, and the summary, PCV, and plotting helpers all
behave as above.

One thing to watch with counts: that latent-scale residual variance is an
*approximation*, and which one you use matters when the outcome is rare.
`fit_maihda(count_approximation = )` offers the three of Nakagawa, Johnson &
Schielzeth (2017) -- `"lognormal"` (the default), `"delta"` and `"trigamma"`.
They agree above a mean count of about 2 and diverge sharply below it: at a mean
count of 0.34 they give level-1 variances of 1.37, 2.94 and 9.76, so the same
model reports a VPC six times larger under one than another. `summary()`
therefore names the approximation and the mean count it was evaluated at, and
warns when that count is at or below 2. Report which one you used.

The VPC can also be computed on the count scale itself. For a Poisson or
negative-binomial model `maihda_vpc_response()` evaluates the exact expressions
of Leckie et al. (2020): the share of the variance of the observed counts that
lies between strata, with no level-1 approximation to choose.
`summary(response_vpc = TRUE)` prints it below the latent-scale VPC. The two
answer different questions, so say which one you report.

## References

- Merlo, J. (2018). Multilevel analysis of individual heterogeneity and discriminatory accuracy (MAIHDA) within an intersectional framework. *Social Science & Medicine*, 203, 74-80.

- Evans, C. R., Williams, D. R., Onnela, J. P., & Subramanian, S. V. (2018). A multilevel approach to modeling health inequalities at the intersection of multiple social identities. *Social Science & Medicine*, 203, 64-73.

- 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. doi:10.1016/j.ssmph.2024.101664

- Merlo, J., Wagner, P., Ghith, N., & Leckie, G. (2016). An original stepwise multilevel logistic regression analysis of discriminatory accuracy: the case of neighbourhoods and health. *PLOS ONE*, 11(4), e0153778.

- Leckie, G., Browne, W. J., Goldstein, H., Merlo, J., & Austin, P. C. (2020). Partitioning variation in multilevel models for count data. *Psychological Methods*, 25(6), 787-801. doi:10.1037/met0000265