yrosseel/lavaan

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...