R/inspect_da.R

In modsem: Latent Interaction (and Moderation) Analysis in Structural Equation Models (SEM)

inspectDA_Matrices <- c("lambda", "theta", "wmat", "tmat",
                        "psi", "beta", "nu", "alpha")


inspectDA_Optim <- c("coefficients.free", "vcov.free", "information",
                     "loglik", "iterations", "convergence")


modsem_inspect_da <- function(model,
                              what = "default",
                              standardized = FALSE,
                              is.public = TRUE) {
  # for brevity...
  .modsemVector <- \(...) modsemVector(..., is.public = is.public)
  .modsemMatrix <- \(...) modsemMatrix(..., is.public = is.public)

  mod_stopif(!length(what), "`what` is of length zero!")

  if (standardized)
    model <- standardize_model(model)

  finalModel <- model$model
  groupModels <- finalModel$models
  ovs      <- c(finalModel$info$allIndsXis, finalModel$info$allIndsEtas)
  n.groups <- length(groupModels)
  is.multi <- n.groups > 1L

  group.names <- names(groupModels)
  if (!length(group.names) || any(!nzchar(group.names))) {
    level.names <- finalModel$info$group.levels

    if (!is.null(level.names) && length(level.names) == n.groups)
      group.names <- level.names
    else
      group.names <- paste0("Group", seq_len(n.groups))
  }

  expected.raw <- model$expected.matrices
  expected.by.group <- vector("list", length = n.groups)
  names(expected.by.group) <- group.names

  if (is.list(expected.raw) && length(expected.raw)) {
    names.raw <- names(expected.raw)
    for (g in seq_len(n.groups)) {
      idx <- integer()
      if (length(names.raw)) {
        candidates <- unique(c(as.character(g), group.names[[g]]))
        idx <- match(candidates, names.raw, nomatch = NA_integer_)
        idx <- idx[!is.na(idx)]
      }
      if (!length(idx) && g <= length(expected.raw)) {
        idx <- g
      }
      if (length(idx)) expected.by.group[[g]] <- expected.raw[[idx[[1]]]]
    }
  } else if (!is.null(expected.raw)) {
    expected.by.group <- replicate(n.groups, expected.raw, simplify = FALSE)
    names(expected.by.group) <- group.names
  }

  buildExpectedPayload <- function(expected) {
    if (is.null(expected)) {
      return(list(
        cov.ov  = NULL,
        cov.lv  = NULL,
        cov.all = NULL,
        cor.ov  = NULL,
        cor.lv  = NULL,
        cor.all = NULL,
        mean.lv = NULL,
        mean.ov = NULL,
        mean.all = NULL,
        r2.all  = NULL,
        r2.lv   = NULL,
        r2.ov   = NULL,
        res.all = NULL,
        res.lv  = NULL,
        res.ov  = NULL
      ))
    }

    cov.ov  <- expected$sigma.ov
    cov.lv  <- expected$sigma.lv
    cov.all <- expected$sigma.all

    cor.ov  <- if (!is.null(cov.ov))  cov2cor(cov.ov)   else NULL
    cor.lv  <- if (!is.null(cov.lv))  cov2cor(cov.lv)   else NULL
    cor.all <- if (!is.null(cov.all)) cov2cor(cov.all)  else NULL

    list(
      cov.ov  = .modsemMatrix(cov.ov,  symmetric = TRUE),
      cov.lv  = .modsemMatrix(cov.lv,  symmetric = TRUE),
      cov.all = .modsemMatrix(cov.all, symmetric = TRUE),

      cor.ov  = .modsemMatrix(cor.ov,  symmetric = TRUE),
      cor.lv  = .modsemMatrix(cor.lv,  symmetric = TRUE),
      cor.all = .modsemMatrix(cor.all, symmetric = TRUE),

      mean.lv  = .modsemMatrix(expected$mu.lv),
      mean.ov  = .modsemMatrix(expected$mu.ov),
      mean.all = .modsemMatrix(expected$mu.all),

      r2.all  = .modsemVector(expected$r2.all),
      r2.lv   = .modsemVector(expected$r2.lv),
      r2.ov   = .modsemVector(expected$r2.ov),

      res.all = .modsemVector(expected$res.all),
      res.lv  = .modsemVector(expected$res.lv),
      res.ov  = .modsemVector(expected$res.ov)
    )
  }

  buildGroupPayload <- function(submodel, expected) {
    # matrices in lavaan notation
    matrices <- submodel$expected.matrices$matrices

    lambda <- matrices$lambda
    theta  <- matrices$theta
    wmat   <- matrices$wmat
    tmat   <- matrices$tmat
    beta   <- matrices$beta
    psi    <- matrices$psi
    alpha  <- matrices$alpha
    nu     <- matrices$nu

    if (!NROW(wmat)) {
      wmat <- lambda
      wmat[TRUE] <- 0
    }

    if (!NROW(tmat)) {
      tmat <- theta 
      tmat[TRUE] <- 0
    }

    if (!is.null(nu))    colnames(nu)   <- "~1"
    if (!is.null(alpha)) colnames(alpha) <- "~1"

    out <- c(
      list(
        N      = submodel$data$n,
        data   = submodel$data$data.full,
        lambda = .modsemMatrix(lambda),
        wmat   = .modsemMatrix(wmat),
        theta  = .modsemMatrix(theta, symmetric = TRUE),
        tmat   = .modsemMatrix(tmat),
        nu     = .modsemMatrix(nu),
        alpha  = .modsemMatrix(alpha),
        beta   = .modsemMatrix(beta),
        psi    = .modsemMatrix(psi, symmetric = TRUE)
      ),
      buildExpectedPayload(expected)
    )
  }

  group.payloads <- Map(buildGroupPayload, groupModels, expected.by.group)
  names(group.payloads) <- group.names

  collapseField <- function(field) {
    values <- lapply(group.payloads, `[[`, field)

    if (!is.multi)
      return(values[[1]])

    names(values) <- group.names
    values
  }

  N.val            <- collapseField("N")
  data.val         <- collapseField("data")
  lambda.val       <- collapseField("lambda")
  wmat.val         <- collapseField("wmat")
  tmat.val         <- collapseField("tmat")
  nu.val           <- collapseField("nu")
  theta.val        <- collapseField("theta")
  beta.val         <- collapseField("beta")
  psi.val          <- collapseField("psi")
  alpha.val        <- collapseField("alpha")
  cov.ov.val       <- collapseField("cov.ov")
  cov.lv.val       <- collapseField("cov.lv")
  cov.all.val      <- collapseField("cov.all")
  cor.ov.val       <- collapseField("cor.ov")
  cor.lv.val       <- collapseField("cor.lv")
  cor.all.val      <- collapseField("cor.all")
  mean.lv.val      <- collapseField("mean.lv")
  mean.ov.val      <- collapseField("mean.ov")
  mean.all.val     <- collapseField("mean.all")
  r2.all.val       <- collapseField("r2.all")
  r2.lv.val        <- collapseField("r2.lv")
  r2.ov.val        <- collapseField("r2.ov")
  res.all.val      <- collapseField("res.all")
  res.lv.val       <- collapseField("res.lv")
  res.ov.val       <- collapseField("res.ov")

  info <- list(
     N                 = N.val,
     vcov.all          = .modsemMatrix(model$vcov.all, symmetric = TRUE),
     vcov.free         = .modsemMatrix(model$vcov.free, symmetric = TRUE),
     information       = .modsemMatrix(model$FIM, symmetric = TRUE),
     data              = data.val,
     coefficients.all  = .modsemVector(model$coefs.all),
     coefficients.free = .modsemVector(model$coefs.free),
     partable          = parameter_estimates(model, is.public = is.public),
     partable.input    = model$originalParTable,
     loglik            = model$logLik,
     iterations        = model$iterations,
     convergence       = model$convergence,
     ovs               = ovs,

     ngroups     = model$model$info$n.groups,
     group       = model$args$group,
     group.label = if (!is.null(model$args$group)) group.names else NULL,

     lambda  = lambda.val,
     theta   = theta.val,
     wmat    = wmat.val,
     tmat    = tmat.val,
     nu      = nu.val,
     alpha   = alpha.val,
     beta    = beta.val,
     psi     = psi.val,

     cov.ov  = cov.ov.val,
     cov.lv  = cov.lv.val,
     cov.all = cov.all.val,

     cor.ov  = cor.ov.val,
     cor.lv  = cor.lv.val,
     cor.all = cor.all.val,

     mean.lv  = mean.lv.val,
     mean.ov  = mean.ov.val,
     mean.all = mean.all.val,

     r2.all  = r2.all.val,
     r2.lv   = r2.lv.val,
     r2.ov   = r2.ov.val,

     res.all = res.all.val,
     res.lv  = res.lv.val,
     res.ov  = res.ov.val
  )

  FIT <- \() {
    h0 <- estimate_h0(model, calc.se = FALSE)
    list(
      fit.h0 = fit_modsem_da(h0, chisq = TRUE),
      fit.h1 = fit_modsem_da(model, chisq = FALSE),
      comparative.fit = compare_fit(est_h1 = model, est_h0 = h0)
    )
  }

  if (length(what) > 1) {
    fields <- info[what]

  } else {
    fields <- switch(
      EXPR      = what,
      coef.all  = info[c("vcov.all", "coefficients.all")],
      coef      = info[c("vcov.all", "coefficients.all")],
      coef.free = info[c("vcov.free", "coefficients.free")],
      default   = info[names(info) != "data"],
      all       = info,
      matrices  = info[inspectDA_Matrices],
      optim     = info[inspectDA_Optim],
      fit       = FIT(),
      info[[what]]
    )
  }

  nullvalues <- vapply(fields, FUN.VALUE = logical(1L), FUN = is.null)
  okifnull   <- names(fields) %in% c("group", "group.label")

  mod_warnif(any(nullvalues & !okifnull),
         "Some fields in `modsem_inspect()` could not be retrieved!")

  fields
}