R/lav_samplestats_gamma.R

In lavaan: Latent Variable Analysis

# YR 21 March 2015
# new approach to compute 'Gamma': the asymptotic variance matrix of
#                                  sqrt{N} times the
#                                  observed sample statistics (means + varcov)
#
# Gamma = N x ACOV[ ybar, vech(S) ]
#       = NACOV[ ybar, vech(S) ]
#
# - one single function for mean + cov
# - handle 'fixed.x' exogenous covariates


# - YR  3 Dec 2015: allow for conditional.x = TRUE
# - YR 22 Jan 2023: add model.based= argument (if object is lavaan object)

# generic public function (not exported yet)
# input for lavGamma can be lavobject, lavdata, data.frame, or matrix
lavGamma <- function(object, group = NULL, missing = "listwise",
                     ov.names.x = NULL, fixed.x = FALSE, conditional.x = FALSE,
                     meanstructure = FALSE, slopestructure = FALSE,
                     gamma.n.minus.one = FALSE, gamma.unbiased = FALSE,
                     ADF = TRUE, model.based = FALSE, NT.rescale = FALSE,
                     Mplus.WLS = FALSE, add.labels = FALSE) {

    # check object

    # 1. object is lavaan object
    if(inherits(object, "lavaan")) {
        lav_object_gamma(lavobject = object, ADF = ADF,
                         model.based = model.based, Mplus.WLS = Mplus.WLS)
    }
    else if(inherits(object, "lavData")) {
        lavdata <- object
        model.based <- FALSE
    } else if(inherits(object, "data.frame") ||
              inherits(object, "matrix")) {
        model.based <- FALSE
        NAMES <- names(object)
        if(!is.null(NAMES) && !is.null(group)) {
            NAMES <- NAMES[- match(group, NAMES)]
        }
        lavdata <- lavData(data = object, group = group,
                           ov.names = NAMES, ordered = NULL,
                           ov.names.x = ov.names.x,
                           lavoptions = list(warn = FALSE,
                                             missing = missing))
    } else {
        stop("lavaan ERROR: lavGamma can not handle objects of class ",
             paste(class(object), collapse= " "))
    }

    # extract data
    Y <- lavdata@X
    if(conditional.x) {
        eXo <- lavdata@eXo
        for(g in seq_len(lavdata@ngroups)) {
            Y[[g]] <- cbind(Y[[g]], eXo[[g]])
        }
    }

    # x.idx
    x.idx <- lapply(seq_len(lavdata@ngroups),
                    function(g) match(lavdata@ov.names.x[[g]],
                                      lavdata@ov.names[[g]]) )

    OUT <- lapply(seq_len(lavdata@ngroups), function(g) {
        if(length(lavdata@cluster) > 0L) {
            cluster.idx <- lavdata@Lp[[g]]$cluster.idx[[2]]
        } else {
            cluster.idx <- NULL
        }
        if(ADF) {
            out <- lav_samplestats_Gamma(Y          = Y[[g]],
                                     Mu             = NULL,
                                     Sigma          = NULL,
                                     x.idx          = x.idx[[g]],
                                     cluster.idx    = cluster.idx,
                                     fixed.x        = fixed.x,
                                     conditional.x  = conditional.x,
                                     meanstructure  = meanstructure,
                                     slopestructure = conditional.x,
                                     gamma.n.minus.one = gamma.n.minus.one,
                                     unbiased       = gamma.unbiased,
                                     Mplus.WLS      = Mplus.WLS)
        } else {
            out <- lav_samplestats_Gamma_NT(Y   = Y[[g]],
                                 wt             = NULL, # for now
                                 rescale        = NT.rescale,
                                 x.idx          = x.idx[[g]],
                                 fixed.x        = fixed.x,
                                 conditional.x  = conditional.x,
                                 meanstructure  = meanstructure,
                                 slopestructure = conditional.x)
        }
        out})

    # todo: labels

    OUT
}



