Skip to contents

This function generates a Monte Carlo method sampling distribution of the standardized direct effect centrality at a particular time interval \(\Delta t\) using the first-order stochastic differential equation model drift matrix \(\boldsymbol{\Phi}\) and process noise covariance matrix \(\boldsymbol{\Sigma}\).

Usage

MCDirectCentralStd(
  phi,
  sigma,
  vcov_theta,
  delta_t,
  R,
  test_phi = TRUE,
  sigma_diag = FALSE,
  ncores = NULL,
  seed = NULL,
  tol = 0.001
)

Arguments

phi

Numeric matrix. The drift matrix (\(\boldsymbol{\Phi}\)). phi should have row and column names pertaining to the variables in the system.

sigma

Numeric matrix. The process noise covariance matrix (\(\boldsymbol{\Sigma}\)).

vcov_theta

Numeric matrix. The sampling variance-covariance matrix of \(\mathrm{vec} \left( \boldsymbol{\Phi} \right)\) and \(\mathrm{vech} \left( \boldsymbol{\Sigma} \right)\)

delta_t

Numeric. Time interval (\(\Delta t\)).

R

Positive integer. Number of replications.

test_phi

Logical. If test_phi = TRUE, the function tests the stability of the generated drift matrix \(\boldsymbol{\Phi}\). If the test returns FALSE, the function generates a new drift matrix \(\boldsymbol{\Phi}\) and runs the test recursively until the test returns TRUE.

sigma_diag

Logical. If sigma_diag = TRUE, treat \(\boldsymbol{\Sigma}\) as a diagonal matrix.

ncores

Positive integer. Number of cores to use. If ncores = NULL, use a single core. Consider using multiple cores when number of replications R is a large value.

seed

Random seed.

tol

Numeric. Smallest possible time interval to allow.

Value

Returns an object of class ctmedmc which is a list with the following elements:

call

Function call.

args

Function arguments.

fun

Function used ("MCDirectCentralStd").

output

A list of length length(delta_t).

Each element in the output list has the following elements:

est

A vector of standardized direct effect centrality.

thetahatstar

A matrix of Monte Carlo standardized direct effect centrality.

Details

See DirectCentralStd() for more details.

Monte Carlo Method

