Total, Direct, and Indirect Effects in Continuous-Time Mediation Model (Bootstrap)
Ivan Jacob Agaloos Pesigan
2025-01-19
Source:vignettes/med-boot.Rmd
med-boot.RmdThe cTMed package provides a bootstrap approach, in
addition to the delta and Monte Carlo methods, for estimating and
quantifying uncertainty in total, direct, and indirect effects within
continuous-time mediation models across different time intervals.
In this example, we will use the fitted model from Fit
the Continuous-Time Vector Autoregressive Model Using the dynr
Package. The object fit represents a fitted CT-VAR
model created using the dynr package.
summary(fit)
#> Coefficients:
#> Estimate Std. Error t value ci.lower ci.upper Pr(>|t|)
#> phi_1_1 -0.293807 0.082574 -3.558 -0.455648 -0.131966 0.0002 ***
#> phi_2_1 0.824396 0.056914 14.485 0.712846 0.935946 <2e-16 ***
#> phi_3_1 -0.454316 0.057647 -7.881 -0.567303 -0.341329 <2e-16 ***
#> phi_1_2 -0.071916 0.068787 -1.045 -0.206735 0.062904 0.1479
#> phi_2_2 -0.568999 0.048532 -11.724 -0.664119 -0.473878 <2e-16 ***
#> phi_3_2 0.707647 0.048636 14.550 0.612322 0.802972 <2e-16 ***
#> phi_1_3 0.025791 0.062589 0.412 -0.096882 0.148463 0.3402
#> phi_2_3 0.082668 0.043529 1.899 -0.002647 0.167983 0.0288 *
#> phi_3_3 -0.687502 0.044092 -15.593 -0.773920 -0.601084 <2e-16 ***
#> sigma_1_1 0.235806 0.022759 10.361 0.191198 0.280413 <2e-16 ***
#> sigma_2_1 0.018668 0.010045 1.858 -0.001020 0.038356 0.0316 *
#> sigma_3_1 -0.064687 0.010768 -6.007 -0.085793 -0.043582 <2e-16 ***
#> sigma_2_2 0.077088 0.009447 8.160 0.058572 0.095603 <2e-16 ***
#> sigma_3_2 0.008278 0.006655 1.244 -0.004765 0.021321 0.1068
#> sigma_3_3 0.080654 0.010196 7.910 0.060670 0.100639 <2e-16 ***
#> theta_1_1 0.201033 0.005107 39.363 0.191023 0.211043 <2e-16 ***
#> theta_2_2 0.191640 0.004373 43.823 0.183069 0.200211 <2e-16 ***
#> theta_3_3 0.197449 0.004521 43.677 0.188588 0.206309 <2e-16 ***
#> mu0_1_1 -0.068547 0.142210 -0.482 -0.347273 0.210180 0.3149
#> mu0_2_1 0.084573 0.176685 0.479 -0.261722 0.430869 0.3161
#> mu0_3_1 0.115721 0.141485 0.818 -0.161585 0.393026 0.2067
#> sigma0_1_1 0.942679 0.205369 4.590 0.540164 1.345195 <2e-16 ***
#> sigma0_2_1 0.161788 0.179880 0.899 -0.190769 0.514346 0.1842
#> sigma0_3_1 0.103257 0.143856 0.718 -0.178696 0.385209 0.2365
#> sigma0_2_2 1.489721 0.311204 4.787 0.879772 2.099670 <2e-16 ***
#> sigma0_3_2 -0.100464 0.178078 -0.564 -0.449491 0.248563 0.2863
#> sigma0_3_3 0.946144 0.204219 4.633 0.545882 1.346406 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log-likelihood value at convergence = 21626.57
#> AIC = 21680.57
#> BIC = 21856.53We need to extract the estimated parameters from the fitted object, which will be used to generate bootstrap samples.
est <- coef(fit)
n
#> [1] 50
time
#> [1] 100
delta_t
#> [1] 0.1
lambda
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 1 0
#> [3,] 0 0 1
nu
#> [1] 0 0 0
mu
#> [1] 0 0 0
mu0 <- est[
c(
"mu0_1_1",
"mu0_2_1",
"mu0_3_1"
)
]
mu0
#> mu0_1_1 mu0_2_1 mu0_3_1
#> -0.06854665 0.08457349 0.11572075
sigma0 <- matrix(
data = est[
c(
"sigma0_1_1",
"sigma0_2_1",
"sigma0_3_1",
"sigma0_2_1",
"sigma0_2_2",
"sigma0_3_2",
"sigma0_3_1",
"sigma0_3_2",
"sigma0_3_3"
)
],
nrow = 3,
ncol = 3
)
sigma0
#> [,1] [,2] [,3]
#> [1,] 0.9426793 0.1617881 0.1032567
#> [2,] 0.1617881 1.4897213 -0.1004635
#> [3,] 0.1032567 -0.1004635 0.9461439
sigma0_l <- t(chol(sigma0))
phi <- matrix(
data = est[
c(
"phi_1_1",
"phi_2_1",
"phi_3_1",
"phi_1_2",
"phi_2_2",
"phi_3_2",
"phi_1_3",
"phi_2_3",
"phi_3_3"
)
],
nrow = 3,
ncol = 3
)
phi
#> [,1] [,2] [,3]
#> [1,] -0.2938068 -0.07191554 0.02579061
#> [2,] 0.8243958 -0.56899879 0.08266777
#> [3,] -0.4543161 0.70764697 -0.68750180
sigma <- matrix(
data = est[
c(
"sigma_1_1", "sigma_2_1", "sigma_3_1",
"sigma_2_1", "sigma_2_2", "sigma_3_2",
"sigma_3_1", "sigma_3_2", "sigma_3_3"
)
],
nrow = 3,
ncol = 3
)
sigma
#> [,1] [,2] [,3]
#> [1,] 0.23580555 0.018668186 -0.064687273
#> [2,] 0.01866819 0.077087752 0.008277734
#> [3,] -0.06468727 0.008277734 0.080654489
sigma_l <- t(chol(sigma))
theta <- diag(3)
diag(theta) <- est[
c(
"theta_1_1",
"theta_2_2",
"theta_3_3"
)
]
theta
#> [,1] [,2] [,3]
#> [1,] 0.2010327 0.0000000 0.0000000
#> [2,] 0.0000000 0.1916402 0.0000000
#> [3,] 0.0000000 0.0000000 0.1974485
theta_l <- t(chol(theta))
R <- 1000L # use at least 1000 in actual research
path <- getwd()
prefix <- "ou"The estimated parameters are then passed as arguments to the
PBSSMOUFixed function from the bootStateSpace
package, which generates a parametric bootstrap sampling distribution of
the parameter estimates. The argument R specifies the
number of bootstrap replications. The generated data and model estimates
are stored in path using the specified prefix
for the file names. The ncores = parallel::detectCores()
argument instructs the function to use all available CPU cores in the
system.
NOTE: Fitting the CT-VAR model multiple times is computationally intensive.
library(bootStateSpace)
boot <- PBSSMOUFixed(
R = R,
path = path,
prefix = prefix,
n = n,
time = time,
delta_t = delta_t,
mu0 = mu0,
sigma0_l = sigma0_l,
mu = mu,
phi = phi,
sigma_l = sigma_l,
nu = nu,
lambda = lambda,
theta_l = theta_l,
ncores = parallel::detectCores(),
seed = 42
)The extract function from the
bootStateSpace package is used to extract the bootstrap phi
matrices as well as the sigma matrices.
phi <- extract(object = boot, what = "phi")
sigma <- extract(object = boot, what = "sigma")In this example, we aim to calculate the total, direct, and indirect
effects of x on y, mediated through
m, over time intervals ranging from 0 to 10.
# time intervals
delta_t <- seq(from = 0, to = 10, length.out = 1000)We also need the estimated drift matrix from the original sample.
# estimated drift matrix
phi_hat <- matrix(
data = est[
c(
"phi_1_1",
"phi_2_1",
"phi_3_1",
"phi_1_2",
"phi_2_2",
"phi_3_2",
"phi_1_3",
"phi_2_3",
"phi_3_3"
)
],
nrow = 3,
ncol = 3
)
colnames(phi_hat) <- rownames(phi_hat) <- c("x", "m", "y")For the standardized effects, the estimated process noise covariance matrix from the original sample is also needed.
# estimated process noise covariance matrix
sigma_hat <- matrix(
data = est[
c(
"sigma_1_1", "sigma_2_1", "sigma_3_1",
"sigma_2_1", "sigma_2_2", "sigma_3_2",
"sigma_3_1", "sigma_3_2", "sigma_3_3"
)
],
nrow = 3,
ncol = 3
)Bootstrap Method
library(cTMed)
boot <- BootMed(
phi = phi,
phi_hat = phi_hat,
delta_t = delta_t,
from = "x",
to = "y",
med = "m",
ncores = parallel::detectCores() # use multiple cores
)
plot(boot)
plot(boot, type = "bc")
The following generates bootstrap confidence intervals for the standardized effects.
boot <- BootMedStd(
phi = phi,
sigma = sigma,
phi_hat = phi_hat,
sigma_hat = sigma_hat,
delta_t = delta_t,
from = "x",
to = "y",
med = "m",
ncores = parallel::detectCores() # use multiple cores
)
plot(boot)
plot(boot, type = "bc")
