# MAIHDA 0.2.1

## New Features

* **`maihda()` and `compare_maihda_groups()` now accept `group` and `context` together (stratified × contextual MAIHDA).** Supplying both now runs the stratified comparison where **each per-group fit is itself a contextual cross-classified model** — `context` is forwarded into every per-group `fit_maihda()` call.

* **`plot(type = "predicted")` gains an `order_by` argument and now orders the strata by predicted value by default (a ranked caterpillar plot).** The predicted-values plot previously drew the strata in native stratum order, while `maihda_table()$strata` already ranked them by predicted outcome.

* **`plot(<maihda_analysis>, type = "vpc")` gains a `model` argument to show the null VPC, the adjusted VPC, or both as one change plot.** A `maihda()` analysis carries both a null and an adjusted model, but `plot(a, type = "vpc")` only ever drew the null-model VPC, with no way to see the adjusted partition and no cue as to which model was shown.

* **`stepwise_pcv()` gains a `context` argument for a contextual cross-classified stepwise PCV.**

## API Changes

* **`calculate_pvc()` has been renamed to `calculate_pcv()`.** `calculate_pvc()` remains as a deprecated alias that warns and forwards, and will be removed in a future release. The result object carries the estimate in `$pcv` and is classed `c("pcv_result", "pvc_result")`; the old `$pvc` element is kept as a deprecated duplicate, so existing code and objects saved by earlier versions keep working (including `print()`). All internal consumers (`maihda()`, `compare_maihda_groups()`, `maihda_table()`, the tidiers, and the Shiny app), the documentation, and the vignettes now use the corrected name.

## Bug Fixes

* The brms engine's linear-predictor reads were broken under current brms (>= ~2.16): `posterior_linpred()` ignores a `summary = TRUE` argument and returns the raw draws matrix, so several brms paths errored.

* A single-row `newdata` errored for response-scale predictions from a brms cumulative (ordinal) fit.

* Longitudinal (growth-curve) fits on a time axis that does not start at 0 could silently converge to a wrong solution; the growth terms are now fit on internally centered time.

* The Poisson / negative-binomial latent-scale level-1 variance was evaluated at the conditional fitted means (BLUPs included); it is now evaluated at the marginal expected count, matching the cited references.

* The longitudinal PCV compared REML variances across models with different fixed effects; it is now computed from ML-refitted growth models.

* Stratum display labels are no longer whitespace-padded for mixed-type dimensions.

* The longitudinal slope PCV read the wrong covariance cell for `time_degree >= 2`

## Improvements

* `maihda()` now warns when a numeric stratum-defining dimension would enter the adjusted model as a linear fixed effect instead of categorical main effects.

