R/lav_modification.R

In lavaan: Latent Variable Analysis

# univariate modification indices
#

modindices <- function(object,
                       standardized = TRUE,
                       cov.std = TRUE,
                       information = "expected",

                       # power statistics?
                       power = FALSE,
                       delta = 0.1,
                       alpha = 0.05,
                       high.power = 0.75,

                       # customize output
                       sort. = FALSE,
                       minimum.value = 0.0,
                       maximum.number = nrow(LIST),
                       free.remove = TRUE,
                       na.remove = TRUE,
                       op = NULL) {

    # check if model has converged
    if(object@optim$npar > 0L && !object@optim$converged) {
        warning("lavaan WARNING: model did not converge")
    }

    # not ready for estimator = "PML"
    if(object@Options$estimator == "PML") {
        stop("lavaan WARNING: modification indices for estimator PML are not implemented yet.")
    }

    # sanity check
    if(power) {
        standardized <- TRUE
    }

    # extended list (fixed-to-zero parameters)
    strict.exo <- FALSE
    if(object@Model@conditional.x) {
        strict.exo <- TRUE
    }
    FULL <- lav_partable_full(partable = object@ParTable,
                              lavpta = object@pta,
                              free = TRUE, start = TRUE,
                              strict.exo = strict.exo)
    FULL$free <- rep(1L, nrow(FULL))
    FULL$user <- rep(10L, nrow(FULL))

    FIT <- lav_object_extended(object, add = FULL, all.free = TRUE)
    LIST <- FIT@ParTable


    # compute information matrix 'extended model'
    # ALWAYS use *expected* information (for now)
    Information <- lavTech(FIT, paste("information", information, sep = "."))

    # compute gradient 'extended model'
    score <- lavTech(FIT, "gradient.logl")

    # Saris, Satorra & Sorbom 1987
    # partition Q into Q_11, Q_22 and Q_12/Q_21
    # which elements of Q correspond with 'free' and 'nonfree' parameters?
    model.idx <- LIST$free[ LIST$free > 0L & LIST$user != 10L ]
    extra.idx <- LIST$free[ LIST$free > 0L & LIST$user == 10L ]

    # catch empty extra.idx (no modification indices!)
    if(length(extra.idx) == 0L) {
        # 2 possibilities: either model is saturated, or we have constraints
        if(object@test[[1]]$df == 0) {
            warning("lavaan WARNING: list with extra parameters is empty; model is saturated")
        } else {
            warning("lavaan WARNING: list with extra parameters is empty; to release equality\n                  constraints, use lavTestScore()")
        }
        LIST <- data.frame(lhs = character(0), op = character(0),
                           rhs = character(0), group = integer(0),
                           mi = numeric(0), epc = numeric(0),
                           sepc.lv = numeric(0), sepc.all = numeric(0),
                           sepc.nox = numeric(0))
        return(LIST)
    }

    # partition
    I11 <- Information[extra.idx, extra.idx, drop = FALSE]
    I12 <- Information[extra.idx, model.idx, drop = FALSE]
    I21 <- Information[model.idx, extra.idx, drop = FALSE]
    I22 <- Information[model.idx, model.idx, drop = FALSE]

    # ALWAYS use *expected* information (for now)
    I22.inv <- try(lavTech(object, paste("inverted.information",
                                         information, sep = ".")),
                   silent = TRUE)
    # just in case...
    if(inherits(I22.inv, "try-error")) {
        stop("lavaan ERROR: could not compute modification indices; information matrix is singular")
    }

    V <- I11 - I12 %*% I22.inv %*% I21
    V.diag <- diag(V)
    # dirty hack: catch very small or negative values in diag(V)
    # this is needed eg when parameters are not identified if freed-up;
    idx <- which(V.diag < .Machine$double.eps^(1/3)) # was 1/2 <0.6-14
    if(length(idx) > 0L) {
        V.diag[idx] <- as.numeric(NA)
    }

    # create and fill in mi
    if(object@Data@nlevels == 1L) {
        N <- object@SampleStats@ntotal
        if(object@Model@estimator %in% ("ML")) {
            score <- -1 * score # due to gradient.logl
        }
    } else {
        # total number of clusters (over groups)
        N <- 0
        for(g in 1:object@SampleStats@ngroups) {
            N <- N + object@Data@Lp[[g]]$nclusters[[2]]
        }
        #score <- score * (2 * object@SampleStats@ntotal) / N
        score <- score / 2 # -2 * LRT
    }
    mi <- numeric( length(score) )
    mi[extra.idx] <- N * (score[extra.idx]*score[extra.idx]) / V.diag
    if(length(model.idx) > 0L) {
        mi[model.idx] <- N * (score[model.idx]*score[model.idx]) / diag(I22)
    }

    LIST$mi <- rep(as.numeric(NA), length(LIST$lhs))
    LIST$mi[ LIST$free > 0 ] <- mi

    # handle equality constraints (if any)
    #eq.idx <- which(LIST$op == "==")
    #if(length(eq.idx) > 0L) {
    #    OUT <- lavTestScore(object, warn = FALSE)
    #    LIST$mi[ eq.idx ] <- OUT$uni$X2
    #}

    # scaled?
    #if(length(object@test) > 1L) {
    #    LIST$mi.scaled <- LIST$mi / object@test[[2]]$scaling.factor
    #}

    # EPC
    d <- (-1 * N) * score
    # needed? probably not; just in case
    d[which(abs(d) < 1e-15)] <- 1.0
    LIST$epc[ LIST$free > 0 ] <- mi/d


    # standardize?
    if(standardized) {

        EPC <- LIST$epc

        if(cov.std) {
            # replace epc values for variances by est values
            var.idx <- which(LIST$op == "~~" & LIST$lhs == LIST$rhs &
                             LIST$exo == 0L)
            EPC[ var.idx ] <- LIST$est[ var.idx ]
        }

        # two problems:
        #   - EPC of variances can be negative, and that is
        #     perfectly legal
        #   - EPC (of variances) can be tiny (near-zero), and we should
        #     not divide by tiny variables
        small.idx <- which(LIST$op == "~~" &
                           LIST$lhs == LIST$rhs &
                           abs(EPC) < sqrt( .Machine$double.eps ) )
        if(length(small.idx) > 0L) {
            EPC[ small.idx ] <- as.numeric(NA)
        }

        # get the sign
        EPC.sign <- sign(LIST$epc)

        LIST$sepc.lv <- EPC.sign * lav_standardize_lv(object,
                                                      partable = LIST,
                                                      est = abs(EPC),
                                                      cov.std = cov.std)
        if(length(small.idx) > 0L) {
            LIST$sepc.lv[small.idx] <- 0
        }
        LIST$sepc.all <- EPC.sign * lav_standardize_all(object,
                                                        partable = LIST,
                                                        est = abs(EPC),
                                                        cov.std = cov.std)
        if(length(small.idx) > 0L) {
            LIST$sepc.all[small.idx] <- 0
        }
        LIST$sepc.nox <- EPC.sign * lav_standardize_all_nox(object,
                                                            partable = LIST,
                                                            est = abs(EPC),
                                                            cov.std = cov.std)
        if(length(small.idx) > 0L) {
            LIST$sepc.nox[small.idx] <- 0
        }

    }

    # power?
    if(power) {
        LIST$delta <- delta
        # FIXME: this is using epc in unstandardized metric
        #        this would be much more useful in standardized metric
        #        we need a lav_standardize_all.reverse function...
        LIST$ncp <- (LIST$mi / (LIST$epc*LIST$epc)) * (delta*delta)
        LIST$power <- 1 - pchisq(qchisq((1.0 - alpha), df=1),
                                 df=1, ncp=LIST$ncp)
        LIST$decision <- character( length(LIST$power) )

        # five possibilities (Table 6 in Saris, Satorra, van der Veld, 2009)
        mi.significant <- ifelse( 1 - pchisq(LIST$mi, df=1) < alpha,
                                  TRUE, FALSE )
        high.power <- LIST$power > high.power
        # FIXME: sepc.all or epc??
        #epc.high <- abs(LIST$sepc.all) > LIST$delta
        epc.high <- abs(LIST$epc) > LIST$delta

        LIST$decision[ which(!mi.significant & !high.power)] <- "(i)"
        LIST$decision[ which( mi.significant & !high.power)] <- "**(m)**"
        LIST$decision[ which(!mi.significant &  high.power)] <- "(nm)"
        LIST$decision[ which( mi.significant &  high.power &
                             !epc.high)] <-  "epc:nm"
        LIST$decision[ which( mi.significant &  high.power &
                              epc.high)] <-  "*epc:m*"

        #LIST$decision[ which(mi.significant &  high.power) ] <- "epc"
        #LIST$decision[ which(mi.significant & !high.power) ] <- "***"
        #LIST$decision[ which(!mi.significant & !high.power) ] <- "(i)"
    }

    # remove rows corresponding to 'fixed.x' exogenous parameters
    #exo.idx <- which(LIST$exo == 1L & nchar(LIST$plabel) > 0L)
    #if(length(exo.idx) > 0L) {
    #    LIST <- LIST[-exo.idx,]
    #}

    # remove some columns
    LIST$id <- LIST$ustart <- LIST$exo <- LIST$label <- LIST$plabel <- NULL
    LIST$start <- LIST$free <- LIST$est <- LIST$se <- LIST$prior <- NULL
    LIST$upper <- LIST$lower <- NULL

    if(power) {
        LIST$sepc.lv <- LIST$sepc.nox <- NULL
    }

    # create data.frame
    LIST <- as.data.frame(LIST, stringsAsFactors = FALSE)
    class(LIST) <- c("lavaan.data.frame", "data.frame")

    # remove rows corresponding to 'old' free parameters
    if(free.remove) {
        old.idx <- which(LIST$user != 10L)
        if(length(old.idx) > 0L) {
            LIST <- LIST[-old.idx,]
        }
    }

    # remove rows corresponding to 'equality' constraints
    eq.idx <- which(LIST$op == "==")
    if(length(eq.idx) > 0L) {
        LIST <- LIST[-eq.idx,]
    }

    # remove even more columns
    LIST$user <- NULL

    # remove block/group/level is only single block
    if(lav_partable_nblocks(LIST) == 1L) {
        LIST$block <- NULL
        LIST$group <- NULL
        LIST$level <- NULL
    }

    # sort?
    if(sort.) {
        LIST <- LIST[order(LIST$mi, decreasing = TRUE),]
    }
    if(minimum.value > 0.0) {
        LIST <- LIST[!is.na(LIST$mi) & LIST$mi > minimum.value,]
    }
    if(maximum.number < nrow(LIST)) {
        LIST <- LIST[seq_len(maximum.number),]
    }
    if(na.remove) {
        idx <- which(is.na(LIST$mi))
        if(length(idx) > 0) {
            LIST <- LIST[-idx,]
        }
    }
    if(!is.null(op)) {
        idx <- LIST$op %in% op
        if(length(idx) > 0) {
            LIST <- LIST[idx,]
        }
    }

    # add header
    # TODO: small explanation of the columns in the header?
#    attr(LIST, "header") <-
# c("modification indices for newly added parameters only; to\n",
#   "see the effects of releasing equality constraints, use the\n",
#   "lavTestScore() function")

    LIST
}

# aliases
modificationIndices <- modificationindices <- modindices