psindex: Parameter Stability Index

# here, we compute various versions of the `information' matrix
# NOTE:
# 1) we ALWAYS compute the UNIT information (not the total information)
# 
# 2) by default, we ignore the constraints (we deal with this when we
#    take the inverse later on)

lav_model_information <- function(lavmodel       = NULL,
                                  lavsamplestats = NULL,
                                  lavdata        = NULL,
                                  lavimplied     = NULL,
                                  lavh1          = NULL,
                                  Delta          = NULL,
                                  lavcache       = NULL,
                                  lavoptions     = NULL,
                                  extra          = FALSE,
                                  augmented      = FALSE,
                                  inverted       = FALSE,
                                  use.ginv       = FALSE) {

    estimator   <- lavmodel@estimator
    information <- lavoptions$information

    # compute information matrix
    if(information == "observed") {
        if(lavsamplestats@missing.flag || lavdata@nlevels > 1L) {
            group.weight <- FALSE
        } else {
            group.weight <- TRUE
        }
        E <- lav_model_information_observed(lavmodel = lavmodel,
            lavsamplestats = lavsamplestats, lavdata = lavdata,
            lavimplied = lavimplied, lavh1 = lavh1,
            lavcache = lavcache, group.weight = group.weight,
            lavoptions = lavoptions, extra = extra,
            augmented = augmented, inverted = inverted, use.ginv = use.ginv)
    } else if(information == "expected") {
        E <- lav_model_information_expected(lavmodel = lavmodel,
            lavsamplestats = lavsamplestats, lavdata = lavdata,
            lavimplied = lavimplied, lavh1 = lavh1,
            lavcache = lavcache, lavoptions = lavoptions, extra = extra, 
            augmented = augmented, inverted = inverted, use.ginv = use.ginv)
    } else if(information == "first.order") {
        E <- lav_model_information_firstorder(lavmodel = lavmodel,
            lavsamplestats = lavsamplestats, lavdata = lavdata,
            lavimplied = lavimplied, lavh1 = lavh1,
            lavcache = lavcache, lavoptions = lavoptions, #extra = extra,
            check.pd = FALSE,
            augmented = augmented, inverted = inverted, use.ginv = use.ginv)
    }

    # information, augmented information, or inverted information
    E
}

# fisher/expected information
#
# information = Delta' I1 Delta, where I1 is the unit information of
# the saturated model (evaluated either at the structured or unstructured
# estimates)
lav_model_information_expected <- function(lavmodel       = NULL,
                                           lavsamplestats = NULL,
                                           lavdata        = NULL,
                                           lavoptions     = NULL,
                                           lavimplied     = NULL,
                                           lavh1          = NULL,
                                           Delta          = NULL,
                                           lavcache       = NULL,
                                           extra          = FALSE,
                                           augmented      = FALSE,
                                           inverted       = FALSE,
                                           use.ginv       = FALSE) {

    if(inverted) {
        augmented <- TRUE
    }

    # 1. Delta
    if(is.null(Delta)) {
        Delta <- computeDelta(lavmodel = lavmodel)
    }


    # 2. H1 information (single level)
    if(lavdata@nlevels == 1L) {
        A1 <- lav_model_h1_information_expected(lavmodel       = lavmodel,
                                                lavsamplestats = lavsamplestats,
                                                lavdata        = lavdata,
                                                lavoptions     = lavoptions,
                                                lavimplied     = lavimplied,
                                                lavh1          = lavh1,
                                                lavcache       = lavcache)
    }

    # 3. compute Information per group
    Info.group  <- vector("list", length=lavsamplestats@ngroups)
    for(g in 1:lavsamplestats@ngroups) {
        # note LISREL documentation suggests (Ng - 1) instead of Ng...
        fg <- lavsamplestats@nobs[[g]]/lavsamplestats@ntotal

        # multilevel
        if(lavdata@nlevels > 1L) {
            # here, we assume only 2 levels, at [[1]] and [[2]]
            if(lavoptions$h1.information == "structured") {
                Sigma.W <- lavimplied$cov[[  (g-1)*2 + 1]]
                Mu.W    <- lavimplied$mean[[ (g-1)*2 + 1]]
                Sigma.B <- lavimplied$cov[[  (g-1)*2 + 2]]
                Mu.B    <- lavimplied$mean[[ (g-1)*2 + 2]]
            } else {
                Sigma.W <- lavh1$implied$cov[[  (g-1)*2 + 1]]
                Mu.W    <- lavh1$implied$mean[[ (g-1)*2 + 1]]
                Sigma.B <- lavh1$implied$cov[[  (g-1)*2 + 2]]
                Mu.B    <- lavh1$implied$mean[[ (g-1)*2 + 2]]
            }
            Lp      <- lavdata@Lp[[g]]

            Info.g <- 
                lav_mvnorm_cluster_information_expected_delta(Lp  = Lp,
                                                      Delta       = Delta[[g]],
                                                      Mu.W        = Mu.W,
                                                      Sigma.W     = Sigma.W,
                                                      Mu.B        = Mu.B,
                                                      Sigma.B     = Sigma.B,
                                                      Sinv.method = "eigen")
            Info.group[[g]] <- fg * Info.g
        } else {
            # compute information for this group
            if(lavmodel@estimator %in% c("DWLS", "ULS")) {
                # diagonal weight matrix
                Delta2 <- sqrt(A1[[g]]) * Delta[[g]]
                Info.group[[g]] <- fg * crossprod(Delta2)
            } else {
                # full weight matrix
                Info.group[[g]] <- 
                    fg * ( crossprod(Delta[[g]], A1[[g]]) %*% Delta[[g]] )
            }
        }
    } # g

    # 4. assemble over groups
    Information <- Info.group[[1]]
    if(lavsamplestats@ngroups > 1) {
        for(g in 2:lavsamplestats@ngroups) {
            Information <- Information + Info.group[[g]]
        }
    }

    # 5. augmented information?
    if(augmented) {
        Information <- 
            lav_model_information_augment_invert(lavmodel    = lavmodel,
                                                 information = Information,
                                                 inverted    = inverted,
                                                 use.ginv    = use.ginv)
    }

    if(extra) {
        attr(Information, "Delta") <- Delta
        attr(Information, "WLS.V") <- A1 # unweighted
    }

    # possibly augmented/inverted
    Information
}

