R/lav_mvreg_cluster.R

In lavaan: Latent Variable Analysis

# loglikelihood clustered/twolevel data -- conditional.x = TRUE

# YR: first version around Sept 2021

# take model-implied mean+variance matrices, and reorder/augment them
# to facilitate computing of (log)likelihood in the two-level case

# when conditional.x = TRUE:
# - sigma.w and sigma.b: same dimensions, level-1 'Y' variables only
# - sigma.zz: level-2 variables only
# - sigma.yz: cov(level-1, level-2)
# - beta.w: beta y within part
# - beta.b: beta y between part
# - beta.z: beta z (between-only)
lav_mvreg_cluster_implied22l <- function(Lp           = NULL,
                                         implied      = NULL,
                                         Res.Int.W    = NULL,
                                         Res.Int.B    = NULL,
                                         Res.Pi.W     = NULL,
                                         Res.Pi.B     = NULL,
                                         Res.Sigma.W  = NULL,
                                         Res.Sigma.B  = NULL) {
    if(!is.null(implied)) {
        # FIXME: only for single-group analysis!
        Res.Sigma.W <- implied$res.cov[[1]]
        Res.Int.W   <- implied$res.int[[1]]
        Res.Pi.W    <- implied$res.slopes[[1]]

        Res.Sigma.B <- implied$res.cov[[2]]
        Res.Int.B   <- implied$res.int[[2]]
        Res.Pi.B    <- implied$res.slopes[[2]]
    }

    # within/between idx
    within.x.idx  <- Lp$within.x.idx[[1]]
    between.y.idx <- Lp$between.y.idx[[2]]
    between.x.idx <- Lp$between.x.idx[[2]]

    # ov.idx per level
    ov.idx      <- Lp$ov.idx

    # 'tilde' matrices: ALL variables within and between
    p.tilde <- length( unique(c(ov.idx[[1]], ov.idx[[2]])) )

    # only 'y'
    ov.y.idx <- Lp$ov.y.idx

    # two levels only (for now)
    ov.y.idx1 <- ov.y.idx[[1]]
    ov.y.idx2 <- ov.y.idx[[2]]

    # Sigma.W.tilde
    Sigma.W.tilde <- matrix(0, p.tilde, p.tilde)
    Sigma.W.tilde[ ov.y.idx1, ov.y.idx1 ] <- Res.Sigma.W

    # INT.W.tilde
    INT.W.tilde <- matrix(0, p.tilde, 1L)
    INT.W.tilde[ ov.y.idx1, 1L ] <- Res.Int.W

    # PI.W.tilde
    PI.W.tilde <- matrix(0, p.tilde, ncol(Res.Pi.W))
    PI.W.tilde[ ov.y.idx1, ] <- Res.Pi.W

    BETA.W.tilde <- rbind(t(INT.W.tilde), t(PI.W.tilde))



    # Sigma.B.tilde
    Sigma.B.tilde <- matrix(0, p.tilde, p.tilde)
    Sigma.B.tilde[ ov.y.idx2, ov.y.idx2 ] <- Res.Sigma.B

    # INT.B.tilde
    INT.B.tilde <- matrix(0, p.tilde, 1L)
    INT.B.tilde[ ov.y.idx2, 1L ] <- Res.Int.B

    # PI.B.tilde
    PI.B.tilde <- matrix(0, p.tilde, ncol(Res.Pi.B))
    PI.B.tilde[ ov.y.idx2, ] <- Res.Pi.B

    BETA.B.tilde <- rbind(t(INT.B.tilde), t(PI.B.tilde))

    if(length(between.y.idx) > 0L) {
        rm.idx <- c(within.x.idx, between.x.idx, between.y.idx) # between AND x
        beta.z   <- BETA.B.tilde[,  between.y.idx, drop = FALSE]
        beta.b   <- BETA.B.tilde[, -rm.idx, drop = FALSE]
        beta.w   <- BETA.W.tilde[, -rm.idx, drop = FALSE]
        sigma.zz <- Sigma.B.tilde[  between.y.idx, between.y.idx, drop = FALSE]
        sigma.yz <- Sigma.B.tilde[ -rm.idx, between.y.idx, drop = FALSE]
        sigma.b  <- Sigma.B.tilde[ -rm.idx, -rm.idx, drop = FALSE]
        sigma.w  <- Sigma.W.tilde[ -rm.idx, -rm.idx, drop = FALSE]
    } else {
        rm.idx <- c(within.x.idx, between.x.idx) # all 'x'
        beta.z   <- matrix(0, 0L, 0L)
        sigma.zz <- matrix(0, 0L, 0L)
        beta.b   <- BETA.B.tilde[, -rm.idx, drop = FALSE]
        beta.w   <- BETA.W.tilde[, -rm.idx, drop = FALSE]
        sigma.b  <- Sigma.B.tilde[ -rm.idx, -rm.idx, drop = FALSE]
        sigma.w  <- Sigma.W.tilde[ -rm.idx, -rm.idx, drop = FALSE]
        sigma.yz <- matrix(0, nrow(sigma.w), 0L)
    }


    # beta.wb # FIXme: not correct if some 'x' are splitted (overlap)
    # but because we ALWAYS treat splitted-x as 'y', this is not a problem
    beta.wb <- rbind(beta.w, beta.b[-1,,drop = FALSE])
    beta.wb[1,] <- beta.wb[1,,drop = FALSE] + beta.b[1,,drop = FALSE]

    list(sigma.w = sigma.w, sigma.b = sigma.b, sigma.zz = sigma.zz,
         sigma.yz = sigma.yz, beta.w = beta.w, beta.b = beta.b, beta.z = beta.z,
         beta.wb = beta.wb)
}


