R/stdmod_lavaan.R

In stdmod: Standardized Moderation Effect and Its Confidence Interval

#' @title Standardized Moderation Effect and Its Bootstrap CI in 'lavaan'
#'
#' @description Compute the standardized moderation effect in a structural
#'              equation model fitted by [lavaan::lavaan()] or its wrappers and
#'              form the nonparametric bootstrap confidence interval.
#'
#' @details
#' ## Important Notes
#'
#' Starting from Version 0.2.7.5, of [stdmod_lavaan()] adopts an approach
#' to bootstrapping different from that in the previous
#' versions (0.2.7.4 and before),
#' yielding bootstrapping results different from those in
#' previous versions (for reasons explained later).
#'
#' To reproduce results from the older version of this function,
#' set `use_old_version` to `TRUE`.
#'
#' ## How it works
#'
#' [stdmod_lavaan()] accepts a [lavaan::lavaan-class] object, the
#' structural equation model output returned
#' by [lavaan::lavaan()] and its wrappers (e.g, [lavaan::sem()]) and computes
#' the standardized moderation effect using the formula in the appendix of
#' Cheung, Cheung, Lau, Hui, and Vong (2022).
#'
#' The standard deviations of the focal variable (the variable with its effect
#' on the outcome variable being moderated), moderator, and outcome
#' variable (dependent variable) are computed from the implied
#' covariance matrix returned by
#' [lavaan::lavInspect()]. Therefore, models fitted to data sets with missing
#' data (e.g., with `missing = "fiml"`) are also supported.
#'
#' Partial standardization can also be requested. For example, standardization
#' can be requested for only the focal variable and the outcome variable.
#'
#' There are two ways to request nonparametric bootstrap
#' confidence interval. First, the model is fitted with `se = "boot"`
#' or `se = "bootstrap"` in `lavaan`. The stored bootstrap
#' estimates will then be retrieved automatically to compute the
#' standardized moderation effect. This is the most efficient approach
#' if the bootstrap confidence intervals are also needed for
#' other parameters in the model. Bootstrapping needs to be
#' done only once.
#'
#' Second, bootstrap estimates can be generated by [manymome::do_boot()].
#' The output is then supplied through the argument `boot_out`.
#' Bootstrapping also only needs to be done once. This approach
#' is appropriate when bootstrapping confidence intervals are
#' not needed for other model parameters, or another type
#' of confidence interval is needed when fitting the model.
#' Please refer to the help page of [manymome::do_boot()]
#' on how to use this function.
#'
#' In both approaches, the standard deviations are also computed
#' in each bootstrap samples. This ensures that the sampling
#' variability of the standard deviations is also taken into
#' account in computing the bootstrap confidence interval of
#' the standardized moderation effect.
#'
#' ## Note on the differences between the current version (Version 0.2.7.5 or later) and previous versions (0.2.7.4 and before)
#'
#' In older versions, [stdmod_lavaan()] does not allow for
#' partial standardization. Moreover,
#' it uses [boot::boot()] to do the bootstrapping. Even with
#' the same seed, the results from [boot::boot()] are not
#' identical to those of `lavaan` with `se = "boot"`
#' because they differ in the way the indices
#' of resampling are generated. Both approaches are correct,
#' They just use the generated random numbers differently.
#' To have results consistent with those from `lavaan`,
#' the current version of [stdmod_lavaan()] adopts a
#' resampling algorithm identical to that of `lavaan`.
#' Last, in older versions, [stdmod_lavaan()] does
#' bootstrapping every time it is called. This is
#' inefficient.
#'
#' The bootstrapping results in the current version are not
#' identical to those in older versions due to the use
#' of different resampling algorithms, To reproduce
#' previous results, set `use_old_version` to `TRUE`
#'
#' @return
#' A list of class `stdmod_lavaan` with these elements:
#'
#'  - `stdmod`: The standardized moderation effect.
#'
#'  - `ci`: The nonparametric bootstrap confidence interval. `NA` if
#'            confidence interval not requested.
#'
#'  - `boot_out`: The raw output from [boot::boot()]. `NA` if
#'            confidence interval not requested.
#'
#'  - `fit`: The original fit object.
#'
#' @param fit The SEM output by [lavaan::lavaan()] or its wrappers.
#' @param x The name of the focal variable in the model, the variable with
#'          its effect on the outcome variable being moderated.
#' @param y The name of the outcome variable (dependent variable) in the model.
#' @param w The name of the moderator in the model.
#' @param x_w The name of the product term (x * w) in the model. It can be
#'            the variable generated by the colon operator, e.g., `"x:w"`,
#'            which is only in the model and not in the original data set.
#' @param standardized_x If `TRUE`, the default, `x` is standardized when
#' computing the standardized moderation effect.
#' @param standardized_y If `TRUE`, the default, `y` is standardized when
#' computing the standardized moderation effect.
#' @param standardized_w If `TRUE`, the default, `w` is standardized when
#' computing the standardized moderation effect.
#' @param boot_ci Boolean. Whether nonparametric bootstrapping will be
#'                conducted. Default is `FALSE`.
#' @param boot_out If set to the output of [manymome::do_boot()], the stored
#' bootstrap estimates will be retrieved to form the bootstrap confidence
#' interval. If set, bootstrap estimates stored in `fit`, if any, will not
#' be used. Default is `NULL`.
#' @param R (Not used in the current version. Used when `use_old_version`
#' is set to `TRUE`.)
#' The number of nonparametric bootstrapping samples. Default is 100.
#' Set this to at least 2000 in actual use.
#' @param conf The level of confidence. Default is .95, i.e., 95%.
#' @param use_old_version If set to `TRUE`, it will use the bootstrapping
#' method used in 0.2.7.4 or before. Included only for reproducing
#' previous results if necessary. Default is `FALSE`.
#' @param ... (Not used in the current version. Used when `use_old_version`
#' is set to `TRUE`.) Optional arguments to be passed to [boot::boot()]. Parallel
#'            processing can be used by adding the appropriate arguments in
#'            [boot::boot()].
#'
#'
#' @author Shu Fai Cheung <https://orcid.org/0000-0002-9871-9448>
#'
#'
#' @references
#' Cheung, S. F., Cheung, S.-H., Lau, E. Y. Y., Hui, C. H., & Vong, W. N.
#' (2022) Improving an old way to measure moderation effect in standardized
#' units. *Health Psychology*, *41*(7), 502-505.
#' \doi{10.1037/hea0001188}
#'
#' @examples
#'
#' #Load a test data of 500 cases
#'
#' dat <- test_mod1
#' library(lavaan)
#' mod <-
#' "
#' med ~ iv + mod + iv:mod + cov1
#' dv ~ med + cov2
#' "
#'
#' fit <- sem(mod, dat)
#'
#' # Compute the standardized moderation effect
#' out_noboot <- stdmod_lavaan(fit = fit,
#'                             x = "iv",
#'                             y = "med",
#'                             w = "mod",
#'                             x_w = "iv:mod")
#' out_noboot
#'
#' # Compute the standardized moderation effect and
#' # its percentile confidence interval using
#' # nonparametric bootstrapping
#' # Fit the model with bootstrap confidence intervals
#' # At least 2000 bootstrap samples should be used
#' # in real research. 50 is used here only for
#' # illustration.
#' fit <- sem(mod, dat, se = "boot", bootstrap = 50,
#'            iseed = 89574)
#' out_boot <- stdmod_lavaan(fit = fit,
#'                           x = "iv",
#'                           y = "med",
#'                           w = "mod",
#'                           x_w = "iv:mod",
#'                           boot_ci = TRUE)
#' out_boot
#'
#' @export

