For didactic purposes, I am simulating data consistent with a fixed effects model and then running the model in SEM using lavaan. I don't understand what is going on with the underlying matrices lavaan is giving me using the inspect() function.
Here is my covariance matrix for reproducibility purposes:
obs <- matrix( c( 2.1531960, 0.9219589, 1.1247727, 2.1399405, 1.6009631, 1.7877342,
0.9219589, 2.0183384, 1.1263892, 1.6586093, 2.0530063, 1.7664361,
1.1247727, 1.1263892, 1.9152931, 1.8183510, 1.7904083, 1.9676634,
2.1399405, 1.6586093, 1.8183510, 4.2225644, 3.0380404, 3.1018654,
1.6009631, 2.0530063, 1.7904083, 3.0380404, 4.0023862, 2.9819462,
1.7877342, 1.7664361, 1.9676634, 3.1018654, 2.9819462, 3.8469132),
nrow = 6, ncol = 6)
and here is my lavaan code:
fe_sem <- '
# Define individual effects variable
n =~ 1*y1 + 1*y2 + 1*y3
# Regressions, constrain coefficient to be equal over time
y1 ~ b*x1
y2 ~ b*x2
y3 ~ b*x3
# Allow unrestricted correlation between eta and covariates
n ~~ x1 + x2 + x3
x1 ~~ x2 + x3
x2 ~~ x3
# Constrain residual variances to be equal over time
y1 ~~ e*y1
y2 ~~ e*y2
y3 ~~ e*y3
'
fe_sem.fit <- sem( model = fe_sem, sample.cov = obs, sample.nobs = 300)
Everything works fine. The estimated coefficient is correct given how I simulated the data. What I don't understand is how lavaan is getting the answer. When I use
inspect( fe_sem.fit)
the lambda matrix is a matrix of zeros:
## $lambda
## n y1 y2 y3 x1 x2 x3
## y1 0 0 0 0 0 0 0
## y2 0 0 0 0 0 0 0
## y3 0 0 0 0 0 0 0
## x1 0 0 0 0 0 0 0
## x2 0 0 0 0 0 0 0
## x3 0 0 0 0 0 0 0
But I define the latent individual effects variable as n =~ 1*y1 + 1*y2 + 1*y3, so as far as I know, the lambda matrix should look like:
## $lambda
## n y1 y2 y3 x1 x2 x3
## y1 1 1 0 0 0 0 0
## y2 1 0 1 0 0 0 0
## y3 1 0 0 1 0 0 0
## x1 0 0 0 0 1 0 0
## x2 0 0 0 0 0 1 0
## x3 0 0 0 0 0 0 1
I have double-checked this by hand. Having lambda as a matrix of zeros messes everything up. Does anyone know what is going on here? Is there something I don't understand about the matrix notation used by lavaan?
Note that inspect is an older function and now used as a shortcut to lavInspect. Further note the dimnames for the matrix.
library(lavaan)
obs <- matrix( c( 2.1531960, 0.9219589, 1.1247727, 2.1399405, 1.6009631, 1.7877342,
0.9219589, 2.0183384, 1.1263892, 1.6586093, 2.0530063, 1.7664361,
1.1247727, 1.1263892, 1.9152931, 1.8183510, 1.7904083, 1.9676634,
2.1399405, 1.6586093, 1.8183510, 4.2225644, 3.0380404, 3.1018654,
1.6009631, 2.0530063, 1.7904083, 3.0380404, 4.0023862, 2.9819462,
1.7877342, 1.7664361, 1.9676634, 3.1018654, 2.9819462, 3.8469132),
nrow = 6, ncol = 6, dimnames = list(NULL, c(paste0(rep(c("x", "y"), each = 3), 1:3)))
)
The output of lavInspect can be regulated with the argument what which defaults to "free". That is, by default, lavInspect gives you the free parameters -- since the loadings are constrained to 1, there are no free parameters. Using the most basic version of your MWE you'd get
fe_sem <- "
# Define individual effects variable
n =~ 1*y1 + 1*y2 + 1*y3
"
fe_sem.fit <- sem( model = fe_sem, sample.cov = obs, sample.nobs = 300)
lavInspect(fe_sem.fit)$lambda
resulting in
## n
## y1 0
## y2 0
## y3 0
You can access the estimated parameters using what = "est".
lavInspect(fe_sem.fit, what = "est")$lambda
## n
## y1 1
## y2 1
## y3 1
You can cross check that behaviour using n =~ y1 + y2 + y3. Finally, with the introduction of the regression statements the y variables became endogenous rather than measurement loadings. That is, for "loadings" on n you'd need to look at the beta matrix.
fe_sem <- "
# Define individual effects variable
n =~ 1*y1 + 1*y2 + 1*y3
# Regressions, constrain coefficient to be equal over time
y1 ~ b*x1
y2 ~ b*x2
y3 ~ b*x3
# Allow unrestricted correlation between eta and covariates
n ~~ x1 + x2 + x3
x1 ~~ x2 + x3
x2 ~~ x3
# Constrain residual variances to be equal over time
y1 ~~ e*y1
y2 ~~ e*y2
y3 ~~ e*y3
"
fe_sem.fit <- sem( model = fe_sem, sample.cov = obs, sample.nobs = 300)
lavInspect(fe_sem.fit, what = "est")$lambda
giving you
## n y1 y2 y3 x1 x2 x3
## y1 0 1 0 0 0 0 0
## y2 0 0 1 0 0 0 0
## y3 0 0 0 1 0 0 0
## x1 0 0 0 0 1 0 0
## x2 0 0 0 0 0 1 0
## x3 0 0 0 0 0 0 1
and
lavInspect(fe_sem.fit, what = "est")$beta
giving you
## n y1 y2 y3 x1 x2 x3
## n 0 0 0 0 0.000 0.000 0.000
## y1 1 0 0 0 0.326 0.000 0.000
## y2 1 0 0 0 0.000 0.326 0.000
## y3 1 0 0 0 0.000 0.000 0.326
## x1 0 0 0 0 0.000 0.000 0.000
## x2 0 0 0 0 0.000 0.000 0.000
## x3 0 0 0 0 0.000 0.000 0.000