Let \(\boldsymbol{\theta}\) be a vector that combines \(\mathrm{vec} \left( \boldsymbol{\Phi} \right)\), that is, the elements of the \(\boldsymbol{\Phi}\) matrix in vector form sorted column-wise and \(\mathrm{vech} \left( \boldsymbol{\Sigma} \right)\), that is, the unique elements of the \(\boldsymbol{\Sigma}\) matrix in vector form sorted column-wise. Let \(\hat{\boldsymbol{\theta}}\) be a vector that combines \(\mathrm{vec} \left( \hat{\boldsymbol{\Phi}} \right)\) and \(\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. $$ \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 \(\hat{\boldsymbol{\theta}}\) which we refer to as \(\hat{\boldsymbol{\theta}}^{\ast}\) can be generated by replacing the population parameters with sample estimates, that is, $$ \hat{\boldsymbol{\theta}}^{\ast} \sim \mathcal{N} \left( \hat{\boldsymbol{\theta}}, \hat{\mathbb{V}} \left( \hat{\boldsymbol{\theta}} \right) \right) . $$ Let \(\mathbf{g} \left( \hat{\boldsymbol{\theta}} \right)\) be a parameter that is a function of the estimated parameters. A sampling distribution of \(\mathbf{g} \left( \hat{\boldsymbol{\theta}} \right)\) , which we refer to as \(\mathbf{g} \left( \hat{\boldsymbol{\theta}}^{\ast} \right)\) , can be generated by using the simulated estimates to calculate \(\mathbf{g}\). The standard deviations of the simulated estimates are the standard errors. Percentiles corresponding to \(100 \left( 1 - \alpha \right) \%\) are the confidence intervals.

References

Bollen, K. A. (1987). Total, direct, and indirect effects in structural equation models. Sociological Methodology, 17, 37. doi:10.2307/271028

Deboeck, P. R., & Preacher, K. J. (2015). No need to be discrete: A method for continuous time mediation analysis. Structural Equation Modeling: A Multidisciplinary Journal, 23 (1), 61-75. doi:10.1080/10705511.2014.973960

Pesigan, I. J. A., Russell, M. A., & Chow, S.-M. (2025). Inferences and effect sizes for direct, indirect, and total effects in continuous-time mediation models. Psychological Methods. doi:10.1037/met0000779

Ryan, O., & Hamaker, E. L. (2021). Time to intervene: A continuous-time approach to network analysis and centrality. Psychometrika, 87 (1), 214-252. doi:10.1007/s11336-021-09767-0

Author

Ivan Jacob Agaloos Pesigan

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 ----------------------------------------------------
MCDirectCentralStd(
  phi = phi,
  sigma = sigma,
  vcov_theta = vcov_theta,
  delta_t = 1,
  R = 100L # use a large value for R in actual research
)
#> Call:
#> MCDirectCentralStd(phi = phi, sigma = sigma, vcov_theta = vcov_theta, 
#>     delta_t = 1, R = 100L)
#> 
#> Direct Effect Centrality
#>   variable interval     est     se   R    2.5%   97.5%
#> 1        x        1  0.5494 0.0373 100  0.4769  0.6287
#> 2        m        1 -0.2858 0.0567 100 -0.3855 -0.1583
#> 3        y        1  0.3888 0.0538 100  0.3004  0.5011

# Range of time intervals ---------------------------------------------------
mc <- MCDirectCentralStd(
  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 -------------------------------------------------------------------
# MCDirectCentralStd has a number of methods including
# print, summary, confint, and plot
print(mc)
#> Call:
#> MCDirectCentralStd(phi = phi, sigma = sigma, vcov_theta = vcov_theta, 
#>     delta_t = 1:5, R = 100L)
#> 
#> Direct Effect Centrality
#>    variable interval     est     se   R    2.5%   97.5%
#> 1         x        1  0.5494 0.0370 100  0.4747  0.6175
#> 2         m        1 -0.2858 0.0552 100 -0.3954 -0.1659
#> 3         y        1  0.3888 0.0608 100  0.2662  0.5115
#> 4         x        2  0.6044 0.0416 100  0.5152  0.6869
#> 5         m        2 -0.3429 0.0708 100 -0.5115 -0.2076
#> 6         y        2  0.5053 0.0738 100  0.3491  0.6443
#> 7         x        3  0.4999 0.0441 100  0.4102  0.5885
#> 8         m        3 -0.3114 0.0741 100 -0.4926 -0.1929
#> 9         y        3  0.4936 0.0798 100  0.3356  0.6446
#> 10        x        4  0.3686 0.0442 100  0.2825  0.4648
#> 11        m        4 -0.2537 0.0725 100 -0.4323 -0.1501
#> 12        y        4  0.4293 0.0856 100  0.2679  0.6023
#> 13        x        5  0.2555 0.0411 100  0.1789  0.3445
#> 14        m        5 -0.1954 0.0681 100 -0.3739 -0.1060
#> 15        y        5  0.3508 0.0896 100  0.1924  0.5406
summary(mc)
#> Call:
#> MCDirectCentralStd(phi = phi, sigma = sigma, vcov_theta = vcov_theta, 
#>     delta_t = 1:5, R = 100L)
#> 
#> Direct Effect Centrality
#>    variable interval     est     se   R    2.5%   97.5%
#> 1         x        1  0.5494 0.0370 100  0.4747  0.6175
#> 2         m        1 -0.2858 0.0552 100 -0.3954 -0.1659
#> 3         y        1  0.3888 0.0608 100  0.2662  0.5115
#> 4         x        2  0.6044 0.0416 100  0.5152  0.6869
#> 5         m        2 -0.3429 0.0708 100 -0.5115 -0.2076
#> 6         y        2  0.5053 0.0738 100  0.3491  0.6443
#> 7         x        3  0.4999 0.0441 100  0.4102  0.5885
#> 8         m        3 -0.3114 0.0741 100 -0.4926 -0.1929
#> 9         y        3  0.4936 0.0798 100  0.3356  0.6446
#> 10        x        4  0.3686 0.0442 100  0.2825  0.4648
#> 11        m        4 -0.2537 0.0725 100 -0.4323 -0.1501
#> 12        y        4  0.4293 0.0856 100  0.2679  0.6023
#> 13        x        5  0.2555 0.0411 100  0.1789  0.3445
#> 14        m        5 -0.1954 0.0681 100 -0.3739 -0.1060
#> 15        y        5  0.3508 0.0896 100  0.1924  0.5406
confint(mc, level = 0.95)
#>    variable interval      2.5 %     97.5 %
#> 1         x        1  0.4747164  0.6174741
#> 2         m        1 -0.3954265 -0.1658916
#> 3         y        1  0.2661921  0.5114527
#> 4         x        2  0.5151977  0.6868680
#> 5         m        2 -0.5115370 -0.2075773
#> 6         y        2  0.3490723  0.6442781
#> 7         x        3  0.4101949  0.5885318
#> 8         m        3 -0.4926175 -0.1928545
#> 9         y        3  0.3355859  0.6446168
#> 10        x        4  0.2824721  0.4648143
#> 11        m        4 -0.4322992 -0.1500715
#> 12        y        4  0.2679152  0.6023399
#> 13        x        5  0.1788787  0.3444951
#> 14        m        5 -0.3738805 -0.1059606
#> 15        y        5  0.1924405  0.5405595
plot(mc)