R/ci_bound_ur_i.R

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

#' @title Likelihood-Based Confidence
#' Bound By Root Finding
#'
#' @description Using root finding
#' 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.
#'
#' For consistency in the interface,
#' most of the arguments in
#' [ci_bound_wn_i()] are also included
#' in this function, even those not used
#' internally.
#'
#' ## Algorithm
#'
#' This function, unlike
#' [ci_bound_wn_i()], use a simple root
#' finding algorithm. Basically, it tries
#' fixing the target parameter to
#' different values until the likelihood
#' ratio test *p*-value, or the
#' corresponding chi-square difference,
#' is equal to the
#' value corresponding to the desired
#' level of confidence. (Internally,
#' the difference between the *p*-value
#' and the target *p*-value, that for
#' the chi-square difference, is the
#' function value.)
#'
#' For finding the bound, this algorithm
#' can be inefficient compared to the
#' one proposed by Wu and Neale (2012).
#' The difference can be less than
#' one second versus 10 seconds.
#' It is included as a backup algorithm
#' for parameters which are difficult
#' for the method by Wu and Neale.
#'
#' Internally, it uses [uniroot()] to
#' find the root.
#'
#' ## Limitation(s)
#'
#' This function does not handle an
#' estimate close to an attainable bound
#' using the method proposed by Wu and
#' Neale (2012). Use it for such
#' parameters with cautions.
#'
#' @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 a lavaan object,
#' generated by
#' [lavaan::parameterTable()].
#'
#' @param npar Ignored by this function.
#' Included consistency in the
#' interface.
#'
#' @param sem_out The fit object.
#' Currently supports
#' [lavaan::lavaan-class] objects only.
#'
#' @param f_constr Ignored by this
#' function. Included consistency in the
#' interface.
#'
#' @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 Ignored by
#' this function. Included consistency
#' in the interface.
#'
#' @param lb_var Ignored by this
#' function. Included consistency in the
#' interface.
#'
#' @param wald_ci_start Ignored by this
#' function. Included consistency in the
#' interface.
#'
#' @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
#' [stats::uniroot()]. 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 Ignored by this function.
#' Included consistency in the
#' interface.
#'
#' @param sf2 Ignored by this function.
#' Included consistency in the
#' interface.
#'
#' @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 Ignored by this
#' function. Included consistency in the
#' interface.
#'
#' @param xtol_rel_factor Ignored by
#' this function. Included consistency
#' in the interface.
#'
#' @param ftol_rel_factor Ignored by
#' this function. Included consistency
#' in the interface.
#'
#' @param lb_prop Ignored by this
#' function. Included consistency in the
#' interface.
#'
#' @param lb_se_k Ignored by this
#' function. Included consistency in the
#' interface.
#'
#' @param d A value used to determine
#' the width of the interval in the
#' initial search. Larger this value,
#' *narrow* the interval. Default is 5.
#' Used by [ci_bound_ur()].
#'
#' @param ... Optional arguments. Not
#' used.
#'
#' @references
#'
#' 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()]; see
#' [ci_bound_wn_i()] on the version
#' for the method by Wu and Neale
#' (2012).
#'
#' @examples
#'
#' library(lavaan)
#' data(simple_med)
#' dat <- simple_med
#' mod <-
#' "
#' m ~ x
#' y ~ m
#' "
#' fit_med <- sem(mod, simple_med, fixed.x = FALSE)
#' # Remove `opts` in real cases.
#' # The options are added just to speed up the example
#' out1l <- ci_bound_ur_i(i = 1,
#'                        sem_out = fit_med,
#'                        which = "lbound",
#'                        opts = list(use_callr = FALSE,
#'                                    interval = c(0.8277, 0.8278)))
#' out1l
#' @export

