# Help moving forward after modIndices() in lavaan

I’m running a model in SEM using lavaan. I’ve really run into a puzzle I can’t quite seem to solve and I would love to tap into everyone’s expertise to help provide some direction. Below I’ve added the steps I’ve taken so far to help shed light on my attempt at Sherlock Holms-ing this situation. Apologies in advance for the length of the post, but hopefully there’s some other brains out there that thrive on solving problems like these.

Please see the image below for my original model. The model is pretty basic with latent variable RL (responsive leadership) predicting latent variable PI (personal initiative) while controlling for manifest gender, manifest ethnicity, latent Control, and latent Cplex (complexity).

Additionally, here is the code I used for the original model:

SEMmodel2 <- '# Latent variables
RL =~ 1*RL_1 + RL_2 + RL_3 + RL_4 + RL_5 + RL_6 + RL_7 + RL_8 + RL_9 + RL_10 + RL_11 + RL_12
PI =~ 1*PI_1 + PI_2 + PI_3 + PI_4 + PI_5 + PI_6 + PI_7
Cplex =~ 1*Cplex_1 + Cplex_2 + Cplex_3 + Cplex_4
Control =~ 1*Cont_1 + Cont_2 + Cont_3 + Cont_4

#regressions
PI ~ RL + Age + Cplex + Control + Gen2 + Eth1 + Eth2 + Eth4 + Eth5'

SEMmodel2 <- lavaan(SEMmodel2, data = dat2, auto.var = TRUE, fixed.x = FALSE, int.ov.free = TRUE, estimator = "dwls")
summary(SEMmodel2, fit.measures = TRUE, standardized = TRUE)


When I ran this model, the model fit wasn’t the best (it wasn’t horrible, but it didn’t meet any cutoffs). For character limit sake, here are a quick summary of some important pieces of the output instead of the full output: CFI = .747, RMSEA = .157, SRMR = .169. Factor loadings between .56 and .86. PI to RL path had a beta of .54 and p <.001.

So, I decided to check the modification indices to see if there were any suggestions that theoretically made sense to modify. I ran the following code:

MOD <- modificationIndices(SEMmodel2, free.remove = TRUE, na.remove = TRUE, sort. = TRUE)
subset(MOD, mi >1000)


And this is where I get tripped up. The top modification indices that are suggested are some paths that are already in my model but just reversed. For instance, PI ~ RL is the main IV to DV path that is in the model. The modification indices are telling me to add RL ~ PI on top of that path, though. It does this with a couple of the covariates as well (control, cplex). See output below:

subset(MOD, mi >1000)
lhs op     rhs       mi   epc sepc.lv sepc.all sepc.nox
RL  ~      PI 4918.437 1.468   1.007    1.007    1.007
Control  ~      PI 4118.074 1.563   1.177    1.177    1.177
Control  ~      RL 3241.769 0.573   0.629    0.629    0.629
RL  ~ Control 3241.769 0.690   0.629    0.629    0.629
RL ~~ Control 3241.769 0.775   0.629    0.629    0.629
Cplex  ~      PI 2525.688 1.197   0.973    0.973    0.973
RL ~~   Cplex 1667.477 0.561   0.491    0.491    0.491
Cplex  ~      RL 1667.477 0.414   0.491    0.491    0.491
RL  ~   Cplex 1667.477 0.583   0.491    0.491    0.491


From these suggestions I tried two different model fixes – both have issues and I’m not sure how to move forward. Any advice or suggestions are welcome.

Fix 1: Add in path RL ~ PI

My inclination is that adding this arrow into the model – like the following image – would make it non-recursive and basically means that RL predicts PI and PI in turn predicts RL. If I add in that path to the model, the model fit is great but there are multiple betas that are above 1 which is concerning (see output below). I’m not as clear with the underlying mathematics that go into this type of analysis. Can someone explain the implications of this and what it would mean about steps forward?

Here is the output for adding in the non-recursive path:

summary(SEMmodel2_3, fit.measures = TRUE, standardized = TRUE)
lavaan 0.6-3 ended normally after 150 iterations

Optimization method                           NLMINB
Number of free parameters                         85

Number of observations                           506

Estimator                                       DWLS
Model Fit Test Statistic                    1330.226
Degrees of freedom                               476
P-value (Chi-square)                           0.000

Model test baseline model:

Minimum Function Test Statistic            24118.401
Degrees of freedom                               513
P-value                                        0.000

User model versus baseline model:

