Multiple mediation analysis for interval-valued data

Abstract

Mediation analysis is an important statistical approach to evaluate the relationships among observed variables. The most commonly used models for mediation analysis handle single-valued variables. However, there are several circumstances (e.g., dimensionality reduction of large datasets, clinical patient courses, repeated measures, masked data, uncertain data) in which the collected information can be represented more naturally by means of intervals. In these cases, standard mediation analyses can be ill-suited. Although interval-valued variables can be transformed into standard single-valued variables, such procedures may mask some relevant information provided by intervals. In this article, we present a novel and simple model (IMedA) to perform mediation analysis on interval-valued variables which is based on both the symbolic regression approach and the regression based mediation framework. We also generalize Stolzenberg’s decomposition of effects to cope with interval-valued data. We further introduce a specific variance based decomposition procedure to descriptively evaluate the sizes of such effects. Finally, to better highlight the IMedA features we apply our model to a real case study from behavioral contexts.

Appendices

Appendix A: Solutions for IMedA model

$$\begin{aligned} \begin{array}{ll}&{{\mathrm{vec}}}(\widehat{\mathbf A}^r) = (\mathbf {I}_k \otimes \mathbf {1}^T \mathbf {1})^{-1} \cdot (\mathbf {I}_k \otimes \mathbf {1})^T {{\mathrm{vec}}}({{\mathbf{M}^\mathbf{r}}} - \mathbf {X}\mathbf {\Xi } {\varvec{\Pi }} - \mathbf {1}\mathbf A^c {\varvec{\Pi }}); \end{array} \end{aligned}$$

(A1)

$$\begin{aligned} \begin{array}{ll} &{}{{\mathrm{vec}}}(\widehat{{\varvec{\Pi }}}) = \begin{bmatrix} (\mathbf {I}_k \otimes \mathbf {\Xi }^T \mathbf {X}^T \mathbf {X} \mathbf {\Xi }) + \\ (\mathbf {I}_k \otimes \mathbf {\Xi }^T \mathbf {X}^T \mathbf {1} \mathbf A^c) + \\ (\mathbf {I}_k \otimes \mathbf A^{cT} \mathbf {1}^T \mathbf {1} \mathbf A^c)\\ \end{bmatrix}^{-1} \cdot \begin{bmatrix} (\mathbf {I}_k \otimes \mathbf {X}\mathbf {\Xi })^T {{\mathrm{vec}}}(\mathbf {M^r - \mathbf {1}\mathbf A^r)} + \\ (\mathbf {I}_k \otimes \mathbf {1}\mathbf A^c)^T {{\mathrm{vec}}}(\mathbf {M^r - \mathbf {1}\mathbf A^r)} \end{bmatrix}; \end{array} \end{aligned}$$

(A2)

$$\begin{aligned} \begin{array}{ll} &{} {{\mathrm{vec}}}(\widehat{\mathbf \Xi }) = \begin{bmatrix} (\mathbf {I}_k \otimes \mathbf {X}^T\mathbf {X}) + \\ ({\varvec{\Pi }} \otimes \mathbf {X}^T\mathbf {X}) \end{bmatrix}^{-1} \cdot \begin{bmatrix} (\mathbf {I}_k \otimes \mathbf {X})^T {{\mathrm{vec}}}({\mathbf{M}^\mathbf{c}} - \mathbf {1}\mathbf A^c) +\\ ({\varvec{\Pi }} \otimes \mathbf {X})^T {{\mathrm{vec}}}({{\mathbf{M}^\mathbf{r}}} - \mathbf {1}\mathbf A^r - \mathbf {1}\mathbf A^c {\varvec{\Pi }}) \\ \end{bmatrix}; \end{array} \end{aligned}$$

(A3)

