R/confint_indirect.R

In manymome: Mediation, Moderation and Moderated-Mediation After Model Fitting

#' @title Confidence Interval of
#' Indirect Effect or Conditional
#' Indirect Effect
#'
#' @description Return the
#' confidence interval of the indirect
#' effect or conditional indirect effect
#' stored in the output of
#' [indirect_effect()] or
#' [cond_indirect()].
#'
#' @details It extracts and returns the
#' stored confidence interval
#' if available.
#'
#' The type of confidence interval
#' depends on the call used to
#' compute the effect. This function
#' merely retrieves the stored estimates,
#' which could be generated by
#' nonparametric bootstrapping,
#' Monte Carlo simulation, or other
#' methods to be supported in
#' the future, and uses them to form the
#' percentile confidence interval.
#'
#' If the following conditions are met, the
#' stored standard errors, if available,
#' will be used to test an effect and
#' form its confidence interval:
#'
#' - Confidence intervals have not been
#'  formed (e.g., by bootstrapping or
#'  Monte Carlo).
#'
#' - The path has no mediators.
#'
#' - The model has only one group.
#'
#' - The path is moderated by one or
#'  more moderators.
#'
#' - Both the `x`-variable and the
#'  `y`-variable are not standardized.
#'
#' If the model is fitted by OLS
#' regression (e.g., using [stats::lm()]),
#' then the variance-covariance matrix
#' of the coefficient estimates will be
#' used, and confidence
#' intervals are computed from the *t*
#' statistic.
#'
#' If the model is fitted by structural
#' equation modeling using `lavaan`, then
#' the variance-covariance computed by
#' `lavaan` will be used,
#' and confidence intervals are computed
#' from the *z* statistic.
#'
#' ## Caution
#'
#' If the model is fitted by structural
#' equation modeling and has moderators,
#' the standard errors, *p*-values,
#' and confidence intervals computed
#' from the variance-covariance matrices
#' for conditional effects
#' can only be trusted if all covariances
#' involving the product terms are free.
#' If any of them are fixed, for example,
#' fixed to zero, it is possible
#' that the model is not invariant to
#' linear transformation of the variables.
#'
#' @param object The output of
#' [indirect_effect()] or
#' [cond_indirect()].
#'
#' @param parm Ignored because the
#' stored object always has only one
#' parameter.
#'
#' @param level If set to `NULL`, the default,
#' then the
#' level of confidence used to generate
#' `object`
#' is used. If set to a value,
#' this value will be used
#' to recompute the confidence intervals.
#' If the confidence interval is to be
#' computed from the standard error,
#' and so `level` is not set in `object`,
#' then the default value is .95.
#' (This new behavior applies to 0.3.6.15
#' and later versions.)
#'
#' @param boot_type If bootstrap
#' confidence interval is to be formed,
#' the type of bootstrap confidence
#' interval. The supported types
#' are `"perc"` (percentile bootstrap
#' confidence interval, the recommended
#' method) and `"bc"`
#' (bias-corrected, or BC, bootstrap
#' confidence interval). If not supplied,
#' the stored `boot_type` will be used.
#'
#' @param ...  Additional arguments.
#' Ignored by the function.
#'
#' @return A numeric vector of
#' two elements, the limits of
#' the confidence interval.
#'
#' @seealso [indirect_effect()] and
#' [cond_indirect()]
#'
#' @examples
#'
#' dat <- modmed_x1m3w4y1
#'
#' # Indirect Effect
#'
#' library(lavaan)
#' mod1 <-
#' "
#' m1 ~ x
#' m2 ~ m1
#' y  ~ m2 + x
#' "
#' fit <- sem(mod1, dat,
#'            meanstructure = TRUE, fixed.x = FALSE,
#'            se = "none", baseline = FALSE)
#' # R should be at least 2000 or 5000 in real research.
#' out1 <- indirect_effect(x = "x", y = "y",
#'                         m = c("m1", "m2"),
#'                         fit = fit,
#'                         boot_ci = TRUE, R = 40, seed = 2345,
#'                         parallel = FALSE,
#'                         progress = FALSE)
#' out1
#' confint(out1)
#'
#'
#' @export


