Package {MAIHDA}

MAIHDA: Multilevel Analysis of Individual Heterogeneity and Discriminatory Accuracy

Description

Tools for Multilevel Analysis of Individual Heterogeneity and Discriminatory Accuracy (MAIHDA) for intersectional inequality research. Methods are described in Merlo (2018) doi:10.1016/j.socscimed.2017.12.026 and Evans et al. (2018) doi:10.1016/j.socscimed.2017.11.011. The package creates intersectional strata, fits multilevel MAIHDA models, estimates variance partition coefficients, proportional change in variance, stratum effects, and discriminatory-accuracy summaries, and provides diagnostic and presentation plots.

Author(s)

Maintainer: Hamid Bulut me@hamidbulut.com

See Also

Useful links:

• https://github.com/hdbt/MAIHDA

• https://hdbt.github.io/MAIHDA/

• Report bugs at https://github.com/hdbt/MAIHDA/issues

Read an element of a PCV result

Description

$ and [[ read a calculate_pcv result as a list. Reading its deprecated pvc element, a copy of pcv, warns; read pcv instead.

Usage

## S3 method for class 'pcv_result'
x$name

## S3 method for class 'pcv_result'
x[[i, ...]]

Arguments

Value

The element.

Subset a MAIHDA group comparison

Description

Indexing method for compare_maihda_groups results that preserves the comparison's metadata attributes (group variable, engine, family, ...) so the print() header stays populated after a row/column subset. Plain [.data.frame keeps the S3 class but drops those attributes, which would otherwise blank the printed header.

Usage

## S3 method for class 'maihda_group_comparison'
x[...]

Arguments

Value

A maihda_group_comparison with the metadata attributes carried over, or a bare vector when the selection drops to a single column.

Subset MAIHDA information criteria

Description

Indexing method for maihda_ic results that preserves the ic_primary metadata attribute (which names the criterion the delta column is computed on) after a row/column subset. Plain [.data.frame keeps the maihda_ic class but drops that attribute, which would leave print.maihda_ic() testing a dropped (NULL) attribute as a scalar condition and erroring after the table prints.

Usage

## S3 method for class 'maihda_ic'
x[...]

Arguments

Value

A maihda_ic with ic_primary carried over, or a bare vector when the selection drops to a single column.

Add Stratum Labels to Estimates

Description

Internal helper function to merge stratum labels into stratum estimates.

Usage

add_stratum_labels(stratum_estimates, strata_info)

Arguments

Value

Data frame with labels merged in

Bootstrap a crossed-dimensions MAIHDA partition (lme4)

Description

Parametric bootstrap (simulate from the fitted model, refit) of the crossed-dimensions VPC and the additive / interaction shares, returning a percentile interval for each via maihda_bootstrap_ci. lme4 only – brms returns posterior credible intervals directly. When ctx_vars names contextual random intercepts, their variance enters each refit's VPC denominator and a context_vpc interval (the contexts' total share) is returned too.

Usage

bootstrap_cc(
  model,
  cc,
  n_boot,
  conf_level,
  ctx_vars = character(0),
  approximation = "lognormal"
)

Arguments

Value

A list with vpc, additive_share, interaction_share (and context_vpc when ctx_vars is non-empty), each a length-2 interval carrying n_ok/mc_se attributes.

Bootstrap a contextual cross-classified MAIHDA partition (lme4)

Description

Parametric bootstrap (simulate from the fitted model, refit) of the between-stratum VPC and the contexts' total share for a contextual cross-classified fit, returning a percentile interval for each via maihda_bootstrap_ci. lme4 only – brms returns posterior credible intervals directly.

Usage

bootstrap_context(
  model,
  ctx_vars,
  n_boot,
  conf_level,
  approximation = "lognormal"
)

Arguments

Value

A list with vpc (between-stratum share) and context_vpc (contexts' total share), each a length-2 interval carrying n_ok/mc_se attributes.

Bootstrap PCV

Description

Internal function to compute bootstrap confidence intervals for PCV.

Usage

bootstrap_pcv(model1, model2, n_boot, conf_level)

Arguments

Value

A vector with lower and upper confidence bounds

Bootstrap VPC/ICC

Description

Internal function to compute bootstrap confidence intervals for VPC.

Usage

bootstrap_vpc(
  model,
  data,
  formula,
  n_boot,
  conf_level,
  approximation = "lognormal"
)

Arguments

Value

A vector with lower and upper confidence bounds

Calculate Proportional Change in Between-Stratum Variance (PCV)

Description

Calculates the proportional change in between-stratum variance (PCV) between two MAIHDA models. The PCV measures how much the between-stratum variance changes when moving from one model to another, and is calculated as: PCV = (Var_model1 - Var_model2) / Var_model1.

Usage

calculate_pcv(
  model1,
  model2,
  bootstrap = FALSE,
  n_boot = 1000,
  conf_level = 0.95,
  estimation = c("fitted", "ML", "REML")
)

Arguments

Details

The PCV is the proportional change in between-stratum variance when moving from model1 to model2: a positive value means model2 has lower between-stratum variance, a negative value means higher. It is the share of model1's between-stratum variance explained by model2 only in the canonical nested case, where model2 adds fixed-effect predictors to model1 on the same outcome, analytic sample and strata. The function does not require nesting, so for non-nested models the PCV is simply a model-dependent difference in variance, not an explained proportion.

REML vs ML (the estimation argument). estimation sets which between-stratum variances are compared:

"fitted" (default)

each model's own estimate – the REML estimate for a Gaussian lmer fit, the variance summary.maihda_model reports.

"ML"

REML lmer fits are refitted by maximum likelihood (refitML) first, and the bootstrap interval is computed on the same basis. Use it to match an analysis fitted by maximum likelihood.

"REML"

REML-type variances. REML lmer fits are used as fitted. For a binomial or Poisson glmer fit whose only random effect is the stratum intercept, the between-stratum variance is re-estimated with the fixed effects integrated out, as REML does for a Gaussian model (see below).

The PCV is a ratio of two variance estimates and uses no likelihood value, so it does not need maximum-likelihood fits (likelihood-ratio tests and information criteria do; see maihda_ic). REML estimates each model's between-stratum variance with little bias. Maximum likelihood makes no allowance for the degrees of freedom spent on the fixed effects, so it underestimates that variance, and by more in the adjusted model than in the null model: the PCV is overstated, most of all with few strata. In a balanced design whose fixed effects are constant within strata, ML multiplies REML's estimate of the variance of the stratum means, \sigma^2_u + \sigma^2_e / n, by (J - p) / J (J strata, p fixed-effect coefficients, neither estimate at zero): 23/24 for the null model of a 2 x 3 x 4 design but 17/24 for the adjusted one. The bootstrap interval inherits the bias, because it simulates from and refits on the same basis.

"fitted" and "ML" differ for Gaussian lmer fits only. glmer and the wemix and ordinal engines estimate variances by maximum likelihood, so under either setting their PCV is overstated in the same way: most with few strata, and more when strata are small or uneven in size, because the shortfall applies to the between-stratum variance plus the sampling variance of a stratum's own estimate. print() notes when the compared variances are maximum-likelihood estimates.

estimation = "REML" corrects that shortfall for binomial and Poisson glmer fits. Each variance minimises the Laplace-approximated deviance evaluated at the joint mode of the fixed and random effects (lme4's nAGQ = 0 deviance) plus the log-determinant of the fixed-effect information, which integrates the fixed effects out; for the logit and log links this is the criterion glmmTMB(REML = TRUE) minimises. The models are not refitted, and the as-fitted values are returned alongside (pcv_fitted, var_model1_fitted, var_model2_fitted).

Use it when strata are few but each holds enough observations for its own estimate to be precise: the PCV then has little bias, though with very few strata it stays imprecise. When the strata are small as well, neither basis is dependable. A variance estimate cannot fall below zero, so with that much noise the REML-type estimate averages too high, most in the adjusted model: the PCV then varies widely from sample to sample and can come out far too low, even negative, where the maximum-likelihood PCV stays too high.

estimation = "REML" is an error for the wemix and ordinal engines, for other glmer families, for a fit with a further random effect such as a context, for an lmer fit made by maximum likelihood, and with bootstrap = TRUE on a glmer fit.

A brms fit is a Bayesian posterior: "ML" and "REML" perform no refit and estimation_used is "posterior". Single-model VPC/ICC summaries always use the model as fitted.

Latent-scale families and rescaling. For families whose level-1 variance is a fixed latent-scale constant – binomial/Bernoulli (\pi^2/3 logit, 1 probit) and the cumulative (ordinal) model – the linear predictor is identified only up to scale. Adding predictors that explain within-stratum (individual-level) variation cannot shrink that fixed level-1 variance; the latent scale stretches instead, inflating the coefficients and the between-stratum variance alike (Bauer 2009; Mood 2010). Part of a null-vs-adjusted change in the between-stratum variance is then rescaling rather than genuinely explained variance, so latent-scale PCVs tend to be understated and can turn negative on this account alone. The canonical MAIHDA adjusted model – which adds the stratum dimensions' main effects, constant within each stratum – is largely unaffected, but the caveat is first-order whenever an added predictor varies within strata (an individual-level covariate, as in the stepwise_pcv steps that add one). The count families' level-1 variance is not a fixed constant, but as with any non-identity link the same non-collapsibility logic applies in attenuated form. Gaussian identity-link PCVs are not subject to this.

When bootstrap = TRUE, the function uses a parametric bootstrap: it simulates new responses from model2 and refits both models with lme4::refit() for each simulated response to obtain confidence intervals for the PCV estimate. For negative-binomial models (glmer.nb) refit() holds the dispersion parameter theta fixed at its original estimate, so the interval is conditional on the estimated theta.

A bootstrap draw whose null-model between-stratum variance lands on the zero boundary has no defined PCV (the denominator is zero); such draws are excluded, so the percentile interval is conditional on estimating a positive null variance. Whenever any draws hit the boundary the function warns, reports the count as n_boot_boundary on the result, and print() repeats the caveat – a sizeable boundary share signals weak between-stratum variation, and the PCV itself is then fragile.

Bootstrap draws whose refit optimiser reports non-convergence are kept in the interval (they still have a defined PCV) and counted in n_boot_nonconverged, which print() shows. When more than half the contributing draws did not converge the interval is still returned, but with a warning and interval_reliable = FALSE; treat it as indicative only and check for singular or failing fits.

The bootstrap is available for the lme4 engine only. For the other engines the PCV is a point estimate: a brms fit's posterior credible interval (reported by summary.maihda_model) covers a single fit's VPC/ICC, not the PCV, which compares two separately fitted models – no posterior interval for the PCV itself is computed – and a design-based (wemix) interval would require replicate weights.

brms fits. For two brms fits the variances compared are the posterior medians of each model's between-stratum variance, the values the variance-components table of summary.maihda_model shows. The posterior of a between-stratum variance is right-skewed when strata are few, and more so in the adjusted model, which spends degrees of freedom on the main effects. The PCV is therefore understated with few strata – the opposite direction to maximum likelihood – and with very few strata it can be negative where the true value is positive. The understatement shrinks as strata are added.

Value

A list containing:

References

Bauer, D. J. (2009). A note on comparing the estimates of models for cluster-correlated or longitudinal data with binary or ordinal outcomes. Psychometrika, 74(1), 97-105.

Mood, C. (2010). Logistic regression: why we cannot do what we think we can do, and what we can do about it. European Sociological Review, 26(1), 67-82.

See Also

stepwise_pcv for the sequential (one-variable-at-a-time) PCV, and maihda which computes the canonical null-vs-adjusted PCV automatically.

Examples

# Create strata and fit two models
strata_result <- make_strata(maihda_sim_data, c("gender", "race"))
model1 <- fit_maihda(health_outcome ~ age + (1 | stratum), data = strata_result$data)
model2 <- fit_maihda(health_outcome ~ age + gender + (1 | stratum), data = strata_result$data)

# Calculate PCV without bootstrap
pcv_result <- calculate_pcv(model1, model2)
print(pcv_result$pcv)

# Calculate PCV with bootstrap CI
# pcv_boot <- calculate_pcv(model1, model2, bootstrap = TRUE, n_boot = 500)
# print(pcv_boot)

Compare MAIHDA Models

Description

Compares variance partition coefficients (VPC/ICC) across multiple MAIHDA models, with optional bootstrap confidence intervals.

Usage

compare_maihda(
  ...,
  model_names = NULL,
  bootstrap = FALSE,
  n_boot = 1000,
  conf_level = 0.95,
  ic = TRUE
)

Arguments

Details

VPCs are only directly comparable when the models share an outcome, family/link, analytic sample, and strata – the canonical use is nested models (e.g. null vs covariate-adjusted) on the same data and strata, to show how the VPC attenuates. If the supplied models differ in any of these, compare_maihda() still returns the table but issues a single warning, because the VPCs are then not directly comparable. The appended information criteria follow the rules of maihda_ic instead: they need the same outcome, analytic sample and weights but not the same strata, and Poisson and negative-binomial fits of the same counts, or binomial or cumulative fits differing in their link, have comparable criteria although their VPCs are not comparable. In addition, when the appended criteria mix scales – likelihood AIC/BIC (lme4/ordinal) shown alongside Bayesian WAIC/LOOIC (brms), which can happen for a same-family lme4-vs-brms comparison that the family/link check does not flag – compare_maihda() warns, because those criteria are on different scales and are not comparable to each other.

Value

A maihda_comparison data frame of VPC/ICC by model. Interval columns (ci_lower/ci_upper) are included when any model supplies an interval – an lme4 bootstrap CI or a brms posterior credible interval. When ic = TRUE, information-criteria columns (AIC/BIC or WAIC/LOOIC, whichever apply) are appended.

Examples

# Canonical use: nested models on the SAME data and strata (null vs adjusted)
strata <- make_strata(maihda_sim_data, vars = c("gender", "race"))

null_model <- fit_maihda(health_outcome ~ 1 + (1 | stratum), data = strata$data)
adj_model  <- fit_maihda(health_outcome ~ age + (1 | stratum), data = strata$data)

# Compare without bootstrap
comparison <- compare_maihda(null_model, adj_model,
                            model_names = c("Null", "Adjusted"))

# Compare with bootstrap CI
comparison_boot <- compare_maihda(null_model, adj_model,
                                 model_names = c("Null", "Adjusted"),
                                 bootstrap = TRUE, n_boot = 500)

Compare MAIHDA Metrics Across Levels of a Grouping Variable

Description

Fits a separate random-intercept MAIHDA model (intercept-only random effects; any fixed-effect covariates in formula are still used) within each level of a higher-level grouping variable (for example country, region, or survey wave) and reports how the variance partition coefficient (VPC/ICC) and the between-/within-stratum variance components differ across those groups. When the strata are defined by at least two dimensions it also fits the adjusted model (the dimensions' additive main effects) within each group and reports the per-group pcv – the proportional change in between-stratum variance, i.e. the additive share of that group's intersectional inequality.

Usage

compare_maihda_groups(
  formula,
  data,
  group,
  engine = "lme4",
  family = "gaussian",
  shared_strata = TRUE,
  min_group_n = 30,
  bootstrap = FALSE,
  n_boot = 1000,
  conf_level = 0.95,
  autobin = TRUE,
  decomposition = c("two-model", "crossed-dimensions"),
  context = NULL,
  sampling_weights = NULL,
  warn_linear = TRUE,
  estimation = c("fitted", "ML"),
  ...
)

Arguments

Details

It estimates one VPC per group as a stratified analysis: each group is modelled independently – the grouping variable itself is not modelled as a crossed random effect, and the strata are not adjusted for it. (A context argument, if supplied, does add a crossed contextual random intercept within each group's model, making each per-group fit contextual cross-classified; the groups themselves remain independently fitted.)

The VPC is the share of the unexplained variance that lies between strata, not the absolute magnitude of intersectional inequality. Because it is a ratio, a group's VPC can differ from another's because the between-stratum variance differs, because the within-stratum (residual) variance differs, or both – two groups with the same between-stratum variance can have very different VPCs. To compare the absolute amount of between-stratum (intersectional) variation across groups, read the returned var_between column alongside the VPC rather than treating a higher VPC as "more inequality".

It is descriptive: it reports each group's VPC (with an interval when available – an lme4 bootstrap CI or a brms credible interval) for side-by-side comparison, but does not test whether the VPCs differ between groups. The per-group intervals describe each group's own uncertainty; whether two intervals overlap is not a valid test of the difference between their VPCs, which would require modelling that difference directly.

Robustness: a group whose analytic sample (rows surviving the model frame) has fewer than min_group_n observations is always skipped with a warning. A group with fewer than two populated strata is also skipped (VPC is undefined with a single stratum) when the stratum membership is known before fitting – that is, when shared_strata = TRUE or data already carries a stratum column. Under shared_strata = FALSE strata are rebuilt inside each group, so a degenerate single-stratum group is instead reported with a "fit failed" status rather than a pre-fit skip. A singular fit yields a VPC of 0 rather than an error (unlike calculate_pcv). A hard fit failure in one group records NA and a status note without aborting the whole comparison; one warning names the groups whose fit failed and gives each distinct reason once.

Fit-quality diagnostics: on every engine, one warning names the groups with a fit problem, and such a group keeps a status of "ok" (the fit did complete). A group is named for a singular fit when its between-stratum variance is estimated at zero, so that var_between and the VPC are 0. A zero in another variance component – a context, a further random effect, the interaction or a dimension of the crossed-dimensions model – leaves them as estimated and does not name the group. A group is also named when its model did not converge and, in a two-model decomposition, when its adjusted model did not. The warning says which: an adjusted model affects the group's PCV, not its VPC, and pcv_status does not record it. A brms model is named when its largest R-hat is above 1.01 or its sampler had a divergent transition, and is never reported as singular.

Value

A data.frame of class maihda_group_comparison with one row per group and columns group, n, n_strata, vpc, var_between, var_other, var_residual, status (and ci_lower/ci_upper whenever a group's summary supplies an interval: an lme4 bootstrap CI with bootstrap = TRUE, or a brms posterior credible interval, which is returned without bootstrapping). When the strata are defined by at least two dimensions, two further columns report the per-group null -> adjusted decomposition: pcv (proportional change in between-stratum variance when the dimensions' additive main effects are added; computed on the variance basis set by estimation – see calculate_pcv), var_between_adjusted (a derived coherence quantity, reported as var_between * (1 - pcv) so it shares the scale of the REML var_between/vpc and the table satisfies pcv = (var_between - var_between_adjusted) / var_between exactly – it is not the adjusted fit's own variance), and var_between_adjusted_ml (the adjusted model's actual between-stratum variance, read straight off the adjusted fit on the same basis as the PCV; under the default estimation = "fitted" this is the REML variance, coinciding with var_between_adjusted, and under estimation = "ML" it is the ML variance, differing from var_between_adjusted only by the small REML-vs-ML gap in the null variance; despite its name it is not always a maximum-likelihood estimate; for a brms fit it is the posterior median of the adjusted variance, which var_between_adjusted equals as well). All three are NA for a group whose adjusted fit failed. A fourth column, pcv_status, records the decomposition outcome per group: "ok" when the PCV was computed, "failed" when the adjusted model or PCV errored (pcv is then NA and the group is named in a warning – the group's own status can still be "ok" because its null VPC model succeeded), and "singular" when the adjusted fit was singular and the PCV saturated near 100% – not warned about, since a singular adjusted fit is indistinguishable from genuinely additive strata (a legitimate, common result), so it is flagged only in this column. These four columns are omitted entirely when the strata have a single dimension. When context is supplied, two further columns report each group's contextual partition: var_context (the between-context variance, summed over contexts) and vpc_context (the contexts' share of the group's unexplained variance); the per-context split is on the "context_per" attribute and the context name(s) on "context_var". These are dropped when no context is supplied. n is the analytic sample size used by the model (after dropping rows with a missing outcome/covariate) for both fitted and skipped groups, falling back to the raw row count only when the model frame cannot be built. var_other is the variance of any additional random effects and is 0 for the canonical single-stratum model. Groups that were skipped or failed have NA metrics and an explanatory status; for a failed fit it is "fit failed: " followed by the first line of the error, whose whole message is in a warning.

See Also

compare_maihda for comparing different models on the same data; plot.maihda_group_comparison for visualising the result.

Examples

data(maihda_country_data)
# How does gender x SES inequality in PISA math scores differ across countries?
cmp <- compare_maihda_groups(
  math ~ 1 + (1 | gender:ses),
  data = maihda_country_data,
  group = "country"
)
print(cmp)
plot(cmp, type = "vpc")

Extract Between-Stratum Variance

Description

Internal function to extract between-stratum variance from a MAIHDA model. For a brms fit it is the posterior median of the variance draws.

Usage

extract_between_variance(model)

Arguments

Value

Numeric value of between-stratum variance

Fit MAIHDA Model

Description

Fits a multilevel model for MAIHDA (Multilevel Analysis of Individual Heterogeneity and Discriminatory Accuracy) using lme4, brms, WeMix (for design-weighted (survey) data), or – for an ordered-factor outcome – a cumulative link mixed model via ordinal::clmm().

Usage

fit_maihda(
  formula,
  data,
  engine = "lme4",
  family = "gaussian",
  autobin = TRUE,
  context = NULL,
  sampling_weights = NULL,
  id = NULL,
  time = NULL,
  time_degree = 1,
  stratum_slope = TRUE,
  interactions = FALSE,
  count_approximation = c("lognormal", "delta", "trigamma"),
  ...
)

Arguments

Value

A maihda_model object containing:

Examples

# Standard approach: manually create strata first
strata_result <- make_strata(maihda_sim_data, vars = c("gender", "race", "education"))
model <- fit_maihda(health_outcome ~ age + (1 | stratum),
                    data = strata_result$data,
                    engine = "lme4")

# Simplified approach: specify stratifying variables directly in the grouping structure
# The function internally calls make_strata() to create intersectionals
model2 <- fit_maihda(health_outcome ~ age + (1 | gender:race:education),
                     data = maihda_sim_data,
                     engine = "lme4")

# Contextual cross-classified MAIHDA: strata crossed with a higher-level
# context (here country) -- the literature's cross-classified MAIHDA.
data(maihda_country_data)
model3 <- fit_maihda(math ~ 1 + (1 | gender:ses),
                     data = maihda_country_data,
                     context = "country")
summary(model3)  # between-stratum vs. between-country vs. residual

Run a Complete MAIHDA Analysis

Description

A single high-level entry point that runs the standard two-model MAIHDA workflow and returns one bundled object. It fits the null model (covariates plus the intersectional random intercept, excluding the stratum dimensions' main effects) and the adjusted model (the null plus the additive main effects of the stratum-defining dimensions), summarises the variance partition (VPC/ICC) of the null model, and reports the PCV – the proportional change in between-stratum variance from the null to the adjusted model, i.e. the additive share of the intersectional inequality. When a higher-level grouping variable is supplied it also compares this decomposition across that variable's levels.

Usage

maihda(
  formula,
  data,
  group = NULL,
  context = NULL,
  engine = "lme4",
  family = "gaussian",
  decomposition = c("two-model", "crossed-dimensions", "longitudinal"),
  autobin = TRUE,
  shared_strata = TRUE,
  min_group_n = 30,
  bootstrap = FALSE,
  n_boot = 1000,
  conf_level = 0.95,
  response_vpc = FALSE,
  seed = NULL,
  sampling_weights = NULL,
  id = NULL,
  time = NULL,
  time_degree = 1,
  stratum_slope = TRUE,
  interactions = TRUE,
  estimation = c("fitted", "ML", "REML"),
  ...
)

Arguments

Details

Binomial companions. For a binary outcome the model summaries also carry the discriminatory accuracy (AUC / C-statistic and Median Odds Ratio) – the "DA" in MAIHDA – automatically, so the null model's strata-only AUC sits alongside its VPC; set response_vpc = TRUE to add the (simulation-based) response-scale VPC as well. These are read from summary(x) and the attached summary_adjusted, and the headline print() shows the null-model AUC.

This is a convenience wrapper around fit_maihda, calculate_pcv, summary.maihda_model and compare_maihda_groups. It always runs a decomposition, never a single plain fit: by default the two-model null/adjusted pair, or – with decomposition = "crossed-dimensions" – one crossed model whose additive dimension random effects carry the additive share. For a single undecomposed fit (e.g. just the null-model VPC / discriminatory accuracy), call fit_maihda directly.

The dimensions' additive main effects. You may write them in the formula – the fully-specified, lme4-native adjusted model outcome ~ covars + var1 + var2 + (1 | var1:var2) – or omit them. Either way the null excludes the dimension main effects and the adjusted includes them: when the formula already lists them it is taken as the adjusted model and the null is derived by dropping them; when they are missing maihda() adds them to the adjusted model and emits a message() so the decomposition stays explicit. Only the additive main effects belong here: a fixed interaction among the stratum dimensions – var1 * var2, which R expands to var1 + var2 + var1:var2 – duplicates the intersectional stratum random intercept (it absorbs the between-stratum variance into fixed cell means, which makes the PCV invalid), so it is rejected with an error. Write var1 + var2; the intersection is estimated by the stratum random effect (and quantified by decomposition = "crossed-dimensions" or maihda_interactions). The dimensions themselves are read from the random term: the shorthand (1 | var1:var2) and make_strata both record them, and a numeric dimension that make_strata() auto-binned enters the adjusted model as its reconstructed tertile factor (matching the strata), not as a linear term. Because maihda() is intrinsically a decomposition, it errors (rather than returning a null-only result) when it cannot build the adjusted model – when the dimensions cannot be recovered (a hand-built stratum column records none) or there is only one dimension (no intersection to decompose). Use fit_maihda for those single-model fits.

Value

An object of class maihda_analysis: a list with

Numeric stratum dimensions

The stratum-defining dimensions are categorical. A factor or character dimension enters the adjusted model as a categorical main effect, and a numeric dimension with more than 10 unique values is tertile-binned when autobin = TRUE (see make_strata). A numeric dimension that is not auto-binned – a category code such as education = 1, 2, 3 (few unique values), or a many-valued numeric fitted with autobin = FALSE – enters the adjusted model as a single linear fixed effect rather than categorical main effects, which changes the PCV/decomposition interpretation. maihda() therefore warns in that case (for three or more distinct values – a binary numeric dimension is exempt, since with two levels a linear term and a factor are equivalent): wrap category codes in factor() before fitting, or, for a genuinely continuous variable, use it as a covariate rather than a stratum dimension.

Realised strata

For a brms crossed-dimensions fit, decomposition$realised and the "Realised strata" block of print() summarise, over the posterior draws, the effects of the J strata in the data. In each draw stratum j has the effect \zeta_j, the sum of its dimensions' level effects and its own (interaction) effect. These are finite-population summaries (Gelman 2005): they describe the strata and levels in the data, where the variance components describe the distributions those were drawn from. A dimension with two levels says a lot about the former and little about the latter, so the realised values depend little on the prior on the standard deviations when each stratum has enough data to pin down its own effect. With few observations per stratum each draw also carries the uncertainty about every stratum's effect: the realised variances then come out too high and move with the prior. Compare the intervals of the stratum estimates with the spread between the strata. Per draw:

between_var

\sum_j (\zeta_j - \bar\zeta)^2 / (J - 1), half the mean squared difference between two distinct strata.

vpc

between_var over itself plus the within-stratum variance: \pi^2/3 for a logit link, 1 for a probit link, \pi^2/6 for a cloglog link, the residual variance of a Gaussian fit. Not given for a Poisson or negative-binomial fit, whose within-stratum variance is computed from the variance components. With a context, the realised variance of the context effects is added to the denominator (see the section Realised strata and contexts).

interaction_var

the residual sum of squares of the least-squares fit of the \zeta_j on an intercept and the dimensions' levels, divided by its residual degrees of freedom, df_interaction. Strata missing from the grid of dimension levels are allowed.

additive_share

1 minus interaction_var / between_var, and interaction_share the rest: a mean-square split, like a PCV computed from REML variances. It is near zero when the dimensions explain none of the differences between the strata, and can then fall below zero.

per_dim

the variance of each dimension's level effects, with its number of levels minus one as divisor.

mor

for a logit link, the median odds ratio between two distinct strata of the data: the exponential of the median of |\zeta_j - \zeta_{j'}| over the pairs. With many strata it is computed from an evenly spaced subset of at least 1000 draws.

Each is the posterior median, with the equal-tailed credible interval at conf_level in the matching *_ci element. Every stratum counts once, whatever its size. n_strata is J. An lme4 crossed-dimensions fit reports its variance components only.

Realised strata and contexts

For a brms fit with a context, summary()$context$realised and the "Realised strata and contexts" block of print() give the realised partition: over the posterior draws, the finite-population variances (Gelman 2005) of the effects of the strata and of each context's levels in the data. A context's variance component describes the distribution its levels were drawn from, and with few levels it follows the prior on its standard deviation; the variance of the effects of the levels in the data describes those levels, which the data pin down when each level has enough observations. With two contexts it is half the squared difference between their effects. Per draw, for a grouping factor with C levels in the data and effects u_k:

var_stratum

the realised variance of the strata: of the stratum intercepts, or of the \zeta_j of a crossed-dimensions fit (its between_var).

per_context

each context's \sum_k (u_k - \bar u)^2 / (C - 1), and context_var_total their sum.

other_var

the same, summed over any other random effect in the formula; given only when there is one.

vpc_stratum

var_stratum over the sum of the realised variances and the within-stratum variance (see vpc in Realised strata), and vpc_context and vpc_context_total each context's and all contexts' share of that sum. Given for a Gaussian identity fit and for a binomial, Bernoulli or cumulative fit on its latent scale, not for a Poisson or negative-binomial fit. For a crossed-dimensions fit vpc_stratum is decomposition$realised$vpc.

Each is the posterior median, with the equal-tailed credible interval at conf_level in the matching *_ci element; for per_context and vpc_context a matrix with a row per context. Every level counts once, whatever its size. n_strata is the number of strata, and the partition's context_levels the number of each context's levels. print() of a two-model analysis shows the realised VPC and context share of the null model. These values depend little on the prior on the standard deviations when each context and stratum has enough data. With few observations per context each draw also carries the uncertainty about every context's effect, and they then move with the prior. An lme4 fit reports its variance components only.

References

Gelman, A. (2005). Analysis of variance – why it is more important than ever. The Annals of Statistics, 33(1), 1-53. doi:10.1214/009053604000001048

See Also

fit_maihda for the single-model fitter, compare_maihda_groups for the group comparison, and summary.maihda_model for the variance summary.

Examples

data(maihda_health_data)

# One call: null + adjusted fit, VPC summary, and PCV decomposition. Writing the
# dimensions' additive main effects (Gender + Race) gives the fully-specified
# adjusted model; maihda() derives the null by dropping them.
a <- maihda(BMI ~ Age + Gender + Race + (1 | Gender:Race), data = maihda_health_data)
a                              # VPC (null) and PCV (null -> adjusted)
a$pcv                          # proportional change in between-stratum variance
a$formula                      # null:     BMI ~ Age + (1 | stratum)
a$adjusted_formula             # adjusted: null + Gender + Race main effects

# Omitting them is equivalent -- maihda() adds them to the adjusted model and
# emits a message; the null and PCV are identical to the explicit form above.
a0 <- maihda(BMI ~ Age + (1 | Gender:Race), data = maihda_health_data)

plot(a, type = "vpc")          # null model
plot(a, type = "effect_decomp")# adjusted model (additive vs intersectional)

# Crossed-dimensions decomposition: one model, the dimensions' main effects entered
# as RANDOM intercepts. The additive and interaction shares of the between-strata
# variance are read directly from the single fit (no null/adjusted pair).
cc <- maihda(BMI ~ Age + (1 | Gender:Race), data = maihda_health_data,
             decomposition = "crossed-dimensions")
cc                                    # VPC and additive/interaction shares
cc$decomposition$additive_share       # crossed-dimensions analogue of the PCV
cc$formula                            # BMI ~ Age + (1|Gender) + (1|Race) + (1|stratum)

# Add a higher-level grouping variable to also compare across its levels.
# maihda_country_data has a real country grouping (PISA achievement data):
data(maihda_country_data)
a2 <- maihda(math ~ 1 + (1 | gender:ses), data = maihda_country_data,
             group = "country")
a2
plot(a2, type = "group_vpc")
plot(a2, type = "group_pcv")

# Contextual cross-classified MAIHDA: instead of one model per country (group=),
# model the strata CROSSED with country in a single fit. The summary partitions
# the unexplained variance into between-stratum vs. between-country vs. residual,
# and the PCV is computed with country partialled out.
a3 <- maihda(math ~ 1 + (1 | gender:ses), data = maihda_country_data,
             context = "country")
a3
a3$summary$context$vpc_context_total  # the country (general contextual) share
plot(a3, type = "context_vpc")

Build the adjusted-model formula and data for a MAIHDA decomposition

Description

Given a fitted null model's formula (in (1 | stratum) form) and its stored strata metadata, returns the adjusted formula – the null formula plus the additive main effects of the stratum dimensions – and the data carrying any reconstructed binned factors. Returns NULL when fewer than two dimensions are available, because there is no intersection to decompose and the single dimension's main effect would render the stratum random intercept redundant (singular).

Usage

maihda_adjusted_formula(null_formula, strata_vars, autobin_info, data)

Arguments

Value

A list with formula and data, or NULL if fewer than two dimensions are available.

Reconstruct the adjusted-model main-effect terms for a MAIHDA decomposition

Description

For each stratum-defining variable, returns the model term to add as an additive fixed main effect in the adjusted model, plus the data augmented with any reconstructed binned factors. A categorical dimension is used directly; a numeric dimension that make_strata tertile-binned is reconstructed as the same binned factor (using the stored breaks/labels), because make_strata() bins only a temporary copy and leaves the original numeric column intact – adding the raw numeric column would wrongly enter a linear term instead of the factor that defines the strata.

Usage

maihda_adjusted_terms(strata_vars, autobin_info, data)

Arguments

Value

A list with terms (character vector of RHS term names) and data (the input augmented with any .maihda_dim_* binned columns).

Area under the ROC curve (C-statistic), rank-based

Description

Computes the AUC / C-statistic as the Mann-Whitney U statistic: the probability that a randomly chosen case (y == 1) is assigned a higher predicted value than a randomly chosen non-case (y == 0), with ties counting as one half. This needs no external package. An AUC of 0.5 is chance; 1 is perfect separation.

Usage

maihda_auc(prob, y)

Arguments

Value

A single number in [0, 1], or NA_real_ if either class is absent.

References

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.

Examples

maihda_auc(c(0.1, 0.4, 0.35, 0.8), c(0, 0, 1, 1))

Null-restricted parametric-bootstrap fixed effects for an lme4 fit

Description

Internal helper. For each fixed-effect coefficient, refits the model with that coefficient constrained to zero, simulates n_boot responses from the restricted fit, refits the full model on each, and refers the observed Wald statistic to the distribution of |t^*| across those refits. The null is the restricted fit itself, its remaining parameters estimated from the same data, so the reference is a parametric-bootstrap approximation rather than the exact null distribution. Returns the shape maihda_fixed_effects_table produces, with df NA: the reference is an empirical distribution, not a t.

Usage

maihda_bootstrap_fixef(model, n_boot, conf_level, longitudinal = NULL)

Arguments

Details

One coefficient at a time, not one term at a time. A term spanning several design columns – a factor with three or more levels, a polynomial, an interaction between factors – holds several coefficients, each reported on its own row; constraining the whole term would test every one of them to zero while the row claims only its own. The siblings are nuisance parameters and are re-estimated, not deleted: dropping them pushes their effect into the variance components, so the draws come from a model the data do not describe. Which coefficients share a term is a property of how the design was spelled – y ~ f and y ~ fb + fc fit the identical model – and a reference distribution may not turn on that.

For a one-column term the restriction is imposed by dropping the term from the formula and then verifying that the refitted design no longer spans the full model's column space. It usually does not, but R's marginality rules recode a surviving higher-order term to absorb a removed marginal one – . ~ . - x applied to y ~ x * f gives f + x:f, whose full dummy expansion spans exactly the original space – and such a "reduction" constrains nothing. Where the spans agree, and for every coefficient of a multi-column term (no formula drops one contrast column), the constraint is imposed on the fitted design columns directly (maihda_restrict_fixef).

The intercept has no reduced model to simulate from – a MAIHDA intercept is a reference-category level rather than a term that can be dropped – so its p-value and interval are NA, not a Wald z that is miscalibrated in exactly the regime this function exists for.

Value

A data frame with term, estimate, se, statistic, df, p_value, lower and upper, carrying n_boot_ok, n_boot_nonconverged and interval_reliable attributes – the last FALSE when more than half the contributing draws failed to converge, as for the VPC and PCV bootstraps (maihda_report_nonconvergence).

Rows complete on every variable a brms model will use

Description

Mirrors the rows brms retains after its own NA exclusion: complete on the response and its addition-term variables (y | trials(n)), on the fixed-effect terms – evaluated through the model frame, so a transformed predictor whose transformation yields NA/NaN counts as incomplete – and on every variable in the random-effect terms.

Usage

maihda_brms_complete_rows(formula, data)

Arguments

Value

A logical vector over the rows of data.

Cumulative brmsfamily for the brms ordinal engine

Description

Turns whatever the caller supplied for an ordinal engine = "brms" fit into a proper brms::cumulative() family, preserving the cumulative-specific options a user set.

Usage

maihda_brms_cumulative_family(family)

Arguments

Details

A rebuild is needed because family may still be the plain list(family = "cumulative", link = ) marker that the family-string path and maihda_cumulative produce, which brm() cannot use. It must be a rebuild rather than a pass-through for the same reason. Rebuilding from the link ALONE, however, silently reset every other option to brms's defaults: brms::cumulative(threshold = "equidistant") fitted the flexible model, with no warning that a different model had been fitted from the one asked for. threshold and link_disc are therefore carried across when the supplied family actually has them; the marker lists carry neither, so those paths still get brms's defaults.

Note the contrast with the clmm engine, where a structured threshold also has to be handled on the way OUT (maihda_clmm_cutpoints). brms needs no such treatment: its generated quantities expand the constrained thresholds into the full b_Intercept[1..K-1] vector, so brms::fixef() reports all K-1 cut points under every threshold structure and maihda_brms_ordinal_thresholds reads them unchanged.

Value

A brms::cumulative() family object.

Posterior-mean cumulative thresholds of a brms cumulative fit

Description

The brms analogue of clmm's object$model$alpha: the ordered cut points \alpha_k of a brms::cumulative() fit, read as posterior means from brms::fixef(), where they appear as Intercept[1], Intercept[2], ... (any location predictors are separate rows and are dropped here). Returned in threshold order so they pair with the brms::posterior_linpred(re_formula = NA) location – which excludes the thresholds, exactly the latent \eta that maihda_ordinal_category_probs expects (P(Y \le k) = g^{-1}(\alpha_k - \eta)).

Usage

maihda_brms_ordinal_thresholds(model)

Arguments

Details

Unlike clmm's $alpha, this needs no expansion under a structured threshold. brms keeps the constrained parameterisation in the model block only – an equidistant fit samples first_Intercept and delta – and its generated quantities write the full b_Intercept[1..K-1] vector, which is what fixef() reports. Verified on a fitted equidistant model: fixef() carries all K-1 Intercept[k] rows (the raw delta is not among them) and the cut points come back equally spaced. Contrast maihda_clmm_cutpoints, which must do that expansion itself.

Value

A numeric vector of thresholds (length K-1 for K categories).

Inject sampling weights into a brms formula

Description

Rewrites y ~ ... as y | weights(w) ~ .... An existing addition term (e.g. an aggregated-binomial y | trials(n)) is extended with + weights(w); a formula that already carries a weights() addition term is rejected (the two weight specifications would conflict).

Usage

maihda_brms_weights_formula(formula, wcol)

Arguments

Value

The rewritten formula (same environment).

Crossed-dimensions variance summary (brms)

Description

brms counterpart of maihda_cc_summary_lme4: computes the additive / interaction partition per posterior draw and returns posterior point estimates with credible intervals for the VPC and the shares (no bootstrap – the posterior already supplies the interval). A contextual random intercept (object$context_info set) enters the per-draw VPC denominator and is reported as its own component row and context element. Each variance component is the posterior median of its draws, and the additive, between-strata and context totals are sums of those components. When the draws hold the level effects, decomposition$realised carries the realised strata and context$realised the realised partition (see the sections Realised strata and Realised strata and contexts of maihda).

Usage

maihda_cc_summary_brms(object, cc, conf_level, point = c("median", "mean"))

Arguments

Value

A list with variance_components, vpc_result, decomposition and context (NULL without a context).

Crossed-dimensions variance summary (lme4)

Description

Internal helper for summary.maihda_model when the model is a crossed-dimensions MAIHDA fit (object$cc_info set). Partitions the crossed random-effect variances into the additive (sum of the dimension REs) and interaction (intersection RE) components, builds the variance-components table and the VPC, and – when bootstrap = TRUE – adds parametric-bootstrap intervals for the VPC and the additive/interaction shares. When the fit also carries a contextual random intercept (object$context_info set), the context variance enters the VPC denominator and is reported as its own component row and context element.

Usage

maihda_cc_summary_lme4(object, cc, vc, bootstrap, n_boot, conf_level)

Arguments

Value

A list with variance_components, vpc_result, decomposition and context (NULL without a context).

Guard against longitudinal ids reused across strata

Description

(time | id) treats every row sharing an id value as the SAME person, so ids that are only unique within a site or group (person "1" in every site) silently merge different people's trajectories into one level-2 unit. An id appearing in more than one stratum is the observable symptom of that (a person's intersectional stratum is a person-level attribute, so a genuinely unique id maps to exactly one stratum); reject it with guidance rather than guessing a nesting. Called after strata resolution, when the stratum column exists. Rows missing id or stratum are ignored (the engines drop them).

Usage

maihda_check_longitudinal_ids(data, id)

Arguments

Value

NULL, invisibly; stops on ambiguous ids.

Guard against unidentified growth terms (distinct-times checks)

Description

A degree-d growth curve is only estimable when the observed times actually span it. Two conditions, both required: the polynomial basis needs at least d + 1 distinct observed time values (on a single time even the linear slope does not exist; on two times a quadratic is collinear with the linear term), and the person-level slope block (time | id) needs at least one id measured at two distinct times (ids repeating only at a single time are replicate measurements, not a growth design). Data failing either check can still FIT: lme4 drops rank-deficient fixed-effect columns but keeps the random-effect columns, leaving the likelihood flat (or near-flat) in the affected covariance parameters, so the "slope variances" it returns are arbitrary – observed: optimizer starting values reported back verbatim on all-equal times, and a ~1300-fold optimizer-dependent swing in the between-stratum variance at unobserved times for a quadratic on two waves – with at most a "boundary (singular) fit" note. Reject with guidance instead.

Usage

maihda_check_longitudinal_times(data, id, time, time_degree)

Arguments

Details

Called on the full input by maihda_validate_longitudinal() and again on the analytic sample by fit_maihda(): dropped rows (subset, missing outcome/covariate/weight) can only REDUCE the distinct times, so the early check never rejects a fittable model, and the re-check catches data whose identifying occasions sit on rows the fit drops.

Value

NULL, invisibly; stops when the growth terms are unidentified.

Expanded cut points of a cumulative (clmm) fit

Description

The K-1 cut points of a K-category cumulative model, whatever threshold structure was fitted.

Usage

maihda_clmm_cutpoints(model)

Arguments

Details

clmm's $alpha holds the free threshold parameters, which equal the cut points only under the default threshold = "flexible". A structured threshold request ("equidistant", "symmetric", "symmetric2", all part of clmm's documented API and reachable through fit_maihda's ...) leaves $alpha holding a shorter, differently-parameterised vector – threshold.1 and spacing for the equidistant case – so reading it as the cut points understates the number of categories and mislocates every one of them. The expanded cut points are stored on the fit as $Theta (a 1 \times (K-1) matrix), equivalently $tJac %*% $alpha; $tJac is the identity under "flexible", so this returns $alpha unchanged there.

Value

A numeric vector of K-1 increasing cut points, named "1|2", "2|3", ... when the fit supplies those labels.

Location linear predictor of a cumulative (clmm) fit

Description

predict.clmm does not exist, so the location part \eta = x'\beta (+ u) is built directly: the fixed design matrix is constructed with the fit's factor coding – the levels it fitted and the contrast matrix it applied to each factor, whatever options(contrasts = ) says now – AND its transformation basis (so a data-dependent term such as scale(x) uses the fit's centre and scale rather than recomputing them from newdata) and multiplied by the location coefficients beta (a clmm has no intercept column – it is absorbed by the thresholds – so beta's names select the right columns), any formula offset term is evaluated on newdata and added (so an offset-only null model still predicts its offset), and include_re adds each row's stratum conditional mode (an unseen stratum contributes 0 – the zero-effect fallback that predict_maihda only reaches when allow_new_levels = TRUE, having otherwise rejected unseen strata upstream). Everything is on the latent (link) scale; map through maihda_ordinal_eta_to_score for the response-scale expected category score.

Usage

maihda_clmm_linpred(object, newdata = NULL, include_re = TRUE)

Arguments

Value

A numeric vector of latent-scale location predictions.

Stratum random-effect table for a cumulative (clmm) fit

Description

Mirrors maihda_stratum_ranef_lme4(): one row per stratum with the conditional mode, its conditional standard error (from ordinal::condVar(), which returns conditional variances), and a 95% interval. At a boundary fit (zero between-stratum variance) the conditional distribution collapses on 0, so the SE is 0.

Usage

maihda_clmm_stratum_ranef(object)

Arguments

Value

A data frame with stratum, stratum_id, random_effect, se, lower_95, upper_95.

Threshold (cut-point) estimates of a cumulative (clmm) MAIHDA fit

Description

The thresholds \alpha_k take the place of the intercept in a cumulative model: P(Y \le k) = g^{-1}(\alpha_k - \eta). Reported as the K-1 expanded cut points (maihda_clmm_cutpoints), so the table means the same thing under every threshold structure rather than echoing the free parameters of a constrained fit. Standard errors come from the Hessian-based vcov() (hence Hess = TRUE at fit time) and degrade to NA when unavailable; under a structured threshold they are the delta-method SEs of the cut points, J V J' for the fit's threshold Jacobian $tJac. That Jacobian is the identity under the default threshold = "flexible", where the transform reproduces the direct vcov() SEs exactly.

Usage

maihda_clmm_thresholds(object)

Arguments

Value

A data frame with term, estimate, se.

Variance components of a cumulative (clmm) MAIHDA fit

Description

Reads the between-stratum variance from ordinal::VarCorr() and pairs it with the latent-scale level-1 variance (\pi^2/3 for logit, 1 for probit), matching the latent treatment of binomial models in the other engines.

Usage

maihda_clmm_variances(object)

Arguments

Value

A list with stratum and residual variances.

Join non-empty caption lines for a plot

Description

Pastes its arguments into a single newline-separated caption, dropping any that are NULL, NA, or empty. Returns NULL when nothing remains, so a plot with no caption content keeps a clean (caption-free) look.

Usage

maihda_compose_caption(...)

Arguments

Value

A single newline-joined string, or NULL.

Containment degrees of freedom for a Gaussian mixed-model fixed effect

Description

Internal helper. Returns between-within (containment) denominator degrees of freedom for each fixed-effect column of a Gaussian lmerMod: a column absorbed by a random-effect term is tested against that term's units minus the columns it absorbs, a column absorbed by none against n minus the random-effect levels. A column counts as absorbed when, within every level of the grouping factor, it is a linear combination of that term's random-effect covariates – constant within the group for a random intercept, any group-specific linear function of time for a random slope on time.

Usage

maihda_containment_df(model)

Arguments

Details

Where several terms absorb a column the smallest applies, taken over the terms that materially contribute to that coefficient's variance (maihda_variance_shares); a term whose variance component is estimated at essentially zero cannot limit precision. That filter can narrow the set but never empty it, so a singular fit keeps its stratum degrees of freedom.

The per-source degrees of freedom are deliberately not Satterthwaite-combined: in the balanced one-way case the stratum and residual shares of a stratum-level contrast come from the same mean square, so combining them as independent chi-squares inflates the result several-fold. Cross-checked against pbkrtest in tests/testthat/test-fixed-effect-df.R; pbkrtest is a suggested package and those checks skip without it.

Value

A named numeric vector of degrees of freedom, or NULL.

Contextual cross-classified variance summary (brms)

Description

brms counterpart of maihda_context_summary_lme4: computes the stratum / context / residual partition per posterior draw and returns posterior point estimates with credible intervals for the between-stratum VPC and the contexts' total share (no bootstrap – the posterior supplies the interval). Each variance component is the posterior median of its draws, and the contexts' total variance is the sum of the context components. When the draws hold the level effects, context$realised carries the realised partition (see the section Realised strata and contexts of maihda).

Usage

maihda_context_summary_brms(
  object,
  ctx,
  conf_level,
  point = c("median", "mean")
)

Arguments

Value

A list with variance_components, vpc_result, context.

Contextual cross-classified variance summary (lme4)

Description

Internal helper for summary.maihda_model when the model carries a contextual random intercept (object$context_info set, fit_maihda(context = )) without the crossed-dimensions decomposition. Partitions the unexplained variance into between-stratum vs. between-context (one share per context variable) vs. residual. The headline VPC stays the between-stratum share of all unexplained variance – numerically identical to the generic single-stratum summary, which folds the context into "Other random effects" – but the context is now named, given its own component row(s), and returned as a context element.

Usage

maihda_context_summary_lme4(object, ctx, vc, bootstrap, n_boot, conf_level)

Arguments

Value

A list with variance_components, vpc_result, context.

Cross-National Educational Achievement Data for MAIHDA

Description

A cross-national dataset for demonstrating how Multilevel Analysis of Individual Heterogeneity and Discriminatory Accuracy (MAIHDA) can be used to compare intersectional inequality across a higher-level grouping variable (here, country) with compare_maihda_groups and maihda. Each row is a 15-year-old student; the intersectional strata are formed by gender and socioeconomic status (ses), and the outcome is the PISA mathematics score.

Usage

maihda_country_data

Format

A data frame with 3,600 rows (600 students in each of 6 countries) and 7 variables:

country

Factor; one of Finland, Germany, United Kingdom, Italy, Japan, Mexico. The higher-level grouping variable.

gender

Factor; student gender (female/male). A stratum dimension.

ses

Factor; socioeconomic status as global tertiles (Low/Medium/High) of escs, computed on the pooled sample so a band means the same in every country. A stratum dimension.

escs

Numeric; the PISA index of economic, social and cultural status (the continuous measure underlying ses).

math

Numeric; PISA mathematics score (first plausible value). The primary outcome.

reading

Numeric; PISA reading score (first plausible value).

low_math

Factor; "Yes" if math is below 420 (PISA proficiency Level 2 baseline), else "No". A binary outcome for logistic examples.

Details

Intersectional inequality (the between-stratum share of variance, VPC/ICC) in mathematics achievement differs across the six countries, which is what makes the dataset a useful showcase for the group-comparison workflow.

The intersectional strata are gender:ses (2 x 3 = 6 strata). A canonical MAIHDA "null" model is math ~ 1 + (1 | gender:ses); comparing its VPC across countries quantifies how much joint gender-by-class inequality in achievement varies between countries.

Note

This is a teaching/illustration dataset only. It uses a single PISA plausible value for each score and does not carry the PISA survey weights or complex sampling design, so results are not survey-representative and should not be used for substantive cross-national inference. (For your own survey data, the package supports design-weighted MAIHDA via the sampling_weights argument of fit_maihda / maihda.)

Source

Derived from the OECD Programme for International Student Assessment (PISA) 2018 student questionnaire data (OECD (2019), PISA 2018 Database), accessed and cleaned via the learningtower R package (MIT licensed), https://CRAN.R-project.org/package=learningtower. A balanced random subsample of 600 complete-case students per country was taken (seed 2026). The data preparation script is in data-raw/maihda_country_data.R.

Examples

data(maihda_country_data)

# Compare intersectional (gender x SES) inequality across countries
analysis <- maihda(
  math ~ 1 + (1 | gender:ses),
  data = maihda_country_data,
  group = "country"
)
analysis
plot(analysis, type = "group_vpc")

Build the crossed-dimensions-model formula and data for a MAIHDA decomposition

Description

The crossed-dimensions alternative to the two-model (fixed-effects PCV) decomposition (the function name keeps the historical "cross_classified" spelling). Given a null model's formula (in (1 | stratum) form, carrying only the covariates) and the stratum metadata, returns the single crossed formula – the covariates plus an additive random intercept for each stratum dimension plus the intersection (stratum) random intercept – together with the data carrying any reconstructed binned factors. In the fitted model each dimension's RE variance is that dimension's additive main-effect variance and the stratum RE variance is the interaction beyond additive; see maihda.

Usage

maihda_cross_classified_formula(
  null_formula,
  strata_vars,
  autobin_info,
  data,
  interaction_group = "stratum",
  context = NULL
)

Arguments

Details

Returns NULL when fewer than two dimensions are available (there is no intersection to decompose). The dimension grouping factor reuses the dimension's own column for a categorical dimension and the reconstructed .maihda_dim_* tertile factor for an auto-binned numeric dimension (via maihda_adjusted_terms), so the additive REs are crossed on exactly the levels that define the strata.

Value

A list with formula, data, dim_groups (a named character vector mapping each strata_var to its random-effect grouping-factor name) and interaction_group ("stratum"); or NULL if fewer than two dimensions are available.

Cumulative (ordinal) family marker for MAIHDA models

Description

Specifies a cumulative (proportional-odds) model for an ordinal outcome in fit_maihda / maihda, with a choice of link: maihda_cumulative("logit") (the default, equivalent to family = "ordinal") or maihda_cumulative("probit"). It plays the role a stats family object plays for the other families – there is no cumulative family constructor in stats, and using brms::cumulative() would require brms for a frequentist fit.

Usage

maihda_cumulative(link = c("logit", "probit"))

Arguments

Value

A family marker list with elements family = "cumulative" and link.

See Also

fit_maihda

Examples

maihda_cumulative()
maihda_cumulative("probit")

Describe the MAIHDA sample before fitting

Description

Builds the standard “Table 1” descriptives of a MAIHDA analytic sample before any model is fitted: the total and complete-case sample sizes, the observed vs. expected intersectional strata (and which are empty or small), the distribution of each stratum-defining dimension, family-aware per-stratum outcome summaries, missing-data accounting, and – for a contextual cross-classified design – the context units. It is the pre-model counterpart of maihda_table.

The description comes from the same machinery as the model: the strata are built by make_strata via the same formula-shorthand resolution as fit_maihda, the analytic sample reproduces the engines' own row selection (missing outcome, covariates, stratum dimensions, context, and – for a weighted design – missing/non-positive sampling weights), and the outcome family uses the same auto-detection as the fitters. Stratum IDs, labels, and counts therefore match a subsequent make_strata() / fit_maihda() / maihda on the same formula and data exactly.

Usage

maihda_describe(
  x,
  data = NULL,
  context = NULL,
  family = NULL,
  sampling_weights = NULL,
  flag_stratum_n = 20,
  include_empty_strata = TRUE,
  autobin = TRUE,
  digits = 3,
  weights = NULL
)

Arguments

Details

Family-aware outcome summaries. The per-stratum (and overall) outcome summary depends on the resolved family, and the applied summary type is recorded in $outcome_summary – a 0/1 outcome is never described with a Gaussian mean/SD as if continuous:

• gaussian (and other continuous families): mean, SD, median, min, max;

• binomial (including an aggregated outcome written either as cbind(success, failure) or, for brms, as success | trials(n)): event count, trials, and the observed proportion. Trials are the binomial denominator, so a row with a missing or non-positive trial count has no observed outcome;

• poisson / negative binomial: the same numeric summaries, read as the observed count mean (rate);

• cumulative (ordinal): the mean and (lower) median category score (categories scored 1..K in level order, the same convention as the package's plots), with the full category distribution in $outcome_levels.

Intersectional sparsity. $overview contrasts n_strata_observed with n_strata_expected (the full Cartesian product of the observed dimension levels) and n_empty_strata; with include_empty_strata = TRUE the empty combinations are enumerated as rows of $strata (with stratum = NA, so real stratum IDs stay identical to a subsequent fit). With many small strata the stratum effects are shrunk toward the mean and the variance components are imprecise with either engine.

Missing data. $overview separates the rows lost to a missing stratum dimension (make_strata() routes them to the NA stratum) from rows with a missing outcome, and reports the resulting complete-case n_analytic – computed with the same row mask the engines apply, so it equals nrow(fit_maihda(...)$data).

Weights. With sampling_weights, weighted counts and weighted outcome means/proportions are reported alongside the unweighted ones, and rows with a missing or non-positive weight are excluded from the analytic sample (as the weighted engines do). As elsewhere in the package, these are population-weighted point summaries – not a complex survey design; no design-based variances are computed.

Units. maihda_describe() treats rows as the unit of description. For longitudinal (long-format) data a row is a measurement occasion, not a person; describing a fitted longitudinal model (fit_maihda(id =, time =)) additionally reports n_individuals and the median occasions per individual in $overview.

Fitted-model input. A maihda_model (or the maihda_analysis bundle from maihda) can be described post hoc; the formula, data, family, context, and sampling weights are taken from the fit, so the description covers the exact analytic sample the model used. For the wemix/ordinal engines the stored data already is the pre-filtered analytic sample, so total and analytic counts coincide there.

Value

An object of class maihda_describe: a list of export-ready data frames (pass to write.csv() or knitr::kable()) plus metadata:

See Also

maihda_table for the post-model tables, make_strata, fit_maihda, plot for the descriptive plots.

Examples

data(maihda_health_data)
desc <- maihda_describe(BMI ~ Age + (1 | Gender:Race:Education),
                        data = maihda_health_data, flag_stratum_n = 20)
desc
desc$overview
head(desc$strata)

# Binary outcome: auto-detected binomial, proportions instead of means
desc_bin <- maihda_describe(Obese ~ (1 | Gender:Race:Education),
                            data = maihda_health_data)
desc_bin$outcome_levels

# Contextual cross-classified design (country as context)
data(maihda_country_data)
desc_ctx <- maihda_describe(math ~ escs + (1 | gender:ses),
                            data = maihda_country_data, context = "country")
desc_ctx$context

# Export-ready tables
# write.csv(desc$strata, "table1_strata.csv", row.names = FALSE)

Discriminatory accuracy of a binary MAIHDA model

Description

Bundles the individual-level discriminatory-accuracy summaries for a binomial MAIHDA model: the AUC / C-statistic (how well the model's predicted probabilities separate cases from non-cases) and the Median Odds Ratio. Applied to a strata-only (null) model, the AUC is the discriminatory accuracy of the intersectional strata themselves – Merlo's central quantity; comparing it with an adjusted model shows whether individual covariates beyond stratum membership sharpen classification. The AUC is computed for any binomial link; the Median Odds Ratio is reported only for the logit link and is NA otherwise (e.g. for a probit fit), since the MOR is an odds-ratio-scale quantity.

The AUC is apparent (in-sample). It scores the same observations used to estimate the model – the fixed effects, the variance components, and the shrunken stratum BLUPs – so it is an apparent (resubstitution) AUC and is optimistically biased, the more so with small or sparse strata. It is reported as the conventional descriptive MAIHDA discriminatory accuracy (Merlo 2018), not as a cross-validated estimate of out-of-sample predictive discrimination; the returned object carries apparent = TRUE and print() says so. For genuine predictive accuracy, validate with an out-of-fold (group-aware) scheme or an optimism correction.

Aggregated-binomial fits are supported in all three forms R accepts – an lme4 cbind(success, failure) response, an lme4 response in [0, 1] whose trial counts are supplied as weights = (see glm), and a brms y | trials(n) response: the AUC is the count-weighted C-statistic over the implied individual-level 0/1 data, and n_case / n_control are the total successes / failures. All three spellings of the same data give the same AUC. This includes individual 0/1 records collapsed to frequency cells (every row all-success or all-failure), whose weights are trial counts like any others; see binomial_weights to override that reading.

Scope. When the model carries random effects beyond the intersectional partition – a contextual (1 | school) from fit_maihda(context = ) or an explicit extra grouping such as (1 | site) – the headline auc is the intersectional-scope concordance: it excludes those other random effects but keeps the fixed effects plus the stratum random effect (and, for a crossed-dimensions fit, the additive dimension effects). The concordance of the full model including the other random effects is reported separately as auc_full (auc_scope = "intersectional"). Caveat – this is not strata-only discrimination. The intersectional-scope score retains the entire fixed-effects predictor, so when the model is adjusted for individual-level covariates (e.g. age, income that vary within strata) those covariates enter this AUC too; it is the concordance of the adjusted fixed effects plus the intersectional random effect(s), not of the strata alone, and it then matches the between-stratum MOR's scope only when the fixed part is intercept-only. For a strata-only discriminatory accuracy, score the null (strata-only) model. For the canonical single-(1 | stratum) model with no other random effects, the full and intersectional scopes coincide and only auc is reported (auc_scope = "model", auc_full absent).

Usage

maihda_discriminatory_accuracy(
  model,
  binomial_weights = c("auto", "trials", "analytic")
)

Arguments

Value

An object of class maihda_da: a list with auc, auc_scope, auc_full, mor, n_case, n_control, family, link, engine, weighted, weight_type, precision_weights_ignored and apparent (always TRUE – the AUC is in-sample; see Description). mor is NA for a non-logit binomial link, where the AUC is still reported. For an aggregated-binomial fit n_case / n_control are the total successes / failures. weighted is TRUE only when the AUC is a genuinely weighted (population-mass) concordance, with weight_type "sampling" for a design-weighted fit (each observation contributes its sampling weight as case/control mass, estimating the population discriminatory accuracy); n_case / n_control stay unweighted observation counts. weight_type is NULL for an unweighted AUC. An aggregated-binomial fit is reported weighted = FALSE with weight_type = NULL: its count-weighted AUC equals the ordinary individual-level concordance over the implied 0/1 data (the trial counts are real observations, not sampling weights), so it is not a design-weighted population quantity. Aggregation is recognised from an lme4 cbind(success, failure) matrix response, from non-unit integral lme4 weights= on a response in [0, 1] (the trial counts, see glm), or from a brms y | trials(n) term. The weights-based form covers a response of exactly 0/1 – individual records collapsed to frequency cells – because for the binomial family a prior weight is a trial count: the weighted log-likelihood of one 0/1 row is exactly that of w trials sharing its outcome, so the fit equals the row-expanded one and its AUC should too. Pass binomial_weights to override. When a response times its trial counts is not a whole number of successes – a malformed binomial, which glmer itself warns about as "non-integer #successes" – the fractional case/control mass is carried into the AUC as it stands, with a warning, rather than rounded into whole observations that were never collected; n_case / n_control are then fractional. Weights not read as trial counts – non-integral weights=, or any weights under binomial_weights = "analytic" – are ignored by the AUC, which is then the ordinary observation-level concordance (weighted = FALSE), with precision_weights_ignored = TRUE flagging that such weights were present.

References

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

See Also

maihda_auc, maihda_mor

Examples

# Obese (Yes/No) by intersectional strata of Gender x Race
strata <- make_strata(maihda_health_data, vars = c("Gender", "Race"))
d <- maihda_health_data
d$stratum <- strata$data$stratum
m <- fit_maihda(Obese ~ (1 | stratum), data = d, family = "binomial")
maihda_discriminatory_accuracy(m)

Drop fixed terms from a formula without dropping the grand mean

Description

The single place a MAIHDA derives a reduced model by removing the stratum dimensions' fixed main effects: the two-model null, the crossed-dimensions base formula, the longitudinal null and the per-group formulas all come through here.

Usage

maihda_drop_fixed_terms(formula, terms, data, notify = FALSE, fn = "maihda")

Arguments

Details

update(f, . ~ . - a - b) removes the named terms, but on a no-intercept formula it removes the grand mean with them. y ~ x + a + b + (1 | a:b) and y ~ 0 + x + a + b + (1 | a:b) are the SAME adjusted model – with no intercept R codes a as cell means, so the two designs span one space and fit identical values – yet the plain reduction turns the first into y ~ x + (1 | stratum) and the second into y ~ x + (1 | stratum) - 1. The second null has no way to represent the outcome's mean, so the stratum random intercept absorbs it: the between-stratum variance, and with it the VPC, MOR and PCV, become functions of the arbitrary origin of the response rather than of the strata.

The grand mean is therefore restored whenever the reduction loses it AND the source model carried it. Conditioning on the source keeps the reduced model nested: a fit that genuinely has no intercept in its span (y ~ 0 + x + a + b with numeric dimensions – a regression through the origin, not an equivalent recoding) keeps its parameterization, because adding a grand mean there would put a column in the null that the adjusted model does not have. A reduction that still spans the intercept some other way – a surviving factor covariate coded as cell means, or a fixed part that collapses to 1 once the bars are stripped – is left exactly as it was.

Value

The reduced formula, in the same environment as formula.

Fit a cumulative MAIHDA model via ordinal::clmm

Description

Internal engine call for fit_maihda(engine = "ordinal"). Builds the analytic sample (complete cases on the model variables) so the stored data matches the rows clmm fits, then calls ordinal::clmm() with Hess = TRUE (needed for the threshold standard errors). The analytic frame is passed by NAME (bound in a private environment) so the call clmm stores stays one line: ordinal's print/summary methods deparse call$data, and embedding the frame there made printing an ordinal fit dump the whole data set. Any formula offset() term is moved ahead of the (1 | stratum) term first, because clmm reads the offset by position in the bar-free variables list (see maihda_offset_before_bars).

Usage

maihda_fit_clmm(formula, data, family, dot_vals)

Arguments

Details

A Laplace fit that ends with every stratum effect at exactly zero, or with a between-stratum SD below 1% of the latent residual SD, is compared with the boundary fit: clmm()'s general routine (clmm.control(useMatrix = TRUE)) started from the first fit's thresholds and coefficients and an SD of zero. The boundary fit is the one returned unless its log-likelihood is below the first fit's by more than 1e-9 of its size.

The model is fitted on its fixed-effects design with every column centred and scaled, the boundary fit included, and returned as the clmm() fit of the formula at that maximum. A covariate far from zero or on a large or tiny scale can leave clmm()'s optimiser at its starting values or short of the maximum, and clmm() takes the gradient and Hessian by differences with a fixed step, too coarse for the coefficient of such a covariate. The returned fit carries the gradient and Hessian of the rescaled fit, mapped exactly, so its standard errors are those of that fit. Under threshold = "symmetric2" the thresholds cannot take up a shift, and the columns are scaled only. A fit given a start is made as the formula stands, since the start is in the formula's own terms.

Value

A list with model (the clmm fit), data (the analytic data frame actually fitted) and coding (the fit's factor levels and contrast matrices, taken while the options that shaped the design are in force; NULL if they could not be read).

Fit a design-weighted MAIHDA model via WeMix

Description

Internal engine call for fit_maihda(engine = "wemix"). Builds the analytic sample (complete cases on the model variables and the weight column, positive weights only) so the stored data matches the rows WeMix fits, attaches the constant level-2 weight column (strata are exhaustive population cells, sampled with certainty), and calls WeMix::mix() with the unconditional weights c(level1, level2).

Usage

maihda_fit_wemix(formula, data, family, sampling_weights, dot_vals)

Arguments

Value

A list with model (the WeMixResults), data (the analytic data frame actually fitted, including the weight columns) and coding (the fit's factor levels and contrast matrices, taken while the options that shaped the design are in force; NULL if they could not be read).

Fixed-effect table with Wald statistics

Description

Internal helper that assembles the fixed_effects slot of a maihda_summary for the frequentist engines: the Wald statistic (estimate / se), its two-sided p-value and the matching Wald interval at conf_level. With df the reference is a t on those degrees of freedom, otherwise the normal.

Usage

maihda_fixed_effects_table(term, estimate, se, conf_level = 0.95, df = NULL)

Arguments

Value

A data frame with term, estimate, se, statistic, df, p_value, lower and upper. df is NA on the normal path.

Does a formula's fixed part carry the grand mean?

Description

TRUE when the fixed-effect design the formula implies – built on data the way the engines build it, model.matrix() on the bar-free formula – carries the all-ones (intercept) vector in its COLUMN SPACE, whether or not the formula writes 1. Writing 0 + does not settle it: with no intercept R codes the FIRST factor term as cell means, so y ~ 0 + f + x still spans the grand mean, while y ~ 0 + x with a numeric x does not. Membership is tested by QR rank – the all-ones vector lies in colspace(X) exactly when rank([X, 1]) == rank(X) – mirroring maihda_re_lhs_spans_intercept() for the random part – but only for the 0 + / - 1 spellings: a formula that writes an intercept returns TRUE without building a design, because model.matrix() always emits the all-ones column for one. Returns NA – never a guess – when the design cannot be built or ranked (an absent or all-NA column, a parse error), so neither caller rewrites a formula or warns about one on the strength of the written 0 + alone.

Usage

maihda_fixed_spans_intercept(formula, data)

Arguments

Value

TRUE if the fixed design spans the intercept, FALSE if it provably does not, NA if the design could not be built.

Glance at a MAIHDA model or analysis

Description

glance methods that return the MAIHDA headline as a one-row tibble: the variance partition coefficient (VPC/ICC), and – for a maihda_analysis – the proportional change in variance (PCV), plus the additive/interaction shares and the discriminatory accuracy (AUC, MOR) for a binomial outcome. The layout is uniform across the lme4, brms, WeMix and ordinal engines. No other package emits this row from the underlying fit, because PCV needs the null+adjusted pair that only a maihda_analysis holds.

Usage

## S3 method for class 'maihda_summary'
glance(x, ...)

## S3 method for class 'maihda_model'
glance(x, ...)

## S3 method for class 'maihda_analysis'
glance(x, ...)

Arguments

Value

A one-row tibble. glance.maihda_analysis adds pcv (with pcv.conf.low/pcv.conf.high when bootstrapped; a brms PCV is a point estimate), the adjusted-model auc.adjusted, nobs, family and mode to the columns produced for a single summary.

See Also

maihda_tidiers for the per-estimate tidy() methods.

Examples

data("maihda_health_data")
a <- maihda(BMI ~ Age + Gender + Race + Education + (1 | Gender:Race:Education),
            data = maihda_health_data)
glance(a)

NHANES Health Data Subset for MAIHDA Use

Description

A pedagogical subset of the National Health and Nutrition Examination Survey (NHANES) dataset, serving as a real-world example for Multilevel Analysis of Individual Heterogeneity and Discriminatory Accuracy (MAIHDA). Contains selected records demonstrating intersectional demographic health inequalities.

Usage

maihda_health_data

Format

A data frame with 3,000 rows and 7 variables:

BMI

Body Mass Index (kg/m^2), a continuous outcome variable.

Obese

Factor indicating obesity status (No/Yes).

Age

Age in years at screening, a continuous covariate.

Gender

Gender of the participant (male/female).

Race

Self-reported race/ethnicity.

Education

Educational attainment level.

Poverty

Poverty to income ratio, a continuous covariate. Some values may be missing.

Note

This is a teaching/illustration dataset only. It is a non-random subsample and does not carry the NHANES survey weights or complex sampling design, so results are not survey-representative and should not be used for substantive population inference. (For your own survey data, the package supports design-weighted MAIHDA via the sampling_weights argument of fit_maihda / maihda.)

Source

Derived from the NHANES R package. Original data collected by the Centers for Disease Control and Prevention (CDC).

Examples

data(maihda_health_data)

# Example usage:
# strata_result <- make_strata(maihda_health_data, vars = c("Gender", "Race", "Education"))
# model <- fit_maihda(BMI ~ Age + (1 | stratum), data = strata_result$data)

Information criteria for MAIHDA models

Description

Reports the relative-fit information criteria for one or more MAIHDA models, to help choose between model structures (different covariate sets, strata definitions, or families) – a question the VPC/ICC and PCV do not address. The criteria reported depend on the engine: AIC and BIC for the likelihood engines (lme4, and ordinal::clmm), and the Bayesian WAIC and LOOIC (leave-one-out information criterion) for brms. Lower is better for all four.

Usage

maihda_ic(..., model_names = NULL)

Arguments

Details

REML vs ML. lmer fits Gaussian models by REML by default, and a REML log-likelihood (hence its AIC/BIC) is not comparable across models with different fixed effects – exactly the canonical MAIHDA null-vs-adjusted comparison. When more than one model is supplied, maihda_ic() therefore refits any REML lmer model with maximum likelihood (refitML) before computing AIC/BIC, matching the behaviour of anova() on lme4 models; the estimator column records when this happened. For a single model the criterion is reported as fitted (the estimator column then reads "REML").

Comparability. Like the VPC, information criteria are only comparable across models fitted to the same analytic sample (same rows and outcome) with the same weights – prior (precision) weights and sampling (design) weights each change which likelihood, or pseudo-likelihood, is being maximised, so the criteria of a weighted and an unweighted fit of the identical model are not on a common scale. Unlike the VPC and PCV, they do not always need the same family: a criterion compares log-likelihoods of the same observations, so Poisson and negative-binomial fits of the same counts are comparable, as are binomial or cumulative (ordinal) fits that differ in their link, and a delta is reported for them. Any other change of family or link withholds it – a Gaussian vs a Poisson fit compares a density with a probability, and the continuous families are compared only within the same family and link. Criteria are never comparable between the likelihood engines and brms (AIC/BIC vs WAIC/LOOIC are different scales). When the supplied models differ in any of these respects maihda_ic() warns and omits the delta column, still reporting each model's own criteria.

Parameter count. df counts the estimated parameters. A negative-binomial lme4 fit with a fixed theta (family = MASS::negative.binomial(theta)) therefore has one fewer than lme4's own count, which includes theta for every negative-binomial family, so its AIC is 2 and its BIC log(n) below AIC() and BIC() on the fitted model; glm() counts that family the same way. A family = "negbinomial" fit estimates theta and counts it.

lme4 GLMM likelihoods. The criteria of a quadrature fit (nAGQ > 1) are the complete likelihood and compare with a Laplace fit's: where lme4's logLik() of such a fit leaves out the saturated log-likelihood, maihda_ic() adds it back. A glmer fit of a family with a scale parameter – a Gaussian with a non-identity link, Gamma, inverse Gaussian – reports NA criteria (estimator "ML (glmer scale family: no likelihood)"): lme4's logLik() for those families is not the marginal likelihood, and on test fits it exceeded the largest value that likelihood attains.

Predictive target of the Bayesian criteria. brms::waic() and brms::loo() are computed from pointwise log-likelihoods conditional on the fitted random effects, so WAIC/LOOIC assess prediction of new observations within the strata (and persons or contexts) already represented in the data – not performance for a new, unseen intersectional stratum. Choosing between strata definitions on LOOIC therefore compares conditional predictive fit; generalisation to new strata is a leave-one-group-out cross-validation question (e.g. brms::kfold() with group = "stratum"), which this package does not wrap. AIC/BIC for the likelihood engines are instead computed from the marginal likelihood (random effects integrated out) – a further reason the likelihood and Bayesian criteria are never comparable with each other.

Design-weighted fits. For the wemix (design-weighted) engine the criteria are reported as NA: a pseudo-likelihood with sampling weights does not define a standard AIC/BIC. A brms fit with sampling_weights is treated the same way: its sampling weights enter as likelihood weights, giving a pseudo-posterior whose weighted pointwise log-likelihoods are not log predictive densities, so WAIC/LOOIC are likewise reported as NA (the estimator column reads "Bayesian (weighted pseudo-posterior)").

Value

A data.frame of class maihda_ic with one row per model and the columns that apply: model, n (analytic sample size), estimator, df (number of parameters; likelihood engines), logLik, AIC, BIC (likelihood engines), WAIC, LOOIC (brms), and – when more than one model is supplied – delta (the difference from the best model on the primary criterion: AIC for the likelihood engines, LOOIC for brms). Columns that are entirely NA across the supplied models are dropped.

See Also

compare_maihda, which reports these criteria alongside the VPC/ICC, and calculate_pcv for the variance decomposition.

Examples

strata <- make_strata(maihda_sim_data, vars = c("gender", "race"))
null_model <- fit_maihda(health_outcome ~ 1 + (1 | stratum), data = strata$data)
adj_model  <- fit_maihda(health_outcome ~ age + (1 | stratum), data = strata$data)

# AIC/BIC for two nested structures (REML lmer fits are ML-refitted first)
maihda_ic(null_model, adj_model, model_names = c("Null", "Adjusted"))

# Or straight from a one-call maihda() analysis (null + adjusted rows)
a <- maihda(health_outcome ~ age + gender + race + (1 | gender:race),
            data = maihda_sim_data)
maihda_ic(a)

Information criteria for a single MAIHDA model

Description

Internal worker for maihda_ic: returns a one-row data frame of the fit criteria for one maihda_model, dispatched on the fitted object's class (mirroring maihda_fit_diagnostics).

Usage

maihda_ic_one(model, ml = FALSE)

Arguments

Value

A one-row data frame with n, estimator, df, logLik, AIC, BIC, WAIC, LOOIC (NA where not applicable to the engine).

Flag strata with credibly non-zero intersectional interaction

Description

Reports, for each intersectional stratum, the interaction component of its outcome – the stratum random effect (BLUP) of an adjusted MAIHDA model, i.e. how far the stratum departs from the additive main-effects prediction of its defining dimensions – and flags the strata whose interaction is credibly different from zero. This is the heart of "where is there genuine intersectionality": a flagged stratum is one whose outcome departs credibly from what the additive parts predict – a descriptive statement about the stratum's outcome, not a causal claim about identity.

Usage

maihda_interactions(
  object,
  conf_level = 0.95,
  adjust = "BH",
  rope = NULL,
  scale = c("link", "response"),
  ...
)

Arguments

Details

It must be read off the adjusted model. Only when the dimensions' additive main effects are in the model (the adjusted model of the two-model decomposition, or the crossed-dimensions model) does the stratum random effect isolate the pure interaction. On a null model the stratum random effect is the total between-stratum deviation (additive + interaction), so passing one is flagged with a warning. The opposite mis-specification is flagged too: a model that adds a fixed interaction among the dimensions (e.g. var1 * var2) absorbs the intersectional effect into fixed cell means, so the stratum random effect is no longer the pure interaction. Passing a maihda result uses the right model automatically.

Frequentist vs. Bayesian evidence. For the frequentist engines (lme4, wemix, ordinal) the flag comes from the BLUP's conditional standard error: a Wald interval at conf_level and a two-sided p-value, with an optional multiplicity correction (adjust). For brms the full posterior is already available, so the exact posterior tail is used – a credible interval at conf_level and the probability of direction pd = max(P(BLUP > 0), P(BLUP < 0)) (in [0.5, 1]; the sign is in direction) – and adjust is not applied: the posterior is already partially pooled, the hierarchical-shrinkage answer to multiplicity (Gelman, Hill & Yajima 2012). That is a per-stratum answer; the marginal intervals carry no formal joint error-rate guarantee for the collective claim "an interaction exists somewhere" (the disjunction reading below applies to Bayesian flags too).

Multiplicity: partial pooling and a correction are different things, and the experts disagree.

• Shrinkage (magnitude/sign). The stratum BLUP is partially pooled, so extreme values are regularised toward the grand mean, attenuating exaggerated-magnitude and wrong-sign (Type M/S) error (Gelman & Carlin 2014). Gelman, Hill & Yajima (2012) argue this shrinkage usually substitutes for a classical multiple-comparisons correction (the problem can "disappear entirely" in the hierarchical model); on that view the flag/no-flag step itself is what to avoid – the null of an exactly zero interaction is rarely the question (McShane, Gelman et al. 2019) – so report the estimate and its interval.

• Whether to correct. If you do want an error-rate screen, whether a correction is warranted depends on the inferential structure of the claim. Each stratum as its own pre-specified hypothesis ("does this stratum interact?") is individual testing and needs none – only if you do not also read the flags collectively. Once the question is "is there an interaction somewhere?" – which an automated all-strata scan effectively is – it is disjunction testing and a correction applies.

adjust = "BH" is the default: fitting and flagging every stratum in one call is the disjunction/screening case, where controlling the expected proportion of false discoveries (FDR) is the appropriate goal. Pass adjust = "none" only when each stratum is a genuine, pre-specified individual hypothesis. The FDR choice (over family-wise "bonferroni"/"holm") is this package's, matching that screening goal. The flag itself is a Wald test on a shrunken BLUP whose conditional SE treats the variance components as known, so it (and any adjust on it) is an explicit, conservative screen – strict, not liberal. Partial pooling deflates a truly-null stratum's BLUP more than its conditional SE, so the null z-statistic is sub-normal (variance about the shrinkage fraction, below 1) and the screen under-flags rather than over-flags, degenerating to no flags at the singular/zero-variance boundary. Lead with the interval (and, for brms, the probability of direction); the substantive question is often not whether an interaction differs from zero but whether it exceeds a smallest interaction of interest (an equivalence reading; Schuirmann 1987; Kruschke 2018), read from the interval.

The interaction is reported on the model's link (latent) scale – a log-odds deviation for a logistic model, etc. – because the additive/interaction split is only exact there. On a non-identity link that makes it a multiplicative departure, which is a different quantity from a shift in the outcome itself; see the section below before reporting it.

Value

An object of class maihda_interactions (a data frame), one row per stratum, sorted flagged-first then by abs(interaction). Columns common to every engine: stratum, label, n (stratum size), interaction (the BLUP, under the default scale), lower/upper (the interval), flagged (logical), and direction ("above"/"below" the additive expectation). Under scale = "response" the interaction/lower/upper columns hold \pi^B_j and its interval, and se is dropped: the conditional standard error is a link-scale quantity, and the response-scale interval is deliberately not symmetric about the estimate, so no single SE would reproduce it. Frequentist fits add se (link scale only) and p_value (and p_adjusted when adjust != "none"). p_value is a conditional screening statistic – a Wald tail on the shrunken BLUP's conditional SE, with the variance components treated as known – not a calibrated frequentist p-value; it is conservative (stochastically large under a true null), so the BH-adjusted flag under-flags truly-null strata rather than delivering an exact error-rate guarantee (see Details). brms instead adds pd (probability of direction, max(P(>0), P(<0)) in [0.5, 1]). When rope is set, a decision column ("relevant"/"negligible"/"inconclusive") is added. Attributes record conf_level, adjust, rope, engine, model_type, n_strata, n_flagged, scale, link (the model's link name, so a caller can tell a log-odds departure from one in the outcome's own units), response_kind (what the model's response scale is – "probability", "count", "score" or "response" – which the link alone does not settle, since a cumulative fit is logit-linked but scores categories), singular (a variance of the fit is at zero, in any term), interaction_at_boundary (the variance of the term the interaction is read from is at zero, so its standard errors collapse; print() then says so), and – on a non-identity link under the default scale – response_interaction, the outcome-scale estimates keyed by stratum. That last attribute exists so print() can show the outcome-scale size beside the link-scale one; the columns are the same either way, and scale = "response" is how to obtain those numbers, with their interval, as data.

What the interaction means on a non-identity link

For a Gaussian (identity-link) fit the stratum BLUP u_j is the part of the stratum's mean outcome attributable to interaction, in the outcome's own units. For a logistic fit it is not: u_j is a deviation in log-odds, so what this function reports is the multiplicative, odds-scale interaction. (What follows is written for the logistic case, the one the literature works in. The same logic holds on any non-identity link, in its own units: a Poisson fit's u_j is a log-rate departure and scale = "response" gives it as a difference in expected counts, and a cumulative (ordinal) fit's is a latent shift reported as a difference in expected category score. The printed column is named for whichever applies.) Evans et al. (2024, section 2.5.1) put it directly – in a logistic MAIHDA one can "no longer directly interpret u_j as the change in mean outcomes (i.e., shift in probabilities) attributable to interaction effects".

Three quantities, of which this function reports two. Keeping them apart is the whole of the difficulty:

1. The multiplicative interaction, u_j – the departure from additivity in log-odds, i.e. whether a dimension multiplies the odds by the same factor at every level of the others. This is what the default scale = "link" reports and what flagged is about.

2. The same departure in outcome units, \pi^B_j – what scale = "response" returns. Evans et al. (2024) define it as

   \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),

   the gap between a stratum's total predicted probability and the probability implied by the additive main effects alone, and rank-plot \pi^B_j where the linear case plots u_j. This is quantity 1 re-expressed, not a second finding: it is zero exactly when u_j is zero and it flags the same strata, so it changes the units you report, not what you may conclude.

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 of the above reports this, and it is generally non-zero even where u_j is exactly zero, because the logistic curve is steeper in the middle than in the tails. With a -2 baseline and +0.7 for each of two dimensions and no interaction at all (u_j = 0 throughout), the second dimension still adds 9.5 percentage points of risk at one level of the first and 14.0 at the other. That excess is a property of the link, not a finding.

So "no strata flagged" supports "no credible multiplicative interaction", not "no interaction": additivity in log-odds does not carry over to probabilities, so quantity 3 is generally non-zero regardless, and it is often the one a policy audience cares about. (Generally, not always. The risk-difference interaction is positive below the curve's midpoint and negative above it, so it passes through zero: with two dimensions of effect a and b it vanishes exactly at an intercept of -(a + b)/2, where the two comparisons straddle the midpoint symmetrically. That is a single configuration, not the general case, and the flags say nothing about which one you are in either way.)

Both scales flag the same strata. Writing g(u) for the map from a stratum's BLUP to its \pi^B_j, g is strictly increasing with g(0) = 0, so estimate, interval endpoints and zero all carry across together. flagged, direction, p_value / p_adjusted and pd are therefore identical under either scale, and the response-scale interval is the exact image of the link-scale one rather than a delta-method approximation – inheriting its conditionality (fixed effects and variance components held at their point estimates), which is why no simulation step is needed. What does change is the size and the ranking: the same log-odds departure is worth more probability near \pi = 0.5 than in the tail, so strata are ordered by the quantity actually reported.

If quantity 3 is what you need, no argument here will give it to you – fit a linear probability model instead: fit_maihda(..., family = "gaussian") on a 0/1 outcome, which fit_maihda signposts 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).

References

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.

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., 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

Gelman, A., Hill, J., & Yajima, M. (2012). Why we (usually) don't have to worry about multiple comparisons. Journal of Research on Educational Effectiveness, 5(2), 189-211.

Gelman, A., & Carlin, J. (2014). Beyond power calculations: assessing Type S (sign) and Type M (magnitude) errors. Perspectives on Psychological Science, 9(6), 641-651.

Rubin, M. (2021). When to adjust alpha during multiple testing: a consideration of disjunction, conjunction, and individual testing. Synthese, 199(3-4), 10969-11000. doi:10.1007/s11229-021-03276-4

McShane, B. B., Gal, D., Gelman, A., Robert, C., & Tackett, J. L. (2019). Abandon statistical significance. The American Statistician, 73(sup1), 235-245.

Schuirmann, D. J. (1987). A comparison of the two one-sided tests procedure and the power approach for assessing the equivalence of average bioavailability. Journal of Pharmacokinetics and Biopharmaceutics, 15(6), 657-680.

Kruschke, J. K. (2018). Rejecting or accepting parameter values in Bayesian estimation. Advances in Methods and Practices in Psychological Science, 1(2), 270-280.

See Also

maihda, calculate_pcv, summary.maihda_model; and plot(..., highlight_interactions = TRUE) to mark the flagged strata on the effect-decomposition / predicted / shrinkage plots.

Examples

data(maihda_health_data)
a <- maihda(BMI ~ Age + Gender + Race + (1 | Gender:Race),
            data = maihda_health_data)
maihda_interactions(a)                  # FDR-screened (default adjust = "BH")
maihda_interactions(a, adjust = "none") # uncorrected per-stratum individual view
maihda_interactions(a, rope = 0.1)      # equivalence: |interaction| within 0.1?

# On a logistic fit the BLUP is a log-odds departure; scale = "response" reports
# Evans et al.'s (2024) pi_j - pi^A_j, the same interaction in probability points.
b <- maihda(Obese ~ Gender + Race + (1 | Gender:Race),
            data = maihda_health_data, family = "binomial")
maihda_interactions(b, scale = "response")

Simulated Longitudinal Data for MAIHDA

Description

A simulated long-format (one row per person-occasion) panel for demonstrating longitudinal / growth-curve MAIHDA (fit_maihda(id =, time =) and maihda(decomposition = "longitudinal")). 600 individuals are each measured over five waves on a continuous wellbeing score, within 12 intersectional strata defined by gender x ethnicity x education. The between-stratum trajectory differences are constructed to be mostly additive (each dimension's main effect on the baseline level and on the rate of change) with one genuine multiplicative interaction, so the longitudinal PCV – a high but sub-1 PCV_slope – is demonstrable.

Usage

maihda_long_data

Format

A data frame with 3000 rows (600 persons x 5 waves) and 8 variables:

id

Person identifier (level 2); repeated across waves.

wave

Measurement occasion, 0 to 4 (the numeric time variable).

gender

Gender (Women/Men); a stratum dimension.

ethnicity

Ethnicity (EthA/EthB/EthC); a stratum dimension.

education

Education (Low/High); a stratum dimension.

age

Baseline age in years, a time-invariant covariate.

wellbeing

Continuous wellbeing outcome (the growth-curve response).

low_wellbeing

Binary companion outcome (1 = bottom 40% of wellbeing), for the logistic longitudinal path.

Source

Simulated for the purpose of the MAIHDA package. The growth structure follows the longitudinal MAIHDA of Bell, Evans, Holman & Leckie (2024) doi:10.1016/j.socscimed.2024.116955.

Examples

data(maihda_long_data)

# Time-varying VPC from a 3-level growth model:
m <- fit_maihda(wellbeing ~ wave + (1 | gender:ethnicity:education),
                data = maihda_long_data, id = "id", time = "wave")
summary(m)

# Additive-vs-multiplicative PCV (null vs adjusted growth model):
a <- maihda(wellbeing ~ wave + (1 | gender:ethnicity:education),
            data = maihda_long_data, id = "id", time = "wave",
            decomposition = "longitudinal")
a$pcv

Build the adjusted-model formula for a longitudinal MAIHDA decomposition

Description

The longitudinal analogue of maihda_adjusted_formula: the null growth model plus the stratum dimensions' additive main effects AND their interactions with the time polynomial (dim:time), so the remaining stratum-level intercept/slope variance is the intersectional interaction beyond additive. Auto-binned numeric dimensions reuse their reconstructed tertile factor (.maihda_dim_*, via maihda_adjusted_terms).

Usage

maihda_longitudinal_adjusted_formula(
  null_formula,
  strata_vars,
  autobin_info,
  data,
  time,
  time_degree
)

Arguments

Value

A list with formula and data, or NULL if fewer than two dimensions are available.

Build the 3-level growth formula for a longitudinal MAIHDA

Description

Given a base formula already carrying the covariates and the resolved stratum grouping (y ~ covars + (1 | stratum)), returns the growth formula y ~ covars + time(+ I(time^2)...) + (time... | id) + (time... | stratum): the time polynomial enters the fixed part (if absent) and a random intercept+slope block is placed at both the individual and stratum levels. Any random effects in the base formula are replaced by this canonical structure.

Usage

maihda_longitudinal_formula(
  base_formula,
  id,
  time,
  time_degree,
  orig_time = time,
  stratum_slope = TRUE
)

Arguments

Details

When the growth terms are built on internally centered time (time is the derived .maihda_ctime column; see the file header), any bare raw-time polynomial the user wrote in the fixed part (orig_time, I(orig_time^2), ...) is replaced by the centered terms rather than kept alongside them, which would be perfectly collinear.

Value

The growth formula (same environment as base_formula).

Longitudinal MAIHDA proportional change in variance (PCV)

Description

Compares the stratum-level random-effect covariance block of the null growth model with that of the adjusted model (null + dimension main effects + their dim:time interactions). Reports the PCV in the baseline variance and in the instantaneous-slope variance at baseline – the additive-vs-multiplicative split of the intersectional trajectory inequality (Bell, Evans, Holman & Leckie 2024) – and the time-specific PCV over the supplied times. Both are evaluated at the observed baseline time (ref_time), so they are invariant to how the time axis is coded (for linear growth the slope variance is the same at every time, so this reduces to the slope-variance cell).

Usage

maihda_longitudinal_pcv(
  null_model,
  adjusted_model,
  times = NULL,
  estimation = c("fitted", "ML")
)

Arguments

Details

As in calculate_pcv, the estimation argument selects the variance-estimation basis. With estimation = "fitted" (the default) each fit's own REML covariance block is used, matching the single-model summaries. With estimation = "ML", REML lmer growth fits are refitted with maximum likelihood (refitML) before the comparison, which overstates both PCVs when strata are few. The stored models (and the single-model summaries computed from them, e.g. the time-varying VPC) always keep their REML fit; ml_refit on the result records whether an ML refit fully applied, and estimation_used records the basis ACTUALLY used – "fitted", "ML", or "mixed" when an ML refit was requested but a model kept its REML fit (a boundary skip or a failed refit), so a mixed REML/ML comparison is not mistaken for a clean one. glmer (GLMM) and brms fits are already on the ML / posterior scale, so the choice does not affect them.

A variance estimate cannot fall below zero. When the strata differ little after adjustment next to the precision of one stratum's own baseline or slope, the adjusted model's variance is too large on average and the null model's is not, so the PCV is understated on average and single values can be far too low or negative. This affects the slope first, where strata usually differ least, and eases with more strata, larger strata and more measurement occasions.

For a brms pair every variance compared is the posterior median of its per-draw value: the between-stratum variance a(t)' \Sigma a(t) at the baseline and at each reporting time, and the slope variance. The stored blocks Sigma_stratum_null and Sigma_stratum_adjusted are the posterior means summary() tabulates, so those variances are not quadratic forms of them. The posterior of a between-stratum variance is right-skewed when strata are few, and more so in the adjusted model, so the PCV of a brms pair is understated on average with few strata – by more when the additive share is moderate, when the strata are small, or, for the slope, when the strata differ little in slope next to the precision of a stratum's slope – and can be negative where the true value is positive. The lme4 (REML) decomposition is less biased then; neither is precise with so few strata.

Value

An object of class maihda_long_pcv. variance_summary names what var_baseline_*, var_slope_* and the variances in pcv_t are ("estimate" for lme4, "posterior median" for brms), and Sigma_summary what the two stored blocks are ("estimate" or "posterior mean"). ml_variances is TRUE when the compared variances are maximum-likelihood estimates: a glmer growth pair, or lmer growth fits made or refitted by maximum likelihood.

Per-stratum trajectory parameters for a longitudinal MAIHDA

Description

The stratum-level random-effect estimates as a wide table, one row per stratum: the stratum's deviation at the baseline time (baseline, the longitudinal analogue of a cross-sectional stratum BLUP), the random intercept (intercept; the deviation at the model's coefficient origin – the internal centering offset, i.e. the observed baseline, when centering applied, or raw time 0 otherwise) and the random slope(s) on time (slope, ...). This is the longitudinal shape of predict_maihda(type = "strata") – a stratum is now a trajectory, not a single value.

Usage

maihda_longitudinal_strata_predictions(object)

Arguments

Details

baseline is a(t_0 - c)' coef with a(t) = (1, t, t^2, ...), t_0 = the reference (baseline) time ref_time = min(time) and c the internal centering offset; it equals intercept whenever ref_time coincides with the centering origin (the usual case for a centered fit, and the zero-anchored case for a raw fit).

Value

A data frame: stratum, stratum_id, optional label, baseline, intercept, slope(, slope2, ...).

Time-varying VPC summary for a longitudinal MAIHDA (brms, linear growth)

Description

Time-varying VPC summary for a longitudinal MAIHDA (brms, linear growth)

Usage

maihda_longitudinal_summary_brms(object, conf_level = 0.95)

Arguments

Value

As maihda_longitudinal_summary_lme4, with posterior bands.

Time-varying VPC summary for a longitudinal MAIHDA (lme4)

Description

Time-varying VPC summary for a longitudinal MAIHDA (lme4)

Usage

maihda_longitudinal_summary_lme4(
  object,
  bootstrap = FALSE,
  n_boot = 1000,
  conf_level = 0.95
)

Arguments

Value

A list with vpc_result (the reference-time VPC, for the headline print), variance_components, and longitudinal (the trajectory).

Median Odds Ratio (MOR) for a logistic MAIHDA model

Description

The Median Odds Ratio translates the between-stratum variance of a logistic MAIHDA model onto the odds-ratio scale: the median relative change in the odds of the outcome when comparing two individuals from randomly chosen strata (higher- vs lower-risk). MOR = exp(sqrt(2 * V_A) * qnorm(0.75)), where V_A is the between-stratum (latent, logit-scale) variance. An MOR of 1 indicates no between-stratum heterogeneity. The MOR is defined only for the logit link (it is the median odds ratio); a non-logit binomial fit such as probit is rejected, because its latent variance is on a different scale and the exp(...) above would not be an odds ratio.

For a cumulative-logit (ordinal) MAIHDA model the same formula applies to the latent logit-scale between-stratum variance and is the median cumulative odds ratio: the median relative change in the odds of being at or below any given outcome category between two randomly chosen strata (under the model's proportional-odds assumption it is the same for every category split).

Scope. V_A is the between-stratum variance. Contextual (context = ) and other non-stratum random effects are never included: the MOR compares two individuals from randomly chosen strata within the same context.

Crossed-dimensions fits use a different calculation. The closed form above assumes the two strata's random effects are independent, which is what makes their difference N(0, 2 V_A). In a crossed-dimensions fit (from maihda(decomposition = "crossed-dimensions")) a stratum's effect is the sum of its dimension effects plus the intersection effect, so two strata sharing a dimension – say two strata that are both "female" – share that dimension's random effect and are correlated. Their difference is then a mixture of normals, one component per pattern of shared dimensions, and applying the closed form to the summed variance overstates the MOR (substantially so when the variance sits mainly in the additive dimensions, and not at all when it sits entirely in the interaction).

The MOR reported for such a fit is therefore computed from that mixture, under an explicit sampling scheme: two distinct strata drawn uniformly at random from the strata present in the analytic sample. Writing \tau^2_d for dimension d's variance and \tau^2_I for the intersection variance, a pair differing on the dimension set D^* has difference variance v = 2(\tau^2_I + \sum_{d \in D^*} \tau^2_d), and the MOR is exp(x) for the x solving \sum_{pairs} (2\Phi(x/\sqrt{v}) - 1) / n_{pairs} = 0.5. For a canonical single-stratum fit the two calculations coincide, and that closed form is used.

brms fits. For a brms fit V_A is the posterior median of the between-stratum variance, the value in the variance-components table of summary.maihda_model. The closed form is monotone in V_A, so the MOR is then the posterior median of the Median Odds Ratio itself. A crossed-dimensions fit evaluates the mixture at the posterior median of each variance. With few strata the posterior of a variance is right-skewed and its median overstates the variance, so the MOR runs high; the overstatement shrinks as strata are added.

Usage

maihda_mor(model)

Arguments

Value

A single number (the MOR, \ge 1), or NA_real_ if the between-stratum variance is unavailable – which for a crossed-dimensions fit also covers the case where the stratum grid needed for the mixture cannot be resolved (fewer than two strata, absent dimension columns, or more than 12 dimensions).

References

Larsen, K., & Merlo, J. (2005). Appropriate assessment of neighborhood effects on individual health: integrating random and fixed effects in multilevel logistic regression. American Journal of Epidemiology, 161(1), 81-88.

See Also

maihda_discriminatory_accuracy

Category probabilities of a cumulative model

Description

Pure function (no fit object): given latent locations eta, ordered thresholds alpha and the link, returns the category probability matrix via P(Y \le k) = g^{-1}(\alpha_k - \eta) and differencing. Rows are observations, columns categories 1..K (K = length(alpha) + 1).

Usage

maihda_ordinal_category_probs(eta, thresholds, link = "logit")

Arguments

Value

A numeric matrix with length(eta) rows that sum to 1.

Latent location to expected category score

Description

Convenience composition of maihda_ordinal_category_probs and maihda_ordinal_expected_score.

Usage

maihda_ordinal_eta_to_score(eta, thresholds, link = "logit")

Arguments

Value

A numeric vector of expected category scores.

Expected category score from a probability matrix

Description

The response-scale summary of a cumulative model used throughout the package (the plot layer's "Average Expected Category Score"): \sum_k k\, p_k, with categories scored 1..K in order.

Usage

maihda_ordinal_expected_score(probs)

Arguments

Value

A numeric vector of expected scores in [1, K].

Prepare data and formula for a sampling-weighted brms fit

Description

Drops rows with missing or non-positive sampling weights and rows incomplete on the model variables (each with a warning), normalizes the remaining weights to mean 1 – likelihood weights scale the effective sample size, so unnormalized expansion weights (summing to the population) would massively overstate the information in the data – and rewrites the formula with a weights() addition term.

Usage

maihda_prepare_brms_sampling_weights(data, formula, sampling_weights)

Arguments

Details

Incomplete rows must go BEFORE the normalization: brms silently excludes them after receiving the data, so normalizing over all weight-valid rows would hand the sampler surviving weights with an arbitrary mean, scaling the pseudo-posterior's effective sample size by that mean.

Value

A list with data (weights column .maihda_sw added), formula (rewritten), and keep (a logical mask over the input rows marking those retained, so the caller can re-slice any row-aligned forwarded arguments to the same rows).

Resolve shared strata and the fitting formula for group comparison

Description

Resolve shared strata and the fitting formula for group comparison

Usage

maihda_prepare_group_strata(
  formula,
  data,
  shared_strata,
  autobin = TRUE,
  bin_rows = NULL
)

Arguments

Value

list(data, formula, strata_vars).

Print the additive vs. intersectional decomposition of a crossed-dimensions summary

Description

Print the additive vs. intersectional decomposition of a crossed-dimensions summary

Usage

maihda_print_cc_decomposition(d)

Arguments

Value

No return value, called for side effects.

Print the stratum vs. context partition of a contextual cross-classified summary

Description

Print the stratum vs. context partition of a contextual cross-classified summary

Usage

maihda_print_context_partition(ctx)

Arguments

Value

No return value, called for side effects.

Parametric-bootstrap proportional-odds test for a cumulative MAIHDA fit

Description

Tests the proportional-odds (parallel-lines) assumption of a fitted cumulative clmm MAIHDA model by calibrating the omnibus nominal-effects likelihood-ratio statistic against its own null distribution under the fitted model.

Usage

maihda_proportional_odds_test(object, n_sim = 199, seed = NULL)

Arguments

Details

The statistic is the omnibus nominal-effects likelihood-ratio statistic: the fixed-effect part of the model is refitted twice with ordinal::clm() – once with all covariates proportional, once with every covariate entering as a threshold-specific (nominal) effect – and twice the log-likelihood difference is taken.

Its usual chi-squared reference (ordinal::nominal_test()) is not valid here. After the stratum random intercept is integrated out, the marginal cumulative-logit slopes generally differ across thresholds once the stratum variance is non-zero. The chi-squared reference also treats the observations as independent, whereas observations that share a stratum share its random effect. A correctly specified model is therefore rejected too often, and more often as the sample grows.

The null distribution is instead simulated under the fitted clmm: each replicate redraws the stratum random effects from N(0, \tau^2) at the fitted variance, forms the conditional category probabilities from the fitted thresholds and location predictor, redraws the response, and recomputes the statistic. The p-value is (1 + \#\{T_b \ge T_{obs}\}) / (1 + B) over the B replicates that refitted. n_sim sets only its Monte Carlo resolution: the smallest attainable value is 1 / (B + 1).

This is a parametric-bootstrap approximation, not an exact test. The replicates are drawn at the fitted stratum variance, which is estimated from only as many units as there are strata, so the test is least reliable when the strata are few, whatever n_sim is.

The test is expensive – every replicate refits two clm() models – and is not run at fit time.

Value

An object of class maihda_po_test: a list with lrt, df, n_terms, p_value (the bootstrap p-value), p_chisq, n_sim (replicates that produced a usable statistic), n_failed, and null_lrt (the simulated null statistics). p_value is the one to report and the only one printed. p_chisq is the chi-squared p-value of the same statistic; it is not calibrated for a mixed model (see Details) and is returned for comparison only.

See Also

fit_maihda, maihda_cumulative

Examples

strata <- make_strata(maihda_sim_data, vars = c("gender", "race"))
d <- strata$data
d$y <- factor(cut(d$health_outcome, 3), labels = 1:3, ordered = TRUE)
m <- fit_maihda(y ~ age + (1 | stratum), data = d, family = "ordinal")
maihda_proportional_odds_test(m, n_sim = 99, seed = 1)

Simulated Health Data for MAIHDA Use

Description

A simulated dataset for demonstrating Multilevel Analysis of Individual Heterogeneity and Discriminatory Accuracy (MAIHDA).

Usage

maihda_sim_data

Format

A data frame with 500 rows and 7 variables:

id

Unique participant identifier.

gender

Gender of the participant.

race

Simulated race/ethnicity category.

education

Educational attainment level.

age

Age in years, a continuous covariate.

health_outcome

A continuous simulated health outcome.

binary_outcome

A binary version of the health outcome.

Source

Simulated for the purpose of the MAIHDA package.

Examples

data(maihda_sim_data)
strata_result <- make_strata(maihda_sim_data, vars = c("gender", "race", "education"))

Sparse Intersectional Data for Bayesian MAIHDA

Description

A simulated cross-sectional dataset in which many intersectional strata hold only a handful of individuals, built with a known intersectional interaction. With strata this small the split of the between-stratum variance into additive and interaction parts is weakly identified. An lme4 fit of the crossed-dimensions decomposition can be singular; a brms fit of it (engine = "brms") gives the split a posterior median and a credible interval, both of which depend on the prior on the standard deviations.

Usage

maihda_sparse_data

Format

A data frame with 240 rows and 6 variables:

gender

Strata dimension (Women/Men).

ethnicity

Strata dimension (White/Black/Asian).

education

Strata dimension (Low/High).

age_group

Strata dimension (Young/Mid/Older).

y

A continuous (Gaussian) outcome with residual SD 1. The 36 stratum effects have a variance of 0.31 (a VPC of 0.24), of which 41% is the intersectional interaction.

event

A binary outcome (No/Yes), ~46% "Yes". On the latent scale its stratum effects have a variance of 1.32 (a VPC of 0.29), again 41% interaction.

attr(maihda_sparse_data, "truth") holds, for each outcome, the intercept and the 36 stratum effects with their additive, interaction and between-stratum variances, the interaction share and the VPC.

Details

Four dimensions form 36 intersectional strata with deliberately skewed sizes (median 6 individuals, 12 of 36 cells below 5, two singletons). The stratum effects were drawn with the interaction aimed at 40% of the between-stratum variance; the 36 effects drawn carry 41% on both outcomes (on the latent scale for event), by the mean-square split of the section Realised strata of maihda. The interaction share a crossed-dimensions fit reports is a ratio of variance parameters, and the four dimensions' variances, each estimated from two or three levels, add up to more than the spread of the main effects among the 36 strata: fitted to these stratum effects with 1,000 individuals in every stratum, lme4 puts the share at about 28% on y and 35% on event, the gap between the two being that of REML and maximum likelihood (the two sets of effects differ only in scale). In vignette("bayesian_sparse_maihda") both lme4 fits are singular through a dimension's variance, not the interaction's, and put the share at about 23% on y and 3% on event. The brms posterior medians are below 41% as well; under the priors the vignette sets on the standard deviations, their credible intervals run from near zero to about 41% or more.

Note

A purely illustrative dataset. The dimension labels are arbitrary and the interaction is constructed, not estimated from any real population – its only purpose is to make the sparse-cell behaviour of the ML and Bayesian estimators visible against a known answer.

Source

Simulated; see data-raw/maihda_sparse_data.R.

Examples

data(maihda_sparse_data)
attr(maihda_sparse_data, "truth")$gaussian$interaction_share  # 0.41

# lme4: the interaction share as a point estimate, from a fit that is singular here:
# m_lme4 <- maihda(y ~ 1 + (1 | gender:ethnicity:education:age_group),
#                  data = maihda_sparse_data, decomposition = "crossed-dimensions")
#
# brms: a posterior median and a credible interval, both resting on the prior:
# m_brms <- maihda(y ~ 1 + (1 | gender:ethnicity:education:age_group),
#                  data = maihda_sparse_data, decomposition = "crossed-dimensions",
#                  engine = "brms",
#                  prior = brms::set_prior("normal(0, 0.5)", class = "sd"))
# See vignette("bayesian_sparse_maihda").

Per-stratum predictions for a cumulative (clmm) fit

Description

Ordinal counterpart of maihda_stratum_predictions_wemix(): per-stratum aggregates of the location prediction plus the stratum effect, on the latent (link) scale or as the expected category score (response scale).

Usage

maihda_stratum_predictions_ordinal(
  object,
  summary_obj,
  scale = c("response", "link")
)

Arguments

Value

A data frame as from maihda_weighted_stratum_aggregate().

Per-stratum predictions for a wemix fit

Description

wemix counterpart of maihda_stratum_predictions_lme4(): per-stratum means of the fixed-part prediction plus the stratum effect, aggregated with the SAMPLING weights so the stratum-level summaries are design-weighted (population-representative under the weights), unlike the lme4 prior-weight aggregation.

Usage

maihda_stratum_predictions_wemix(
  object,
  summary_obj,
  scale = c("response", "link")
)

Arguments

Value

A data frame as from maihda_weighted_stratum_aggregate().

Remove the package-internal weights(.maihda_sw) addition term from a formula

Description

The inverse of the injection maihda_brms_weights_formula() performs for the reserved .maihda_sw column. maihda() and compare_maihda_groups() derive the null / adjusted models from a sampling-weighted fit's stored formula, which already carries the injected weights(.maihda_sw) term; that term must be stripped before the derived model is re-prepared, or the reserved-column guard (.maihda_sw both in the formula and in the carried-over data) and the weights-formula rewrite (which rejects an existing weights() term) abort the package's own refit. Only the internal weights(.maihda_sw) call is removed – a user's own weights() term (a genuine conflict) is left in place so it is still rejected – and any other addition term (e.g. trials(n)) is preserved.

Usage

maihda_strip_brms_weights_term(formula)

Arguments

Value

The formula with the internal weights(.maihda_sw) term removed, or the input unchanged when it carries no such term.

Canonical MAIHDA results table and ranked-strata table

Description

Assembles the two standard MAIHDA write-up deliverables from a fitted analysis in one call: (a) a model-results table contrasting the null and adjusted models (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 the intersectional strata by their predicted outcome, so the best- and worst-off strata can be read directly. It introduces no new estimator: the model-results table reuses the quantities from summary() (calling summary() itself for a bare fit_maihda model), and the ranked-strata table reuses the same stratum predictions as plot(type = "predicted"), so the table agrees exactly with summary() and plot().

Usage

maihda_table(
  x,
  n_strata = 10L,
  scale = c("response", "link"),
  which = c("null", "adjusted"),
  digits = 3,
  ...
)

Arguments

Details

The model-results table is mostly numeric and export-ready (e.g. write.csv(maihda_table(a)$models, ...) or pass it to knitr::kable()): statistics are rows, models are columns, and each estimate has accompanying *_lower/*_upper columns that hold the confidence/credible interval when one is available (the intercept's interval, the VPC bootstrap or posterior interval, and the bootstrap PCV interval) and NA otherwise. The intercept carries the same interval summary() reports for it – a Wald interval for the likelihood engines, a credible interval for brms. The variance and SD rows are point estimates. The print() method renders the same table in the familiar “estimate [low, high]” layout.

For a "crossed-dimensions" analysis (one model, no null/adjusted pair) the results table has a single estimate column and gains “Additive share” / “Interaction share” rows instead of the PCV. For a contextual cross-classified analysis (maihda(context = )) it gains a “Context share (VPC)” row. A bare fit_maihda model is also accepted and yields a single-model table (no PCV).

The ranked-strata table ranks every stratum by its model-predicted outcome (on the scale requested), using the same stratum predictions as plot(type = "predicted"): the predicted value carries the conditional (random-effect) interval, and the stratum random effect (BLUP) is reported alongside it. By default the ranking uses the null model – the headline intersectional inequality (which strata fare best/worst overall); set which = "adjusted" to rank by the adjusted model instead. The full ranked table is returned in $strata; print() shows the top and bottom n_strata.

Value

An object of class maihda_table: a list with

See Also

maihda, summary.maihda_model, calculate_pcv, maihda_discriminatory_accuracy.

Examples

data(maihda_health_data)
a <- maihda(BMI ~ Age + Gender + Race + (1 | Gender:Race), data = maihda_health_data)

tab <- maihda_table(a)
tab                 # printed: model-results table + top/bottom strata
tab$models          # the numeric, export-ready results table
tab$strata          # all strata ranked by predicted BMI

# write.csv(tab$models, "results.csv", row.names = FALSE)

Tag a summary with its role in a two-model analysis

Description

Internal helper. The null and adjusted summaries of a maihda analysis are indistinguishable once pulled out of the analysis object, which makes a printed summary easy to misread – the null model's fixed effects are the intercept and covariates only, because the strata dimensions are its random-effect grouping rather than fixed effects. Stamping the role lets print.maihda_summary say which model it is showing.

Usage

maihda_tag_role(s, role)

Arguments

Value

s with a "maihda_role" attribute (NULL in, NULL out).

Tidy a MAIHDA summary, model, or analysis

Description

tidy methods that return the MAIHDA estimates as a tidy tibble, ready for downstream tables (gt, flextable) and ggplot2. They read the slots that summary.maihda_model already computes and add no new statistics.

Usage

## S3 method for class 'maihda_summary'
tidy(x, component = c("strata", "variance", "fixed"), ...)

## S3 method for class 'maihda_model'
tidy(x, component = c("strata", "variance", "fixed"), ...)

## S3 method for class 'maihda_analysis'
tidy(
  x,
  component = c("strata", "variance", "fixed"),
  which = c("null", "adjusted"),
  ...
)

Arguments

Value

A tibble. For component = "strata": columns stratum, label, estimate, std.error, conf.low, conf.high. For "variance": component, variance, sd, proportion. For "fixed": term, estimate, std.error, statistic, df, p.value, conf.low, conf.high.

Fixed-effect statistics

For the lme4, WeMix and ordinal engines statistic is the Wald statistic (estimate / std.error), p.value its two-sided p-value, conf.low/conf.high the matching Wald interval at the summary's conf_level, and df the denominator degrees of freedom they were built on.

A Gaussian lme4 fit reports containment (between-within) df and a t; every other engine leaves df NA and uses a z. Under summary(df_method = "bootstrap") p.value and the interval come from the parametric bootstrap and df is NA. See summary.maihda_model. For brms the estimate is the posterior mean with its Est.Error and credible interval, and statistic/df/p.value are NA. Standard errors are NA where they are undefined (a boundary fit).

See Also

glance.maihda_analysis for the one-row model summary.

Examples

data("maihda_health_data")
m <- fit_maihda(BMI ~ Age + (1 | Gender:Race:Education), data = maihda_health_data)
tidy(m)                       # stratum estimates
tidy(m, component = "variance")
tidy(m, component = "fixed")

Recommended Figure Size for the UpSet Stratum Plot

Description

Computes sensible width and height (in inches) for plot(object, type = "upset"), so a knitr chunk or ggsave() call can size the figure to its content. The UpSet composite grows taller with the number of matrix rows (one per binary 0/1 dimension, one per level of a multi-level factor) and wider with the number of strata columns shown, so a single fixed size tends to crop or stretch it – particularly for multi-level designs (many rows) or a large n_strata (many columns; UpSet is an inherently wide format).

Usage

maihda_upset_size(object, n_strata = 50)

Arguments

Value

A list with numeric width and height (inches) plus the rows (matrix rows) and cols (strata shown) they derive from.

See Also

plot.maihda_model

Examples

strata <- make_strata(maihda_sim_data, vars = c("gender", "race"))
model <- fit_maihda(health_outcome ~ (1 | stratum), data = strata$data)
sz <- maihda_upset_size(model, n_strata = 30)
ggplot2::ggsave(
  tempfile(fileext = ".png"),
  plot(model, type = "upset", n_strata = 30),
  width = sz$width, height = sz$height)

Validate a longitudinal (id / time) specification

Description

Validate a longitudinal (id / time) specification

Usage

maihda_validate_longitudinal(
  id,
  time,
  time_degree,
  data,
  engine = "lme4",
  sampling_weights = NULL,
  context = NULL,
  formula = NULL,
  stratum_slope = TRUE,
  family = "gaussian"
)

Arguments

Value

A list list(id, time, time_degree, stratum_slope).

Validate a sampling-weights specification

Description

Validate a sampling-weights specification

Usage

maihda_validate_sampling_weights(sampling_weights, data)

Arguments

Value

The validated column name.

Response-scale VPC for a binomial or count MAIHDA model

Description

Computes the variance partition coefficient on the scale of the observed response: the probability scale for a binomial MAIHDA model, the count scale for a Poisson or negative-binomial one.

For a binomial model it uses the simulation method of Goldstein, Browne & Rasbash (2002). Stratum random effects u \sim N(0, \sigma^2_u) are simulated and converted to predicted probabilities p = g^{-1}(\eta + u) (with \eta the fixed-part linear predictor); the VPC is then the between-stratum variance of p as a share of the total (between + the binomial within-stratum variance \overline{p(1-p)}). Unlike the latent-scale VPC (fixed level-1 variance \pi^2/3 for the logit), the response-scale VPC depends on the overall outcome prevalence, so report it as a complement to – not a replacement for – the latent-scale value.

For a Poisson or negative-binomial model with the log link it evaluates the exact expressions of Leckie, Browne, Goldstein, Merlo & Austin (2020). With \lambda = \exp(\eta + \sigma^2_u / 2) the marginal expected count, the between-stratum variance of the counts is \lambda^2 (\exp(\sigma^2_u) - 1) and their within-stratum variance is \lambda + \alpha \lambda^2 \exp(\sigma^2_u), where \alpha = 1 / \theta is the negative-binomial overdispersion (0 for a Poisson model); the VPC is the between-stratum variance as a share of their sum. It is the share of the variance of the counts themselves that lies between strata, and it involves no approximation of the level-1 variance. The latent-scale VPC of summary.maihda_model refers to the log scale instead and depends on the count_approximation of fit_maihda, so the two answer different questions. Report the count-scale VPC when the partition should describe the observed counts, or when the expected count is low and the latent-scale approximations disagree.

Usage

maihda_vpc_response(model, n_sim = 10000, seed = NULL)

Arguments

Details

Binomial models. The fixed part \eta is collapsed to a single value – the mean linear predictor \bar\eta over the analytic sample – before the random effect is simulated around it. The result is therefore a VPC evaluated at the mean covariate profile (a conditional-at-mean estimate), not one marginalised over the empirical covariate distribution. For the canonical strata-only (null) model \eta is constant (the intercept), so the two coincide and the value is exact. For an adjusted model (one with covariates) they can differ, because the inverse link is nonlinear and g^{-1}(\bar\eta) \neq \overline{g^{-1}(\eta)}: read the response-scale VPC from the null model, or interpret an adjusted value as conditional on the average covariate profile rather than as a covariate-averaged (marginal) VPC.

The method is binomial-link agnostic: it maps the simulated stratum effects through whichever inverse link the model uses (logit, probit, cloglog, ...), so a non-logit binomial fit is computed on its own scale rather than rejected. Only the family is required to be binomial.

When the model carries random intercepts beyond the stratum (a contextual (1 | school) or an explicit (1 | site)), the simulation integrates over them: the reported estimate is the stratum share Var(E[p \mid u_{stratum}]) of the total response-scale variance, where the total includes the variation the other random effects induce in p plus the binomial within-variance. Simulating the stratum effect alone would overstate the stratum share. The inner integral over the combined non-stratum effect is evaluated by deterministic Gauss-Hermite quadrature, not an inner Monte-Carlo sample, so it contributes no simulation error of its own and increasing n_sim converges the whole estimate (the outer stratum draw is then the only Monte-Carlo dimension). The number of quadrature nodes used is reported in inner_nodes.

Count models. \lambda is computed for each observation of the analytic sample from its fixed-part linear predictor, any offset included, so with covariates or an offset the VPC differs between observations. The estimate is the mean of the observation-level VPCs; their median, quartiles and range are returned with it, as are the mean and range of \lambda. Prior weights of the fit enter these summaries as frequencies. A null model without an offset has a single \lambda and a single VPC. n_sim and seed are not used.

The expressions are those of a model whose only random effect is the stratum random intercept. A count model with a context or any other additional random effect, a crossed-dimensions fit, a longitudinal fit, and a link other than the log are refused.

Value

An object of class maihda_vpc_response: a list with estimate, scale = "response", method ("simulation" for a binomial model, "closed-form" for a count model), var_between (the latent-scale between-stratum variance) and var_other (the summed latent-scale variance of any non-stratum random intercepts, 0 when there are none).

For a binomial model it also holds n_sim, lp_fixed (the mean fixed-part linear predictor), inner_method ("gauss-hermite" when non-stratum effects are integrated out, else "none") and inner_nodes (the number of Gauss-Hermite quadrature nodes used for that inner integral, NA when there are no non-stratum effects).

For a count model it also holds family, the median, quartiles and range of the observation-level VPCs (all equal to estimate when \lambda is the same for every observation), lambda (the mean marginal expected count), lambda_range and theta (the negative-binomial dispersion, Inf for a Poisson model).

References

Goldstein, H., Browne, W., & Rasbash, J. (2002). Partitioning variation in multilevel models. Understanding Statistics, 1(4), 223-231.

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

See Also

maihda_discriminatory_accuracy, summary.maihda_model

Examples

strata <- make_strata(maihda_health_data, vars = c("Gender", "Race"))
d <- maihda_health_data
d$stratum <- strata$data$stratum
m <- fit_maihda(Obese ~ (1 | stratum), data = d, family = "binomial")
maihda_vpc_response(m, seed = 1)

# A count outcome: the VPC on the count scale, in closed form
set.seed(1)
u <- rnorm(max(d$stratum), sd = 0.4)
d$visits <- rpois(nrow(d), exp(0.5 + u[d$stratum]))
m_count <- fit_maihda(visits ~ (1 | stratum), data = d, family = "poisson")
maihda_vpc_response(m_count)

Warn when a numeric stratum dimension enters the adjusted model as a linear term

Description

In MAIHDA the stratum-defining dimensions are categorical, but the adjusted / decomposition model adds an un-binned numeric dimension by its raw column name – i.e. as a single linear slope, not categorical main effects (see the else branch of maihda_adjusted_terms). That silently changes the PCV/decomposition interpretation. This helper flags the offending dimensions and warns; it does not alter the model (the term is still entered as-is). The flagged case is a stratum variable that is numeric in the data, not a key in autobin_info, and has three or more distinct values – either a category code with few unique values, or a many-valued numeric fitted with autobin = FALSE. Categorical (factor/character) dimensions and auto-binned numerics – reconstructed as the .maihda_dim_* tertile factors – are excluded, and so is a binary (two-level) numeric dimension: with only two levels a linear term and the two-level factor span the same design space, so the fit and the PCV coincide and there is nothing to warn about.

Usage

maihda_warn_linear_strata_dims(strata_vars, autobin_info, data, fn = "maihda")

Arguments

Details

Called once per decomposition at each public entry point (not inside maihda_adjusted_terms(), which runs several times per decomposition), so the warning fires exactly once.

Value

Invisibly, the character vector of flagged dimension names (empty if none).

Warn when a MAIHDA fit's fixed part cannot represent the outcome's mean

Description

A MAIHDA reads the stratum random intercept's variance as between-stratum inequality. When the fixed design does not span the intercept – a 0 + or - 1 formula whose surviving terms are all numeric, so no factor is coded as cell means – there is nothing in the fixed part to carry the outcome's mean and the random intercept absorbs it instead. The between-stratum variance, and with it the VPC, the MOR and any PCV computed from the fit, then measure the location of the outcome rather than the strata.

Usage

maihda_warn_no_grand_mean(formula, data_list)

Arguments

Details

Silent unless the span test can be run and provably fails, so a formula whose design could not be built is never warned about on the strength of a written 0 +.

Value

Invisibly, TRUE if a warning was emitted.

Fixed-part (and optionally full) linear predictor of a wemix fit

Description

WeMix's own predict() method needs the grouping structure re-resolved and offers no fixed-only form, so predictions are built directly from the coefficient vector and the stored stratum effects: the fixed design matrix is constructed with the fit's factor coding – the levels it fitted and the contrast matrix it applied to each factor, whatever options(contrasts = ) says now – AND its transformation basis (so a data-dependent term such as scale(x) uses the fit's centre and scale rather than recomputing them from newdata) and multiplied by coef, any formula offset term is evaluated on newdata and added (only a fit saved before fit_maihda() refused offsets for this engine carries one; WeMix fitted it without the offset), and include_re adds each row's stratum effect (conditional mode; an unseen stratum contributes 0 – the zero-effect fallback that predict_maihda only reaches when allow_new_levels = TRUE, having otherwise rejected unseen strata upstream). Everything is on the link scale. A fit made with WeMix's internal centring (center_grand or center_group, which fit_maihda() now rejects) is refused: its coefficients are on a centred scale whose constants WeMix does not keep.

Usage

maihda_wemix_linpred(object, newdata = NULL, include_re = TRUE)

Arguments

Value

A numeric vector of link-scale predictions.

Stratum random-effect table for a wemix fit

Description

Mirrors maihda_stratum_ranef_lme4(): one row per stratum with the random-effect estimate (conditional mode), a conditional standard error, and a 95% interval. WeMix reports no conditional variances, so the SE is computed analytically from the weighted pseudo-likelihood: for a Gaussian model the conditional precision of u_j is 1/\tau^2 + \sum_j w_{ij}/\sigma^2 (the design-weighted analogue of lme4's condVar, to which it reduces at unit weights), and for a binomial-logit model the Laplace curvature at the conditional mode, 1/\tau^2 + \sum_j w_{ij}\,\hat p_{ij}(1-\hat p_{ij}). These are model-based approximations, not design-based (replicate-weight) uncertainty.

Usage

maihda_wemix_stratum_ranef(object)

Arguments

Value

A data frame with stratum, stratum_id, random_effect, se, lower_95, upper_95.

Variance components of a wemix MAIHDA fit

Description

Reads the between-stratum variance (and, for a linear model, the residual variance) from the WeMixResults variance table. For a binomial-logit model the level-1 variance is the usual latent-scale \pi^2/3, matching the lme4/brms summaries.

Usage

maihda_wemix_variances(object)

Arguments

Value

A list with stratum and residual variances.

Create Strata from Multiple Variables

Description

This function creates strata (intersectional categories) from multiple categorical variables in a dataset.

Usage

make_strata(
  data,
  vars,
  sep = " × ",
  min_n = 1,
  autobin = TRUE,
  bin_rows = NULL
)

Arguments

Details

If any of the specified variables has a missing value (NA) for a given observation, that observation will be assigned to the NA stratum (stratum = NA), rather than creating a stratum that includes the missing value.

The strata_info data frame is also attached as an attribute to the data, which allows fit_maihda() to automatically capture stratum labels for use in plots and summaries.

When autobin discretises a numeric grouping variable v, the adjusted-model and prediction machinery later add an internal factor column named .maihda_dim_<v>; the .maihda_dim_ prefix is therefore reserved. make_strata() errors if data already holds the .maihda_dim_<v> column for a variable it is about to auto-bin, so an existing user column is never silently overwritten (rename it, or pass autobin = FALSE).

A numeric grouping variable with 10 or fewer unique values (or any numeric when autobin = FALSE) is kept as its raw values: distinct values still define distinct strata, but the variable is not binned. If such a variable is then used to build an adjusted MAIHDA (maihda), it enters the adjusted model as a linear fixed effect rather than categorical main effects, and maihda() warns (when it has three or more distinct values). Wrap category codes in factor() before creating strata if you want them treated categorically.

Value

An object of class maihda_strata: a list with elements

Examples

# Create strata from gender and race variables
result <- make_strata(maihda_sim_data, vars = c("gender", "race"))
print(result$strata_info)

Order-Invariant PCV Attribution Across Predictors (Shapley / Dominance)

Description

Apportions the between-stratum variance reduction (the PCV) among a set of predictors. Where stepwise_pcv adds the variables one at a time – so each variable's contribution depends on its entry order – pcv_importance() treats the PCV as a value function over variable subsets and attributes it fairly: the flagship "shapley" method averages each variable's marginal PCV over every possible entry order, and "dominance" adds Budescu's pairwise dominance detail (its general dominance weights coincide with the Shapley values). Both methods satisfy the efficiency identity: the contributions sum exactly to the full model's total PCV.

Two attribution targets are useful in a MAIHDA. Passing the stratum dimensions (e.g. c("gender", "race", "education")) splits the additive share – the canonical null-to-adjusted PCV – fairly across the dimensions: which social dimension drives the additive between-stratum inequality. Passing individual-level covariates gives an order-invariant counterpart to stepwise_pcv(): how much each covariate explains of the between-stratum variance (subject to the latent-scale caveat below for non-Gaussian families).

Usage

pcv_importance(
  data,
  outcome,
  vars,
  method = c("shapley", "dominance"),
  approx = c("exact", "montecarlo"),
  n_perm = 2000,
  engine = "lme4",
  family = "gaussian",
  context = NULL,
  sampling_weights = NULL,
  bootstrap = FALSE,
  n_boot = 1000,
  conf_level = 0.95,
  estimation = c("fitted", "ML")
)

Arguments

Details

Value function and efficiency. Write V_0 for the null model's between-stratum variance and V(S) for the between-stratum variance after adding the variable subset S as fixed effects (each put on the estimation basis – the default "fitted" keeps each lmer fit's own REML variance, "ML" refits it with maximum likelihood – exactly as in calculate_pcv and stepwise_pcv, so all attributions live on the same scale as the rest of the package). The value function is v(S) = (V_0 - V(S)) / V_0 – the total PCV of the model that adds S – and the Shapley contribution of variable i averages its marginal PCV v(S \cup \{i\}) - v(S) over all subsets with the usual Shapley weights. Because every permutation's marginals telescope to v(N), the contributions of every method here (including the Monte-Carlo approximation, for any n_perm) sum exactly to the full-model total PCV. This is the multilevel-PCV analogue of the LMG / Shorrocks-Shapley decomposition of R^2 (Groemping 2006; Shorrocks 2013).

Exact vs. Monte-Carlo cost. The exact attribution fits 2^k - 1 non-empty subset models plus the null – 2^k fits in total, with every fit cached and reused across all marginal differences: 256 fits at k = 8, 1024 at k = 10. That is feasible for lme4 but quickly infeasible beyond, and each of those fits is a separate Stan run under engine = "brms" – exact attribution on brms is therefore strongly discouraged except for very small k. The Monte-Carlo route samples n_perm entry orders (unbiased for the Shapley values, since a uniform random permutation reproduces the Shapley weights), reports a per-variable Monte-Carlo standard error, and warns when the largest MC_SE exceeds 0.01 on the PCV scale (increase n_perm).

Latent-scale families and rescaling. For binomial/ordinal (and, in attenuated form, count) families, adding a predictor that varies within strata rescales the latent metric, so its marginal PCV mixes explained variance with rescaling (see the latent-scale note in calculate_pcv); Shapley values of the PCV inherit this. Attribution among the stratum dimensions – constant within each stratum – is largely unaffected.

Suppression and negative contributions. The PCV can be negative, so a contribution (and its Share) can be negative or exceed 100%. Efficiency still holds; contributions are reported signed and no non-negative normalisation is applied. A negative Shapley contribution flags a suppressor-style variable whose inclusion tends to raise the between-stratum variance.

Bootstrap. With bootstrap = TRUE the whole attribution is bootstrapped: responses are simulated from the full model and every subset model is refit per draw (lme4::refit), giving percentile intervals per contribution – n_boot times n_fits refits, so gate it by cost. Available for engine = "lme4" with exact attribution (any method); the Monte-Carlo approximation already carries its own sampling error and is not combined with the bootstrap. As in calculate_pcv, refit() holds a glmer.nb dispersion parameter fixed at its original estimate.

Value

An object of class maihda_pcv_importance: a list with

Reproducibility

approx = "montecarlo" draws random permutations with the session RNG; call set.seed() first for reproducible results. The exact methods give the same contributions on every run of a likelihood engine, but the bootstrap = TRUE intervals come from simulated responses and an engine = "brms" fit is an MCMC draw, so call set.seed() first for those as well.

References

Budescu, D. V. (1993). Dominance analysis: a new approach to the problem of relative importance of predictors in multiple regression. Psychological Bulletin, 114(3), 542-551.

Groemping, U. (2006). Relative importance for linear regression in R: the package relaimpo. Journal of Statistical Software, 17(1), 1-27.

Shorrocks, A. F. (2013). Decomposition procedures for distributional analysis: a unified framework based on the Shapley value. Journal of Economic Inequality, 11, 99-126.

See Also

stepwise_pcv for the sequential path with the discriminatory-accuracy trajectory, calculate_pcv for the two-model PCV and the latent-scale caveat, and maihda for the canonical null-vs-adjusted decomposition.

Examples

strata <- make_strata(maihda_sim_data, c("gender", "race"))

# Order-invariant split of the PCV across the two dimensions and age
imp <- pcv_importance(strata$data, "health_outcome", c("gender", "race", "age"))
print(imp)
plot(imp)

# Dominance analysis: adds conditional / complete dominance detail
pcv_importance(strata$data, "health_outcome", c("gender", "race", "age"),
               method = "dominance")

# Monte-Carlo approximation (set a seed for reproducibility)
set.seed(42)
pcv_importance(strata$data, "health_outcome", c("gender", "race", "age"),
               approx = "montecarlo", n_perm = 500)

Plot a MAIHDA Analysis

Description

Dispatches each type to the model it is valid on. The VPC and shrinkage views ("vpc", "obs_vs_shrunken", "predicted") use the null model by default; for type = "vpc" the model argument selects the null model, the adjusted model, or "both" – one combined plot that shows the null-to-adjusted change with the PCV annotated. The additive-vs-intersectional views ("effect_decomp", "prediction_deviation") use the adjusted model, whose fixed effects carry the dimensions' additive part so the stratum random effect is the pure interaction; with fewer than two dimensions (no adjusted model) they fall back to the null model. Group types ("group_vpc", "group_components", "group_between_variance", "group_pcv", "group_additive_share") use the group comparison when maihda was called with a group.

Usage

## S3 method for class 'maihda_analysis'
plot(
  x,
  type = "all",
  highlight_interactions = FALSE,
  only_flagged = FALSE,
  highlight_by = c("flag", "rope"),
  rope = NULL,
  select = c("order", "deviation"),
  order_by = c("predicted_desc", "stratum", "predicted_asc", "deviation", "size"),
  model = c("null", "adjusted", "both"),
  ...
)

Arguments

Value

A ggplot2 object – including type = "vpc", model = "both", which is a single combined change plot – or (for type = "all") an invisible list of them.

Plot a MAIHDA Model Comparison

Description

Plots VPC/ICC across the models compared by compare_maihda. Dispatched via plot() on the classed result.

Usage

## S3 method for class 'maihda_comparison'
plot(x, ...)

Arguments

Value

A ggplot2 object.

Examples

strata <- make_strata(maihda_sim_data, vars = c("gender", "race"))

null_model <- fit_maihda(health_outcome ~ 1 + (1 | stratum), data = strata$data)
adj_model  <- fit_maihda(health_outcome ~ age + (1 | stratum), data = strata$data)

comparison <- compare_maihda(null_model, adj_model, bootstrap = TRUE)
plot(comparison)

Plot a MAIHDA sample description

Description

Descriptive plots for a maihda_describe object: the stratum-size distribution (the sparsity of the strata space), the family-aware outcome distribution, and per-variable missingness.

Usage

## S3 method for class 'maihda_describe'
plot(x, type = c("stratum_size", "outcome", "missingness"), ...)

Arguments

Value

A ggplot2 object, extendable with the usual + grammar.

Examples

desc <- maihda_describe(BMI ~ Age + (1 | Gender:Race:Education),
                        data = maihda_health_data)
plot(desc, type = "stratum_size")
plot(desc, type = "outcome")
plot(desc, type = "missingness")

Plot a MAIHDA Group Comparison

Description

Visualises the output of compare_maihda_groups as a point/forest plot of the VPC/ICC by group, as stacked variance-composition bars (between- vs within-stratum share) by group, as bars of the absolute between-stratum (intersectional) variance by group, or as bars of the additive share (PCV) by group. Dispatched via plot() on the classed result.

Usage

## S3 method for class 'maihda_group_comparison'
plot(
  x,
  type = c("vpc", "components", "between_variance", "pcv", "additive_share"),
  ...
)

Arguments

Value

A ggplot2 object.

Examples

data(maihda_health_data)
cmp <- compare_maihda_groups(BMI ~ Age + (1 | Gender:Race),
                             data = maihda_health_data, group = "Education")
plot(cmp, type = "vpc")
plot(cmp, type = "components")
plot(cmp, type = "between_variance")

Plot MAIHDA Model Results

Description

Creates various plots for visualizing MAIHDA model results including variance partition coefficient comparisons, observed vs. shrunken estimates, and predicted subgroup values with conditional intervals.

Usage

## S3 method for class 'maihda_model'
plot(
  x,
  type = c("all", "vpc", "obs_vs_shrunken", "predicted", "upset", "effect_decomp",
    "prediction_deviation", "context_vpc", "vpc_trajectory", "trajectories"),
  summary_obj = NULL,
  n_strata = 50,
  highlight_interactions = FALSE,
  only_flagged = FALSE,
  highlight_by = c("flag", "rope"),
  rope = NULL,
  select = c("order", "deviation"),
  order_by = c("predicted_desc", "stratum", "predicted_asc", "deviation", "size"),
  quantity = c("predicted", "interaction"),
  ...
)

Arguments

Value

For a single type, a ggplot2 object that you can extend with the usual + grammar (themes, labs(), added layers, or a replacement fill/colour scale). Some types return a richer object: "prediction_deviation" returns a patchwork of two panels and "upset" a patchwork of three panels (theme every panel at once with & theme_*()). type = "all" returns a named list of ggplot2 objects.

Examples

strata_result <- make_strata(maihda_sim_data, vars = c("gender", "race"))
model <- fit_maihda(health_outcome ~ age + (1 | stratum), data = strata_result$data)

# VPC plot
plot(model, type = "vpc")

# Single-type plots are ggplot objects -- restyle them with ggplot2:
plot(model, type = "vpc") +
  ggplot2::theme_classic() +
  ggplot2::labs(title = "Variance partition, restyled")

# Generate all plots (a named list); pick one out to restyle it:
plots <- plot(model)
plots$predicted + ggplot2::theme_bw()

Plot a PCV importance attribution

Description

Bar chart of the per-variable contributions to the total PCV, ordered by size, with a dashed zero line (negative bars flag suppression) and, when available, parametric-bootstrap confidence intervals. The bars sum to the full-model Total PCV (the efficiency property).

Usage

## S3 method for class 'maihda_pcv_importance'
plot(x, ...)

Arguments

Value

A ggplot object.

Stratum vs. Context Variance Plot (contextual cross-classified MAIHDA)

Description

One bar per variance component – the between-stratum (intersectional) variance, each context's variance, any other random effects, and the residual – on the variance scale, with each component's share of the total printed above its bar. Complements plot_vpc()'s stacked proportion bar by showing the magnitudes the shares are computed from. For a crossed-dimensions fit with a context the between-stratum variance is drawn as its additive dimension and interaction components.

Usage

plot_context_vpc(summary_obj)

Arguments

Value

A ggplot2 object.

Effect Decomposition Plot

Description

Decomposes the total deviation from the overall mean into the additive (fixed) component and the intersectional (random) component for each stratum.

Usage

plot_effect_decomposition(
  object,
  summary_obj,
  top_n_labels = 10,
  highlight = NULL
)

Arguments

Details

Only the intersectional random effects enter the deviations and the overall mean: the stratum effect and, for a crossed-dimensions fit, the dimension effects. A contextual random effect, or any other grouping, is part of neither component and is excluded.

Value

A ggplot2 object

Observed vs. Shrunken Estimates Plot

Description

Observed vs. Shrunken Estimates Plot

Usage

plot_obs_vs_shrunken(
  object,
  summary_obj,
  highlight = NULL,
  only_flagged = FALSE
)

Arguments

Details

The x-axis is each stratum's raw observed mean; the y-axis is the model-based stratum estimate, which includes the fixed-effect contribution. Both are on the scale the response was fitted on, which is what makes the y = x diagonal meaningful: for a transformed response such as log(y) ~ x the observed values are stratum means of log(y), not of y. To read the panel on the original scale, fit the model on that scale; back-transforming only one axis would not give comparable quantities. For an intercept-only (null) model the vertical distance from the diagonal is pure shrinkage toward the grand mean. For a covariate-adjusted model the model estimate also moves with the stratum's covariate profile, so distance from the diagonal reflects both shrinkage and covariate adjustment and should not be read as shrinkage alone.

Value

A ggplot2 object

Time-specific PCV plot (longitudinal MAIHDA)

Description

The additive share – the proportional change in the between-stratum (trajectory) variance from the null to the adjusted model – as a function of time. A high, flat curve means intersectional trajectory inequalities are "mostly additive".

Usage

plot_pcv_trajectory(pcv)

Arguments

Value

A ggplot2 object.

Plot Predicted Stratum Values with Confidence Intervals

Description

Plot Predicted Stratum Values with Confidence Intervals

Usage

plot_predicted_strata(
  object,
  summary_obj,
  n_strata,
  scale = c("response", "link"),
  highlight = NULL,
  only_flagged = FALSE,
  select = c("order", "deviation"),
  order_by = c("predicted_desc", "stratum", "predicted_asc", "deviation", "size")
)

Arguments

Value

A ggplot2 object

Plot Prediction Deviation Panels

Description

Creates an advanced, publication-ready two-panel dashboard for visualizing predicted values and highlighting the most notable cases or strata. What "notable" means depends on the model type, and the labelled points are not statistical outliers in the regression-diagnostic sense:

• Gaussian and Poisson (and the ordinal "expected_score" mode): the cases/strata whose prediction sits furthest from the mean prediction (largest deviation), ranked by absolute deviation.

• Binomial: the cases/strata with the largest absolute deviance residual, i.e. where the observed 0/1 outcome – for an aggregated binomial (cbind(successes, failures), a proportion with its trial counts as prior weights, or a brms y | trials(n) fit), the observed successes out of trials – is least consistent with the fitted probability (worst-fit points), ranked by |deviance residual|. As in R's binomial family, a prior weight counts as that many trials, so a 0/1 outcome with weights is scored as that many successes or failures. Every row's outcome is coded as the model coded it – by label, whatever the order of a factor's levels in data – and a factor with more than two levels, as in R's binomial family, as its first level against all the others. Predictions are per-trial probabilities for every binomial fit. A row whose outcome or trial count is unknown has no residual: it is left out of the stratum residual means, is never labelled, and is drawn as a hollow point.

• Ordinal "surprise" mode: the cases/strata with the highest surprise -\log P(\text{observed category}), i.e. the least probable observations under the model. It needs the observed response in data: a row whose response is missing is left out, and one whose category is not among the model's fitted categories is left out with a warning.

Usage

plot_prediction_deviation_panels(
  model,
  data = NULL,
  type = c("auto", "gaussian", "poisson", "binomial", "ordinal"),
  ordinal_mode = c("surprise", "expected_score"),
  top_n_labels = 5,
  strata_info = NULL
)

Arguments

Details

The vertical bars are 95% intervals for the plotted values. On a stratum panel each bar is an interval for the stratum's weighted mean prediction, the point drawn: posterior quantiles of that mean for a brms fit, and for an lme4 fit a delta-method interval from the joint covariance of the fixed and random effects (Henderson's mixed-model equations), formed on the logit scale for a probability and the log scale for a count, so that a bar stays inside the range of the response. A case-level panel of an lme4 fit draws the same interval for each row. The lme4 interval treats the variance parameters as known: it runs narrow when strata are few or small or an outcome is rare, and a random effect whose variance is estimated at zero adds no uncertainty to it. Other models draw case-level bars from the standard errors of their predict() method and no stratum bars; WeMix and ordinal fits draw none.

Value

A patchwork object containing two ggplot2 panels.

References

Henderson, C. R. (1975). Best linear unbiased estimation and prediction under a selection model. Biometrics, 31(2), 423-447. doi:10.2307/2529430

Bates, D., Maechler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01

Kackar, R. N., & Harville, D. A. (1984). Approximations for standard errors of estimators of fixed and random effects in mixed linear models. Journal of the American Statistical Association, 79(388), 853-862. doi:10.1080/01621459.1984.10477102

Stratum mean-trajectory plot (longitudinal MAIHDA)

Description

One predicted line per stratum over time – that stratum's fixed-part trajectory plus its own random intercept and slope (BLUPs) – the longitudinal analogue of the predicted-strata caterpillar. Shows how the intersectional groups fan out (or converge) over time.

Usage

plot_stratum_trajectories(
  object,
  summary_obj,
  n_strata = 50,
  select = c("order", "deviation")
)

Arguments

Details

The fixed part varies with the stratum drawn: the stratum-defining dimensions take that stratum's values, so an adjusted growth model's dimension main effects and dim:time interactions are included, while every other covariate is held at a shared reference profile (mean for a numeric, modal level for a factor). The lines are therefore predicted stratum trajectories net of covariate composition, excluding the individual-level random effects. For a null growth model the fixed part carries no dimension terms and every line shares the one population trajectory.

Value

A ggplot2 object.

UpSet-style Predicted Stratum Plot

Description

Composite alternative to the text-labelled "predicted" view that replaces the long intersectional axis labels with an UpSet-style category matrix. Three panels share one column order: a top bar of intersection (stratum) sizes, a middle matrix encoding each stratum's category on every dimension, and a bottom panel of predicted values with conditional intervals. Columns are ordered by intersection size (largest first) by default, or by any other order_by rule – "predicted_desc" turns the view into the ranked caterpillar of plot_predicted_strata() drawn against the category matrix instead of long text labels. Binary 0/1 (or logical) dimensions collapse to a single present/absent row; multi-level factors get one row per level, and each column lights exactly one dot per dimension.

Usage

plot_upset_strata(
  object,
  summary_obj,
  n_strata,
  scale = c("response", "link"),
  highlight = NULL,
  only_flagged = FALSE,
  select = c("order", "deviation"),
  order_by = c("size", "predicted_desc", "stratum", "predicted_asc", "deviation"),
  quantity = c("predicted", "interaction")
)

Arguments

Value

A patchwork object stacking the three panels.

VPC Visualization Plot

Description

VPC Visualization Plot

Usage

plot_vpc(summary_obj)

Arguments

Value

A ggplot2 object

Null-vs-Adjusted Variance-Partition Change Plot

Description

Draws the null and adjusted models' variance partitions as two stacked bars on a shared axis and annotates the PCV, so the drop in the between-stratum share after accounting for the dimensions' additive main effects is visible in one figure. Backs plot(<maihda_analysis>, type = "vpc", model = "both") for a cross-sectional analysis.

Usage

plot_vpc_change(null_summary, adjusted_summary, pcv = NULL)

Arguments

Value

A ggplot2 object.

Time-varying VPC trajectory plot (longitudinal MAIHDA)

Description

The between-stratum share of variance (VPC/ICC) as a function of time, with a confidence/credible ribbon when available. The headline reference-time VPC is marked. For a longitudinal MAIHDA the VPC is not a single number – the between-stratum variance is a random intercept + slope on time – so this curve replaces the cross-sectional VPC bar.

Usage

plot_vpc_trajectory(summary_obj)

Arguments

Value

A ggplot2 object.

Null-vs-Adjusted VPC-over-time Change Plot (longitudinal MAIHDA)

Description

Overlays the null and adjusted models' time-varying VPC/ICC curves on one axis, so the reduction in the between-stratum share after adjustment is visible over the whole time range. Backs plot(<maihda_analysis>, type = "vpc", model = "both") for a longitudinal analysis; complements plot_pcv_trajectory (the additive-share PCV(t)).

Usage

plot_vpc_trajectory_change(null_summary, adjusted_summary)

Arguments

Value

A ggplot2 object.

Predict method for maihda_model objects

Description

S3 predict method for MAIHDA models: a thin alias of predict_maihda, so predict(model, type = "strata") and predict_maihda(model, type = "strata") are interchangeable. See predict_maihda for the full documentation of the arguments and return value.

Usage

## S3 method for class 'maihda_model'
predict(
  object,
  newdata = NULL,
  type = c("individual", "strata", "response", "link"),
  scale = c("response", "link"),
  allow_new_levels = FALSE,
  ...
)

Arguments

Value

Depending on type:

• For "individual": A numeric vector of predicted values on the requested scale, one per row of newdata – or, when newdata is NULL, one per ANALYTIC row, matching nrow(object$data). A model fitted with na.action = na.exclude is not padded back out to the original input rows the way predict() on the underlying fit is, so the result can always be bound to object$data.

• For "strata": A data frame with stratum ID and predicted random effect. When newdata is supplied, the result is restricted to the strata present in newdata (and a stratum the model never saw is an error, as for "individual"); when newdata is NULL, every training stratum is returned.

See Also

predict_maihda

Predict from MAIHDA Model

Description

Makes predictions from a fitted MAIHDA model, either at the stratum level or individual level.

Usage

predict_maihda(
  object,
  newdata = NULL,
  type = c("individual", "strata", "response", "link"),
  scale = c("response", "link"),
  allow_new_levels = FALSE,
  ...
)

Arguments

Value

Depending on type:

• For "individual": A numeric vector of predicted values on the requested scale, one per row of newdata – or, when newdata is NULL, one per ANALYTIC row, matching nrow(object$data). A model fitted with na.action = na.exclude is not padded back out to the original input rows the way predict() on the underlying fit is, so the result can always be bound to object$data.

• For "strata": A data frame with stratum ID and predicted random effect. When newdata is supplied, the result is restricted to the strata present in newdata (and a stratum the model never saw is an error, as for "individual"); when newdata is NULL, every training stratum is returned.

Examples

strata_result <- make_strata(maihda_sim_data, vars = c("gender", "race"))
model <- fit_maihda(health_outcome ~ age + (1 | stratum), data = strata_result$data)

# Individual predictions
pred_ind <- predict_maihda(model, type = "individual")

# Stratum predictions
pred_strata <- predict_maihda(model, type = "strata")

Print a MAIHDA Analysis

Description

Print a MAIHDA Analysis

Usage

## S3 method for class 'maihda_analysis'
print(x, ...)

Arguments

Value

No return value, called for side effects.

Print a MAIHDA sample description

Description

Print a MAIHDA sample description

Usage

## S3 method for class 'maihda_describe'
print(x, digits = x$digits, ...)

Arguments

Value

Invisibly, x.

Print method for MAIHDA group comparisons

Description

Print method for MAIHDA group comparisons

Usage

## S3 method for class 'maihda_group_comparison'
print(x, ...)

Arguments

Value

No return value, called for side effects.

Print MAIHDA information criteria

Description

Print MAIHDA information criteria

Usage

## S3 method for class 'maihda_ic'
print(x, ...)

Arguments

Value

No return value, called for side effects.

Print a MAIHDA interaction diagnostic

Description

Print a MAIHDA interaction diagnostic

Usage

## S3 method for class 'maihda_interactions'
print(x, ...)

Arguments

Value

No return value, called for side effects.

Print a longitudinal MAIHDA PCV

Description

Print a longitudinal MAIHDA PCV

Usage

## S3 method for class 'maihda_long_pcv'
print(x, ...)

Arguments

Value

The object, invisibly.

Print method for maihda_model

Description

Print method for maihda_model

Usage

## S3 method for class 'maihda_model'
print(x, ...)

Arguments

Value

No return value, called for side effects.

Print a PCV importance attribution

Description

Print a PCV importance attribution

Usage

## S3 method for class 'maihda_pcv_importance'
print(x, digits = 4, ...)

Arguments

Value

Invisibly, x.

Print a stepwise MAIHDA table

Description

Print a stepwise MAIHDA table

Usage

## S3 method for class 'maihda_stepwise'
print(x, ...)

Arguments

Value

Invisibly, x.

Print method for maihda_strata objects

Description

Print method for maihda_strata objects

Usage

## S3 method for class 'maihda_strata'
print(x, ...)

Arguments

Value

No return value, called for side effects.

Print method for maihda_summary objects

Description

Print method for maihda_summary objects

Usage

## S3 method for class 'maihda_summary'
print(x, ...)

Arguments

Value

No return value, called for side effects.

Print a MAIHDA results table

Description

Print a MAIHDA results table

Usage

## S3 method for class 'maihda_table'
print(x, digits = x$digits, ...)

Arguments

Value

Invisibly, x.

Print method for PCV results

Description

Print method for PCV results

Usage

## S3 method for class 'pcv_result'
print(x, ...)

## S3 method for class 'pvc_result'
print(x, ...)