Skip to contents

Longitudinal MAIHDA: intersectional inequalities over time

Hamid Bulut

2026-09-28

Source: vignettes/longitudinal.Rmd
longitudinal.Rmd

From a snapshot to a trajectory

Standard (cross-sectional) MAIHDA answers “do intersectional strata differ?” at a single point in time, summarised by the between-stratum VPC. With repeated measurements we can ask a richer question: “do strata differ in how they change over time?” – and decompose those trajectory differences into an additive part (the dimensions’ main effects and their interactions with time) and a multiplicative part (an intersectional trajectory beyond additive).

This is the longitudinal MAIHDA of Bell, Evans, Holman & Leckie (2024). It is a 3-level growth-curve model with measurement occasions (level 1) within individuals (level 2) within intersectional strata (level 3), with a random intercept and a slope on time at both the individual and stratum levels:

ytij=β0+β1t+(u0j+u1jt)⏟stratum+(v0ij+v1ijt)⏟individual+etij. y_{tij} = \beta_0 + \beta_1 t + \underbrace{(u_{0j} + u_{1j}\,t)}_{\text{stratum}} + \underbrace{(v_{0ij} + v_{1ij}\,t)}_{\text{individual}} + e_{tij}.

Because the stratum random effect now has a slope, the between-stratum variance becomes a function of time:

Var⁡S(t)=σu02+2tσu01+t2σu12,VPC⁡S(t)=Var⁡S(t)Var⁡S(t)+Var⁡I(t)+σe2. \operatorname{Var}_S(t) = \sigma^2_{u0} + 2t\,\sigma_{u01} + t^2 \sigma^2_{u1}, \qquad \operatorname{VPC}_S(t) = \frac{\operatorname{Var}_S(t)} {\operatorname{Var}_S(t) + \operatorname{Var}_I(t) + \sigma^2_e}.

The data

The bundled maihda_long_data is a simulated panel: 600 people, each measured over five waves, within 12 strata of gender ×\times ethnicity ×\times education. The trajectory differences are constructed to be mostly additive with one genuine interaction, so the decomposition below is interpretable.

data(maihda_long_data)
head(maihda_long_data)
#>      id wave gender ethnicity education age wellbeing low_wellbeing
#> 1 P0001    0    Men      EthB       Low  50     3.732             1
#> 2 P0001    1    Men      EthB       Low  50     3.696             1
#> 3 P0001    2    Men      EthB       Low  50     1.638             1
#> 4 P0001    3    Men      EthB       Low  50     3.585             1
#> 5 P0001    4    Men      EthB       Low  50     3.594             1
#> 6 P0002    0  Women      EthB       Low  56     5.631             0

It is long format, one row per person-occasion, with a person id (id) and a numeric time (wave).

Fitting and the time-varying VPC

Supply id and time to fit_maihda(). You only need to write the strata shorthand (1 | var1:var2); the growth slopes on time are added for you.

m <- fit_maihda(wellbeing ~ wave + (1 | gender:ethnicity:education),
                data = maihda_long_data, id = "id", time = "wave")