ci_bound_ur_i <- function(i = NULL,
                          npar = NULL, # Not used by ur
                          sem_out = NULL,
                          f_constr = NULL, # Not used by ur
                          which = NULL,
                          history = FALSE,
                          perturbation_factor = .9, # Can be used by sem_out_userp()
                          lb_var = -Inf, # Not used by ur. Set in sem_out
                          standardized = FALSE,
                          wald_ci_start = !standardized, # Not used by ur
                          opts = list(),
                          ciperc = .95,
                          ci_limit_ratio_tol = 1.5,
                          verbose = FALSE,
                          sf = 1, # Not used by ur
                          sf2 = 0, # Not used by ur
                          p_tol = 5e-4,
                          std_method = "internal", # Not used by ur
                          bounds = "none", # Not used by ur
                          xtol_rel_factor = 1, # Not used by ur
                          ftol_rel_factor = 1, # Not used by ur
                          lb_prop = .05, # Not used by ur
                          lb_se_k = 3, # Not used by ur
                          d = 5,
                          ...) {

    # Basic info

    # 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_limit <- switch(which,
                        lbound = p_est[i, "ci.lower"],
                        ubound = p_est[i, "ci.upper"])
      } 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_limit <- switch(which,
                                lbound = p_est[i, "ci.lower"],
                                ubound = p_est[i, "ci.upper"])
      }

    # Initial interval to search
    a <- abs(i_org_ci_limit - i_est) / d
    interval0 <- c(i_org_ci_limit - a, i_org_ci_limit + a)
    attr(interval0, "bound_start") <- i_org_ci_limit
    attr(interval0, "user_est") <- i_est

    npar <- lavaan::lavTech(sem_out, "npar")

    # Generate the user-parameter function
    est_i_func <- gen_est_i(i = i,
                            sem_out = sem_out,
                            standardized = standardized)

    # Find the root
    on.exit(try(rs$kill(), silent = TRUE), add = TRUE)
    rs <- callr::r_session$new()
    ur_opts <- list(sem_out = sem_out,
                    func = est_i_func,
                    level = ciperc,
                    which = which,
                    interval = interval0,
                    uniroot_maxiter = 500,
                    use_callr = FALSE,
                    rs = rs)
    ur_opts <- utils::modifyList(ur_opts,
                                 opts)
    out <- do.call(ci_bound_ur,
                       ur_opts)

    # Process the results

    bound <- as.numeric(out$bound)
    bound_unchecked <- bound

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

    # Check convergence
    # NOTE:
    # - Not necessary to check the code
    #   because the validity of the bound
    #   can be determined by the LR test.
    # - Find the status code of uniroot
    # - Need to be done in the call to uniroot()
    # 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

    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

    fit_final <- out$sem_out_bound

    fit_post_check <- lavaan::lavTech(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

    ciperc_final <- 1 - out$lrt[2, "Pr(>Chisq)"]
    if (abs(ciperc_final - ciperc) > p_tol) {
        status <- 1
        check_level_of_confidence <- FALSE
        # bound <- NA
      } else {
        check_level_of_confidence <- TRUE
      }

    if (!is.na(bound)) {
        check_direction <- switch(which,
                             lbound = (bound < i_est),
                             ubound = (bound > i_est))
        if (!check_direction) {
            status <- 1
          }
      } else {
        check_direction <- NA
      }

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

    coef0 <- methods::getMethod("coef",
              signature = "lavaan",
              where = asNamespace("lavaan"))(sem_out)
    coef_final <- methods::getMethod("coef",
              signature = "lavaan",
              where = asNamespace("lavaan"))(fit_final)
    xstart <- coef0
    xstart[] <- 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 = unname(xstart),
                 bound_unchecked = bound_unchecked,
                 org_values = coef0,
                 final_values = coef_final,
                 ciperc = ciperc,
                 ciperc_final = ciperc_final,
                 i = i,
                 i_lor = get_lhs_op_rhs(i, sem_out, more =  TRUE),
                 optim_message = "Nil", # out$message
                 optim_iterations = out$optimize_out$iter,
                 optim_status = NA, # out$status
                 optim_termination_conditions = "Nil", # out$termination_conditions,
                 method = "ur",
                 which = which,
                 standardized = standardized,
                 check_optimization = check_optimization,
                 check_post_check = check_post_check,
                 check_level_of_confidence = check_level_of_confidence,
                 check_direction = check_direction,
                 interval0 = interval0
                 )
    if (verbose) {
        out0 <- out
        out0$sem_out_bound <- NULL
        diag$history <- out0
        diag$fit_final <- fit_final
      } else {
        diag$history <- NULL
        diag$fit_final <- NULL
      }
    if (status != 0) {
        # If status != 0, override verbose
        diag$history <- out
      }
    out <- list(bound = bound,
                diag = diag,
                call = match.call(),
                method_name = "Root Finding")
    class(out) <- c("cibound", class(out))
    out
  }