* calculate_pcv(bootstrap = TRUE)` on a non-lme4 engine now explains what is available.

* Documented the latent-scale rescaling caveat on the PCV.

## Documentation

* Fixed a mislabelled figure in the reporting vignette: the `tidy()` caterpillar plotted the null model's stratum estimates, the total between-stratum deviations, but the text called them "interaction estimates" (the pure interactions require `tidy(a, which = "adjusted")`; the vignette now says so).
* The plot-interpretation vignette no longer claims a bare `fit_maihda()` fit is "equivalent" to the `maihda()` analysis for its plots: `plot()` on the analysis routes the VPC/predicted/shrinkage views to the null model and the effect-decomposition views to the adjusted model, which a single fit cannot do.
* The group-comparison vignette now explains its own `group_pcv` figure: with the real PISA data every country's PCV is essentially 1 (the gender-by-SES gaps are purely additive, so the per-country adjusted fits are singular-boundary), and the flat row of 100% bars is itself the substantive result, not a rendering problem.
* Restored the accidentally truncated opening paragraph of the "Planning a MAIHDA analysis" vignette.

# MAIHDA 0.2.0

## New Features

* **Removed the `"risk_vs_effect"` plot type (and the internal `plot_risk_vs_effect()`).** The quadrant view of each stratum's mean prediction against its stratum random effect largely duplicated what `"effect_decomp"` and `"predicted"` already show more directly, so it has been dropped from `plot()` for both `maihda_model` and `maihda_analysis` objects, from `type = "all"`, and from the Shiny app's plot menu.
* Removed the `"ternary"` plot type and the `compute_maihda_ternary_data()`, `maihda_ternary_plot()`, and `plot_maihda_ternary()` functions. The ternary diagnostic normalised each stratum's additive, intersection-specific, and uncertainty signals to sum to 1, which discards effect magnitude and puts two effect components and an estimation-error term on a single triangle. Its content is covered more directly by `"effect_decomp"` (the additive-vs-intersectional split with magnitudes intact) and the `"predicted"` / `"upset"` views (uncertainty shown as intervals), so it has been dropped from `plot()` for both `maihda_model` and `maihda_analysis` objects, from `type = "all"`, and from the Shiny app, and the optional `ggtern` dependency removed.
* **Added the `"upset"` plot type and `maihda_upset_size()`.** An UpSet-style alternative to `"predicted"` that replaces the long intersectional axis labels with a category matrix (an intersection-size bar, a dot matrix encoding each stratum's level on every dimension, and a predicted-value panel, all sharing one column order). Binary 0/1 dimensions collapse to a present/absent row; multi-level factors get one row per level. A `quantity` argument switches the bottom panel between the predicted value (fixed + random) and the stratum random effect (interaction). `maihda_upset_size()` returns content-scaled figure dimensions.
* **`theme_maihda()` now is set as standard theme for ggplot objects**
* **The `predicted` (and longitudinal `trajectories`) plot can now keep the most extreme strata when it truncates, instead of the first by stratum order.** When there are more strata than the `n_strata` cap, the view dropped strata in stratum order, effectively arbitrary with respect to how far a subgroup sits from the population mean, so the most striking strata could be the ones omitted. `plot()` gains an opt-in `select` argument: `"order"` (default, unchanged) keeps the first n_strata in stratum order; `"deviation"` keeps the n_strata furthest from the reference line (largest `|predicted - reference|`, so the most extreme strata in both  directions, not just one tail). For a longitudinal `trajectories` plot it keeps the strata whose trajectories swing furthest from the population curve (peak `|random deviation|` over the time grid). Selection and display are kept separate: `select` changes which strata appear, but the x-axis stays in stratum order. It composes with the flag-aware cap (flagged strata are always kept; `select` governs the fill) and the caption names the rule used ("the 12 strata furthest from the reference, of 200").
* **The BLUP plots can now show only the flagged strata, and never hide a flagged stratum behind the display cap.** When there are many strata the `"predicted"` view caps the number drawn (`n_strata`, default 50) and kept them in stratum order, so a stratum flagged by `maihda_interactions()` that fell past the cap was dropped from the figure entirely; the highlight could not reach it. `plot()` gains `only_flagged`: `TRUE` restricts the `"predicted"` and `"obs_vs_shrunken"` views to the strata carrying a credibly non-zero interaction (and turns the highlight on with the stored diagnostic if it was off), with a caption naming the screen (e.g. "Showing the 7 flagged strata (95% interval, BH-adjusted) of 200 total") and a graceful captioned empty panel when none are flagged. Independently, whenever interactions are highlighted the `n_strata` cap on `"predicted"` is now flag-aware: every flagged stratum is kept and the remaining slots are filled in stratum order, so the signal the highlight exists to surface is never silently truncated away. The across-strata reference line is still computed from the full set, so filtering never shifts it. `"effect_decomp"` deliberately ignores `only_flagged` (its waterfall exists to show each flagged stratum's place in the full distribution) and says so; the framing stays on "flagged", not "significant", consistent with the diagnostic's exploratory, partial-pooling reading.
* **The interaction diagnostic now defaults to FDR control and gains an equivalence (ROPE) reading.** `maihda_interactions()` now defaults to `adjust = "BH"` (Benjamini-Hochberg) rather than `"none"`: fitting and flagging every stratum in one call is a screening question, so the flags should control the false-discovery rate by default. Pass `adjust = "none"` (or `interactions = "none"` to `maihda()`/`fit_maihda()`) for the uncorrected, per-stratum individual-testing view. A new `rope` argument adds an equivalence / smallest-interaction-of-interest reading (Schuirmann 1987; Kruschke 2018): supply a region of practical equivalence (a half-width `d` for `c(-d, d)`, or `c(lower, upper)`, on the link scale) and each stratum gains a `decision` of `"relevant"` (interval entirely outside the region), `"negligible"` (entirely inside), or `"inconclusive"` (straddling a bound), reported by `print()`. The documentation now also keeps two ideas the literature distinguishes apart -- partial pooling regularises magnitude/sign (Gelman, Hill & Yajima 2012) while whether to correct depends on the inferential structure of the claim (Rubin 2021) -- rather than treating shrinkage as itself a multiplicity correction.

## Bug Fixes

* **`compare_maihda(ic = TRUE)` could show information criteria on incompatible scales without a warning.** The comparability check only warned when the models differed in outcome, weights, family/link, analytic sample, or strata, then appended whichever information-criterion columns existed. Two same-family models fitted on different engines, e.g. a Gaussian `lme4` fit (reporting AIC/BIC) and a Gaussian `brms` fit (reporting WAIC/LOOIC) -- agree on all the checked aspects, so no warning fired, yet the appended table placed the likelihood AIC/BIC and the Bayesian WAIC/LOOIC side by side. Those are different scales and are not comparable to each other (as `maihda_ic()`'s documentation already notes). `compare_maihda()` now emits an explicit warning whenever the appended criteria span both the likelihood (AIC/BIC) and the Bayesian (WAIC/LOOIC) scales, so a cross-engine comparison can no longer be read as if the four numbers were on one ruler. An all-`lme4`/all-`ordinal` (AIC/BIC only) or all-`brms` (WAIC/LOOIC only) comparison is unaffected.
* **A fixed interaction among the stratum dimensions (e.g. `Gender * Race`) silently corrupted the decomposition.** R expands `Gender * Race` to `Gender + Race + Gender:Race`, but the workflow only looked for the dimensions' additive main effects when deciding whether the supplied formula was the fully-specified adjusted model. It found both `Gender` and `Race`, treated the fit as the adjusted model, and derived the null by removing only those main effects -- leaving the fixed `Gender:Race` interaction in **both** the null and the adjusted formulas. That fixed cell-means term duplicates the intersectional `(1 | Gender:Race)` random intercept: it absorbs the between-stratum variance into fixed effects, pinning the stratum variance at a singular boundary so the PCV came back `NULL`/invalid (and every per-group PCV came back `NA` in `compare_maihda_groups()` / `maihda(group = )`). `maihda()` and `compare_maihda_groups()` now **reject** such a formula up front with an actionable error -- the MAIHDA adjusted model is defined to carry only the dimensions' additive main effects, with the intersection estimated by the stratum random effect (write `Gender + Race`, and use `decomposition = "crossed-dimensions"` or `maihda_interactions()` to quantify the interaction). `maihda_interactions()` likewise now warns when handed a bare model that carries such a fixed dimension interaction (its stratum random effects are no longer the pure interaction the diagnostic reports). Detection is robust to variable order within the interaction and to non-syntactic names, and a legitimate covariate-by-dimension interaction (e.g. `age * Gender`) is left untouched.
* **Discriminatory accuracy silently dropped the AUC/MOR for a `brms` aggregated-binomial (`y | trials(n)`) fit.** The model layer fits such responses, and `summary()` attaches the discriminatory accuracy automatically for any binomial/Bernoulli fit, but `maihda_discriminatory_accuracy()` only recognised the aggregated structure for `lme4` `cbind(success, failure)` fits. A `brms` `y | trials(n)` fit therefore reached `maihda_auc()` with the per-row success counts* as the response, errored on the non-0/1 values, and the summary quietly omitted the DA. It now detects a `brms` aggregated binomial via the existing trial-count extraction path (`maihda_brms_trial_counts()`, which parses the `trials()` addition term off the stored formula -- `brms` exposes no `weights.brmsfit`) and computes the same count-weighted C-statistic the `lme4` `cbind()` path uses. The rows are ranked by the per-trial probability `predict_maihda()` returns (now normalised to a probability on both engines -- see below), so the AUC ranks by probability and not by a count that confounds the probability with the number of trials; `n_case` / `n_control` are the total successes / failures. Bernoulli and `lme4` fits are unaffected.
* `calculate_pvc()` (and therefore `maihda()`, `stepwise_pcv()`, and `compare_maihda_groups()`) could **silently report a REML-based PCV for a singular multi-random-effect model**. The ML refit (`lme4::refitML()`, needed because REML between-stratum variances are not comparable across models with different fixed effects) was skipped whenever the fit was  globally singular, but a fit can be singular because a **non-stratum** random effect (e.g. an extra `(1 | site)` grouping factor) is on the boundary while the stratum variance is comfortably nonzero. In that case the package returned the REML PCV instead of the intended ML PCV. The ML refit is now skipped only when the **stratum** variance itself is on the boundary (effectively zero, judged on `lme4::isSingular()`'s own relative tolerance), or when `refitML()` fails; a non-stratum boundary no longer suppresses it. The boundary skip for a zero-variance stratum is preserved, so the zero-variance guard in `calculate_pvc()` is not masked.

* **Aggregated-binomial stratum predictions were row-weighted instead of trial-weighted.** For a `cbind(success, failure)` (or `y | trials(n)`) fit the rows of an intersectional stratum can carry very different numbers of binomial trials, but the per-stratum  prediction means averaged the fitted probabilities / linear predictors with unit weights, so a 5-trial row counted as much as a 200-trial row. The unit weighting is correct for the  observed  stratum summaries (which sum successes over summed trials, the trials are already in the denominator), but wrong for the model predictions. Prediction aggregation now uses a dedicated `maihda_prediction_weights()` that weights each row by its binomial trial count -- read from `stats::weights(model, type = "prior")` for an `lme4` `cbind()` fit, or parsed from the formula's `trials()` addition term for a `brms` `y | trials(n)` fit (which exposes `model.frame.brmsfit` but no `weights.brmsfit`, so the prior-weights route is unavailable and the counts would otherwise fall back to unit weights) -- while the observed numerator/denominator path is unchanged. This corrects the predicted stratum means and the `w_sum` weights feeding `maihda_table()`, the predicted-strata / effect-decomposition plots, and the prediction-deviation panels; unweighted and non-binomial fits are unaffected (the weights reduce exactly to the previous values).
* **`brms` cumulative (ordinal) stratum predictions were on the wrong response scale.** The stratum-level prediction helper applied the **scalar inverse link** for the `engine = "brms"`, `family = "ordinal"` response scale, returning a single cumulative probability in `[0, 1]` rather than the expected category score `sum_k k * P(Y = k)` in `[1, K]` that the rest of the package documents and computes -- the individual `predict_maihda()` path already collapses the `fitted()` category-probability array to that score, and the `clmm` (`engine = "ordinal"`) stratum path already used it. This silently mis-scaled (and mis-ranked) `maihda_table()` and the stratum plots for a Bayesian cumulative model. The brms stratum helper now maps the latent location to the expected category score with the same shared cumulative helpers (using the posterior-mean thresholds), so the brms and `clmm` cumulative paths agree and match the documented response scale; the link scale is unchanged.
* **`maihda_interactions()` returned a meaningless scalar diagnostic for a longitudinal fit.** A direct call on a longitudinal `maihda()` analysis (or a bare longitudinal model) fell through to the crossed-dimensions branch and reported a single per-stratum BLUP -- which drops the random slope and so describes a cross-section, not the trajectory. It now errors with a pointer to the trajectory tools, matching the automatic-attachment path that already skips longitudinal objects.
* **The Bayesian `pd` column was mislabelled.** `maihda_interactions()` reported `pd = P(BLUP > 0)` under the heading "probability of direction", so a strongly **negative** interaction (whose direction is near-certain) showed `pd` near 0 rather than near 1. `pd` is now the conventional probability of direction `max(P(BLUP > 0), P(BLUP < 0))` (in `[0.5, 1]`, à la `bayestestR::p_direction`), with the sign reported separately in the existing `direction` column.
* **`predict_maihda()` silently accepted an unseen stratum on the `wemix` and `ordinal` paths.** When `newdata` already carried a `stratum` column, the preparation step returned early and skipped the unseen-stratum check, so a misspelled or genuinely new stratum flowed through to the WeMix/`clmm` linear-predictor helpers, which map a missing random effect to 0, yielding a fixed-only prediction that  looked  valid and contradicted both the documented contract (unseen strata are an error, as for `type = "strata"`) and `lme4`'s default behaviour. `predict_maihda()` now rejects an unseen stratum by default for  every  engine and prediction type, whether the stratum is supplied directly or rebuilt from the grouping variables. A new  `allow_new_levels`  argument (default `FALSE`) opts into the previous behaviour  explicitly: for `type = "individual"` it returns a population-average(fixed-effects-only) prediction for unseen strata, dropping the stratum random effect (forwarded as `allow.new.levels` to lme4 and `allow_new_levels` to brms). Stratum-level predictions have no random effect to report for an unseen stratum, so they remain an error regardless.
* **`predict_maihda(scale = "response")` returned expected counts, not probabilities, for a `brms` aggregated-binomial (`y | trials(n)`) fit.** `brms`'s `fitted()` / `posterior_epred()` return the expected number of successes (`trials * p`) for a binomial fit with a `trials()` term, whereas `lme4`'s `type = "response"` for a `cbind(success, failure)` fit returns the per-trial probability. The same `predict_maihda(scale = "response")` call therefore produced two different scales depending on the engine, so anyone ranking or comparing predicted risks across `brms` strata was comparing expected counts (which confound the probability with the per-row number of trials) rather than probabilities. The `brms` response-scale prediction is now normalised by the per-row trial counts (evaluated on the prediction `newdata`), so both engines return a per-trial probability in `[0, 1]`; the discriminatory-accuracy AUC, the Shiny app's absolute-prediction panel, and any downstream ranking now read a true probability. Bernoulli and non-binomial fits are unaffected (no `trials()` term -> no normalisation).
* **`predict_maihda(allow_new_levels = TRUE)` did not give the documented zero-effect prediction for an unseen stratum on the `brms` engine.** The contract is to drop the stratum random effect (treat it as zero) for a stratum the model never saw. The `brms` path merely forwarded `allow_new_levels = TRUE`, but `brms`'s default `sample_new_levels = "uncertainty"` draws a new-stratum effect from the estimated random-effects distribution rather than zeroing it, so the prediction silently carried a random per-call group effect. Unseen-stratum rows are now split off and predicted with an `re_formula` that excludes the stratum term, while seen-stratum rows keep their estimated effect (a blanket exclusion would have dropped the seen strata's effects too). Any other  random effect the unseen row participates in -- a contextual `(1 | school)` intercept from `fit_maihda(context = )`, a longitudinal `(time | id)` growth term -- is  kept , exactly as `lme4`'s `allow.new.levels` zeroes only the unseen level's effect and retains seen ones; for the usual single-stratum model the excluding `re_formula` is `NA`, the fixed-effects-only population average. `lme4` and the `wemix`/`ordinal` paths already behaved this way and are unchanged, so the engines now agree on the same model.
* **The crossed-dimensions decomposition reported `NaN` additive / interaction shares when there was no between-stratum variance.** The shares split the between-strata variance, so for a degenerate fit whose additive  and  interaction variances are both estimated at exactly zero they are `0 / 0 = NaN`, which leaked through `summary()` and the comparison/tidier outputs. `maihda_cc_partition()` now returns `NA` (a defined "undefined") for the shares when `between == 0`, and for the VPC when the total variance is `0`; `print()` shows `NA (no between-strata variance to split)`. Non-degenerate fits are unchanged.
* **A formula `offset()` term was silently dropped from `wemix` and `ordinal` predictions.** `WeMix::mix()` and `ordinal::clmm()` both honour an `offset(.)` term in the model formula when fitting, but the manual linear-predictor helpers (no `predict.clmm` exists, and WeMix's own `predict()` offers no fixed-only/scale control) rebuilt the linear predictor as `X %*% beta (+ u)` from the design matrix only -- and `model.matrix()` never produces a column for an offset term, so the offset was omitted. Every link-scale prediction (and the response-scale value derived from it) was therefore off by exactly the per-row offset: `predict_maihda()` for both engines, and the ordinal `plot_prediction_deviation_panels()` probabilities, which rebuild the `clmm` location independently. The helpers now evaluate the formula's offset on the prediction data (`stats::model.offset()` of the rebuilt model frame) and add it back,  including for an offset-only null model  whose location is otherwise identically zero. Fits with no offset are unaffected.
* **`maihda()` and `compare_maihda_groups()` chose the engine from the raw outcome, before `subset`/weights.** The wrappers pass `engine` explicitly to every sub-fit, so they pre-selected `engine = "ordinal"` from an ordered-factor outcome -- but they checked the  raw  column, while `fit_maihda()` detects the family on the  analytic sample (after `subset` and weight filtering). An ordered 3-level outcome subset to two observed levels is a binary analytic sample: a direct `fit_maihda()` fits it as `binomial`/`lme4`, but the wrappers pinned `engine = "ordinal"` and then errored (`maihda()`) or failed every group (`compare_maihda_groups()`) with an engine/family contradiction. The ordinal pre-check now runs on the same analytic sample `fit_maihda()` uses (the forwarded `subset`/weights are resolved against the data mask first), so the engine choice can no longer contradict the resolved family. Relatedly, `maihda()` now forwards the evaluated  `subset`/`weights` (not the raw expressions) to its derived null/adjusted fits, mirroring `compare_maihda_groups()`: a response-referencing `subset` (e.g. `subset = y %in% c("lo", "mid")`) was otherwise re-evaluated by each derived fit against `$original_data`, whose response has already been recoded to 0/1, selecting zero rows.

* **Auto-binning a numeric stratum dimension could silently overwrite a user column.** When `make_strata()` (or the `(1 | a:b)` shorthand) auto-bins a numeric dimension `v`, the adjusted-model and prediction machinery add an internal factor column named `.maihda_dim_<v>` (referenced by the adjusted / crossed-dimensions formulae and rebuilt for prediction `newdata`). Neither writer checked whether the user's data already carried a column of that name, so an existing `.maihda_dim_<v>` variable was clobbered -- and, because that column then enters the model, the fit or prediction changed silently. The `.maihda_dim_` prefix is now **reserved**: `make_strata()` rejects a pre-existing `.maihda_dim_<v>` column for any dimension it is about to auto-bin, with an actionable rename hint (or set `autobin = FALSE`), mirroring the existing reserved-weight-column guard. The internal name is now generated through one shared helper so the fit-time and predict-time writers cannot drift, and predicting on a fitted model's own data (which legitimately carries the package's copy) is unaffected.
* **`maihda_table()` could show REML variance rows alongside an ML-based PCV without saying so.** `calculate_pvc()` refits a Gaussian `lme4` fit with maximum likelihood before differencing the between-stratum variances (REML variances are not comparable across models with different fixed effects), but the table's *Between-stratum variance* / *SD* and *VPC/ICC* rows are each model's own (REML) estimate -- the quantities `summary()` reports. The two are on different variance bases by design, so `PCV != (var_null - var_adj) / var_null` read off the table, which looks like an inconsistency. `maihda_table()` now attaches a `models_note` (shown by `print()`) explaining this whenever the PCV's variance basis actually differs from the displayed rows; it stays silent for already-ML engines (`glmer`/`brms`/`wemix`/`ordinal`) and for a boundary null where `calculate_pvc()` keeps the REML fit. The displayed numbers are unchanged, so the table still agrees exactly with `summary()`.

## Improvements

* **Console output is now colour-coded** (via `cli`). The print methods across the package -- `maihda()` analyses, model/summary, PCV, discriminatory accuracy, the interaction diagnostic, group comparison, tables, information criteria, and the rest -- bold section titles, accent the headline values, and dim secondary notes, reusing the **plot accent colour** so the console and the figures agree. The palette is deliberately **valence-neutral**: colour marks emphasis (a finding to look at, a neutral conclusion, a de-emphasised note), never good-vs-bad -- so, e.g., a `negligible` equivalence result is shown in a neutral colour, not green. It **degrades to plain text automatically** wherever ANSI is unsupported (knitr/vignettes, `R CMD check`, `testthat`, `NO_COLOR`), so rendered and captured output is byte-for-byte unchanged. `cli` joins `Imports`.
* `compare_maihda_groups()` (and `maihda(group = )`) gain a **`var_between_adjusted_ml`** column: the adjusted model's *actual* between-stratum variance, read straight off the adjusted fit on the maximum-likelihood scale the PCV is computed on. The existing `var_between_adjusted` column is, by design, a derived coherence quantity (`var_between * (1 - pcv)`, on the REML `var_between`/`vpc` scale so the table satisfies `pcv = (var_between - var_between_adjusted) / var_between` exactly) and is **not** the adjusted fit's own variance; the documentation now states this explicitly and `var_between_adjusted_ml` reports the literal adjusted variance for anyone who needs it. The two differ only by the small REML-vs-ML gap in the null variance.

# MAIHDA 0.1.11

## New Features

* **`maihda()` now reports the intersectional interactions by default.** "Which strata genuinely interact" is the scientific payoff of MAIHDA, so it no longer needs a separate call: `maihda()` computes `maihda_interactions()` on the adjusted (or crossed-dimensions) model, stores it as the `interactions` slot, and surfaces a one-line headline in `print()` (flagged count, the strongest stratum, and -- crucially -- the multiplicity stance actually used, so an uncorrected scan is never silently read as corrected). The computation is essentially free (it reads the stratum estimates the summary already holds; no refit), and it is skipped for a longitudinal decomposition (whose interaction is a trajectory, not a scalar). Control it with the new `interactions` argument: `TRUE` (default) computes it with the diagnostic's own default correction; `FALSE` skips it; or pass a `p.adjust` method name to choose the correction. `fit_maihda()` gains the same argument as an **opt-in** (`interactions = FALSE` by default, since a single fit is often a null model where the diagnostic does not apply). `plot(..., highlight_interactions = TRUE)` now reuses the stored diagnostic, so the plot highlights and the printed headline can no longer disagree.
* Added **`broom`-style `tidy()` and `glance()` methods** for `maihda_model`, `maihda_summary`, and `maihda_analysis`, so MAIHDA results drop straight into tidy data for `ggplot2`, `gt`/`flextable` tables, and other downstream tooling. `tidy()` returns the stratum (intersection) random-effect estimates by default (with the human-readable intersectional label), or the variance-components table (`component = "variance"`) or fixed effects (`component = "fixed"`), all in broom's `estimate`/`std.error`/`conf.low`/`conf.high` shape; `tidy()` on a `maihda()` analysis takes `which = c("null", "adjusted")`. `glance()` returns the **MAIHDA headline as a one-row tibble** -- the VPC/ICC (with its bootstrap/posterior interval), and for a `maihda()` analysis the **PCV**, the adjusted-model AUC, the additive/interaction shares, AUC and MOR for a binary outcome, plus `n_strata`/`nobs`/`engine`/`family` -- a row no generic mixed-model tool emits, since the PCV needs the null+adjusted pair. The layout is uniform across the `lme4`, `brms`, `wemix`, and `ordinal` engines. The generics come from the lightweight **`generics`** package (the same ones `broom`/`broom.mixed` re-export), so `tidy()`/`glance()` dispatch whether the user has `broom`, `generics`, or just `MAIHDA` loaded, with no hard `broom` dependency; the raw fixed-effect/per-row tidying `broom.mixed` already does on the underlying fit is intentionally not duplicated.
* Added **longitudinal (growth-curve) MAIHDA** -- the life-course extension of Bell, Evans, Holman & Leckie (2024, *Soc Sci Med* 351:116955) -- via new `id`, `time`, and `time_degree` arguments on `fit_maihda()` and `maihda()`. Supplying `id` (the person/unit measured repeatedly) and `time` (a numeric measurement-time column) fits a **3-level growth curve** -- occasions within individuals within intersectional strata -- with a random intercept *and* slope on time at *both* the individual and stratum levels (`outcome ~ time + (time | id) + (time | stratum)`; the growth slopes are added automatically, so you write only the strata shorthand `(1 | var1:var2)`). The between-stratum variance, and therefore the **VPC, is then a function of time** (`VarS(t) = a(t)' Sigma_s a(t)`): `summary()` reports the baseline (reference-time) VPC, the full time-varying VPC trajectory, and the stratum/individual intercept-slope-covariance blocks, and a new plot type `"vpc_trajectory"` draws the VPC-over-time curve (with `"trajectories"` for the predicted per-stratum mean trajectories; `plot(type = "vpc")`/`"all"` route there for a longitudinal fit). `maihda(decomposition = "longitudinal")` (selected automatically when `id`/`time` are supplied) fits a null and an adjusted growth model -- the adjusted model adding the dimensions' main effects **and their `dim:time` interactions** -- and reports the **PCV separately for the baseline (intercept) and the slope variance**: the additive-vs-multiplicative split of the intersectional *trajectory* inequality (the paper's "partly multiplicative but mostly additive" finding), returned as a `maihda_long_pcv` object with `pcv_intercept`, `pcv_slope`, and a time-specific `pcv_t` (and a `"pcv_trajectory"` plot). All families with a defined level-1 variance are supported (`gaussian`/`binomial`/`poisson`/`negbinomial`, latent scale for non-Gaussian) on `engine = "lme4"` (any degree) and `engine = "brms"` (linear growth; posterior credible bands on the VPC trajectory). `predict_maihda(type = "strata")` returns each stratum's deviation at the baseline (reference) time plus its raw intercept and slope (a stratum is now a *trajectory*; the baseline deviation, evaluated at `ref_time = min(time)`, is the longitudinal analogue of a cross-sectional stratum BLUP and differs from the raw time-0 intercept whenever time is not zero-referenced). The intercept-only guards elsewhere are untouched -- a longitudinal fit is tagged and routed to the time-varying path, while every other model still rejects random slopes -- and `extract_between_variance()`/`calculate_pvc()` reject a longitudinal model with a pointer to the time-varying decomposition. Design-weighted, contextual, `wemix`/`ordinal`, and group-comparison longitudinal models are out of scope (each errors). A new bundled dataset **`maihda_long_data`** (600 persons x 5 waves, gender x ethnicity x education strata, with constructed mostly-additive trajectory differences) demonstrates it.
* Added **ordinal (cumulative) MAIHDA** for ordered-factor outcomes -- the multicategorical extension the MAIHDA tutorials flag as priority future work. `family = "ordinal"` (alias `"cumulative"`; probit via the new exported `maihda_cumulative("probit")` or `brms::cumulative("probit")`) fits a cumulative (proportional-odds) model: a new **`engine = "ordinal"`** wraps `ordinal::clmm()` for the frequentist path (selected automatically for an ordinal family under the default engine, with a message) and `engine = "brms"` uses `brms::cumulative()`. An ordered-factor outcome with 3+ levels under the all-default family/engine selects the model automatically (with a warning; a 2-level ordered factor still takes the binomial auto-detect path), and an unordered factor response is coerced to ordered in its declared level order with a message. The VPC/ICC lives on the latent scale -- level-1 variance pi^2/3 (logit) or 1 (probit), the same latent treatment as binomial models, so cumulative VPCs are comparable to logistic ones -- and `summary()` gains a `thresholds` slot (the cut points that take the intercept's place, with Hessian SEs) printed alongside the location fixed effects. Because `predict.clmm` does not exist, predictions are built from the stored coefficients: the link scale is the latent location eta = x'beta + u, and the response scale is the **expected category score** sum_k k*P(Y = k) (the quantity the plots label "Average Expected Category Score"), with P(Y <= k) = g^-1(alpha_k - eta) differenced by pure, unit-tested helpers shared with the brms path (whose `fitted()` probability array `predict_maihda()` now collapses to the same score). The two-model `maihda()` decomposition and PCV, `stepwise_pcv()`, `compare_maihda_groups()` (engine handshakes mirror the `sampling_weights` precedent), the stratum plots, `obs_vs_shrunken` (observed mean category score), `risk_vs_effect`, `effect_decomp` (exact on the latent scale), and the deviation panels all work; `maihda_mor()` now also accepts cumulative-logit fits, returning the **median cumulative odds ratio**. Explicit limits with targeted errors: logit/probit links only, the canonical single `(1 | stratum)` structure (no `context`/crossed-dimensions -- use brms), no `sampling_weights` on the clmm path, no parametric bootstrap (clmm has no simulate/refit; brms provides credible intervals), AUC/`maihda_vpc_response()` stay binomial-only, and the ternary diagnostic is not yet supported. `ordinal` joins `Suggests`.
* Added the **negative-binomial family** for overdispersed count outcomes: `family = "negbinomial"` on `fit_maihda()` (and therefore `maihda()`, `stepwise_pcv()`, and `compare_maihda_groups()`). The dispersion parameter theta is *estimated* from the data -- `lme4::glmer.nb()` under the default engine, the `shape` parameter under `engine = "brms"` -- and a fixed-theta `MASS::negative.binomial(theta)` family object is also accepted with lme4, honouring the supplied theta. The VPC/ICC uses the lognormal latent-scale level-1 variance `log(1 + 1/mu + 1/theta)` (Nakagawa, Johnson & Schielzeth 2017, *J R Soc Interface*), the negative-binomial analogue of the Stryhn et al. (2006) Poisson approximation already used by the package (it reduces to it as theta grows); for brms the `shape` draws are propagated into the VPC credible interval. The two-model `maihda()` decomposition, `calculate_pvc()`, parametric-bootstrap intervals (via `lme4::refit()`, which holds theta fixed at its estimate -- the interval is conditional on theta, as documented), prediction, and the stratum plots (routed to the count branch) all work; the log link is required, and the wemix engine, `maihda_discriminatory_accuracy()`, and `maihda_vpc_response()` reject the family with targeted messages. Internally, engine-specific family labels are now canonicalised so the family/link comparability checks in `calculate_pvc()` and `compare_maihda()` no longer depend on raw label strings: two fits whose theta is *estimated* (lme4's theta-embedding `"Negative Binomial(<theta>)"`, brms's `shape`) compare equal, since each label otherwise carried its own theta estimate and could never match. A *fixed*, user-specified theta (`MASS::negative.binomial(theta)`) is treated differently -- it is part of the model specification, so it stays in the comparability key and two fits with different fixed thetas are (correctly) rejected as incomparable.
* Added **design-weighted MAIHDA** (survey/sampling weights) via a new `sampling_weights` argument on `fit_maihda()`, `maihda()`, `stepwise_pcv()`, and `compare_maihda_groups()`. Sampling weights are *not* the same thing as lme4's `weights=` (precision weights, which rescale the residual variance), so feeding survey weights to lmer/glmer maximises the wrong objective and gives invalid population estimates; supplying `sampling_weights` with `engine = "lme4"` is therefore an **error**, and supplying it with the default engine switches (with a message) to the new **`engine = "wemix"`**: weighted pseudo-maximum-likelihood via `WeMix::mix()` (Rabe-Hesketh & Skrondal 2006), the estimator built for NAEP/PISA-style survey analysis. The individual weights enter at level 1 unchanged and the level-2 (stratum) weights are 1, because intersectional strata are exhaustive population cells sampled with certainty. The wemix engine supports the canonical MAIHDA structure -- `gaussian(identity)` or `binomial(logit)` with the single `(1 | stratum)` intercept -- and flows through the whole toolkit: `summary()` reports the design-weighted VPC/ICC (latent-scale pi^2/3 level-1 variance for logistic fits, matching the other engines) and design-consistent (sandwich) fixed-effect standard errors; stratum estimates carry analytic conditional SEs (the design-weighted analogue of lme4's `condVar`, reducing to it at unit weights); `predict_maihda()`, the stratum plots (aggregated with the sampling weights, so stratum summaries are population-representative), `calculate_pvc()` (which now also refuses to compare fits with *different* sampling weights), the `maihda()` two-model decomposition, `stepwise_pcv()`, and `compare_maihda_groups()` all work. For a binomial fit, `maihda_discriminatory_accuracy()` computes the **design-weighted AUC** (each observation contributes its weight as case/control mass) and flags it in the print method. Alternatively `engine = "brms"` accepts `sampling_weights` as likelihood weights (`y | weights(w)`, normalized to mean 1), giving a *pseudo-posterior*: point estimates are design-consistent but credible intervals are not design-based (a message says so). Limitations are explicit rather than silent: no parametric bootstrap for wemix (a design-based interval would need replicate weights -- a possible future extension), and no crossed random effects, so `context =` and `decomposition = "crossed-dimensions"` require lme4/brms. Rows with missing or non-positive weights are dropped with a warning. A unit-weight wemix fit reproduces the lme4 ML fit to numerical precision. `WeMix` joins `Suggests`.
* Added the **contextual cross-classified MAIHDA** -- the "cross-classified MAIHDA" of the literature (e.g. patients cross-classified by intersectional stratum and *hospital*, or students by stratum and *school*) -- via a new `context` argument on `fit_maihda()` and `maihda()`. `context = "school"` (one or more column names) enters each context as a crossed random intercept alongside the intersectional stratum effect, `outcome ~ covars + (1 | stratum) + (1 | school)`, and the summaries then partition the unexplained variance into **between-stratum vs. between-context vs. residual**: the headline VPC/ICC becomes the between-stratum share *net of* the context, each context gets its own `Context: <name>` row in the variance-components table, and a new `$context` summary element carries the per-context variances and shares (with parametric-bootstrap intervals for lme4 via `bootstrap = TRUE`, and posterior credible intervals for brms). In `maihda()` the context random intercept is carried by *both* the null and the adjusted model, so the PCV is computed with the context partialled out; `context` also composes with `decomposition = "crossed-dimensions"` (the context variance then enters the single fit and its VPC denominator). A new plot type `"context_vpc"` (on both `maihda_model` and `maihda_analysis`) bars the stratum vs. context variances with their shares, and `plot(type = "vpc")` renders the contextual split automatically. `context` cannot be combined with `group` (a *stratified* per-level comparison vs. a *joint* crossed model are different designs; supplying both errors), a context variable may not be a stratum dimension or appear as a fixed effect, and a context with few levels weakly identifies its variance (consider `engine = "brms"`). A manually written `... + (1 | school)` still fits and is summarised generically as "Other random effects"; only `context =` activates the labelled partition.
* **Renamed** the `decomposition` value `"cross-classified"` to **`"crossed-dimensions"`** in `maihda()` and `compare_maihda_groups()` (and in the Shiny app's decomposition toggle), because that mode crosses the stratum *dimensions'* main effects as random intercepts -- whereas "cross-classified MAIHDA" in the literature means the contextual stratum-by-place model now fitted via `context` (see above). The old value is accepted as a **deprecated alias** that warns and maps to `"crossed-dimensions"`. Note for code that inspects results: a `maihda_analysis` from this mode now has `mode = "crossed-dimensions"` (was `"cross-classified"`), and `compare_maihda_groups()` results carry `attr(, "decomposition") == "crossed-dimensions"`; printed output and plot titles say "crossed-dimensions" accordingly.
* Added `maihda()`, a single high-level entry point that runs the standard two-model MAIHDA workflow in one call. It fits the **null** model (the formula you supply) and the **adjusted** model (the same model plus the additive main effects of the stratum dimensions), summarises the null model's VPC/ICC, and reports the **PCV** -- the proportional change in between-stratum variance from null to adjusted (the additive share of the intersectional inequality). When a `group` is supplied it also runs this decomposition within each group (the `compare_maihda_groups()` result gains a `pcv` column). It returns one consistent `maihda_analysis` object with new `model_adjusted`/`summary_adjusted`/`pcv` slots (alongside the unchanged null-model `model`/`summary`), and `print`, `summary`, and `plot` methods; `plot()` routes the VPC/shrinkage views to the null model and the additive-vs-intersectional views (`effect_decomp`, `risk_vs_effect`, `ternary`) to the adjusted model, and gains `type = "group_pcv"`. `maihda()` is intrinsically a decomposition and has no single-model mode -- use `fit_maihda()` for a single fit. A numeric dimension auto-binned for the strata enters the adjusted model as its reconstructed tertile factor (not a linear term). Because it cannot decompose without an intersection, `maihda()` errors (rather than returning a half-result) when the stratum-defining variables are unidentifiable (e.g. a hand-built `stratum` column) or define only one dimension; the shorthand `(1 | var1:var2)` and `make_strata()` paths both record the dimensions and decompose normally.
* Added the `maihda_country_data` dataset (OECD PISA 2018, accessed via the `learningtower` package): 3,600 students across six countries with gender x socioeconomic-status strata and mathematics-achievement outcomes. It showcases `compare_maihda_groups()` / `maihda(group = "country")`, since intersectional inequality (VPC/ICC) genuinely differs across the countries.
* Added the `maihda_sparse_data` dataset and a new vignette, **"Bayesian MAIHDA for sparse intersections"**, showing why `engine = "brms"` is the safer estimator when intersectional cells are small. The (simulated) data carry a *known* interaction -- 40% of the between-stratum variance, on a Gaussian outcome `y` and a binary outcome `event` -- across 36 strata with deliberately skewed sizes (median 6, half the cells below 5). Under that sparsity the maximum-likelihood (lme4) crossed-dimensions fit is singular and over-shrinks the interaction (to ~23% Gaussian, ~3% binary -- i.e. a real interaction read as "fully additive"), with no uncertainty attached; a weakly-informative prior on the random-effect SDs regularises the variance off the boundary and returns a calibrated credible interval that covers the truth. The vignette also documents the precompute-and-cache pattern for shipping brms results in a build (Stan cannot run on CRAN's / pkgdown's builders).
* Added `maihda_table()`, a one-call **publication-ready results export** that assembles the two canonical MAIHDA write-up deliverables (cf. Evans et al. 2024's tutorial) from a fitted `maihda()` analysis (or a single `fit_maihda()` model): (a) a **model-results table** contrasting the null (Model 1) and adjusted (Model 2) fits -- intercept, between-stratum variance and SD, VPC/ICC, the PCV, and, for a binary outcome, the AUC and Median Odds Ratio -- and (b) a **ranked-strata table** ordering every intersectional stratum by its model-predicted outcome (their Table 4), with the predicted value's conditional interval, the stratum size, and the stratum random effect. It introduces no new estimator -- the model-results table reuses the quantities from `summary()` and the ranked-strata table reuses `plot(type = "predicted")`'s stratum predictions, so the table agrees exactly with `summary()`/`plot()`. The `$models` data frame is numeric and export-ready (statistics in rows, an estimate + `*_lower`/`*_upper` interval columns per model: the VPC bootstrap/posterior interval and the bootstrap PCV interval are carried, other rows are point estimates), and the `print()` method renders the familiar "estimate [low, high]" layout plus the top/bottom strata (`n_strata`). It adapts to every fit type: a crossed-dimensions analysis gets "Additive share"/"Interaction share" rows instead of the PCV, a contextual cross-classified analysis (`context =`) gets a "Context share (VPC)" row, an ordinal fit leaves the intercept row `NA` (its category thresholds are reported by `summary()`, not the table), and `which = "adjusted"` ranks the strata by the adjusted rather than the null model. Works across the lme4, brms, wemix, and ordinal engines.
* `summary()` of a binomial/Bernoulli MAIHDA model -- and therefore `maihda()`, whose summaries it builds -- now reports the **discriminatory accuracy** automatically: the AUC / C-statistic and Median Odds Ratio (the "DA" in MAIHDA), as a new `discriminatory_accuracy` slot shown by the `print()` methods, so the strata's discriminatory accuracy sits alongside the VPC without a separate call. For `maihda()`, the headline `print()` shows the null model's AUC/MOR with the adjusted model's AUC for comparison. The response-scale (probability-scale) VPC is attached on request via `summary(model, response_vpc = TRUE, seed = )` or `maihda(..., response_vpc = TRUE, seed = )` (it is simulation-based, hence opt-in and seeded). Both are computed from the already-fitted model with no refit, and are skipped for non-binomial outcomes and for the cross-classified fit (whose single-stratum between-variance the MOR/response-VPC require is not defined across crossed random effects). The standalone `maihda_discriminatory_accuracy()`, `maihda_auc()`, `maihda_mor()`, and `maihda_vpc_response()` are unchanged.
* Added `maihda_ic()`, which **surfaces relative-fit information criteria** for choosing between model *structures* (different covariate sets, strata definitions, or families) -- the question the VPC/PCV do not answer. It reports `AIC`/`BIC` for the likelihood engines (lme4, ordinal `clmm`) and the Bayesian `WAIC`/`LOOIC` for brms, takes one or more `maihda_model`s (or a `maihda_analysis`, which expands into its null and adjusted rows), and adds a `delta` column (gap from the best model on the primary criterion). Crucially it handles the **REML pitfall**: `lmer` fits Gaussian models by REML, whose AIC/BIC are *not* comparable across models with different fixed effects (the canonical null-vs-adjusted MAIHDA case), so when more than one model is supplied any REML fit is refitted with maximum likelihood via `lme4::refitML()` before the criteria are read -- matching what `anova()` does on `lme4` models -- and the `estimator` column records it. Design-weighted (`wemix`) pseudo-likelihood fits report `NA` (no standard AIC/BIC is defined). `compare_maihda()` now appends these criteria alongside the VPC/ICC by default (`ic = TRUE`); set `ic = FALSE` for the lean VPC-only table.
* `stepwise_pcv()` now also reports the **discriminatory-accuracy trajectory** for a binary (binomial/Bernoulli) outcome, alongside the between-stratum-variance / PCV trajectory it already produced -- the stepwise discriminatory-accuracy analysis of Merlo et al. (2016). Each step gains an `AUC` (that model's C-statistic; step 0 is the strata-only discriminatory accuracy), `Step_AUC` and `Total_AUC` (the **absolute** change in AUC -- delta-AUC -- versus the previous step and versus the null, in contrast to the **proportional** `Step_PCV`/`Total_PCV`), and `MOR` (the Median Odds Ratio, logit link only) column. No extra models are fit: the columns are read off each step's already-fitted model via `maihda_discriminatory_accuracy()`, so a design-weighted stepwise (`sampling_weights`) reports the design-weighted AUC, and the columns are simply absent for non-binary outcomes (the gaussian/poisson/ordinal table is unchanged). A new `print()` method for the `maihda_stepwise` table notes the proportional-vs-absolute distinction.
* Added `maihda_interactions()`, a diagnostic that flags which strata carry a credibly non-zero **intersectional interaction** -- the heart of "where is there genuine intersectionality". It reads each stratum's interaction BLUP (the stratum random effect of the **adjusted** / crossed-dimensions model, i.e. the departure from the additive main-effects prediction) and flags the strata whose effect is credibly different from zero, returning a ranked `maihda_interactions` table (flagged strata first, by interaction magnitude). For the frequentist engines (`lme4`/`wemix`/`ordinal`) it uses the BLUP's conditional standard error with an optional multiplicity correction (`adjust`, default `"none"`; the docs steer to FDR `"BH"` for screening many strata, with the full set of `p.adjust` methods available for a reviewer who needs a specific one); for `brms` it uses the **exact posterior tail** -- a credible interval and the probability of direction `pd = P(BLUP > 0)` -- rather than a normal approximation, and `adjust` is inert (the Bayesian answer is multiplicity-free). Passing a `maihda()` analysis uses the right (adjusted) model automatically; passing a bare null model is caught with a warning, since its stratum effects mix the additive and interaction parts. The help explains why a correction is optional on already-shrunken BLUP estimates (Gelman, Hill & Yajima 2012; Gelman & Carlin 2014).
* `plot()` gains a `highlight_interactions` argument (on both a `maihda_model` and a `maihda()` analysis) that focuses and stars the `maihda_interactions()`-flagged strata on the BLUP-based views (`effect_decomp`, `predicted`, `obs_vs_shrunken`). Pass `TRUE` (flags computed with defaults), a multiple-testing method such as `"BH"`, or a precomputed `maihda_interactions` object to reuse a specific `conf_level`/`adjust`; for an analysis the flags are computed once from the adjusted model and reused across views. In `effect_decomp`, labels follow the selected flagged set, so a BH screen labels only strata that survive the BH adjustment. `FALSE` (default) leaves every plot unchanged.

## Improvements

* `maihda_mor()` now requires the **logit** link, not merely the binomial family. The Median Odds Ratio is an odds-ratio-scale quantity derived from the logistic latent variance, so applying its `exp(sqrt(2 * V_A) * qnorm(0.75))` formula to a `binomial(link = "probit")` fit returned a number that was not on the model's scale. `maihda_mor()` now errors for a non-logit binomial link, and `maihda_discriminatory_accuracy()` reports the (link-agnostic) AUC with `mor = NA` for such fits, recording the link and explaining the `NA` in its `print()` method.
* Clarified that `maihda_vpc_response()` collapses the fixed part to its mean linear predictor before simulating, so for an adjusted (covariate) model the response-scale VPC is a conditional-at-mean estimate (evaluated at the average covariate profile), not one marginalised over the covariate distribution. It is exact for the canonical null/strata-only model; the documentation now states this and recommends reading it from the null model.
* Updated the "MAIHDA for binary outcomes" vignette, which still claimed the package shipped no AUC/MOR helper and defined a local one. It now uses the exported `maihda_discriminatory_accuracy()`, `maihda_auc()`, and `maihda_mor()`, and points to `maihda_vpc_response()` for the response-scale VPC.
* Added the missing `shinycssloaders` dependency to the README's interactive-dashboard list (it is in `DESCRIPTION` Suggests and used by `run_maihda_app()`).
* Plotting is now unified under the base `plot()` generic. `compare_maihda()` and `compute_maihda_ternary_data()` now return classed objects, so `plot()` dispatches automatically:
  * `plot()` on a `compare_maihda()` result (was `plot_comparison()`)
  * `plot()` on a `compare_maihda_groups()` result, with `type = "vpc"`/"components" (was `plot_group_comparison()`)
  * `plot()` on a `compute_maihda_ternary_data()` result (was `plot_maihda_ternary()`)
* The old `plot_comparison()`, `plot_group_comparison()`, and `plot_maihda_ternary()` functions still work but are **deprecated** and emit a one-time warning pointing to `plot()`.
* `fit_maihda()` now records lme4 fit-quality diagnostics (singular fit and convergence warnings) on the model object; `print()` on the model and `summary()` surface a "Fit diagnostics" note so a boundary/non-converged fit -- which makes the VPC/PCV unreliable -- is no longer silent.
* `compare_maihda_groups()` now raises a single aggregated warning naming any group whose lme4 fit was singular or failed to converge (per-group fits on small strata are where this is most likely), so an unreliable per-group VPC -- often pinned at 0 by a boundary fit -- is no longer silent.
* The `"risk_vs_effect"` plot no longer frames the outcome axis as "risk". A higher predicted value is not universally "bad" (it depends on the outcome), so the axis and title now read as a neutral "mean predicted value/probability", with a note that the direction depends on the outcome.
* Clarified the documentation of `plot_prediction_deviation_panels()` to match the implementation: the labelled points use a per-type metric -- deviation from the mean prediction for Gaussian/Poisson (and the ordinal expected-score mode), the absolute deviance residual for binomial, and surprise for the ordinal surprise mode -- and are not statistical outliers or model-misfit "deviants".
* Clarified that the per-group VPC/ICC in `compare_maihda_groups()` is the between-stratum *share* of variance, which can differ across groups because of the residual variance as well as the between-stratum variance. The documentation now points to the `var_between` column for comparing the absolute amount of intersectional variation, and notes that overlap of separate per-group intervals is not a valid test of whether two groups' VPCs differ.
* `plot()` on a `compare_maihda_groups()` result gains `type = "between_variance"` (and `plot()` on a `maihda(group = )` analysis gains `type = "group_between_variance"`), which bars the absolute between-stratum variance by group -- the *magnitude* of intersectional variation that the VPC *share* cannot convey. All group plots now name any groups omitted because their VPC was not estimable instead of dropping them silently.
* Clarified the PCV documentation (`calculate_pvc()`, `stepwise_pcv()`, and the print method): the PCV is a model-dependent change in between-stratum variance and equals variance "explained" only when the second model nests the first; the stepwise PCV is order-dependent and not a variable's unique contribution. The vignette and Shiny app no longer describe PCV as variance causally "explained" by main effects or treat a negative PCV as evidence of hidden structural inequality.
* Corrected the `summary()` VPC/ICC documentation: the denominator is the total unexplained variance, which includes the variance of any additional random effects (not just between-stratum + residual), and a note on the weighted-Gaussian level-1 variance was added.
* Documented which families the MAIHDA variance summaries support (`gaussian("identity")`, binomial/Bernoulli with logit/probit, `poisson("log")`); other families such as `Gamma(link = "log")` will fit but `summary()`/VPC/PCV will stop with a clear "not implemented" error rather than returning an invalid number.
* README clarifications: a note that the CRAN release can lag this repository (so the newest features may require the GitHub development version), and the interactive-dashboard dependency list now includes `future`, `promises`, and `haven` (for SPSS/Stata uploads).
* Completed the PCV wording cleanup in the remaining vignette and Shiny-app text (the prior pass missed several spots), and softened over-strong app labels: the quadrant plot is "Mean Prediction vs. Stratum Effect" (not "Risk vs. Intersectional Effect"), the cumulative-PCV chart is "Change in Between-Stratum Variance" (not "Variance Explained"), and "deviant strata" is now "most extreme strata".
* Removed the stale checked-in rendered vignette HTML (`vignettes/*.html`); these are build artifacts generated from the `.Rmd` sources and had drifted out of sync with the corrected text. They are now git-ignored and added to `.Rbuildignore`, so a locally rendered HTML is never shipped in the tarball and R CMD build regenerates `inst/doc` from the `.Rmd`.
* Fixed an invalid documented URL flagged by `R CMD check --as-cran`: the `maihda_country_data` `@source` linked to `https://www.oecd.org/pisa/data/`, which returns HTTP 403. It now links to the CRAN page of the `learningtower` package (the reproducible access route used to build the dataset), keeping the OECD PISA 2018 attribution in the text.
* The `data-raw/maihda_health_data.R` regeneration script no longer calls `install.packages("NHANES")` as a side effect; like the country-data script it now stops with a clear message asking the developer to install the dependency.

## Performance

* `make_strata()` (and the prediction-time stratum lookup) now matches rows to strata with a single vectorised, collision-safe key match instead of an O(rows x strata x variables) row-by-row scan, so it scales to large survey datasets. Behaviour is unchanged, including the correct handling of values that contain the stratum-label separator.

## Bug Fixes

* `calculate_pvc()` and `compare_maihda()` now apply their analytic-sample identity checks to **design-weighted (`wemix`) fits**. A `WeMixResults` object exposes no `nobs()`/`model.frame()`, so the row-count, row-identity and outcome-fingerprint guards previously fell through to a silent pass: two wemix fits sharing a formula, `n` and strata but fit to **completely different outcome values** compared as if they were the same analytic sample (`calculate_pvc()` returned a PCV, `compare_maihda()` reported VPCs, neither warned). The guards now fall back to the wrapper's stored analytic `$data` (the rows the engine actually fit), so a mismatched outcome is caught -- `calculate_pvc()` errors and `compare_maihda()` warns -- exactly as for the lme4 path. Genuinely matched wemix fits are unaffected.
* The **longitudinal PCV baseline** is now the between-stratum variance at the *observed* baseline time (`ref_time = min(time)`), not the raw time-0 random-intercept variance. `maihda_longitudinal_pcv()` (and the `maihda_long_pcv` print method) evaluated `Sigma[1, 1]` directly, which is the variance at `time = 0` -- correct only when time is centred so the baseline is 0. For a model whose time does not start at zero (e.g. waves 10:12) this reported the PCV of an extrapolated time-0 variance, which is meaningless and can even be negative; the baseline PCV is now `a(t)'Sigma a(t)` evaluated at `ref_time`, matching how `summary()` reports the baseline VPC. The slope-variance PCV is unchanged (it is invariant to where time is zeroed).
* Cross-sectional, single-value-per-stratum summaries are now **refused for longitudinal (growth-curve) models** rather than silently returning a misleading scalar. A growth model's stratum estimand is a *trajectory* (random intercept + slope), so the scalar BLUP that `maihda_table()`'s ranked-strata table and `plot(type = "predicted" / "obs_vs_shrunken" / "risk_vs_effect" / "effect_decomp" / "prediction_deviation" / "ternary")` build (which adds only the random intercept and drops the slope) misrepresents it. These paths now error -- or, for `maihda_table()`, omit the ranking with an explanatory note -- and point to the trajectory tools (`predict(type = "strata")`, `plot(type = "trajectories")`, `plot(type = "vpc_trajectory")`). `summary()`'s stratum-estimates print is relabelled "baseline (intercept) deviations" for a longitudinal fit, with the same pointer.
* `maihda_ic()` now reports the analytic sample size `n` for a **design-weighted (`wemix`) fit** instead of `NA`. `WeMixResults` has no `nobs()` method, so the IC table fell back to `NA`; it now uses `nrow(model$data)` (the rows the fit used), matching the `nobs` that `glance()` already reports for the same model.
* **Corrected the PCV for Gaussian `lmer` models (the REML pitfall).** The proportional change in between-stratum variance compares the stratum variance of a *null* and an *adjusted* model that differ in their fixed effects (the adjusted model adds the dimensions' additive main effects). `lmer` fits Gaussian models by REML, and a REML variance estimate is **not** comparable across models with different fixed effects, so the PCV was biased -- it overstated the adjusted model's residual between-stratum variance and therefore **understated the additive share (the PCV)**. `calculate_pvc()` -- and hence `maihda()`'s PCV, `stepwise_pcv()`, and the per-group PCV in `compare_maihda_groups()` -- now refits any REML `lmer` model with maximum likelihood (`lme4::refitML()`) before reading the variances (and before the parametric bootstrap, so the interval matches), exactly as `maihda_ic()` already does for AIC/BIC and as `anova()` does on `lme4` models. In simulations with a known 60% additive share (40% interaction), the reported PCV moved from ~50% (biased, REML) to ~58--60% (ML), recovering the truth. `glmer` (GLMM) fits, the brms/wemix/ordinal engines, single-model VPC/ICC summaries (which correctly stay REML, being comparison-free), and singular/boundary fits are unaffected. In `compare_maihda_groups()`, `var_between_adjusted` is now reported as `var_between * (1 - pcv)` so it shares the REML scale of the VPC's `var_between` and the table stays internally coherent.
* VPC/ICC is no longer reported for a Gaussian model fit with a non-identity link (e.g. `gaussian(link = "log")`). The residual variance is on the response scale while the between-stratum variance is on the link scale, so `summary()`, `compare_maihda()`, and `calculate_pvc()` now raise a clear error instead of silently returning an invalid variance partition.
* Binary-outcome auto-detection now keys off the analytic model frame -- after applying covariate transformations (e.g. `log(x)`), dropping rows with missing values, dropping rows with a missing prior weight, and applying any `subset=` (including negative/character row indices) -- instead of the raw outcome column. An outcome that is only 0/1 once excluded rows are removed is now correctly fit with `family = "binomial"`.
* Data-masked engine arguments forwarded through `...` (e.g. `weights=`, `subset=`, `offset=`) now work at any nesting depth, including through the `maihda()` and `compare_maihda_groups()` wrappers. Arguments are captured as quosures and resolved with the data mask (`rlang::eval_tidy`), fixing the previous "object not found" / "..1 used in an incorrect context" failures (a direct `fit_maihda()` call worked, but the wrappers did not).
* A binary outcome is now recoded to 0/1 in a way that no longer breaks a `subset=` expression referencing the original response labels (e.g. `subset = y %in% c("no", "yes")`): the subset is evaluated against the original response before recoding.
* `compare_maihda_groups()` now slices forwarded `weights=`/`subset=`/`offset=` to each group's rows before fitting, not just for the row-count guard. An external (non-column) `weights` vector previously failed every group with a length mismatch, and an external `subset` vector could be recycled onto the wrong rows of later groups; both now align correctly per group.
* The Gaussian VPC/ICC now accounts for prior `weights`. With weights the per-observation residual variance is `sigma^2 / w_i`, so the level-1 variance reported is the mean conditional residual variance `mean(sigma^2 / w_i)` rather than a single `sigma^2`; unweighted models are unchanged.
* `plot_effect_decomposition()` now uses the stratum random effect (BLUP) itself as the intersectional component instead of (full prediction - fixed prediction). With additional random effects such as `(1 | site)` the latter wrongly absorbed those effects into the stratum component.
* The stratum-level "surprise" in `plot_prediction_deviation_panels()` (ordinal, surprise mode) is now the average per-observation surprise, `mean(-log(p))` (log loss), instead of `-log(mean(p))`, which could change stratum rankings.
* The Shiny app (`run_maihda_app()`) now also auto-detects a numeric 0/1 outcome and fits it as binomial under the default family, matching the core API, instead of silently fitting a linear probability model.
* `compare_maihda_groups()` now warns when groups end up with different populated strata even under `shared_strata = TRUE`, since their VPCs are then estimated over different stratum support and are not strictly directly comparable.
* `plot_prediction_deviation_panels()` now plots Poisson/count models on the response (expected-count) scale with count labels, rather than routing them through the Gaussian link-scale branch.
* `compare_maihda_groups(min_group_n = ...)` now guards the analytic sample size (the rows the model actually fits) rather than the raw group row count, so a group with enough raw rows but a tiny usable sample is skipped instead of being fit on a handful of observations.
* `predict_maihda(type = "strata")` now respects `newdata`: it returns only the strata present in `newdata` (and errors on a stratum the model never saw, as `type = "individual"` does) instead of always returning every training stratum. With `newdata = NULL` it still returns all strata.
* `compare_maihda_groups()` now counts populated strata for the pre-fit guard on the analytic model frame, not the raw subgroup. A group with two raw strata but only one stratum left after missing-row removal is now cleanly skipped as VPC-undefined instead of failing during fitting with "grouping factors must have > 1 sampled level".
* `n_boot` for bootstrap intervals must now be at least 10 (the minimum number of successful refits an interval requires); an unusably small `n_boot` fails immediately with a clear message instead of only erroring after the bootstrap runs.
* `calculate_pvc()` and `compare_maihda()` now detect differing prior `weights`: previously two models fit to the same rows/outcome/strata but with different weights compared as if equivalent (PCV returned silently, no warning), even though weights change the variance estimates. `calculate_pvc()` now errors and `compare_maihda()` warns; an unweighted fit and an explicit unit-weight fit are still treated as equal.
* When a two-level non-0/1 outcome is recoded to 0/1, the chosen mapping (which level becomes the modeled event = 1) is now reported via a `message()` and stored on the model as `$response_recoding`. Previously the mapping followed alphabetical (character) or declared (factor) level order silently, with no signal at all when `family = "binomial"` was passed explicitly, so the modeled event could be inverted unnoticed. The 0/1 assignment rule is unchanged (it matches base `glm`).
* Stratum-level plot aggregations (`plot(type = "predicted")`, `"risk_vs_effect"`, `"obs_vs_shrunken"`, `"effect_decomp"`) now honour prior `weights`, collapsing per-stratum predictions with a prior-weight-weighted mean (and weighting the reference lines by the summed weights). This makes the plots consistent with the weighted Gaussian VPC for survey/weighted fits; unweighted models are unaffected, and aggregated-binomial (`cbind`) fits, whose `weights(type = "prior")` are the trials, are left unweighted to avoid double-counting.
* The Shiny app's "Compute Bootstrap CIs" control now actually produces the VPC/ICC bootstrap intervals it (and the vignette) advertise. The bootstrap was previously applied only to the PCV, so the headline VPC/ICC was shown as a point estimate with no interval. The interval is now computed in the background worker (keeping the UI responsive) and displayed alongside the VPC/ICC in the Model Summary and Interactive Explorer tabs; the PCV interval in the PCV Results tab is unchanged.
* The Shiny app (`run_maihda_app()`) no longer aborts the whole analysis when the baseline (null) model has zero or negative between-stratum variance (a singular / no-between-variation fit), which makes `calculate_pvc()` error by design. The fitted model, VPC/ICC, summaries, and plots are now shown as usual and the PCV is reported as unavailable (with the underlying reason) instead of failing the entire workflow.
* `compare_maihda_groups()` now accepts a bare family *function* (e.g. `family = stats::gaussian`) -- one of the documented forms ("as in `fit_maihda`"), alongside a family name (`"gaussian"`) and a family object (`gaussian()`). The per-group fits already handled it; only the family metadata recorded on the result did not, so the call fit every group and then failed with "object of type 'closure' is not subsettable". The family function is now resolved to a family object up front, exactly as `fit_maihda()` does.

## Diagnostics

* `brms` fits now record MCMC convergence diagnostics (maximum Rhat and the number of divergent transitions) alongside the engine, surfaced in the "Fit diagnostics" block of `print()`/`summary()` so a non-converged or divergent Bayesian fit is no longer silent.
* Bootstrap VPC/PCV intervals now report Monte Carlo error: `summary()` and `calculate_pvc()` print the number of successful bootstrap draws and the Monte Carlo standard error so the precision of an interval can be judged.


# MAIHDA 0.1.10

## New Features

* Added `compare_maihda_groups()` to compare intersectional inequality (VPC/ICC and between-/within-stratum variance) across levels of a higher-level grouping variable such as country, region, or survey wave. It fits a stratified MAIHDA model per group, by default using shared/global strata so VPCs are directly comparable, with optional per-group bootstrap confidence intervals.
* Added `plot_group_comparison()` to visualize the result either as a VPC-by-group forest plot or as stacked variance-composition bars.

## Bug Fixes

* Fixed parametric-bootstrap confidence intervals for VPC (`summary()`) and PVC (`calculate_pvc()`): failed `refit()` iterations were silently recorded as `0` instead of being dropped, biasing intervals toward zero and suppressing the high-failure-rate warning. Failed iterations are now excluded correctly.
* Corrected the Poisson VPC residual-variance approximation to `log(1 + 1/mu)` (Stryhn et al. 2006); the previous `1/mu` linearization biased the VPC downward for low-count outcomes.
* `plot_prediction_deviation_panels()` no longer draws zero-width "95% CI" bars when the underlying model does not supply standard errors (e.g. `lme4::merMod`); intervals are omitted instead of collapsed.


# MAIHDA 0.1.9

## Bug Fixes

* Clarified the Shiny dashboard PVC/HUD interpretation so negative PVC values are shown as variance unmasking rather than as unexplained interaction variance.
* Fixed the coverage workflow failure-artifact upload configuration.


# MAIHDA 0.1.8

## General Updates & New Features

* Added `plot_prediction_deviation_panels()` function for visualizing predicted values and identifying deviant cases.
* Added `plot_risk_vs_effect()` function to create a quadrant scatterplot comparing overall marginal predicted risk against pure intersectional effects.
* Added `plot_effect_decomposition()` function to visually decompose the total deviation from the overall mean into additive and intersectional components.
* Replaced the redundant "caterpillar" plot with the "predicted" plot in `plot()` and the interactive dashboard.
* Added automatic tertile binning (via an `autobin` parameter) for numeric grouping variables with more than 10 unique values in `make_strata()`.
* Updated the interactive Shiny Dashboard (`run_maihda_app()`) to include the new visualizations and a toggle for auto-binning continuous strata variables.
* Added detection for binomial data. `fit_maihda()` will now automatically detect binomial outcomes and switch to the appropriate family.

## Bug Fixes

* **VPC/ICC Calculation Fix**: Corrected the residual variance estimation for binomial and ordinal models. The package now accurately applies the theoretical level-1 variance ($\pi^2 / 3$ for `"logit"` links and $1$ for `"probit"` links) internally when summarizing models or bootstrapping the variance partition coefficient, avoiding deflated VPC/ICC metrics.

# MAIHDA 0.1.7

## General Updates & New Features

* Added `stepwise_pcv()` function to sequentially estimate proportional change in variance (PCV) by adding predictors one-by-one.
* Added a fully-featured interactive Shiny Dashboard (via `run_maihda_app()`) for visual data exploration, model fitting, and performance visualization.
* Improved bootstrap methods for more efficient confidence interval estimation.
* Added missing documentation block for the `maihda_sim_data` dataset to resolve `R CMD check` warnings.
* Updated test suite setup: `tests/testthat.R` was modified to correctly use `test_check("MAIHDA")` instead of `shinytest2`.
* Added `importFrom(stats, as.formula)` for the `stepwise_pcv` function to prevent undefined warnings.
* Updated `introduction.Rmd` vignette: added standard CRAN installation instructions, and improved text clarity.

# MAIHDA 0.1.0

## Initial Release

* Initial CRAN submission
* Added `make_strata()` function for creating intersectional strata
* Added `fit_maihda()` function for fitting multilevel models with lme4 (default) or brms engines
* Added `summary()` function for variance partition and stratum estimates
* Added `predict_maihda()` function for individual and stratum-level predictions
* Added `plot()` function with three plot types:
  * Caterpillar plots of stratum random effects
  * Variance partition coefficient visualization
  * Observed vs. shrunken estimates comparison
* Added `compare_maihda()` function for comparing models with bootstrap confidence intervals
* Added comprehensive documentation and vignettes
* Added unit tests for core functionality

## Bug Fixes and Improvements

* Enhanced `make_strata()` to properly handle missing values (NA) in input variables:
  * Observations with missing values in any stratum variable are now assigned NA stratum
  * Missing values are no longer included as valid stratum categories
  * Added comprehensive tests for missing value handling
