R/cTMed-mc-beta-std.R

In cTMed: Continuous-Time Mediation

#' Monte Carlo Sampling Distribution
#' for the Elements of the Standardized Matrix of Lagged Coefficients
#' Over a Specific Time Interval
#' or a Range of Time Intervals
#'
#' This function generates a Monte Carlo method
#' sampling distribution
#' for the elements of the standardized matrix of lagged coefficients
#' \eqn{\boldsymbol{\beta}}
#' over a specific time interval \eqn{\Delta t}
#' or a range of time intervals
#' using the first-order stochastic differential equation model
#' drift matrix \eqn{\boldsymbol{\Phi}}
#' and process noise covariance matrix \eqn{\boldsymbol{\Sigma}}.
#'
#' @details See [TotalStd()].
#'
#' ## 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) .
#'   }
#'   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 MCMedStd
#' @inherit IndirectStd 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 ("MCBetaStd").}
#'     \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}{Estimated elements of the standardized matrix
#'       of lagged coefficients.}
#'     \item{thetahatstar}{A matrix of Monte Carlo
#'       elements of the standardized matrix of lagged coefficients.}
#'   }
#'
#' @examples
#' 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
#' )
#' vcov_theta <- matrix(
#'   data = c(
#'     0.00843, 0.00040, -0.00151, -0.00600, -0.00033,
#'     0.00110, 0.00324, 0.00020, -0.00061, -0.00115,
#'     0.00011, 0.00015, 0.00001, -0.00002, -0.00001,
#'     0.00040, 0.00374, 0.00016, -0.00022, -0.00273,
#'     -0.00016, 0.00009, 0.00150, 0.00012, -0.00010,
#'     -0.00026, 0.00002, 0.00012, 0.00004, -0.00001,
#'     -0.00151, 0.00016, 0.00389, 0.00103, -0.00007,
#'     -0.00283, -0.00050, 0.00000, 0.00156, 0.00021,
#'     -0.00005, -0.00031, 0.00001, 0.00007, 0.00006,
#'     -0.00600, -0.00022, 0.00103, 0.00644, 0.00031,
#'     -0.00119, -0.00374, -0.00021, 0.00070, 0.00064,
#'     -0.00015, -0.00005, 0.00000, 0.00003, -0.00001,
#'     -0.00033, -0.00273, -0.00007, 0.00031, 0.00287,
#'     0.00013, -0.00014, -0.00170, -0.00012, 0.00006,
#'     0.00014, -0.00001, -0.00015, 0.00000, 0.00001,
#'     0.00110, -0.00016, -0.00283, -0.00119, 0.00013,
#'     0.00297, 0.00063, -0.00004, -0.00177, -0.00013,
#'     0.00005, 0.00017, -0.00002, -0.00008, 0.00001,
#'     0.00324, 0.00009, -0.00050, -0.00374, -0.00014,
#'     0.00063, 0.00495, 0.00024, -0.00093, -0.00020,
#'     0.00006, -0.00010, 0.00000, -0.00001, 0.00004,
#'     0.00020, 0.00150, 0.00000, -0.00021, -0.00170,
#'     -0.00004, 0.00024, 0.00214, 0.00012, -0.00002,
#'     -0.00004, 0.00000, 0.00006, -0.00005, -0.00001,
#'     -0.00061, 0.00012, 0.00156, 0.00070, -0.00012,
#'     -0.00177, -0.00093, 0.00012, 0.00223, 0.00004,
#'     -0.00002, -0.00003, 0.00001, 0.00003, -0.00013,
#'     -0.00115, -0.00010, 0.00021, 0.00064, 0.00006,
#'     -0.00013, -0.00020, -0.00002, 0.00004, 0.00057,
#'     0.00001, -0.00009, 0.00000, 0.00000, 0.00001,
#'     0.00011, -0.00026, -0.00005, -0.00015, 0.00014,
#'     0.00005, 0.00006, -0.00004, -0.00002, 0.00001,
#'     0.00012, 0.00001, 0.00000, -0.00002, 0.00000,
#'     0.00015, 0.00002, -0.00031, -0.00005, -0.00001,
#'     0.00017, -0.00010, 0.00000, -0.00003, -0.00009,
#'     0.00001, 0.00014, 0.00000, 0.00000, -0.00005,
#'     0.00001, 0.00012, 0.00001, 0.00000, -0.00015,
#'     -0.00002, 0.00000, 0.00006, 0.00001, 0.00000,
#'     0.00000, 0.00000, 0.00010, 0.00001, 0.00000,
#'     -0.00002, 0.00004, 0.00007, 0.00003, 0.00000,
#'     -0.00008, -0.00001, -0.00005, 0.00003, 0.00000,
#'     -0.00002, 0.00000, 0.00001, 0.00005, 0.00001,
#'     -0.00001, -0.00001, 0.00006, -0.00001, 0.00001,
#'     0.00001, 0.00004, -0.00001, -0.00013, 0.00001,
#'     0.00000, -0.00005, 0.00000, 0.00001, 0.00012
#'   ),
#'   nrow = 15
#' )
#'
#' # Specific time interval ----------------------------------------------------
#' MCBetaStd(
#'   phi = phi,
#'   sigma = sigma,
#'   vcov_theta = vcov_theta,
#'   delta_t = 1,
#'   R = 100L # use a large value for R in actual research
#' )
#'
#' # Range of time intervals ---------------------------------------------------
#' mc <- MCBetaStd(
#'   phi = phi,
#'   sigma = sigma,
#'   vcov_theta = vcov_theta,
#'   delta_t = 1:5,
#'   R = 100L # use a large value for R in actual research
#' )
#' plot(mc)
#'
#' # Methods -------------------------------------------------------------------
#' # MCBetaStd 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 beta mc
#' @export
MCBetaStd <- function(phi,
                      sigma,
                      vcov_theta,
                      delta_t,
                      R,
                      test_phi = TRUE,
                      ncores = NULL,
                      seed = NULL,
                      tol = 0.01) {
  idx <- rownames(phi)
  stopifnot(
    idx == colnames(phi)
  )
  args <- list(
    phi = phi,
    sigma = sigma,
    vcov_theta = vcov_theta,
    delta_t = delta_t,
    R = R,
    test_phi = test_phi,
    ncores = ncores,
    seed = seed,
    method = "mc",
    network = FALSE
  )
  delta_t <- sort(
    ifelse(
      test = delta_t < tol,
      yes = tol, # .Machine$double.xmin
      no = delta_t
    )
  )
  output <- .MCBetaStd(
    phi = phi,
    sigma = sigma,
    vcov_theta = vcov_theta,
    delta_t = delta_t,
    R = R,
    test_phi = test_phi,
    ncores = ncores,
    seed = seed
  )
  names(output) <- delta_t
  out <- list(
    call = match.call(),
    args = args,
    fun = "MCBetaStd",
    output = output
  )
  class(out) <- c(
    "ctmedmc",
    class(out)
  )
  out
}