# for internal use -- lavobject or internal slots
lav_object_gamma <- function(lavobject = NULL,
                             # or individual slots
                             lavdata        = NULL,
                             lavoptions     = NULL,
                             lavsamplestats = NULL,
                             lavh1          = NULL,
                             lavimplied     = NULL,
                             # other options
                             ADF = TRUE, model.based = FALSE,
                             Mplus.WLS = FALSE) {

    # extract slots
    if(!is.null(lavobject)) {
        lavdata        <- lavobject@Data
        lavoptions     <- lavobject@Options
        lavsamplestats <- lavobject@SampleStats
        if(.hasSlot(lavobject, "h1")) {
            lavh1          <- lavobject@h1
        } else {
            # only for <0.6
            out <- lav_h1_implied_logl(lavdata        = lavdata,
                                       lavsamplestats = lavsamplestats,
                                       lavpta         = lavobject@pta,
                                       lavoptions     = lavoptions)
            h1.implied      <- out$implied
            h1.loglik       <- out$logl$loglik
            h1.loglik.group <- out$logl$loglik.group
            lavh1 <- list(implied      = h1.implied,
                          loglik       = h1.loglik,
                          loglik.group = h1.loglik.group)
            lavoptions$gamma.n.minus.one <- FALSE
            lavoptions$gamma.unbiased <- FALSE
        }
        lavimplied     <- lavobject@implied


    }

    missing <- lavoptions$missing
    if(!missing %in% c("listwise", "pairwise")) {
            model.based <- TRUE
    }
    fixed.x           <- lavoptions$fixed.x
    conditional.x     <- lavoptions$conditional.x
    meanstructure     <- lavoptions$meanstructure
    gamma.n.minus.one <- lavoptions$gamma.n.minus.one
    gamma.unbiased    <- lavoptions$gamma.unbiased

    if(ADF && model.based && conditional.x) {
        stop("lavaan ERROR: ADF + model.based + conditional.x is not supported yet.")
    }

    # output container
    OUT <- vector("list", length = lavdata@ngroups)

    # compute Gamma matrix for each group
    for(g in seq_len(lavdata@ngroups)) {
        x.idx <- lavsamplestats@x.idx[[g]]
        COV <- MEAN <- NULL
        if(!ADF || model.based) {
            implied <- lavh1$implied   # saturated/unstructured
            if(model.based) {
                implied <- lavimplied   # model-based/structured
            }
            if(conditional.x) {
                # convert to joint COV/MEAN
                res.S      <- implied$res.cov[[g]]
                res.slopes <- implied$res.slopes[[g]]
                res.int    <- implied$res.int[[g]]
                S.xx       <- implied$cov.x[[g]]
                M.x        <- implied$mean.x[[g]]

                S.yy <- res.S + res.slopes %*% S.xx %*% t(res.slopes)
                S.yx <- res.slopes %*% S.xx
                S.xy <- S.xx %*% t(res.slopes)
                M.y <- res.int + res.slopes %*% M.x

                COV <- rbind( cbind(S.yy, S.yx), cbind(S.xy, S.xx) )
                MEAN <- c(M.y, M.x)
            } else {
                # not conditional.x
                COV  <- implied$cov[[g]]
                MEAN <- implied$mean[[g]]
            }
        } # COV/MEAN

        if(ADF) {
            if(conditional.x) {
                Y <- cbind(lavdata@X[[g]], lavdata@eXo[[g]])
            } else {
                Y <- lavdata@X[[g]]
            }
            if(length(lavdata@cluster) > 0L) {
                cluster.idx <- lavdata@Lp[[g]]$cluster.idx[[2]]
            } else {
                cluster.idx <- NULL
            }
            OUT[[g]] <- lav_samplestats_Gamma(Y = Y,
                                 Mu             = MEAN,
                                 Sigma          = COV,
                                 x.idx          = x.idx,
                                 cluster.idx    = cluster.idx,
                                 fixed.x        = fixed.x,
                                 conditional.x  = conditional.x,
                                 meanstructure  = meanstructure,
                                 slopestructure = conditional.x,
                                 gamma.n.minus.one = gamma.n.minus.one,
                                 unbiased       = gamma.unbiased,
                                 Mplus.WLS      = Mplus.WLS)
        } else {
            OUT[[g]] <- lav_samplestats_Gamma_NT(COV = COV,  # joint!
                                      MEAN           = MEAN, # joint!
                                      x.idx          = x.idx,
                                      fixed.x        = fixed.x,
                                      conditional.x  = conditional.x,
                                      meanstructure  = meanstructure,
                                      slopestructure = conditional.x)
        }
    } # g

    OUT
}