# recreate implied matrices from 2L matrices
lav_mvreg_cluster_2l2implied <- function(Lp,
                                         sigma.w  = NULL,
                                         sigma.b  = NULL,
                                         sigma.zz = NULL,
                                         sigma.yz = NULL,
                                         beta.w   = NULL,
                                         beta.b   = NULL,
                                         beta.z   = NULL) {

    # within/between idx
    within.x.idx  <- Lp$within.x.idx[[1]]
    between.y.idx <- Lp$between.y.idx[[2]]
    between.x.idx <- Lp$between.x.idx[[2]]

    # ov.idx per level
    ov.idx      <- Lp$ov.idx

    # 'tilde' matrices: ALL variables within and between
    p.tilde <- length( unique(c(ov.idx[[1]], ov.idx[[2]])) )

    # only 'y'
    ov.y.idx <- Lp$ov.y.idx

    # two levels only (for now)
    ov.y.idx1 <- ov.y.idx[[1]]
    ov.y.idx2 <- ov.y.idx[[2]]

    # Sigma.W.tilde
    Sigma.W.tilde <- matrix(0, p.tilde, p.tilde)
    Sigma.W.tilde[ ov.y.idx1, ov.y.idx1 ] <- sigma.w

    # INT.W.tilde
    INT.W.tilde <- matrix(0, p.tilde, 1L)
    INT.W.tilde[ ov.y.idx1, 1L ] <- beta.w[1L, ]

    # PI.W.tilde
    PI.W.tilde <- matrix(0, p.tilde, nrow(beta.w) - 1L)
    PI.W.tilde[ ov.y.idx1, ] <- t(beta.w[-1L, ])

    # Sigma.B.tilde
    Sigma.B.tilde <- matrix(0, p.tilde, p.tilde)
    Sigma.B.tilde[ ov.y.idx1, ov.y.idx1 ] <- sigma.b

    # INT.B.tilde
    INT.B.tilde <- matrix(0, p.tilde, 1L)
    INT.B.tilde[ ov.y.idx1, 1L ] <- beta.b[1L, ]

    # PI.B.tilde
    PI.B.tilde <- matrix(0, p.tilde, nrow(beta.b) - 1L)
    PI.B.tilde[ ov.y.idx1, ] <- t(beta.b[-1L, ])

    if(length(between.y.idx) > 0L) {
        INT.B.tilde[ between.y.idx, 1L ] <- beta.z[1L, ]
         PI.B.tilde[ between.y.idx,    ] <- t(beta.z[-1L, ])
        Sigma.B.tilde[ between.y.idx, between.y.idx ] <- sigma.zz
        Sigma.B.tilde[ ov.y.idx1,     between.y.idx ] <- sigma.yz
        Sigma.B.tilde[ between.y.idx, ov.y.idx1     ] <- t(sigma.yz)
    }

    Res.Sigma.W  <- Sigma.W.tilde[ ov.y.idx1, ov.y.idx1, drop = FALSE]
    Res.Int.W    <-   INT.W.tilde[ ov.y.idx1,          , drop = FALSE]
    Res.Pi.W     <-    PI.W.tilde[ ov.y.idx1,          , drop = FALSE]

    Res.Sigma.B  <- Sigma.B.tilde[ ov.y.idx2, ov.y.idx2, drop = FALSE]
    Res.Int.B    <-   INT.B.tilde[ ov.y.idx2,          , drop = FALSE]
    Res.Pi.B     <-    PI.B.tilde[ ov.y.idx2,          , drop = FALSE]

    implied <- list(res.cov    = list(Res.Sigma.W, Res.Sigma.B),
                    res.int    = list(Res.Int.W, Res.Int.B),
                    res.slopes = list(Res.Pi.W, Res.Pi.B))

    # Note: cov.x and mean.x must be added by the caller
    implied
}

lav_mvreg_cluster_loglik_samplestats_2l <- function(YLp          = NULL,
                                                    Lp           = NULL,
                                                    Res.Sigma.W  = NULL,
                                                    Res.Int.W    = NULL,
                                                    Res.Pi.W     = NULL,
                                                    Res.Sigma.B  = NULL,
                                                    Res.Int.B    = NULL,
                                                    Res.Pi.B     = NULL,
                                                    out          = NULL, # 2l
                                                    Sinv.method  = "eigen",
                                                    log2pi       = FALSE,
                                                    minus.two    = TRUE) {

    # map implied to 2l matrices
    if(is.null(out)) {
        out <- lav_mvreg_cluster_implied22l(Lp = Lp, implied = NULL,
                   Res.Sigma.W = Res.Sigma.W,
                   Res.Int.W = Res.Int.W, Res.Pi.W = Res.Pi.W,
                   Res.Sigma.B = Res.Sigma.B,
                   Res.Int.B = Res.Int.B, Res.Pi.B = Res.Pi.B)
    }
    sigma.w  <- out$sigma.w
    sigma.b  <- out$sigma.b
    sigma.zz <- out$sigma.zz
    sigma.yz <- out$sigma.yz
    beta.w   <- out$beta.w
    beta.b   <- out$beta.b
    beta.z   <- out$beta.z
    beta.wb  <- out$beta.wb

    # check for beta.wb
    if(is.null(out$beta.wb)) {
        beta.wb <- rbind(beta.w, beta.b[-1,,drop = FALSE])
        beta.wb[1,] <- beta.wb[1,,drop = FALSE] + beta.b[1,,drop = FALSE]
    }

    # log 2*pi
    LOG.2PI <- log(2 * pi)

    # Lp
    nclusters      <- Lp$nclusters[[2]]
    cluster.size   <- Lp$cluster.size[[2]]
    cluster.sizes  <- Lp$cluster.sizes[[2]]
    ncluster.sizes <- Lp$ncluster.sizes[[2]]
    n.s            <- Lp$cluster.size.ns[[2]]

    # dependent 'y' level-2 ('Z') only variables?
    between.y.idx  <- Lp$between.y.idx[[2]]

    # extract (the many) sample statistics from YLp
    sample.wb         <- YLp[[2]]$sample.wb
    sample.YYres.wb1  <- YLp[[2]]$sample.YYres.wb1
    sample.XX.wb1     <- YLp[[2]]$sample.XX.wb1
    sample.wb2        <- YLp[[2]]$sample.wb2
    sample.YYres.wb2  <- YLp[[2]]$sample.YYres.wb2
    sample.YresX.wb2  <- YLp[[2]]$sample.YresX.wb2
    sample.XX.wb2     <- YLp[[2]]$sample.XX.wb2
    sample.clz.Y2.res <- YLp[[2]]$sample.clz.Y2.res
    sample.clz.Y2.XX  <- YLp[[2]]$sample.clz.Y2.XX
    sample.clz.Y2.B   <- YLp[[2]]$sample.clz.Y2.B
    if(length(between.y.idx) > 0L) {
        sample.clz.ZZ.res <- YLp[[2]]$sample.clz.ZZ.res
        sample.clz.ZZ.XX  <- YLp[[2]]$sample.clz.ZZ.XX
        sample.clz.ZZ.B   <- YLp[[2]]$sample.clz.ZZ.B
        sample.clz.YZ.res <- YLp[[2]]$sample.clz.YZ.res
        sample.clz.YZ.XX  <- YLp[[2]]$sample.clz.YZ.XX
        sample.clz.YresXZ <- YLp[[2]]$sample.clz.YresXZ  # zero?
        sample.clz.XWZres <- YLp[[2]]$sample.clz.XWZres
    }

    # reconstruct S.PW
    wb1.diff <- sample.wb - beta.wb
    Y1Y1.wb.res <- ( sample.YYres.wb1 +
                     t(wb1.diff) %*% sample.XX.wb1 %*% (wb1.diff) )

    # this one is weighted -- not the same as crossprod(Y2w.res)
    wb2.diff <- sample.wb2 - beta.wb
    Y2Y2w.res <- ( sample.YYres.wb2 +
                   sample.YresX.wb2 %*% (wb2.diff) +
                   t(wb2.diff) %*% t(sample.YresX.wb2) +
                   t(wb2.diff) %*% sample.XX.wb2 %*% (wb2.diff) )
    S.PW <- (Y1Y1.wb.res - Y2Y2w.res) / sum(cluster.size - 1)

    # common parts:
    sigma.w.inv <- lav_matrix_symmetric_inverse(S = sigma.w,
                       logdet = TRUE)
    sigma.w.logdet <- attr(sigma.w.inv, "logdet")
    if(length(between.y.idx) > 0L) {
        sigma.zz.inv <- lav_matrix_symmetric_inverse(S = sigma.zz,
                           logdet = TRUE)
        sigma.zz.logdet <- attr(sigma.zz.inv, "logdet")
        sigma.yz.zi <- sigma.yz %*% sigma.zz.inv
        sigma.zi.zy <- t(sigma.yz.zi)
        sigma.b.z <- sigma.b - sigma.yz %*% sigma.zi.zy
    } else {
        sigma.b.z <- sigma.b
    }

    # min 2* logliklihood
    DIST   <- numeric(ncluster.sizes)
    LOGDET <- numeric(ncluster.sizes)
    CONST  <- numeric(ncluster.sizes)
    for(clz in seq_len(ncluster.sizes)) {
        # cluster size
        nj <- cluster.sizes[clz]

        # data between
        nj.idx <- which(cluster.size == nj)
        y2.diff <- sample.clz.Y2.B[[clz]] - beta.wb
        Y2Yc.yy <- ( sample.clz.Y2.res[[clz]] +
                     t(y2.diff) %*% sample.clz.Y2.XX[[clz]] %*% (y2.diff) )
        if(length(between.y.idx) > 0L) {
            zz.diff <- sample.clz.ZZ.B[[clz]] - beta.z
            Y2Yc.zz <- ( sample.clz.ZZ.res[[clz]] +
                         t(zz.diff) %*% sample.clz.ZZ.XX[[clz]] %*% (zz.diff) )
            Y2Yc.yz <- ( sample.clz.YZ.res[[clz]] +
                         sample.clz.YresXZ[[clz]] %*% zz.diff + # zero?
                         t(y2.diff) %*% sample.clz.XWZres[[clz]] +
                         t(y2.diff) %*% sample.clz.YZ.XX[[clz]] %*% zz.diff )
        }

        # construct sigma.j
        sigma.j <- (nj * sigma.b.z) + sigma.w
        sigma.j.inv <- lav_matrix_symmetric_inverse(S = sigma.j,
                       logdet = TRUE)
        sigma.j.logdet <- attr(sigma.j.inv, "logdet")

        if(length(between.y.idx) > 0L) {
            sigma.ji.yz.zi <- sigma.j.inv %*% sigma.yz.zi

            # part 1 -- zz
            Vinv.11 <- sigma.zz.inv + nj*(sigma.zi.zy %*% sigma.ji.yz.zi)
            q.zz <- sum(Vinv.11 * Y2Yc.zz)

            # part 2 -- yz
            q.yz <- - nj * sum(sigma.ji.yz.zi * Y2Yc.yz)
        } else {
            q.zz <- q.yz <- sigma.zz.logdet <- 0
        }

        # part 5 -- yyc
        q.yyc <- -nj * sum(sigma.j.inv * Y2Yc.yy )

        if(log2pi) {
            P <- nj*nrow(sigma.w) + nrow(sigma.zz)
            CONST[clz] <- P * LOG.2PI
        }
        LOGDET[clz] <- sigma.zz.logdet + sigma.j.logdet
        DIST[clz] <-  q.zz + 2*q.yz - q.yyc
    }
    # q.yya + q.yyb
    q.W  <- sum(cluster.size - 1) * sum(sigma.w.inv * S.PW)
    # logdet within part
    L.W  <- sum(cluster.size - 1) * sigma.w.logdet

    # -2*times logl (without the constant) (for optimization)
    loglik <- sum(LOGDET*n.s) + sum(DIST) + q.W + L.W

    if(log2pi) {
        loglik <- loglik + sum(CONST*n.s)
    }

    # functions below compute -2 * logl
    if(!minus.two) {
        loglik <- loglik / (-2)
    }

    loglik
}

