R/lav_h1_logl.R

In lavaan: Latent Variable Analysis

# compute logl for the unrestricted (h1) model -- per group
lav_h1_logl <- function(lavdata        = NULL,
                        lavsamplestats = NULL,
                        lavoptions     = NULL) {

    # number of groups
    ngroups <- lavdata@ngroups

    logl.group <- rep(as.numeric(NA), ngroups)

    # should compute logl, or return NA?
    logl.ok <- FALSE
    if(lavoptions$estimator %in% c("ML", "MML")) {
        # check if everything is numeric, OR if we have exogenous
        # factor with 2 levels only
        if(all(lavdata@ov$type == "numeric")) {
            logl.ok <- TRUE
        } else {
            not.idx <- which(lavdata@ov$type != "numeric")
            for(i in not.idx) {
                if(lavdata@ov$type[i] == "factor" &&
                   lavdata@ov$exo[i] == 1L &&
                   lavdata@ov$nlev[i] == 2L) {
                    logl.ok <- TRUE
                } else {
                    logl.ok <- FALSE
                    break
                }
            }
        }
    }

    # new in 0.6-9 (so SAM can handle N<P)
    if(!is.null(lavoptions$sample.icov) && !lavoptions$sample.icov) {
        logl.ok <- FALSE
    }

    if(logl.ok) {
        for(g in seq_len(ngroups) ) {
            if(lavdata@nlevels > 1L) {
                OUT <- lav_mvnorm_cluster_em_sat(YLp  = lavsamplestats@YLp[[g]],
                                             Lp       = lavdata@Lp[[g]],
                                             verbose  = FALSE,
                                             tol      = 1e-04,     # option?
                                             min.variance = 1e-05, # option?
                                             max.iter = 5000L)     # option?
                # store logl per group
                logl.group[g] <- OUT$logl
            } else if(lavsamplestats@missing.flag) {
                logl.group[g] <-
                    lav_mvnorm_missing_loglik_samplestats(
                        Yp     = lavsamplestats@missing[[g]],
                        Mu     = lavsamplestats@missing.h1[[g]]$mu,
                        Sigma  = lavsamplestats@missing.h1[[g]]$sigma,
                        x.idx  = lavsamplestats@x.idx[[g]],
                        x.mean = lavsamplestats@mean.x[[g]],
                        x.cov  = lavsamplestats@cov.x[[g]])
            } else { # single-level, complete data
                # all we need is: logdet of covariance matrix, nobs and nvar
                if(lavoptions$conditional.x) {
                    logl.group[g] <-
                        lav_mvnorm_h1_loglik_samplestats(
                            sample.cov.logdet =
                                lavsamplestats@res.cov.log.det[[g]],
                            sample.nvar       =
                                NCOL(lavsamplestats@res.cov[[g]]),
                            sample.nobs       = lavsamplestats@nobs[[g]])
                } else {
                    logl.group[g] <-
                        lav_mvnorm_h1_loglik_samplestats(
                            sample.cov.logdet = lavsamplestats@cov.log.det[[g]],
                            sample.nvar       = NCOL(lavsamplestats@cov[[g]]),
                            sample.nobs       = lavsamplestats@nobs[[g]],
                            x.idx             = lavsamplestats@x.idx[[g]],
                            x.cov             = lavsamplestats@cov.x[[g]])
                }
            } # complete
        } # g
    } # logl.ok is TRUE

    out <- list(loglik       = sum(logl.group),
                loglik.group = logl.group)

    out
}