Multilevel path anaylysis in lavaan: Testing for mediation at both Level 1 and Level 2

Multilevel-path-anaylysis-in-lavaan-Testing-for-mediation-at-both-Level-1-and-Level-2-web.knit

Overview

For this example, we are testing a model that incorporates mediation at both Level 1 and Level 2. We will include compositional variables (i.g., group means of Level 1) variables as predictors of the mediator at Level 2 and (indirectly) the outcome variable (also at Level 2). This discussion picks up from previous demonstrations on multilevel path analysis with lavaan found at https://mikesquanthub.blogspot.com/2024/04/multilevel-path-analysis-in-lavaan.html and https://mikesquanthub.blogspot.com/2024/04/multilevel-path-model-using-lavaan.html.

The example is adapted from a study by Xie et al. (2015), “Linking Colleague Support to Employees’ Promotive Voice: A Moderated Mediation Model” published in PLOS One. Although we will not be testing the authors’ model, we will rely on the dataset they published in the Supplementary Materials on the journal site for their article (https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0132123).

I have constructed a smaller dataset from the original data and made some changes in naming, etc. You can download a copy of the .csv file containing the data here. The .csv file is named prvoiceV5. Once you have downloaded the data to your computer, you will need to import it into RStudio for use.

For our analysis, we are modeling data based on \(n_1\)=162 employees (level 1 units) nested within \(n_2\)=51 work teams (level 2 units). At the employee level, we are predicting promotive voice as a function of felt obligation, psychological safety, and team identification. Felt obligation is treated as a mediator of the effects of psychological safety and team identification on promotive voice.

At level 2, we are modeling the effects of the compositional effects of psychological safety and team identification on felt obligation and promotive voice at the team-level.

Here is a figure of the model we will be specifying and testing using lavaan.



Let’s take a look at the names of the variables in our dataset.

names(prvoiceV5)
##  [1] "X"                           "teamcode"                   
##  [3] "colleaguesupport"            "feltobl"                    
##  [5] "gender"                      "age"                        
##  [7] "eduation"                    "tenure"                     
##  [9] "teamsize"                    "teamID"                     
## [11] "subgroupformationindividual" "psysafe"                    
## [13] "prvoice"

The variables we will be using when running our analysis include:

  1. prvoice (level 1 outcome; promotive voice)
  2. feltobl (level 1 mediator; felt obligation)
  3. teamID (level 1 antecedent; team identification)
  4. psysafe (level 1 antecedent; psychological safety)

All variables are treated as continuous in the analysis. [Keep in mind too that lavaan currently does not allow for categorical endogenous variables]


Before carrying out our main analyses, we will want to generate cluster means (i.e., team-level means) for our exogenous variables to serve as predictors at Level 2.

library(mitml)
## *** This is beta software. Please report any bugs!
## *** See the NEWS file for recent changes.
prvoiceV5$teamID_M <- clusterMeans(prvoiceV5$teamID, cluster = prvoiceV5$teamcode)
prvoiceV5$psysafe_M <- clusterMeans(prvoiceV5$psysafe, cluster = prvoiceV5$teamcode)

names(prvoiceV5) 
##  [1] "X"                           "teamcode"                   
##  [3] "colleaguesupport"            "feltobl"                    
##  [5] "gender"                      "age"                        
##  [7] "eduation"                    "tenure"                     
##  [9] "teamsize"                    "teamID"                     
## [11] "subgroupformationindividual" "psysafe"                    
## [13] "prvoice"                     "teamID_M"                   
## [15] "psysafe_M"
head(prvoiceV5)
##   X teamcode colleaguesupport feltobl gender age eduation tenure teamsize
## 1 1        2             4.50    5.00      1   2        4     28       10
## 2 2        2             4.50    4.20      1   2        4     12       10
## 3 3        2             3.50    4.45      1   3        3     23       10
## 4 4        4             3.75    4.20      0   2        2      3        4
## 5 5        4             5.00    4.20      0   2        3     18        4
## 6 6        4             2.25    3.00      1   2        3     18        4
##   teamID subgroupformationindividual psysafe prvoice teamID_M psysafe_M
## 1   6.50                        2.00    6.57     4.2 6.553333  6.046667
## 2   6.33                        3.50    6.00     3.8 6.553333  6.046667
## 3   6.83                        3.00    5.57     5.4 6.553333  6.046667
## 4   5.50                        5.50    5.71     4.0 5.113333  4.236667
## 5   4.67                        4.50    4.14     5.2 5.113333  4.236667
## 6   5.17                        2.75    2.86     6.0 5.113333  4.236667