$$\begin{aligned} \begin{array}{ll} &{} {{\mathrm{vec}}}(\widehat{\mathbf A}^c) = \begin{bmatrix} (\mathbf {I}_k \otimes \mathbf {1}^T\mathbf {1}) + \\ ({\varvec{\Pi }} \otimes \mathbf {1}^T\mathbf {1}) \end{bmatrix}^{-1} \cdot \begin{bmatrix} (\mathbf {I}_k \otimes \mathbf {1})^T {{\mathrm{vec}}}({\mathbf{M}^\mathbf{c}} - \mathbf {X}\mathbf {\Xi }) +\\ ({\varvec{\Pi }} \otimes \mathbf {1})^T {{\mathrm{vec}}}({{\mathbf{M}^\mathbf{r}}} - \mathbf {1}\mathbf A^r - \mathbf {X}\mathbf {\Xi } {\varvec{\Pi }}) \\ \end{bmatrix}; \end{array} \end{aligned}$$

(A4)

$$\begin{aligned} \begin{array}{ll} &{} \widehat{\delta } = \begin{bmatrix} \alpha ^c \mathbf {1}^T\mathbf {1}\alpha ^c + 2\alpha ^c \mathbf {1}^T\mathbf {X}{\varvec{\beta }} + 2\alpha ^c \mathbf {1}^T\mathbf {M}{\varvec{\gamma }} +\\ {\varvec{\beta }}^T \mathbf {X}^T\mathbf {1}\alpha ^c + 2{\varvec{\beta }}^T \mathbf {X}^T\mathbf {X}{\varvec{\beta }} + 2{\varvec{\beta }}^T \mathbf {X}^T\mathbf {M}{\varvec{\gamma }} +\\ {\varvec{\gamma }}^T \mathbf {M}^T\mathbf {1}\alpha ^c + 2{\varvec{\gamma }}^T \mathbf {M}^T\mathbf {X}{\varvec{\beta }} + 2{\varvec{\gamma }}^T \mathbf {M}^T\mathbf {M}{\varvec{\gamma }} \end{bmatrix}^{-1} \cdot \begin{bmatrix} \alpha ^c \mathbf {1}^T (\mathbf {y}^r - \mathbf {1}\alpha ^r ) + \\ {\varvec{\beta }}^T \mathbf {X}^T (\mathbf {y}^r - \mathbf {1}\alpha ^r ) + \\ {\varvec{\gamma }}^T \mathbf {M}^T (\mathbf {y}^r - \mathbf {1}\alpha ^r ) \end{bmatrix}; \end{array}\nonumber \\ \end{aligned}$$

(A5)

$$\begin{aligned} \widehat{{\varvec{\beta }}} = \begin{bmatrix} \mathbf {X}^T \mathbf {X} + \\ \delta \mathbf {X}^T \mathbf {X} \delta \end{bmatrix}^{-1} \cdot \begin{bmatrix} \mathbf {X}^T (\mathbf {y}^c - \mathbf {1}\alpha ^c - \mathbf {M}{\varvec{\gamma }}) + \\ \mathbf {X}^T [~ \mathbf {y}^r - \mathbf {1}\alpha ^r + (-\mathbf {1}\alpha ^c - \mathbf {M}{\varvec{\gamma }})\delta ~]\delta \end{bmatrix}; \end{aligned}$$

(A6)

$$\begin{aligned} \widehat{{\varvec{\gamma }}^c} = \begin{bmatrix} {\mathbf{M}^\mathbf{c}}^T {\mathbf{M}^\mathbf{c}} + \\ \delta {\mathbf{M}^\mathbf{c}}^T {\mathbf{M}^\mathbf{c}} \delta \end{bmatrix}^{-1} \cdot \begin{bmatrix} {\mathbf{M}^\mathbf{c}}^T (\mathbf {y}^c - \mathbf {1}\alpha ^c - \mathbf {X}{\varvec{\beta }} - {{\mathbf{M}^\mathbf{r}}}{\varvec{\gamma }}^{\varvec{r}}) + \\ {\mathbf{M}^\mathbf{c}}^T [~ \mathbf {y}^r - \mathbf {1}\alpha ^r + (-\mathbf {1}\alpha ^c - \mathbf {X}{\varvec{\beta }} - {{\mathbf{M}^\mathbf{r}}}{\varvec{\gamma }}^{\varvec{r}})\delta ~]\delta \end{bmatrix}; \end{aligned}$$

(A7)