Comparative Fit Index (CFI)                    0.964
Tucker-Lewis Index (TLI)                       0.961

Root Mean Square Error of Approximation:

RMSEA                                          0.060
90 Percent Confidence Interval          0.056  0.063
P-value RMSEA <= 0.05                          0.000

Standardized Root Mean Square Residual:

SRMR                                           0.076

Parameter Estimates:

Information                                 Expected
Information saturated (h1) model        Unstructured
Standard Errors                             Standard

Latent Variables:
Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
RL =~
RL_1              1.000                               1.176    0.817
RL_2              0.968    0.029   32.879    0.000    1.139    0.858
RL_3              0.894    0.028   31.569    0.000    1.052    0.796
RL_4              1.000    0.031   32.592    0.000    1.177    0.812
RL_5              0.943    0.029   32.449    0.000    1.109    0.831
RL_6              0.992    0.031   32.417    0.000    1.168    0.844
RL_7              0.978    0.030   32.764    0.000    1.150    0.846
RL_8              0.877    0.028   31.560    0.000    1.032    0.801
RL_9              1.007    0.031   32.855    0.000    1.185    0.867
RL_10             0.781    0.025   31.113    0.000    0.919    0.782
RL_11             0.935    0.029   32.153    0.000    1.099    0.800
RL_12             0.891    0.028   31.675    0.000    1.048    0.804
PI =~
PI_1              1.000                               0.794    0.726
PI_2              0.897    0.036   25.161    0.000    0.712    0.632
PI_3              1.052    0.040   26.039    0.000    0.836    0.713
PI_4              1.226    0.044   27.748    0.000    0.974    0.790
PI_5              1.154    0.043   26.887    0.000    0.917    0.728
PI_6              1.002    0.041   24.271    0.000    0.795    0.625
PI_7              1.071    0.040   26.505    0.000    0.850    0.711
Cplex =~
Cplex_1           1.000                               0.807    0.484
Cplex_2           0.866    0.060   14.429    0.000    0.699    0.407
Cplex_3           1.654    0.090   18.335    0.000    1.335    0.864
Cplex_4           1.464    0.080   18.266    0.000    1.182    0.892
Control =~
Cont_1            1.000                               1.044    0.775
Cont_2            0.984    0.037   26.470    0.000    1.027    0.772
Cont_3            0.948    0.038   25.247    0.000    0.990    0.727
Cont_4            1.136    0.041   27.373    0.000    1.186    0.904

Regressions:
Estimate  Std.Err  z-value  P(>|z|)   Std.lv Std.all
PI ~
RL               -0.782    0.081   -9.631    0.000   -1.159   -1.159
Age               0.021    0.003    5.922    0.000    0.026    0.234
Cplex             1.098    0.094   11.672    0.000    1.117    1.117
Control           1.009    0.072   14.013    0.000    1.327    1.327
Gen2             -0.113    0.059   -1.915    0.056   -0.142   -0.068
Eth1              0.727    0.122    5.979    0.000    0.915    0.294
Eth2              0.025    0.089    0.278    0.781    0.031    0.009
Eth4              0.294    0.119    2.474    0.013    0.371    0.085
Eth5              0.878    0.405    2.166    0.030    1.106    0.120
RL ~
PI                1.548    0.063   24.687    0.000    1.045    1.045

Covariances:
Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
Age ~~
Gen2              0.481    0.197    2.440    0.015    0.481    0.112
Eth1             -0.490    0.098   -4.986    0.000   -0.490   -0.170
Eth2             -0.065    0.120   -0.541    0.589   -0.065   -0.025
Eth4             -0.219    0.076   -2.889    0.004   -0.219   -0.107
Eth5              0.005    0.045    0.112    0.911    0.005    0.005
Gen2 ~~
Eth1              0.011    0.007    1.562    0.118    0.011    0.073
Eth2             -0.004    0.006   -0.653    0.514   -0.004   -0.028
Eth4             -0.010    0.004   -2.344    0.019   -0.010   -0.092
Eth5             -0.002    0.002   -1.143    0.253   -0.002   -0.044
Eth1 ~~
Eth2             -0.011    0.002   -5.749    0.000   -0.011   -0.116
Eth4             -0.006    0.001   -4.731    0.000   -0.006   -0.088
Eth5             -0.001    0.001   -2.388    0.017   -0.001   -0.040
Eth2 ~~
Eth4             -0.005    0.001   -4.518    0.000   -0.005   -0.077
Eth5             -0.001    0.000   -2.359    0.018   -0.001   -0.035
Eth4 ~~
Eth5             -0.001    0.000   -2.270    0.023   -0.001   -0.027

