High-dimensional Multivariate Mediation: with Application to Neuroimaging Data

1 Introduction

• High-dimensional mediation has received scarce attention, despite growing applications measuring massive numbers of variables and a need for methods when mediators greatly outnumber observations.The paper notes brain imaging, genetics, epidemiology, and public health as examples of such applications.Source: §1 · pp. 3–4
• In the thermal-pain fMRI setting, the question is whether brain response mediates the relationship between stimulus temperature and reported pain.Brain activity is measured across many voxels or regions.Source: §1 · p. 5
• Unlike voxel-wise mediation, which ignores relationships between voxels, DM seeks multivariate linear combinations and bases dimension reduction on the complete treatment–mediator–outcome relationship.The paper contrasts DM with partial least squares, which bases dimension reduction on the correlation between a response and explanatory variables.Source: §1 · pp. 5–6
• DMs combine high-dimensional mediator activity into orthogonal components ranked by the proportion of the LSEM likelihood each accounts for.The components may provide more interpretable mediators as linear combinations of voxels and reduce the number needed to estimate the indirect effect.Source: §1 · p. 6
• The paper develops an estimation algorithm and inference method for DMs, then illustrates the approach with simulations and a thermal-pain fMRI application.The proposed framework is presented as extending beyond LSEMs to other model-based mediation approaches.Source: §1 · pp. 4, 6

2 A Mutivariate Causal Mediation Model

• The total effect is decomposed into direct and indirect effects, separating the outcome contrast not through mediators from the contrast operating through them.The section defines controlled and natural direct effects and the natural indirect effect using potential outcomes.Source: §2 · p. 7
• The observed-data LSEM assumes linear relationships among treatment, mediators, and outcome, with no treatment–mediator interaction in the outcome regression.The controlled and natural direct effects are equivalent under the assumed absence of treatment–mediator interaction.Source: §2 · p. 8
• Under the stated assumptions, regression quantities identify causal mediation effects; when assumptions fail or counterfactuals are undefined, they may still support exploratory identification of potential mediators.The assumptions rule out confounding for treatment–outcome, mediator–outcome, and treatment–mediator relationships, including mediator–outcome confounding affected by treatment.Source: §2 · pp. 7–8
• When the mediator count grows, estimates become unstable and can greatly exceed the sample size, motivating an orthogonal transformation that simplifies and stabilizes estimation.After transformation, a series of LSEMs can estimate direct and indirect effects, one for each transformed mediator.Source: §2 · p. 9
• The transformation is chosen to rank transformed mediators by their share of the full LSEM likelihood, potentially yielding interpretable voxel combinations and fewer mediators for estimating the indirect effect.The figure depicts correlated mediators before transformation and independent mediators afterward.Source: §2 · p. 9, Figure 2 · p. 9

3.1 Model Formulation

• The model transforms the mediator vector into orthogonal variables so its linear structural equation model can be estimated separately for each transformed mediator.The first direction is estimated jointly with the model parameters.Source: §3.1 · pp. 10–11
• DMs are orthogonal, unit-norm mediator combinations selected to maximize the conditional joint likelihood, with each later direction constrained by earlier ones.The signs of the DMs are unidentifiable.Source: §3.1 · pp. 11–12

3.2 Estimation

• An iterative procedure jointly estimates the first DM and path parameters; the authors establish consistency and asymptotic normality for the estimated parameters.The asymptotic results are given in Supplemental Material Theorems 1 and 2.Source: §3.2 · p. 13

3.3 Higher Order Directions of Mediation

• For higher-order DMs, subtraction with Gram–Schmidt projections is computationally more efficient than adding penalty parameters and gives a good approximation.The penalty-parameter approach is described as more likely to achieve global maxima.Source: §3.3 · pp. 13–14

3.4 High-dimensional Directions of Mediation

• As the number of mediators grows, the low-dimensional estimation procedure becomes cumbersome, motivating a matrix decomposition to reduce the problem’s dimension.The method uses a generalized version of Population Value Decomposition (PVD).Source: §3.4 · p. 14

3.4.1 Generalized PVD

• GPVD allows subjects to have different numbers of trials while retaining the original PVD’s dimension-reduction benefits.It represents subject-specific data using subject-specific coefficients and a population-specific matrix; the reduced dimension B can be chosen using a criterion such as total variance explained.Source: §3.4.1 · pp. 14–15

3.4.2 Estimation using GPVD

• Estimation using GPVD fits the reduced mediator representation and retrieves the high-dimensional DM estimate through the generalized inverse.The reduced direction is defined by projecting the original direction through D.Source: §3.4.2 · p. 16

• In low-dimensional settings, variance estimates for the first DM and path coefficients can be obtained using Theorems 1 and 2 from the Supplemental Material.Source: §4 · p. 16
• In high-dimensional settings, generalized-inverse variance estimation is underestimated because D is random, while covariance estimation for D is computationally infeasible; the authors therefore use bootstrap inference.The cited dimensions are D of size B × p, with B ≪ p.Source: §4 · p. 16
• The bootstrap procedure resamples rows of the reduced mediator matrix, maps each bootstrap DM back to the high-dimensional space, and repeats this to form a matrix of bootstrap estimates.The stacked estimates are J × p, with columns giving bootstrap values for individual voxels.Source: §4 · pp. 16–17
• Because DM signs are unidentifiable, voxel-wise bootstrap distributions can be unimodal or bimodal and are modeled as mixed normals for p-value calculation.The passage identifies the mixtools package as an EM-based option for estimating mixture parameters.Source: §4 · p. 17

