Path Analysis

This week we are temporarily putting aside the idea of latent variables, and looking in more depth at the framework for fitting and testing models as a system of variables connected by arrows (or “paths”).

Next week, we bring this technique together with what we have learned about latent variable modelling (i.e., the idea of a latent “factor”), to enable us to build some really quite sophisticated models.

Motivating Path Analysis

Over the course of USMR and the first block of this course, you have hopefully become pretty comfortable with the regression world, and can see how it is extended to lots of different types of outcome and data structures. So far in MSMR our discussion of regressions have all focused on how this can be extended to multiple levels. This has brought so many more potential study designs that we can now consider modelling - pretty much any study where we are interested in explaining some outcome variable, and where we have sampled clusters of observations (or clusters of clusters of clusters of … etc.).

But we are still restricted to thinking, similar to how we thought in USMR, about one single outcome variable. In fact, if we think about the structure of the fixed effects part of a model (i.e., the bit we’re specifically interested in), then we’re still limited to thinking of the world in terms of “this is my outcome variable, everything else predicts it”.

Regression as a path diagram

1. Imagine writing the names of all your variables on a whiteboard
2. Specify which one is your dependent (or “outcome” or “response”) variable.
3. Sit back and relax, you’re done!

In terms of a theoretical model of the world, there’s not really much to it. We have few choices in the model we construct beyond specifying which is our outcome variable.
We can visualise our multiple regression model like this:

Figure 1: In multiple regression, we decide which variable is our outcome variable, and then everything else is done for us

Of course, there are a few other things that are included (an intercept term, the residual error, and the fact that our predictors can be correlated with one another), but the idea remains pretty much the same:

Figure 2: Multiple regression with intercept, error, predictor covariances

Theories and Models

What if I my theoretical model of the world doesn’t fit this structure?

Let’s suppose I have 5 variables: Age, Parental Income, Income, Autonomy, and Job Satisfaction. I draw them up on my whiteboard:

Figure 3: My variables

My theoretical understanding of how these things fit together leads me to link my variables to end up with something like that in Figure 4.

Figure 4: My theory about my system of variables

In this diagram, a persons income is influenced by their age, their parental income, and their level of autonomy, and in turn their income predicts their job satisfaction. Job satisfaction is also predicted by a persons age directly, and by their level of autonomy, which is also predicted by age. It’s complicated to look at, but in isolation each bit of this makes theoretical sense.

Take each arrow in turn and think about what it represents:

If we think about trying to fit this “model” with the tools that we have, then we might end up wanting to fit three separate regression models, which between them specify all the different arrows in the diagram:

\[ \begin{align} \textrm{Job Satisfaction} & = \beta_0 + \beta_1(\textrm{Age}) + \beta_2(\textrm{Autonomy}) + \beta_3(\textrm{Income}) + \varepsilon \\ \textrm{Income} & = \beta_0 + \beta_1(\textrm{Age}) + \beta_2(\textrm{Autonomy}) + \beta_3(\textrm{Parental Income}) + \varepsilon \\ \textrm{Autonomy} & = \beta_0 + \beta_1(\textrm{Age}) + \varepsilon \\ \end{align} \]

In this case, our theoretical model involves having multiple endogenous (think “outcome”) variables. So what do we do if we want to talk about how well the entire model (Figure 4) fits the data we observed? This is where path analysis techniques come in handy.

New terminology!

• Exogenous variables are a bit like what we have been describing with words like “independent variable” or “predictor”. In a SEM diagram, they have no paths coming from other variables in the system, but have paths going to other variables.
• Endogenous variables are more like the “outcome”/“dependent”/“response” variables we are used to. They have some path coming from another variable in the system (and may also - but not necessarily - have paths going out from them).

Introducing Path Analysis

It might help to think of the starting point of path analysis as drawing the variables on a whiteboard and drawing arrows to reflect your theory about the system of variables that you observed (much like in Figure 3 and 4 above).

