Longitudinal MAIHDA: intersectional inequalities over time

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:

\[ 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:

\[ \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 $\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:

\[\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):

\[\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 $\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.

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 $(\text{outcome}/\text{time})^2$ while $\sigma^2_e$ is in $\text{outcome}^2$, so the sum would be dimensionally meaningless. (The headline VPC is well formed because $a(t)'\Sigma a(t)$ returns to $\text{outcome}^2$.) The intercept VPC drops it for symmetry with the slope.

Both are evaluated at the baseline $t_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).
#> Variance basis: as fitted (REML for Gaussian lmer, matching summary())
#> 
#> 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

Both PCVs compare REML variance estimates, and such an estimate cannot fall below zero: the stratum block must remain a covariance matrix. Where the strata differ little after adjustment next to how precisely one stratum’s own baseline or slope is estimated, the adjusted model’s variance is therefore too large on average while the null model’s is not, and the PCV is understated on average; single values scatter widely and can be negative. PCV_slope is the more exposed, because strata usually differ less in slope than in level. The yardstick is the sampling variance of one stratum’s slope, about $(\sigma^2_{v1} + \sigma^2_e / \sum_k (t_k - \bar t)^2) / n_j$ for a stratum of $n_j$ individuals observed at times $t_k$, with $\sigma^2_{v1}$ the slope variance between individuals. Here it is about 0.0014 (0.0006 to 0.0025 across the twelve strata) against a null between-stratum slope variance of 0.043, and the shortfall is negligible; with a dozen strata it shows from about a third of the null variance, and when the two are of similar size the data do not support a decomposition of the slope. More strata, larger strata and more waves reduce it.

With engine = "brms" each variance in this decomposition is the posterior median of its per-draw value: the between-stratum variance at the baseline, the slope variance, and the variance at each time behind the trajectory plot. The covariance blocks in summary() are posterior means, so they do not reproduce those numbers. The posterior of a between-stratum variance is right-skewed when the strata are few, and more so in the adjusted model: its variance is smaller next to the sampling noise of the stratum estimates, and it spends degrees of freedom on the main effects and their time interactions. A brms PCV is therefore understated on average with few strata – by more when the additive share is moderate, when the strata are small, or, for PCV_slope, when the strata differ little in slope next to how precisely a stratum’s slope is estimated – and can come out negative where the true share is positive. With around a dozen strata the lme4 decomposition is less biased, though neither is precise with so few strata; the understatement shrinks as strata are added.

Scope and cautions

Time coding. You can code time however is convenient for your workflow (waves 0, 1, 2, …, waves 10, 11, …, age, calendar year, or days since baseline). 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. With engine = "lme4", time whose span is not between 1 and 10 is also divided by a power of ten for the fit (days over five years by 1000, with a message): counted in fine units, lme4’s optimisers stop far from the optimum, mostly without a convergence warning. The growth model of a count or binary outcome is also refitted from more starting values, because glmer’s optimiser often stops with a variance at zero, short of the maximum. All results (the VPC trajectory, the PCV, the fixed effects, 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 (engine = "brms" reports an interval instead, very wide when the data say little). 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 are still reported. A rank 1 block has no between-stratum variance at one time, which print() names: the further that time lies from the observed times, the more both are carried by the intercept variance, which is still estimable; close to it the VPC is near zero and a PCV is not interpretable.

How often this happens depends on the number of strata and on how large the stratum slope variance is next to the sampling variance of one stratum’s slope (the yardstick in the decomposition section above). Where the strata do not diverge it is common in any design. Where they do, it is common with few strata unless their slope variance is well above the yardstick, and the slope variance of such a fit tends to be too small; more strata, more people per stratum and more measurement occasions make it rarer. With few strata, a rank-deficient block is therefore weak 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
```
With few strata this comparison has little power: it often favours dropping a slope variance that is there.

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. Where the strata do diverge, the constraint moves their slope differences into the individual-level slope variance, so the VPC understates the between-stratum share at later times and can fall where the true VPC rises; PCV_intercept can then also take in part of that divergence.
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