Variances:
Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
.RL_1              0.690    0.149    4.624    0.000    0.690    0.333
.RL_2              0.466    0.135    3.439    0.001    0.466    0.264
.RL_3              0.640    0.129    4.965    0.000    0.640    0.366
.RL_4              0.714    0.146    4.905    0.000    0.714    0.340
.RL_5              0.552    0.134    4.112    0.000    0.552    0.310
.RL_6              0.552    0.151    3.648    0.000    0.552    0.288
.RL_7              0.527    0.133    3.970    0.000    0.527    0.285
.RL_8              0.593    0.134    4.415    0.000    0.593    0.358
.RL_9              0.462    0.140    3.310    0.001    0.462    0.248
.RL_10             0.538    0.117    4.580    0.000    0.538    0.389
.RL_11             0.681    0.126    5.413    0.000    0.681    0.360
.RL_12             0.602    0.135    4.471    0.000    0.602    0.354
.PI_1              0.565    0.079    7.134    0.000    0.565    0.472
.PI_2              0.761    0.112    6.803    0.000    0.761    0.600
.PI_3              0.676    0.116    5.816    0.000    0.676    0.492
.PI_4              0.573    0.109    5.272    0.000    0.573    0.376
.PI_5              0.746    0.101    7.367    0.000    0.746    0.470
.PI_6              0.985    0.127    7.738    0.000    0.985    0.609
.PI_7              0.706    0.101    6.998    0.000    0.706    0.494
.Cplex_1           2.132    0.153   13.931    0.000    2.132    0.766
.Cplex_2           2.457    0.154   16.004    0.000    2.457    0.834
.Cplex_3           0.605    0.195    3.106    0.002    0.605    0.254
.Cplex_4           0.360    0.154    2.332    0.020    0.360    0.205
.Cont_1            0.726    0.131    5.553    0.000    0.726    0.400
.Cont_2            0.716    0.129    5.535    0.000    0.716    0.404
.Cont_3            0.876    0.140    6.253    0.000    0.876    0.472
.Cont_4            0.313    0.138    2.265    0.024    0.313    0.182
.RL                1.343    0.080   16.793    0.000    0.970    0.970
.PI                0.285    0.090    3.166    0.002    0.452    0.452
Cplex             0.652    0.062   10.503    0.000    1.000    1.000
Control           1.090    0.063   17.203    0.000    1.000    1.000
Age              80.342    5.238   15.338    0.000   80.342    1.000
Gen2              0.231    0.006   39.102    0.000    0.231    1.000
Eth1              0.103    0.011    9.424    0.000    0.103    1.000
Eth2              0.084    0.011    8.027    0.000    0.084    1.000
Eth4              0.052    0.009    5.789    0.000    0.052    1.000
Eth5              0.012    0.005    2.497    0.013    0.012    1.000


Fix 2: Add in RL ~ Control path

The variable “Control” theoretically makes sense to use as a control variable for RL, so that would be a good path to enter. When I try adding the covariate “Control” to the IV (“RL”), the model fit increases but it is still not superb (i.e., still below the cutoffs; see output below).

summary(SEMmodel2_2, fit.measures = TRUE, standardized = TRUE)
lavaan 0.6-3 ended normally after 135 iterations

Optimization method                           NLMINB
Number of free parameters                         85

Number of observations                           506

Estimator                                       DWLS
Model Fit Test Statistic                    3153.444
Degrees of freedom                               476
P-value (Chi-square)                           0.000

Model test baseline model:

Minimum Function Test Statistic            24118.401
Degrees of freedom                               513
P-value                                        0.000

User model versus baseline model:

Comparative Fit Index (CFI)                    0.887
Tucker-Lewis Index (TLI)                       0.878

Root Mean Square Error of Approximation:

RMSEA                                          0.106
90 Percent Confidence Interval          0.102  0.109
P-value RMSEA <= 0.05                          0.000

Standardized Root Mean Square Residual:

SRMR                                           0.116

Parameter Estimates:

Information                                 Expected
Information saturated (h1) model        Unstructured
Standard Errors                             Standard

