Error in reliability()

[perlus]

Hello together, Thank you so much for all your work on semTools and lavaan.

Since I had an issue with runMI when trying to inspect fit scales I switched to the dev Version of semTools, as proposed by Terrence in another post. Then most worked fine. But when trying to get the reliability measures from my large lavaan.mi object (run as seen below) I get an error. When inspecting the lavInspect as well as the lavListInspect functions I indeed cannot find a what option 'est'. Unfortunatelly I am not witty enough in coding to correct the problem myself.

reliability(KFAfin_impRun)
Error in lavListInspect(object = object, what = what, add.labels = add.labels,  : 
  unknown `what' argument in inspect function: `est'

Thanks for your help in advance.

Kind regards from Germany, Alex

[TDJorgensen]

est is only a lavInspect() option for lavaan objects, not lavaanList objects (because there would be a different set of estimates for each data set to which the model was fitted).

The reliability() function is not designed to work with lavaan.mi objects. I'm sure I can allow it to do so, but I can't say when I will be able to make time for that on my semTools to-do list. Currently, you could save reliability estimates from each imputation:

HSMiss <- HolzingerSwineford1939[ , c(paste("x", 1:9, sep = ""),
                                      "ageyr","agemo","school")]
set.seed(12345)
HSMiss$x5 <- ifelse(HSMiss$x5 <= quantile(HSMiss$x5, .3), NA, HSMiss$x5)
age <- HSMiss$ageyr + HSMiss$agemo/12
HSMiss$x9 <- ifelse(age <= quantile(age, .3), NA, HSMiss$x9)

## specify CFA model from lavaan's ?cfa help page
HS.model <- '
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
'

## impute data within runMI...
out1 <- cfa.mi(HS.model, data = HSMiss, m = 3, seed = 12345,
               miArgs = list(noms = "school"), FUN = function(fit) {
                   REL <- reliability(fit)
                   list(alpha = REL["alpha",],
                        omega = REL["omega",])
               })
out1@funList

Then take the average to naively pool them:

## stack imputations in rows
(omegas <- do.call(rbind, lapply(out1@funList, function(x) x$omega)))
## means
colMeans(omegas)

I'll keep this issue open until I add some syntax to the reliability() source code to grab pooled est output from lavaan.mi objects.

[perlus]

Hi Terrence, thank you so much for your quick and helpful answer.

I will also need the AVE from the original reliabilities function. Unfortunately

REL <- reliability(fit) list(alpha = REL["alpha",], omega = REL["omega",], ave = REL["ave",])

provoked the error "Index out of bounds"...

[TDJorgensen]

Because the row is labeled "avevar", not "ave"

[perlus]

Thank you so much for this quick answer. Safed my day!

[perlus]

Last question: Is it possible to also retrieve the rsquare from the inspect function that way?

I tried:

{
  INS <- inspect(fit)
  list(rsquare = INS["rsquare",])
}

Which throws the error wrong number of Dimensions.

[TDJorgensen]

Yes, you can add any output you want to the list. You get an error because you need to ask inspect() for "rsquare" rather than the default "free" output.

out1 <- cfa.mi(HS.model, data = HSMiss, m = 3, seed = 12345,
               miArgs = list(noms = "school"), FUN = function(fit) {
                   REL <- reliability(fit)
                   list(alpha = REL["alpha",],
                        omega = REL["omega",],
                        rsq = lavInspect(fit, "rsquare"))
               })

[perlus]

Thank you. Terrence, you cannot imagine how much this helped.

But I need to come back on another possible related issue:

In a very small sem model (three factors, nine indicators) I get a rather high AIC (>20000), whilst all other fit measures seem reasonable and very good.

[TDJorgensen]

AIC is calculated from the -2 log-likelihood, so (holding model fit constant) it increases with sample size, not model size. AIC is useful only for model comparison (lowest is preferred), not for judging a model in isolation.

[perlus]

Thank you, yes I want to compare several models. Just was surprised because another model showed AIC of around 700 or so, with n = 250, in AMOS though. Here estimated a model with n = 850 in lavaan wit runMI.

[TDJorgensen]

reliability(), reliabilityL2(), and maximalRelia() now all accept lavaan.mi objects, and will calculate reliability using pooled estimates. The development version is currently needed for these features:

devtools::install_github("simsem/semTools/semTools")

I expect to submit a semTools update to CRAN by the end of the summer.

[TDJorgensen]

FYI, not sure why it didn't occur to me at the time you asked about rsquare, but you can already obtain that from pooled estimates using summary(fit, rsquare = TRUE).