$$\begin{aligned} \widehat{{\varvec{\gamma }}^r} = \begin{bmatrix} {{\mathbf{M}^\mathbf{r}}}^T {{\mathbf{M}^\mathbf{r}}} + \\ \delta {{\mathbf{M}^\mathbf{r}}}^T {{\mathbf{M}^\mathbf{r}}} \delta \end{bmatrix}^{-1} \cdot \begin{bmatrix} {{\mathbf{M}^\mathbf{r}}}^T (\mathbf {y}^r - \mathbf {1}\alpha ^r - \mathbf {X}{\varvec{\beta }} - {\mathbf{M}^\mathbf{c}}{\varvec{\gamma }}^{\varvec{c}}) + \\ {{\mathbf{M}^\mathbf{r}}}^T [~ \mathbf {y}^r - \mathbf {1}\alpha ^r + (-\mathbf {1}\alpha ^r - \mathbf {X}{\varvec{\beta }} - {\mathbf{M}^\mathbf{c}}{\varvec{\gamma }}^{\varvec{c}})\delta ~]\delta \end{bmatrix}; \end{aligned}$$

(A8)

$$\begin{aligned} \widehat{\alpha }^c = \frac{1}{n(1 + {\delta }^2)} \cdot \begin{bmatrix} \mathbf {1}^T (\mathbf {y}^c - \mathbf {X}{\varvec{\beta }} - \mathbf {M}{\varvec{\gamma }}) + \\ \mathbf {1}^T [~ \mathbf {y}^r - \mathbf {1}\alpha ^r + (-\mathbf {X}{\varvec{\beta }} - \mathbf {M}{\varvec{\gamma }})\delta ~]\delta \end{bmatrix}; \end{aligned}$$

(A9)

$$\begin{aligned} \widehat{\alpha }^r = \frac{1}{n} \cdot \mathbf {1}^T [ ~ {{\mathbf{y}^\mathbf{r}}} - (\mathbf {1}\alpha ^c - \mathbf {X}{\varvec{\beta }} - \mathbf {M}{\varvec{\gamma }})\delta ]; \end{aligned}$$

(A10)

where vec(.) is the linear operator that converts a $n \times k$ matrix into a $kn \times 1$ vector, $\otimes$ denotes the Kronecker product, $\mathbf {I}_k$ is a $k \times k$ identity matrix whereas $\mathbf {1}$ is a $n\times k$ matrix of all ones.

Appendix B: Decomposition of effects for IMedA model

In order to derive direct and indirect effects for the IMedA model, we proceed as follows. Consider the regression systems $\mathcal {S}_1$ and $\mathcal {S}_2$ shown in Eq. 1:

$$\begin{aligned} \mathcal {S}_1: {\left\{ \begin{array}{ll} {\mathbf{M}^\mathbf{c}} = \mathbf {1} \mathbf {A}^c + \mathbf {X} \mathbf {\Xi } + \mathbf {E}^c\\ {{\mathbf{M}^\mathbf{r}}} = \mathbf {1} \mathbf {A}^r + (\mathbf {1} \mathbf {A}^c + \mathbf {X} \mathbf {\Xi }) {\varvec{\Pi }} + \mathbf {E^r} \end{array}\right. } \end{aligned}$$

$$\begin{aligned} \mathcal {S}_2: {\left\{ \begin{array}{ll} {{\mathbf{y}^\mathbf{c}}} = \mathbf {1} {\alpha }^c + \mathbf {X} {\varvec{\beta }} + {\mathbf{M}^\mathbf{c}} {\varvec{\gamma }}^{\varvec{c}} + {{\mathbf{M}^\mathbf{r}}} {\varvec{\gamma }}^{\varvec{r}} + {\varvec{\epsilon }}^{\varvec{c}}\\ {{\mathbf{y}^\mathbf{r}}} = \mathbf {1} {\alpha }^r + (\mathbf {1} {\alpha }^c + \mathbf {X} {\varvec{\beta }} + {\mathbf{M}^\mathbf{c}} {\varvec{\gamma }}^{\varvec{c}} + {{\mathbf{M}^\mathbf{r}}} {\varvec{\gamma }^{\varvec{r}}})\delta + {\varvec{\epsilon }}^{\varvec{r}}\\ \end{array}\right. } \end{aligned}$$

Firstly, substitute the equations of ${\mathbf{M}^\mathbf{c}}$ and ${{\mathbf{M}^\mathbf{r}}}$ into ${{\mathbf{y}^\mathbf{c}}}$ and ${{\mathbf{y}^\mathbf{r}}}$, as follows:

$$\begin{aligned} \mathcal {S'}_2: {\left\{ \begin{array}{ll} {{\mathbf{y}^\mathbf{c}}} = \mathbf {1} {\alpha }^c + \mathbf {X} {\varvec{\beta }} + \, [\mathbf {1} \mathbf {A}^c + \mathbf {X} \mathbf {\Xi } + \mathbf {E}^c]{\varvec{\gamma }}^c + \\ \quad \quad +\, [\mathbf {1} \mathbf {A}^r + (\mathbf {1} \mathbf {A}^c + \mathbf {X} \mathbf {\Xi }) {\varvec{\Pi }} + \mathbf {E^r}]\gamma ^r + {\varvec{\epsilon }}^c\\ {{\mathbf{y}^\mathbf{r}}} = \mathbf {1} {\alpha }^r + (\mathbf {1} {\alpha }^c + \mathbf {X} {\varvec{\beta }} + [\mathbf {1} \mathbf {A}^c + \mathbf {X} \mathbf {\Xi } + \mathbf {E}^c]{\varvec{\gamma }}^c + \\ \quad \quad +\, [\mathbf {1} \mathbf {A}^r + (\mathbf {1} \mathbf {A}^c + \mathbf {X} \mathbf {\Xi }) {\varvec{\Pi }} + \mathbf {E^r}]\gamma ^r)\delta ^r_y + {\varvec{\epsilon }}^{\varvec{r}} \end{array}\right. } \end{aligned}$$

(B1)

Multiplying through and expanding terms, using a little algebra, we obtain the following reduced form system $\mathcal {S'}_2$:

$$\begin{aligned} \mathcal {S'}_2: {\left\{ \begin{array}{ll} \begin{array}{ll} {{\mathbf{y}^\mathbf{c}}} &{} = \mathbf {1} {\alpha }^c + \mathbf {1} \mathbf {A}^c {\varvec{\gamma }}^c + \mathbf {1} \mathbf {A}^r {\varvec{\gamma }}^r + \mathbf {1} \mathbf {A}^c {\varvec{\Pi }} {\varvec{\gamma }}^r + {{\mathbf{x}^\mathbf{c}}} [{\varvec{\beta }}^c + {\varvec{\xi }}^c ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r)]\\ &{}\quad +\, {{\mathbf{x}^\mathbf{r}}} [{\varvec{\beta }}^r + {\varvec{\xi }}^r ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r)] + \mathbf {E}^c {\varvec{\gamma }}^c + \mathbf {E}^r {\varvec{\gamma }}^r + \varvec{\epsilon }^c \end{array}\\ \begin{array}{ll} {{\mathbf{y}^\mathbf{r}}} &{} = \mathbf {1} {\alpha }^r + \mathbf {1} {\alpha }^c\delta + \mathbf {1} \mathbf {A}^c {\varvec{\gamma }}^c\delta + \mathbf {1} \mathbf {A}^r {\varvec{\gamma }}^r\delta + \mathbf {1} \mathbf {A}^c {\varvec{\Pi }} {\varvec{\gamma }}^r\delta +\\ &{} \quad +\, {{\mathbf{x}^\mathbf{c}}} [{\varvec{\beta }}^c + {\varvec{\xi }}^c ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r)]\delta + {{\mathbf{x}^\mathbf{r}}} [{\varvec{\beta }}^r + {\varvec{\xi }}^r ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r)]\delta \\ &{} \quad +\, \mathbf {E}^c {\varvec{\gamma }}^c\delta + \mathbf {E}^r {\varvec{\gamma }}^r\delta + \varvec{\epsilon }^r \end{array} \end{array}\right. } \end{aligned}$$

(B2)

Next, taking the partial derivatives of ${{\mathbf{y}^\mathbf{c}}}$ and ${{\mathbf{y}^\mathbf{r}}}$ with respect to ${{\mathbf{x}^\mathbf{c}}}$ and ${{\mathbf{x}^\mathbf{r}}}$ we have the equations for the total effect (TE) of the model, as follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial ~ {{\mathbf{y}^\mathbf{c}}}}{\partial ~ {{\mathbf{x}^\mathbf{c}}}} = \delta \beta ^c + {\varvec{\xi }}^c ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r) \frac{\partial ~ {{\mathbf{y}^\mathbf{c}}}}{\partial ~ {{\mathbf{x}^\mathbf{r}}}} = \delta \beta ^r + {\varvec{\xi }}^r ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r)\\ \frac{\partial ~ {{\mathbf{y}^\mathbf{r}}}}{\partial ~ {{\mathbf{x}^\mathbf{c}}}} = \delta \beta ^c + {\varvec{\xi }}^c ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r)\delta \frac{\partial ~ {{\mathbf{y}^\mathbf{r}}}}{\partial ~ {{\mathbf{x}^\mathbf{r}}}} = \delta \beta ^r + {\varvec{\xi }}^r ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r)\delta \\ \end{array}\right. } \end{aligned}$$