#' @title Find a Likelihood-Based
#' Confidence Bound By Root Finding
#'
#' @description 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()]
#' using [uniroot()].
#'
#' @details
#' This function is called xby
#' [ci_bound_ur_i()]. This function is
#' exported because it is a
#' stand-alone function that can be used
#' directly for any function that
#' receives a lavaan object and returns
#' a scalar.
#'
#' The function [ci_bound_ur_i()] is a
#' wrapper of this function, with an
#' interface similar to that of
#' [ci_bound_wn_i()] and returns a
#' `cibound`-class object. The
#' user-parameter function is generated
#' internally by [ci_bound_wn_i()].
#'
#' This function, on the other hand,
#' requires users to supply the function
#' directly through the `func` argument.
#' This provides the flexibility to find
#' the bound for any function of the
#' model parameter, even one that cannot
#' be easily coded in `lavaan` model
#' syntax.
#'
#' @param sem_out The fit object.
#' Currently supports
#' [lavaan::lavaan-class] objects only.
#'
#' @param func A function that receives
#' a lavaan object and returns a scalar.
#' This function is to be used by
#' [gen_userp()] and so there are
#' special requirements on it.
#' Alternatively, it can be the output
#' of [gen_est_i()].
#'
#' @param ... Optional arguments to be
#' passed to `func`. Usually not used
#' but included in case the function
#' has such arguments.
#'
#' @param level The level of confidence
#' of the confidence interval. Default
#' is .95, or 95%.
#'
#' @param which Whether the lower bound
#' or the upper bound is to be found.
#' Must be `"lbound"` or `"ubound"`.
#'
#' @param interval A numeric vector of
#' two values, which is the initial
#' interval to be searched. If `NULL`,
#' the default, it will be determined
#' internally using Wald or delta
#' method confidence interval, if
#' available.
#'
#' @param progress Whether progress will
#' be reported on screen during the
#' search. Default is
#' `FALSE`.
#'
#' @param method The actual function to
#' be used in the search. which can only
#' be `"uniroot"`, the default, for now.
#' May include other function in the
#' future.
#'
#' @param lrt_method The method used in
#' [lavaan::lavTestLRT()]. Default is
#' `"default"`. It is automatically set
#'  to `"satorra.2000"` and cannot be
#' overridden if a scaled test statistic
#'  is requested in `sem_out`.
#'
#' @param tol The tolerance used in
#' [uniroot()], default is .005.
#'
#' @param root_target Whether the
#' chi-square difference (`"chisq"`),
#' the default, or its *p*-value
#' (`"pvalue"`) is used as the function
#' value in finding the root. Should have
#' little impact on the results.
#'
#' @param d A value used to determine
#' the width of the interval in the
#' initial search. Larger this value,
#' *narrow* the interval. Default is 5.
#'
#' @param uniroot_extendInt To be passed
#' to the argument `extendInt` of
#' [uniroot()]. Whether the interval
#' should be extended if the root is not
#' found. Default value depends on
#' the bound to be searched. Refer to
#' the help page of [uniroot()] for
#' possible values.
#'
#' @param uniroot_trace To be passed to
#' the argument `trace` of [uniroot()].
#' How much information is printed
#' during the search. Default is 0, and
#' no information is printed during the
#' search. Refer to the help page of
#' [uniroot()] for possible values.
#'
#' @param uniroot_maxiter The maximum
#' number of iteration in the search.
#' Default is 1000.
#'
#' @param use_callr Whether the
#' `callr` package will be used to
#' do the search in a separate R
#' process. Default is `TRUE`. Should
#' not set to `FALSE` if used in an
#' interactive environment unless this is
#' intentional.
#'
#' @param rs Optional. If set to
#' a persistent R process created by
#' `callr`, it will be used instead of
#' starting a new one, and it will not
#' be terminated on exit.
#'
#' @return
#' The function [ci_bound_ur()] returns
#' a list with the following elements:
#'
#' - `bound`: The bound found.
#'
#' - `optimize_out`: THe output of the
#' root finding function, [uniroot()]
#' for now. (Called `optimize_out`
#' because an earlier version of this
#' function also uses [optimize()]).
#'
#' - `sem_out_bound`: The `lavaan` model
#' with the user-defined parameter fixed
#' to the bound.
#'
#' - `lrt`: The output of
#' [lavaan::lavTestLRT()] comparing
#' `sem_out` and `sem_out_bound`.
#'
#' - `bound_start`: The Wald or delta
#' method confidence bound returned when
#' determining the interval internally.
#'
#' - `user_est`: The estimate of the
#' user-defined parameter when
#' determining the interval internally.
#'
#' @examples
#'
#' library(lavaan)
#' data(simple_med)
#' dat <- simple_med
#' mod <-
#' "
#' m ~ x
#' y ~ m
#' "
#' fit_med <- lavaan::sem(mod, simple_med, fixed.x = FALSE)
#' parameterTable(fit_med)
#' # Create a function to get the second parameter
#' est_i <- gen_est_i(i = 2, sem_out = fit_med)
#' # Find the lower bound of the likelihood-based confidence interval
#' # of the second parameter.
#' # user_callr should be TRUE or omitted in read research.
#' # Remove interval in read research. It is added to speed up the example.
#' out1l <- ci_bound_ur(sem_out = fit_med,
#'                      func = est_i,
#'                      which = "lbound",
#'                      use_callr = FALSE,
#'                      interval = c(.39070, .39075))
#' out1l
#'
#' @rdname ci_bound_ur
#'
#' @export

