lavInspect(fit, "gamma") returns error for ML-SEMs
TDJorgensen opened this issue · 2 comments
I'm trying to get an idea about the required structure of @SampleStats@NACOV
for a multilevel SEM, related to the lavaan.mi
package. The unstructured model's NACOV is generally stored in a list, but I am unsure whether the list should be length(ngroups)
or length(nblocks)
. The vcov()
for a multigroup SEM is a single matrix because parameters can be constrained across groups, but in the unstructured h1 model, the group blocks of that matrix are independent (thus, separate in the NACOV
list). I expect same would hold for multilevel SEMs (the level-specific covariance-matrix components are independent), but maybe there is additional complexity.
To see how lavaan
provides it, I tried:
example(Demo.twolevel)
lavInspect(fit, "gamma")
Which returned an error from lav_object_inspect_sampstat_gamma()
:
Error in dimnames(x) <- dn :
length of 'dimnames' [1] not equal to array extent
When lav_object_inspect_sampstat_gamma()
runs lav_object_gamma()
, it yields a single list but with a 42x42 matrix, which I did not expect. In a simpler model without exogenous covariates:
model <- '
level: 1
fw =~ y1 + y2 + y3
level: 2
fb =~ y1 + y2 + y3
'
fit2 <- sem(model, data = Demo.twolevel, cluster = "cluster")
lavaan:::lav_object_gamma(fit2)
The matrix is only 9x9, when I would expect 6 (co)variances per level + 3 means at Level 2 (so 15, or 18 if 3 means are included at Level 1). I think this is leading to the error when lav_object_inspect_sampstat_gamma()
tries to assign names.
for a multilevel SEM, related to the lavaan.mi package
Specifically, the NACOV would be needed for lavResiduals()
. But I notice that lavResiduals(fit2, summary = TRUE)
is not available. So perhaps the reason is related to the structure of Gamma not being as obvious as I expect.
A partial explanation:
- in the single-group, two-level setting, there is only one 'Gamma' matrix
- the 'Gamma' matrix is here defined as A1^{-1} B1 A1^{-1}, where A1 is the expected or observed information, and B1 is the first-order information, in both cases for the unrestricted (h1) model
- the 'order' of the elements is as follows: mu.w within, vech(sigma.w) within, mu.b between, vech(sigma.b) between
- the code to compute those ingredients (A1, B1) can be found in lav_mvnorm_cluster.R
- an internal convenience function to compute this 'Gamma' is lav_model_h1_omega()
- lav_object_inspect_gamma() is not calling this lav_model_h1_omega() function (yet); this is for historical reasons: Originally, 'Gamma' was only supposed to be sample-based, and it was computed using lav_samplestats_gamma(); the latter can handle 'clustered' data, but that Gamma is not the same as the two-level 'Omega' version
- the problem is that the term 'Gamma' is overloaded: it means different things in different contexts, and therefore, it is somewhat unclear what lavInspect(, "gamma") should return in this setting...