R/lav_samplestats_igamma.R

In lavaan: Latent Variable Analysis

# YR 18 Dec 2015
# - functions to (directly) compute the inverse of 'Gamma' (the asymptotic
#   variance matrix of the sample statistics)
# - often used as 'WLS.V' (the weight matrix in WLS estimation)
#   and when computing the expected information matrix

# NOTE:
#  - three types:
#       1) plain         (conditional.x = FALSE, fixed.x = FALSE)
#       2) fixed.x       (conditional.x = FALSE, fixed.x = TRUE)
#       3) conditional.x (conditional.x = TRUE)
#  - if conditional.x = TRUE, we ignore fixed.x (can be TRUE or FALSE)

# NORMAL-THEORY
lav_samplestats_Gamma_inverse_NT <- function(Y              = NULL,
                                             COV            = NULL,
                                             ICOV           = NULL,
                                             MEAN           = NULL,
                                             rescale        = TRUE,
                                             x.idx          = integer(0L),
                                             fixed.x        = FALSE,
                                             conditional.x  = FALSE,
                                             meanstructure  = FALSE,
                                             slopestructure = FALSE) {

    # check arguments
    if(length(x.idx) == 0L) {
        conditional.x <- FALSE
        fixed.x       <- FALSE
    }

    if(is.null(ICOV)) {
        if(is.null(COV)) {
            stopifnot(!is.null(Y))

            # coerce to matrix
            Y <- unname(as.matrix(Y)); N <- nrow(Y)
            COV <- cov(Y)
            if(rescale) {
                COV <- COV * (N-1) / N # ML version
            }
        }

        ICOV <- solve(COV)
    }

    # if conditional.x, we may also need COV and MEAN
    if(conditional.x && length(x.idx) > 0L &&
       (meanstructure || slopestructure)) {

        if(is.null(COV)) {
            stopifnot(!is.null(Y))

            # coerce to matrix
            Y <- unname(as.matrix(Y)); N <- nrow(Y)
            COV <- cov(Y)

            if(rescale) {
                COV <- COV * (N-1) / N # ML version
            }
        }

        if(is.null(MEAN)) {
            stopifnot(!is.null(Y))
            MEAN <- unname(colMeans(Y))
        }
    }

    # rename
    S.inv <- ICOV
    S     <- COV
    M     <- MEAN

    # unconditional
    if(!conditional.x) {

        # unconditional - stochastic x
        if(!fixed.x) {
            Gamma.inv <- 0.5*lav_matrix_duplication_pre_post(S.inv %x% S.inv)
            if(meanstructure) {
                Gamma.inv <- lav_matrix_bdiag(S.inv, Gamma.inv)
            }

        # unconditional - fixed x
        } else {
            # handle fixed.x = TRUE
            Gamma.inv <- 0.5*lav_matrix_duplication_pre_post(S.inv %x% S.inv)

            # zero rows/cols corresponding with x/x combinations
            nvar <- NROW(ICOV); pstar <- nvar*(nvar+1)/2
            M <- matrix(0, nvar, nvar)
            M[ lav_matrix_vech_idx(nvar) ] <- seq_len(pstar)
            zero.idx <- lav_matrix_vech(M[x.idx, x.idx, drop = FALSE])
            Gamma.inv[zero.idx,] <- 0
            Gamma.inv[,zero.idx] <- 0

            if(meanstructure) {
                S.inv.nox <- S.inv
                S.inv.nox[x.idx,]  <- 0; S.inv.nox[,x.idx] <- 0
                Gamma.inv <- lav_matrix_bdiag(S.inv.nox, Gamma.inv)
            }
        }

    } else {
        # conditional.x

        # 4 possibilities:
        # - no meanstructure, no slopes
        # -    meanstructure, no slopes
        # - no meanstructure, slopes
        # -    meanstructure, slopes

        S11 <- S.inv[-x.idx, -x.idx, drop = FALSE]

        Gamma.inv <- 0.5*lav_matrix_duplication_pre_post(S11 %x% S11)

        if(meanstructure || slopestructure) {
            C <- S[ x.idx,  x.idx, drop=FALSE]
            MY <- M[-x.idx]; MX <- M[x.idx]
            C3 <- rbind(c(1,MX),
                        cbind(MX, C + tcrossprod(MX)))
        }

        if(meanstructure) {
            if(slopestructure) {
                A11 <- C3 %x% S11
            } else {
                c11 <- 1 / solve(C3)[1, 1, drop=FALSE]
                A11 <- c11 %x% S11
            }
        } else {
            if(slopestructure) {
                A11 <- C %x% S11
            } else {
                A11 <- matrix(0,0,0)
            }
        }

        if(meanstructure || slopestructure) {
            Gamma.inv <- lav_matrix_bdiag(A11, Gamma.inv)
        }
    }

    Gamma.inv
}