R/lav_model_h1_information.R

In lavaan: Latent Variable Analysis

# the information matrix of the unrestricted (H1) model
# taking into account:
#   - the estimator (ML or (D)WLS/ULS)
#   - missing or not
#   - fixed.x = TRUE or FALSE
#   - conditional.x = TRUE or FALSE
#   - h1.information is "structured" or "unstructured"
#
# Note: this replaces the (old) lav_model_wls_v() function
#
# - YR 22 Okt 2017: initial version
# - YR 03 Dec 2017: add lavh1, implied is either lavimplied or lavh1
#                   add support for clustered data: first.order
# - YR 03 Jan 2018: add support for clustered data: expected
# - YR 23 Aug 2018: lav_model_h1_acov (0.6-3)

lav_model_h1_information <- function(lavobject      = NULL,
                                     lavmodel       = NULL,
                                     lavsamplestats = NULL,
                                     lavdata        = NULL,
                                     lavimplied     = NULL,
                                     lavh1          = NULL,
                                     lavcache       = NULL,
                                     lavoptions     = NULL) {

    if(!is.null(lavobject) && inherits(lavobject, "lavaan")) {
        lavmodel       <- lavobject@Model
        lavsamplestats <- lavobject@SampleStats
        lavdata        <- lavobject@Data
        lavimplied     <- lavobject@implied
        if(.hasSlot(lavobject, "h1")) {
            lavh1      <- lavobject@h1
        } else {
            lavh1      <- lav_h1_implied_logl(lavdata = lavobject@Data,
                                      lavsamplestats = lavobject@SampleStats,
                                      lavoptions = lavobject@Options)
        }
        lavcache       <- lavobject@Cache
        lavoptions     <- lavobject@Options
    }

    # sanity check
    if(length(lavh1) == 0L) {
        lavh1 <- lav_h1_implied_logl(lavdata = lavdata,
                                     lavsamplestats = lavsamplestats,
                                     lavoptions = lavoptions)
    }
    if(length(lavimplied) == 0L) {
        lavimplied <- lav_model_implied(lavmodel = lavmodel)
    }

    # information
    information <- lavoptions$information[1] # ALWAYS take the first one
                                             # the caller must control it


    # compute information matrix
    if(information == "observed") {
        I1 <- lav_model_h1_information_observed(lavmodel = lavmodel,
            lavsamplestats = lavsamplestats, lavdata = lavdata,
            lavimplied = lavimplied, lavh1 = lavh1,
            lavcache = lavcache, lavoptions = lavoptions)
    } else if(information == "expected") {
        I1 <- lav_model_h1_information_expected(lavmodel = lavmodel,
            lavsamplestats = lavsamplestats, lavdata = lavdata,
            lavimplied = lavimplied, lavh1 = lavh1,
            lavcache = lavcache, lavoptions = lavoptions)
    } else if(information == "first.order") {
        I1 <- lav_model_h1_information_firstorder(lavmodel = lavmodel,
            lavsamplestats = lavsamplestats, lavdata = lavdata,
            lavimplied = lavimplied, lavh1 = lavh1,
            lavcache = lavcache, lavoptions = lavoptions)
    }

    # I1 information, as a list per group
    I1
}