We see here that both compositional variables have been added to the data frame.

The main analyses

library(lavaan)
model1a <- '

level: 1  # specification of level 1 model

prvoice ~ feltobl  
feltobl ~ psysafe + teamID

level: 2  # specification of level 2 model

prvoice ~ feltobl   
feltobl ~ psysafe_M + teamID_M 

'


The following code will fit the model and generate our output. Make sure you include the cluster argument (see below) and set it equal to the name of the clustering variable (enclosed in quotes). [The output will contain warning messages about lack of variance within clusters for certain variables. This is simply a warning. But I do not reproduce it here.]

fit1 <- sem(model1a, data = prvoiceV5, cluster = "teamcode")

summary(fit1, fit.measures=T, rsquare=T)
## lavaan 0.6.13 ended normally after 51 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        12
## 
##   Number of observations                           162
##   Number of clusters [teamcode]                     51
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 2.502
##   Degrees of freedom                                 4
##   P-value (Chi-square)                           0.644
## 
## Model Test Baseline Model:
## 
##   Test statistic                               127.043
##   Degrees of freedom                                10
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000
##   Tucker-Lewis Index (TLI)                       1.032
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)               -343.478
##   Loglikelihood unrestricted model (H1)       -342.227
##                                                       
##   Akaike (AIC)                                 710.956
##   Bayesian (BIC)                               748.008
##   Sample-size adjusted Bayesian (SABIC)        710.018
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.095
##   P-value H_0: RMSEA <= 0.050                    0.788
##   P-value H_0: RMSEA >= 0.080                    0.091
## 
## Standardized Root Mean Square Residual (corr metric):
## 
##   SRMR (within covariance matrix)                0.026
##   SRMR (between covariance matrix)               0.075
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Observed
##   Observed information based on                Hessian
## 
## 
## Level 1 [within]:
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   prvoice ~                                           
##     feltobl           0.430    0.109    3.926    0.000
##   feltobl ~                                           
##     psysafe           0.365    0.051    7.179    0.000
##     teamID            0.285    0.051    5.559    0.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .prvoice           0.000                           
##    .feltobl           0.000                           
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .prvoice           0.625    0.084    7.434    0.000
##    .feltobl           0.229    0.029    7.933    0.000
## 
## R-Square:
##                    Estimate
##     prvoice           0.123
##     feltobl           0.518
## 
## 
## Level 2 [teamcode]:
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   prvoice ~                                           
##     feltobl           3.945    2.372    1.663    0.096
##   feltobl ~                                           
##     psysafe_M        -0.041    0.071   -0.575    0.565
##     teamID_M         -0.077    0.092   -0.838    0.402
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .prvoice           0.067    2.699    0.025    0.980
##    .feltobl           1.457    0.452    3.223    0.001
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .prvoice           0.328    0.291    1.130    0.258
##    .feltobl           0.024    0.017    1.423    0.155
## 
## R-Square:
##                    Estimate
##     prvoice           0.569
##     feltobl           0.131

The chi-square goodness of fit test was not statistically significant, \(\chi^2\)(4)=2.502, p=.644 (lead to a failure to reject the exact fit hypothesis). The CFI=1 and TLI = 1.032 were above .95, indicating a very good fitting model. Similarly, the RMSEA of approximately 0 indicates a very good fitting model. The SRMR within (.026) was below .05 suggesting good fit at Level 1. The SRMR between (.075) was rather high indicating a fit that was not particularly strong.

At Level 1, team identification (b=.285, p<.001) and psychological safety (b=.365, p<.001) were both significant positive predictors of felt obligation. Felt obligation, in turn, was a positive and significant (b=.430, p<.001) predictor of promotive voice.

At Level 2, neither the compositional effects of team identification or psychological safety predicted felt obligation (at the team level). Moreover, felt obligation did not significantly predict variation between teams in promotive voice.

If you would like to report on the intraclass correlations for the endogenous level 1 variables use the lavInspect() function.

lavInspect(fit1, "icc")
## prvoice feltobl psysafe  teamID 
##   0.523   0.051   0.000   0.000

