R/mbco.R

In quantPsych/RMediation: Mediation Analysis Confidence Intervals

#' Model-based Constrained Optimization (MBCO) Chi-squared Test
#'
#' This function computes asymptotic MBCO chi-squared test for a smooth function of model parameters including a function of indirect effects.
#'
#' @param h0 An \code{OpenMx} model estimated under a null hypothesis, which is a more constrained model
#' @param h1 An \code{OpenMx} model estimated under an alternative hypothesis, which is a less constrained model. This is usually a model hypothesized by a researcher.
#' @param R The number of bootstrap draws.
#' @param type If 'asymp' (default), the asymptotic MBCO chi-squares test comparing fit of h0 and h1. If 'parametric', the parametric bootstrap MBCO chi-squared test is computed. If 'semi', the semi-parametric MBCO chi-squared is computed.
#' @param alpha Significance level with the default value of .05
#' @param checkHess If 'No' (default), the Hessian matrix would not be calculated.
#' @param checkSE if 'No' (default), the standard errors would not be calculated.
#' @param optim Choose optimizer available in OpenMx. The default optimizer is "SLSQP". Other optimizer choices are available. See \link{mxOption} for more details.
#' @param precision Functional precision. The default value is set to 1e-9. See \link{mxOption} for more details.
#' @return A \link{list} that contains \item{chisq}{asymptotic chi-squared test statistic value} \item{\code{df}}{chi-squared df} \item{p}{chi-squared p-value computed based on the method specified by the argument \code{type}}
#' @author Davood Tofighi \email{dtofighi@@gmail.com}
#' @import OpenMx
#' @importFrom modelr resample_bootstrap
#' @export
#' @references  Tofighi, D., & Kelley, K. (2020). Indirect effects in sequential mediation models: Evaluating methods for hypothesis testing and confidence interval formation. \emph{Multivariate Behavioral Research}, \bold{55}, 188–210. \doi{https://doi.org/10.1080/00273171.2019.1618545}
#'
#'  Tofighi, D. (2020). Bootstrap Model-Based Constrained Optimization Tests of Indirect Effects. \emph{Frontiers in Psychology}, \bold{10}, 2989. \doi{https://doi.org/10.3389/fpsyg.2019.02989}

#' @examples
#' library(OpenMx)
#' data(memory_exp)
#' memory_exp$x <- as.numeric(memory_exp$x) - 1 # manually creating dummy codes
#' endVar <- c("x", "repetition", "imagery", "recall")
#' manifests <- c("x", "repetition", "imagery", "recall")
#' full_model <- mxModel(
#'   "memory_example",
#'   type = "RAM",
#'   manifestVars = manifests,
#'   mxPath(
#'     from = "x",
#'     to = endVar,
#'     arrows = 1,
#'     free = TRUE,
#'     values = .2,
#'     labels = c("a1", "a2", "cp")
#'   ),
#'   mxPath(
#'     from = "repetition",
#'     to = "recall",
#'     arrows = 1,
#'     free = TRUE,
#'     values = .2,
#'     labels = "b1"
#'   ),
#'   mxPath(
#'     from = "imagery",
#'     to = "recall",
#'     arrows = 1,
#'     free = TRUE,
#'     values = .2,
#'     labels = "b2"
#'   ),
#'   mxPath(
#'     from = manifests,
#'     arrows = 2,
#'     free = TRUE,
#'     values = .8
#'   ),
#'   mxPath(
#'     from = "one",
#'     to = endVar,
#'     arrows = 1,
#'     free = TRUE,
#'     values = .1
#'   ),
#'   mxAlgebra(a1 * b1, name = "ind1"),
#'   mxAlgebra(a2 * b2, name = "ind2"),
#'   mxCI("ind1", type = "both"),
#'   mxCI("ind2", type = "both"),
#'   mxData(observed = memory_exp, type = "raw")
#' )
#' ## Reduced  Model for indirect effect: a1*b1
#' null_model1 <- mxModel(
#'   model = full_model,
#'   name = "Null Model 1",
#'   mxConstraint(ind1 == 0, name = "ind1_eq0_constr")
#' )
#' full_model <- mxTryHard(full_model, checkHess = FALSE, silent = TRUE)
#' null_model1 <- mxTryHard(null_model1, checkHess = FALSE, silent = TRUE)
#' mbco(null_model1, full_model)
mbco <- function(h0 = NULL,
                 h1 = NULL,
                 R = 10L,
                 type = "asymp",
                 alpha = .05,
                 checkHess = "No",
                 checkSE = "No",
                 optim = "SLSQP",
                 precision = 1e-9) {
  if (missing(h0)) {
    stop("'h0' argument be a MxModel object")
  }

  if (missing(h1)) {
    stop("'h1' argument be a MxModel object")
  }

  if (!all(sapply(c(h0, h1), is, "MxModel"))) {
    stop("The 'h0' and 'h1' argument must be MxModel objects")
  }

  type <-
    match.arg(type, c("asymp", "parametric", "semi")) # checks if one of the types of partially matches
  optim <- match.arg(optim, c("NPSOL", "CSOLNP", "SLSQP"))

  if (!OpenMx::imxHasNPSOL()) {
    optim <- "SLSQP"
  } # if NPSOL is not available, use SLSQP

  res <-
    if (type == "asymp") {
      # Asymptotic MBCO LRT
      mbco_asymp(h0 = h0, h1 = h1, alpha = alpha)
    } else if (type == "parametric") {
      # Parametric bootstrap MBCO LRT
      mbco_parametric(
        h0 = h0,
        h1 = h1,
        R = R,
        alpha = alpha,
        checkHess = checkHess,
        checkSE = checkSE,
        optim = optim,
        precision = precision
      )
    } else if (type == "semi") {
      # Semiparametric bootstrap MBCO LRT
      mbco_semi(
        h0 = h0,
        h1 = h1,
        R = R,
        alpha = alpha,
        checkHess = checkHess,
        checkSE = checkSE,
        optim = optim,
        precision = precision
      )
    }

  return(res)
}