# fisher/expected information of H1
lav_model_h1_information_expected <- function(lavobject      = NULL,
                                              lavmodel       = NULL,
                                              lavsamplestats = NULL,
                                              lavdata        = NULL,
                                              lavoptions     = NULL,
                                              lavimplied     = NULL,
                                              lavh1          = NULL,
                                              lavcache       = NULL) {

    if(!is.null(lavobject) && inherits(lavobject, "lavaan")) {
        lavmodel       <- lavobject@Model
        lavsamplestats <- lavobject@SampleStats
        lavdata        <- lavobject@Data
        lavimplied     <- lavobject@implied
        if(.hasSlot(lavobject, "h1")) {
            lavh1      <- lavobject@h1
        } else {
            lavh1      <- lav_h1_implied_logl(lavdata = lavobject@Data,
                                      lavsamplestats = lavobject@SampleStats,
                                      lavoptions = lavobject@Options)
        }
        lavcache       <- lavobject@Cache
        lavoptions     <- lavobject@Options
    }

    # sanity check
    if(length(lavh1) == 0L) {
        lavh1 <- lav_h1_implied_logl(lavdata = lavdata,
                                     lavsamplestats = lavsamplestats,
                                     lavoptions = lavoptions)
    }
    if(length(lavimplied) == 0L) {
        lavimplied <- lav_model_implied(lavmodel = lavmodel)
    }

    estimator <- lavmodel@estimator

    # structured of unstructured? (since 0.5-23)
    if(!is.null(lavoptions) &&
       !is.null(lavoptions$h1.information[1]) &&
       lavoptions$h1.information[1] == "unstructured") {
        structured <- FALSE
    } else {
        structured <- TRUE
    }


    # 1. WLS.V (=A1) for GLS/WLS
    if(lavmodel@estimator == "GLS"  || lavmodel@estimator == "WLS") {
        A1 <- lavsamplestats@WLS.V
    }

    # 1b.
    else if(lavmodel@estimator == "DLS") {
        if(lavmodel@estimator.args$dls.GammaNT == "sample") {
            A1 <- lavsamplestats@WLS.V
        } else {
            A1 <- vector("list", length = lavsamplestats@ngroups)
            for(g in seq_len(lavsamplestats@ngroups)) {
                dls.a <- lavmodel@estimator.args$dls.a
                GammaNT <- lav_samplestats_Gamma_NT(
                               COV            = lavimplied$cov[[g]],
                               MEAN           = lavimplied$mean[[g]],
                               rescale        = FALSE,
                               x.idx          = lavsamplestats@x.idx[[g]],
                               fixed.x        = lavmodel@fixed.x,
                               conditional.x  = lavmodel@conditional.x,
                               meanstructure  = lavmodel@meanstructure,
                               slopestructure = lavmodel@conditional.x)
                W.DLS <- (1 - dls.a)*lavsamplestats@NACOV[[g]] + dls.a*GammaNT
                A1[[g]] <- lav_matrix_symmetric_inverse(W.DLS)
            }
        }
    }

    # 2. DWLS/ULS diagonal @WLS.VD slot
    else if(lavmodel@estimator == "DWLS"  || lavmodel@estimator == "ULS") {
        # diagonal only!!
        A1 <- lavsamplestats@WLS.VD
    }

    # 3a. ML single level
    else if( lavmodel@estimator %in% c("ML", "NTRLS", "DLS", "catML")
             && lavdata@nlevels == 1L ) {
        A1 <- vector("list", length=lavsamplestats@ngroups)

        # structured? compute model-implied statistics
        if(structured && length(lavimplied) == 0L) {
            lavimplied <- lav_model_implied(lavmodel)
        }

        for(g in 1:lavsamplestats@ngroups) {

            if(.hasSlot(lavdata, "weights")) {
                WT <- lavdata@weights[[g]]
            } else {
                WT <- NULL
            }

            if(lavsamplestats@missing.flag) {
                # mvnorm
                # FIXME: allow for meanstructure = FALSE
                # FIXME: allow for conditional.x = TRUE
                if(lavmodel@meanstructure && structured) {
                    MEAN <- lavimplied$mean[[g]]
                } else {
                    MEAN <- lavsamplestats@missing.h1[[g]]$mu
                }

                if(structured) {
                    A1[[g]] <-
                      lav_mvnorm_missing_information_expected(
                          Y = lavdata@X[[g]],
                          Mp = lavdata@Mp[[g]],
                          wt = WT,
                          Mu = MEAN,
                          # meanstructure = lavmodel@meanstructure,
                          Sigma = lavimplied$cov[[g]],
                          x.idx = lavsamplestats@x.idx[[g]])
                } else {
                    A1[[g]] <-
                      lav_mvnorm_missing_information_expected(
                          Y = lavdata@X[[g]],
                          Mp = lavdata@Mp[[g]],
                          wt = WT,
                          Mu = MEAN,
                          # meanstructure = lavmodel@meanstructure,
                          Sigma = lavsamplestats@missing.h1[[g]]$sigma,
                          x.idx = lavsamplestats@x.idx[[g]])
                }
            } else {
                if(lavmodel@conditional.x) {
                    # mvreg
                    if(lavmodel@meanstructure && structured) {
                        RES.INT    <- lavimplied$res.int[[g]]
                        RES.SLOPES <- lavimplied$res.slopes[[g]]
                    } else {
                        RES.INT    <- lavsamplestats@res.int[[g]]
                        RES.SLOPES <- lavsamplestats@res.slopes[[g]]
                    }

                    if(structured) {
                        A1[[g]] <- lav_mvreg_information_expected(
                            sample.mean.x     = lavsamplestats@mean.x[[g]],
                            sample.cov.x      = lavsamplestats@cov.x[[g]],
                            sample.nobs       = lavsamplestats@nobs[[g]],
                            res.int           = RES.INT,
                            res.slopes        = RES.SLOPES,
                            #wt               = WT,
                            #meanstructure    = lavmodel@meanstructure,
                            res.cov           = lavimplied$res.cov[[g]])
                    } else {
                        A1[[g]] <- lav_mvreg_information_expected(
                            sample.mean.x     = lavsamplestats@mean.x[[g]],
                            sample.cov.x      = lavsamplestats@cov.x[[g]],
                            sample.nobs       = lavsamplestats@nobs[[g]],
                            res.int           = lavsamplestats@res.int[[g]],
                            res.slopes        = lavsamplestats@res.slopes[[g]],
                            #wt               = WT,
                            #meanstructure    = lavmodel@meanstructure,
                            res.cov           = lavsamplestats@res.cov[[g]])
                    }

                } else {
                    # conditional.x = FALSE
                    # mvnorm
                    if(lavmodel@meanstructure && structured) {
                        MEAN <- lavimplied$mean[[g]]
                    } else {
                        MEAN <- lavsamplestats@mean[[g]]
                    }

                    correlation.flag <- FALSE
                    if(.hasSlot(lavmodel, "correlation")) {
                        correlation.flag <- lavmodel@correlation
                    }
                    if(structured) {
                        A1[[g]] <- lav_mvnorm_information_expected(
                                  Sigma         = lavimplied$cov[[g]],
                                  #wt = WT, # not needed
                                  x.idx         = lavsamplestats@x.idx[[g]],
                                  meanstructure = lavmodel@meanstructure,
                                  correlation   = correlation.flag)
                    } else {
                        A1[[g]] <- lav_mvnorm_h1_information_expected(
                                  sample.cov.inv = lavsamplestats@icov[[g]],
                                  #wt = WT, not needed
                                  x.idx          = lavsamplestats@x.idx[[g]],
                                  meanstructure  = lavmodel@meanstructure,
                                  correlation    = correlation.flag)
                    }
                } # conditional.x
            } # missing

            # stochastic group weight
            if(lavmodel@group.w.free) {
                # unweight!! (as otherwise, we would 'weight' again)
                a <- exp(lavimplied$group.w[[g]]) / lavsamplestats@nobs[[g]]
                A1[[g]] <- lav_matrix_bdiag( matrix(a, 1L, 1L), A1[[g]])
            }

        } # g
    } # ML

    # 3b. ML + multilevel
    else if(lavmodel@estimator == "ML" && lavdata@nlevels > 1L) {

        A1 <- vector("list", length = lavsamplestats@ngroups)

        # structured? compute model-implied statistics
        if(structured && length(lavimplied) == 0L) {
            lavimplied <- lav_model_implied(lavmodel)
        }

        # structured? lavimplied vs lavh1
        if(structured) {
            implied <- lavimplied
        } else {
            implied <- lavh1$implied
        }

        for(g in 1:lavsamplestats@ngroups) {

            MU.W    <- implied$mean[[ (g-1)*lavdata@nlevels + 1L ]]
            MU.B    <- implied$mean[[ (g-1)*lavdata@nlevels + 2L ]]
            SIGMA.W <- implied$cov[[  (g-1)*lavdata@nlevels + 1L ]]
            SIGMA.B <- implied$cov[[  (g-1)*lavdata@nlevels + 2L ]]

            # clustered data
            A1[[g]] <- lav_mvnorm_cluster_information_expected(
                           Lp           = lavdata@Lp[[g]],
                           Mu.W         = MU.W,
                           Sigma.W      = SIGMA.W,
                           Mu.B         = MU.B,
                           Sigma.B      = SIGMA.B,
                           x.idx        = lavsamplestats@x.idx[[g]])
        } # g
    } # ML + multilevel


    A1
}