# first derivative -2*logl wrt Beta.W, Beta.B, Sigma.W, Sigma.B
lav_mvreg_cluster_dlogl_2l_samplestats <- function(YLp          = NULL,
                                                   Lp           = NULL,
                                                   Res.Sigma.W  = NULL,
                                                   Res.Int.W    = NULL,
                                                   Res.Pi.W     = NULL,
                                                   Res.Sigma.B  = NULL,
                                                   Res.Int.B    = NULL,
                                                   Res.Pi.B     = NULL,
                                                   out          = NULL, # 2l
                                                   return.list  = FALSE,
                                                   Sinv.method  = "eigen") {

    # map implied to 2l matrices
    if(is.null(out)) {
        out <- lav_mvreg_cluster_implied22l(Lp = Lp, implied = NULL,
                   Res.Sigma.W = Res.Sigma.W,
                   Res.Int.W = Res.Int.W, Res.Pi.W = Res.Pi.W,
                   Res.Sigma.B = Res.Sigma.B,
                   Res.Int.B = Res.Int.B, Res.Pi.B = Res.Pi.B)
    }
    sigma.w  <- out$sigma.w
    sigma.b  <- out$sigma.b
    sigma.zz <- out$sigma.zz
    sigma.yz <- out$sigma.yz
    beta.w   <- out$beta.w
    beta.b   <- out$beta.b
    beta.z   <- out$beta.z
    beta.wb  <- out$beta.wb

    # check for beta.wb
    if(is.null(out$beta.wb)) {
        beta.wb <- rbind(beta.w, beta.b[-1,,drop = FALSE])
        beta.wb[1,] <- beta.wb[1,,drop = FALSE] + beta.b[1,,drop = FALSE]
    }

    # Lp
    nclusters      <- Lp$nclusters[[2]]
    cluster.size   <- Lp$cluster.size[[2]]
    cluster.sizes  <- Lp$cluster.sizes[[2]]
    ncluster.sizes <- Lp$ncluster.sizes[[2]]
    n.s            <- Lp$cluster.size.ns[[2]]

    within.x.idx  <- Lp$within.x.idx[[1]]
    between.y.idx <- Lp$between.y.idx[[2]]

    w1.idx <- seq_len(length(within.x.idx) + 1L)
    b1.idx <- c(1L, seq_len(nrow(beta.wb))[-w1.idx])

    # extract (the many) sample statistics from YLp
    sample.wb         <- YLp[[2]]$sample.wb
    sample.YYres.wb1  <- YLp[[2]]$sample.YYres.wb1
    sample.XX.wb1     <- YLp[[2]]$sample.XX.wb1
    sample.wb2        <- YLp[[2]]$sample.wb2
    sample.YYres.wb2  <- YLp[[2]]$sample.YYres.wb2
    sample.YresX.wb2  <- YLp[[2]]$sample.YresX.wb2
    sample.XX.wb2     <- YLp[[2]]$sample.XX.wb2
    sample.clz.Y2.res <- YLp[[2]]$sample.clz.Y2.res
    sample.clz.Y2.XX  <- YLp[[2]]$sample.clz.Y2.XX
    sample.clz.Y2.B   <- YLp[[2]]$sample.clz.Y2.B
    if(length(between.y.idx) > 0L) {
        sample.clz.ZZ.res <- YLp[[2]]$sample.clz.ZZ.res
        sample.clz.ZZ.XX  <- YLp[[2]]$sample.clz.ZZ.XX
        sample.clz.ZZ.B   <- YLp[[2]]$sample.clz.ZZ.B
        sample.clz.YZ.res <- YLp[[2]]$sample.clz.YZ.res
        sample.clz.YZ.XX  <- YLp[[2]]$sample.clz.YZ.XX
        sample.clz.YresXZ <- YLp[[2]]$sample.clz.YresXZ  # zero?
        sample.clz.XWZres <- YLp[[2]]$sample.clz.XWZres
    }

    # reconstruct S.PW
    wb1.diff <- sample.wb - beta.wb
    Y1Y1.wb.res <- ( sample.YYres.wb1 +
                     t(wb1.diff) %*% sample.XX.wb1 %*% (wb1.diff) )

    # this one is weighted -- not the same as crossprod(Y2w.res)
    wb2.diff <- sample.wb2 - beta.wb
    Y2Y2w.res <- ( sample.YYres.wb2 +
                   sample.YresX.wb2 %*% (wb2.diff) +
                   t(wb2.diff) %*% t(sample.YresX.wb2) +
                   t(wb2.diff) %*% sample.XX.wb2 %*% (wb2.diff) )
    S.PW <- (Y1Y1.wb.res - Y2Y2w.res) / sum(cluster.size - 1)

    # common parts:
    sigma.w.inv <- lav_matrix_symmetric_inverse(S = sigma.w)

    G.beta.w   <- matrix(0, ncluster.sizes, length(beta.w))
    G.beta.b   <- matrix(0, ncluster.sizes, length(beta.b))
    G.beta.wb  <- matrix(0, ncluster.sizes, length(beta.wb))
    G.sigma.w1 <- matrix(0, ncluster.sizes, length(lav_matrix_vech(sigma.w)))
    G.sigma.b  <- matrix(0, ncluster.sizes, length(lav_matrix_vech(sigma.b)))

    if(length(between.y.idx) > 0L) {

        G.beta.z   <- matrix(0, ncluster.sizes, length(beta.z))
        G.sigma.zz <- matrix(0, ncluster.sizes, length(lav_matrix_vech(sigma.zz)))
        G.sigma.yz <- matrix(0, ncluster.sizes, length(sigma.yz))

        sigma.zz.inv <- lav_matrix_symmetric_inverse(S = sigma.zz)
        sigma.yz.zi <- sigma.yz %*% sigma.zz.inv
        sigma.zi.zy <- t(sigma.yz.zi)
        sigma.b.z <- sigma.b - sigma.yz %*% sigma.zi.zy

        for(clz in seq_len(ncluster.sizes)) {
            # cluster size
            nj <- cluster.sizes[clz]

            y2.diff <- sample.clz.Y2.B[[clz]] - beta.wb
            XX.y2.diff <- sample.clz.Y2.XX[[clz]] %*% y2.diff
            Y2Yc.yy <- sample.clz.Y2.res[[clz]] + crossprod(y2.diff, XX.y2.diff)

            zz.diff <- sample.clz.ZZ.B[[clz]] - beta.z
            Y2Yc.zz <- ( sample.clz.ZZ.res[[clz]] +
                         t(zz.diff) %*% sample.clz.ZZ.XX[[clz]] %*% (zz.diff) )
            Y2Yc.yz <- ( sample.clz.YZ.res[[clz]] +
                         sample.clz.YresXZ[[clz]] %*% zz.diff + # zero?
                         t(y2.diff) %*% sample.clz.XWZres[[clz]] +
                         t(y2.diff) %*% sample.clz.YZ.XX[[clz]] %*% zz.diff )

            # construct sigma.j
            sigma.j <- (nj * sigma.b.z) + sigma.w
            sigma.j.inv <- lav_matrix_symmetric_inverse(S = sigma.j)
            sigma.ji.yz.zi    <- sigma.j.inv %*% sigma.yz.zi
            sigma.zi.zy.ji    <- t(sigma.ji.yz.zi)
            sigma.ji.yz       <- sigma.j.inv %*% sigma.yz
            ns.sigma.j.inv    <- n.s[clz] * sigma.j.inv
            ns.sigma.zz.inv   <- n.s[clz] * sigma.zz.inv
            ns.sigma.yz       <- n.s[clz] * sigma.yz
            ns.sigma.ji.yz.zi <- n.s[clz] * sigma.ji.yz.zi

            # common parts
            ZZ.zi.yz.ji <- Y2Yc.zz %*% sigma.zi.zy.ji
            ji.YZ.zi    <-  sigma.j.inv %*% Y2Yc.yz %*% sigma.zz.inv

            jYZj.yy <- sigma.j.inv %*% Y2Yc.yy %*% sigma.j.inv
            jYZj.yz <- tcrossprod(ji.YZ.zi, sigma.ji.yz)
            jYZj.zz <- sigma.ji.yz.zi %*% ZZ.zi.yz.ji

            jYZj <- nj * (jYZj.yy + jYZj.zz - jYZj.yz - t(jYZj.yz))

            # SIGMA.W (between part)
            g.sigma.w1 <- ns.sigma.j.inv - jYZj
            tmp <- g.sigma.w1*2; diag(tmp) <- diag(g.sigma.w1)
            G.sigma.w1[clz,] <- lav_matrix_vech(tmp)

            # SIGMA.B
            g.sigma.b <- nj * g.sigma.w1
            tmp <- g.sigma.b*2; diag(tmp) <- diag(g.sigma.b)
            G.sigma.b[clz,] <- lav_matrix_vech(tmp)

            # SIGMA.ZZ
            YZ1 <- ZZ.zi.yz.ji %*% sigma.yz
            YZ2 <- crossprod(Y2Yc.yz, sigma.ji.yz)
            tmp <- (  t(sigma.yz) %*% g.sigma.w1 %*% sigma.yz
                     - 1/nj * Y2Yc.zz - t(YZ1) - YZ1 + t(YZ2) + YZ2 )
            g.sigma.zz <- ( ns.sigma.zz.inv +
                            nj * sigma.zz.inv %*% tmp %*% sigma.zz.inv )
            tmp <- g.sigma.zz*2; diag(tmp) <- diag(g.sigma.zz)
            G.sigma.zz[clz,] <- lav_matrix_vech(tmp)

            # SIGMA.ZY
            tmp1 <- crossprod(ZZ.zi.yz.ji, sigma.zz.inv)
            tmp2 <- ns.sigma.ji.yz.zi
            tmp3 <- ji.YZ.zi
            tmp4 <- jYZj %*% sigma.yz.zi
            g.sigma.yz <- 2 * nj * (tmp1 - tmp2 - tmp3 + tmp4)
            G.sigma.yz[clz,] <- lav_matrix_vec(g.sigma.yz)

            # BETA.Z
            A <- (sigma.zz.inv + nj*(sigma.zi.zy.ji %*% sigma.yz.zi)) # symm!
            B <- nj*(sigma.zi.zy.ji)
            tmp.z <- ( sample.clz.ZZ.XX[[clz]] %*% zz.diff %*% A -
                       (t(sample.clz.YresXZ[[clz]]) +
                        t(sample.clz.YZ.XX[[clz]]) %*% y2.diff) %*% t(B) )
            G.beta.z[clz, ] <- as.vector(-2 * tmp.z)

            # BETA.W (between part only) + BETA.B
            tmp <- ( sample.clz.XWZres[[clz]] +
                     sample.clz.YZ.XX[[clz]] %*% zz.diff )
            out.b <- tmp %*% sigma.zi.zy.ji - XX.y2.diff %*% sigma.j.inv
            out.w <- out.b + XX.y2.diff %*% sigma.w.inv
            tmp.b <- out.b[b1.idx,,drop = FALSE]
            tmp.w <- out.w[w1.idx,,drop = FALSE]
            G.beta.b[clz,] <- as.vector(2 * nj * tmp.b)
            G.beta.w[clz,] <- as.vector(2 * nj * tmp.w)

        } # clz

        d.beta.w1  <- matrix(colSums(G.beta.w), nrow(beta.w), ncol(beta.w))
        d.beta.b   <- matrix(colSums(G.beta.b), nrow(beta.b), ncol(beta.b))
        d.sigma.w1 <- lav_matrix_vech_reverse(colSums(G.sigma.w1))
        d.sigma.b  <- lav_matrix_vech_reverse(colSums(G.sigma.b))

        # z
        d.beta.z   <- matrix(colSums(G.beta.z), nrow(beta.z), ncol(beta.z))
        d.sigma.zz <- lav_matrix_vech_reverse(colSums(G.sigma.zz))
        d.sigma.yz <- matrix(colSums(G.sigma.yz), nrow(sigma.yz),ncol(sigma.yz))

    } # between.y.idx

    else { # no beween.y.idx

        for(clz in seq_len(ncluster.sizes)) {
            # cluster size
            nj <- cluster.sizes[clz]

            y2.diff    <- sample.clz.Y2.B[[clz]] - beta.wb
            XX.y2.diff <- sample.clz.Y2.XX[[clz]] %*% y2.diff
            Y2Yc.yy <- sample.clz.Y2.res[[clz]] + crossprod(y2.diff, XX.y2.diff)

            # construct sigma.j
            sigma.j <- (nj * sigma.b) + sigma.w
            sigma.j.inv <- lav_matrix_symmetric_inverse(S = sigma.j)

            # common part
            jYYj <- nj * sigma.j.inv %*% Y2Yc.yy %*% sigma.j.inv

            # SIGMA.W (between part)
            g.sigma.w1 <- (n.s[clz]*sigma.j.inv) - jYYj
            tmp <- g.sigma.w1*2; diag(tmp) <- diag(g.sigma.w1)
            G.sigma.w1[clz,] <- lav_matrix_vech(tmp)

            # SIGMA.B
            g.sigma.b <- nj * g.sigma.w1
            tmp <- g.sigma.b*2; diag(tmp) <- diag(g.sigma.b)
            G.sigma.b[clz,] <- lav_matrix_vech(tmp)

            # BETA.W (between part only) + BETA.B
            out.b <- -1 * XX.y2.diff %*% sigma.j.inv
            out.w <- out.b + XX.y2.diff %*% sigma.w.inv
            tmp.b <- out.b[b1.idx,,drop = FALSE]
            tmp.w <- out.w[w1.idx,,drop = FALSE]
            G.beta.b[clz,] <- as.vector(2 * nj * tmp.b)
            G.beta.w[clz,] <- as.vector(2 * nj * tmp.w)

        } # cl

        d.beta.w1  <- matrix(colSums(G.beta.w), nrow(beta.w), ncol(beta.w))
        d.beta.b   <- matrix(colSums(G.beta.b), nrow(beta.b), ncol(beta.b))
        d.sigma.w1 <- lav_matrix_vech_reverse(colSums(G.sigma.w1))
        d.sigma.b  <- lav_matrix_vech_reverse(colSums(G.sigma.b))

        # z
        d.beta.z   <- matrix(0, 0L, 0L)
        d.sigma.zz <- matrix(0, 0L, 0L)
        d.sigma.yz <- matrix(0, 0L, 0L)

    } # no-between-y

    # Sigma.W (bis)
    d.sigma.w2 <- sum(cluster.size - 1) * ( sigma.w.inv
                   - sigma.w.inv %*% S.PW %*% sigma.w.inv )
    tmp <- d.sigma.w2*2; diag(tmp) <- diag(d.sigma.w2)
    d.sigma.w2 <- tmp

    d.sigma.w <- d.sigma.w1 + d.sigma.w2

    # beta.w (bis)
    d.beta.w2 <- -2 * (sample.XX.wb1 %*% (sample.wb - beta.wb))[w1.idx,,drop = FALSE] %*% sigma.w.inv

    d.beta.w <- d.beta.w1 + d.beta.w2

    # rearrange
    dimplied <- lav_mvreg_cluster_2l2implied(Lp,
                    sigma.w = d.sigma.w, sigma.b = d.sigma.b,
                    sigma.zz = d.sigma.zz, sigma.yz = d.sigma.yz,
                    beta.w = d.beta.w, beta.b = d.beta.b, beta.z = d.beta.z)

    if(return.list) {
        return(dimplied)
    }

    # as a single vector
    out <- c(drop(dimplied$res.int[[1]]),
             lav_matrix_vec(dimplied$res.slopes[[1]]),
             lav_matrix_vech(dimplied$res.cov[[1]]),
             drop(dimplied$res.int[[2]]),
             lav_matrix_vec(dimplied$res.slopes[[2]]),
             lav_matrix_vech(dimplied$res.cov[[2]]))
    out
}

