1.2.1 Single mediator

grViz("digraph regression {
graph [rankdir = LR bgcolor=transparent]

forcelabels=true;

node [shape = box, fontcolor=gray25 color=gray80]

node [fontname='Helvetica']
X; M;

node [fillcolor=gray90 style=filled]
Y;

edge [color=gray50 style=filled]
X -> Y [label='Direct']
X -> M
M -> Y
}")

1.2.2 Multiple mediators

grViz("digraph regression {
graph [rankdir = LR, bgcolor=transparent, layout=dot]

forcelabels=true;

node [shape = box, fontcolor=gray25 color=gray80]

node [fontname='Helvetica']
X; M1; M2;

node [fillcolor=gray90 style=filled]
Y;

edge [color=gray50 style=filled]
X -> Y [xlabel='e']
X -> M1 [label='a']
M1 -> M2 [label='b']
M2 -> Y [label='c']
M1 -> Y [label='d']
}")

In a multiple mediation model, one considers the direct effect (\(e\)) of a predictor \(X\) on a criterion \(Y\) after accounting for all mediating variables (here, \(M1\) and \(M2\)).

1.2.3 Specific indirect effect

Multiple mediation also exposes the new concept of specific indirect effects, which summarize the association between a predictor, \(X\) and outcome \(Y\) via one or more specific intervening pathways. For example, in the above model, we could measure how much of the relationship between \(X\) and \(Y\) is transmitted via the \(X \rightarrow M1 \rightarrow Y\), which would be \(a \cdot d\).

1.2.4 Total indirect effect

The total indirect effect is the sum of all indirect pathways between \(X\) and \(Y\) (i.e., excluding the \(X \rightarrow Y\) path).

1.2.5 Testing indirect effects

When we think about the magnitude of an indirect effect, it is best to conceptualize it as a product of the pathways between the predictor and outcome. Thus, the specific indirect effect of \(X\) on \(Y\) via the \(X \rightarrow M1 \rightarrow M2 \rightarrow Y\) pathway is the product of the corresponding parameter estimates \(a \cdot b \cdot c\).

2.2.1 Total effect

Adopting the notation of Andrew Hayes on conditional process analysis (e.g., Hayes, 2018, Introduction to mediation, moderation, and conditional process analysis: A regression-based approach), the magnitude of the relationship between the predictor \(X\) and outcome \(Y\) can be quantified by:

\[ Y = i_Y + cX + e_Y \]

where \(c\) is the regression coefficient that summarizes the total effect of \(X\) on \(Y\). If \(X\) is a 0/1 group variable, such as our CBT versus TAU example, then \(c\) quantifies the mean difference between the treatment (1) and control (0) group.

2.2.2 Effect of \(X\) on \(M\)

One key part of mediation is the strength of the relationship between the predictor, \(X\), and the mediating variable \(M\). If \(X\) does not influence \(M\), it is hard to imagine how it \(M\) could explain the relationship the between \(X\) and \(Y\) (Shrout and Bolger, 2002). For example, if CBT does not reduce the frequency of negative automatic thoughts, it is hard to imagine how automatic thoughts are a mediator of depression symptom change in CBT.

Thus, we wish to quantify the relationship:

\[ M = i_M + aX + e_M \]

where \(a\) is the strength of association between the predictor and mediator.

2.2.3 Effect of \(M\) on \(Y\)

For mediation to be plausible, we also expect that there must be a meaningful relationship between the mediator and the outcome. Otherwise, how could it be a proximate cause of the outcome?