summary(m)
#> MAIHDA Model Summary
#> ====================
#> 
#> Variance Partition Coefficient (VPC/ICC) at baseline (wave = 0):
#>   Estimate: 0.3985
#> 
#> Variance Components:
#>                                            component variance     sd
#>                Between-stratum: intercept (time = 0)  0.40774 0.6385
#>                        Between-stratum: slope (wave)  0.04284 0.2070
#>          Between-stratum: intercept-slope covariance  0.08848     NA
#>        Between-individual (id): intercept (time = 0)  0.25558 0.5055
#>                Between-individual (id): slope (wave)  0.01936 0.1391
#>  Between-individual (id): intercept-slope covariance -0.00110     NA
#>                                    Within (residual)  0.35982 0.5999
#> 
#> Time-varying VPC/ICC (between-stratum share over wave):
#>   range 0.3985 to 0.6628 across wave in [0, 4].
#>   The between-stratum variance is a function of time (random intercept +
#>   slope), so the VPC varies; it depends on where time is zeroed. See
#>   plot(type = "vpc_trajectory") for the full curve.
#> 
#> Trajectory VPCs (Bell et al. 2024, eq. 5; occasion-level variance excluded):
#>   Intercept (wave = 0): 0.6147    Slope: 0.6888
#>   These ask how intersectionally patterned trajectories are; the VPC
#>   above asks how much of an observed measurement a stratum explains.
#>   Intercept exceeds it at the reference time; Slope is on a different
#>   scale and is not comparable to it. See ?summary.maihda_model.
#> 
#> Fixed Effects (Wald t, 95% intervals):
#>         term  estimate     se statistic df  p_value   lower  upper
#>  (Intercept)  4.986465 0.1870   26.6719 10 1.27e-10  4.5699 5.4030
#>         wave -0.006306 0.0607   -0.1039 10    0.919 -0.1415 0.1289
#>   df: containment (between-within).
#> 
#> Stratum baseline (intercept) deviations (first 10):
#>  stratum stratum_id               label random_effect      se lower_95 upper_95
#>        1          1    Men × EthB × Low      -0.63793 0.08953 -0.81340  -0.4625
#>        2          2  Women × EthB × Low      -0.06461 0.08594 -0.23306   0.1038
#>        3          3   Men × EthA × High       0.86073 0.07782  0.70820   1.0133
#>        4          4   Men × EthB × High       0.09123 0.10259 -0.10985   0.2923
#>        5          5 Women × EthA × High       0.93600 0.08804  0.76344   1.1086
#>        6          6 Women × EthC × High       0.35847 0.13544  0.09300   0.6239
#>        7          7   Men × EthC × High      -0.11459 0.13286 -0.37500   0.1458
#>        8          8    Men × EthA × Low      -0.32848 0.07050 -0.46665  -0.1903
#>        9          9  Women × EthA × Low       0.15991 0.07326  0.01633   0.3035
#>       10         10    Men × EthC × Low      -0.82706 0.13286 -1.08747  -0.5667
#>   ... and 2 more strata
#>   (random slope not shown; use predict(type = "strata") for the per-stratum intercept and slope, or plot(type = "trajectories")).

summary() reports the VPC at the baseline (reference time, the earliest wave) and the full trajectory of the VPC across the observed times. Plot it:

plot(m, type = "vpc_trajectory")   # VPC(t), with the reference time marked

plot(m, type = "trajectories")     # predicted per-stratum mean trajectories

A rising VPC trajectory means the strata fan out over time (intersectional inequality grows).

Two VPCs, and which one to report

summary() reports two different variance partitions for a growth model, and they are not interchangeable. They differ in a single denominator term: the level-1 (occasion) residual σe2\sigma^2_e, which here is within-individual volatility – how far one measurement falls from that person’s own smooth trajectory. It mixes measurement error, real short-term fluctuation, and any misfit of the functional form. A cross-sectional MAIHDA cannot separate it from between-individual variance at all; repeated measures are what split the two.

The headline VPC keeps the term, and therefore answers the discriminatory-accuracy question of how much of an observed measurement a stratum accounts for:

VPC⁡S(t)=Var⁡S(t)Var⁡S(t)+Var⁡I(t)+σe2.\operatorname{VPC}_S(t) = \frac{\operatorname{Var}_S(t)} {\operatorname{Var}_S(t) + \operatorname{Var}_I(t) + \sigma^2_e}.

The trajectory VPCs drop this term, following Bell et al. (2024), equation (5):

VPC⁡intercept=Var⁡S(t0)Var⁡S(t0)+Var⁡I(t0),VPC⁡slope=SlopeVar⁡S(t0)SlopeVar⁡S(t0)+SlopeVar⁡I(t0).\operatorname{VPC}_{\text{intercept}} = \frac{\operatorname{Var}_S(t_0)} {\operatorname{Var}_S(t_0) + \operatorname{Var}_I(t_0)}, \qquad \operatorname{VPC}_{\text{slope}} = \frac{\operatorname{SlopeVar}_S(t_0)} {\operatorname{SlopeVar}_S(t_0) + \operatorname{SlopeVar}_I(t_0)}.

