lavResidualsY (a re-imagining of "casewise" residuals)

[jebyrnes]

Given the existence of lavPredictY() we can easily write a function lavResidualsY() to get residuals of observed variables.

lavResidualsY <- function(object,
                          ynames = lavNames(object, "ov.nox"),
                          xnames = lavNames(object, "ov.x")){
  pred <- lavPredictY(object,
                      ynames = ynames,
                      xnames = xnames) |> 
    as.data.frame()
  
  d <- inspect(object, "data") |> 
    as.data.frame()

  r <- lapply(names(pred), function(x){
    pred[[x]] - d[[x]]
  })

  res <- do.call(cbind, r) |> 
    as.data.frame()
  names(res) <- names(pred)
  
  res
}

mod <- "
rich ~ cover + firesev 
firesev ~ age
"

f <- sem(mod, data = piecewiseSEM::keeley) 

#let's get residuals and then do a test for MV normality
r <- lavResidualsY(f)

mvnormtest::mshapiro.test(t(r))

This should also work with LV models, but might take a bit more - e.g.,

lavResidualsY(fit, ynames = c("y5", "y6", "y7", "y8"),
            xnames = c("x1", "x2", "x3", "y1", "y2", "y3", "y4"))

[jebyrnes]

(and you can un-pipe it if you want. I'm just a pipe fan.)