# cluster-wise scores -2*logl wrt Beta.W, Beta.B, Sigma.W, Sigma.B
lav_mvreg_cluster_scores_2l <- function(Y1           = NULL,
                                        YLp          = NULL,
                                        Lp           = NULL,
                                        Res.Sigma.W  = NULL,
                                        Res.Int.W    = NULL,
                                        Res.Pi.W     = NULL,
                                        Res.Sigma.B  = NULL,
                                        Res.Int.B    = NULL,
                                        Res.Pi.B     = NULL,
                                        out          = NULL, # 2l
                                        Sinv.method  = "eigen") {

    # map implied to 2l matrices
    if(is.null(out)) {
        out <- lav_mvreg_cluster_implied22l(Lp = Lp, implied = NULL,
                   Res.Sigma.W = Res.Sigma.W,
                   Res.Int.W = Res.Int.W, Res.Pi.W = Res.Pi.W,
                   Res.Sigma.B = Res.Sigma.B,
                   Res.Int.B = Res.Int.B, Res.Pi.B = Res.Pi.B)
    }
    sigma.w  <- out$sigma.w
    sigma.b  <- out$sigma.b
    sigma.zz <- out$sigma.zz
    sigma.yz <- out$sigma.yz
    beta.w   <- out$beta.w
    beta.b   <- out$beta.b
    beta.z   <- out$beta.z
    beta.wb  <- out$beta.wb

    # check for beta.wb
    if(is.null(out$beta.wb)) {
        beta.wb <- rbind(beta.w, beta.b[-1,,drop = FALSE])
        beta.wb[1,] <- beta.wb[1,,drop = FALSE] + beta.b[1,,drop = FALSE]
    }

    # Lp
    nclusters      <- Lp$nclusters[[2]]
    cluster.size   <- Lp$cluster.size[[2]]
    cluster.idx     <- Lp$cluster.idx[[2]]

    within.x.idx  <- Lp$within.x.idx[[1]]
    between.idx   <- Lp$between.idx[[2]]
    between.y.idx <- Lp$between.y.idx[[2]]
    between.x.idx <- Lp$between.x.idx[[2]]

    y1.idx <- Lp$ov.y.idx[[1]]
    x1.idx <- c(within.x.idx, between.x.idx) # in that order

    # residuals for 'Y'
    Y1.wb <- Y1[, y1.idx, drop = FALSE]

    if(length(x1.idx) > 0L) {
        EXO.wb <- cbind(1, Y1[, x1.idx, drop = FALSE])
        Y1.wb.hat <- EXO.wb %*% beta.wb
        Y1.wb.res <- Y1.wb - Y1.wb.hat
    } else {
        Y1.wb.res <- Y1.wb
    }

    # residuals 'Y' (level 2)
    Y2 <- YLp[[2]]$Y2
    if(length(x1.idx) > 0L) {
        EXO.wb2 <- cbind(1, Y2[, x1.idx, drop = FALSE])
        Y2w.res <- Y2[, y1.idx, drop = FALSE] - EXO.wb2 %*% beta.wb
    } else {
        EXO.wb2 <- matrix(1, nrow(Y2), 1L)
        Y2w.res <- Y2[, y1.idx, drop = FALSE]
    }

    # residual 'Z' (level 2)
    if(length(between.y.idx) > 0L) {
        if(length(between.x.idx) > 0L) {
            EXO.z <- cbind(1, Y2[, between.x.idx, drop = FALSE])
            Y2.z <- Y2[, between.y.idx, drop = FALSE]
            Y2z.res <-  Y2.z - EXO.z %*% beta.z
            # sample.z
            #XX.z <- crossprod(EXO.z)
            #sample.z <- try(solve(XX.z, crossprod(EXO.z, Y2.z)))
            #if(inherits(sample.z, "try-error")) {
            #    sample.z <- MASS::ginv(XX.z) %*% crossprod(EXO.z, Y2.z)
            #}

            # sample.wb2
            #sample.wb2 <- YLp[[2]]$sample.wb2
        } else {
            Y2z.res <- Y2[, between.y.idx, drop = FALSE]
        }
    }

    # common parts:
    sigma.w.inv <- lav_matrix_symmetric_inverse(S = sigma.w)

    G.beta.w1  <- matrix(0, nclusters, length(beta.w))
    G.beta.b   <- matrix(0, nclusters, length(beta.b))
    G.beta.wb  <- matrix(0, nclusters, length(beta.wb))
    G.sigma.w1 <- matrix(0, nclusters, length(lav_matrix_vech(sigma.w)))
    G.sigma.b  <- matrix(0, nclusters, length(lav_matrix_vech(sigma.b)))


    if(length(between.y.idx) > 0L) {

        G.beta.z   <- matrix(0, nclusters, length(beta.z))
        G.sigma.zz <- matrix(0, nclusters, length(lav_matrix_vech(sigma.zz)))
        G.sigma.yz <- matrix(0, nclusters, length(sigma.yz))

        sigma.zz.inv <- lav_matrix_symmetric_inverse(S = sigma.zz)
        sigma.yz.zi <- sigma.yz %*% sigma.zz.inv
        sigma.zi.zy <- t(sigma.yz.zi)
        sigma.b.z <- sigma.b - sigma.yz %*% sigma.zi.zy

        for(cl in seq_len(nclusters)) {
            # cluster size
            nj <- cluster.size[cl]

            # data within for the cluster (centered)
            Y1m <- Y1.wb.res[cluster.idx == cl,, drop = FALSE]
            yc  <- Y2w.res[cl,]

            # data between
            zc  <- Y2z.res[cl,]
            Y2Yc.yy <- tcrossprod(Y2w.res[cl,])
            Y2Yc.zz <- tcrossprod(Y2z.res[cl,])
            Y2Yc.yz <- tcrossprod(Y2w.res[cl,], Y2z.res[cl,])

            # construct sigma.j
            sigma.j <- (nj * sigma.b.z) + sigma.w
            sigma.j.inv <- lav_matrix_symmetric_inverse(S = sigma.j)
            sigma.ji.yz.zi    <- sigma.j.inv %*% sigma.yz.zi
            sigma.zi.zy.ji    <- t(sigma.ji.yz.zi)
            sigma.ji.yz       <- sigma.j.inv %*% sigma.yz

            # common parts
            ZZ.zi.yz.ji <- Y2Yc.zz %*% sigma.zi.zy.ji
            ji.YZ.zi    <-  sigma.j.inv %*% Y2Yc.yz %*% sigma.zz.inv

            jYZj.yy <- sigma.j.inv %*% Y2Yc.yy %*% sigma.j.inv
            jYZj.yz <- tcrossprod(ji.YZ.zi, sigma.ji.yz)
            jYZj.zz <- sigma.ji.yz.zi %*% ZZ.zi.yz.ji

            jYZj <- nj * (jYZj.yy + jYZj.zz - jYZj.yz - t(jYZj.yz))

            # SIGMA.W (between part)
            g.sigma.w1 <- sigma.j.inv - jYZj
            tmp <- g.sigma.w1*2; diag(tmp) <- diag(g.sigma.w1)
            G.sigma.w1[cl,] <- lav_matrix_vech(tmp)

            # SIGMA.W (within part)
            #g.sigma.w2 <- ( (nj-1) * sigma.w.inv
            #    - sigma.w.inv %*% (crossprod(Y1m) - nj*Y2Yc.yy) %*% sigma.w.inv )
            #tmp <- g.sigma.w2*2; diag(tmp) <- diag(g.sigma.w2)
            #G.sigma.w2[cl,] <- lav_matrix_vech(tmp)
            #G.sigma.w[cl,] <- G.sigma.w1[cl,] + G.sigma.w2[cl,]

            # SIGMA.B
            g.sigma.b <- nj * g.sigma.w1
            tmp <- g.sigma.b*2; diag(tmp) <- diag(g.sigma.b)
            G.sigma.b[cl,] <- lav_matrix_vech(tmp)

            # SIGMA.ZZ
            YZ1 <- ZZ.zi.yz.ji %*% sigma.yz
            YZ2 <- crossprod(Y2Yc.yz, sigma.ji.yz)
            tmp <- ( t(sigma.yz) %*% g.sigma.w1 %*% sigma.yz
                     -(1/nj * Y2Yc.zz + t(YZ1) + YZ1 - t(YZ2) - YZ2) )
            g.sigma.zz <- ( sigma.zz.inv +
                            nj * sigma.zz.inv %*% tmp %*% sigma.zz.inv )
            tmp <- g.sigma.zz*2; diag(tmp) <- diag(g.sigma.zz)
            G.sigma.zz[cl,] <- lav_matrix_vech(tmp)

            # SIGMA.ZY
            #g.sigma.yz <- 2 * nj * (
            #              (sigma.j.inv %*%
            #                  (sigma.yz.zi %*% Y2Yc.zz - sigma.yz - Y2Yc.yz)
            #                   + jYZj %*% sigma.yz) %*% sigma.zz.inv )
            tmp1 <- crossprod(ZZ.zi.yz.ji, sigma.zz.inv)
            tmp2 <- sigma.ji.yz.zi
            tmp3 <- ji.YZ.zi
            tmp4 <- jYZj %*% sigma.yz.zi
            g.sigma.yz <- 2 * nj * (tmp1 - tmp2 - tmp3 + tmp4)
            G.sigma.yz[cl,] <- lav_matrix_vec(g.sigma.yz)

            # BETA.Z
            # here, we avoid the (sample.z - beta.z) approach
            exo.z <- cbind(1, Y2[cl, between.x.idx, drop = FALSE])
            tmp1 <- (sigma.zz.inv + nj*(sigma.zi.zy.ji %*% sigma.yz.zi)) %*% zc
            tmp2 <- nj * (sigma.zi.zy.ji) %*% yc
            tmp.z <- crossprod(exo.z, drop(tmp1 - tmp2))
            G.beta.z[cl, ] <- as.vector(-2 * tmp.z)

            # BETA.W
#            exo.w <- cbind(1,
#                           Y1[cluster.idx == cl, within.x.idx, drop = FALSE])
#            G.beta.w[cl,] <- as.vector( 2 * t(exo.w) %*% (
#                                 matrix(1, nj, 1) %x% (zc %*% sigma.zi.zy.ji -
#                                                       yc %*% sigma.j.inv +
#                                                       yc %*% sigma.w.inv) -
#                                 Y1m %*% sigma.w.inv) )

            # BETA.W (between part only)
            exo2.w <- cbind(1, Y2[cl, within.x.idx, drop = FALSE])
            tmp2 <- (zc %*% sigma.zi.zy.ji -
                     yc %*% sigma.j.inv +
                     yc %*% sigma.w.inv)
            G.beta.w1[cl,] <- as.vector( 2 * nj * crossprod(exo2.w, tmp2) )

            # BETA.W (within part only)
            #exo.w <- cbind(1,
            #               Y1[cluster.idx == cl, within.x.idx, drop = FALSE])
            #tmp1 <- - Y1m %*% sigma.w.inv
            #G.beta.ww <- as.vector( 2 * crossprod(exo.w, tmp1) )
            #G.beta.w[cl,] <- G.beta.w1 + G.beta.ww
            #G.beta.w2[cl,] <- G.beta.ww

            # BETA.B
            exo2.b <- cbind(1, Y2[cl, between.x.idx, drop = FALSE])
            tmp <- ( zc %*% sigma.zi.zy.ji - yc %*% sigma.j.inv )
            G.beta.b[cl,] <- as.vector(  2 * nj * crossprod(exo2.b, tmp) )
        } # cl

    } # between.y.idx

    else { # no beween.y.idx

        for(cl in seq_len(nclusters)) {
            # cluster size
            nj <- cluster.size[cl]

            # data within for the cluster (centered)
            Y1m <- Y1.wb.res[cluster.idx == cl,, drop = FALSE]
            yc  <- Y2w.res[cl,]

            # data between
            Y2Yc.yy <- tcrossprod(Y2w.res[cl,])

            # construct sigma.j
            sigma.j <- (nj * sigma.b) + sigma.w
            sigma.j.inv <- lav_matrix_symmetric_inverse(S = sigma.j)

            # common part
            jYYj <- nj * sigma.j.inv %*% Y2Yc.yy %*% sigma.j.inv

            # SIGMA.W
            #g.sigma.w <- ( (nj-1) * sigma.w.inv
            #    - sigma.w.inv %*% (crossprod(Y1m) - nj*Y2Yc.yy) %*% sigma.w.inv
            #    + sigma.j.inv - jYYj )
            #tmp <- g.sigma.w*2; diag(tmp) <- diag(g.sigma.w)
            #G.sigma.w[cl,] <- lav_matrix_vech(tmp)

            # SIGMA.W (between part)
            g.sigma.w1 <- sigma.j.inv - jYYj
            tmp <- g.sigma.w1*2; diag(tmp) <- diag(g.sigma.w1)
            G.sigma.w1[cl,] <- lav_matrix_vech(tmp)

            # SIGMA.B
            g.sigma.b <- nj * (sigma.j.inv - jYYj)
            tmp <- g.sigma.b*2; diag(tmp) <- diag(g.sigma.b)
            G.sigma.b[cl,] <- lav_matrix_vech(tmp)

            # BETA.W (between part only)
            exo2.w <- cbind(1, Y2[cl, within.x.idx, drop = FALSE])
            tmp2 <- ( - yc %*% sigma.j.inv + yc %*% sigma.w.inv )
            G.beta.w1[cl,] <- as.vector( 2 * nj * crossprod(exo2.w, tmp2) )

            # BETA.B
            exo2.b <- cbind(1, Y2[cl, between.x.idx, drop = FALSE])
            tmp <- - yc %*% sigma.j.inv
            G.beta.b[cl,] <- as.vector(  2 * nj * crossprod(exo2.b, tmp) )
        } # cl

    } # no-between-y

    # beta.w (bis)

#    d.beta.w2 <- -2 * t(EXO.wb[,1:(length(within.x.idx) + 1L), drop = FALSE]) %*% Y1.wb.res %*% sigma.w.inv

    Y1.wb.res.i <- Y1.wb.res %*% sigma.w.inv
    w1.idx <- seq_len(length(within.x.idx) + 1L)
    a1.idx <- rep(w1.idx, times = ncol(Y1.wb.res.i))
    b1.idx <- rep(seq_len(ncol(Y1.wb.res.i)), each = length(w1.idx))
    TMP <- EXO.wb[,a1.idx,drop = FALSE] * Y1.wb.res.i[,b1.idx,drop = FALSE]
    G.beta.w2 <- -2 * rowsum.default(TMP, cluster.idx, reorder = FALSE,
                                     na.rm = TRUE)
    G.beta.w <- G.beta.w1 + G.beta.w2

    # Sigma.W (bis)
    #d.sigma.w2 <- sum(cluster.size - 1) * ( sigma.w.inv
    #               - sigma.w.inv %*% S.PW %*% sigma.w.inv )
    #tmp <- d.sigma.w2*2; diag(tmp) <- diag(d.sigma.w2)
    #d.sigma.w2 <- tmp

    #g.sigma.w2 <- ( (nj-1) * sigma.w.inv
    #        - sigma.w.inv %*% (crossprod(Y1m) - nj*Y2Yc.yy) %*% sigma.w.inv )

    Y1a.res <- Y1.wb.res - Y2w.res[cluster.idx, , drop = FALSE]
    Y1a.res.i <- Y1a.res %*% sigma.w.inv
    idx1 <- lav_matrix_vech_col_idx(nrow(sigma.w))
    idx2 <- lav_matrix_vech_row_idx(nrow(sigma.w))
    SW2 <- matrix(lav_matrix_vech(sigma.w.inv), nrow = nclusters,
                  length(lav_matrix_vech(sigma.w.inv)), byrow = TRUE)
    SW2 <- SW2 * (cluster.size - 1)
    TMP <- Y1a.res.i[, idx1, drop = FALSE] * Y1a.res.i[, idx2, drop = FALSE]
    TMP2 <- rowsum.default(TMP, cluster.idx, reorder = FALSE, na.rm = TRUE)
    G.sigma.w2 <- 2 * (SW2 - TMP2)
    diagh.idx <- lav_matrix_diagh_idx(nrow(sigma.w))
    G.sigma.w2[,diagh.idx] <- G.sigma.w2[,diagh.idx,drop = FALSE]/2
    G.sigma.w <- G.sigma.w1 + G.sigma.w2



    # rearrange columns to Res.Int.W, Res.Pi.W, Res.Sigma.W,
    #                      Res.Int.B, Res.Pi.B, Res.Sigma.B

    # ov.idx per level
    ov.idx      <- Lp$ov.idx

    # 'tilde' matrices: ALL variables within and between
    p.tilde <- length( unique(c(ov.idx[[1]], ov.idx[[2]])) )
    p.tilde.star <- p.tilde * (p.tilde + 1) / 2
    B.tilde <- lav_matrix_vech_reverse(seq_len(p.tilde.star))

    # only 'y'
    ov.y.idx <- Lp$ov.y.idx

    # two levels only (for now)
    ov.y.idx1 <- ov.y.idx[[1]]
    ov.y.idx2 <- ov.y.idx[[2]]

    # WITHIN (is easy)
    BETA.W.idx <- matrix(seq_len(length(beta.w)), nrow(beta.w), ncol(beta.w))
    BETA.B.idx <- matrix(seq_len(length(beta.b)), nrow(beta.b), ncol(beta.b))

    Res.Int.W   <- G.beta.w[,BETA.W.idx[1L,],drop = FALSE]
    Res.Pi.W    <- G.beta.w[,lav_matrix_vecr(BETA.W.idx[-1L,]),drop = FALSE]
    Res.Sigma.W <- G.sigma.w

    # Sigma.B
    Sigma.B.tilde <- matrix(0, nclusters, p.tilde.star)
    col.idx <- lav_matrix_vech(B.tilde[ov.y.idx1, ov.y.idx1, drop = FALSE])
    Sigma.B.tilde[ , col.idx ] <- G.sigma.b

    # Int.B
    BETA.B.tilde <- matrix(seq_len(nrow(beta.b)*p.tilde), nrow(beta.b), p.tilde)
    Int.B <- matrix(0, nclusters, p.tilde)
    Int.B[,ov.y.idx1] <- G.beta.b[,BETA.B.idx[1L,]]

    # Pi.B
    Pi.B <- matrix(0, nclusters, p.tilde * (nrow(beta.b) - 1L))
    col.idx <- lav_matrix_vecr(BETA.B.tilde[-1L, ov.y.idx1, drop = FALSE])
    Pi.B[, col.idx] <- G.beta.b[, lav_matrix_vecr(BETA.B.idx[-1L,]),drop=FALSE]

    if(length(between.y.idx) > 0L) {
        # Sigma.B: add yz/zz parts
        col.idx <- lav_matrix_vec(B.tilde[ov.y.idx1,between.y.idx,drop = FALSE])
        Sigma.B.tilde[ , col.idx ] <- G.sigma.yz
        col.idx <- lav_matrix_vech(B.tilde[between.y.idx, between.y.idx,
                                           drop = FALSE ])
        Sigma.B.tilde[ , col.idx ] <- G.sigma.zz

        # Int.B: add z-part
        BETA.Z.idx <- matrix(seq_len(length(beta.z)),nrow(beta.z),ncol(beta.z))
        Int.B[,between.y.idx] <- G.beta.z[,BETA.Z.idx[1L,],drop = FALSE]

        # Pi.B: add beta.z
        col.idx <- lav_matrix_vecr(BETA.B.tilde[-1L,between.y.idx,drop = FALSE])
        Pi.B[, col.idx] <-
            G.beta.z[, lav_matrix_vecr(BETA.Z.idx[-1L,]), drop = FALSE]
    }

    # only extract ov.y.idx2 for BETWEEN
    col.idx <- lav_matrix_vech(B.tilde[ov.y.idx2, ov.y.idx2, drop = FALSE])
    Res.Sigma.B <- Sigma.B.tilde[ , col.idx, drop = FALSE ]

    Res.Int.B <- Int.B[, ov.y.idx2, drop = FALSE]

    col.idx <- lav_matrix_vecr(BETA.B.tilde[-1, ov.y.idx2])
    Res.Pi.B <- Pi.B[, col.idx, drop = FALSE]

    SCORES <- cbind(Res.Int.W, Res.Pi.W, Res.Sigma.W,
                    Res.Int.B, Res.Pi.B, Res.Sigma.B)

    SCORES
}