Finally, collecting and simplifying the ensuing terms, we obtain the following equations for TE:

$$\begin{aligned} \mathrm{TE}^{y^c}= & {} {\beta }^c_y + {\beta }^r_y + ({{\varvec{\beta }}}^c_m \circ {{\varvec{\gamma }}}^{c~T}_m)\mathbf {1}_m + ({{\varvec{\beta }}}^c_m \circ {{\varvec{\gamma }}}^{r~T}_m \circ {{\varvec{\Pi }}})\mathbf {1}_m + ({{\varvec{\beta }}}^r_m \circ {{\varvec{\gamma }}}^{c~T}_m )\mathbf {1}_m + \nonumber \\&+\, ({{\varvec{\beta }}}^r_m \circ {{\varvec{\gamma }}}^{r~T}_m \circ {{\varvec{\Pi }}})\mathbf {1}_m \nonumber \\ \mathrm{TE}^{y^r}= & {} \delta ^r_y [{\beta }^c_y + {\beta }^r_y + ({{\varvec{\beta }}}^c_m \circ {{\varvec{\gamma }}}^{c~T}_m)\mathbf {1}_m + ({{\varvec{\beta }}}^c_m \circ {{\varvec{\gamma }}}^{r~T}_m \circ {{\varvec{\Pi }}})\mathbf {1}_m + \nonumber \\&+\, ({{\varvec{\beta }}}^r_m \circ {{\varvec{\gamma }}}^{c~T}_m )\mathbf {1}_m + ({{\varvec{\beta }}}^r_m \circ {{\varvec{\gamma }}}^{r~T}_m \circ {{\varvec{\Pi }}})\mathbf {1}_m] \end{aligned}$$

(B3)

which are in the general form of $\mathrm{TE} = \mathrm{DE}_c + \mathrm{DE}_r + \mathrm{IE}_{c/c} + \mathrm{IE}_{c/r} + \mathrm{IE}_{r/c} + \mathrm{IE}_{r/r}$. Note that the equation $\mathrm{TE}^{y^r}$ for ${{\mathbf{y}^\mathbf{r}}}$ is obtained as linear combination of $\mathrm{TE}^{y^c}$ through the parameter $\delta$.

Appendix C: Decomposition of variance for IMedA model