s <- summary(m)
c(intercept = s$longitudinal$vpc_intercept, slope = s$longitudinal$vpc_slope)
#> intercept     slope 
#> 0.6146997 0.6887613

These ask how intersectionally patterned people’s trajectories are – what share of the between-individual variation in where a trajectory starts, and in how fast it changes, lies between strata. Because σe2\sigma^2_e is absent, neither shrinks when the outcome is measured noisily, which makes them comparable across studies using different instruments where the headline VPC is not.

Both are shares of a between-stratum variance, which a MAIHDA estimates from few strata by construction, so report them with their interval – vpc_intercept_ci and vpc_slope_ci, filled in by summary(bootstrap = TRUE) here (not run in the vignette, because it refits the model many times) and by the posterior for a brms fit, which needs no refitting:

s_ci <- summary(m, bootstrap = TRUE, n_boot = 500)
s_ci$longitudinal$vpc_intercept_ci
s_ci$longitudinal$vpc_slope_ci

The point estimates differ in kind between the two engines, and the difference is not cosmetic. For lme4 they are the plug-in from the fitted covariance blocks – the two numbers printed above. For brms they are the posterior median of the per-draw ratio, not the ratio of the posterior-mean variance components: the between-stratum variance has a right-skewed posterior whenever the strata are few, so a plug-in built from its posterior mean runs high. Refitting these same twelve strata with brms, over 150 of the individuals, the plug-in gives 0.6145 where the posterior median is 0.5765 – a 6.6% overstatement, and one that widens as the strata get fewer. The interval matters more than the shift: on that fit it runs from 0.34 to 0.82, half the unit interval, and it was previously not reported at all.

Report the headline VPC unless you specifically mean the trajectory question. It is the quantity comparable to published cross-sectional MAIHDA VPCs, and it is the one that speaks to discriminatory accuracy. The trajectory VPCs should be reported next to PCV_intercept and PCV_slope, and vpc_slope is the only scale-free way to describe how intersectional the rates of change are.

