R/lav_model_lik.R

In nietsnel/psindex: Parameter Stability Index

# casewise likelihoods

# closed-form marginal likelihood
# - classic SEM models, continous observed variables only
lav_model_lik_ml <- function(lavmodel       = NULL,
                             GLIST          = NULL,
                             lavdata        = NULL,
                             lavsamplestats = NULL) {


}


# marginal ML
lav_model_lik_mml <- function(lavmodel      = NULL,
                              THETA         = NULL,
                              TH            = NULL,
                              GLIST         = NULL,
                              group         = 1L,
                              lavdata       = NULL,
                              sample.mean   = NULL,
                              sample.mean.x = NULL,
                              lavcache      = NULL) {

    conditional.x <- lavmodel@conditional.x

    # data for this group
    X <- lavdata@X[[group]]; nobs <- nrow(X); nvar <- ncol(X)
    eXo <- lavdata@eXo[[group]]

    # MLIST (for veta and yhat)
    mm.in.group <- 1:lavmodel@nmat[group] + cumsum(c(0,lavmodel@nmat))[group]
    MLIST <- GLIST[ mm.in.group ]

    # quadrature points
    GH <- lavcache[[group]]$GH; nGH <- nrow(GH$x)
    nfac <- ncol(GH$x)

    # compute VETAx (latent lv only)
    lv.dummy.idx <- c(lavmodel@ov.y.dummy.lv.idx[[group]],
                      lavmodel@ov.x.dummy.lv.idx[[group]])
    VETAx <- computeVETAx.LISREL(MLIST = MLIST,
                                 lv.dummy.idx = lv.dummy.idx)
    #VETAx <- computeVETAx.LISREL(MLIST = MLIST)
    # check for negative values?
    if(any(diag(VETAx) < 0)) {
        warning("psindex WARNING: --- VETAx contains negative values")
        print(VETAx)
        return(0)
    }

    # cholesky?
    #if(is.null(lavmodel@control$cholesky)) {
        CHOLESKY <- TRUE    
    #} else {
    #    CHOLESKY <- as.logical(lavmodel@control$cholesky)
        #if(nfac > 1L && !CHOLESKY) {
        #    warning("psindex WARNING: CHOLESKY is OFF but nfac > 1L")
        #}
    #}

    if(!CHOLESKY) {
        # we should still 'scale' the factors, if std.lv=FALSE
        ETA.sd <- sqrt( diag(VETAx) )
    } else {
        # cholesky takes care of scaling
        tchol.VETA <- try(chol(VETAx), silent = TRUE)
        if(inherits(tchol.VETA, "try-error")) {
            warning("psindex WARNING: --- VETAx not positive definite")
            print(VETAx)
            return(0)
        }
        if(!is.null(MLIST$alpha) || !is.null(MLIST$gamma)) {
            if(conditional.x) {
                EETAx <- computeEETAx.LISREL(MLIST = MLIST, eXo = eXo, N = nobs,
                        sample.mean = sample.mean,
                        ov.y.dummy.ov.idx = lavmodel@ov.y.dummy.ov.idx[[group]],
                        ov.x.dummy.ov.idx = lavmodel@ov.x.dummy.ov.idx[[group]],
                        ov.y.dummy.lv.idx = lavmodel@ov.y.dummy.lv.idx[[group]],
                        ov.x.dummy.lv.idx = lavmodel@ov.x.dummy.lv.idx[[group]])
            } else {
                EETA <- computeEETA.LISREL(MLIST = MLIST, 
                        mean.x = sample.mean.x,
                        sample.mean = sample.mean,
                        ov.y.dummy.ov.idx = lavmodel@ov.y.dummy.ov.idx[[group]],
                        ov.x.dummy.ov.idx = lavmodel@ov.x.dummy.ov.idx[[group]],
                        ov.y.dummy.lv.idx = lavmodel@ov.y.dummy.lv.idx[[group]],
                        ov.x.dummy.lv.idx = lavmodel@ov.x.dummy.lv.idx[[group]])
            }
            #if(length(lv.dummy.idx) > 0L) {
            #    EETAx <- EETAx[,-lv.dummy.idx,drop=FALSE]
            #}
        }
    }

    # compute (log)lik for each node, for each observation
    SUM.LOG.FY <- matrix(0, nrow=nGH, ncol=nobs)
    for(q in 1:nGH) {

        # current value(s) for ETA
        #eta <- matrix(0, nrow = 1, ncol = ncol(MLIST$lambda))
   
        # non-dummy elements -> quadrature points
        #eta[1L, -lv.dummy.idx] <- GH$x[q,,drop=FALSE]
        XQ <- GH$x[q,,drop=FALSE]
     
        # rescale/unwhiten
        if(CHOLESKY) {
            # un-orthogonalize
            XQ <- XQ %*% tchol.VETA
        } else {
            # no unit scale? (un-standardize)
            XQ <- sweep(XQ, MARGIN=2, STATS=ETA.sd, FUN="*")
        }

        eta <- matrix(0, nrow = 1, ncol = ncol(MLIST$lambda))
        if(length(lv.dummy.idx) > 0L) {
            eta[, -lv.dummy.idx] <- XQ
        } else {
            eta <- XQ
        }

        # eta_i = alpha + BETA eta_i + GAMMA eta_i + error
        #
        # - direct effect of BETA is already in VETAx, and hence tchol.VETA
        # - need to add alpha, and GAMMA eta_i
        if(!is.null(MLIST$alpha) || !is.null(MLIST$gamma)) {
            if(conditional.x) {
                eta <- sweep(EETAx, MARGIN=2, STATS=eta, FUN="+")
            } else {
                eta <- eta + EETA
            }
        }

        # compute yhat for this node (eta)
        if(lavmodel@conditional.x) {
            yhat <- computeEYetax.LISREL(MLIST = MLIST, eXo = eXo,
                    ETA = eta, sample.mean = sample.mean,
                    ov.y.dummy.ov.idx = lavmodel@ov.y.dummy.ov.idx[[group]],
                    ov.x.dummy.ov.idx = lavmodel@ov.x.dummy.ov.idx[[group]],
                    ov.y.dummy.lv.idx = lavmodel@ov.y.dummy.lv.idx[[group]],
                    ov.x.dummy.lv.idx = lavmodel@ov.x.dummy.lv.idx[[group]])
        } else {
            yhat <- computeEYetax3.LISREL(MLIST = MLIST,
                    ETA = eta, sample.mean = sample.mean,
                    mean.x = sample.mean.x,
                    ov.y.dummy.ov.idx = lavmodel@ov.y.dummy.ov.idx[[group]],
                    ov.x.dummy.ov.idx = lavmodel@ov.x.dummy.ov.idx[[group]],
                    ov.y.dummy.lv.idx = lavmodel@ov.y.dummy.lv.idx[[group]],
                    ov.x.dummy.lv.idx = lavmodel@ov.x.dummy.lv.idx[[group]])
        }
 
        # compute fy.var, for this node (eta): P(Y_i =  y_i | eta_i, x_i)
        log.fy.var <- lav_predict_fy_internal(X = X, yhat = yhat,
                          TH = TH, THETA = THETA,
                          num.idx = lavmodel@num.idx[[group]],
                          th.idx  = lavmodel@th.idx[[group]],
                          link    = lavmodel@link, log. = TRUE)

        # if log, fy is just the sum of log.fy.var
        log.fy <- apply(log.fy.var, 1L, sum)

        # store log likelihoods for this node
        SUM.LOG.FY[q,] <- log.fy
    }

    # integration
    lik <- as.numeric( t(GH$w) %*% exp(SUM.LOG.FY) )

    # avoid underflow
    idx <- which(lik < exp(-600))
    if(length(idx) > 0L) {
        lik[idx] <- exp(-600)
    }

    lik
}