stdmod_lavaan <- function(fit,
                          x,
                          y,
                          w,
                          x_w,
                          standardized_x = TRUE,
                          standardized_y = TRUE,
                          standardized_w = TRUE,
                          boot_ci = FALSE,
                          boot_out = NULL,
                          R = 100,
                          conf = 0.95,
                          use_old_version = FALSE,
                          ...) {
    if (use_old_version) {
        out <- (stdmod_lavaan_old(
                  fit = fit,
                  x = x,
                  y = y,
                  w = w,
                  x_w = x_w,
                  boot_ci = boot_ci,
                  R = R,
                  conf = conf,
                  ...
                ))
      } else {
        out <- (stdmod_lavaan2(
                  fit = fit,
                  x = x,
                  y = y,
                  w = w,
                  x_w = x_w,
                  standardized_x = standardized_x,
                  standardized_y = standardized_y,
                  standardized_w = standardized_w,
                  boot_ci = boot_ci,
                  conf = conf,
                  boot_out = boot_out
                ))
      }
    out$call <- match.call()
    out$old_version <- use_old_version
    return(out)
  }

#' @noRd
stdmod_lavaan_old <- function(fit,
                          x,
                          y,
                          w,
                          x_w,
                          boot_ci = FALSE,
                          R = 100,
                          conf = 0.95, ...) {
    if (!requireNamespace("lavaan", quietly = TRUE)) {
        stop(paste("lavaan needs to be installed to run this function."))
      }
    boot_i <- boot_i_gen(fit = fit, x = x,
                                    y = y,
                                    w = w,
                                    x_w = x_w)
    dat_org <- lavaan::lavInspect(fit, "data")
    if (!boot_ci) {
        stdmod <- stdmod_from_fit(fit = fit,
                                  x = x,
                                  y = y,
                                  w = w,
                                  x_w = x_w)
        boot_out <- NA
        stdmod_ci <- NA
      } else {
        boot_out <- boot::boot(dat_org, boot_i, R = R, ...)
        stdmod <- boot_out$t0
        if (any(is.na(boot_out$t))) {
            tmp <- paste0("Fit not safe in at least one bootstrap ",
                          "samples, either with an error or with ",
                          "a warning message. Please check the original ",
                          "fit object, or set 'warn' to FALSE ",
                          "in lavaan if the warning can be ignored.")
            warning(tmp)
            cat(paste0("\nNumber of valid bootstrap estimates: ",
                           sum(!is.na(boot_out$t)), "\n"))
          }
        stdmod_ci <- boot::boot.ci(boot_out,
                                   type = "perc",
                                   conf = conf)$percent[4:5]
        names(stdmod_ci) <- c("lower_bound", "upper_bound")
      }
    out <- list(stdmod = stdmod,
                ci     = stdmod_ci,
                boot_out = boot_out)
    out$call <- match.call()
    out$fit <- fit
    class(out) <- c("stdmod_lavaan", class(out))
    return(out)
  }

