R/ci_bound_wn_i.R

In semlbci: Likelihood-Based Confidence Interval in Structural Equation Models

#' @title Likelihood-based Confidence Bound By Wu-Neale-2012
#'
#' @description User the method proposed by Wu and Neale
#' (2012) to find the lower or upper bound of the likelihood-based
#'  confidence interval (LBCI) for one parameter in a structural
#'  equation model fitted in [lavaan::lavaan()].
#'
#' @details
#'
#' ## Important Notice
#'
#' This function is not supposed to be used directly by users in
#' typical scenarios. Its interface is user-*unfriendly* because it
#' should be used through [semlbci()]. It is exported such that
#' interested users can examine how a confidence bound is found, or
#' use it for experiments or simulations.
#'
#' ## Usage
#'
#' This function is the lowest level function used by [semlbci()].
#' [semlbci()] calls this function once for each bound of each
#' parameter. To use it, [set_constraint()] needs to be called first
#' to create the equality constraint required by the algorithm
#' proposed by Wu and Neale (2012).
#'
#' ## Algorithm
#'
#' This function implements the algorithm presented in Wu and Neale
#' (2012; see also Pek & Wu, 2015, Equation 12) that estimates all
#' free parameters in the optimization.
#'
#' ## Limitation(s)
#'
#' This function does not yet implement the method by Wu and Neale
#' (2012) for an estimate close to an attainable bound.
#'
#' @return A `cibound`-class object which is a list with three elements:
#'
#' - `bound`: A single number. The value of the bound located. `NA` is
#'            the search failed for various reasons.
#'
#' - `diag`: A list of diagnostic information.
#'
#' - `call`: The original call.
#'
#' A detailed and organized output can be printed by the default print
#' method ([print.cibound()]).
#'
#' @param i The position of the target parameter as appeared in the
#'  parameter table of an lavaan object, generated by
#'  [lavaan::parameterTable()].
#'
#' @param npar The number of free parameters, including those
#'  constrained to be equal.
#'
#' @param sem_out The fit object. Currently supports
#'  [lavaan::lavaan-class] objects only.
#'
#' @param f_constr The constraint function generated by
#' [set_constraint()].
#'
#' @param which Whether the lower bound or the upper bound is to be
#'  found. Must be `"lbound"` or `"ubound"`.
#'
#' @param history Not used. Kept for backward compatibility.
#'
#' @param perturbation_factor A number multiplied to the parameter
#'  estimates in `sem_out`. Using the parameter estimates as starting
#'  values may lead to errors in the first few iterations. Default is
#'  .90. This argument is ignored if `wald_ci_start` is `TRUE.
#'
#' @param lb_var The lower bound for free parameters that are
#'  variances. If equal to `-Inf`, the default, `lb_prop` and
#'  `lb_se_k` will be used to set the lower bounds for free variances.
#'  If it is a number, it will be used to set the lower bounds for all
#'  free variances.
#'
#' @param wald_ci_start If `TRUE`, there are no equality
#'  constraints in the model, and the target parameter is not a
#'  user-defined parameter, the Wald confidence bounds will be used as
#'  the starting value.
#'
#' @param standardized If `TRUE`, the LBCI is for the requested
#'  estimate in the standardized solution. Default is `FALSE`.
#'
#' @param opts Options to be passed to [nloptr::nloptr()], the current
#'   optimizer. Default is `list()`.
#'
#' @param ciperc The intended coverage probability for the confidence
#'  interval. Default is .95, and the bound for a 95% confidence
#'  interval will be sought.
#'
#' @param ci_limit_ratio_tol The tolerance for the ratio of `a` to
#'  `b`, where `a` is the distance between an LBCI limit and the point
#'  estimate, and the `b` is the distance between the original
#'  confidence limit (by default the Wald CI in [lavaan::lavaan()])
#'  and the point estimate. If the ratio is larger than this value or
#' smaller than the reciprocal of this value, a warning is set in the
#'  status code. Default is 1.5.
#'
#' @param verbose If `TRUE`, the function will store more diagnostic
#'   information in the attribute `diag`. Default is `FALSE`.
#'
#' @param sf A scaling factor. Used for robust confidence bounds.
#'  Default is 1. Computed by an internal function called by
#'  [semlbci()] when `robust = "satorra.2000"`.
#'
#' @param sf2 A shift factor. Used for robust confidence bounds.
#'  Default is 0. Computed by an internal function called by
#'  [semlbci()] when `robust = "satorra.2000"`.
#'
#' @param p_tol Tolerance for checking the achieved level of
#'  confidence. If the absolute difference between the achieved level
#'  and `ciperc` is greater than this amount, a warning is set in the
#'  status code and the bound is set to `NA`. Default is 5e-4.
#'
#' @param std_method The method used to find the standardized
#'  solution. If equal to `"lavaan"`,
#'  [lavaan::standardizedSolution()] will be used. If equal to
#'  `"internal"`, an internal function will be used. The `"lavaan"`
#'  method should work in all situations, but the `"internal"` method
#'  is usually much faster. Default is `"internal"`.
#'
#' @param bounds Default is `""` and this function will set the lower
#'  bounds to `lb_var` for variances. Other valid values are those
#'  accepted by [lavaan::lavaan()]. Ignored for now.
#'
#' @param xtol_rel_factor Multiply the default `xtol_rel` by a number,
#'   usually a positive number equal to or less than 1, to change the
#'   default termination criterion. Default is 1.
#'
#' @param ftol_rel_factor Multiply the default `ftol_rel` by a number,
#'   usually a positive number equal to or less than 1, to change the
#'   default termination criterion. Default is 1.
#'
#' @param lb_prop Used by an internal function to set the lower bound
#'   for free variances. Default is .05, setting the lower bound to
#'   .05 * estimate. Used only if the lower bound set by `lb_se_k` is
#'   negative.
#'
#' @param lb_se_k Used by an internal function to set the lower bound
#'   for free variances. Default is 3, the estimate minus 3 standard
#'   error. If negative, the lower bound is set using `lb_prop`.
#'
#' @param try_harder If error occurred
#' in the optimization, how many more
#' times to try. In each new attempt,
#' the starting values will be randomly
#' jittered. Default is 0.
#'
#' @param fit_lb The vector of lower
#' bounds of parameters. Default is
#' `-Inf`, setting the lower bounds to
#' `-Inf` for all parameters except for
#' free variances which are controlled
#' by `lb_var`.
#'
#' @param fit_ub The vector of upper
#' bounds of parameters. Default is
#' `+Inf`, setting the lower bounds to
#' `+Inf` for all parameters.
#'
#' @param timeout The approximate
#' maximum time for the search, in
#' second. Default is 300 seconds
#' (5 minutes).
#'
#' @param ... Optional arguments. Not used.
#'
#' @references
#'
#' Pek, J., & Wu, H. (2015). Profile likelihood-based confidence
#'  intervals and regions for structural equation models.
#'  *Psychometrika, 80*(4), 1123-1145.
#'  \doi{10.1007/s11336-015-9461-1}
#'
#' Wu, H., & Neale, M. C. (2012). Adjusted confidence intervals for a
#'  bounded parameter. *Behavior Genetics, 42*(6), 886-898.
#'  \doi{10.1007/s10519-012-9560-z}
#'
#'
#' @seealso [print.cibound()], [semlbci()], [ci_i_one()]
#'
#' @examples
#'
#' data(simple_med)
#' dat <- simple_med
#'
#' mod <-
#' "
#' m ~ x
#' y ~ m
#' "
#'
#' fit_med <- lavaan::sem(mod, simple_med, fixed.x = FALSE)
#'
#' fn_constr0 <- set_constraint(fit_med)
#'
#' out1l <- ci_bound_wn_i(i = 1,
#'                        npar = 5,
#'                        sem_out = fit_med,
#'                        f_constr = fn_constr0,
#'                        which = "lbound")
#' out1l
#'
#' @export