# 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_NT <- function(Y              = NULL, # should include
                                                            # eXo if
                                                            #conditional.x=TRUE
                                     wt             = NULL,
                                     COV            = NULL, # joint!
                                     MEAN           = NULL, # joint!
                                     rescale        = FALSE,
                                     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
    }

    # compute COV from Y
    if(is.null(COV)) {
        stopifnot(!is.null(Y))

        # coerce to matrix
        Y <- unname(as.matrix(Y)); N <- nrow(Y)
        if(is.null(wt)) {
            COV <- cov(Y)
            if(rescale) {
                COV <- COV * (N-1) / N # (normal) ML version
            }
        } else {
            out <- stats::cov.wt(Y, wt = wt, method = "ML")
            COV <- out$cov
        }
    } else {
        if(!missing(rescale)) {
            warning("lavaan WARNING: rescale= argument has no effect if COV is given")
        }
        if(!missing(wt)) {
            warning("lavaan WARNING: wt= argument has no effect if COV is given")
        }
    }

    # if needed, compute MEAN from Y
    if(conditional.x && length(x.idx) > 0L && is.null(MEAN) &&
       (meanstructure || slopestructure)) {
       stopifnot(!is.null(Y))
       if(is.null(wt)) {
           MEAN <- colMeans(Y, na.rm = TRUE)
       } else {
           MEAN <- out$center
       }
    }

    # rename
    S <- COV
    M <- MEAN

    # unconditional
    if(!conditional.x) {

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

        # unconditional - fixed x
        } else {
            # handle fixed.x = TRUE

            # cov(Y|X) = A - B C^{-1} B'
            # where A = cov(Y), B = cov(Y,X), C = cov(X)
            A <- S[-x.idx, -x.idx, drop=FALSE]
            B <- S[-x.idx,  x.idx, drop=FALSE]
            C <- S[ x.idx,  x.idx, drop=FALSE]
            YbarX <- A - B %*% solve(C, t(B))

            # reinsert YbarX in Y+X (residual) covariance matrix
            YbarX.aug <- matrix(0, nrow = NROW(S), ncol = NCOL(S))
            YbarX.aug[ -x.idx, -x.idx ] <- YbarX

            # take difference
            R <- S - YbarX.aug

            Gamma.S <- 2*lav_matrix_duplication_ginv_pre_post(S %x% S)
            Gamma.R <- 2*lav_matrix_duplication_ginv_pre_post(R %x% R)
            Gamma <- Gamma.S - Gamma.R

            if(meanstructure) {
                Gamma <- lav_matrix_bdiag(YbarX.aug, Gamma)
            }
        }

    } else {
        # conditional.x

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

        # regress Y on X, and compute covariance of residuals 'R'
        A <- S[-x.idx, -x.idx, drop=FALSE]
        B <- S[-x.idx,  x.idx, drop=FALSE]
        C <- S[ x.idx,  x.idx, drop=FALSE]
        Cov.YbarX <- A - B %*% solve(C) %*% t(B)
        Gamma <- 2*lav_matrix_duplication_ginv_pre_post(Cov.YbarX %x% Cov.YbarX)

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

        if(meanstructure) {
            if(slopestructure) {
                #A11 <- solve(C3) %x% Cov.YbarX
                A11 <- Cov.YbarX %x% solve(C3)
            } else {
                #A11 <- solve(C3)[1, 1, drop=FALSE] %x% Cov.YbarX
                A11 <- Cov.YbarX %x% solve(C3)[1, 1, drop=FALSE]
            }
        } else {
            if(slopestructure) {
                #A11 <- solve(C3)[-1, -1, drop=FALSE] %x% Cov.YbarX
                A11 <- Cov.YbarX %x% solve(C3)[-1, -1, drop=FALSE]
            } else {
                A11 <- matrix(0,0,0)
            }
        }

        Gamma <- lav_matrix_bdiag(A11, Gamma)
    }

    Gamma
}

# 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)
#
#  - new in 0.6-1: if Mu/Sigma is provided, compute 'model-based' Gamma
#                  (only if conditional.x = FALSE, for now)
#  - new in 0.6-2: if cluster.idx is not NULL, correct for clustering
#  - new in 0.6-13: add unbiased = TRUE (for the 'plain' setting only)

