--- title: "The-Structural-Equation-Model" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{The-Structural-Equation-Model} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```r library(lessSEM) ``` The lessSEM package comes with a custom implementation of structural equation models (SEM). This implementation supports full-information-maximum-likelihood computation in case of missing data and could also be used by other packages. Identical to [regsem](https://github.com/Rjacobucci/regsem) lessSEM also builds on [lavaan](https://github.com/yrosseel/lavaan) to set up the model. That is, if you are already familiar with lavaan, setting up models with lessSEM should be relatively easy. We will use the political democracy example from the sem documentation of lavaan in the following: ```r library(lavaan) # see ?lavaan::sem model <- ' # latent variable definitions ind60 =~ x1 + x2 + x3 dem60 =~ y1 + a*y2 + b*y3 + c*y4 dem65 =~ y5 + a*y6 + b*y7 + c*y8 # regressions dem60 ~ ind60 dem65 ~ ind60 + dem60 # residual correlations y1 ~~ y5 y2 ~~ y4 + y6 y3 ~~ y7 y4 ~~ y8 y6 ~~ y8 ' lavaanModel <- sem(model, data = PoliticalDemocracy) ``` ## From lavaan to lessSEM To translate the model from lavaan to lessSEM, we have to use the `lessSEM:::.SEMFromLavaan` function. Importantly, this function is not exported by lessSEM. That is, you _must_ use the three colons as shown above to access this function! ```r library(lessSEM) # won't work: mySEM <- .SEMFromLavaan(lavaanModel = lavaanModel) # will work: mySEM <- lessSEM:::.SEMFromLavaan(lavaanModel = lavaanModel) ``` ```r show(mySEM) #> Internal C++ model representation of lessSEM #> Parameters: #> ind60=~x2 ind60=~x3 a b c dem60~ind60 dem65~ind60 dem65~dem60 y1~~y5 #> 2.1796566 1.8182100 1.1907820 1.1745407 1.2509789 1.4713302 0.6004746 0.8650430 0.5825389 #> y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1 x2~~x2 x3~~x3 y1~~y1 #> 1.4402477 2.1829448 0.7115901 0.3627964 1.3717741 0.0813878 0.1204271 0.4666596 1.8546417 #> y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60 #> 7.5813926 4.9556766 3.2245521 2.3130404 4.9681408 3.5600367 3.3076854 0.4485989 3.8753039 #> dem65~~dem65 x1~1 x2~1 x3~1 y1~1 y2~1 y3~1 y4~1 y5~1 #> 0.1644633 5.0543838 4.7921946 3.5576898 5.4646667 4.2564429 6.5631103 4.4525330 5.1362519 #> y6~1 y7~1 y8~1 #> 2.9780741 6.1962639 4.0433897 #> #> Objective value: 3097.6361581071 ``` The `lessSEM:::.SEMFromLavaan` function comes with some additional arguments to fine tune the initialization of the model. 1. `whichPars`: with the `whichPars` arguments, we can change which parameters are used in the mySEM created above. By default, we will use the estimates (`whichPars = "est"`) of the lavaan model, but we could also use the starting values (`whichPars = "start"`) or supply custom parameter values 2. `fit`: When `fit = TRUE`, lessSEM will fit the model once and compare the fitting function value to that of the lavaanModel. If you supplied parameters other than "est", this should be set to `fit = FALSE` 3. `addMeans`: Should a mean structure be added? It is currenlty recomended to set this to `TRUE` 4. `activeSet`: This allows for only using part of the data set. This can be useful for cross-validation. 5. `dataSet`: This allows for passing a different data set to mySEM. This can be useful for cross-validation. In most cases, we recommend setting up the model as shown above, with none of the additional arguments being used. ## Working with the Rcpp_SEMCpp class The mySEM object is implemented in C++ to make everything run faster. The underlying class is `Rcpp_SEMCpp` and was created using the wonderful [Rcpp](https://github.com/RcppCore/Rcpp) and [RcppArmadillo](https://github.com/RcppCore/RcppArmadillo) packages. ```r class(mySEM) #> [1] "Rcpp_SEMCpp" #> attr(,"package") #> [1] "lessSEM" ``` You can access its elements using the dollar-operator: ```r mySEM$A #> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] #> [1,] 0.0000000 0.000000 0.000000 0 0 0 0 0 0 0 0 0 0 0 #> [2,] 1.4713302 0.000000 0.000000 0 0 0 0 0 0 0 0 0 0 0 #> [3,] 0.6004746 0.865043 0.000000 0 0 0 0 0 0 0 0 0 0 0 #> [4,] 1.0000000 0.000000 0.000000 0 0 0 0 0 0 0 0 0 0 0 #> [5,] 2.1796566 0.000000 0.000000 0 0 0 0 0 0 0 0 0 0 0 #> [6,] 1.8182100 0.000000 0.000000 0 0 0 0 0 0 0 0 0 0 0 #> [7,] 0.0000000 1.000000 0.000000 0 0 0 0 0 0 0 0 0 0 0 #> [8,] 0.0000000 1.190782 0.000000 0 0 0 0 0 0 0 0 0 0 0 #> [9,] 0.0000000 1.174541 0.000000 0 0 0 0 0 0 0 0 0 0 0 #> [10,] 0.0000000 1.250979 0.000000 0 0 0 0 0 0 0 0 0 0 0 #> [11,] 0.0000000 0.000000 1.000000 0 0 0 0 0 0 0 0 0 0 0 #> [12,] 0.0000000 0.000000 1.190782 0 0 0 0 0 0 0 0 0 0 0 #> [13,] 0.0000000 0.000000 1.174541 0 0 0 0 0 0 0 0 0 0 0 #> [14,] 0.0000000 0.000000 1.250979 0 0 0 0 0 0 0 0 0 0 0 ``` Note that, identical to regsem, the model is implemented with the RAM notation (McArdle & McDonald, 1984). If you are not familiar with this notation, Fox (2006) provides a short introduction. However, you won't need to know the details for the time being. Instead, we will focus on how to get and set the parameters, fit the model, get its gradients, etc. ### Accessing the Parameters The parameters of the model can be accessed with the `lessSEM:::.getParameters` function: ```r (myParameters <- lessSEM:::.getParameters(mySEM)) #> ind60=~x2 ind60=~x3 a b c dem60~ind60 dem65~ind60 dem65~dem60 y1~~y5 #> 2.1796566 1.8182100 1.1907820 1.1745407 1.2509789 1.4713302 0.6004746 0.8650430 0.5825389 #> y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1 x2~~x2 x3~~x3 y1~~y1 #> 1.4402477 2.1829448 0.7115901 0.3627964 1.3717741 0.0813878 0.1204271 0.4666596 1.8546417 #> y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60 #> 7.5813926 4.9556766 3.2245521 2.3130404 4.9681408 3.5600367 3.3076854 0.4485989 3.8753039 #> dem65~~dem65 x1~1 x2~1 x3~1 y1~1 y2~1 y3~1 y4~1 y5~1 #> 0.1644633 5.0543838 4.7921946 3.5576898 5.4646667 4.2564429 6.5631103 4.4525330 5.1362519 #> y6~1 y7~1 y8~1 #> 2.9780741 6.1962639 4.0433897 ``` The naming is identical to that of the lavaanModel. By default, the parameters are returned in the transformed format. This requires some more explanation: In lessSEM we assume that negative variances are outside of the parameter space. That is, negative variances are **not** allowed (this is different from lavaan!). To ensure that all variances are positive, we use a transformation: Say we are interested in the variance `ind60~~ind60`. Internally, there is a parameter called `x1~~x1` and this parameter has a `rawValue` and a transformed value (called just `value`). We can access these values with: ```r mySEM$getParameters() #> label value rawValue location isTransformation #> 1 ind60=~x2 2.1796566 2.1796566 Amatrix FALSE #> 2 ind60=~x3 1.8182100 1.8182100 Amatrix FALSE #> 3 a 1.1907820 1.1907820 Amatrix FALSE #> 4 b 1.1745407 1.1745407 Amatrix FALSE #> 5 c 1.2509789 1.2509789 Amatrix FALSE #> 6 dem60~ind60 1.4713302 1.4713302 Amatrix FALSE #> 7 dem65~ind60 0.6004746 0.6004746 Amatrix FALSE #> 8 dem65~dem60 0.8650430 0.8650430 Amatrix FALSE #> 9 y1~~y5 0.5825389 0.5825389 Smatrix FALSE #> 10 y2~~y4 1.4402477 1.4402477 Smatrix FALSE #> 11 y2~~y6 2.1829448 2.1829448 Smatrix FALSE #> 12 y3~~y7 0.7115901 0.7115901 Smatrix FALSE #> 13 y4~~y8 0.3627964 0.3627964 Smatrix FALSE #> 14 y6~~y8 1.3717741 1.3717741 Smatrix FALSE #> 15 x1~~x1 0.0813878 -2.5085299 Smatrix FALSE #> 16 x2~~x2 0.1204271 -2.1167106 Smatrix FALSE #> 17 x3~~x3 0.4666596 -0.7621551 Smatrix FALSE #> 18 y1~~y1 1.8546417 0.6176915 Smatrix FALSE #> 19 y2~~y2 7.5813926 2.0256969 Smatrix FALSE #> 20 y3~~y3 4.9556766 1.6005337 Smatrix FALSE #> 21 y4~~y4 3.2245521 1.1707941 Smatrix FALSE #> 22 y5~~y5 2.3130404 0.8385629 Smatrix FALSE #> 23 y6~~y6 4.9681408 1.6030457 Smatrix FALSE #> 24 y7~~y7 3.5600367 1.2697708 Smatrix FALSE #> 25 y8~~y8 3.3076854 1.1962487 Smatrix FALSE #> 26 ind60~~ind60 0.4485989 -0.8016262 Smatrix FALSE #> 27 dem60~~dem60 3.8753039 1.3546241 Smatrix FALSE #> 28 dem65~~dem65 0.1644633 -1.8050678 Smatrix FALSE #> 29 x1~1 5.0543838 5.0543838 Mvector FALSE #> 30 x2~1 4.7921946 4.7921946 Mvector FALSE #> 31 x3~1 3.5576898 3.5576898 Mvector FALSE #> 32 y1~1 5.4646667 5.4646667 Mvector FALSE #> 33 y2~1 4.2564429 4.2564429 Mvector FALSE #> 34 y3~1 6.5631103 6.5631103 Mvector FALSE #> 35 y4~1 4.4525330 4.4525330 Mvector FALSE #> 36 y5~1 5.1362519 5.1362519 Mvector FALSE #> 37 y6~1 2.9780741 2.9780741 Mvector FALSE #> 38 y7~1 6.1962639 6.1962639 Mvector FALSE #> 39 y8~1 4.0433897 4.0433897 Mvector FALSE ``` For all parameters which are _not_ variances, the `rawValue` will be identical to the `value`. For variances, the `rawValue` can be any real value. The `value` itself is then computed as $e^{\text{rawValue}}$; this ensures that the `value` is always positive. You can access the raw values as follows: ```r lessSEM:::.getParameters(mySEM, raw = TRUE) #> ind60=~x2 ind60=~x3 a b c dem60~ind60 dem65~ind60 dem65~dem60 y1~~y5 #> 2.1796566 1.8182100 1.1907820 1.1745407 1.2509789 1.4713302 0.6004746 0.8650430 0.5825389 #> y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1 x2~~x2 x3~~x3 y1~~y1 #> 1.4402477 2.1829448 0.7115901 0.3627964 1.3717741 -2.5085299 -2.1167106 -0.7621551 0.6176915 #> y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60 #> 2.0256969 1.6005337 1.1707941 0.8385629 1.6030457 1.2697708 1.1962487 -0.8016262 1.3546241 #> dem65~~dem65 x1~1 x2~1 x3~1 y1~1 y2~1 y3~1 y4~1 y5~1 #> -1.8050678 5.0543838 4.7921946 3.5576898 5.4646667 4.2564429 6.5631103 4.4525330 5.1362519 #> y6~1 y7~1 y8~1 #> 2.9780741 6.1962639 4.0433897 ``` Note that the raw value for `ind60~~ind60` is negative while the transformed value is positive. ### Changing the Parameters Being able to change the parameters is essential for fitting a model. In lessSEM, this is facilitated by the `lessSEM:::.setParameters` function: ```r # first, let's change one of the parameters: myParameters["a"] <- 1 # now, let's change the parameters of the model mySEM <- lessSEM:::.setParameters(SEM = mySEM, # the model labels = names(myParameters), # names of the parameters values = myParameters, # values of the parameters raw = FALSE) ``` Note that we had to specify if the parameters in `myParameters` are given in raw format. Here, we already used the transformed parameters, so we set `raw = FALSE`. Using the raw parameters instead would look as follows: ```r myParameters <- lessSEM:::.getParameters(mySEM, raw = TRUE) # first, let's change one of the parameters: myParameters["a"] <- 1 # now, let's change the parameters of the model mySEM <- lessSEM:::.setParameters(SEM = mySEM, # the model labels = names(myParameters), # names of the parameters values = myParameters, # values of the parameters raw = TRUE) ``` Let's check the parameters: ```r lessSEM:::.getParameters(mySEM) #> ind60=~x2 ind60=~x3 a b c dem60~ind60 dem65~ind60 dem65~dem60 y1~~y5 #> 2.1796566 1.8182100 1.0000000 1.1745407 1.2509789 1.4713302 0.6004746 0.8650430 0.5825389 #> y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1 x2~~x2 x3~~x3 y1~~y1 #> 1.4402477 2.1829448 0.7115901 0.3627964 1.3717741 0.0813878 0.1204271 0.4666596 1.8546417 #> y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60 #> 7.5813926 4.9556766 3.2245521 2.3130404 4.9681408 3.5600367 3.3076854 0.4485989 3.8753039 #> dem65~~dem65 x1~1 x2~1 x3~1 y1~1 y2~1 y3~1 y4~1 y5~1 #> 0.1644633 5.0543838 4.7921946 3.5576898 5.4646667 4.2564429 6.5631103 4.4525330 5.1362519 #> y6~1 y7~1 y8~1 #> 2.9780741 6.1962639 4.0433897 ``` Note that `a` now has the value `1`. ## Fitting the model To compute the -2-log-likelihood of the model, we use the `lessSEM:::.fit` function: ```r mySEM <- lessSEM:::.fit(SEM = mySEM) ``` The -2-log-likelihood can be accessed with: ```r mySEM$objectiveValue #> [1] 3100.741 ``` ## Computing the gradients To compute the gradients, use the `lessSEM:::.getGradients` function. Gradients can be computed for the transformed parameters ```r lessSEM:::.getGradients(mySEM, raw = FALSE) #> ind60=~x2 ind60=~x3 a b c dem60~ind60 dem65~ind60 dem65~dem60 y1~~y5 #> 0.361097622 0.105564095 -32.837359814 2.453232158 17.076222049 -0.450648533 -1.078272542 -5.357920036 0.004650747 #> y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1 x2~~x2 x3~~x3 y1~~y1 #> -0.157923486 -0.567076177 0.161163293 0.266869495 -0.271533827 0.230054158 0.306685898 -0.093753843 -0.015516535 #> y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60 #> -0.343353554 0.099359383 0.137753880 0.131454473 -0.330083693 0.073331567 0.148964628 -0.103960291 -0.252921392 #> dem65~~dem65 x1~1 x2~1 x3~1 y1~1 y2~1 y3~1 y4~1 y5~1 #> -1.955349241 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 #> y6~1 y7~1 y8~1 #> 0.000000000 0.000000000 0.000000000 ``` or for the raw parameters ```r lessSEM:::.getGradients(mySEM, raw = TRUE) #> ind60=~x2 ind60=~x3 a b c dem60~ind60 dem65~ind60 dem65~dem60 y1~~y5 #> 0.361097622 0.105564095 -32.837359814 2.453232158 17.076222049 -0.450648533 -1.078272542 -5.357920036 0.004650747 #> y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1 x2~~x2 x3~~x3 y1~~y1 #> -0.157923486 -0.567076177 0.161163293 0.266869495 -0.271533827 0.018723602 0.036933297 -0.043751134 -0.028777612 #> y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60 #> -2.603098086 0.492392971 0.444194569 0.304059508 -1.639902248 0.261063066 0.492728130 -0.046636468 -0.980147244 #> dem65~~dem65 x1~1 x2~1 x3~1 y1~1 y2~1 y3~1 y4~1 y5~1 #> -0.321583201 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 #> y6~1 y7~1 y8~1 #> 0.000000000 0.000000000 0.000000000 ``` ## Computing the Hessian To compute the Hessian, use the `lessSEM:::.getHessian` function. The Hessian can be computed for the transformed parameters ```r lessSEM:::.getHessian(mySEM, raw = FALSE) ``` or for the raw parameters ```r lessSEM:::.getHessian(mySEM, raw = TRUE) ``` ## Computing the Scores To compute the scores (derivative of the -2-log-likelihood for each person), use the `lessSEM:::.getScores` function. The scores can be computed for the transformed parameters ```r lessSEM:::.getScores(mySEM, raw = FALSE) ``` or for the raw parameters ```r lessSEM:::.getScores(mySEM, raw = TRUE) ``` ## Using lessSEM with general purpose optimizers The most important part about the whole SEM implementation mentioned above is that we can use it flexibly with different optimizers. For instance, we may want to try out the BFGS optimizer from optim. **Important**: We highly recommend that you use the raw parameters for any optimization. Using the non-raw parameters can cause errors and unnecessary headaches! Let's have a look at the `optim` function: ```r args(optim) #> function (par, fn, gr = NULL, ..., method = c("Nelder-Mead", #> "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"), lower = -Inf, #> upper = Inf, control = list(), hessian = FALSE) #> NULL ``` Note that the function requires a `par` argument - the parameter estimates - a `fn` argument - the fitting function - and also allows for the gradients to be passed to the function using the `gr` argument. We could build such functions based on the `lessSEM:::.fit` and `lessSEM:::.getGradients` functions shown above, however for convenience such wrappers are already implemented in lessSEM. The fitting function is called with `lessSEM:::.fitFunction` and the gradient function is called `lessSEM:::.gradientFunction`. Both expect a vector with parameters, a SEM, and an argument specifying if the parameters are in raw format. We can use this in optim as follows: ```r # let's get the starting values: par <- lessSEM:::.getParameters(mySEM, raw = TRUE) # important: Use raw = TRUE! print(par) #> ind60=~x2 ind60=~x3 a b c dem60~ind60 dem65~ind60 dem65~dem60 y1~~y5 #> 2.1796566 1.8182100 1.0000000 1.1745407 1.2509789 1.4713302 0.6004746 0.8650430 0.5825389 #> y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1 x2~~x2 x3~~x3 y1~~y1 #> 1.4402477 2.1829448 0.7115901 0.3627964 1.3717741 -2.5085299 -2.1167106 -0.7621551 0.6176915 #> y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60 #> 2.0256969 1.6005337 1.1707941 0.8385629 1.6030457 1.2697708 1.1962487 -0.8016262 1.3546241 #> dem65~~dem65 x1~1 x2~1 x3~1 y1~1 y2~1 y3~1 y4~1 y5~1 #> -1.8050678 5.0543838 4.7921946 3.5576898 5.4646667 4.2564429 6.5631103 4.4525330 5.1362519 #> y6~1 y7~1 y8~1 #> 2.9780741 6.1962639 4.0433897 opt <- optim(par = par, fn = lessSEM:::.fitFunction, # use the fitting function wrapper gr = lessSEM:::.gradientFunction, # use the gradient function wrapper SEM = mySEM, # use the SEM we created above raw = TRUE, # make sure to tell the functions that we are using raw parameters method = "BFGS" # use the BFGS optimizer ) print(opt$par) #> ind60=~x2 ind60=~x3 a b c dem60~ind60 dem65~ind60 dem65~dem60 y1~~y5 #> 2.1791276 1.8180458 1.1909397 1.1740909 1.2511328 1.4725867 0.6007137 0.8649836 0.5817910 #> y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1 x2~~x2 x3~~x3 y1~~y1 #> 1.4336940 2.1828828 0.7229781 0.3605874 1.3774602 -2.5093044 -2.1126718 -0.7633056 0.6178333 #> y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60 #> 2.0246995 1.6022529 1.1690474 0.8385775 1.6039528 1.2711405 1.1971146 -0.8013687 1.3542390 #> dem65~~dem65 x1~1 x2~1 x3~1 y1~1 y2~1 y3~1 y4~1 y5~1 #> -1.8057174 5.0543838 4.7921946 3.5576898 5.4646667 4.2564429 6.5631103 4.4525330 5.1362519 #> y6~1 y7~1 y8~1 #> 2.9780741 6.1962639 4.0433897 ``` Note that the parameter `a` is now back at the maximum likelihood estimate from before. However, all parameters are still in raw format. To get the transformed parameters, let's take one more step: ```r mySEM <- lessSEM:::.setParameters(SEM = mySEM, # the model labels = names(opt$par), # names of the parameters values = opt$par, # values of the parameters raw = TRUE) print(lessSEM:::.getParameters(mySEM, raw = FALSE)) #> ind60=~x2 ind60=~x3 a b c dem60~ind60 dem65~ind60 dem65~dem60 y1~~y5 #> 2.17912764 1.81804575 1.19093968 1.17409087 1.25113276 1.47258673 0.60071368 0.86498364 0.58179103 #> y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1 x2~~x2 x3~~x3 y1~~y1 #> 1.43369401 2.18288278 0.72297808 0.36058737 1.37746019 0.08132479 0.12091448 0.46612308 1.85490471 #> y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60 #> 7.57383439 4.96420359 3.21892497 2.31307420 4.97264951 3.56491599 3.31055101 0.44871438 3.87381195 #> dem65~~dem65 x1~1 x2~1 x3~1 y1~1 y2~1 y3~1 y4~1 y5~1 #> 0.16435650 5.05438384 4.79219463 3.55768979 5.46466667 4.25644288 6.56311025 4.45253304 5.13625192 #> y6~1 y7~1 y8~1 #> 2.97807408 6.19626389 4.04338968 ``` Compare those to the parameter estimates from **lavaan**: ```r coef(lavaanModel) #> ind60=~x2 ind60=~x3 a b c a b c dem60~ind60 #> 2.180 1.818 1.191 1.175 1.251 1.191 1.175 1.251 1.471 #> dem65~ind60 dem65~dem60 y1~~y5 y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1 #> 0.600 0.865 0.583 1.440 2.183 0.712 0.363 1.372 0.081 #> x2~~x2 x3~~x3 y1~~y1 y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 #> 0.120 0.467 1.855 7.581 4.956 3.225 2.313 4.968 3.560 #> y8~~y8 ind60~~ind60 dem60~~dem60 dem65~~dem65 #> 3.308 0.449 3.875 0.164 ``` Finally, we can compute the standard errors: ```r lessSEM:::.standardErrors(SEM = mySEM, raw = FALSE) #> ind60=~x2 ind60=~x3 a b c dem60~ind60 dem65~ind60 dem65~dem60 y1~~y5 #> 0.13885220 0.15204330 0.14166120 0.11987057 0.12295637 0.39139697 0.23828914 0.07567860 0.36462027 #> y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1 x2~~x2 x3~~x3 y1~~y1 #> 0.68977247 0.73096919 0.62119517 0.46062832 0.57969390 0.01968652 0.06991196 0.08897395 0.45717112 #> y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60 #> 1.34332170 0.96373267 0.74092220 0.48364101 0.89600779 0.73922562 0.71425331 0.08675480 0.88802932 #> dem65~~dem65 x1~1 x2~1 x3~1 y1~1 y2~1 y3~1 y4~1 y5~1 #> 0.23331748 0.08406657 0.17326967 0.16121433 0.29892606 0.43891242 0.39404806 0.37957637 0.30446534 #> y6~1 y7~1 y8~1 #> 0.39247640 0.36442149 0.37545879 ``` Let's compare this to **lavaan** again: ```r parameterEstimates(lavaanModel)[,1:6] #> lhs op rhs label est se #> 1 ind60 =~ x1 1.000 0.000 #> 2 ind60 =~ x2 2.180 0.138 #> 3 ind60 =~ x3 1.818 0.152 #> 4 dem60 =~ y1 1.000 0.000 #> 5 dem60 =~ y2 a 1.191 0.139 #> 6 dem60 =~ y3 b 1.175 0.120 #> 7 dem60 =~ y4 c 1.251 0.117 #> 8 dem65 =~ y5 1.000 0.000 #> 9 dem65 =~ y6 a 1.191 0.139 #> 10 dem65 =~ y7 b 1.175 0.120 #> 11 dem65 =~ y8 c 1.251 0.117 #> 12 dem60 ~ ind60 1.471 0.392 #> 13 dem65 ~ ind60 0.600 0.226 #> 14 dem65 ~ dem60 0.865 0.075 #> 15 y1 ~~ y5 0.583 0.356 #> 16 y2 ~~ y4 1.440 0.689 #> 17 y2 ~~ y6 2.183 0.737 #> 18 y3 ~~ y7 0.712 0.611 #> 19 y4 ~~ y8 0.363 0.444 #> 20 y6 ~~ y8 1.372 0.577 #> 21 x1 ~~ x1 0.081 0.019 #> 22 x2 ~~ x2 0.120 0.070 #> 23 x3 ~~ x3 0.467 0.090 #> 24 y1 ~~ y1 1.855 0.433 #> 25 y2 ~~ y2 7.581 1.366 #> 26 y3 ~~ y3 4.956 0.956 #> 27 y4 ~~ y4 3.225 0.723 #> 28 y5 ~~ y5 2.313 0.479 #> 29 y6 ~~ y6 4.968 0.921 #> 30 y7 ~~ y7 3.560 0.710 #> 31 y8 ~~ y8 3.308 0.704 #> 32 ind60 ~~ ind60 0.449 0.087 #> 33 dem60 ~~ dem60 3.875 0.866 #> 34 dem65 ~~ dem65 0.164 0.227 ``` ## References * Fox, J. (2006). Teacher’s corner: Structural equation modeling with the sem package in R. Structural Equation Modeling: A Multidisciplinary Journal, 13(3), 465–486. https://doi.org/10.1207/s15328007sem1303_7 * McArdle, J. J., & McDonald, R. P. (1984). Some algebraic properties of the Reticular Action Model for moment structures. British Journal of Mathematical and Statistical Psychology, 37(2), 234–251. https://doi.org/10.1111/j.2044-8317.1984.tb00802.x