ci_bound_wn_i <- function(i = NULL,
                          npar = NULL,
                          sem_out = NULL,
                          f_constr = NULL,
                          which = NULL,
                          history = FALSE,
                          perturbation_factor = .9,
                          lb_var = -Inf,
                          standardized = FALSE,
                          wald_ci_start = !standardized,
                          opts = list(),
                          ciperc = .95,
                          ci_limit_ratio_tol = 1.5,
                          verbose = FALSE,
                          sf = 1,
                          sf2 = 0,
                          p_tol = 5e-4,
                          std_method = "internal",
                          bounds = "none",
                          xtol_rel_factor = 1,
                          ftol_rel_factor = 1,
                          lb_prop = .05,
                          lb_se_k = 3,
                          try_harder = 0,
                          fit_lb = -Inf,
                          fit_ub = +Inf,
                          timeout = 300,
                          ...) {
    k <- switch(which,
                lbound = 1,
                ubound = -1)
    # Check if the parameter is a user-defined parameter
    p_table <- lavaan::parameterTable(sem_out)
    i_op <- p_table[i, "op"]
    # The id in the vector of free parameters
    i_in_free <- p_table[i, "free"]

    # Get original point estimate and CI
    if (standardized) {
        ## Standardized solution
        p_est <- lavaan::standardizedSolution(sem_out,
                    type = "std.all",
                    se = TRUE,
                    zstat = FALSE,
                    pvalue = FALSE,
                    ci = TRUE,
                    level = ciperc,
                    remove.eq = FALSE,
                    remove.ineq = FALSE,
                    remove.def = FALSE,
                    output = "data.frame")
        i_est <- p_est[i, "est.std"]
        i_org_ci_lower <- p_est[i, "ci.lower"]
        i_org_ci_upper <- p_est[i, "ci.upper"]
        if (std_method == "internal") {
            # Test whether the internal function can reproduce the
            # solution from lavaan::standardizedSolution().
            p_est_int <- std_lav(lavaan::coef(sem_out),
                                  sem_out)
            p_est_test <- lavaan::parameterEstimates(
                                    sem_out,
                                    standardized = TRUE,
                                    se = FALSE,
                                    zstat = FALSE,
                                    pvalue = FALSE,
                                    ci = FALSE,
                                    cov.std = FALSE,
                                    remove.system.eq = FALSE,
                                    remove.eq = FALSE,
                                    remove.ineq = FALSE,
                                    remove.def = FALSE,
                                    remove.nonfree = FALSE)
            if (lavaan::lavTech(sem_out, "ngroups") > 1) {
                p_est_test <- p_est_test[, c("lhs", "op", "rhs",
                                              "group",
                                              "std.all")]
              } else {
                p_est_test <- p_est_test[, c("lhs", "op", "rhs",
                                              "std.all")]
              }
            p_est_int_tmp <- p_est_test
            p_est_int_tmp$std.all <- p_est_int
            p_est_test <- p_est_test[!(p_est_test$op %in% c("~1", "==")), ]
            p_est_int_tmp <- p_est_int_tmp[!(p_est_int_tmp$op %in% c("~1", "==")), ]
            if (all(p_est_test$std.all != p_est_int_tmp$std.all)) {
                stop(paste0("The internal method cannot reproduce std.all.\n",
                            "Please use std_method = ", dQuote("lavaan"), "."))
              }
          }
      } else {
        ## Unstandardized solution
        p_est <- lavaan::parameterEstimates(sem_out,
                                            se = TRUE,
                                            ci = TRUE,
                                            level = ciperc,
                                            zstat = FALSE,
                                            fmi = FALSE,
                                            rsquare = TRUE,
                                            output = "data.frame")
        i_est <- p_est[i, "est"]
        i_org_ci_lower <- p_est[i, "ci.lower"]
        i_org_ci_upper <- p_est[i, "ci.upper"]
      }

    # Create the objective function and the gradient function
    ## Standardized solution
    if (standardized && (i_op %in% c("=~", "~", "~~", ":="))) {
        if (!(std_method %in% c("lavaan", "internal"))) {
            stop(paste0("std_method must be either ",
                        dQuote("lavaan"), " or ",
                        dQuote("internal")))
          }
        if (std_method == "lavaan") {
            p_std <- lavaan::standardizedSolution(sem_out,
                                                  type = "std.all",
                                                  se = FALSE,
                                                  zstat = FALSE,
                                                  pvalue = FALSE,
                                                  ci = FALSE,
                                                  remove.eq = FALSE,
                                                  remove.ineq = FALSE,
                                                  remove.def = FALSE,
                                                  output = "data.frame")
            p_std$id <- seq_len(nrow(p_std))
            if (lavaan::lavTech(sem_out, "ngroups") > 1) {
                i_lor <- get_lhs_op_rhs(i, sem_out, more = TRUE)
                i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs", "group"))$id
              } else {
                i_lor <- get_lhs_op_rhs(i, sem_out)
                i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs"))$id
              }
            start0 <- lavaan::parameterTable(sem_out)
            # The function to be minimized.
            lbci_b_f <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
                start1 <- start0
                start1[start1$free > 0, "est"] <- param
                sem_out2 <- sem_out
                sem_out2@ParTable <- as.list(start1)
                sem_model <- sem_out2@Model
                sem_model <- lavaan::lav_model_set_parameters(sem_model,
                                          start1[start1$free > 0, "est"])
                sem_out2@Model <- sem_model
                std0 <- lavaan::standardizedSolution(sem_out2,
                                                type = "std.all",
                                                se = FALSE,
                                                zstat = FALSE,
                                                pvalue = FALSE,
                                                ci = FALSE,
                                                remove.eq = FALSE,
                                                remove.ineq = FALSE,
                                                remove.def = FALSE,
                                                output = "data.frame")
                k * std0[i_std, "est.std"]
              }
          }
        if (std_method == "internal") {
            p_std <- lavaan::parameterEstimates(
                                  sem_out,
                                  standardized = TRUE,
                                  se = FALSE,
                                  zstat = FALSE,
                                  pvalue = FALSE,
                                  ci = FALSE,
                                  cov.std = FALSE,
                                  remove.system.eq = FALSE,
                                  remove.eq = FALSE,
                                  remove.ineq = FALSE,
                                  remove.def = FALSE,
                                  remove.nonfree = FALSE)
            p_std$id <- seq_len(nrow(p_std))
            if (lavaan::lavTech(sem_out, "ngroups") > 1) {
                i_lor <- get_lhs_op_rhs(i, sem_out, more = TRUE)
                i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs", "group"))$id
              } else {
                i_lor <- get_lhs_op_rhs(i, sem_out)
                i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs"))$id
              }
            start0 <- lavaan::parameterTable(sem_out)
            # The function to be minimized.
            lbci_b_f <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
                std0 <- tryCatch(std_lav(param, sem_out),
                                 error = function(e) e)
                if (inherits(std0, "error")) {
                    if (stop_on_error) {
                        stop("Error in std_lav().")
                      } else {
                        return(Inf)
                      }
                  } else {
                    return(k * std0[i_std])
                  }
              }
          }
        # The gradient of the function to be minimized
        lbci_b_grad <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
            lavaan::lav_func_gradient_complex(
                                      lbci_b_f, param, sem_out = sem_out,
                                      debug = debug, lav_warn = lav_warn,
                                      sf = sf,
                                      sf2 = sf2
                                    )
          }
      }
    ## Unstandardized solution
    if (i_op == ":=") {
        ## User-defined parameter
        if (standardized) {
            # Not used
          } else {
            # Get the name of the defined parameter
            i_name <- p_table[i, "label"]
            # The function to be minimized.
            lbci_b_f <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
                out <- tryCatch(k * sem_out@Model@def.function(param)[i_name],
                                error = function(e) e)
                if (inherits(out, "error")) {
                    if (stop_on_error) {
                        stop("Error in cmoputing user-parameter(s).")
                      } else {
                        return(Inf)
                      }
                  } else {
                    return(out)
                  }
              }
            # The gradient of the function to be minimized
            lbci_b_grad <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
                lavaan::lav_func_gradient_complex(
                                          lbci_b_f, param, sem_out = sem_out,
                                          debug = debug, lav_warn = lav_warn,
                                          sf = sf,
                                          sf2 = sf2
                                        )
              }
          }
      } else {
        ## Free parameter
        if (standardized) {
            # Not used
          } else {
            # The function to be minimized.
            # force() may not be needed but does not harm to keep them.
            lbci_b_f <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
                force(k)
                force(i_in_free)
                k * param[i_in_free]
              }
            # The gradient of the function to be minimized
            grad_c <- rep(0, npar)
            grad_c[i_in_free] <- k
            lbci_b_grad <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
                    grad_c
              }
          }
      }

    # Set the starting values of the parameters
    if (wald_ci_start) {
        xstart <- set_start_wn(i, sem_out,
                               which = which,
                               standardized = standardized,
                               ciperc = ciperc)
        xstart <- xstart[xstart$free > 0, "est"]
      } else {
        xstart <- perturbation_factor * lavaan::coef(sem_out)
      }

    # Get SEs. For bounds
    p_table <- lavaan::parameterTable(sem_out)
    se_free <- p_table[p_table$free > 0, "se"]
    est_free <- p_table[p_table$free > 0, "est"]

    # Set lower bounds
    if (fit_lb == -Inf) {
        fit_lb <- rep(-Inf, npar)
      } else if (length(fit_lb) == 1) {
        fit_lb <- est_free - se_free * fit_lb
      }
    if (isTRUE(identical(lb_var, -Inf))) {
        # Override lb_var = -Inf for free variances
        fit_lb[find_variance_in_free(sem_out)] <- find_variance_in_free_lb(sem_out,
                                                      prop = lb_prop,
                                                      se_k = lb_se_k)
      } else {
        fit_lb[find_variance_in_free(sem_out)] <- lb_var
      }

    # Set upper bounds
    if (fit_ub == -Inf) {
        fit_ub <- rep(-Inf, npar)
      } else if (length(fit_ub) == 1) {
        fit_ub <- est_free + se_free * fit_ub
      }

    # For pmax() and pmin()
    fit_lb_na <- fit_lb
    fit_ub_na <- fit_ub
    fit_lb_na[fit_lb_na == -Inf] <- NA
    fit_ub_na[fit_ub_na == +Inf] <- NA

    # Need further testing on using lavaan bounds
    # if (!is.null(p_table$lower)) {
    #     fit_lb <- p_table$lower[p_table$free > 0]
    #   } else {
    #     if (bounds == "") {
    #         fit_lb <- rep(-Inf, npar)
    #         fit_lb[find_variance_in_free(sem_out)] <- lb_var
    #       } else {
    #         p_bounds <- set_bounds(sem_out = sem_out, bounds = bounds)
    #         fit_lb <- p_bounds$lower[p_bounds$free > 0]
    #       }
    #   }
    # if (!is.null(p_table$upper)) {
    #     fit_ub <- p_table$upper[p_table$free > 0]
    #   } else {
    #     if (bounds == "") {
    #         fit_ub <- rep(+Inf, npar)
    #       } else {
    #         p_bounds <- set_bounds(sem_out = sem_out, bounds = bounds)
    #         fit_ub <- p_bounds$upper[p_bounds$free > 0]
    #       }
    #   }

    # Optimization

    # Default options in nloptr
    opts_final <- utils::modifyList(list("algorithm" = "NLOPT_LD_SLSQP",
                        # "xtol_rel" = 1.0e-10,
                        "xtol_rel" = 1.0e-5 * xtol_rel_factor,
                        "ftol_rel" = 1.0e-5 * ftol_rel_factor,
                        "maxeval" = 500,
                        "print_level" = 0,
                        "maxtime" = timeout),
                        opts)
    try_harder <- as.integer(max(try_harder, 0))
    try_harder_count <- 0
    # Make sure the starting values are within the bounds
    xstart <- pmin(pmax(xstart,
                        fit_lb_na,
                        na.rm = TRUE),
                   fit_ub_na,
                   na.rm = TRUE)
    xstart_i <- xstart
    while(try_harder_count <= try_harder) {
        out <- tryCatch(nloptr::nloptr(
                            x0 = xstart_i,
                            eval_f = lbci_b_f,
                            lb = fit_lb,
                            ub = fit_ub,
                            eval_grad_f = lbci_b_grad,
                            eval_g_eq = f_constr,
                            opts = opts_final,
                            sem_out = sem_out,
                            lav_warn = FALSE,
                            debug = FALSE,
                            sf = sf,
                            sf2 = sf2,
                            stop_on_error = TRUE),
                        error = function(e) e)
        if (inherits(out, "error")) {
            xstart_i <- stats::runif(length(xstart_i),
                                     min = -2,
                                     max = 2) * xstart_i
            xstart_i <- pmin(pmax(xstart_i,
                                  fit_lb_na,
                                  na.rm = TRUE),
                             fit_ub_na,
                             na.rm = TRUE)
            try_harder_count <- try_harder_count + 1
          } else {
            try_harder_count <- Inf
          }
      }

    # Failed after try harder k times

    if (inherits(out, "error")) {
        search_error <- as.character(out)
        opts_error <- utils::modifyList(opts_final,
                                        list("maxeval" = 1))
        out <- nloptr::nloptr(x0 = xstart,
                              eval_f = lbci_b_f,
                              lb = fit_lb,
                              ub = fit_ub,
                              eval_grad_f = lbci_b_grad,
                              eval_g_eq = f_constr,
                              opts = opts_error,
                              sem_out = sem_out,
                              lav_warn = FALSE,
                              debug = FALSE,
                              sf = sf,
                              sf2 = sf2,
                              stop_on_error = TRUE)
      } else {
        search_error <- NULL
      }


    # Process the results
    ## Need as.numeric() to handle NA
    bound <- as.numeric(k * out$objective)
    bound_unchecked <- bound

    # Initialize the status code
    ## Any values other than 0 denotes a problem
    status <- 0

    # Check convergence
    if (out$status < 0) {
        # Verify the range of nloptr that denotes an error
        status <- 1
        check_optimization <- FALSE
        # bound <- NA
      } else {
        check_optimization <-  TRUE
      }

    # Check the limit
    i_org_ci_limit <- switch(which,
                        lbound = i_org_ci_lower,
                        ubound = i_org_ci_upper)
    if (!is.na(bound)) {
        ci_limit_ratio <- abs((bound - i_est) / (i_org_ci_limit - i_est))
        ci_limit_ratio2 <- 1 / ci_limit_ratio
        if (any(ci_limit_ratio > ci_limit_ratio_tol,
                ci_limit_ratio2 > ci_limit_ratio_tol)) {
            # status <- 1
            # Do not set the bound to NA because the limit may still be valid.
          }
      } else {
        ci_limit_ratio <- NA
      }

    # Check whether admissible
    start0 <- lavaan::parameterTable(sem_out)
    i_free <- find_free(sem_out)
    start0[i_free, "est"] <- out$solution
    start0[i_free, "start"] <- out$solution
    tmp <- sem_out@Options
    tmp$check.start <- TRUE
    tmp$check.post <- TRUE
    tmp$check.vcov <- TRUE
    tmp$warn <- FALSE
    tmp$do.fit <- FALSE
    tmp$start <- start0
    fit_final <- lavaan::lavaan(model = start0,
                                slotOptions = tmp,
                                # slotParTable = tmp2,
                                # slotModel = tmp3,
                                slotSampleStats = sem_out@SampleStats,
                                slotData = sem_out@Data)
    # fit_final <- lavaan::update(sem_out, start = start0, do.fit = FALSE,
    #                             check.start = TRUE,
    #                             check.post = TRUE,
    #                             check.vcov = TRUE)
    fit_post_check <- lavaan::lavInspect(fit_final, "post.check")
    if (!fit_post_check) {
        status <- 1
        check_post_check <- FALSE
        # bound <- NA
        # The warning should be raised by the calling function, not this one.
      } else {
        check_post_check <- TRUE
      }

    # Achieved level of confidence
    if (lavaan::lavTech(sem_out, "ngroups") > 1) {
        ntotal <- lavaan::lavTech(sem_out, "ntotal")
      } else {
        ntotal <- lavaan::lavTech(sem_out, "nobs")
      }
    fmin_org <- lavaan::lavTech(sem_out, "optim")$fx
    fmin_final <- lavaan::lavTech(fit_final, "optim")$fx
    chisq_diff <- (fmin_final - fmin_org) * 2 * ntotal / sf + sf2
    ciperc_final <- stats::pchisq(chisq_diff, 1)
    if (abs(ciperc_final - ciperc) > p_tol) {
        status <- 1
        check_level_of_confidence <- FALSE
        # bound <- NA
      } else {
        check_level_of_confidence <- TRUE
      }

    if (!all(check_optimization,
             check_post_check,
             check_level_of_confidence)) {
        bound <- NA
      }

    ## Collect diagnostic information
    diag <- list(status = status,
                 est_org = i_est,
                 ci_org_limit = i_org_ci_limit,
                 ci_limit_ratio = ci_limit_ratio,
                 fit_post_check = fit_post_check,
                 start_values = xstart,
                 bound_unchecked = bound_unchecked,
                 org_values = lavaan::coef(sem_out),
                 final_values = lavaan::coef(fit_final),
                 ciperc = ciperc,
                 ciperc_final = ciperc_final,
                 i = i,
                 i_lor = get_lhs_op_rhs(i, sem_out, more =  TRUE),
                 optim_message = out$message,
                 optim_iterations = out$iterations,
                 optim_status = out$status,
                 optim_termination_conditions = out$termination_conditions,
                 method = "wn",
                 which = which,
                 standardized = standardized,
                 check_optimization = check_optimization,
                 check_post_check = check_post_check,
                 check_level_of_confidence = check_level_of_confidence,
                 search_error = search_error
                 )
    if (verbose) {
        diag$history <- out
        diag$fit_final <- fit_final
      } else {
        diag$history <- NULL
        diag$fit_final <- NULL
      }
    if (out$status < 0) {
        # If convergence status < 0, override verbose
        diag$history <- out
      }
    out <- list(bound = bound,
                diag = diag,
                call = match.call(),
                method_name = "Wu-Neale-2012")
    class(out) <- c("cibound", class(out))
    out
  }