Considering the reduced form system $\mathcal {S'}_2$:

$$\begin{aligned} \mathcal {S'}_2: {\left\{ \begin{array}{ll} \begin{array}{ll} {{\mathbf{y}^\mathbf{c}}} &{} = \mathbf {1} {\alpha }^c + \mathbf {1} \mathbf {A}^c {\varvec{\gamma }}^c + \mathbf {1} \mathbf {A}^r {\varvec{\gamma }}^r + \mathbf {1} \mathbf {A}^c {\varvec{\Pi }} {\varvec{\gamma }}^r + {{\mathbf{x}^\mathbf{c}}} [{\varvec{\beta }}^c + {\varvec{\xi }}^c ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r)]\\ &{}\quad +\, {{\mathbf{x}^\mathbf{r}}} [{\varvec{\beta }}^r + {\varvec{\xi }}^r ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r)] + \mathbf {E}^c {\varvec{\gamma }}^c + \mathbf {E}^r {\varvec{\gamma }}^r + {\varvec{\epsilon }}^c \end{array}\\ \begin{array}{ll} {{\mathbf{y}^\mathbf{r}}} &{} = \mathbf {1} {\alpha }^r + \mathbf {1} {\alpha }^c\delta + \mathbf {1} \mathbf {A}^c {\varvec{\gamma }}^c\delta + \mathbf {1} \mathbf {A}^r {\varvec{\gamma }}^r\delta + \mathbf {1} \mathbf {A}^c {\varvec{\Pi }} {\varvec{\gamma }}^r\delta +\\ &{}\quad +\,{{\mathbf{x}^\mathbf{c}}} [{\varvec{\beta }}^c + {\varvec{\xi }}^c ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r)]\delta + {{\mathbf{x}^\mathbf{r}}} [{\varvec{\beta }}^r + {\varvec{\xi }}^r ({\varvec{\gamma }}^c + {\varvec{\Pi }} {\varvec{\gamma }}^r)]\delta \\ &{}\quad +\, \mathbf {E}^c {\varvec{\gamma }}^c\delta + \mathbf {E}^r {\varvec{\gamma }}^r\delta + {\varvec{\epsilon }}^r \end{array} \end{array}\right. } \end{aligned}$$

(C1)

the following identities hold:

$$\begin{aligned} {{\mathrm{var}}}({{\mathbf{y}^\mathbf{c}}})= & {} {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{c}}}{\beta }^c) +\nonumber \\&+ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{r}}}{\beta }^r) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{cT})\mathbf {1}_k) + \nonumber \\&+ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k)+ \nonumber \\&+ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{cT} )\mathbf {1}_k) +{{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k) + \nonumber \\&+ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},\mathbf {E}^c {\varvec{\gamma }}^c) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},\mathbf {E}^r {\varvec{\gamma }}^r) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{\varvec{\epsilon }}^c)\nonumber \\ {{\mathrm{var}}}({{\mathbf{y}^\mathbf{r}}})= & {} {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},{\delta }{{\mathbf{x}^\mathbf{c}}}{\beta }^c) +{{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},{\delta }{{\mathbf{x}^\mathbf{r}}}{\beta }^r) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},{\delta }{{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{cT})\mathbf {1}_k) + \nonumber \\&+ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},{\delta }{{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},{\delta }{{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{cT} )\mathbf {1}_k) +\nonumber \\&+ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},{\delta }{{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k) + \nonumber \\&+ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\mathbf {E}^c {\varvec{\gamma }}^c \delta ) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\mathbf {E}^r {\varvec{\gamma }}^r \delta ) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{\varvec{\epsilon }}^r) \end{aligned}$$

(C2)

after noticing that:

$$\begin{aligned}&\displaystyle {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},\mathbf {1} {\alpha }^c + \mathbf {1} \mathbf {A}^c {\varvec{\gamma }}^c + \mathbf {1} \mathbf {A}^r {\varvec{\gamma }}^r + \mathbf {1} \mathbf {A}^c {\varvec{\Pi }} {\varvec{\gamma }}^r) = 0 \\&\displaystyle {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},\mathbf {1} {\alpha }^r + \mathbf {1} {\alpha }^c\delta + \mathbf {1} \mathbf {A}^c {\varvec{\gamma }}^c\delta + \mathbf {1} \mathbf {A}^r {\varvec{\gamma }}^r\delta + \mathbf {1} \mathbf {A}^c {\varvec{\Pi }} {\varvec{\gamma }}^r\delta ) = 0 \end{aligned}$$

where ${{\mathrm{cov}}}(.)$ and ${{\mathrm{var}}}(.)$ indicates the covariance and variance operators, $\circ$ denotes the Hadamard product whereas $\mathbf {1}_k$ is a $k\times 1$ vector of all ones. The following properties hold:

$$\begin{aligned} \begin{array}{ll} &{} \begin{bmatrix} {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{c}}}{\beta }^c) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{r}}}{\beta }^r) + \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{cT})\mathbf {1}_k) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k)+ \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{cT} )\mathbf {1}_k) +{{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k) + \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},\mathbf {E}^c {\varvec{\gamma }}^c) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},\mathbf {E}^r {\varvec{\gamma }}^r) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{\varvec{\epsilon }}^c) \end{bmatrix} \begin{bmatrix} {{\mathrm{var}}}({{\mathbf{y}^\mathbf{c}}}) \end{bmatrix}^{-1} = 1\\ &{} \begin{bmatrix} {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{c}}}{\beta }^c) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{r}}}{\beta }^r) + \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{cT})\mathbf {1}_k) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k)+ \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{cT} )\mathbf {1}_k) +{{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k) + \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\mathbf {E}^c {\varvec{\gamma }}^c\delta ) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\mathbf {E}^r {\varvec{\gamma }}^r\delta ) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},{\varvec{\epsilon }}^r) \end{bmatrix} \cdot \begin{bmatrix} {{\mathrm{var}}}({{\mathbf{y}^\mathbf{r}}}) \end{bmatrix}^{-1} = 1 \end{array}\nonumber \\ \end{aligned}$$

(C3)

