library(modsem)6 LMS and QML approaches
7 The Latent Moderated Structural Equations (LMS) and the Quasi Maximum Likelihood (QML) approach
Both the LMS- and QML approach works on most models, but interaction effects with endogenous can be a bit tricky to estimate (see the vignette. Both approaches (particularily the LMS approach) are quite computationally intensive, and are thus partly implemented in C++ (using Rcpp and RcppArmadillo). Additionally starting parameters are estimated using the double centering approach (and the means of the observed variables) are used to generate good starting parameters for faster convergence. If you want to see the progress of the estimation process you can use ´verbose = TRUE´.
7.1 A Simple Example
Here you can see an example of the LMS approach for a simple model. By default the summary function calculates fit measures compared to a null model (i.e., the same model without an interaction term).
library(modsem)
m1 <- '
# Outer Model
X =~ x1
X =~ x2 + x3
Z =~ z1 + z2 + z3
Y =~ y1 + y2 + y3
# Inner model
Y ~ X + Z
Y ~ X:Z
'
lms1 <- modsem(m1, oneInt, method = "lms")
summary(lms1, standardized = TRUE) # standardized estimates
#> Estimating null model
#> EM: Iteration = 1, LogLik = -17831.87, Change = -17831.875
#> EM: Iteration = 2, LogLik = -17831.87, Change = 0.000
#>
#> modsem (version 1.0.1):
#> Estimator LMS
#> Optimization method EM-NLMINB
#> Number of observations 2000
#> Number of iterations 92
#> Loglikelihood -14687.85
#> Akaike (AIC) 29437.71
#> Bayesian (BIC) 29611.34
#>
#> Numerical Integration:
#> Points of integration (per dim) 24
#> Dimensions 1
#> Total points of integration 24
#>
#> Fit Measures for H0:
#> Loglikelihood -17832
#> Akaike (AIC) 35723.75
#> Bayesian (BIC) 35891.78
#> Chi-square 17.52
#> Degrees of Freedom (Chi-square) 24
#> P-value (Chi-square) 0.826
#> RMSEA 0.000
#>
#> Comparative fit to H0 (no interaction effect)
#> Loglikelihood change 3144.02
#> Difference test (D) 6288.04
#> Degrees of freedom (D) 1
#> P-value (D) 0.000
#>
#> R-Squared:
#> Y 0.596
#> R-Squared Null-Model (H0):
#> Y 0.395
#> R-Squared Change:
#> Y 0.201
#>
#> Parameter Estimates:
#> Coefficients standardized
#> Information expected
#> Standard errors standard
#>
#> Latent Variables:
#> Estimate Std.Error z.value Pr(>|z|)
#> X =~
#> x1 0.926
#> x2 0.891 0.014 63.7 0.000
#> x3 0.912 0.014 67.4 0.000
#> Z =~
#> z1 0.927
#> z2 0.898 0.014 66.2 0.000
#> z3 0.913 0.013 67.7 0.000
#> Y =~
#> y1 0.969
#> y2 0.954 0.009 105.3 0.000
#> y3 0.961 0.009 111.2 0.000
#>
#> Regressions:
#> Estimate Std.Error z.value Pr(>|z|)
#> Y ~
#> X 0.427 0.020 21.7 0.000
#> Z 0.370 0.019 20.0 0.000
#> X:Z 0.454 0.018 25.9 0.000
#>
#> Covariances:
#> Estimate Std.Error z.value Pr(>|z|)
#> X ~~
#> Z 0.199 0.024 8.3 0.000
#>
#> Variances:
#> Estimate Std.Error z.value Pr(>|z|)
#> x1 0.142 0.007 19.3 0.000
#> x2 0.206 0.009 23.5 0.000
#> x3 0.169 0.008 21.5 0.000
#> z1 0.141 0.008 18.4 0.000
#> z2 0.193 0.009 22.6 0.000
#> z3 0.167 0.008 20.9 0.000
#> y1 0.061 0.003 17.8 0.000
#> y2 0.090 0.004 22.7 0.000
#> y3 0.077 0.004 20.7 0.000
#> X 1.000 0.017 60.6 0.000
#> Z 1.000 0.018 54.9 0.000
#> Y 0.404 0.015 26.3 0.000Here you can see the same example using the QML approach.
qml1 <- modsem(m1, oneInt, method = "qml")
summary(qml1)
#> Estimating null model
#> Starting M-step
#>
#> modsem (version 1.0.1):
#> Estimator QML
#> Optimization method NLMINB
#> Number of observations 2000
#> Number of iterations 110
#> Loglikelihood -17496.22
#> Akaike (AIC) 35054.43
#> Bayesian (BIC) 35228.06
#>
#> Fit Measures for H0:
#> Loglikelihood -17832
#> Akaike (AIC) 35723.75
#> Bayesian (BIC) 35891.78
#> Chi-square 17.52
#> Degrees of Freedom (Chi-square) 24
#> P-value (Chi-square) 0.826
#> RMSEA 0.000
#>
#> Comparative fit to H0 (no interaction effect)
#> Loglikelihood change 335.66
#> Difference test (D) 671.32
#> Degrees of freedom (D) 1
#> P-value (D) 0.000
#>
#> R-Squared:
#> Y 0.607
#> R-Squared Null-Model (H0):
#> Y 0.395
#> R-Squared Change:
#> Y 0.211
#>
#> Parameter Estimates:
#> Coefficients unstandardized
#> Information observed
#> Standard errors standard
#>
#> Latent Variables:
#> Estimate Std.Error z.value Pr(>|z|)
#> X =~
#> x1 1.000
#> x2 0.803 0.013 63.97 0.000
#> x3 0.914 0.013 67.80 0.000
#> Z =~
#> z1 1.000
#> z2 0.810 0.012 65.12 0.000
#> z3 0.881 0.013 67.62 0.000
#> Y =~
#> y1 1.000
#> y2 0.798 0.007 107.57 0.000
#> y3 0.899 0.008 112.55 0.000
#>
#> Regressions:
#> Estimate Std.Error z.value Pr(>|z|)
#> Y ~
#> X 0.674 0.032 20.94 0.000
#> Z 0.566 0.030 18.95 0.000
#> X:Z 0.712 0.028 25.46 0.000
#>
#> Intercepts:
#> Estimate Std.Error z.value Pr(>|z|)
#> x1 1.023 0.024 42.89 0.000
#> x2 1.215 0.020 60.99 0.000
#> x3 0.919 0.022 41.48 0.000
#> z1 1.012 0.024 41.57 0.000
#> z2 1.206 0.020 59.27 0.000
#> z3 0.916 0.022 42.06 0.000
#> y1 1.038 0.033 31.45 0.000
#> y2 1.221 0.027 45.49 0.000
#> y3 0.954 0.030 31.86 0.000
#> Y 0.000
#>
#> Covariances:
#> Estimate Std.Error z.value Pr(>|z|)
#> X ~~
#> Z 0.200 0.024 8.24 0.000
#>
#> Variances:
#> Estimate Std.Error z.value Pr(>|z|)
#> x1 0.158 0.009 18.14 0.000
#> x2 0.162 0.007 23.19 0.000
#> x3 0.165 0.008 20.82 0.000
#> z1 0.166 0.009 18.34 0.000
#> z2 0.159 0.007 22.62 0.000
#> z3 0.158 0.008 20.71 0.000
#> y1 0.159 0.009 17.98 0.000
#> y2 0.154 0.007 22.67 0.000
#> y3 0.164 0.008 20.71 0.000
#> X 0.983 0.036 27.00 0.000
#> Z 1.019 0.038 26.95 0.000
#> Y 0.943 0.038 24.87 0.0007.2 A more complicated example
Here you can see an example of a more complicated example using the model from the theory of planned behaviour (TPB), where there are two endogenous variables, where there is an interaction between an endogenous and exogenous variable. When estimating more complicated models with the LMS-approach, it is recommended that you increase the number of nodes used for numerical integration. By default the number of nodes is set to 16, and can be increased using the nodes argument. The argument has no effect on the QML approach. You can also get robust standard errors by setting robust.se = TRUE in the modsem() function.
Note: If you want the lms-approach to give as similar results as possible to mplus, you would have to increase the number of nodes (e.g., nodes = 100), and lower the convergence criterion (e.g., conv_crit = 1e-5).
# ATT = Attitude,
# PBC = Perceived Behavioural Control
# INT = Intention
# SN = Subjective Norms
# BEH = Behaviour
tpb <- '
# Outer Model (Based on Hagger et al., 2007)
ATT =~ att1 + att2 + att3 + att4 + att5
SN =~ sn1 + sn2
PBC =~ pbc1 + pbc2 + pbc3
INT =~ int1 + int2 + int3
BEH =~ b1 + b2
# Inner Model (Based on Steinmetz et al., 2011)
INT ~ ATT + SN + PBC
BEH ~ INT + PBC
BEH ~ INT:PBC
'
lms2 <- modsem(tpb, TPB, method = "lms", nodes = 32)
summary(lms2)
#> Estimating null model
#> EM: Iteration = 1, LogLik = -26393.22, Change = -26393.223
#> EM: Iteration = 2, LogLik = -26393.22, Change = 0.000
#>
#> modsem (version 1.0.1):
#> Estimator LMS
#> Optimization method EM-NLMINB
#> Number of observations 2000
#> Number of iterations 55
#> Loglikelihood -23439.03
#> Akaike (AIC) 46986.05
#> Bayesian (BIC) 47288.5
#>
#> Numerical Integration:
#> Points of integration (per dim) 32
#> Dimensions 1
#> Total points of integration 32
#>
#> Fit Measures for H0:
#> Loglikelihood -26393
#> Akaike (AIC) 52892.45
#> Bayesian (BIC) 53189.29
#> Chi-square 66.27
#> Degrees of Freedom (Chi-square) 82
#> P-value (Chi-square) 0.897
#> RMSEA 0.000
#>
#> Comparative fit to H0 (no interaction effect)
#> Loglikelihood change 2954.20
#> Difference test (D) 5908.39
#> Degrees of freedom (D) 1
#> P-value (D) 0.000
#>
#> R-Squared:
#> INT 0.364
#> BEH 0.259
#> R-Squared Null-Model (H0):
#> INT 0.367
#> BEH 0.210
#> R-Squared Change:
#> INT -0.003
#> BEH 0.049
#>
#> Parameter Estimates:
#> Coefficients unstandardized
#> Information expected
#> Standard errors standard
#>
#> Latent Variables:
#> Estimate Std.Error z.value Pr(>|z|)
#> PBC =~
#> pbc1 1.000
#> pbc2 0.914 0.013 69.72 0.000
#> pbc3 0.802 0.012 65.65 0.000
#> ATT =~
#> att1 1.000
#> att2 0.878 0.013 70.13 0.000
#> att3 0.789 0.012 65.48 0.000
#> att4 0.695 0.011 60.51 0.000
#> att5 0.887 0.013 69.99 0.000
#> SN =~
#> sn1 1.000
#> sn2 0.889 0.017 52.23 0.000
#> INT =~
#> int1 1.000
#> int2 0.913 0.016 58.89 0.000
#> int3 0.807 0.014 55.85 0.000
#> BEH =~
#> b1 1.000
#> b2 0.959 0.032 30.44 0.000
#>
#> Regressions:
#> Estimate Std.Error z.value Pr(>|z|)
#> INT ~
#> PBC 0.217 0.029 7.44 0.000
#> ATT 0.214 0.026 8.15 0.000
#> SN 0.176 0.028 6.40 0.000
#> BEH ~
#> PBC 0.233 0.022 10.40 0.000
#> INT 0.188 0.025 7.59 0.000
#> PBC:INT 0.205 0.019 10.95 0.000
#>
#> Intercepts:
#> Estimate Std.Error z.value Pr(>|z|)
#> pbc1 0.991 0.022 45.69 0.000
#> pbc2 0.979 0.020 48.32 0.000
#> pbc3 0.986 0.018 54.20 0.000
#> att1 1.009 0.023 43.43 0.000
#> att2 1.003 0.021 48.38 0.000
#> att3 1.013 0.019 53.05 0.000
#> att4 0.996 0.017 57.18 0.000
#> att5 0.988 0.021 47.20 0.000
#> sn1 1.001 0.023 42.84 0.000
#> sn2 1.006 0.021 47.87 0.000
#> int1 1.011 0.021 48.20 0.000
#> int2 1.009 0.020 51.35 0.000
#> int3 1.003 0.018 55.80 0.000
#> b1 0.999 0.021 47.19 0.000
#> b2 1.017 0.020 51.24 0.000
#> INT 0.000
#> BEH 0.000
#>
#> Covariances:
#> Estimate Std.Error z.value Pr(>|z|)
#> PBC ~~
#> ATT 0.668 0.021 31.97 0.000
#> SN 0.668 0.022 30.93 0.000
#> ATT ~~
#> SN 0.623 0.019 33.42 0.000
#>
#> Variances:
#> Estimate Std.Error z.value Pr(>|z|)
#> pbc1 0.148 0.008 18.78 0.000
#> pbc2 0.159 0.007 21.62 0.000
#> pbc3 0.155 0.006 23.83 0.000
#> att1 0.167 0.007 23.51 0.000
#> att2 0.150 0.006 24.33 0.000
#> att3 0.160 0.006 26.53 0.000
#> att4 0.163 0.006 27.59 0.000
#> att5 0.159 0.006 25.12 0.000
#> sn1 0.178 0.015 12.11 0.000
#> sn2 0.156 0.012 13.15 0.000
#> int1 0.157 0.009 18.25 0.000
#> int2 0.160 0.008 20.69 0.000
#> int3 0.168 0.007 23.73 0.000
#> b1 0.185 0.019 9.81 0.000
#> b2 0.136 0.017 7.96 0.000
#> PBC 0.947 0.017 55.36 0.000
#> ATT 0.992 0.014 69.43 0.000
#> SN 0.981 0.015 64.02 0.000
#> INT 0.491 0.020 24.97 0.000
#> BEH 0.456 0.023 19.84 0.000
qml2 <- modsem(tpb, TPB, method = "qml")
summary(qml2, standardized = TRUE) # standardized estimates
#> Estimating null model
#> Starting M-step
#>
#> modsem (version 1.0.1):
#> Estimator QML
#> Optimization method NLMINB
#> Number of observations 2000
#> Number of iterations 74
#> Loglikelihood -26326.25
#> Akaike (AIC) 52760.5
#> Bayesian (BIC) 53062.95
#>
#> Fit Measures for H0:
#> Loglikelihood -26393
#> Akaike (AIC) 52892.45
#> Bayesian (BIC) 53189.29
#> Chi-square 66.27
#> Degrees of Freedom (Chi-square) 82
#> P-value (Chi-square) 0.897
#> RMSEA 0.000
#>
#> Comparative fit to H0 (no interaction effect)
#> Loglikelihood change 66.97
#> Difference test (D) 133.95
#> Degrees of freedom (D) 1
#> P-value (D) 0.000
#>
#> R-Squared:
#> INT 0.366
#> BEH 0.263
#> R-Squared Null-Model (H0):
#> INT 0.367
#> BEH 0.210
#> R-Squared Change:
#> INT 0.000
#> BEH 0.053
#>
#> Parameter Estimates:
#> Coefficients standardized
#> Information observed
#> Standard errors standard
#>
#> Latent Variables:
#> Estimate Std.Error z.value Pr(>|z|)
#> PBC =~
#> pbc1 0.933
#> pbc2 0.913 0.013 69.47 0.000
#> pbc3 0.894 0.014 66.10 0.000
#> ATT =~
#> att1 0.925
#> att2 0.915 0.013 71.56 0.000
#> att3 0.892 0.013 66.38 0.000
#> att4 0.865 0.014 61.00 0.000
#> att5 0.912 0.013 70.85 0.000
#> SN =~
#> sn1 0.921
#> sn2 0.913 0.017 52.61 0.000
#> INT =~
#> int1 0.912
#> int2 0.895 0.015 59.05 0.000
#> int3 0.866 0.016 55.73 0.000
#> BEH =~
#> b1 0.877
#> b2 0.900 0.028 31.71 0.000
#>
#> Regressions:
#> Estimate Std.Error z.value Pr(>|z|)
#> INT ~
#> PBC 0.243 0.033 7.35 0.000
#> ATT 0.242 0.030 8.16 0.000
#> SN 0.199 0.031 6.37 0.000
#> BEH ~
#> PBC 0.289 0.028 10.37 0.000
#> INT 0.212 0.028 7.69 0.000
#> PBC:INT 0.227 0.020 11.33 0.000
#>
#> Covariances:
#> Estimate Std.Error z.value Pr(>|z|)
#> PBC ~~
#> ATT 0.692 0.030 23.45 0.000
#> SN 0.695 0.030 23.08 0.000
#> ATT ~~
#> SN 0.634 0.029 21.70 0.000
#>
#> Variances:
#> Estimate Std.Error z.value Pr(>|z|)
#> pbc1 0.130 0.007 18.39 0.000
#> pbc2 0.166 0.008 21.43 0.000
#> pbc3 0.201 0.008 23.89 0.000
#> att1 0.144 0.006 23.53 0.000
#> att2 0.164 0.007 24.71 0.000
#> att3 0.204 0.008 26.38 0.000
#> att4 0.252 0.009 27.65 0.000
#> att5 0.168 0.007 24.93 0.000
#> sn1 0.153 0.013 12.09 0.000
#> sn2 0.167 0.013 13.26 0.000
#> int1 0.168 0.009 18.11 0.000
#> int2 0.199 0.010 20.41 0.000
#> int3 0.249 0.011 23.55 0.000
#> b1 0.231 0.023 10.12 0.000
#> b2 0.191 0.024 8.10 0.000
#> PBC 1.000 0.037 27.07 0.000
#> ATT 1.000 0.037 26.93 0.000
#> SN 1.000 0.040 25.22 0.000
#> INT 0.634 0.026 24.64 0.000
#> BEH 0.737 0.037 20.17 0.000