\[ Y = i_Y + c'X + bM + e_Y \]

Notice that in this equation, the effect of the predictor/IV \(X\) on the outcome is typically included, representing the idea that the total effect reflects both a direct component (\(c'\)) and indirect component (\(a \cdot b\)). The \(c'\) coefficient represents the expected change in the outcome for \(Y\) a unit change in the predictor \(X\) at a given level of the mediator \(M\). In the two group case (1/0 \(X\)), the direct effect \(c'\) represents the expected mean difference between groups at the same (candidate) level of the mediator.

Likewise, the coefficient \(b\) represents the effect of the mediator \(M\) on the outcome \(Y\) at a given level of the predictor (aka controlling for \(X\)). In the case of a 0/1 \(X\), the \(b\) coefficient captures the effect of the mediator on the outcome in the control (0) group.

2.3.1 Sidebar: Nonparametric bootstrap resampling

Nonparametric bootstrap resampling is a powerful and general tool to infer the distribution of one or more parameters drawn from an unknown sampling distribution. The intuition is that if our sample represents a random subset of possible cases drawn from the population of interest, then we could potentially obtain a plausible replication sample by resampling the data with replacement. That is, if we have \(N = 200\) cases that represent a random sample of the population of interest, then if we drew observations from our sample one at a time with equal probability until we reached \(N = 200\), we would have a replication sample. This process necessarily means that we will select some observations several times and others not at all.

We repeat the resampling process many times (typically \(k = 1,000-10,000\)), estimating our model in each bootstrap sample. This gives us \(k\) estimates of our parameter(s) of interest. These \(k\) estimates form an empirical sampling distribution that allow us to infer the probability that the parameter falls in a given region (e.g., between .1 and .4) or differs from some value (e.g., > 0). Sidebar to sidebar: this approach has formal similarities to Bayesian statistics where we setup a sampler to draw plausible values for all parameters from the joint posterior distribution. In the same way, draws from this posterior can be used to make inferences on the value of parameters from unknown statistical distributions.

2.3.2 Sidebar: no need for X and Y to be correlated for mediation to exist

I won’t belabor the details (if you want to know more, see Cerin & MacKinnon, 2009; Hayes, 2009; Shrout & Bolger, 2002), but please know that the old Barron and Kenny approach to mediation required that X and Y be associated. Subsequent work has demonstrated that mediation can exist even when X and Y are not correlated. This typically happens when the direct (\(c'\)) and indirect effects (\(ab\)) are roughly equal, but opposite, in magnitude. In this scenario, the total effect, \(c\), will be around zero.

2.4.1 Plausible equivalent models

The problem of equivalent models can also rear its head when considering mediation. In particular, when considering full mediation, we can reverse the direction of the effects and achieve identical fit

grViz("digraph regression {
graph [rankdir = LR bgcolor=transparent]

forcelabels=true;

node [shape = box, fontcolor=gray25 color=gray80]

node [fontname='Helvetica']
X; M; Y;

edge [color=gray50 style=filled]
X -> M -> Y
}")

is equivalent to

grViz("digraph regression {
graph [rankdir = LR bgcolor=transparent]

forcelabels=true;

node [shape = box, fontcolor=gray25 color=gray80]

node [fontname='Helvetica']
X; M; Y;

edge [color=gray50 style=filled]
Y -> M -> X
}")

Let’s try it out!

medmodel <- '
Y ~ 0.3*1 + 0.3*X + 0.5*M
M ~ 0.5*1 + 0.4*X
X ~ 0.9*1
'

data <- lavaan::simulateData(medmodel, sample.nobs=500, meanstructure=TRUE)
psych::describe(data)

cor(data)

##       Y     M     X
## Y 1.000 0.552 0.382
## M 0.552 1.000 0.324
## X 0.382 0.324 1.000

m1_syntax <- '
  M ~ X
  Y ~ M
'

m1 <- sem(m1_syntax, data=data, fixed.x=FALSE, conditional.x=FALSE)
summary(m1)

## lavaan 0.6-3 ended normally after 10 iterations
## 
##   Optimization method                           NLMINB
##   Number of free parameters                          5
## 
##   Number of observations                           500
## 
##   Estimator                                         ML
##   Model Fit Test Statistic                      34.439
##   Degrees of freedom                                 1
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard Errors                             Standard
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   M ~                                                 
##     X                 0.359    0.047    7.665    0.000
##   Y ~                                                 
##     M                 0.616    0.042   14.795    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .M                 1.097    0.069   15.811    0.000
##    .Y                 1.063    0.067   15.811    0.000
##     X                 1.000    0.063   15.811    0.000

semPaths_default(m1)

versus:

m2_syntax <- '
  M ~ Y
  X ~ M
'

m2 <- sem(m2_syntax, data=data, fixed.x=FALSE, conditional.x=FALSE)
summary(m2)

## lavaan 0.6-3 ended normally after 11 iterations
## 
##   Optimization method                           NLMINB
##   Number of free parameters                          5
## 
##   Number of observations                           500
## 
##   Estimator                                         ML
##   Model Fit Test Statistic                      34.439
##   Degrees of freedom                                 1
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard Errors                             Standard
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   M ~                                                 
##     Y                 0.494    0.033   14.795    0.000
##   X ~                                                 
##     M                 0.293    0.038    7.665    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .M                 0.852    0.054   15.811    0.000
##    .X                 0.895    0.057   15.811    0.000
##     Y                 1.529    0.097   15.811    0.000

semPaths_default(m2)

If you’re unsure of equivalence or nesting, remember NET!

semTools::net(m1, m2)

## 
##         If cell [R, C] is TRUE, the model in row R is nested within column C.
## 
##         If the models also have the same degrees of freedom, they are equivalent.
## 
##         NA indicates the model in column C did not converge when fit to the
##         implied means and covariance matrix from the model in row R.
## 
##         The hidden diagonal is TRUE because any model is equivalent to itself.
##         The upper triangle is hidden because for models with the same degrees
##         of freedom, cell [C, R] == cell [R, C].  For all models with different
##         degrees of freedom, the upper diagonal is all FALSE because models with
##         fewer degrees of freedom (i.e., more parameters) cannot be nested
##         within models with more degrees of freedom (i.e., fewer parameters).
##         
##             m1   m2
## m1 (df = 1)        
## m2 (df = 1) TRUE

Also note that models with the \(c'\) path included are saturated and therefore fit the data equally, and perfectly, well.

But: if we switch what is the predictor/IV \(X\) versus the mediator \(M\), these are not equivalent models.

m3_syntax <- '
Y ~ X
X ~ M
'

m3 <- sem(m3_syntax, data, conditional.x=FALSE, fixed.x=FALSE)
summary(m3, fit.measures=TRUE)

## lavaan 0.6-3 ended normally after 12 iterations
## 
##   Optimization method                           NLMINB
##   Number of free parameters                          5
## 
##   Number of observations                           500
## 
##   Estimator                                         ML
##   Model Fit Test Statistic                     136.919
##   Degrees of freedom                                 1
##   P-value (Chi-square)                           0.000
## 
## Model test baseline model:
## 
##   Minimum Function Test Statistic              271.549
##   Degrees of freedom                                 3
##   P-value                                        0.000
## 
## User model versus baseline model:
## 
##   Comparative Fit Index (CFI)                    0.494
##   Tucker-Lewis Index (TLI)                      -0.518
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -2218.033
##   Loglikelihood unrestricted model (H1)      -2149.573
## 
##   Number of free parameters                          5
##   Akaike (AIC)                                4446.065
##   Bayesian (BIC)                              4467.138
##   Sample-size adjusted Bayesian (BIC)         4451.268
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.521
##   90 Percent Confidence Interval          0.450  0.597
##   P-value RMSEA <= 0.05                          0.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.175
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard Errors                             Standard
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   Y ~                                                 
##     X                 0.473    0.051    9.256    0.000
##   X ~                                                 
##     M                 0.293    0.038    7.665    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .Y                 1.305    0.083   15.811    0.000
##    .X                 0.895    0.057   15.811    0.000
##     M                 1.226    0.078   15.811    0.000

semPaths_default(m3)

This model is not nested in the mediation models above!

#there's a small bug in net for some models where it objects to variable name ordering
semTools::net(m1, m2, m3)

## 
##         If cell [R, C] is TRUE, the model in row R is nested within column C.
## 
##         If the models also have the same degrees of freedom, they are equivalent.
## 
##         NA indicates the model in column C did not converge when fit to the
##         implied means and covariance matrix from the model in row R.
## 
##         The hidden diagonal is TRUE because any model is equivalent to itself.
##         The upper triangle is hidden because for models with the same degrees
##         of freedom, cell [C, R] == cell [R, C].  For all models with different
##         degrees of freedom, the upper diagonal is all FALSE because models with
##         fewer degrees of freedom (i.e., more parameters) cannot be nested
##         within models with more degrees of freedom (i.e., fewer parameters).
##         
##             m1    m2    m3
## m1 (df = 1)               
## m2 (df = 1) TRUE          
## m3 (df = 1) FALSE FALSE

Remember that a good theory of the role of each variable, and their likely temporal precedence, should guide your understanding of plausible models.

2.4.2 No unmodeled third variable confounders

Random assignment to experimental conditions is one of the few ways to handle this problem. Random assignment should nullify the relationship between \(X\) and possible confounders. This does not handle the problem of confounders undermining the relationship between \(M\) and \(Y\), which should be handled by including these as covariates of the \(Y\) regression equation.

More generally, any unmodeled variables that a correlated with \(M\) and/or \(Y\) could undermine the validity of our mediation tests because they could provide a plausible alternative explanation.

2.4.3 Measured variables are proxies for true causal variables

Proxy causal variables. An independent variable/predictor may be highly correlated with a true, but unmeasured, cause of the DV (\(Y\)). For example, imagine that we believe that leaves changing colors and falling off of trees causes decreases in temperature by altering available oxygen (mediator) in the atmosphere. We might observe data consistent with this mediation hypothesis (although oxygen is probably distributed more quickly by wind), but in actuality, changes in temperature reflect the orientation of the local land mass relative to the sun.

Proxy mediator variables. If we measure a broad mediator (e.g., SES), but it is simply highly correlated with something more specific and causal (e.g., food insecurity), we may obtain mediation results spuriously.

Proxy dependent variables. Lastly, one’s DV may also be a proxy and/or there may be a variable upstream in the causal chain:

grViz("digraph regression {
graph [rankdir = LR bgcolor=transparent]

forcelabels=true;

node [shape = box, fontcolor=gray25 color=gray80]
Y1 [style=filled, fontname='Helvetica']; 

node [fontname='Helvetica']
X; M; Y2

edge [color=gray50 style=filled]
X -> M -> Y1 -> Y2
}")

If we measure Y2, but not Y1, then our conclusion about how \(M\) influences the relationship between \(X\) and \(Y2\) may be misleading.

2.4.4 Differential reliability of measurement

Lastly, if measures have different reliabilities, this will have a differential on attenuation of the mediated effects. As Little notes, this is most problematic when a measure is reliable and meaningful in one group, but not another. For example, if we are measuring social media knowledge in young versus old people, there is a chance that the reliability will be poorer in older people who have less exposure to Snapchat, Tumblr, Tindr, etc. This could lead mediation to appear more as between-groups moderation. This reminds us of the importance of measurement invariance and, when possible, the value of including a measurement model of our key constructs to account for unreliability. As we talked about in week 1, in standard regression and path analysis, we assume that predictors are measured with perfect reliability.

2.5.1 Parallel mediators

When there is no (hypothesized) relationship between mediators, but one believes that the total effect of \(X\) on \(Y\) reflects multiple mediators, this is called a parallel mediator model. The mediators are permitted to be correlated to account for their covariance, but one does not include directed paths among mediators since this alters the presumed causal chain.

Recall from the discussion above that when we have parallel mediators, we can partition indirect effects into specific versus total. Moreover, in a parallel mediator model, the scale of the mediators does not influence the magnitude of the specific indirect effects since those are in the units of the \(X - Y\) association. Therefore, once can test whether the specific indirect effects are different in magnitude between competing mediators. This is done using the conventional equality constraint approach with an LRT or Wald test.

2.5.2 Sequential mediators

There can also be chains among mediators in which one mediator is thought to cause another. These increases the number of specific indirect effects and also makes a stronger statement about causal ordering. Nevertheless, sequential mediation can be a powerful test of a hypothesis that involves multiple steps (e.g., in development). Notice also how indirect (mediated) effects can be the product of more than two coefficients.

3.1.1 Interpreting effects of individual predictors in a moderation model

Darlington and Hayes (2017) warn that the coefficients for two predictors in a regression model that includes their interaction should not typically be thought of as ‘main effects’ (unless \(X\) and \(W\) are categorical). Rather, \(b1\) and \(b2\) represent something closer to simple effects in this model. When the two-way \(XW\) interaction is included, \(b1\) is the relationship between \(X\) and \(Y\) when \(W = 0\). Likewise, \(b2\) is the effect of \(W\) on \(Y\) when \(X = 0\). One challenge that inheres in this representation is that \(X\) or \(W\) may not have meaningful or at least interpretable zero points.

3.1.2 The importance of centering in moderation analysis

Mean-centering \(X\) and \(W\) before computing their interaction is almost always the appropriate decision for two reasons. First, this eliminates what West calls non-essential collinearity. In other words, the interaction can be collinear with either predictor simply because of the scaling of the variables, which is largely eliminated by centering (esp. when the variables are bivariate normal). This collinearity affects the individual predictor effects, \(b1\) and \(b2\), but not the interaction, \(b3\). Second, centering improves the interpretability of \(b1\) and \(b2\) by shifting their interpretation: under mean centering, \(b1\) reflects the relationship between \(X\) and \(Y\) at the mean of \(W\).

3.2.1 Binary moderator

When the moderator is binary (e.g., two groups), the \(b3\) coefficient (i.e., the interaction) quantifies the difference between the reference group (0) and treatment group (1).

3.2.2 Continuous moderator

When moderators are continuous, the most common approach is to plot the relationship between \(X\) and \(Y\) at different levels of \(W\). One effective convention is to plot the predicted \(X - Y\) association at the mean of \(W\) +/- 1 SD (i.e., three plots).

Interaction plot

states <- as.data.frame(state.x77)
states$HSGrad <- states$`HS Grad`
fit <- lm(Income ~ HSGrad + Murder*Illiteracy, data = states)

interact_plot(model = fit, pred = Murder, modx = Illiteracy)

3.2.3 Understanding the regions of significance in a moderated effect

Visually examining the association of \(X\) and \(Y\) at different levels of the moderator \(W\) helps to give an intuition of the values of \(W\) where moderation is likely to exist (i.e., the \(X-Y\) slope is non-zero). But to test this formally, one needs to estimate the \(X-Y\) relationship at different candidate values of \(W\) and test the significance of the relationship (i.e., at what level of \(W\) are \(X\) and \(Y\) significantly related?).

One way to handle this is the ‘pick a point’ method in which one computes the predicted \(X-Y\) relationship at a given value of \(W\) and tests the significance of the slope (a kind of simple slopes analysis). This can be accomplished by centering \(W\) at the candidate value. Under this approach, \(b1\) provides the simple slope, standard error, and appropriate p-value for the relationship of \(X\) and \(Y\) at the candidate value of \(W\).

The generalization of this approach is called the Johnson-Neyman technique, which identifies the values of \(W\) at which the relationship between \(X\) and \(Y\) is significant (Bauer & Curran, 2005). Although I’m not aware of a lavaan-friendly package for this, you might check out the interactions package, which implements this approach in standard regression models. https://cran.r-project.org/web/packages/interactions/index.html.

Johnson-Neyman plot

johnson_neyman(model = fit, pred = Murder, modx = Illiteracy)

## JOHNSON-NEYMAN INTERVAL 
## 
## When Illiteracy is OUTSIDE the interval [0.80, 2.67], the slope of
## Murder is p < .05.
## 
## Note: The range of observed values of Illiteracy is [0.50, 2.80]