R/gen_userp.R

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

#' @title Create a Wrapper To Be Used
#' in 'lavaan' Models
#'
#' @description Make a function on
#' `lavaan` object usable in a `lavaan`
#' model syntax.
#'
#' @details
#'
#' ## `gen_userp`
#'
#' There are cases in which we want to
#' create a user parameter which is a
#' function of other free parameters,
#' computed by a function. However such
#' a function may work only on a
#' `lavaan` object.
#'
#' If the target function works by
#' extracting parameter estimates stored
#' in the `Model` slot and/or the
#' `implied` slot, then [gen_userp()]
#' can be used to convert it to a
#' function that retrieves the parameter
#' estimates when being called by
#' [lavaan::lavaan()] or its wrappers,
#' modifies the stored
#' `lavaan` object using
#' [lavaan::lav_model_set_parameters()]
#' and [lavaan::lav_model_implied()] to
#' change the estimates, and call the
#' target function.
#'
#' Note that this is an unconventional
#' way to define a user parameter and
#' the generated function should always
#' be checked to see whether it works as
#' expected.
#'
#' As shown in the examples, the
#' parameter computed this may not
#' have standard error nor *p*-value.
#'
#' The main purpose is for the point
#' estimate, for searching the
#' likelihood-based confidence bound
#' using [ci_bound_ur()] and
#' [ci_bound_ur_i()].
#'
#' Note that the target function
#' specified in `func` should work
#' directly on the parameter estimates
#' stored in the `Model` slot and then
#' get the estimates using
#' [lavaan::lav_model_get_parameters()].
#' Functions that work on the unmodified
#' output generated by
#' [lavaan::lavaan()] usually do not
#' work.
#'
#' Users are not recommended to use
#' [gen_userp()] and [gen_sem_out_userp()]
#' directly because they require
#' unconventional way to extract
#' parameter estimates from a lavaan
#' model. However, developers may use
#' them to include functions
#' they wrote in a lavaan model. This
#' is the technique used by
#' [ci_bound_ur_i()] to constrain any
#' parameter in a model to an arbitrary
#' value.
#'
#' @param func A function that receives a
#' `lavaan`-object and returns a scalar.
#' See Details on the restriction on this
#' function.
#'
#' @param sem_out A `lavaan`-class object
#' to be modified.
#'
#' @return
#'
#' ## `gen_userp`
#'
#' It returns a function that
#' accepts a numeric vector of length
#' equals to the number of free parameters
#' in `sem_out`, and returns a scalar
#' which is the output of `func`. If this
#' vector is not supplied, it will try to
#' find it in the `parent.frame()`. This
#' is how it works inside a `lavaan`
#' model.
#'
#' @examples
#'
#' library(lavaan)
#'
#' data(simple_med)
#' dat <- simple_med
#' mod <-
#' "
#' m ~ a*x
#' y ~ b*m
#' ab := a*b
#' "
#' fit_med <- sem(mod, simple_med, fixed.x = FALSE)
#' parameterEstimates(fit_med)
#'
#' # A trivial example for verifying the results
#' my_ab <- function(object) {
#'     # Need to use lav_model_get_parameters()
#'     # because the object is only a modified
#'     # lavaan-object, not one directly
#'     # generated by lavaan function
#'     est <- lavaan::lav_model_get_parameters(object@Model, type = "user")
#'     unname(est[1] * est[2])
#'   }
#'
#' # Check the function
#' my_ab(fit_med)
#' coef(fit_med, type = "user")["ab"]
#'
#' # Create the function
#' my_userp <- gen_userp(func = my_ab,
#'                       sem_out = fit_med)
#' # Try it on the vector of free parameters
#' my_userp(coef(fit_med))
#'
#' # Generate a modified lavaan model
#' fit_userp <- gen_sem_out_userp(userp = my_userp,
#'                                userp_name = "my_userp",
#'                                sem_out = fit_med)
#'
#' # This function can then be used in the model syntax.
#'
#' # Note that the following example only work when called inside the
#' # workspace or inside other functions such as ci_bound_ur()`
#' # and `ci_bound_ur_i()` because `lavaan::sem()` will
#' # search `my_userp()` in the global environment.
#'
#' # Therefore, the following lines are commented out.
#' # They should be run only in a "TRUE" interactive
#' # session.
#'
#' # mod2 <-
#' # "
#' # m ~ x
#' # y ~ m
#' # ab := my_userp()
#' # "
#' # fit_med2 <- sem(mod2, simple_med, fixed.x = FALSE)
#' # parameterEstimates(fit_med2)
#' #
#' # # Fit the model with the output of the function, a*b
#' # # fixed to .50
#' #
#' # fit_new <- fit_userp(.50)
#' #
#' # # Check if the parameter ab is fixed to .50
#' # parameterEstimates(fit_new)
#'
#'
#' @rdname gen_userp
#'
#' @export