$$\begin{aligned} \begin{array}{ll} &{} \begin{bmatrix} {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{c}}}{\beta }^c) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{r}}}{\beta }^r) + \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{cT})\mathbf {1}_k) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k)+ \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{cT} )\mathbf {1}_k) +{{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},{{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k) + \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},\mathbf {E}^c {\varvec{\gamma }}^c) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{c}}},\mathbf {E}^r {\varvec{\gamma }}^r) \end{bmatrix} \cdot \begin{bmatrix} {{\mathrm{var}}}({{\mathbf{y}^\mathbf{c}}}) \end{bmatrix}^{-1} \approx \omega _1 \\ &{} \begin{bmatrix} {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{c}}}{\beta }^c) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{r}}}{\beta }^r) + \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{cT})\mathbf {1}_k) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{c}}}({{\varvec{\xi }}}^c \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k)+ \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{cT} )\mathbf {1}_k) +{{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\delta {{\mathbf{x}^\mathbf{r}}}({{\varvec{\xi }}}^r \circ {{\varvec{\gamma }}}^{rT} \circ {{\varvec{\Pi }}}^{T})\mathbf {1}_k) + \\ {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\mathbf {E}^c {\varvec{\gamma }}^c\delta ) + {{\mathrm{cov}}}({{\mathbf{y}^\mathbf{r}}},\mathbf {E}^r {\varvec{\gamma }}^r\delta ) \end{bmatrix} \cdot \begin{bmatrix} {{\mathrm{var}}}({{\mathbf{y}^\mathbf{r}}}) \end{bmatrix}^{-1} \approx \omega _2 \end{array}\nonumber \\ \end{aligned}$$

(C4)

where $\omega _1 = \Vert {{\mathbf{y}^\mathbf{c}}} - {{\mathbf{y}^\mathbf{c}}}^* \Vert ^2 /~ \Vert {{\mathbf{y}^\mathbf{c}}} - {\overline{\mathbf{y}^\mathbf{c}}} \Vert ^2$ whereas $\omega _2 = \Vert {{\mathbf{y}^\mathbf{r}}} - {{\mathbf{y}^\mathbf{r}}}^* \Vert ^2 /~ \Vert {{\mathbf{y}^\mathbf{r}}} - {\overline{\mathbf{y}^\mathbf{r}}} \Vert ^2$. Note that ${\bar{\mathbf{y^c}}}$ and ${\bar{\mathbf{y^r}}}$ are $n\times 1$ vectors containing the mean values of $\mathbf {{y^c}}$ and $\mathbf {{y^r}}$ whereas ${{\mathbf{y}^\mathbf{c}}}^*$ and ${{\mathbf{y}^\mathbf{r}}}^*$ refers to the estimated reduced equations in $\mathcal {S'}_2$ without considering the residual terms ${\varvec{\epsilon }}^c$ and ${\varvec{\epsilon }}^r$. Note also that the terms in the right side of Eq. C4 denotes the variance explained by the reduced system $\mathcal {S'}_2$.

Appendix D: Bias-corrected and accelerated (BCa) bootstrap procedure

Bias-corrected and accelerated (BCa) is a powerful bootstrap procedure usually adopted in mediation analysis. More precisely, Q samples (with $Q\ge 1000$) of size n are row-wise randomly drawn (with replacement) from the original matrices ${\mathbf{M}^\mathbf{c}}$, ${{\mathbf{M}^\mathbf{r}}}$ and original vectors ${{\mathbf{y}^\mathbf{c}}}$, ${{\mathbf{y}^\mathbf{r}}}$. For each $q-$th sample, the mediation parameters are estimated by applying the IMedA procedure on the sample matrices ${\mathbf{M}^\mathbf{c}}_q$, ${{\mathbf{M}^\mathbf{r}}}_q$ and vectors ${{\mathbf{y}^\mathbf{c}}}_q$,${{\mathbf{y}^\mathbf{r}}}_q$. These steps are repeated for Q times. The ensuing sample parameter distributions are then used for computing the standard errors or BCa based confidence intervals (95% CIs) for every estimated parameter in the model. In what follows we briefly describe how the BCa based CIs can be obtained. For the sake of generality, considering the i-th parameter of $\widehat{{\varvec{\gamma }}^c}_i$ with $\widehat{{\varvec{\gamma }}^c}_i^* = \widehat{ {\varvec{\gamma }}^c}_{i1}, \widehat{ {\varvec{\gamma }}^c}_{i2}, \ldots , \widehat{ {\varvec{\gamma }}^c}_{iQ}$ denoting its empirical distribution. The 95% BCa based confidence interval for such parameter takes the form of $[~ \widehat{ {\varvec{\gamma }}^c}^*_{i_g}, \widehat{ {\varvec{\gamma }}^c}^*_{i_v}~]$, where g and v are, in turn, computed as follows:

$$\begin{aligned} g = n \cdot \Phi _\mathcal {N} \biggl ( z_{\phi _2} + \frac{z_{\phi _2} + z_{\alpha /2}}{1-\phi _1 (z_{\phi _2} + z_{\alpha /2})} \biggl ) \quad \quad v = n \cdot \Phi _\mathcal {N} \biggl ( z_{\phi _2} + \frac{z_{\phi _2} + z_{1-\alpha /2}}{1-\phi _1 (z_{\phi _2} + z_{1-\alpha /2})} \biggl ) \end{aligned}$$

where $\Phi _\mathcal {N}$ is the cumulative normal distribution function, $z_{\alpha /2} = -1.96$ and $z_{1-\alpha /2} = 1.96$ for $\alpha =0.05$ whereas $z_{\phi _2} = \Phi _\mathcal {N}^{-1}(\phi _1,0,1)$, with $\Phi _\mathcal {N}^{-1}$ being the inverse cumulative normal distribution function. Note that the term $z_{\phi _2}$ measures the median bias of the bootstrap distribution $\widehat{ \gamma }^*_i$ where $\phi _2$ is computed as follows:

$$\begin{aligned} \phi _2 = \frac{1}{Q}\sum _{\widehat{{\varvec{\gamma }}^c}^*_i ~<~ \underline{\widehat{{\varvec{\gamma }}^c}}} \widehat{{\varvec{\gamma }}^c}^*_i \quad \quad \quad \mathrm{where } \ \underline{\widehat{{\varvec{\gamma }}^c}} = \frac{\sum _{i=1}^Q \widehat{{\varvec{\gamma }}^c}^*_i}{Q} \end{aligned}$$

whereas, on the contrary, the acceleration term $\phi _1$, which measures the rate of change of the standard deviation of $\widehat{ \gamma }^*_i$, is computed as follows:

$$\begin{aligned} \phi _1 = \sum _{i=1}^n (\widehat{{\varvec{\gamma }}^c}^*_i - \underline{\widehat{{\varvec{\gamma }}^{c}}})^3 \cdot \biggl [ \biggl ( 6 \sum _{i=1}^n (\widehat{{\varvec{\gamma }}^c}^*_i - \underline{\widehat{{\varvec{\gamma }}^c}})^2 \biggl )^{\frac{3}{2}} \biggl ]^{-1} \end{aligned}$$

Note that, when $\phi _1 = \phi _2 = 0$ the BCa based CIs simply reduce to the standard percentile CIs.