lav_model_h1_information_observed <- function(lavobject      = NULL,
                                              lavmodel       = NULL,
                                              lavsamplestats = NULL,
                                              lavdata        = NULL,
                                              lavimplied     = NULL,
                                              lavh1          = NULL,
                                              lavcache       = NULL,
                                              lavoptions     = NULL) {

    if(!is.null(lavobject) && inherits(lavobject, "lavaan")) {
        lavmodel       <- lavobject@Model
        lavsamplestats <- lavobject@SampleStats
        lavdata        <- lavobject@Data
        lavimplied     <- lavobject@implied
        if(.hasSlot(lavobject, "h1")) {
            lavh1      <- lavobject@h1
        } else {
            lavh1      <- lav_h1_implied_logl(lavdata = lavobject@Data,
                                      lavsamplestats = lavobject@SampleStats,
                                      lavoptions = lavobject@Options)
        }
        lavcache       <- lavobject@Cache
        lavoptions     <- lavobject@Options
    }

    # sanity check
    if(length(lavh1) == 0L) {
        lavh1 <- lav_h1_implied_logl(lavdata = lavdata,
                                     lavsamplestats = lavsamplestats,
                                     lavoptions = lavoptions)
    }
    if(length(lavimplied) == 0L) {
        lavimplied <- lav_model_implied(lavmodel = lavmodel)
    }

    estimator <- lavmodel@estimator

    # structured?
    if(!is.null(lavoptions) &&
       !is.null(lavoptions$h1.information[1]) &&
       lavoptions$h1.information[1] == "unstructured") {
        structured <- FALSE
    } else {
        structured <- TRUE
    }

    # 1. WLS.V (=A1) for GLS/WLS
    if(lavmodel@estimator == "GLS"  || lavmodel@estimator == "WLS" ||
       lavmodel@estimator == "DLS") {
        A1 <- lavsamplestats@WLS.V
    }

    # 1b.
    else if(lavmodel@estimator == "DLS") {
        if(lavmodel@estimator.args$dls.GammaNT == "sample") {
            A1 <- lavsamplestats@WLS.V
        } else {
            A1 <- vector("list", length = lavsamplestats@ngroups)
            for(g in seq_len(lavsamplestats@ngroups)) {
                dls.a <- lavmodel@estimator.args$dls.a
                GammaNT <- lav_samplestats_Gamma_NT(
                               COV            = lavimplied$cov[[g]],
                               MEAN           = lavimplied$mean[[g]],
                               rescale        = FALSE,
                               x.idx          = lavsamplestats@x.idx[[g]],
                               fixed.x        = lavmodel@fixed.x,
                               conditional.x  = lavmodel@conditional.x,
                               meanstructure  = lavmodel@meanstructure,
                               slopestructure = lavmodel@conditional.x)
                W.DLS <- (1 - dls.a)*lavsamplestats@NACOV[[g]] + dls.a*GammaNT
                A1[[g]] <- lav_matrix_symmetric_inverse(W.DLS)
            }
        }
    }

    # 2. DWLS/ULS diagonal @WLS.VD slot
    else if(lavmodel@estimator == "DWLS"  || lavmodel@estimator == "ULS") {
        # diagonal only!!
        A1 <- lavsamplestats@WLS.VD
    }

    # 3a. ML single level
    else if(lavmodel@estimator == "ML" && lavdata@nlevels == 1L) {
        A1 <- vector("list", length=lavsamplestats@ngroups)

        # structured? compute model-implied statistics
        if(structured && length(lavimplied) == 0L) {
            lavimplied <- lav_model_implied(lavmodel)
        }

        for(g in 1:lavsamplestats@ngroups) {

            if(lavsamplestats@missing.flag) {
                # mvnorm
                # FIXME: allow for meanstructure = FALSE
                # FIXME: allow for conditional.x = TRUE
                if(lavmodel@meanstructure && structured) {
                    MEAN <- lavimplied$mean[[g]]
                } else {
                    MEAN <- lavsamplestats@missing.h1[[g]]$mu
                }

                if(structured) {
                    A1[[g]] <-
                      lav_mvnorm_missing_information_observed_samplestats(
                          Yp = lavsamplestats@missing[[g]],
                          # wt not needed
                          Mu = MEAN,
                          # meanstructure = lavmodel@meanstructure,
                          Sigma = lavimplied$cov[[g]],
                          x.idx = lavsamplestats@x.idx[[g]])
                } else {
                    A1[[g]] <-
                      lav_mvnorm_missing_information_observed_samplestats(
                          Yp = lavsamplestats@missing[[g]],
                          # wt not needed
                          Mu = MEAN,
                          # meanstructure = lavmodel@meanstructure,
                          Sigma = lavsamplestats@missing.h1[[g]]$sigma,
                          x.idx = lavsamplestats@x.idx[[g]])
                }
            } else {
                if(lavmodel@conditional.x) {
                    # mvreg
                    if(lavmodel@meanstructure && structured) {
                        RES.INT    <- lavimplied$res.int[[g]]
                        RES.SLOPES <- lavimplied$res.slopes[[g]]
                    } else {
                        RES.INT    <- lavsamplestats@res.int[[g]]
                        RES.SLOPES <- lavsamplestats@res.slopes[[g]]
                    }

                    if(structured) {
                        A1[[g]] <- lav_mvreg_information_observed_samplestats(
                            sample.res.int    = lavsamplestats@res.int[[g]],
                            sample.res.slopes = lavsamplestats@res.slopes[[g]],
                            sample.res.cov    = lavsamplestats@res.cov[[g]],
                            sample.mean.x     = lavsamplestats@mean.x[[g]],
                            sample.cov.x      = lavsamplestats@cov.x[[g]],
                            res.int           = RES.INT,
                            res.slopes        = RES.SLOPES,
                            #wt               = WT,
                            #meanstructure    = lavmodel@meanstructure,
                            res.cov           = lavimplied$res.cov[[g]])
                    } else {
                        A1[[g]] <- lav_mvreg_information_observed_samplestats(
                            sample.res.int    = lavsamplestats@res.int[[g]],
                            sample.res.slopes = lavsamplestats@res.slopes[[g]],
                            sample.res.cov    = lavsamplestats@res.cov[[g]],
                            sample.mean.x     = lavsamplestats@mean.x[[g]],
                            sample.cov.x      = lavsamplestats@cov.x[[g]],
                            res.int           = lavsamplestats@res.int[[g]],
                            res.slopes        = lavsamplestats@res.slopes[[g]],
                            #wt               = WT,
                            #meanstructure    = lavmodel@meanstructure,
                            res.cov           = lavsamplestats@res.cov[[g]])
                    }

                } else {
                    # conditional.x = FALSE
                    # mvnorm
                    if(lavmodel@meanstructure && structured) {
                        MEAN <- lavimplied$mean[[g]]
                    } else {
                        MEAN <- lavsamplestats@mean[[g]]
                    }

                    if(structured) {
                        A1[[g]] <- lav_mvnorm_information_observed_samplestats(
                                  sample.mean   = lavsamplestats@mean[[g]],
                                  sample.cov    = lavsamplestats@cov[[g]],
                                  Mu            = MEAN,
                                  Sigma         = lavimplied$cov[[g]],
                                  #wt           = WT, # not needed
                                  x.idx         = lavsamplestats@x.idx[[g]],
                                  meanstructure = lavmodel@meanstructure)
                    } else {
                        A1[[g]] <- lav_mvnorm_h1_information_observed_samplestats(
                                  sample.mean    = lavsamplestats@mean[[g]],
                                  sample.cov     = lavsamplestats@cov[[g]],
                                  sample.cov.inv = lavsamplestats@icov[[g]],
                                  #wt            = WT, not needed
                                  x.idx          = lavsamplestats@x.idx[[g]],
                                  meanstructure  = lavmodel@meanstructure)
                    }
                } # conditional.x
            } # missing

            # stochastic group weight
            if(lavmodel@group.w.free) {
                # unweight!!
                a <- exp(lavimplied$group.w[[g]]) / lavsamplestats@nobs[[g]]
                A1[[g]] <- lav_matrix_bdiag( matrix(a,1,1), A1[[g]])
            }

        } # g
    } # ML

    # 3b. ML + multilevel
    else if(lavmodel@estimator == "ML" && lavdata@nlevels > 1L) {

        A1 <- vector("list", length = lavsamplestats@ngroups)

        # structured? compute model-implied statistics
        if(structured && length(lavimplied) == 0L) {
            lavimplied <- lav_model_implied(lavmodel)
        }

        # structured? lavimplied vs lavh1
        if(structured) {
            implied <- lavimplied
        } else {
            implied <- lavh1$implied
        }

        for(g in 1:lavsamplestats@ngroups) {

            MU.W    <- implied$mean[[ (g-1)*lavdata@nlevels + 1L ]]
            MU.B    <- implied$mean[[ (g-1)*lavdata@nlevels + 2L ]]
            SIGMA.W <- implied$cov[[  (g-1)*lavdata@nlevels + 1L ]]
            SIGMA.B <- implied$cov[[  (g-1)*lavdata@nlevels + 2L ]]

            # clustered data
            A1[[g]] <- lav_mvnorm_cluster_information_observed(
                           Lp           = lavdata@Lp[[g]],
                           YLp          = lavsamplestats@YLp[[g]],
                           Mu.W         = MU.W,
                           Sigma.W      = SIGMA.W,
                           Mu.B         = MU.B,
                           Sigma.B      = SIGMA.B,
                           x.idx        = lavsamplestats@x.idx[[g]])
        } # g
    } # ML + multilevel

    A1
}