Latent Variables:
Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
RL =~
RL_1              1.000                               1.167    0.810
RL_2              0.973    0.031   31.535    0.000    1.135    0.855
RL_3              0.897    0.030   30.246    0.000    1.047    0.793
RL_4              1.001    0.032   31.216    0.000    1.168    0.806
RL_5              0.956    0.031   31.208    0.000    1.116    0.836
RL_6              1.003    0.032   31.135    0.000    1.170    0.846
RL_7              0.994    0.031   31.575    0.000    1.160    0.853
RL_8              0.891    0.029   30.376    0.000    1.040    0.808
RL_9              1.013    0.032   31.526    0.000    1.182    0.865
RL_10             0.793    0.026   29.943    0.000    0.926    0.788
RL_11             0.940    0.031   30.822    0.000    1.098    0.798
RL_12             0.896    0.030   30.240    0.000    1.045    0.802
PI =~
PI_1              1.000                               0.796    0.728
PI_2              0.890    0.035   25.129    0.000    0.709    0.629
PI_3              1.047    0.040   26.039    0.000    0.834    0.711
PI_4              1.226    0.044   27.801    0.000    0.977    0.792
PI_5              1.153    0.043   26.925    0.000    0.918    0.729
PI_6              0.996    0.041   24.261    0.000    0.794    0.624
PI_7              1.069    0.040   26.552    0.000    0.851    0.712
Cplex =~
Cplex_1           1.000                               0.981    0.588
Cplex_2           0.980    0.078   12.530    0.000    0.961    0.560
Cplex_3           1.291    0.094   13.695    0.000    1.266    0.819
Cplex_4           1.032    0.077   13.425    0.000    1.012    0.764
Control =~
Cont_1            1.000                               1.045    0.775
Cont_2            0.981    0.037   26.468    0.000    1.025    0.770
Cont_3            0.946    0.037   25.236    0.000    0.989    0.726
Cont_4            1.136    0.041   27.398    0.000    1.188    0.906

Regressions:
Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
PI ~
RL                0.196    0.020    9.692    0.000    0.288    0.288
Age               0.011    0.003    4.476    0.000    0.014    0.127
Cplex             0.435    0.034   12.887    0.000    0.536    0.536
Control           0.297    0.030    9.870    0.000    0.390    0.390
Gen2             -0.043    0.045   -0.960    0.337   -0.054   -0.026
Eth1              0.488    0.086    5.660    0.000    0.613    0.197
Eth2              0.040    0.072    0.560    0.575    0.051    0.015
Eth4              0.126    0.095    1.328    0.184    0.158    0.036
Eth5              0.494    0.258    1.918    0.055    0.621    0.067
RL ~
Control           0.714    0.036   19.995    0.000    0.639    0.639

Covariances:
Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
Age ~~
Gen2              0.481    0.197    2.440    0.015    0.481    0.112
Eth1             -0.490    0.098   -4.986    0.000   -0.490   -0.170
Eth2             -0.065    0.120   -0.541    0.589   -0.065   -0.025
Eth4             -0.219    0.076   -2.889    0.004   -0.219   -0.107
Eth5              0.005    0.045    0.112    0.911    0.005    0.005
Gen2 ~~
Eth1              0.011    0.007    1.562    0.118    0.011    0.073
Eth2             -0.004    0.006   -0.653    0.514   -0.004   -0.028
Eth4             -0.010    0.004   -2.344    0.019   -0.010   -0.092
Eth5             -0.002    0.002   -1.143    0.253   -0.002   -0.044
Eth1 ~~
Eth2             -0.011    0.002   -5.749    0.000   -0.011   -0.116
Eth4             -0.006    0.001   -4.731    0.000   -0.006   -0.088
Eth5             -0.001    0.001   -2.388    0.017   -0.001   -0.040
Eth2 ~~
Eth4             -0.005    0.001   -4.518    0.000   -0.005   -0.077
Eth5             -0.001    0.000   -2.359    0.018   -0.001   -0.035
Eth4 ~~
Eth5             -0.001    0.000   -2.270    0.023   -0.001   -0.027