#' @noRd

boot_i_gen <- function(fit, x, y, w, x_w) {
  force(fit)
  force(x)
  force(y)
  force(w)
  force(x_w)
  function(dat, i = NULL) {
      if (is.null(i)) {
          dat_i <- dat
        } else {
          dat_i <- dat[i, ]
        }
      fit_test <- tryCatch(fit_i <- lavaan::update(fit,
                                       data = dat_i,
                                       se = "none",
                                       h1 = FALSE,
                                       baseline = FALSE,
                                       test = "standard"
                                       ),
                            error = function(e) e,
                            warning = function(e) e)
      if (inherits(fit_test, "error")) {
          print(fit_test)
          return(NA)
        }
      if (inherits(fit_test, "warning")) {
          # options_old <- options("warn" = -1)
          warning(fit_test)
          # options("warn" = options_old$warn)
          return(NA)
        }
      if (!inherits(fit_i, "lavaan")) {
          warning("Something's wrong. Fit failed.")
          return(NA)
        } else {
          x_w_std_i <- stdmod_from_fit(fit = fit_i,
                                       x = x,
                                       y = y,
                                       w = w,
                                       x_w = x_w)
          return(x_w_std_i)
        }
    }
  }

#' @noRd

stdmod_from_fit <- function(fit, x, y, w, x_w) {
    fit_cov_implied <- lavaan::lavInspect(fit, "implied")
    x_sd <- sqrt(diag(fit_cov_implied$cov)[x])
    w_sd <- sqrt(diag(fit_cov_implied$cov)[w])
    y_sd <- sqrt(diag(fit_cov_implied$cov)[y])
    x_w_i <- lavaan::coef(fit)[paste0(y, "~", x_w)]
    x_w_std_i <- x_w_i * x_sd * w_sd / y_sd
    return(x_w_std_i)
  }

