R/cTMed-mc-phi-sigma.R

In cTMed: Continuous-Time Mediation

#' Generate Random Drift Matrices
#' and Process Noise Covariance Matrices
#' Using the Monte Carlo Method
#'
#' This function generates random
#' drift matrices \eqn{\boldsymbol{\Phi}}
#' and process noise covariabces matrices \eqn{\boldsymbol{\Sigma}}
#' using the Monte Carlo method.
#'
#' @details
#' ## Monte Carlo Method
#'   Let \eqn{\boldsymbol{\theta}} be
#'   a vector that combines
#'   \eqn{\mathrm{vec} \left( \boldsymbol{\Phi} \right)},
#'   that is,
#'   the elements of the \eqn{\boldsymbol{\Phi}} matrix
#'   in vector form sorted column-wise and
#'   \eqn{\mathrm{vech} \left( \boldsymbol{\Sigma} \right)},
#'   that is,
#'   the unique elements of the \eqn{\boldsymbol{\Sigma}} matrix
#'   in vector form sorted column-wise.
#'   Let \eqn{\hat{\boldsymbol{\theta}}} be
#'   a vector that combines
#'   \eqn{\mathrm{vec} \left( \hat{\boldsymbol{\Phi}} \right)} and
#'   \eqn{\mathrm{vech} \left( \hat{\boldsymbol{\Sigma}} \right)}.
#'   Based on the asymptotic properties of maximum likelihood estimators,
#'   we can assume that estimators are normally distributed
#'   around the population parameters.
#'   \deqn{
#'   	\hat{\boldsymbol{\theta}}
#'   	\sim
#'   	\mathcal{N}
#'   	\left(
#'   	\boldsymbol{\theta},
#'   	\mathbb{V} \left( \hat{\boldsymbol{\theta}} \right)
#'   	\right)
#'   }
#'   Using this distributional assumption,
#'   a sampling distribution of \eqn{\hat{\boldsymbol{\theta}}}
#'   which we refer to as \eqn{\hat{\boldsymbol{\theta}}^{\ast}}
#'   can be generated by replacing the population parameters
#'   with sample estimates,
#'   that is,
#'   \deqn{
#'   	\hat{\boldsymbol{\theta}}^{\ast}
#'   	\sim
#'   	\mathcal{N}
#'   	\left(
#'   	\hat{\boldsymbol{\theta}},
#'   	\hat{\mathbb{V}} \left( \hat{\boldsymbol{\theta}} \right)
#'   	\right) .
#'   }
#'
#' @author Ivan Jacob Agaloos Pesigan
#'
#' @inheritParams IndirectStd
#' @inheritParams MCPhi
#' @param vcov_theta Numeric matrix.
#'   The sampling variance-covariance matrix of
#'   \eqn{\mathrm{vec} \left( \boldsymbol{\Phi} \right)} and
#'   \eqn{\mathrm{vech} \left( \boldsymbol{\Sigma} \right)}
#'
#' @return Returns an object
#'   of class `ctmedmc` which is a list with the following elements:
#'   \describe{
#'     \item{call}{Function call.}
#'     \item{args}{Function arguments.}
#'     \item{fun}{Function used ("MCPhiSigma").}
#'     \item{output}{A list simulated drift matrices.}
#'   }
#'
#' @examples
#' set.seed(42)
#' phi <- matrix(
#'   data = c(
#'     -0.357, 0.771, -0.450,
#'     0.0, -0.511, 0.729,
#'     0, 0, -0.693
#'   ),
#'   nrow = 3
#' )
#' colnames(phi) <- rownames(phi) <- c("x", "m", "y")
#' sigma <- matrix(
#'   data = c(
#'     0.24455556, 0.02201587, -0.05004762,
#'     0.02201587, 0.07067800, 0.01539456,
#'     -0.05004762, 0.01539456, 0.07553061
#'   ),
#'   nrow = 3
#' )
#' MCPhiSigma(
#'   phi = phi,
#'   sigma = sigma,
#'   vcov_theta = 0.1 * diag(15),
#'   R = 100L # use a large value for R in actual research
#' )
#'
#' @family Continuous-Time Mediation Functions
#' @keywords cTMed mc
#' @export
MCPhiSigma <- function(phi,
                       sigma,
                       vcov_theta,
                       R,
                       test_phi = TRUE,
                       ncores = NULL,
                       seed = NULL) {
  idx <- rownames(phi)
  stopifnot(
    idx == colnames(phi)
  )
  args <- list(
    phi = phi,
    sigma = sigma,
    vcov_theta = vcov_theta,
    R = R,
    test_phi = test_phi,
    ncores = ncores,
    seed = seed,
    method = "mc"
  )
  # nocov start
  par <- FALSE
  if (!is.null(ncores)) {
    ncores <- as.integer(ncores)
    if (ncores > R) {
      ncores <- R
    }
    if (ncores > 1) {
      par <- TRUE
    }
  }
  if (par) {
    os_type <- Sys.info()["sysname"]
    if (os_type == "Darwin") {
      fork <- TRUE
    } else if (os_type == "Linux") {
      fork <- TRUE
    } else {
      fork <- FALSE
    }
    if (fork) {
      if (!is.null(seed)) {
        set.seed(seed)
      }
      output <- parallel::mclapply(
        X = seq_len(R),
        FUN = function(i) {
          .MCPhiSigmaI(
            theta = c(.Vec(phi), .Vech(sigma)),
            vcov_theta = vcov_theta,
            test_phi = test_phi
          )
        },
        mc.cores = ncores
      )
    } else {
      cl <- parallel::makeCluster(ncores)
      on.exit(
        parallel::stopCluster(cl = cl)
      )
      if (!is.null(seed)) {
        parallel::clusterSetRNGStream(
          cl = cl,
          iseed = seed
        )
      }
      output <- parallel::parLapply(
        cl = cl,
        X = seq_len(R),
        fun = function(i) {
          .MCPhiSigmaI(
            theta = c(.Vec(phi), .Vech(sigma)),
            vcov_theta = vcov_theta,
            test_phi = test_phi
          )
        }
      )
    }
    # nocov end
  } else {
    if (!is.null(seed)) {
      set.seed(seed)
    }
    output <- .MCPhiSigma(
      theta = c(.Vec(phi), .Vech(sigma)),
      vcov_theta = vcov_theta,
      R = R,
      test_phi = test_phi
    )
  }
  output <- lapply(
    X = output,
    FUN = function(x) {
      colnames(x[[1]]) <- rownames(x[[1]]) <- idx
      x
    }
  )
  out <- list(
    call = match.call(),
    args = args,
    fun = "MCPhiSigma",
    output = output
  )
  class(out) <- c(
    "ctmedmcphi",
    class(out)
  )
  out
}