R/lav_model_properties.R

In lavaan: Latent Variable Analysis

# collect information about the model that we can use
# (eg. is theta diagonal or not, is the structurual model recursive or not,
#  is the model just a regression model, etc)
#
# initial version: YR 15 March 2021
# - YR 05 Oct 2021: use det(I - B) to check if B is acyclic
# - YR 11 Nov 2021: if no latents, and conditional.x = TRUE, we may have no
#                   beta matrix


# note: there is no 'lavmodel' yet, because we call this in lav_model.R
lav_model_properties <- function(GLIST, lavpartable = NULL, lavpta = NULL,
                                 nmat = NULL, m.free.idx = NULL) {

    if(is.null(lavpta)) {
        lavpta <- lav_partable_attributes(lavpartable)
    }

    nblocks <- lavpta$nblocks

    # is the model a univariate/multivariate linear multiple regression
    # model (per block)?
    uvreg <- logical(nblocks)
    uvord <- logical(nblocks)
    mvreg <- logical(nblocks)
    acyclic <- rep(as.logical(NA), nblocks)
    bowfree <- rep(as.logical(NA), nblocks)
    nexo <- integer(nblocks)

    for(g in seq_len(nblocks)) {

        # at least 1 regression
        if(length(lavpta$vnames$eqs.y[[g]]) == 0L) {
            next
        }

        # find beta index for this block
        mm.in.block <- 1:nmat[g] + cumsum(c(0L,nmat))[g]
        MLIST <- GLIST[mm.in.block]
        beta.idx <- which(names(MLIST) == "beta") + cumsum(c(0L,nmat))[g]
         psi.idx <- which(names(MLIST) == "psi" ) + cumsum(c(0L,nmat))[g]

        if(length(beta.idx) > 0L) {
            # 1. acyclic?
            B <- GLIST[[beta.idx]]
            # keep fixed values (if any); fill in 1 in all 'free' positions
            B[ m.free.idx[[beta.idx]] ] <- 1
            IminB <- diag(nrow(B)) - B
            # if B is acyclic, we should be able to permute the rows/cols of B
            # so that B is upper/lower triangular, and so det(I-B) = 1
            if(det(IminB) == 1) {
                acyclic[g] <- TRUE
            } else {
                acyclic[g] <- FALSE
            }

            # 2. bow-free?
            B.one <- as.integer(B != 0)
            Psi <- GLIST[[psi.idx]]
            # keep fixed values (if any); fill in 1 in all 'free' positions
            Psi[ m.free.idx[[psi.idx]] ] <- 1
            Psi.one <- as.integer(Psi != 0)
            Both.one <- B.one + Psi.one
            if(any(Both.one > 1)) {
                bowfree[g] <- FALSE
            } else {
                bowfree[g] <- TRUE
            }
        } else {
            # perhaps conditional.x = TRUE?
            # if there is no BETA, then we only have Gamma, and the
            # system must be acyclic
            acyclic[g] <- TRUE
            # and also bowfree
            bowfree[g] <- TRUE
        }


        # no latent variables, at least 1 dependent variable
        if(lavpta$nfac[[g]] > 0L) {
            next
        }

        # no mediators
        if(length(lavpta$vnames$eqs.y[[g]]) !=
           length(lavpta$vnames$ov.y[[g]])) {
            next
        }

        # categorical y?
        if(length(lavpta$vnames$ov.ord[[g]]) > 0L) {
            # we only flag the univariate version
            if(length(lavpta$vnames$ov.ord[[g]]) == 1L &&
               length(lavpta$vnames$ov.y[[g]])   == 1L &&
               lavpta$vnames$ov.ord[[g]][1] == lavpta$vnames$ov.y[[g]][1]) {
                uvord[g] <- TRUE
            }

        # mvreg?
        } else {
            if(length(lavpta$vnames$ov.y[[g]]) > 1L) {
                mvreg[g] <- TRUE
            } else {
                uvreg[g] <- TRUE
            }
        }

        nexo[g] <- length(lavpta$vnames$eqs.x[[g]])

    } # g

    modprop <- list( uvreg = uvreg, uvord = uvord, mvreg = mvreg,
                     nexo = nexo, acyclic = acyclic, bowfree = bowfree )

    modprop
}