As it happens, we have already seen the conventions for how to depict variables and parameters in these type of diagrams in last week’s lab: by using rectangles (observed variables), ovals (latent variables), single headed arrows (regression paths) and double headed arrows (covariances), we can draw various model structures (as mentioned at the top of this page, we are temporarily putting aside the latent variables today).

We could try to fit some of these more complex models by just using lm() many times over and fitting various regression models. Path analysis is a bit like this, but all in one. It involves fitting a set of regression equations simultaneously. One obvious benefit of this is that it can allow us to talk about “model fit” (in the way we discussed in last week’s lab) in relation to our entire theory, rather than to individual regressions that make up only part of the theoretical model.

Refresher on path diagrams

• Observed variables are represented by squares or rectangles. These are the named variables of interest which exist in our dataset - i.e. the ones which we have measured directly.
  
• Latent variables are represented as ovals/ellipses or circles.
  
• Covariances are represented by double-headed arrows. In many diagrams these are curved.
  
• Regressions are shown by single headed arrows (e.g., an arrow from \(x\) to \(y\) for the path \(y \sim x\)). Factor loadings are also regression paths.
  Recall that specifying a factor structure is simply to say that some measured variables \(y_i\) are each regressed onto some unmeasured factor(s) - \(y = \lambda \cdot F + u\) looks an awful lot like \(y = \beta \cdot x + \epsilon\)!!).

Mediation

Another benefit is that we are no longer limited to studying only simple relationships between x and y, we can now study how x might change z which in turn might change y.

As an example, let’s imagine we are interested in peoples’ intention to get vaccinated, and we observe the following variables:

• Intention to vaccinate (scored on a range of 0-100)
  
• Health Locus of Control (HLC) score (average score on a set of items relating to perceived control over ones own health)
  
• Religiosity (average score on a set of items relating to an individual’s religiosity).

We are assuming here that we do not have the individual items, but only the scale scores (if we had the individual items we might be inclined to model religiosity and HLC as latent variables!).
If we draw out our variables, and think about this in the form of a standard regression model with “Intention to vaccinate” as our outcome variable, then all the lines are filled in for us (see Figure 5)

Figure 5: Multiple regression: choose your outcome, sit back and relax

But what if our theory suggests that some other model might be of more relevance? For instance, what if we believe that participants’ religiosity has an effect on their Health Locus of Control score, which in turn affects the intention to vaccinate (see Figure 6)?

In this case, the HLC variable is thought of as a mediator, because is mediates the effect of religiosity on intention to vaccinate. We are specifying presence of a distinct type of effect: direct and indirect.

Direct vs Indirect

In path diagrams:

• Direct effect = one single-headed arrow between the two variables concerned
  
• Indirect effect = An effect transmitted via some other variables

If we have a variable \(X\) that we take to ‘cause’ variable \(Y\), then our path diagram will look like so: In this diagram, path \(c\) is the total effect. This is the unmediated effect of \(X\) on \(Y\).

However, while the effect of \(X\) on \(Y\) could in part be explained by the process of being mediated by some variable \(M\), the variable \(X\) could still affect \(Y\) directly.
Our mediating model is shown below:

In this case, path \(c'\) is the direct effect, and paths \(a\) and \(b\) make up the indirect effect.

You will find in some areas people talk about the ideas of “complete” vs “partial” mediation. “Complete mediation” is when \(X\) no longer affects \(Y\) after \(M\) has been controlled (so path \(c'\) is not significantly different from zero), and “partial mediation” is when the path from \(X\) to \(Y\) is reduced in magnitude when the mediator \(M\) is introduced, but still different from zero.

Figure 6: Mediation as a path model. If you’re interested, you can find the inspiration for this data from the paper here. I haven’t properly read it though!

Now that we’ve seen how path analysis works, we can use that same logic to investigate models which have quite different structures, such as those including mediating variables. So if we can’t fit our theoretical model into a regression framework, let’s just fit it into a framework which is lots of regressions smushed together!

Luckily, we can just get the lavaan package to do all of this for us. So let’s look at fitting the model in Figure 6.

First we read in our data:

vax <- read_csv("https://uoepsy.github.io/data/vaxdat.csv")
summary(vax)

##   religiosity        hlc         intention   
##  Min.   :-1.0   Min.   :0.40   Min.   :39.0  
##  1st Qu.: 1.8   1st Qu.:2.00   1st Qu.:59.0  
##  Median : 2.4   Median :3.00   Median :64.0  
##  Mean   : 2.4   Mean   :2.99   Mean   :65.1  
##  3rd Qu.: 3.0   3rd Qu.:3.60   3rd Qu.:74.0  
##  Max.   : 4.6   Max.   :5.80   Max.   :88.0

Then we specify the relevant paths:

med_model <- " 
    intention ~ religiosity
    intention ~ hlc
    hlc ~ religiosity
"

If we fit this model as it is, we won’t actually be testing the indirect effect, we will simply be fitting a couple of regressions.

To do that, we need to explicitly define the indirect effect in the model, by first creating a label for each of its sub-component paths, and then defining the indirect effect itself as the product of these (click here for a lovely pdf explainer from Aja on why the indirect effect is calculated as the product).

In lavaan, we use a new operator, :=, to create this estimate.

med_model <- " 
    intention ~ religiosity
    intention ~ b*hlc
    hlc ~ a*religiosity
    
    indirect := a*b
"

:=

This operator ‘defines’ new parameters which take on values that are an arbitrary function of the original model parameters. The function, however, must be specified in terms of the parameter labels that are explicitly mentioned in the model syntax.

(the lavaan project)

Note. The labels we use are completely up to us. This would be equivalent:

med_model <- " 
    intention ~ religiosity
    intention ~ peppapig * hlc
    hlc ~ kermit * religiosity
    
    indirect := kermit * peppapig
"

Finally, we estimate our model.
It is common to estimate the indirect effect using bootstrapping (a method of resampling the data with replacement, thousands of times, in order to empirically generate a sampling distribution). We can do this easily in lavaan, using the sem() function:

mm1.est <- sem(med_model, data=vax, se = "bootstrap") 

And we can get out estimates for all our parameters, including the one we created called “indirect”.

summary(mm1.est, ci = TRUE)

## lavaan 0.6.15 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         5
## 
##   Number of observations                           100
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                            Bootstrap
##   Number of requested bootstrap draws             1000
##   Number of successful bootstrap draws            1000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##   intention ~                                                           
##     religiosty        0.270    1.022    0.265    0.791   -1.583    2.382
##     hlc        (b)    5.971    0.933    6.397    0.000    4.044    7.747
##   hlc ~                                                                 
##     religiosty (a)    0.508    0.086    5.916    0.000    0.348    0.680
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##    .intention        62.090    7.926    7.834    0.000   45.650   77.632
##    .hlc               0.753    0.100    7.572    0.000    0.552    0.927
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##     indirect          3.033    0.655    4.633    0.000    1.825    4.344

Mediation Exercises

This week’s lab focuses on the technique of path analysis using the same context as previous weeks: conduct problems in adolescence. In this week’s example, a researcher has collected data on n=557 adolescents and would like to know whether there are associations between conduct problems (both aggressive and non-aggressive) and academic performance and whether the relations are mediated by the quality of relationships with teachers.

The data is available at https://uoepsy.github.io/data/cp_teachacad.csv

Question A1

First, read in the dataset.

Solution

cp_teach<-read_csv("https://uoepsy.github.io/data/cp_teachacad.csv")
summary(cp_teach)

##        ID           Acad           Teach_r         Non_agg           Agg       
##  Min.   :  1   Min.   :-3.042   Min.   :-3.68   Min.   :-3.28   Min.   :-3.27  
##  1st Qu.:140   1st Qu.:-0.736   1st Qu.:-0.89   1st Qu.:-0.69   1st Qu.:-0.73  
##  Median :279   Median :-0.026   Median : 0.01   Median :-0.09   Median :-0.02  
##  Mean   :279   Mean   :-0.061   Mean   :-0.09   Mean   :-0.06   Mean   :-0.02  
##  3rd Qu.:418   3rd Qu.: 0.606   3rd Qu.: 0.62   3rd Qu.: 0.57   3rd Qu.: 0.68  
##  Max.   :557   Max.   : 3.057   Max.   : 3.52   Max.   : 3.44   Max.   : 3.31

Question A2

Use the sem() function in lavaan to specify and estimate a straightforward linear regression model to test whether aggressive and non-aggressive conduct problems significantly predict academic performance.

How do your results compare to those you obtain using the lm() function?

Solution

# we can fit the model in lavaan as follows:
# first we specify the model using lavaan syntax
sr_lavaan<-'Acad ~ Non_agg+Agg'

# next we can estimate the model using the sem() function
sr_lavaan.est<-sem(sr_lavaan, data=cp_teach)

# we can inspect the results using the summary() function
summary(sr_lavaan.est)

## lavaan 0.6.15 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         3
## 
##   Number of observations                           557
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   Acad ~                                              
##     Non_agg           0.182    0.057    3.178    0.001
##     Agg               0.318    0.057    5.599    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .Acad              0.943    0.057   16.688    0.000

# the same model can be fit using lm():
sr_lm<-lm(Acad~Non_agg+Agg, data=cp_teach)
summary(sr_lm)

## 
## Call:
## lm(formula = Acad ~ Non_agg + Agg, data = cp_teach)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8962 -0.5957  0.0073  0.6219  3.1325 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.0420     0.0414   -1.02   0.3097    
## Non_agg       0.1824     0.0576    3.17   0.0016 ** 
## Agg           0.3181     0.0570    5.58  3.7e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.974 on 554 degrees of freedom
## Multiple R-squared:  0.194,  Adjusted R-squared:  0.191 
## F-statistic: 66.5 on 2 and 554 DF,  p-value: <2e-16

We can see that both non-aggressive and aggressive conduct problems significantly predict academic performance. We can also see that we get the same results when we use the sem() function as we do when we use the lm() function. Lavaan will give essentially the same results as lm() for simple and multiple regression problems. However, if we have multiple outcome variables in our model it is advantageous to do this using path mediation model with lavaan. This allows us to include all the regressions in a single model.

Question A3

Now specify a model in which non-aggressive conduct problems have both a direct and indirect effect (via teacher relationships) on academic performance.

Sketch the path diagram using the website suggested in the lecture, https://www.diagrams.net/. To export the diagram to PNG, click File -> Export As -> PNG.

Solution

model1<-'
    #we regress academic performance on non-aggressive conduct problems (the direct effect)
    Acad ~ Non_agg
    
    #we regress academic peformance on teacher relationship quality
    Acad ~ Teach_r
    
    #we regress teacher relationship quality on non-aggressive conduct problems
    Teach_r ~ Non_agg 
'

Question A4

Now define the indirect effect in order to test the hypothesis that non-aggressive conduct problems have both a direct and an indirect effect (via teacher relationships) on academic performance.

Fit the model and examine the 95% CI.

Solution

#model specification
model1 <- '
    Acad ~ Non_agg
    #we label the two parameters that comprise the indirect effect b and a
    Acad ~ b*Teach_r    
    Teach_r ~ a*Non_agg  
    
    # the indirect effect is the product of a and b. We create a new parameter (ind) to estimate the indirect effect
    ind := a*b   
'

#model estimation
model1.est <- sem(model1, data=cp_teach, se='bootstrap') 

# we request bootstrapped standard errors to assess the significance of the indirect effect
summary(model1.est, ci=T)

## lavaan 0.6.15 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         5
## 
##   Number of observations                           557
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                            Bootstrap
##   Number of requested bootstrap draws             1000
##   Number of successful bootstrap draws            1000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##   Acad ~                                                                
##     Non_agg           0.158    0.057    2.793    0.005    0.042    0.271
##     Teach_r    (b)    0.328    0.047    6.994    0.000    0.234    0.418
##   Teach_r ~                                                             
##     Non_agg    (a)    0.769    0.034   22.795    0.000    0.700    0.838
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##    .Acad              0.919    0.057   16.079    0.000    0.804    1.030
##    .Teach_r           0.713    0.042   17.088    0.000    0.632    0.799
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##     ind               0.252    0.037    6.782    0.000    0.179    0.326

Note

The confidence intervals computed using bootstrapping may be slightly different for you. This is normal, as bootstrap is based on random sampling!

We can see that the 95% bootstrapped confidence interval for the indirect effect of non-aggressive conduct problems on academic performance (‘ind’) does not include zero. We can conclude that the indirect effect is significant at \(\alpha = .05\). The direct effect is also statistically significant.

Question A5

Specify a new parameter which is the total (direct+indirect) effect of non-aggressive conduct problems on academic performance

Solution

We can create a new parameter that is the sum of the direct and indirect effect to evaluate the total effect of non-aggressive conduct problems on academic performance.

#model specification
model1 <- '
    # we now also label the direct effect of non-aggressive conduct problems on academic performance
    Acad ~ c*Non_agg    
    Acad ~ b*Teach_r    
    Teach_r ~ a*Non_agg  
    
    ind := a*b
    
    #the total effect is the indirect effect plus the direct effect
    total := c + a*b 
'

#model estimation
model1.est <- sem(model1, data=cp_teach, se='bootstrap') 

# we request bootstrapped standard errors to assess the significance of the indirect effect
summary(model1.est, ci=T)

## lavaan 0.6.15 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         5
## 
##   Number of observations                           557
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                            Bootstrap
##   Number of requested bootstrap draws             1000
##   Number of successful bootstrap draws            1000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##   Acad ~                                                                
##     Non_agg    (c)    0.158    0.057    2.759    0.006    0.053    0.277
##     Teach_r    (b)    0.328    0.049    6.630    0.000    0.230    0.425
##   Teach_r ~                                                             
##     Non_agg    (a)    0.769    0.032   23.797    0.000    0.704    0.835
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##    .Acad              0.919    0.056   16.423    0.000    0.805    1.029
##    .Teach_r           0.713    0.041   17.216    0.000    0.633    0.790
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##     ind               0.252    0.039    6.453    0.000    0.176    0.331
##     total             0.410    0.043    9.498    0.000    0.327    0.497

Question A6

Now visualise the estimated model and its parameters using the semPaths() function from the semPlot package.

Solution

#to include the parameter estimates we set what='est'
semPaths(model1.est, what='est')

#to change label size, add edge.label.cex = 1.5
#you can adjust 1.5 to other values

A more complex model

Question B1

Using the website website suggested in the lecture, https://www.diagrams.net/, sketch the path diagram for a model in which both aggressive and non-aggressive conduct problems have both direct and indirect effects (via teacher relationships) on academic performance. To export the diagram to PNG, click File -> Export As -> PNG.

Now specify the model in R, taking care to also define the parameters for the indirect effects.

Solution

There are many ways in which you could position the variables in your diagram, so the two path diagrams below are equivalent:

Option A

Option B

We now have two predictors, one mediator and one outcome (and two indirect effects, one for each predictor). We can represent this in two lines: one where we specify academic performance as the outcome variable and one where we specify teacher relationships (the mediator) as the outcome variable.

model2 <- '
    Acad ~ Agg + Non_agg + b*Teach_r
    Teach_r ~ a_A*Agg + a_NA*Non_agg
   
    ind_A  := a_A*b  #indirect effect for aggressive conduct problems
    ind_NA := a_NA*b #indirect effect for non-aggressive conduct problems
'

Question B2

Now estimate the model and test the significance of the indirect effects

Solution

model2.est<-sem(model2,  data=cp_teach,se='bootstrap') 
summary(model2.est, ci=T)

## lavaan 0.6.15 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         7
## 
##   Number of observations                           557
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                            Bootstrap
##   Number of requested bootstrap draws             1000
##   Number of successful bootstrap draws            1000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##   Acad ~                                                                
##     Agg               0.171    0.060    2.852    0.004    0.058    0.293
##     Non_agg           0.091    0.061    1.490    0.136   -0.033    0.203
##     Teach_r    (b)    0.256    0.054    4.767    0.000    0.146    0.358
##   Teach_r ~                                                             
##     Agg      (a_A)    0.574    0.044   12.997    0.000    0.488    0.659
##     Non_agg (a_NA)    0.358    0.041    8.648    0.000    0.275    0.440
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##    .Acad              0.908    0.058   15.686    0.000    0.791    1.019
##    .Teach_r           0.540    0.028   19.077    0.000    0.481    0.591
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##     ind_A             0.147    0.032    4.556    0.000    0.081    0.210
##     ind_NA            0.092    0.023    4.051    0.000    0.048    0.137

We can see that the 95% confidence intervals for both indirect effects do not include zero, therefore, we can conclude that they are significant at the \(\alpha = .05\) level.

Question B3

Write a brief paragraph reporting on the results of the model estimates in Question B2. Include a Figure or Table to display the parameter estimates.

Solution

A path mediation model was used to test the direct and indirect effects (via teacher relationship quality) of aggressive and non-aggressive conduct problems on academic performance. In the model, academic performance was regressed on teacher relationship quality, non-aggressive conduct problems and aggressive conduct problems while teacher relationship quality (the mediator) was regressed on aggressive and non-aggressive conduct problems. The indirect effects were tested using the product of the coefficient for the regression of outcome on mediator and the coefficient for the regression of mediator on predictor. The statistical significance of the indirect effects were evaluated using bootstrapped 95% confidence intervals with 1000 bootstrap samples.

Unstandardised parameter estimates are provided in Figure 7. Solid lines indicate that a parameter is significant at the 5% significance level, while dashed lines represent non-significant paths.The indirect effects of both non-aggressive (\(b = 0.09\), 95% CI=0.05-0.14) and aggressive (\(b = 0.15\), 95% CI=0.08-0.22) conduct problems on academic performance were statistically significant.

Figure 7: Effect of conduct problems on academic performance mediated by quality of teacher relationship.

Optional: Mediation the more manual way: back to lm()

Following Baron & Kenny 1986, we can conduct mediation analysis by using three separate regression models.

1. \(y \sim x\)
2. \(m \sim x\)
3. \(y \sim x + m\)

Step 1. Determine the presence of y ~ x:
if x predicts y, then there is possibility to detect mediation

vax <- read_csv("https://uoepsy.github.io/data/vaxdat.csv")

mod1 <- lm(intention ~ religiosity, data = vax)
summary(mod1)$coefficients

##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)     57.2      2.422   23.61 3.17e-42
## religiosity      3.3      0.929    3.55 5.84e-04

Step 2. Determine the presence of m ~ x: if x predicts m, then there is possibility to detect mediation

mod2 <- lm(hlc ~ religiosity, data = vax)
summary(mod2)$coefficients

##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    1.775     0.2229    7.96 3.04e-12
## religiosity    0.508     0.0855    5.94 4.36e-08

Step 3. Examine the effect of y ~ x + m:
If the x no longer predicts y after partialling out effects due to m, then there is full mediation. If the effect of x on y is smaller, then there is partial mediation.

mod3 <- lm(intention ~ religiosity + hlc, data = vax)
summary(mod3)$coefficients

##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    46.58      2.610  17.846 2.10e-32
## religiosity     0.27      0.910   0.297 7.67e-01
## hlc             5.97      0.922   6.478 3.85e-09

Step 4. Test for the mediation.
There are various ways to do this, but the simplest is probably:

library(mediation)
summary(mediate(mod2, mod3, treat='religiosity', mediator='hlc', boot=TRUE, sims=500))

## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value    
## ACME              3.033        1.808         4.37  <2e-16 ***
## ADE               0.270       -1.771         2.35     0.8    
## Total Effect      3.303        1.351         5.30  <2e-16 ***
## Prop. Mediated    0.918        0.505         2.07  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 100 
## 
## 
## Simulations: 500

ACME: Average Causal Mediation Effects
ADE: Average Direct Effects
Total Effect: sum of the mediation (indirect) effect and the direct effect.