ci_bound_ur <- function(sem_out,
                        func,
                        ...,
                        level = .95,
                        which = c("lbound", "ubound"),
                        interval = NULL,
                        progress = FALSE,
                        method = "uniroot",
                        lrt_method = "default",
                        tol = .0005, # This is the tolerance in x, not in y
                        root_target = c("chisq", "pvalue"),
                        d = 5,
                        uniroot_extendInt = switch(which, lbound = "downX", ubound = "upX"),
                        uniroot_trace = 0,
                        uniroot_maxiter = 1000,
                        use_callr = TRUE,
                        rs = NULL) {

    # Internal Workflow
    # - Define the function that works on the lavaan object.
    # - Call add_func() to generate a modified lavaan object.
    # - Call sem_out_userp_run() to fix to a value

    # satorra.2000 is used automatically and cannot be overridden
    scaled <- any(names(lavaan::lavInspect(sem_out, "test")) %in%
                        c("satorra.bentler",
                          "yuan.bentler",
                          "yuan.bentler.mplus",
                          "mean.var.adjusted",
                          "scaled.shifted"))
    if (scaled) {
        lrt_method <- "satorra.2000"
      }

    which <- match.arg(which)
    root_target <- match.arg(root_target)
    # Only support uniroot() for now
    method <- match.arg(method)
    lrt_alpha <- 1 - level
    # Create the modified model
    fit_i <- add_func(func = func,
                      sem_out = sem_out)
    if (is.null(interval)) {
        # Set the interval to search
        if (progress) {
            cat("\nFinding the initial interval ...")
          }
        interval <- uniroot_interval(sem_out = sem_out,
                                     func = func,
                                     which = which,
                                     d = d)
        if (progress) {
            cat("\rInitial interval computed.      ")
          }
      }
    bound_start <- attr(interval, "bound_start")
    user_est <- attr(interval, "user_est")
    bound_start <- ifelse(is.null(bound_start), NA, bound_start)
    user_est <- ifelse(is.null(user_est), NA, user_est)
    if (method == "uniroot") {
        # Define the function for which to find the root
        lrtp_i <- function(x,
                            alpha,
                            root_target = "pvalue",
                            global_ok = FALSE,
                            rs = NULL) {
            sem_out_x <- sem_out_userp_run(target = x,
                                            object = fit_i,
                                            global_ok = global_ok,
                                            rs = rs)
            lrt0 <- lavaan::lavTestLRT(sem_out_x,
                                       sem_out,
                                       method = lrt_method)
            cname <- switch(root_target,
                            pvalue = "Pr(>Chisq)",
                            chisq = "Chisq diff")
            out1 <- lrt0[2, cname]
            y_target <- switch(root_target,
                                pvalue = alpha,
                                chisq = stats::qchisq(level, df = 1))
            out1 - y_target
          }
        if (use_callr) {
            if (is.null(rs)) {
                on.exit(try(rs$kill(), silent = TRUE), add = TRUE)
                rs <- callr::r_session$new()
              }
            # May fix: Cannot pass an R session to another R session
            if (progress) {
                optimize_out <- rs$call(function(f,
                                                interval,
                                                alpha,
                                                root_target,
                                                global_ok,
                                                extendInt,
                                                tol,
                                                trace,
                                                maxiter) {
                                          stats::uniroot(f = f,
                                              interval = interval,
                                              alpha = alpha,
                                              root_target = root_target,
                                              global_ok = global_ok,
                                              extendInt = extendInt,
                                              tol = tol,
                                              trace = trace,
                                              maxiter = maxiter)
                                        },
                                      args = list(f = lrtp_i,
                                                  interval = interval,
                                                  alpha = lrt_alpha,
                                                  root_target = root_target,
                                                  global_ok = TRUE,
                                                  extendInt = uniroot_extendInt,
                                                  tol = tol,
                                                  trace = uniroot_trace,
                                                  maxiter = uniroot_maxiter))
                cat("\nSearch started ...\n")
                while (rs$poll_process(500) != "ready") {
                    rtime <- rs$get_running_time()["total"]
                    # TODO:
                    # - Need to handle time longer than 1 mins.
                    cat("\r", "Time spent: ", round(rtime, 2), " sec.", spe = "")
                  }
                cat("\nSearch finished. Finalize the results ... \n")
                optimize_out <- rs$read()$result
              } else {
                optimize_out <- rs$run(function(f,
                                                interval,
                                                alpha,
                                                root_target,
                                                global_ok,
                                                extendInt,
                                                tol,
                                                trace,
                                                maxiter) {
                                          stats::uniroot(f = f,
                                              interval = interval,
                                              alpha = alpha,
                                              root_target = root_target,
                                              global_ok = global_ok,
                                              extendInt = extendInt,
                                              tol = tol,
                                              trace = trace,
                                              maxiter = maxiter)
                                        },
                                      args = list(f = lrtp_i,
                                                  interval = interval,
                                                  alpha = lrt_alpha,
                                                  root_target = root_target,
                                                  global_ok = TRUE,
                                                  extendInt = uniroot_extendInt,
                                                  tol = tol,
                                                  trace = uniroot_trace,
                                                  maxiter = uniroot_maxiter))
              }
          } else {
            # Run locally. Need global_ok == FALSE.
            if (progress) {
                cat("\nSearch started ...\n")
              }
            if (is.null(rs)) {
                on.exit(try(rs$kill(), silent = TRUE), add = TRUE)
                rs <- callr::r_session$new()
              }
            optimize_out <- stats::uniroot(lrtp_i,
                                           interval = interval,
                                           alpha = lrt_alpha,
                                           root_target = root_target,
                                           global_ok = FALSE,
                                           rs = rs,
                                           extendInt = uniroot_extendInt,
                                           tol = tol,
                                           trace = uniroot_trace,
                                           maxiter = uniroot_maxiter)
            if (progress) {
                cat("\nSearch finished. Finalize the results ... \n")
              }
        }
      }
    out_root <- optimize_out$root
    sem_out_final <- sem_out_userp_run(target = out_root,
                                       object = fit_i,
                                       rs = rs)
    lrt_final <- lavaan::lavTestLRT(sem_out_final,
                                    sem_out,
                                    method = lrt_method)
    out <- list(bound = out_root,
                optimize_out = optimize_out,
                sem_out_bound = sem_out_final,
                lrt = lrt_final,
                bound_start = bound_start,
                user_est = user_est)
    out
  }

