R/lrt.R

In semlrtp: Likelihood Ratio Test P-Values for Structural Equation Models

#' @title Fix a Free Parameter to Zero and Do an LR Test
#'
#' @description Fix the designated
#' free parameter to zero and do a
#' likelihood ratio test.
#'
#' @details It fixes the designated
#' free parameter in a `lavaan` output,
#' refit the model, and do a likelihood
#' ratio test comparing this model with
#' the original model.
#'
#' The model to be fixed is generated
#' by [fix_to_zero()].
#'
#' If the parameter to be fixed is
#' a variance, related covariance(s),
#' if any, will also be fixed to zero.
#'
#' Users should usually call
#' [lrtp()] directly instead of calling
#' this function. It is exported for
#' developers.
#'
#' @return
#' A `lrt`-class object, which is a
#' list with the following elements:
#'
#' - `lrt`: The output of [lavaan::lavTestLRT()].
#' If there is an error message or
#' warning, it is set to `NA`.
#'
#' - `par_id`: The row number of the
#' designated free parameter.
#'
#' - `par_label`: The label of the
#' designated free parameter,
#' generated by [lavaan::lav_partable_labels()].
#'
#' - `fit1`: The original `lavaan`
#' output, if `store_fit` is `TRUE`.
#' `NA` if `store_fit` is `FALSE`,
#' the default.
#'
#' - `fix_to_zero`: The output of
#' `fit_to_zero()`.
#'
#' - `call`: The call to this function.
#'
#' - `lrt_status`: Integer. If 0, then
#' there is no error nor warning
#' in the likelihood ratio test and
#' [lavaan::lavTestLRT()] returns a
#' table (`data.frame`) of the test.
#' If -1, then something is wrong,
#' e.g., an error or warning occurred
#' when doing the likelihood ratio
#' test.
#'
#' - `lrt_msg`: If something went wrong
#' when doing the likelihood ratio
#' test, this is the error or warning
#' message when calling
#' [lavaan::lavTestLRT()].
#' If no error nor warning,
#' this is `NA`.
#'
#' @param fit A `lavaan`-class object.
#'
#' @param par_id It can be an integer.
#' or a string. If it is an integer,
#' it should be the row
#' number of the free parameter in the
#' parameter table of `fit` to be
#' fixed to zero. If it is a string,
#' it must be
#' a valid `lavaan` model syntax for
#' a parameter, or
#' the label of a labelled parameter.
#'
#' @param store_fit Logical. If `TRUE`,
#' `fit` will be stored in the output.
#' Default is `FALSE`.
#'
#' @param group If a model syntax
#' is used in `par_id` and the model
#' is a multigroup model, this should
#' be either the group label or the
#' group number of the parameter.
#'
#' @param se_keep_bootstrap Logical.
#' If `TRUE` and `fit` used
#' bootstrapping standard error
#' (with `se = "bootstrap"`), bootstrapping
#' will also be use in fitting the
#' restricted model. If `FALSE`, the
#' default, then `se` will be set
#' to `"standard"` if it is `"bootstrap"`
#' in `fit`, to speed up the computation.
#'
#' @param LRT_method String. Passed to
#' the `method` argument of
#' [lavaan::lavTestLRT()]. Default is
#' `"default"`, and let
#' [lavaan::lavTestLRT()] decide the
#' method based on `fit`.
#'
#' @param scaled.shifted Logical.
#' Used when the method used in
#' [lavaan::lavTestLRT()] is
#' `"satorra.2000"`. Default is
#' `TRUE` and a scaled and shifted
#' test statistic is used, the same
#' default of [lavaan::lavTestLRT()].
#'
#' @param fallback_method The default
#' method of [lavaan::lavTestLRT()],
#' `"satorra.bentler.2001"`,
#' may sometimes fail. If failed,
#' this function will call
#' [lavaan::lavTestLRT()]
#' again using `fallback_method`. which
#' is `"satorra.2000"` by default.
#'
#' @seealso [print.lrt()] for its
#' print-method, and [lrtp()] for the
#' main function.
#'
#' @author Shu Fai Cheung <https://orcid.org/0000-0002-9871-9448>
#'
#' @examples
#' library(lavaan)
#' data(data_sem16)
#' mod <-
#' "
#' f1 =~ x1 + x2 + x3
#' f2 =~ x4 + x5 + x6
#' "
#' fit <- sem(mod, data_sem16)
#' # Fix the factor covariance to zero
#' out <- lrt(fit, par_id = 15)
#' out$lrt
#' parameterEstimates(fit)[15, ]
#' parameterEstimates(out$fix_to_zero$fit0)[15, ]
#'
#' # Can use model syntax for par_id
#'
#' out <- lrt(fit, par_id = "f1 =~ x3")
#' out$lrt
#'
#' @export

lrt <- function(fit,
                par_id,
                store_fit = FALSE,
                group = NULL,
                se_keep_bootstrap = FALSE,
                LRT_method = "default",
                scaled.shifted = TRUE,
                fallback_method = "satorra.2000") {
    if (isFALSE(inherits(fit, "lavaan"))) {
        stop("The fit object is not a lavaan object.")
      }
    if (missing(par_id)) {
        stop("par_id must be supplied.")
      }
    if (!isTRUE(is.numeric(par_id))) {
        par_id <- lhs_to_par_id(fit = fit,
                                par = par_id,
                                group = group)
      }
    fit_i <- fix_to_zero(fit,
                         par_id = par_id,
                         se_keep_bootstrap = se_keep_bootstrap)
    fit0 <- fit_i$fit0
    lrt_out <- NA
    lrt_msg <- NA
    # fit0 is NA if
    # - Failed to converge (fit0_converged == FALSE), or
    # - VCOV warning (fit0_vcov_ok == FALSE), or
    # - DF change > 1 and the par is not a variance
    if (inherits(fit0, "lavaan")) {
        lrt_msg <- tryCatch(lrt_out <- lavaan::lavTestLRT(fit,
                                                          fit0,
                                                          method = LRT_method,
                                                          scaled.shifted = scaled.shifted),
                               error = function(e) e,
                               warning = function(w) w)
        if (inherits(lrt_msg, "warning")) {
            tmp <- as.character(lrt_msg)
            # Satorra-Bentler 2001 failed. Use Satorra 2000
            if (grepl("scaling factor is negative",
                      tmp, fixed = TRUE)) {
                lrt_msg <- tryCatch(lrt_out <-  lavaan::lavTestLRT(fit,
                                                                   fit0,
                                                                   method = fallback_method,
                                                                   scaled.shifted = scaled.shifted),
                                      error = function(e) e,
                                      warning = function(w) w)
              }
          }
      }
    lrt_status <- c(NotOK = -1)
    # lrt_out is NA if
    # - LRT failed (error or warning), or
    # - LRT is not an ANOVA table for whatever reason
    if (inherits(lrt_out, "anova")) {
        if (inherits(lrt_msg, "error") ||
            inherits(lrt_msg, "warning")) {
              lrt_out <- NA
              lrt_status <- c(NotOK = -1)
            } else {
              lrt_status <- c(OK = 0)
            }
      } else {
        lrt_out <- NA
      }
    pt <- lavaan::parameterTable(fit)
    tmp <- lavaan::lav_partable_labels(pt)
    par_label <- tmp[par_id]
    out <- list(lrt = lrt_out,
                par_id = par_id,
                par_label = par_label,
                fit1 = ifelse(store_fit,
                              fit,
                              NA),
                fix_to_zero = fit_i,
                call = match.call(),
                lrt_status = lrt_status,
                lrt_msg = lrt_msg)
    class(out) <- c("lrt", class(out))
    out
  }