gen_userp <- function(func,
                      sem_out) {
    stopifnot(is.function(func))
    stopifnot(inherits(sem_out, "lavaan"))
    out <- function(.x., ...) {
        if (missing(.x.)) {
            .x. <- get(".x.", parent.frame())
          }
        # Just in case ...
        force(sem_out)

        sem_out@Model <- lavaan::lav_model_set_parameters(sem_out@Model, .x.)
        sem_out@implied <- lavaan::lav_model_implied(sem_out@Model)
        f0 <- func(sem_out)
        if (is.na(f0) || is.nan(f0)) {
            f0 <- Inf
          }
        f0
      }
    out
  }

# @title Generate a Function to Fix a
# User-Defined Parameter
#
# @description Add a function to a
# `lavaan` object as a user-defined
# parameter and return a function which
# can fix this user-defined parameter
# to a value.
#
#' @details
#'
#' ## `gen_sem_out_userp`
#'
#' The function [gen_sem_out_userp()]
#' is to be used internally
#' for generating a function for searching
#' a likelihood-based confidence bound.
#' It is exported because it needs to
#' be run in an fresh external R process,
#' usually created by `callr` in other
#' internal functions.
#'
#' @param userp A function that is
#' generated by [gen_userp()].
#'
#' @param userp_name The name of the
#' function `userp` to be used in the
#' `lavaan` model. It does not have to
#' be the name of the function in
#' `userp`. Should be changed only if it
#' conflicts with another object in the
#' parent environment, which should not
#' happen if the model is always fitted
#' in a clean R session.
#'
#' @param fix If `TRUE`, the default,
#' the function generated is used to fix
#' the value of `userp` to a target
#' value using an equality constraint.
#' If `FALSE`, then the function simply
#' fits the model to the data.
#'
#' @param control_args To be passed to
#' the argument of the same name in
#' [lavaan::lavaan()]. Default is
#' `list()`. Can be used to set the
#' default values of this argument
#' in the generated function.
#'
#' @param iter.max The maximum number of
#' iteration when the generated function
#' fit the model. Default is 10000.
#'
#' @param max_attempts If the initial
#' fit with the equality constraint
#' fails, how many more attempts
#' will be made by the generated
#' function. Default is 5.
#'
#' @return
#'
#' ## `gen_sem_out_userp`
#'
#' If `fix` is `TRUE`, it returns a
#' function with these arguments:
#'
#' - `target`: The value to which the
#' user-defined parameter will be fixed
#' to.
#'
#' - `verbose`: If `TRUE`, additional
#' information will be printed when
#' fitting the model.
#'
#' - `control`: The values to be passed
#' as a list to the argument of the
#' same name in [lavaan::lavaan()].
#'
#' - `seed`: Numeric. If supplied, it
#' will be used in [set.seed()] to
#' initialize the random number
#' generator. Necessary to reproduce
#' some results because random numbers
#' are used in some steps in `lavaan`.
#' If `NULL`, the default, [set.seed()]
#' will not be called.
#'
#' If `fix` is `FALSE, then it returns a
#' function with optional arguments that
#' will be ignored, Calling it will
#' simply fit the modified model to the
#' data. Useful for getting the value of
#' the user-defined parameter.
#'
#' @rdname gen_userp
#'
#' @export