# outer product of the case-wise scores (gradients)
lav_model_h1_information_firstorder <- function(lavobject      = NULL,
                                                lavmodel       = NULL,
                                                lavsamplestats = NULL,
                                                lavdata        = NULL,
                                                lavimplied     = NULL,
                                                lavh1          = NULL,
                                                lavcache       = NULL,
                                                lavoptions     = NULL) {

    if(!is.null(lavobject) && inherits(lavobject, "lavaan")) {
        lavmodel       <- lavobject@Model
        lavsamplestats <- lavobject@SampleStats
        lavdata        <- lavobject@Data
        lavimplied     <- lavobject@implied
        if(.hasSlot(lavobject, "h1")) {
            lavh1      <- lavobject@h1
        } else {
            lavh1      <- lav_h1_implied_logl(lavdata = lavobject@Data,
                                      lavsamplestats = lavobject@SampleStats,
                                      lavoptions = lavobject@Options)
        }
        lavcache       <- lavobject@Cache
        lavoptions     <- lavobject@Options
    }

    # sanity check
    if(length(lavh1) == 0L) {
        lavh1 <- lav_h1_implied_logl(lavdata = lavdata,
                                     lavsamplestats = lavsamplestats,
                                     lavoptions = lavoptions)
    }
    if(length(lavimplied) == 0L) {
        lavimplied <- lav_model_implied(lavmodel = lavmodel)
    }

    estimator <- lavmodel@estimator
    if(!estimator %in% c("ML", "PML")) {
        stop("lavaan ERROR: information = \"first.order\" not available for estimator ", sQuote(estimator))
    }

    # structured?
    if(!is.null(lavoptions) &&
       !is.null(lavoptions$h1.information[1]) &&
       lavoptions$h1.information[1] == "unstructured") {
        structured <- FALSE
    } else {
        structured <- TRUE
    }

    # clustered?
    if(!is.null(lavoptions) &&
       !is.null(lavoptions$.clustered) &&
       lavoptions$.clustered) {
        clustered <- TRUE
        if(is.null(lavdata@Lp[[1]])) {
            stop("lavaan ERROR: lavdata@Lp is empty, while clustered = TRUE")
        }
        if(estimator == "PML") {
            stop("lavaan ERROR: clustered information is not (yet) available when estimator = \"PML\"")
        }
        #if(lavsamplestats@missing.flag) {
        #    stop("lavaan ERROR: clustered information is not (yet) available when missing = \"ML\"")
        #}
        #if(lavmodel@conditional.x) {
        #    stop("lavaan ERROR: clustered information is not (yet) available when conditional.x = TRUE")
        #}
        #if(!structured) {
        #    stop("lavaan ERROR: clustered information is not (yet) available when h1.information = \"unstructured\"")
        #}
    } else {
        clustered <- FALSE
    }

    # structured? compute model-implied statistics
    if(estimator == "PML" || structured) {
        if(length(lavimplied) == 0L) {
            lavimplied <- lav_model_implied(lavmodel)
        }
    }

    # structured? lavimplied vs lavh1
    if(structured) {
        implied <- lavimplied
    } else {
        implied <- lavh1$implied
    }

    B1 <- vector("list", length=lavsamplestats@ngroups)
    for(g in 1:lavdata@ngroups) {

        if(.hasSlot(lavdata, "weights")) {
            WT <- lavdata@weights[[g]]
        } else {
            WT <- NULL
        }

        if(estimator == "PML") {
            # slow approach: compute outer product of case-wise scores

            if(lavmodel@conditional.x) {
                SIGMA <- implied$res.cov[[g]]
                MU    <- implied$res.mean[[g]]
                TH    <- implied$res.th[[g]]
                PI    <- implied$res.slopes[[g]]
                EXO   <- lavdata@eXo[[g]]
            } else {
                SIGMA <- implied$cov[[g]]
                MU    <- implied$mean[[g]]
                TH    <- implied$th[[g]]
                PI    <- NULL
                EXO   <- NULL
            }
            SC <- pml_deriv1(Sigma.hat  = SIGMA,
                             Mu.hat     = MU,
                             TH         = TH,
                             th.idx     = lavmodel@th.idx[[g]],
                             num.idx    = lavmodel@num.idx[[g]],
                             X          = lavdata@X[[g]],
                             eXo        = EXO,
                             PI         = PI,
                             lavcache   = lavcache[[g]],
                             missing    = lavdata@missing,
                             scores     = TRUE,
                             negative   = FALSE)
            # information H1
            B1[[g]] <- lav_matrix_crossprod(SC)

        } else if(estimator == "ML" && lavdata@nlevels > 1L) {

            # if not-structured, we use lavh1, and that is always
            # 'unconditional' (for now)
            if(lavmodel@conditional.x && structured) {
                Res.Sigma.W <- implied$res.cov[[    (g-1)*lavdata@nlevels + 1L]]
                Res.Int.W   <- implied$res.int[[    (g-1)*lavdata@nlevels + 1L]]
                Res.Pi.W    <- implied$res.slopes[[ (g-1)*lavdata@nlevels + 1L]]

                Res.Sigma.B <- implied$res.cov[[    (g-1)*lavdata@nlevels + 2L]]
                Res.Int.B   <- implied$res.int[[    (g-1)*lavdata@nlevels + 2L]]
                Res.Pi.B    <- implied$res.slopes[[ (g-1)*lavdata@nlevels + 2L]]
                B1[[g]] <- lav_mvreg_cluster_information_firstorder(
                               Y1            = lavdata@X[[g]],
                               YLp           = lavsamplestats@YLp[[g]],
                               Lp            = lavdata@Lp[[g]],
                               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,
                               divide.by.two = TRUE)
            } else {
                MU.W    <- implied$mean[[ (g-1)*lavdata@nlevels + 1L ]]
                MU.B    <- implied$mean[[ (g-1)*lavdata@nlevels + 2L ]]
                SIGMA.W <- implied$cov[[  (g-1)*lavdata@nlevels + 1L ]]
                SIGMA.B <- implied$cov[[  (g-1)*lavdata@nlevels + 2L ]]
                B1[[g]] <- lav_mvnorm_cluster_information_firstorder(
                               Y1            = lavdata@X[[g]],
                               YLp           = lavsamplestats@YLp[[g]],
                               Lp            = lavdata@Lp[[g]],
                               Mu.W          = MU.W,
                               Sigma.W       = SIGMA.W,
                               Mu.B          = MU.B,
                               Sigma.B       = SIGMA.B,
                               x.idx         = lavsamplestats@x.idx[[g]],
                               divide.by.two = TRUE)
            }

        } else if(estimator == "ML" && lavdata@nlevels == 1L) {

            if(length(lavdata@cluster) > 0L) {
                cluster.idx <- lavdata@Lp[[g]]$cluster.idx[[2]]
            } else {
                cluster.idx <- NULL
            }

            if(lavsamplestats@missing.flag) {
                # mvnorm
                # FIXME: allow for meanstructure = FALSE
                # FIXME: allow for conditional.x = TRUE
                if(lavmodel@meanstructure && structured) {
                    MEAN <- lavimplied$mean[[g]]
                } else {
                    MEAN <- lavsamplestats@missing.h1[[g]]$mu
                }

                B1[[g]] <- lav_mvnorm_missing_information_firstorder(
                              Y = lavdata@X[[g]],
                              Mp = lavdata@Mp[[g]], wt = WT,
                              cluster.idx = cluster.idx,
                              Mu = MEAN,
                              # meanstructure = lavmodel@meanstructure,
                              Sigma = implied$cov[[g]],
                              x.idx = lavsamplestats@x.idx[[g]])

            } else {
                if(lavmodel@conditional.x) {
                    # mvreg
                    if(lavmodel@meanstructure && structured) {
                        RES.INT    <- lavimplied$res.int[[g]]
                        RES.SLOPES <- lavimplied$res.slopes[[g]]
                    } else {
                        RES.INT    <- lavsamplestats@res.int[[g]]
                        RES.SLOPES <- lavsamplestats@res.slopes[[g]]
                    }

                    B1[[g]] <- lav_mvreg_information_firstorder(
                                  Y              = lavdata@X[[g]],
                                  eXo            = lavdata@eXo[[g]],
                                  res.int        = RES.INT,
                                  res.slopes     = RES.SLOPES,
                                  #wt            = WT,
                                  #meanstructure = lavmodel@meanstructure,
                                  res.cov        = implied$res.cov[[g]])
                } else {
                    # conditional.x = FALSE
                    # mvnorm
                    if(lavmodel@meanstructure && structured) {
                        MEAN <- lavimplied$mean[[g]]
                    } else {
                        # NOTE: the information matrix will be the same (minus
                        # the meanstructure block), but once INVERTED, the
                        # standard errors will be (slightly) smaller!!!
                        # This is only visibile when estimator = "MLF"
                        # (or information = "first.order")
                        MEAN <- lavsamplestats@mean[[g]] # saturated
                    }

                    if(structured) {
                        B1[[g]] <- lav_mvnorm_information_firstorder(
                                  Y = lavdata@X[[g]],
                                  Mu = MEAN, Sigma = lavimplied$cov[[g]],
                                  wt = WT,
                                  cluster.idx = cluster.idx,
                                  x.idx = lavsamplestats@x.idx[[g]],
                                  meanstructure = lavmodel@meanstructure)
                    } else {
                        B1[[g]] <- lav_mvnorm_h1_information_firstorder(
                                  Y = lavdata@X[[g]],
                                  sample.cov.inv = lavsamplestats@icov[[g]],
                                  Gamma = lavsamplestats@NACOV[[g]],
                                  wt = WT,
                                  cluster.idx = cluster.idx, # only if wt
                                  x.idx = lavsamplestats@x.idx[[g]],
                                  meanstructure = lavmodel@meanstructure)
                    }
                } # mvnorm
            } # missing
        } # ML

        # stochastic group weight
        if(lavmodel@group.w.free) {
            # unweight!!
            a <- exp(lavimplied$group.w[[g]]) / lavsamplestats@nobs[[g]]
            B1[[g]] <- lav_matrix_bdiag( matrix(a,1,1), B1[[g]])
        }

    } # g

    B1
}