# only for Mplus MLM
lav_model_information_expected_MLM <- function(lavmodel       = NULL, 
                                               lavsamplestats = NULL, 
                                               Delta          = NULL,
                                               extra          = FALSE,
                                               augmented      = FALSE,
                                               inverted       = FALSE,
                                               use.ginv       = FALSE) {

    if(inverted) {
        augmented <- TRUE
    }

    if(is.null(Delta)) {
        Delta = computeDelta(lavmodel = lavmodel)
    }

    # compute A1
    A1 <- vector("list", length=lavsamplestats@ngroups)
    if(lavmodel@group.w.free) {
        GW <- unlist(computeGW(lavmodel = lavmodel))
    }
    for(g in 1:lavsamplestats@ngroups) {
        A1[[g]] <- lav_mvnorm_h1_information_expected(
                          sample.cov     = lavsamplestats@cov[[g]],
                          sample.cov.inv = lavsamplestats@icov[[g]],
                          x.idx          = lavsamplestats@x.idx[[g]])
        # the same as GLS... (except for the N/N-1 scaling)
        if(lavmodel@group.w.free) {
            # unweight!!
            a <- exp(GW[g]) / lavsamplestats@nobs[[g]]
            # a <- exp(GW[g]) * lavsamplestats@ntotal / lavsamplestats@nobs[[g]]
            A1[[g]] <- lav_matrix_bdiag( matrix(a,1,1), A1[[g]])
        }
    }

    # compute Information per group
    Info.group  <- vector("list", length=lavsamplestats@ngroups)
    for(g in 1:lavsamplestats@ngroups) {
        fg <- lavsamplestats@nobs[[g]]/lavsamplestats@ntotal
        # compute information for this group
        Info.group[[g]] <- fg * (t(Delta[[g]]) %*% A1[[g]] %*% Delta[[g]])
    }

    # assemble over groups
    Information <- Info.group[[1]]
    if(lavsamplestats@ngroups > 1) {
        for(g in 2:lavsamplestats@ngroups) {
            Information <- Information + Info.group[[g]]
        }
    }

    # augmented information?
    if(augmented) {
        Information <-
            lav_model_information_augment_invert(lavmodel    = lavmodel,
                                                 information = Information,
                                                 inverted    = inverted,
                                                 use.ginv    = use.ginv)
    }

    if(extra) {
        attr(Information, "Delta") <- Delta
        attr(Information, "WLS.V") <- A1 # unweighted
    }

    Information
}


