R/lav_representation_ram.R

In lavaan: Latent Variable Analysis

# RAM representation
#
# initial version: YR 2021-10-04

lav_ram <- function(partable = NULL,
                    target   = NULL,
                    extra    = FALSE,
                    remove.nonexisting = TRUE) {

    # prepare target list
    if(is.null(target)) target <- partable

    stopifnot(!is.null(target$block))

    # not for categorical data (yet)
    if(any(partable$op == "|")) {
        stop("lavaan ERROR: RAM representation is not (yet) supported for categorical endogenous variables.")
    }

    # not for conditional.x = TRUE yet
    conditional.x <- any(partable$exo > 0L & partable$op == "~")
    if(conditional.x) {
        stop("lavaan ERROR: RAM representation is not (yet) supported if conditional.x = TRUE")
    }

    # prepare output
    N <- length(target$lhs)
    tmp.mat <- character(N); tmp.row <- integer(N); tmp.col <- integer(N)

    # global settings
    meanstructure <- any(partable$op == "~1")
    categorical   <- any(partable$op == "|")
    group.w.free  <- any(partable$lhs == "group" & partable$op == "%")

    # number of blocks
    nblocks <- lav_partable_nblocks(partable)

    # always return ov.idx
    ov.idx <- vector("list", nblocks)
    ov.dummy.names.nox <- vector("list", nblocks)
    ov.dummy.names.x   <- vector("list", nblocks)
    if(extra) {
        REP.mmNames     <- vector("list", nblocks)
        REP.mmNumber    <- vector("list", nblocks)
        REP.mmRows      <- vector("list", nblocks)
        REP.mmCols      <- vector("list", nblocks)
        REP.mmDimNames  <- vector("list", nblocks)
        REP.mmSymmetric <- vector("list", nblocks)
    }

    for(g in 1:nblocks) {

        # info from user model per block
        ov.names   <- vnames(partable, "ov",  block = g)
        nvar <- length(ov.names)
        ov.idx[[g]] <- seq_len(nvar)
        ov.dummy.names.nox[[g]] <- character(0)
        ov.dummy.names.x[[g]]   <- character(0)

        lv.names   <- vnames(partable, "lv",  block = g)
        both.names <- c(ov.names, lv.names)
        nboth <- length(both.names)

        # 1. "=~" indicators
        idx <- which(target$block == g & target$op == "=~")
        tmp.mat[idx] <- "A"
        tmp.row[idx] <- match(target$rhs[idx], both.names)
        tmp.col[idx] <- match(target$lhs[idx], both.names)

        # 2. "~" regressions
        idx <- which(target$block == g & (target$op == "~" |
                                          target$op == "<~") )
        tmp.mat[idx] <- "A"
        tmp.row[idx] <- match(target$lhs[idx], both.names)
        tmp.col[idx] <- match(target$rhs[idx], both.names)

        # 3. "~~" variances/covariances
        idx <- which(target$block == g & target$op == "~~")
        tmp.mat[idx] <- "S"
        tmp.row[idx] <- match(target$lhs[idx], both.names)
        tmp.col[idx] <- match(target$rhs[idx], both.names)

        # catch lower-elements in theta/psi
        idx.lower <- which(tmp.mat == "S" & tmp.row > tmp.col)
        if(length(idx.lower) > 0L) {
            tmp <- tmp.row[idx.lower]
            tmp.row[idx.lower] <- tmp.col[idx.lower]
            tmp.col[idx.lower] <- tmp
        }

        # 4. "~1" means/intercepts
        idx <- which(target$block == g & target$op == "~1")
        tmp.mat[idx] <- "m"
        tmp.row[idx] <- match(target$lhs[idx], both.names)
        tmp.col[idx] <- 1L

        # 5. "|" th
        # not used yet

        # 6. "~*~" scales
        # not used yet

        # 7. group weights
        idx <- which(target$block == g & target$lhs == "group" &
                     target$op == "%")
        tmp.mat[idx] <- "gw"
        tmp.row[idx] <- 1L
        tmp.col[idx] <- 1L

        if(extra) {
            # mRows
            mmRows <- list(ov.idx = 1L,
                           A      = nboth,
                           S      = nboth,
                           m      = nboth,
                           gw     = 1L)

            # mCols
            mmCols <- list(ov.idx = nvar,
                           A      = nboth,
                           S      = nboth,
                           m      = 1L,
                           gw     = 1L)

            # dimNames for LISREL model matrices
            mmDimNames <- list(ov.idx = list( "ov.idx",     ov.names),
                               A      = list( both.names, both.names),
                               S      = list( both.names, both.names),
                               m      = list( both.names, "intercept"),
                               gw      = list( "group",       "weight"))
            # isSymmetric
            mmSymmetric <- list(ov.idx = FALSE,
                                A      = FALSE,
                                S      = TRUE,
                                m      = FALSE,
                                gw     = FALSE)

            # which mm's do we need? (always include ov.idx, A and S)
            IDX <- which(target$block == g)
            mmNames <- c("ov.idx", "A", "S")
            if(meanstructure) {
                mmNames <- c(mmNames, "m")
            }
            if("gw" %in% tmp.mat[IDX]) {
                mmNames <- c(mmNames, "gw")
            }

            REP.mmNames[[g]]     <- mmNames
            REP.mmNumber[[g]]    <- length(mmNames)
            REP.mmRows[[g]]      <- unlist(mmRows[ mmNames ])
            REP.mmCols[[g]]      <- unlist(mmCols[ mmNames ])
            REP.mmDimNames[[g]]  <- mmDimNames[ mmNames ]
            REP.mmSymmetric[[g]] <- unlist(mmSymmetric[ mmNames ])
        } # extra
    } # nblocks

    REP <- list(mat = tmp.mat,
                row = tmp.row,
                col = tmp.col)

    # always return ov.idx attribute
    attr(REP, "ov.idx")             <- ov.idx
    attr(REP, "ov.dummy.names.nox") <- ov.dummy.names.nox
    attr(REP, "ov.dummy.names.x")   <- ov.dummy.names.x

    if(extra) {
        attr(REP, "mmNames")     <- REP.mmNames
        attr(REP, "mmNumber")    <- REP.mmNumber
        attr(REP, "mmRows")      <- REP.mmRows
        attr(REP, "mmCols")      <- REP.mmCols
        attr(REP, "mmDimNames")  <- REP.mmDimNames
        attr(REP, "mmSymmetric") <- REP.mmSymmetric
    }

    REP
}