5.1 Simulation Set-up

• The simulation process generates mediator and treatment data, sets path coefficients and a first direction, then compares estimated parameters with their true values.The simulations repeatedly generated data under the specified model and recorded parameter estimates.Source: §5.1 · p. 18

5.2 Simulation Results

• For p = 10, estimates become more accurate and increasingly normal with smaller standard deviations as sample size rises; signs also become more consistent.For p = 3, results are shown for n = 10, 100, 500, and 1,000.Source: §5.2 · p. 20
• Turn-out rates increase with n for fixed p and with p for fixed n, with Table 1 covering n,p combinations across 1,000 simulations.The turn-out rate counts the estimating results produced out of a fixed number of simulations.Source: §5.2 · p. 21, Table 1 · p. 21
• At small sample sizes, λ(θ) can be poorly behaved and Newton–Raphson can fail at an intermediary step; performance improves with larger p or sufficiently large n.For p ∼ 3, the algorithm runs better with sufficiently large sample size, e.g. n ∼ 300; finer grid points improve performance at a computational-efficiency cost.Source: §5.2 · pp. 21–22
• Under the null, empirical p-values show adequate false-positive-rate control, though results are somewhat over-conservative.Figure 5 plots empirical against theoretical p-values, with 95% confidence bounds.Source: §5.2 · p. 22, Figure 5 · p. 22
• Across 100 bootstrap repetitions, the estimated first direction of mediation is consistent with the simulated signal in the 1,000 active voxels.The simulated mediator was designed to have strong activations in the center 1,000 voxels.Source: §5.2 · p. 22, Figure 6 · p. 23

6.1 Data Description

• The analysis uses 1,149 trial observations across 206,777 voxels from 33 healthy, right-handed participants.Each trial's temperature and pain rating are represented as vectors of length 1,149, while the activation matrix has one row per trial and one column per voxel.Source: §6.1 · pp. 23, 25
• Across nine runs, seven were passive and involved participants experiencing and rating heat stimuli, while two were regulation runs; stimuli used six temperature levels, and pain intensity was rated on a 0-to-100 scale.Source: §6.1 · pp. 23–24
• Trial-specific GLMs modeled thermal stimulation and rating periods, and their stimulus coefficients supplied the voxelwise activation measures.The functional images were motion- and slice-timing-corrected, normalized, smoothed, and high-pass filtered before analysis.Source: §6.1 · p. 24

6.2 Results

• The first direction has positive weights in ascending nociceptive pathways and negative weights in regions associated with processing pain context.Positive-weight regions include the anterior cingulate, insula, parietal operculum/S2, S1 hand area, and cerebellum; negative-weight regions include lateral prefrontal cortex, parahippocampal cortex, and ventral caudate.Source: §6.2 · p. 27
• The second direction assigns opposite weights to S1 and the parietal operculum/S2 and amygdala, separating nociceptive regions after accounting for the first direction.It also includes non-nociceptive regions that partly overlap with, and partly differ from, those in the first direction.Source: §6.2 · p. 27
• The third direction has positive weights in visual and superior cerebellar regions and parts of the hippocampus, with negative weights in prefrontal regions and the nucleus accumbens/ventral striatum.The authors also identify parahippocampal cortex and anterior insula/VLPFC as involved in this direction.Source: §6.2 · p. 27

7 Discussion

• The first direction is the mediator combination that maximizes the SEM likelihood, and each subsequent direction maximizes it subject to orthogonality to earlier directions.The discussion frames this as an approach to high-dimensional mediation analysis.Source: §7 · p. 27
• Causal interpretation depends on an untestable assumption that mediators behave as if randomly assigned, whose validity may differ across brain regions.The authors say this context-dependent assumption makes causal claims more difficult to assess in neuroimaging.Source: §7 · pp. 27–28
• The approach is intended for exploratory detection of potentially mediating region sets that can be examined in more targeted studies.The authors note that the framework's principles extend beyond LSEMs to other model-based mediation approaches.Source: §7 · p. 28
• The derivation assumes each subject is scanned under one condition, while the motivating fMRI data include multiple conditions; extending the causal model is left for future work.Such an extension could allow single-subject mediation studies when unit direct effects on mediators and unit total effects on outcomes are observed.Source: §7 · p. 28
• The authors identify measurement error in brain activation as a limitation and propose adapting the model to handle systematic mediator measurement errors.Source: §7 · p. 28
• The framework does not identify the signs of its estimates; the authors suggest Bayesian sign constraints or non-negativity constraints as possible responses.Mixed positive and negative mediator weights can cancel, so a zero weighted mediator does not necessarily imply no mediation effect.Source: §7 · pp. 28–29