lav_model_information_observed <- function(lavmodel       = NULL,
                                           lavsamplestats = NULL,
                                           lavdata        = NULL,
                                           lavimplied     = NULL,
                                           lavh1          = NULL,
                                           lavcache       = NULL,
                                           lavoptions     = NULL,
                                           extra          = FALSE,
                                           group.weight   = TRUE,
                                           augmented      = FALSE,
                                           inverted       = FALSE,
                                           use.ginv       = FALSE) {
    if(inverted) {
        augmented <- TRUE
    }

    # observed.information:
    #     - "hessian": second derivative of objective function
    #     - "h1": observed information matrix of saturated (h1) model, 
    #             pre- and post-multiplied by the jacobian of the model
    #             parameters (Delta), usually evaluated at the structured
    #             sample statistics (but this depends on the h1.information 
    #             option)
    if(!is.null(lavoptions) &&
       !is.null(lavoptions$observed.information) &&
       lavoptions$observed.information == "h1") {
        observed.information <- "h1"
    } else {
        observed.information <- "hessian"
    }

    # HESSIAN based
    if(observed.information == "hessian") {
        Hessian <- lav_model_hessian(lavmodel       = lavmodel,
                                     lavsamplestats = lavsamplestats,
                                     lavdata        = lavdata,
                                     lavoptions     = lavoptions,
                                     lavcache       = lavcache,
                                     group.weight   = group.weight)

        # NOTE! What is the relationship between the Hessian of the objective
        # function, and the `information' matrix (unit or total)
      
        # 1. in psindex, we ALWAYS minimize, so the Hessian is already pos def
        # 2. currently, all estimators give unit information, except MML and PML
        #    so, no need to divide by N
        Information <- Hessian

        # divide by 'N' for MML and PML
        if(lavmodel@estimator == "PML" || lavmodel@estimator == "MML") {
            Information <- Information / lavsamplestats@ntotal
        }

        # if multilevel, we should divide by 'J', the number of clusters
        if(lavdata@nlevels > 1L) {
            NC <- 0
            for(g in 1:lavsamplestats@ngroups) {
                NC <- NC + lavdata@Lp[[g]]$nclusters[[2]]
            }
            Information <- Information * lavsamplestats@ntotal / NC
        }
    }

    # using 'observed h1 information'
    # we need DELTA and 'WLS.V' (=A1)

    if(observed.information == "h1" || extra) {
        # 1. Delta
        Delta <- computeDelta(lavmodel = lavmodel)

        # 2. H1 information
    
        A1 <- lav_model_h1_information_observed(lavmodel       = lavmodel,
                                                lavsamplestats = lavsamplestats,
                                                lavdata        = lavdata,
                                                lavoptions     = lavoptions,
                                                lavimplied     = lavimplied,
                                                lavh1          = lavh1,
                                                lavcache       = lavcache)
    }

    if(observed.information == "h1") {
        # compute Information per group
        Info.group  <- vector("list", length=lavsamplestats@ngroups)
        for(g in 1:lavsamplestats@ngroups) {
            fg <- lavsamplestats@nobs[[g]]/lavsamplestats@ntotal
            # compute information for this group
            if(lavmodel@estimator %in% c("DWLS", "ULS")) {
                # diagonal weight matrix
                Delta2 <- sqrt(A1[[g]]) * Delta[[g]]
                Info.group[[g]] <- fg * crossprod(Delta2)
            } else {
                # full weight matrix
                Info.group[[g]] <-
                    fg * ( crossprod(Delta[[g]], A1[[g]]) %*% Delta[[g]] )
            }
        }

        # assemble over groups
        Information <- Info.group[[1]]
        if(lavsamplestats@ngroups > 1) {
            for(g in 2:lavsamplestats@ngroups) {
                Information <- Information + Info.group[[g]]
            }
        }
    }

    # augmented information?
    if(augmented) {
        Information <-
            lav_model_information_augment_invert(lavmodel    = lavmodel,
                                                 information = Information,
                                                 inverted    = inverted,
                                                 use.ginv    = use.ginv)
    }

    if(extra) {
        attr(Information, "Delta") <- Delta
        attr(Information, "WLS.V") <- A1
    }

    Information
}