# the model-implied variance/covariance matrix of the observed variables
lav_ram_sigmahat <- function(MLIST = NULL, delta = NULL) {

    ov.idx <- as.integer(MLIST$ov.idx[1,])
    A      <- MLIST$A
    S      <- MLIST$S

    # get (I-A)^{-1}
    IA.inv <- lav_matrix_inverse_iminus(A)

    # compute Sigma for all ov and lv
    VYeta <- tcrossprod(IA.inv %*% S, IA.inv)

    # select only observed part
    VY <- VYeta[ov.idx, ov.idx, drop = FALSE]

    # if delta, scale
    if(!is.null(MLIST$delta) && delta) {
        nvar <- ncol(VY)
        DELTA <- diag(MLIST$delta[,1L], nrow = nvar, ncol = nvar)
        VY <- DELTA %*% VY %*% DELTA
    }

    VY
}

# VETA: the variance/covariance matrix of the latent variables only
lav_ram_veta <- function(MLIST = NULL) {

    ov.idx <- as.integer(MLIST$ov.idx[1,])
    A      <- MLIST$A
    S      <- MLIST$S

    # get (I-A)^{-1}
    IA.inv <- lav_matrix_inverse_iminus(A)

    # compute Sigma for all ov and lv
    VYeta <- tcrossprod(IA.inv %*% S, IA.inv)

    # select only latent part
    VETA <- VYeta[-ov.idx, -ov.idx, drop = FALSE]

    VETA
}

# MuHat: the model-implied means/intercepts
lav_ram_muhat <- function(MLIST = NULL) {

    ov.idx <- as.integer(MLIST$ov.idx[1,])
    A      <- MLIST$A
    m      <- MLIST$m

    # shortcut
    if(is.null(m)) return(matrix(0, nrow = length(ov.idx), 1L))

    # get (I-A)^{-1}
    IA.inv <- lav_matrix_inverse_iminus(A)

    # all means/intercepts
    EYeta <- IA.inv %*% m

    # select observed only
    muhat <- EYeta[ov.idx, , drop = FALSE]

    muhat
}