For the slope the residual could not have been included even in principle: a slope variance is in (outcome/time)2(\text{outcome}/\text{time})^2 while σe2\sigma^2_e is in outcome2\text{outcome}^2, so the sum would be dimensionally meaningless. (The headline VPC is well formed because a(t)′Σa(t)a(t)'\Sigma a(t) returns to outcome2\text{outcome}^2.) The intercept VPC drops it for symmetry with the slope.

Both are evaluated at the baseline t0t_0; the intercept VPC depends on where time is zeroed and the slope VPC does not, as Bell et al. note. Their own examples centre on mean age rather than the baseline, so an intercept VPC replicated from the paper will differ from the one reported here unless the reference points are aligned.

Decomposing the trajectory: additive vs. multiplicative

maihda(decomposition = "longitudinal") (selected automatically when id/time are supplied) fits a null growth model and an adjusted growth model. The adjusted model adds the dimensions’ main effects and their interactions with time (dim:time), so the stratum-level variance it leaves behind is the interaction beyond additive.

a <- maihda(wellbeing ~ wave + (1 | gender:ethnicity:education),
            data = maihda_long_data, id = "id", time = "wave",
            decomposition = "longitudinal")
a$pcv
#> Longitudinal PCV (additive vs. multiplicative intersectionality)
#> ================================================================
#> 
#> Baseline (wave = 0) variance: 0.4077 (null) -> 0.0332 (adjusted)
#>   PCV_intercept: 91.9% of the baseline between-stratum inequality is additive.
#> Slope (wave) variance at baseline: 0.0428 (null) -> 0.0058 (adjusted)
#>   PCV_slope:     86.6% of the *trajectory* between-stratum inequality is additive
#>                  (the remainder is the multiplicative/interaction part).
#> 
#> The PCV is the share of the null model's between-stratum (trajectory) variance
#> explained by the dimensions' additive main effects and their time interactions;
#> a high PCV_slope means trajectory inequalities are 'mostly additive'.
#>

The PCV is reported separately for the two pieces of the trajectory:

  • PCV_intercept describes the share of the baseline between-stratum inequality explained by the additive main effects.
  • PCV_slope describes the share of the trajectory (slope) between-stratum inequality explained additively. A high PCV_slope means the “trajectories are mostly additive” finding; the remainder is the multiplicative/interaction part.
plot(a, type = "pcv_trajectory")   # the additive share over time

Scope and cautions

  • Time coding. You can code time however is convenient for your workflow ( waves 0, 1, 2, …, waves 10, 11, …, age, or calendar year). When the axis does not start at 0, the growth terms are fit on internally centered time (time - min(time), with a message), because the raw polynomial basis over an offset range is ill-conditioned and can silently converge to a wrong solution. All results ( the VPC trajectory, the PCV, plots, and predictions ) are reported on your original time scale, and the baseline (ref_time) is always the earliest observed time.

  • Identifiability. A stratum random slope needs enough occasions per stratum; sparse strata, few waves, or irregular measurement times can give a singular lme4 fit (the brms engine handles this better). The fit diagnostics name the block that is at the boundary. Which block it is decides whether the result is affected. For example a boundary (time | id) block leaves the between-stratum variance well estimated, a boundary (time | stratum) block does not.

  • Rank-deficient stratum blocks. The usual boundary outcome is a stratum intercept-slope correlation of exactly 1 or -1 (Corr = 1.00 in VarCorr()). The block is then rank 1 (one row), one parameter short of a full covariance matrix, and the stratum slope variance is not separately identified. PCV_slope is NA in that case, and print() on the PCV object says why. PCV_intercept and the VPC trajectory come from the intercept variance, which is still estimable.

    How often this happens depends on the stratum-level slope variance. it is common when the slope variance is near zero, but rare when it is near zero. A rank-deficient block is usually evidence that the strata do not diverge over time.

    Both structures share the same fixed effects, so their fits can be compared with AIC:

    m_slope <- fit_maihda(wellbeing ~ wave + (1 | gender:ethnicity:education),
                      data = maihda_long_data, id = "id", time = "wave")
m_level <- fit_maihda(wellbeing ~ wave + (1 | gender:ethnicity:education),
                      data = maihda_long_data, id = "id", time = "wave",
                      stratum_slope = FALSE)
AIC(m_slope$model) - AIC(m_level$model)   # negative favours keeping the slope

    stratum_slope = FALSE fits (time | id) + (1 | stratum). Individuals keep their growth curves; strata are constrained to differ in level only. Report the constraint, when the absence of trajectory divergence is assumed, and therefore not estimated. The between-stratum variance is then one time-constant number and the decomposition reports PCV_intercept (no slope) only. The VPC still varies with time, because the individual-level slope variance and the residual are in its denominator.

  • Not a single number. The VPC is time-varying, so extract_between_variance() and calculate_pcv() deliberately refuse a longitudinal model; use the longitudinal decomposition above.

  • Which VPC. summary() reports both the discriminatory-accuracy VPC(t) and the Bell et al. trajectory VPCs, which use a different denominator; see Two VPCs, and which one to report.

  • Out of scope (v1). Design-weighted, contextual (stratum ×\times place ×\times time), and wemix/ordinal longitudinal models are not yet supported.

Reference

Bell, A., Evans, C., Holman, D., & Leckie, G. (2024). Extending intersectional multilevel analysis of individual heterogeneity and discriminatory accuracy (MAIHDA) to study individual longitudinal trajectories, with application to mental health in the UK. Social Science & Medicine, 351, 116955. doi:10.1016/j.socscimed.2024.116955