confint.indirect <- function(
  object,
  parm,
  level = NULL,
  boot_type,
  ...
) {

    if (missing(boot_type)) {
        ci_boot_type <- object$boot_type
      } else {
        ci_boot_type <- boot_type
      }
    has_ci <- FALSE
    level_default <- .95
    if (isTRUE(!is.null(object$boot_ci))) {

        # ==== Boot CI found ====

        has_ci <- TRUE
        old_ci <- object$boot_ci
        ci_type <- "boot"
        ind_i <- object$boot_indirect
        if (is.null(level)) {
          level <- object$level
        }
        if ((level == object$level) &&
            (ci_boot_type == object$boot_type)) {
            # Retrieve stored CI if
            # no change to level and boot_type
            new_ci <- FALSE
          } else {
            # Recompute CI if
            # level or boot_type changed
            new_ci <- TRUE
          }
      }
    if (isTRUE(!is.null(object$mc_ci))) {

        # ==== Monte Carlo CI found ====

        has_ci <- TRUE
        old_ci <- object$mc_ci
        ci_type <- "mc"
        ind_i <- object$mc_indirect
        if (is.null(level)) {
          level <- object$level
        }
        if (level == object$level) {
            # Retrieve stored CI if
            # no change to level
            new_ci <- FALSE
          } else {
            # Recompute CI if
            # level changed
            new_ci <- TRUE
          }
      }

    if (!has_ci &&
        is.null(level)) {
      # No stored CI and so no stored level
      # Set default level if level is NULL
      level <- level_default
    }
    # level should have a value up to this point
    # - If set, always used.
    # - If not set, it depends on the type of CI.

    # ==== SE CI: Try ====

    se_out <- cond_effect_original_se(object,
                                      level = level)
    has_original_se <- !is.null(se_out)
    has_m <- isTRUE(!is.null(object$m))
    standardized_x <- object$standardized_x
    standardized_y <- object$standardized_y
    has_groups <- !is.null(object$group_number)
    has_wlevels <- !is.null(object$wvalues)
    if (!has_ci &&
        !has_m &&
        !has_groups &&
        has_wlevels &&
        !standardized_x &&
        !standardized_y &&
        has_original_se) {

        # ==== Use SE CI ====

        ci_type <- "se"
        has_ci <- TRUE
      }

    # ==== has_ci? ====

    if (has_ci) {
        if (ci_type == "se") {
            # Use SE CI
            out0 <- c(se_out$cilo, se_out$cihi)
          } else {
            # Boot CI or MC CI
            if (new_ci) {

                # ==== Recompute CI (new_ci TRUE) ====

                out0 <- boot_ci_internal(t0 = object$indirect,
                                t = ind_i,
                                level = level,
                                boot_type = ifelse(ci_type == "boot",
                                                  ci_boot_type,
                                                  "perc"),
                                add_names = FALSE)
              } else {

                # ==== Retrieved stored CI ====

                out0 <- old_ci
              }
          }
      } else {
        warning("Confidence interval not in the object.")
        out0 <- c(NA, NA)
      }
    # Borrowed from stats::confint()
    probs <- c((1 - level) / 2, 1 - (1 - level) / 2)
    cnames <- paste(format(100 * probs,
                           trim = TRUE,
                           scientific = FALSE,
                           digits = 2), "%")
    if (has_ci) {
        if (ci_type == "boot") {
            tmp <- switch(ci_boot_type,
                          perc = "Percentile: ",
                          bc = "Bias-Corrected: ")
            cnames <- paste0(tmp, cnames)
          }
        if (ci_type == "mc") {
            cnames <- paste0("Monte Carlo: ", cnames)
          }
      }
    rnames <- paste0(object$y, "~", object$x)
    out <- array(data = out0,
                 dim = c(1, 2),
                 dimnames = list(rnames, cnames))
    out
  }

#' @noRd

cond_effect_original_se <- function(object,
                                    level = .95) {
    if (!is.null(object$m)) {
        return(NULL)
      }
    if (is.null(object$original_se)) {
        return(NULL)
      }
    est <- object$indirect
    se <- unname(object$original_se)
    if (is.na(se)) {
        return(NULL)
      }
    dfres <- object$df_residual
    test_stat <- est / se
    p <- 2 * stats::pt(abs(test_stat),
                       df = dfres,
                       lower.tail = FALSE)
    sig <- stats::symnum(p,
                         corr = FALSE,
                         na = FALSE,
                         cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                         symbols = c("***", "**", "*", " "))
    z_crit <- -1 * stats::qt((1 - level) / 2,
                             df = dfres,
                             lower.tail = TRUE)
    cilo <- est - z_crit * se
    cihi <- est + z_crit * se
    out <- list(se = se,
                stat = test_stat,
                pvalue = p,
                sig = sig,
                cilo = cilo,
                cihi = cihi)
    attr(out, "sig_legend") <- attr(sig, "legend")
    attr(out, "original_se_level") <- level
    return(out)
  }