R/lav_cfa_utils.R

In lavaan: Latent Variable Analysis

# utility functions needed for lav_cfa_*

# compute THETA and PSI, given lambda using either ULS or GLS
# this function assumes:
# - THETA is diagonal
# - PSI is unrestricted
# - we assume W = S^{-1}
#
# YR 17 oct 2022: - add lower/upper bounds for theta (only to compute PSI)
#                 - use 'lambda' correction to ensure PSI is positive definite
# YR 02 feb 2023: - add psi.mapping.ML argument
lav_cfa_lambda2thetapsi <- function(lambda = NULL, S = NULL, S.inv = NULL,
                                    GLS = FALSE, psi.mapping.ML = FALSE,
                                    nobs = 20L) {
    LAMBDA <- as.matrix(lambda)
    nvar <- nrow(LAMBDA); nfac <- ncol(LAMBDA)

    if(GLS) {
        # see Browne, 1974 section 4 case II
        if(is.null(S.inv)) {
            W <- solve(S)
        } else {
            W <- S.inv
        }
        tLW <- crossprod(LAMBDA, W)
        M <- solve(tLW %*% LAMBDA, tLW) # GLS mapping
        #D <- W %*% LAMBDA %*% M # symmmetric
        D <- crossprod(M, tLW)
        #theta <- solve(W*W - D*D, diag(W %*% S %*% W - D %*% S %*% D))
        theta  <- try(solve(W*W - D*D, diag(W - D)), # because W == S^{-1}
                      silent = TRUE)
        if(inherits(theta, "try-error")) {
            # what to do?
            warning("lavaan WARNING: problem computing THETA values; trying pace algorithm")
            theta <- lav_efa_pace(S = S, nfactors = nfac, theta.only = TRUE)
        }
    } else {
        # see Hagglund 1982, section 4
        M <- solve(crossprod(LAMBDA), t(LAMBDA)) # ULS mapping function
        D <- LAMBDA %*% M
        theta <-  try(solve(diag(nvar) - D*D, diag(S - (D %*% S %*% D))),
                      silent = TRUE)
        if(inherits(theta, "try-error")) {
            # what to do?
            warning("lavaan WARNING: problem computing THETA values; trying pace algorithm")
            theta <- lav_efa_pace(S = S, nfactors = nfac, theta.only = TRUE)
        }
    }
    theta.nobounds <- theta

    # ALWAYS check bounds for theta (only to to compute PSI)!
    theta.bounds <- TRUE
    if(theta.bounds) {
        diagS <- diag(S)
        # lower bound
        lower.bound <- diagS * 0 # * 0.01
        too.small.idx <- which(theta < lower.bound)
        if(length(too.small.idx) > 0L) {
            theta[ too.small.idx ] <- lower.bound[ too.small.idx ]
        }

        # upper bound
        upper.bound <- diagS * 1 # * 0.99
        too.large.idx <- which(theta > upper.bound)
        if(length(too.large.idx) > 0L) {
            theta[ too.large.idx ] <- upper.bound[ too.large.idx ]
        }
    }

    # psi
    diag.theta <- diag(theta, nvar)
    lambda <- try(lav_matrix_symmetric_diff_smallest_root(S, diag.theta),
                  silent = TRUE)
    if(inherits(lambda, "try-error")) {
        warning("lavaan WARNING: failed to compute lambda")
        SminTheta <- S - diag.theta # and hope for the best
    } else {
        cutoff <- 1 + 1/(nobs - 1)
        if(lambda < cutoff) {
            lambda.star <- lambda - 1/(nobs - 1)
            SminTheta <- S - lambda.star * diag.theta
        } else {
            SminTheta <- S - diag.theta
        }
    }

    # just like local SAM
    if(psi.mapping.ML) {
        Ti <- 1/theta
        zero.theta.idx <- which(abs(theta) < 0.01) # be conservative
        if(length(zero.theta.idx) > 0L) {
            Ti[zero.theta.idx] <- 1
        }
        M <- solve(t(LAMBDA) %*% diag(Ti, nvar) %*% LAMBDA) %*% t(LAMBDA) %*% diag(Ti, nvar)
        PSI <- M %*% SminTheta %*% t(M) # ML
    } else {
        PSI <- M %*% SminTheta %*% t(M) # ULS/GLS
    }

    # we take care of the bounds later!
    list(lambda = LAMBDA, theta = theta.nobounds, psi = PSI)
}

