yrosseel / lavaan

an R package for structural equation modeling and more
http://lavaan.org
429 stars 98 forks source link

Add h1.model= argument to fitMeasures() #328

Closed TDJorgensen closed 5 months ago

TDJorgensen commented 5 months ago

lav_fit_measures() also checks fitted model for @external$h1.model, as does lav_object_summary(). Here is an example to demonstrate:

## simulate data from 2 x 2 design with orthogonal factors
des <- expand.grid(therapy = c(-1,1), drug = c(-1,1))
des$int <- des$therapy*des$drug
set.seed(123)
dat <- do.call(rbind, lapply(1:5, function(n) {
  ## n=5 per group (4 groups)
  yHat <- 1 + 2*des$therapy + 3*des$drug - 4*des$int
  des$DV <- rnorm(4, mean = yHat, sd = sqrt(4*var(yHat))) # Rsq = 20%
  des
}))

## Fit in lavaan with orthogonality constraints
mod <- ' DV ~ therapy + drug + int
## only estimate variances
  DV ~~ DV
  therapy ~~ therapy
  drug ~~ drug
  int ~~ int
'
fit <- lavaan(mod, data = dat, fixed.x = FALSE,# h1 = FALSE,
              test = c("browne.residual.nt","satorra.bentler","scaled.shifted"))
fit # unscaled ChiSq=0 with df=3

## These 3 df will make any more-restricted model seem to fit better
mod0 <- ' DV ~ b1*therapy + b1*drug # equal effects
          DV ~ 0*int                # no interaction
## only estimate variances
  DV ~~ DV
  therapy ~~ therapy
  drug ~~ drug
  int ~~ int
'
fit0 <- lavaan(mod0, data = dat, fixed.x = FALSE, #h1 = FALSE,
               test = c("browne.residual.nt","satorra.bentler","scaled.shifted"))
fit0 # unscaled ChiSq=6.7 with df=5, NOT significant

## compare to a more appropriate h1 model:
lavTestLRT(fit0, fit) # same ChiSq=6.7, but with df=2, IS significant
## we should prefer this as the default h1 model

## this updated allows us to get fit indices using df=2 rather than 5
fitMeasures(fit0, h1.model = fit, output = "text")

## Optionally override default saturated model in @h1 slot
fit@h1  <- list()
fit0@h1 <- list()
fit@external$h1.model <- fit
fit0@external$h1.model <- fit

## because show() calls lav_fit_summary(), these all display df=2
fit0
lavTestLRT(fit0)
fitMeasures(fit0, output = "text")

Tests with unadjusted df=5 are still available from lavTest(fit0, output = "text")

Contrary to what I speculated when we discussed this in person, I no longer expect a problem with RMRs. The issue here is that there are structural constraints in the unrestricted model's implied moments, which should be reflected by df > 0. But the unrestricted moments are the same, so residuals will be the same. Their _SE_s and CIs in lavResiduals() might be improved, if the custom h1.model parameters are estimated more efficiently. But that's for future research.

TDJorgensen commented 5 months ago

Oh, and the first 5 commits in this PR address some other issues we discussed (e.g., new lavTestLRT() options).

TDJorgensen commented 5 months ago

I just updated the documentation, too.

yrosseel commented 5 months ago

This works really well. Tests did not reveal any issues. Thanks a lot for this PR. This was a lot of work.