#' @noRd

stdmod_lavaan2 <- function(fit,
                           x = NULL,
                           y = NULL,
                           w = NULL,
                           x_w = NULL,
                           standardized_x = TRUE,
                           standardized_y = TRUE,
                           standardized_w = TRUE,
                           boot_ci = FALSE,
                           conf = 0.95,
                           boot_out = NULL) {
    if (lavaan::lavInspect(fit, "ngroups") > 1) {
        stop("Multiple-group models are not supported.")
      }
    tmp <- try(lavaan::lavInspect(fit, "boot"),
               silent = TRUE)
    if (boot_ci && inherits(tmp, "try-error") &&
        is.null(boot_out)) {
        stop("Bootstrap confidence interval ",
             "requested but either bootstrapping was ",
             "not done when fitting the model or ",
             "boot_out was not set.")
      }
    stdmod0 <- stdmod_from_fit2(fit = fit,
                               x = x,
                               y = y,
                               w = w,
                               x_w = x_w,
                               standardized_x = standardized_x,
                               standardized_y = standardized_y,
                               standardized_w = standardized_w)
    if (boot_ci) {
        if (is.null(boot_out)) {
            boot_out <- manymome::fit2boot_out(fit)
          }
        est_all <- lapply(boot_out, function(x) x$est)
        implied_cov_all <- lapply(boot_out, function(x) x$implied$cov)
        stdmod_all <- mapply(stdmod_from_fit2,
                             est = est_all,
                             implied_cov = implied_cov_all,
                             MoreArgs = list(x = x,
                                             y = y,
                                             w = w,
                                             x_w = x_w,
                                             standardized_x = standardized_x,
                                             standardized_y = standardized_y,
                                             standardized_w = standardized_w))
        boot_tmp <- list(t0 = stdmod0,
                         t = matrix(stdmod_all, ncol = 1),
                         R = length(stdmod_all))
        stdmod_ci <- boot::boot.ci(boot_tmp,
                            type = "perc",
                            conf = conf)$percent[4:5]
        names(stdmod_ci) <- c("lower_bound", "upper_bound")
      } else {
        stdmod_ci <- NA
        boot_tmp <- NA
      }
    out <- list(stdmod = stdmod0,
                ci     = stdmod_ci,
                boot_out = boot_tmp)
    out$call <- match.call()
    out$fit <- fit
    class(out) <- c("stdmod_lavaan", class(out))
    return(out)
  }

#' @noRd

stdmod_from_fit2 <- function(fit,
                             est,
                             implied_cov,
                             x,
                             y,
                             w,
                             x_w,
                             standardized_x = TRUE,
                             standardized_y = TRUE,
                             standardized_w = TRUE,
                             full = FALSE
                             ) {
    if (!missing(fit)) {
        est <- lavaan::parameterEstimates(fit)
        implied_cov <- lavaan::lavInspect(fit, "implied")$cov
      }
    sd_x <- ifelse(standardized_x,
                   sqrt(implied_cov[x, x]),
                   1)
    sd_y <- ifelse(standardized_y,
                   sqrt(implied_cov[y, y]),
                   1)
    sd_w <- ifelse(standardized_w,
                   sqrt(implied_cov[w, w]),
                   1)
    est_x_w <- est[(est$lhs == y) &
                   (est$rhs == x_w) &
                   (est$op == "~"), "est"]
    out <- est_x_w * sd_x * sd_w / sd_y
    if (full) {
        return(c(out, est_x_w, sd_x, sd_w, sd_y))
      } else {
        return(out)
      }
  }