# compute PSI, given lambda and theta using either ULS, GLS, ML
# this function assumes:
# - THETA is diagonal
# - PSI is unrestricted
#
# YR 08 Mar 2023: - first version
lav_cfa_lambdatheta2psi <- function(lambda = NULL, theta = NULL, # vector!
                                    S = NULL, S.inv = NULL,
                                    mapping = "ML", nobs = 20L) {
    LAMBDA <- as.matrix(lambda)
    nvar <- nrow(LAMBDA); nfac <- ncol(LAMBDA)

    theta.nobounds <- theta

    # ALWAYS check bounds for theta to compute PSI
    diagS <- diag(S)
    # lower bound
    lower.bound <- diagS * 0 # * 0.01
    too.small.idx <- which(theta < lower.bound)
    if(length(too.small.idx) > 0L) {
        theta[ too.small.idx ] <- lower.bound[ too.small.idx ]
    }

    # upper bound
    upper.bound <- diagS * 1 # * 0.99
    too.large.idx <- which(theta > upper.bound)
    if(length(too.large.idx) > 0L) {
        theta[ too.large.idx ] <- upper.bound[ too.large.idx ]
    }

    # psi
    diag.theta <- diag(theta, nvar)
    lambda <- try(lav_matrix_symmetric_diff_smallest_root(S, diag.theta),
                  silent = TRUE)
    if(inherits(lambda, "try-error")) {
        warning("lavaan WARNING: failed to compute lambda")
        SminTheta <- S - diag.theta # and hope for the best
    } else {
        cutoff <- 1 + 1/(nobs - 1)
        if(lambda < cutoff) {
            lambda.star <- lambda - 1/(nobs - 1)
            SminTheta <- S - lambda.star * diag.theta
        } else {
            SminTheta <- S - diag.theta
        }
    }

    # mapping matrix
    if(mapping == "ML") {
        Ti <- 1/theta
        zero.theta.idx <- which(abs(theta) < 0.01) # be conservative
        if(length(zero.theta.idx) > 0L) {
            Ti[zero.theta.idx] <- 1
        }
        M <- solve(t(LAMBDA) %*% diag(Ti, nvar) %*% LAMBDA) %*% t(LAMBDA) %*% diag(Ti, nvar)
    } else if(mapping == "GLS") {
        if(is.null(S.inv)) {
            S.inv <- try(solve(S), silent = TRUE)
        }
        if(inherits(S.inv, "try-error")) {
            M <- tcrossprod(solve(crossprod(LAMBDA)), LAMBDA)
        } else {
            M <- solve(t(LAMBDA) %*% S.inv %*% LAMBDA) %*% t(LAMBDA) %*% S.inv
        }
    } else if(mapping == "ULS") {
        M <- tcrossprod(solve(crossprod(LAMBDA)), LAMBDA)
    }

    # compute PSI
    PSI <- M %*% SminTheta %*% t(M)

    PSI
}

# compute theta elements for a 1-factor model
lav_cfa_theta_spearman <- function(S, bounds = "wide") {
    p <- ncol(S); out <- numeric(p)
    R <- cov2cor(S)
    for(p.idx in seq_len(p)) {
        var.p <- R[p.idx, p.idx]
        x <- R[,p.idx][-p.idx]
        aa <- lav_matrix_vech(tcrossprod(x), diagonal = FALSE)
        ss <- lav_matrix_vech(R[-p.idx, -p.idx, drop = FALSE], diagonal = FALSE)
        h2 <- mean(aa/ss) # communaliteit
        if(bounds == "standard") {
            h2[h2 < 0] <- 0
            h2[h2 > 1] <- 1
        } else if(bounds == "wide") {
            h2[h2 < -0.05] <- -0.05 # correponds to lower bound ov.var "wide"
            h2[h2 > +1.20] <- +1.20 # correponds to upper bound ov.var "wide"
        }
        out[p.idx] <- (1 - h2) * S[p.idx, p.idx]
    }
    out
}