R/cTMed-mc-indirect-central.R

In cTMed: Continuous-Time Mediation

#' Monte Carlo Sampling Distribution
#' of Indirect Effect Centrality
#' Over a Specific Time Interval
#' or a Range of Time Intervals
#'
#' This function generates a Monte Carlo method
#' sampling distribution
#' of the indirect effect centrality
#' at a particular time interval \eqn{\Delta t}
#' using the first-order stochastic differential equation model
#' drift matrix \eqn{\boldsymbol{\Phi}}.
#'
#' @details See [IndirectCentral()] for more details.
#'
#' ## Monte Carlo Method
#'   Let \eqn{\boldsymbol{\theta}} be
#'   \eqn{\mathrm{vec} \left( \boldsymbol{\Phi} \right)},
#'   that is,
#'   the elements of the \eqn{\boldsymbol{\Phi}} matrix
#'   in vector form sorted column-wise.
#'   Let \eqn{\hat{\boldsymbol{\theta}}} be
#'   \eqn{\mathrm{vec} \left( \hat{\boldsymbol{\Phi}} \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) .
#'   }
#'   Let
#'   \eqn{\mathbf{g} \left( \hat{\boldsymbol{\theta}} \right)}
#'   be a parameter that is a function of the estimated parameters.
#'   A sampling distribution of
#'   \eqn{\mathbf{g} \left( \hat{\boldsymbol{\theta}} \right)} ,
#'   which we refer to as
#'   \eqn{\mathbf{g} \left( \hat{\boldsymbol{\theta}}^{\ast} \right)} ,
#'   can be generated by using the simulated estimates
#'   to calculate
#'   \eqn{\mathbf{g}}.
#'   The standard deviations of the simulated estimates
#'   are the standard errors.
#'   Percentiles corresponding to
#'   \eqn{100 \left( 1 - \alpha \right) \%}
#'   are the confidence intervals.
#'
#' @author Ivan Jacob Agaloos Pesigan
#'
#' @inheritParams Indirect
#' @inheritParams MCPhi
#' @inheritParams MCMed
#' @inherit Indirect references
#'
#' @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 ("MCIndirectCentral").}
#'     \item{output}{A list the length of which is equal to
#'         the length of `delta_t`.}
#'   }
#'   Each element in the `output` list has the following elements:
#'   \describe{
#'     \item{est}{A vector of indirect effect centrality.}
#'     \item{thetahatstar}{A matrix of Monte Carlo
#'     indirect effect centrality.}
#'   }
#'
#' @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")
#' vcov_phi_vec <- matrix(
#'   data = c(
#'     0.00843, 0.00040, -0.00151,
#'     -0.00600, -0.00033, 0.00110,
#'     0.00324, 0.00020, -0.00061,
#'     0.00040, 0.00374, 0.00016,
#'     -0.00022, -0.00273, -0.00016,
#'     0.00009, 0.00150, 0.00012,
#'     -0.00151, 0.00016, 0.00389,
#'     0.00103, -0.00007, -0.00283,
#'     -0.00050, 0.00000, 0.00156,
#'     -0.00600, -0.00022, 0.00103,
#'     0.00644, 0.00031, -0.00119,
#'     -0.00374, -0.00021, 0.00070,
#'     -0.00033, -0.00273, -0.00007,
#'     0.00031, 0.00287, 0.00013,
#'     -0.00014, -0.00170, -0.00012,
#'     0.00110, -0.00016, -0.00283,
#'     -0.00119, 0.00013, 0.00297,
#'     0.00063, -0.00004, -0.00177,
#'     0.00324, 0.00009, -0.00050,
#'     -0.00374, -0.00014, 0.00063,
#'     0.00495, 0.00024, -0.00093,
#'     0.00020, 0.00150, 0.00000,
#'     -0.00021, -0.00170, -0.00004,
#'     0.00024, 0.00214, 0.00012,
#'     -0.00061, 0.00012, 0.00156,
#'     0.00070, -0.00012, -0.00177,
#'     -0.00093, 0.00012, 0.00223
#'   ),
#'   nrow = 9
#' )
#'
#' # Specific time interval ----------------------------------------------------
#' MCIndirectCentral(
#'   phi = phi,
#'   vcov_phi_vec = vcov_phi_vec,
#'   delta_t = 1,
#'   R = 100L # use a large value for R in actual research
#' )
#'
#' # Range of time intervals ---------------------------------------------------
#' mc <- MCIndirectCentral(
#'   phi = phi,
#'   vcov_phi_vec = vcov_phi_vec,
#'   delta_t = 1:5,
#'   R = 100L # use a large value for R in actual research
#' )
#' plot(mc)
#'
#' # Methods -------------------------------------------------------------------
#' # MCIndirectCentral has a number of methods including
#' # print, summary, confint, and plot
#' print(mc)
#' summary(mc)
#' confint(mc, level = 0.95)
#' plot(mc)
#'
#' @family Continuous-Time Mediation Functions
#' @keywords cTMed network mc
#' @export
MCIndirectCentral <- function(phi,
                              vcov_phi_vec,
                              delta_t,
                              R,
                              test_phi = TRUE,
                              ncores = NULL,
                              seed = NULL,
                              tol = 0.01) {
  idx <- rownames(phi)
  stopifnot(
    idx == colnames(phi)
  )
  total <- FALSE
  args <- list(
    phi = phi,
    vcov_phi_vec = vcov_phi_vec,
    delta_t = delta_t,
    R = R,
    test_phi = test_phi,
    ncores = ncores,
    seed = seed,
    method = "mc",
    network = TRUE,
    total = total
  )
  delta_t <- sort(
    ifelse(
      test = delta_t < tol,
      yes = tol, # .Machine$double.xmin
      no = delta_t
    )
  )
  output <- .MCCentral(
    phi = phi,
    vcov_phi_vec = vcov_phi_vec,
    delta_t = delta_t,
    total = total,
    R = R,
    test_phi = test_phi,
    ncores = ncores,
    seed = seed
  )
  names(output) <- delta_t
  out <- list(
    call = match.call(),
    args = args,
    fun = "MCIndirectCentral",
    output = output
  )
  class(out) <- c(
    "ctmedmc",
    class(out)
  )
  out
}