gen_sem_out_userp <- function(userp,
                              sem_out,
                              userp_name = "semlbciuserp1234",
                              fix = TRUE,
                              control_args = list(),
                              iter.max = 10000,
                              max_attempts = 5) {

    # Generate a function that:
    # - refits the model with the user-parameter fixed to a target value.
    # CHECKED: Identical to the experimental version
    # For add_fun() using callr.

    iter.max <- max(lavaan::lavInspect(sem_out, "optim")$iterations,
                    iter.max)
    ptable <- lavaan::parameterTable(sem_out)
    # Create a unique label: userp_label
    labels <- unique(ptable$label)
    labels <- labels[nchar(labels) > 0]
    userp_label <- make.unique(c(labels, "user"))
    userp_label <- userp_label[length(userp_label)]
    # Generate the user parameter row using userp()
    user_pt1 <- data.frame(lhs = userp_label[length(userp_label)],
                           op = ":=",
                           rhs = paste0(userp_name, "()"),
                           user = 1,
                           block = 0,
                           group = 0,
                           free = 0,
                           ustart = NA,
                           label = userp_label,
                           exo = 0)
    user_pt1 <- lavaan::lav_partable_complete(partable = user_pt1,
                                              start = TRUE)
    # Add the user parameter to the model
    ptable1 <- lavaan::lav_partable_add(partable = ptable,
                                        add = user_pt1)
    ptable1$se[which(ptable1$label == userp_label)] <- NA

    # Extract slots to be used
    slot_opt <- sem_out@Options
    slot_dat <- sem_out@Data

    if (fix) {
        # Create the row of equality constraint
        user_pt2 <- data.frame(lhs = userp_label,
                               op = "==",
                               rhs = NA,
                               user = 1,
                               block = 0,
                               group = 0,
                               free = 0,
                               ustart = NA,
                               exo = 0)
        user_pt2 <- lavaan::lav_partable_complete(partable = user_pt2,
                                                  start = TRUE)
        # Add the equality constraint to the parameter tahle
        ptable2 <- lavaan::lav_partable_add(partable = ptable1,
                                            add = user_pt2)
        # Store the position of the constraint
        i_eq <- which((ptable2$lhs == userp_label) &
                      (ptable2$op == "=="))

        # Fix the ids of the rows
        ptable2$id <- seq_along(ptable2$id)

        # Set the options
        slot_opt$do.fit <- TRUE
        # The "mplus" version does not take into account the user-parameter
        slot_opt$start <- "simple"
        # Heywood cases are allowed
        slot_opt$check.start <- FALSE
        slot_opt$check.post <- FALSE
        # Disable some options for efficiency
        slot_opt$se <- "none"
        # slot_opt$h1 <- FALSE
        # slot_opt$baseline <- FALSE

        out <- function(target,
                        verbose = FALSE,
                        control = list(),
                        seed = NULL) {
            if (!is.null(seed)) set.seed(seed)
            # Just in case ...
            force(control_args)
            force(slot_opt)
            force(slot_dat)
            force(ptable2)
            force(userp_label)
            force(i_eq)

            ptable3 <- ptable2
            ptable3$rhs[i_eq]  <- target
            control_args <- utils::modifyList(control,
                                              list(iter.max = iter.max))
            slot_opt$control <- control_args

            # Fit once to update the model slot
            attempted <- 0
            fit_new_converged <- FALSE
            if (verbose) {
                cat("Target:", target, "\n")
              }
            if (verbose) {
                cat("First attempt ...\n")
              }
            # Fail when
            # eq.constraints.K is 0x0
            # Succeed when
            # fit_new@Model@eq.constraints, or
            # length(fit_new@Model@ceq.nonlinear.idx) > 0
            slot_opt_tmp <- slot_opt
            slot_opt_tmp$do.fit <- FALSE
            eq_not_ok <- TRUE
            while (eq_not_ok) {
                fit_new <- lavaan::lavaan(model = ptable3,
                                          slotOptions = slot_opt_tmp,
                                          slotData = slot_dat)
                eq_not_ok <- !fit_new@Model@eq.constraints &&
                             !(length(fit_new@Model@ceq.nonlinear.idx) > 0)
              }

            # Extract the updated slots
            slot_pat_new <- fit_new@ParTable
            slot_model_new <- fit_new@Model
            slot_smp_new <- fit_new@SampleStats

            # Fit the model with the user parameter constrained
            # to the target value
            fit_new <- tryCatch(lavaan::lavaan(slotParTable = slot_pat_new,
                                               slotModel = slot_model_new,
                                               slotSampleStats = slot_smp_new,
                                               slotOptions = slot_opt,
                                               slotData = slot_dat),
                                error = function(e) e)
            if (inherits(fit_new, "error")) {
                stop("Something's with the function.",
                     "Error occurred when fitting the modified model.")
              }

            fit_new_converged <- lavaan::lavInspect(fit_new, "converged")

            # Try max_attempts more times
            if (!fit_new_converged) {
                if (verbose) {
                    cat("More attempts ...\n")
                  }
                while (!((attempted > max_attempts) ||
                          fit_new_converged)) {
                    if (verbose) {cat("Attempt", attempted, "\n")}
                    ptable_new <- fit_new@ParTable
                    i_free <- ptable_new$free > 0
                    ptable_new$est[i_free] <- jitter(ptable_new$est[i_free],
                                                     factor = 50)
                    slot_opt$start <- as.data.frame(ptable_new)
                    fit_old <- fit_new
                    fit_new <- tryCatch(lavaan::lavaan(slotParTable = ptable_new,
                                                       slotModel = slot_model_new,
                                                       slotSampleStats = slot_smp_new,
                                                       slotOptions = slot_opt,
                                                       slotData = slot_dat),
                                        error = function(e) e)
                    if (inherits(fit_new, "error")) {
                        fit_new <- fit_old
                        fit_new_converged <- FALSE
                      } else {
                        fit_new_converged <- lavaan::lavInspect(fit_new, "converged")
                      }
                    attempted <- attempted + 1
                  }
              }
            if (!fit_new_converged) {
                if (verbose) {
                    cat("The modified model failed to converge.")
                  }
              }
            if (verbose) {
                cat("Done!\n")
              }
            fit_new
          }
      } else {
        # The arguments are not used in this version
        out <- function(target,
                        verbose = FALSE,
                        control = list(),
                        seed = NULL) {
            # Generate a model with the user parameter added.
            # For debugging
            # Just in case ...
            force(slot_opt)
            force(slot_dat)
            force(ptable1)
            slot_opt$do.fit <- TRUE
            fit_new <- lavaan::lavaan(model = ptable1,
                                      slotOptions = slot_opt,
                                      slotData = slot_dat)
            fit_new
          }
      }
    out
  }