#' @param i The position of the target
#' parameter as appeared in the
#' parameter table of an lavaan object,
#' generated by
#' [lavaan::parameterTable()].
#'
#' @param standardized If `TRUE`, the
#' standardized estimate is to be
#' retrieved. Default is `FALSE`.
#' Only support `"std.all"` for now.
#'
#' @return
#' The function [gen_est_i()] returns
#' a special function can inspects the
#' `Model` slot (and `implied` slot
#' if necessary) of a modified `lavaan`
#' object and return the parameter
#' estimate. This function is to be used
#' by [ci_bound_ur()] or
#' [gen_sem_out_userp()].
#'
#' @rdname ci_bound_ur
#'
#' @export

gen_est_i <- function(i,
                      sem_out,
                      standardized = FALSE) {
    ptable <- lavaan::parameterTable(sem_out)
    # i_user <- (ptable$free > 0) | (ptable$user > 0)
    # i_est <- which(which(i_user) == i)
    # From the help page of lavaan-class object,
    # type = "user" should return *all* parameters
    # in the parameter table, including equality constraints.
    # Therefore, the row number should be equal to the order
    # in the vector of coefficients.
    i_user <- i
    i_est <- i
    is_def <- (ptable$op[i] == ":=")
    i_lhs <- (ptable$lhs[i])
    i_rhs <- (ptable$rhs[i])
    i_op <- ptable$op[i]
    i_gp <- ptable$group[i]
    if (standardized) {
        if (is_def) {
            out <- function(object) {
                # Just in case
                force(i_est)
                std <- lavaan::standardizedSolution(object,
                                est = lavaan::lav_model_get_parameters(object@Model, type = "user"),
                                GLIST = object@Model@GLIST,
                                se = FALSE,
                                zstat = FALSE,
                                pvalue = FALSE,
                                ci = FALSE,
                                remove.eq = FALSE,
                                remove.ineq = FALSE,
                                remove.def = FALSE,
                                type = "std.all")
                unname(std[i_est, "est.std"])
              }
          } else {
            out <- function(object) {
                # Just in case
                force(i_est)
                force(i_gp)
                force(i_op)
                force(i_lhs)
                force(i_rhs)
                object@implied <- lavaan::lav_model_implied(object@Model)
                est <- lavaan::lav_model_get_parameters(object@Model, type = "user")[i_est]
                implied <- lavaan::lavInspect(object, "cov.all", drop.list.single.group = FALSE)
                implied_sd <- sqrt(diag(implied[[i_gp]]))
                est1 <- switch(i_op,
                               `~~` = est / (implied_sd[i_lhs] * implied_sd[i_rhs]),
                               `~` = est * implied_sd[i_rhs] / implied_sd[i_lhs],
                               `=~` = est * implied_sd[i_lhs] / implied_sd[i_rhs])
                unname(est1)
              }
          }
      } else {
        out <- function(object) {
            force(i_est)
            unname(lavaan::lav_model_get_parameters(object@Model, type = "user")[i_est])
          }
      }
    return(out)
  }