We can test for mediation at Level 1 and Level 2 using lavaan. To do this, you will need to (a) assign labels to the paths used in the construction of the indirect effects and then (b) create estimands for the indirect effects. The figure below shows you the paths involved in the construction of the indirect effect estimands.




Specification of our model…

model1a <- '

level: 1  # specification of level 1 model

prvoice ~ l1b*feltobl  
feltobl ~ l1a1*psysafe + l1a2*teamID

# Adding estimands to capture indirect effects

IE1 := l1a1*l1b #psysafe->feltobl->prvoice
IE2 := l1a2*l1b #teamID->feltobl->prvoice

level: 2  # specification of level 2 model

prvoice ~ l2b*feltobl  
feltobl ~ l2a1*psysafe_M + l2a2*teamID_M 

# Adding estimands to capture indirect effects

IE3 := l2a1*l2b #psysafe_M ->feltobl->prvoice
IE4 := l2a2*l2b #teamID_M ->feltobl->prvoice
'

fit1 <- sem(model1a, data = prvoiceV5, cluster = "teamcode")

summary(fit1, fit.measures=T, rsquare=T)
## lavaan 0.6.13 ended normally after 51 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        12
## 
##   Number of observations                           162
##   Number of clusters [teamcode]                     51
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 2.502
##   Degrees of freedom                                 4
##   P-value (Chi-square)                           0.644
## 
## Model Test Baseline Model:
## 
##   Test statistic                               127.043
##   Degrees of freedom                                10
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000
##   Tucker-Lewis Index (TLI)                       1.032
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)               -343.478
##   Loglikelihood unrestricted model (H1)       -342.227
##                                                       
##   Akaike (AIC)                                 710.956
##   Bayesian (BIC)                               748.008
##   Sample-size adjusted Bayesian (SABIC)        710.018
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.095
##   P-value H_0: RMSEA <= 0.050                    0.788
##   P-value H_0: RMSEA >= 0.080                    0.091
## 
## Standardized Root Mean Square Residual (corr metric):
## 
##   SRMR (within covariance matrix)                0.026
##   SRMR (between covariance matrix)               0.075
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Observed
##   Observed information based on                Hessian
## 
## 
## Level 1 [within]:
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   prvoice ~                                           
##     feltobl  (l1b)    0.430    0.109    3.926    0.000
##   feltobl ~                                           
##     psysafe (l1a1)    0.365    0.051    7.179    0.000
##     teamID  (l1a2)    0.285    0.051    5.559    0.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .prvoice           0.000                           
##    .feltobl           0.000                           
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .prvoice           0.625    0.084    7.434    0.000
##    .feltobl           0.229    0.029    7.933    0.000
## 
## R-Square:
##                    Estimate
##     prvoice           0.123
##     feltobl           0.518
## 
## 
## Level 2 [teamcode]:
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   prvoice ~                                           
##     feltobl  (l2b)    3.945    2.372    1.663    0.096
##   feltobl ~                                           
##     psysf_M (l2a1)   -0.041    0.071   -0.575    0.565
##     temID_M (l2a2)   -0.077    0.092   -0.838    0.402
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .prvoice           0.067    2.699    0.025    0.980
##    .feltobl           1.457    0.452    3.223    0.001
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .prvoice           0.328    0.291    1.130    0.258
##    .feltobl           0.024    0.017    1.423    0.155
## 
## R-Square:
##                    Estimate
##     prvoice           0.569
##     feltobl           0.131
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     IE1               0.157    0.047    3.345    0.001
##     IE2               0.122    0.039    3.104    0.002
##     IE3              -0.161    0.290   -0.555    0.579
##     IE4              -0.304    0.268   -1.135    0.256

We see that the both indirect effects at Level 1 were positive and statistically significant. Both indirect effects at Level 2 were not statistically significant.

References

Xie, X.-Y., Ling, C.-D., Mo, S.-J., Luan, K. (2015). Linking colleague support to employees’ promotive voice: A moderated mediation model. PLoS ONE 10(7): e0132123. doi:10.1371/journal.pone.0132123. Downloaded April 21, 2024 from: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0132123

Comments

Popular posts from this blog

Factor analysis of EBI items: Tutorial with RStudio and EFA.dimensions package

Process model 7 using Hayes Process macro with RStudio

A tutorial on exploratory factor analysis with SPSS