# outer product of the case-wise scores (gradients)
lav_model_information_firstorder <- function(lavmodel       = NULL,
                                             lavsamplestats = NULL,
                                             lavdata        = NULL,
                                             lavimplied     = NULL,
                                             lavh1          = NULL,
                                             lavcache       = NULL,
                                             lavoptions     = NULL,
                                             check.pd       = FALSE,
                                             extra          = FALSE,
                                             augmented      = FALSE,
                                             inverted       = FALSE,
                                             use.ginv       = FALSE) {
    if(!lavmodel@estimator %in% c("ML", "PML")) {
        stop("psindex ERROR: information = \"first.order\" not available for estimator ", sQuote(lavmodel@estimator))
    }

    if(inverted) {
        augmented <- TRUE
    }

    B0.group <- vector("list", lavsamplestats@ngroups)

    # 1. Delta
    Delta <- computeDelta(lavmodel = lavmodel)

    # 2. H1 information
    B1 <- lav_model_h1_information_firstorder(lavmodel       = lavmodel,
                                              lavsamplestats = lavsamplestats,
                                              lavdata        = lavdata,
                                              lavoptions     = lavoptions,
                                              lavimplied     = lavimplied,
                                              lavh1          = lavh1,
                                              lavcache       = lavcache)

    # 3. compute Information per group
    Info.group  <- vector("list", length=lavsamplestats@ngroups)
    for(g in 1:lavsamplestats@ngroups) {

        # unweighted (needed in lav_test?)
        B0.group[[g]] <- t(Delta[[g]]) %*% B1[[g]] %*% Delta[[g]] 
       
        fg <- lavsamplestats@nobs[[g]]/lavsamplestats@ntotal
        # compute information for this group
        Info.group[[g]] <- fg * B0.group[[g]]
    }

    # 4. assemble over groups
    Information <- Info.group[[1]]
    if(lavsamplestats@ngroups > 1) {
        for(g in 2:lavsamplestats@ngroups) {
            Information <- Information + Info.group[[g]]
        }
    }


    # NOTE: for MML and PML, we get 'total' information (instead of unit)
    # divide by 'N' for MML and PML
    if(lavmodel@estimator == "PML" || lavmodel@estimator == "MML") {
        Information <- Information / lavsamplestats@ntotal
        for(g in 1:lavsamplestats@ngroups) {
            B0.group[[g]] <- B0.group[[g]] / lavsamplestats@ntotal
        }
    }

    # augmented information?
    if(augmented) {
        Information <-
            lav_model_information_augment_invert(lavmodel    = lavmodel,
                                                 information = Information,
                                                 check.pd    = check.pd,
                                                 inverted    = inverted,
                                                 use.ginv    = use.ginv)
    }

    if(extra) {
        attr(Information, "B0.group") <- B0.group
        attr(Information, "Delta") <- Delta
        attr(Information, "WLS.V") <- B1
    }

    Information
}


# create augmented information matrix (if needed), and take the inverse
# (if inverted = TRUE), returning only the [1:npar, 1:npar] elements
lav_model_information_augment_invert <- function(lavmodel    = NULL,
                                                 information = NULL,
                                                 inverted    = FALSE,
                                                 check.pd    = FALSE,
                                                 use.ginv    = FALSE) {

    npar <- nrow(information)
    is.augmented <- FALSE

    # handle constraints
    if(nrow(lavmodel@con.jac) > 0L) {
        H <- lavmodel@con.jac
        inactive.idx <- attr(H, "inactive.idx")
        lambda <- lavmodel@con.lambda # lagrangean coefs
        if(length(inactive.idx) > 0L) {
            H <- H[-inactive.idx,,drop=FALSE]
            lambda <- lambda[-inactive.idx]
        }
        if(nrow(H) > 0L) {
            is.augmented <- TRUE   
            H0 <- matrix(0,nrow(H),nrow(H))
            H10 <- matrix(0, ncol(information), nrow(H))
            DL <- 2*diag(lambda, nrow(H), nrow(H))
            # FIXME: better include inactive + slacks??
            E3 <- rbind( cbind(     information,  H10, t(H)),
                         cbind(          t(H10),   DL,  H0),
                         cbind(               H,   H0,  H0)  )
            information <- E3   
        }
    }

    if(check.pd) {
        eigvals <- eigen(information, symmetric = TRUE, 
                         only.values = TRUE)$values
        if(any(eigvals < -1 * .Machine$double.eps^(3/4))) {
            warning("psindex WARNING: matrix based on first order outer product of the derivatives is not positive definite; the model may not be identified")
        }
    }

    if(inverted) {
        if(is.augmented) {
            # note: default tol in MASS::ginv is sqrt(.Machine$double.eps)
            #       which seems a bit too conservative
            #       from 0.5-20, we changed this to .Machine$double.eps^(3/4)
            information <- 
                try( MASS::ginv(information, 
                                tol = .Machine$double.eps^(3/4))[1:npar, 
                                                                 1:npar, 
                                                                 drop = FALSE],
                     silent = TRUE )
        } else {
            if(use.ginv) {
                information <- try( MASS::ginv(information,
                                               tol = .Machine$double.eps^(3/4)),
                                    silent = TRUE )
            } else {
                information <- try( solve(information), silent = TRUE )
            }
        }
    }

    # augmented/inverted information
    information
}

lav_model_information_expected_2l <- function(lavmodel       = NULL,
                                              lavsamplestats = NULL,
                                              lavdata        = NULL,
                                              lavoptions     = NULL,
                                              lavimplied     = NULL,
                                              lavh1          = NULL,
                                              g              = 1L) {
    # see Yuan & Bentler (2002), p.549 top line
    # I.j = nj. Delta.mu' sigma.j.inv +
    #       Delta.sigma.j' W.j Delta.sigma.j +
    #       (nj-1) Delta.sigma.w' W.w Delta.sigma.w
    #
    # where 
    # - sigma.j = sigma.w + n.j * sigma.b
    # - W.w = 1/2 * D'(sigma.w.inv %x% sigma.w.inv) D
    # - W.j = 1/2 * D'(sigma.j.inv %x% sigma.j.inv) D
    
    
}