# first-order information: outer crossprod of scores per cluster
lav_mvreg_cluster_information_firstorder <- function(Y1           = NULL,
                                                     YLp          = NULL,
                                                     Lp           = NULL,
                                                     Res.Sigma.W  = NULL,
                                                     Res.Int.W    = NULL,
                                                     Res.Pi.W     = NULL,
                                                     Res.Sigma.B  = NULL,
                                                     Res.Int.B    = NULL,
                                                     Res.Pi.B     = NULL,
                                                     divide.by.two = FALSE,
                                                     Sinv.method  = "eigen") {
    N <- NROW(Y1)

    SCORES <- lav_mvreg_cluster_scores_2l(Y1           = Y1,
                                          YLp          = YLp,
                                          Lp           = Lp,
                                          Res.Sigma.W  = Res.Sigma.W,
                                          Res.Int.W    = Res.Int.W,
                                          Res.Pi.W     = Res.Pi.W,
                                          Res.Sigma.B  = Res.Sigma.B,
                                          Res.Int.B    = Res.Int.B,
                                          Res.Pi.B     = Res.Pi.B,
                                          Sinv.method  = Sinv.method)

    # divide by 2 (if we want scores wrt objective function)
    if(divide.by.two) {
        SCORES <- SCORES / 2
    }

    # unit information
    information <- crossprod(SCORES)/Lp$nclusters[[2]]

    information
}