Variances:
Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
.RL_1              0.712    0.150    4.753    0.000    0.712    0.343
.RL_2              0.474    0.136    3.482    0.000    0.474    0.269
.RL_3              0.649    0.130    5.011    0.000    0.649    0.372
.RL_4              0.735    0.146    5.023    0.000    0.735    0.350
.RL_5              0.537    0.135    3.968    0.000    0.537    0.301
.RL_6              0.546    0.152    3.582    0.000    0.546    0.285
.RL_7              0.505    0.134    3.767    0.000    0.505    0.273
.RL_8              0.576    0.135    4.265    0.000    0.576    0.348
.RL_9              0.469    0.140    3.340    0.001    0.469    0.251
.RL_10             0.525    0.118    4.446    0.000    0.525    0.380
.RL_11             0.685    0.127    5.409    0.000    0.685    0.363
.RL_12             0.607    0.135    4.480    0.000    0.607    0.357
.PI_1              0.562    0.079    7.086    0.000    0.562    0.470
.PI_2              0.766    0.112    6.853    0.000    0.766    0.604
.PI_3              0.679    0.116    5.849    0.000    0.679    0.494
.PI_4              0.568    0.109    5.226    0.000    0.568    0.373
.PI_5              0.744    0.101    7.339    0.000    0.744    0.469
.PI_6              0.988    0.127    7.764    0.000    0.988    0.611
.PI_7              0.705    0.101    6.987    0.000    0.705    0.493
.Cplex_1           1.821    0.175   10.416    0.000    1.821    0.654
.Cplex_2           2.022    0.178   11.387    0.000    2.022    0.686
.Cplex_3           0.785    0.211    3.727    0.000    0.785    0.329
.Cplex_4           0.733    0.156    4.707    0.000    0.733    0.417
.Cont_1            0.725    0.131    5.536    0.000    0.725    0.399
.Cont_2            0.720    0.129    5.571    0.000    0.720    0.407
.Cont_3            0.878    0.140    6.274    0.000    0.878    0.473
.Cont_4            0.309    0.138    2.232    0.026    0.309    0.180
.RL                0.805    0.046   17.395    0.000    0.591    0.591
.PI                0.182    0.025    7.243    0.000    0.287    0.287
Cplex             0.962    0.105    9.168    0.000    1.000    1.000
Control           1.092    0.063   17.211    0.000    1.000    1.000
Age              80.342    5.238   15.338    0.000   80.342    1.000
Gen2              0.231    0.006   39.102    0.000    0.231    1.000
Eth1              0.103    0.011    9.424    0.000    0.103    1.000
Eth2              0.084    0.011    8.027    0.000    0.084    1.000
Eth4              0.052    0.009    5.789    0.000    0.052    1.000
Eth5              0.012    0.005    2.497    0.013    0.012    1.000


Since the model was good but not great, I looked back to the modification indices one more time (see output below).

subset(MOD2_2, mi >1000)
lhs op     rhs       mi   epc sepc.lv sepc.all sepc.nox
Control  ~      PI 2582.365 1.722   1.312    1.312    1.312
Cplex  ~      PI 2521.305 1.196   0.971    0.971    0.971
Control  ~   Cplex 2511.789 0.773   0.725    0.725    0.725
Cplex  ~ Control 2511.789 0.681   0.725    0.725    0.725
Cplex ~~ Control 2511.789 0.744   0.725    0.725    0.725
Cplex  ~      RL 2312.928 0.461   0.548    0.548    0.548
RL  ~      PI 1729.922 1.308   0.893    0.893    0.893
RL  ~   Cplex 1668.640 0.585   0.492    0.492    0.492
RL ~~   Cplex 1668.640 0.563   0.639    0.639    0.639
Control =~ Cplex_4 1121.039 0.863   0.902    0.680    0.680
RL =~ Cplex_4 1092.501 0.612   0.715    0.539    0.539
PI =~ Cplex_4 1062.157 1.407   1.120    0.845    0.845


The modification indices suggest either more non-recursive paths or things that theoretically don't make sense. Indicating I should either go with the first fix of adding the non-recursive (RL ~ PI) path or just stick with fix 2 as my final model with just good but not great model fit. Would that be a correct conclusion?

Any help to my queries in italics or suggestions about how to move forward given this information would be wildly welcomed.

If there’s any other information you need me to add, please let me know and I will add it asap. Additionally, if I somehow missed another post that already answers this question (I scoured the site and couldn’t find anything but I could’ve missed something in my search), feel free to just direct me to that post instead. Thanks!

Try respecifying the model by adding covariances among all the exogenous variables, latent and observed. Check to ensure all the covariances you intend are there by using lavaanify(SEMmodel2, OPTIONS), where OPTIONS stands in for the options you included in the call to lavaan(). Refit the model once you have ensured your specifications are correct, and then check the modification indices again. This time, they should only refer to indicator disturbances (because the only restrictions you're making are on the covariances between the disturbances and the other observed variables (i.e., you're assuming them to be zero in the measurement models).