# ADF THEORY
lav_samplestats_Gamma <- function(Y,                          # Y+X if cond!
                                  Mu                 = NULL,
                                  Sigma              = NULL,
                                  x.idx              = integer(0L),
                                  cluster.idx        = NULL,
                                  fixed.x            = FALSE,
                                  conditional.x      = FALSE,
                                  meanstructure      = FALSE,
                                  slopestructure     = FALSE,
                                  gamma.n.minus.one  = FALSE,
                                  unbiased           = FALSE,
                                  Mplus.WLS          = FALSE) {
    # coerce to matrix
    Y <- unname(as.matrix(Y)); N <- nrow(Y); p <- ncol(Y)

    # unbiased?
    if(unbiased) {
        if(conditional.x || fixed.x || !is.null(Sigma) ||
           !is.null(cluster.idx)) {
            stop("lavaan ERROR: unbiased Gamma only available for the simple (not conditional.x or fixed.x or model-based or clustered) setting.")
        } else {
            COV <- COV.unbiased <- cov(Y)
            COV <- COV * (N-1)/N
            cov.vech <- lav_matrix_vech(COV)
        }
    }

    # model-based?
    if(!is.null(Sigma)) {
        stopifnot(!conditional.x)
        model.based <- TRUE
        if(meanstructure) {
            stopifnot(!is.null(Mu))
            sigma <- c(as.numeric(Mu), lav_matrix_vech(Sigma))
        } else {
            Mu <- colMeans(Y, na.rm = TRUE) # for centering!
            sigma <- lav_matrix_vech(Sigma)
        }
    } else {
        model.based <- FALSE
    }

    # denominator
    if(gamma.n.minus.one) {
        N <- N - 1
    }

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

    if(!conditional.x && !fixed.x) {
        # center, so we can use crossprod instead of cov
        if(model.based) {
            Yc <- t( t(Y) - as.numeric(Mu) )
        } else {
            Yc <- t( t(Y) - colMeans(Y, na.rm = TRUE) )
        }

        # create Z where the rows_i contain the following elements:
        #  - Y_i (if meanstructure is TRUE)
        #  - vech(Yc_i' %*% Yc_i) where Yc_i are the residuals
        idx1 <- lav_matrix_vech_col_idx(p)
        idx2 <- lav_matrix_vech_row_idx(p)
        if(meanstructure) {
            Z <- cbind(Y, Yc[,idx1, drop = FALSE] *
                          Yc[,idx2, drop = FALSE] )
        } else {
            Z <- ( Yc[,idx1, drop = FALSE] *
                   Yc[,idx2, drop = FALSE] )
        }

        if(model.based) {
            if(meanstructure) {
                stopifnot(!is.null(Mu))
                sigma <- c(as.numeric(Mu), lav_matrix_vech(Sigma))
            } else {
                sigma <- lav_matrix_vech(Sigma)
            }
            Zc <- t( t(Z) - sigma )
        } else {
            Zc <- t( t(Z) - colMeans(Z, na.rm = TRUE) )
        }

        # clustered?
        if(length(cluster.idx) > 0L) {
            Zc <- rowsum(Zc, cluster.idx)
        }

        if(anyNA(Zc)) {
            Gamma <- lav_matrix_crossprod(Zc) / N
        } else {
            Gamma <- base::crossprod(Zc) / N
        }
    } else if(!conditional.x && fixed.x) {

        if(model.based) {
            Yc <- t( t(Y) - as.numeric(Mu) )
            Y.bar <- colMeans(Y, na.rm = TRUE)
            res.cov <- ( Sigma[-x.idx, -x.idx, drop = FALSE] -
                         Sigma[-x.idx, x.idx, drop = FALSE] %*%
                          solve(Sigma[x.idx, x.idx, drop = FALSE]) %*%
                         Sigma[x.idx, -x.idx, drop = FALSE] )
            res.slopes <- ( solve(Sigma[x.idx, x.idx, drop = FALSE]) %*%
                            Sigma[x.idx, -x.idx, drop = FALSE] )
            res.int <- ( Y.bar[-x.idx] -
                         as.numeric(colMeans(Y[,x.idx,drop = FALSE],
                                    na.rm = TRUE) %*% res.slopes) )
            x.bar <- Y.bar[x.idx]
            yhat.bar <- as.numeric(res.int + as.numeric(x.bar) %*% res.slopes)
            YHAT.bar <- numeric(p)
            YHAT.bar[-x.idx] <- yhat.bar; YHAT.bar[x.idx] <- x.bar
            YHAT.cov <- Sigma
            YHAT.cov[-x.idx, -x.idx] <- Sigma[-x.idx, -x.idx] - res.cov


            yhat <- cbind(1, Y[,x.idx]) %*% rbind(res.int, res.slopes)
            YHAT <- Y; YHAT[,-x.idx] <- yhat
            #YHAT <- cbind(yhat, Y[,x.idx])
            YHATc <- t( t(YHAT) - YHAT.bar )
            idx1 <- lav_matrix_vech_col_idx(p)
            idx2 <- lav_matrix_vech_row_idx(p)
            if(meanstructure) {
                Z <- ( cbind(Y, Yc[,idx1, drop = FALSE] *
                                Yc[,idx2, drop = FALSE] ) -
                       cbind(YHAT, YHATc[,idx1, drop = FALSE] *
                                   YHATc[,idx2, drop = FALSE]) )
                sigma1 <- c(Mu, lav_matrix_vech(Sigma))
                sigma2 <- c(YHAT.bar, lav_matrix_vech(YHAT.cov))
            } else {
                Z <- ( Yc[,idx1, drop = FALSE] *
                       Yc[,idx2, drop = FALSE]  -
                       YHATc[,idx1, drop = FALSE] *
                       YHATc[,idx2, drop = FALSE] )
                sigma1 <- lav_matrix_vech(Sigma)
                sigma2 <- lav_matrix_vech(YHAT.cov)
            }
            Zc <- t( t(Z) - (sigma1 - sigma2) )
        } else {
            QR <- qr(cbind(1, Y[, x.idx, drop = FALSE]))
            yhat <- qr.fitted(QR, Y[, -x.idx, drop = FALSE])
            # YHAT <- cbind(yhat, Y[,x.idx])
            YHAT <- Y; YHAT[,-x.idx] <- yhat

            Yc <- t( t(Y) - colMeans(Y, na.rm = TRUE) )
            YHATc <- t( t(YHAT) - colMeans(YHAT, na.rm = TRUE) )
            idx1 <- lav_matrix_vech_col_idx(p)
            idx2 <- lav_matrix_vech_row_idx(p)
            if(meanstructure) {
                Z <- ( cbind(Y, Yc[,idx1, drop = FALSE] *
                                Yc[,idx2, drop = FALSE]) -
                       cbind(YHAT, YHATc[,idx1, drop = FALSE] *
                                   YHATc[,idx2, drop = FALSE]) )
            } else {
                Z <- ( Yc[,idx1, drop = FALSE] *
                       Yc[,idx2, drop = FALSE] -
                       YHATc[,idx1, drop = FALSE] *
                       YHATc[,idx2, drop = FALSE] )
            }
            Zc <- t( t(Z) - colMeans(Z, na.rm = TRUE) )
        }

        # clustered?
        if(length(cluster.idx) > 0L) {
            Zc <- rowsum(Zc, cluster.idx)
        }

        if(anyNA(Zc)) {
            Gamma <- lav_matrix_crossprod(Zc) / N
        } else {
            Gamma <- base::crossprod(Zc) / N
        }


    } else {
        # conditional.x

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

        # regress Y on X, and compute residuals
        X    <- cbind(1, Y[,  x.idx, drop = FALSE])
        QR   <- qr(X)
        RES  <- qr.resid(QR, Y[, -x.idx, drop = FALSE])
        p    <- ncol(RES)

        idx1 <- lav_matrix_vech_col_idx(p)
        idx2 <- lav_matrix_vech_row_idx(p)

        if(meanstructure || slopestructure) {
            XtX.inv <- unname(solve(crossprod(X)))
            Xi <- (X %*% XtX.inv) * N ## FIXME, shorter way?
            ncX <- NCOL(X); ncY <- NCOL(RES)
        }

        if(meanstructure) {
            if(slopestructure) {
                # Xi.idx <- rep(seq_len(ncX), each  = ncY)
                #Res.idx <- rep(seq_len(ncY), times = ncX)
                 Xi.idx <- rep(seq_len(ncX), times = ncY)
                Res.idx <- rep(seq_len(ncY), each  = ncX)
                Z <- cbind( Xi[, Xi.idx, drop = FALSE] *
                           RES[,Res.idx, drop = FALSE],
                           RES[,   idx1, drop = FALSE] *
                           RES[,   idx2, drop = FALSE] )
            } else {
                Xi.idx <- rep(1L, each = ncY)
                Z <- cbind( Xi[, Xi.idx ,drop = FALSE] *
                           RES,
                           RES[,   idx1, drop = FALSE] *
                           RES[,   idx2, drop = FALSE] )
            }
        } else {
            if(slopestructure) {
                # Xi.idx <- rep(seq_len(ncX), each  = ncY)
                # Xi.idx <- Xi.idx[ -seq_len(ncY) ]
                Xi.idx <- rep(seq(2,ncX), times = ncY)
                #Res.idx <- rep(seq_len(ncY), times = (ncX - 1L))
                Res.idx <- rep(seq_len(ncY), each = (ncX - 1L))
                Z <- cbind( Xi[, Xi.idx, drop = FALSE] *
                           RES[,Res.idx, drop = FALSE],
                           RES[,   idx1, drop = FALSE] *
                           RES[,   idx2, drop = FALSE] )
            } else {
                Z <- RES[,idx1, drop = FALSE] * RES[,idx2, drop = FALSE]
            }
        }

        if(model.based) {
            Zc <- t( t(Z) - sigma )
        } else {
            Zc <- t( t(Z) - colMeans(Z, na.rm = TRUE) )
        }

        # clustered?
        if(length(cluster.idx) > 0L) {
            Zc <- rowsum(Zc, cluster.idx)
        }

        if(anyNA(Zc)) {
            Gamma <- lav_matrix_crossprod(Zc) / N
        } else {
            Gamma <- base::crossprod(Zc) / N
        }
    }


    # only to mimic Mplus when estimator = "WLS"
    if(Mplus.WLS && !fixed.x && !conditional.x) {
        # adjust G_22 (the varcov part)
        S <- cov(Y, use = "pairwise")
        w <- lav_matrix_vech(S)
        w.biased <- (N-1)/N * w
        diff <- outer(w,w) - outer(w.biased, w.biased)
        if(meanstructure) {
            Gamma[-seq_len(p), -seq_len(p)] <-
                Gamma[-seq_len(p), -seq_len(p), drop = FALSE] - diff
        } else {
            Gamma <- Gamma - diff
        }

        if(meanstructure) {
            # adjust G_12/G_21 (third-order)
            # strange rescaling?
            N1 <- (N - 1) / N
            Gamma[seq_len(p),-seq_len(p)] <- Gamma[seq_len(p),-seq_len(p)] * N1
            Gamma[-seq_len(p),seq_len(p)] <- Gamma[-seq_len(p),seq_len(p)] * N1
        }
    }

    # clustered?
    if(length(cluster.idx) > 0L) {
        nC <- nrow(Zc)
        Gamma <- Gamma * nC / (nC - 1)
    }

    # unbiased?
    if(unbiased) {
        # normal-theory Gamma (cov only)
        GammaNT.cov <- 2 * lav_matrix_duplication_ginv_pre_post(COV %x% COV)

        if(meanstructure) {
            Gamma.cov      <- Gamma[-(1:p), -(1:p), drop = FALSE]
            Gamma.mean.cov <- Gamma[1:p,    -(1:p), drop = FALSE]
        } else {
            Gamma.cov <- Gamma
        }

        # Browne's unbiased DF estimator (COV part)
        Gamma.u <- ( N*(N-1)/(N-2)/(N-3) * Gamma.cov -
                     N/(N-2)/(N-3) * ( GammaNT.cov -
                                       2/(N-1) * tcrossprod(cov.vech) )  )
        if(meanstructure) {
            Gamma <- lav_matrix_bdiag(COV,  Gamma.u)

            # 3-rd order:
            Gamma[1:p,(p+1):ncol(Gamma)] <- Gamma.mean.cov * N/(N-2)
            Gamma[(p+1):ncol(Gamma),1:p] <- t( Gamma.mean.cov * N/(N-2) )

        } else {
            Gamma <- Gamma.u
        }
    } # unbiased

    Gamma
}