# asymptotic variance matrix (=Gamma/N) of the unrestricted (H1)
# sample statistics
#
# FIXME: make this work for categorical/GLS/WLS/...
#
lav_model_h1_acov <- function(lavobject      = NULL,
                              lavmodel       = NULL,
                              lavsamplestats = NULL,
                              lavdata        = NULL,
                              lavoptions     = NULL,
                              lavimplied     = NULL,
                              lavh1          = NULL,
                              lavcache       = NULL,
                              meanstructure  = NULL,  # if specified, use it
                              h1.information = NULL,  # if specified, use it
                              se = NULL) {            # if specified, use it

    if(!is.null(lavobject) && inherits(lavobject, "lavaan")) {
        lavmodel       <- lavobject@Model
        lavsamplestats <- lavobject@SampleStats
        lavdata        <- lavobject@Data
        lavimplied     <- lavobject@implied
        if(.hasSlot(lavobject, "h1")) {
            lavh1      <- lavobject@h1
        } else {
            lavh1      <- lav_h1_implied_logl(lavdata = lavobject@Data,
                                      lavsamplestats = lavobject@SampleStats,
                                      lavoptions = lavobject@Options)
        }
        lavcache       <- lavobject@Cache
        lavoptions     <- lavobject@Options
    }

    # sanity check
    if(length(lavh1) == 0L) {
        lavh1 <- lav_h1_implied_logl(lavdata = lavdata,
                                     lavsamplestats = lavsamplestats,
                                     lavoptions = lavoptions)
    }
    if(length(lavimplied) == 0L) {
        lavimplied <- lav_model_implied(lavmodel = lavmodel)
    }

    # override
    if(!is.null(meanstructure)) {
        lavoptions$meanstructure <- meanstructure
    }
    if(!is.null(h1.information)) {
        lavoptions$h1.information[1] <- h1.information
    }
    if(!is.null(se)) {
        lavoptions$se <- se
    }



    # information
    information <- lavoptions$information[1] # ALWAYS used the first

    # compute information matrix
    if(information == "observed") {
        I1 <- lav_model_h1_information_observed(lavmodel = lavmodel,
            lavsamplestats = lavsamplestats, lavdata = lavdata,
            lavimplied = lavimplied, lavh1 = lavh1,
            lavcache = lavcache, lavoptions = lavoptions)
    } else if(information == "expected") {
        I1 <- lav_model_h1_information_expected(lavmodel = lavmodel,
            lavsamplestats = lavsamplestats, lavdata = lavdata,
            lavimplied = lavimplied, lavh1 = lavh1,
            lavcache = lavcache, lavoptions = lavoptions)
    } else if(information == "first.order") {
        I1 <- lav_model_h1_information_firstorder(lavmodel = lavmodel,
            lavsamplestats = lavsamplestats, lavdata = lavdata,
            lavimplied = lavimplied, lavh1 = lavh1,
            lavcache = lavcache, lavoptions = lavoptions)
    }

    if(lavoptions$se %in% c("robust.huber.white", "robust.sem")) {
        J1 <- lav_model_h1_information_firstorder(lavmodel = lavmodel,
            lavsamplestats = lavsamplestats, lavdata = lavdata,
            lavimplied = lavimplied, lavh1 = lavh1,
            lavcache = lavcache, lavoptions = lavoptions)
    }

    # compute ACOV per group
    ACOV <- vector("list", length = lavdata@ngroups)
    for(g in 1:lavdata@ngroups) {

        # denominator
        if(lavdata@nlevels == 1L) {
            Ng <- lavsamplestats@nobs[[g]]
        } else {
            Ng <- lavdata@Lp[[g]]$nclusters[[2]]
        }

        # invert information
        I1.g.inv <- try(lav_matrix_symmetric_inverse(I1[[g]]), silent = TRUE)
        if(inherits(I1.g.inv, "try-error")) {
            stop("lavaan ERROR: could not invert h1 information matrix in group ", g)
        }

        # which type of se?
        if(lavoptions$se %in% c("standard", "none")) {
            ACOV[[g]] <- 1/Ng * I1.g.inv
        } else if(lavoptions$se %in% c("robust.huber.white", "robust.sem")) {
            ACOV[[g]] <- 1/Ng * (I1.g.inv %*% J1[[g]] %*% I1.g.inv)
        }
    }

    ACOV
}