# derivative of 'Sigma' wrt the (freel) elements in A and/or S
lav_ram_dsigma <- function(m = "A",
                           idx = seq_len(length(MLIST[[m]])),
                           MLIST = NULL,
                           vech  = TRUE) {

    ov.idx <- as.integer(MLIST$ov.idx[1,])
    A      <- MLIST$A
    S      <- MLIST$S

    nvar  <- length(ov.idx)
    nboth <- nrow(A)

    # shortcut for ov.idx, m, ...
    if(!m %in% c("A", "S")) {
        pstar <- nvar*(nvar+1)/2
        return( matrix(0.0, nrow = pstar, ncol = length(idx)) )
    }

    # get (I-A)^{-1}
    IA.inv <- lav_matrix_inverse_iminus(A)

    if(m == "A") {
        L1 <- (IA.inv %*% S %*% t(IA.inv))[ov.idx,,drop = FALSE]
        KOL.idx <- matrix(1:(nboth*nboth), nboth, nboth, byrow = TRUE)[idx]
        DX <- (L1 %x% IA.inv[ov.idx,,drop = FALSE])[,idx, drop = FALSE] +
              (IA.inv[ov.idx,,drop = FALSE] %x% L1)[, KOL.idx, drop = FALSE]
        # this is not really needed (because we select idx=m.el.idx)
        # but just in case we need all elements of beta...
        DX[, which(idx %in% lav_matrix_diag_idx(nboth))] <- 0.0
    } else if(m == "S") {
        DX <- (IA.inv[ov.idx,,drop = FALSE] %x% IA.inv[ov.idx,,drop = FALSE])
        # symmetry correction, but keeping all duplicated elements
        # since we depend on idx=m.el.idx
        lower.idx <- lav_matrix_vech_idx(nboth, diagonal = FALSE)
        upper.idx <- lav_matrix_vechru_idx(nboth, diagonal = FALSE)
        offdiagSum <- DX[,lower.idx] + DX[,upper.idx]
        DX[,c(lower.idx, upper.idx)] <- cbind(offdiagSum, offdiagSum)
        DX <- DX[,idx, drop = FALSE]
    } else {
        stop("wrong model matrix names: ", m, "\n")
    }

    # vech?
    if(vech) {
        v.idx <- lav_matrix_vech_idx(nvar)
        DX <- DX[v.idx,, drop = FALSE]
    }

    DX
}

# derivative of 'Mu' wrt the (free) elements in A and/or m
lav_ram_dmu <- function(m = "A",
                        idx = seq_len(length(MLIST[[m]])),
                        MLIST = NULL,
                        vech  = TRUE) {

    ov.idx <- as.integer(MLIST$ov.idx[1,])
    A      <- MLIST$A
    S      <- MLIST$S

    nvar  <- length(ov.idx)
    nboth <- nrow(A)

    # shortcut for ov.idx, m, ...
    if(!m %in% c("A", "m")) {
        return( matrix(0.0, nrow = nvar, ncol = length(idx)) )
    }

    # get (I-A)^{-1}
    IA.inv <- lav_matrix_inverse_iminus(A)

    if(m == "A") {
        DX <- (t(IA.inv %*% MLIST$m) %x% IA.inv)[ov.idx, idx, drop = FALSE]
    } else if(m == "m") {
        DX <- IA.inv[ov.idx, idx, drop = FALSE]
    } else {
        stop("wrong model matrix names: ", m, "\n")
    }

    DX
}

# derivative of ML/GLS objective function F wrt the free parameters
lav_ram_df <- function(MLIST = NULL, Omega = NULL, Omega.mu = NULL) {

    ov.idx <- as.integer(MLIST$ov.idx[1,])
    A      <- MLIST$A
    S      <- MLIST$S

    nvar  <- length(ov.idx)
    nboth <- nrow(A)

    # get (I-A)^{-1}
    IA.inv <- lav_matrix_inverse_iminus(A)

    # meanstructure?
    meanstructure <- FALSE; if(!is.null(Omega.mu)) meanstructure <- TRUE

    # pre-compute
    tIA.inv <- t(IA.inv)
    Omega..IA.inv..S..tIA.inv <- ( Omega %*% IA.inv[ov.idx,,drop = FALSE] %*% S %*% t(IA.inv) )

    # 1. A
    if(meanstructure) {
        A.deriv <-
            -1.0*(( t(IA.inv)[,ov.idx,drop = FALSE] %*% (Omega.mu %*% t(MLIST$m)) %*% t(IA.inv)) +
                  (tIA.inv[,ov.idx,drop = FALSE] %*% Omega..IA.inv..S..tIA.inv))
    } else {
        A.deriv <- -1.0 * ( tIA.inv[,ov.idx,drop = FALSE] %*% Omega..IA.inv..S..tIA.inv )
    }

    # 2. S
    S.deriv <- -1.0 * ( tIA.inv[,ov.idx,drop = FALSE] %*% Omega %*% IA.inv[ov.idx,,drop = FALSE] )
    diag(S.deriv) <- 0.5 * diag(S.deriv)

    if(meanstructure) {
        m.deriv <- -1.0 * t( t(Omega.mu) %*% IA.inv[ov.idx,,drop = FALSE] )
    } else {
        m.deriv <- NULL
    }

    list(A = A.deriv,
         S = S.deriv,
         m = m.deriv)
}