Skip to contents

Model

The measurement model is given by 𝐲i,t=π›Ž+πš²π›ˆi,t+𝛆i,t,with𝛆i,tβˆΌπ’©(𝟎,𝚯)\begin{equation} \mathbf{y}_{i, t} = \boldsymbol{\nu} + \boldsymbol{\Lambda} \boldsymbol{\eta}_{i, t} + \boldsymbol{\varepsilon}_{i, t}, \quad \mathrm{with} \quad \boldsymbol{\varepsilon}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \boldsymbol{\Theta} \right) \end{equation} where 𝐲i,t\mathbf{y}_{i, t}, π›ˆi,t\boldsymbol{\eta}_{i, t}, and 𝛆i,t\boldsymbol{\varepsilon}_{i, t} are random variables and π›Ž\boldsymbol{\nu}, 𝚲\boldsymbol{\Lambda}, and 𝚯\boldsymbol{\Theta} are model parameters. 𝐲i,t\mathbf{y}_{i, t} represents a vector of observed random variables, π›ˆi,t\boldsymbol{\eta}_{i, t} a vector of latent random variables, and 𝛆i,t\boldsymbol{\varepsilon}_{i, t} a vector of random measurement errors, at time tt and individual ii. π›Ž\boldsymbol{\nu} denotes a vector of intercepts, 𝚲\boldsymbol{\Lambda} a matrix of factor loadings, and 𝚯\boldsymbol{\Theta} the covariance matrix of 𝛆\boldsymbol{\varepsilon}.

An alternative representation of the measurement error is given by 𝛆i,t=𝚯12𝐳i,t,with𝐳i,tβˆΌπ’©(𝟎,𝐈)\begin{equation} \boldsymbol{\varepsilon}_{i, t} = \boldsymbol{\Theta}^{\frac{1}{2}} \mathbf{z}_{i, t}, \quad \mathrm{with} \quad \mathbf{z}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \mathbf{I} \right) \end{equation} where 𝐳i,t\mathbf{z}_{i, t} is a vector of independent standard normal random variables and (𝚯12)(𝚯12)β€²=𝚯\left( \boldsymbol{\Theta}^{\frac{1}{2}} \right) \left( \boldsymbol{\Theta}^{\frac{1}{2}} \right)^{\prime} = \boldsymbol{\Theta} .

The dynamic structure is given by π›ˆi,t=𝛂+π›ƒπ›ˆi,tβˆ’1+𝛇i,t,with𝛇i,tβˆΌπ’©(𝟎,𝚿)\begin{equation} \boldsymbol{\eta}_{i, t} = \boldsymbol{\alpha} + \boldsymbol{\beta} \boldsymbol{\eta}_{i, t - 1} + \boldsymbol{\zeta}_{i, t}, \quad \mathrm{with} \quad \boldsymbol{\zeta}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \boldsymbol{\Psi} \right) \end{equation} where π›ˆi,t\boldsymbol{\eta}_{i, t}, π›ˆi,tβˆ’1\boldsymbol{\eta}_{i, t - 1}, and 𝛇i,t\boldsymbol{\zeta}_{i, t} are random variables, and 𝛂\boldsymbol{\alpha}, 𝛃\boldsymbol{\beta}, and 𝚿\boldsymbol{\Psi} are model parameters. Here, π›ˆi,t\boldsymbol{\eta}_{i, t} is a vector of latent variables at time tt and individual ii, π›ˆi,tβˆ’1\boldsymbol{\eta}_{i, t - 1} represents a vector of latent variables at time tβˆ’1t - 1 and individual ii, and 𝛇i,t\boldsymbol{\zeta}_{i, t} represents a vector of dynamic noise at time tt and individual ii. 𝛂\boldsymbol{\alpha} denotes a vector of intercepts, 𝛃\boldsymbol{\beta} a matrix of autoregression and cross regression coefficients, and 𝚿\boldsymbol{\Psi} the covariance matrix of 𝛇i,t\boldsymbol{\zeta}_{i, t}.

An alternative representation of the dynamic noise is given by 𝛇i,t=𝚿12𝐳i,t,with𝐳i,tβˆΌπ’©(𝟎,𝐈)\begin{equation} \boldsymbol{\zeta}_{i, t} = \boldsymbol{\Psi}^{\frac{1}{2}} \mathbf{z}_{i, t}, \quad \mathrm{with} \quad \mathbf{z}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \mathbf{I} \right) \end{equation} where (𝚿12)(𝚿12)β€²=𝚿\left( \boldsymbol{\Psi}^{\frac{1}{2}} \right) \left( \boldsymbol{\Psi}^{\frac{1}{2}} \right)^{\prime} = \boldsymbol{\Psi} .

Parameters

Notation

Let t=100t = 100 be the number of time points and n=5n = 5 be the number of individuals.

Let the measurement model intecept vector π›Ž\boldsymbol{\nu} be given by

π›Ž=(000).\begin{equation} \boldsymbol{\nu} = \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) . \end{equation}

Let the factor loadings matrix 𝚲\boldsymbol{\Lambda} be given by

𝚲=(100010001).\begin{equation} \boldsymbol{\Lambda} = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right) . \end{equation}

Let the measurement error covariance matrix 𝚯\boldsymbol{\Theta} be given by

𝚯=(0.20000.20000.2).\begin{equation} \boldsymbol{\Theta} = \left( \begin{array}{ccc} 0.2 & 0 & 0 \\ 0 & 0.2 & 0 \\ 0 & 0 & 0.2 \\ \end{array} \right) . \end{equation}

Let the initial condition π›ˆ0\boldsymbol{\eta}_{0} be given by

π›ˆ0βˆΌπ’©(π›π›ˆβˆ£0,πšΊπ›ˆβˆ£0)\begin{equation} \boldsymbol{\eta}_{0} \sim \mathcal{N} \left( \boldsymbol{\mu}_{\boldsymbol{\eta} \mid 0}, \boldsymbol{\Sigma}_{\boldsymbol{\eta} \mid 0} \right) \end{equation}

π›π›ˆβˆ£0=(000)\begin{equation} \boldsymbol{\mu}_{\boldsymbol{\eta} \mid 0} = \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) \end{equation}

πšΊπ›ˆβˆ£0=(10.20.20.210.20.20.21).\begin{equation} \boldsymbol{\Sigma}_{\boldsymbol{\eta} \mid 0} = \left( \begin{array}{ccc} 1 & 0.2 & 0.2 \\ 0.2 & 1 & 0.2 \\ 0.2 & 0.2 & 1 \\ \end{array} \right) . \end{equation}

Let the constant vector 𝛂\boldsymbol{\alpha} be given by

𝛂=(000).\begin{equation} \boldsymbol{\alpha} = \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) . \end{equation}

Let the transition matrix 𝛃\boldsymbol{\beta} be given by

𝛃=(0.7000.50.60βˆ’0.10.40.5).\begin{equation} \boldsymbol{\beta} = \left( \begin{array}{ccc} 0.7 & 0 & 0 \\ 0.5 & 0.6 & 0 \\ -0.1 & 0.4 & 0.5 \\ \end{array} \right) . \end{equation}

Let the dynamic process noise 𝚿\boldsymbol{\Psi} be given by

𝚿=(0.10000.10000.1).\begin{equation} \boldsymbol{\Psi} = \left( \begin{array}{ccc} 0.1 & 0 & 0 \\ 0 & 0.1 & 0 \\ 0 & 0 & 0.1 \\ \end{array} \right) . \end{equation}

R Function Arguments

n
#> [1] 5
time
#> [1] 100
mu0
#> [1] 0 0 0
sigma0
#>      [,1] [,2] [,3]
#> [1,]  1.0  0.2  0.2
#> [2,]  0.2  1.0  0.2
#> [3,]  0.2  0.2  1.0
sigma0_l # sigma0_l <- t(chol(sigma0))
#>      [,1]      [,2]      [,3]
#> [1,]  1.0 0.0000000 0.0000000
#> [2,]  0.2 0.9797959 0.0000000
#> [3,]  0.2 0.1632993 0.9660918
alpha
#> [1] 0 0 0
beta
#>      [,1] [,2] [,3]
#> [1,]  0.7  0.0  0.0
#> [2,]  0.5  0.6  0.0
#> [3,] -0.1  0.4  0.5
psi
#>      [,1] [,2] [,3]
#> [1,]  0.1  0.0  0.0
#> [2,]  0.0  0.1  0.0
#> [3,]  0.0  0.0  0.1
psi_l # psi_l <- t(chol(psi))
#>           [,1]      [,2]      [,3]
#> [1,] 0.3162278 0.0000000 0.0000000
#> [2,] 0.0000000 0.3162278 0.0000000
#> [3,] 0.0000000 0.0000000 0.3162278
nu
#> [1] 0 0 0
lambda
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    0    1    0
#> [3,]    0    0    1
theta
#>      [,1] [,2] [,3]
#> [1,]  0.2  0.0  0.0
#> [2,]  0.0  0.2  0.0
#> [3,]  0.0  0.0  0.2
theta_l # theta_l <- t(chol(theta))
#>           [,1]      [,2]      [,3]
#> [1,] 0.4472136 0.0000000 0.0000000
#> [2,] 0.0000000 0.4472136 0.0000000
#> [3,] 0.0000000 0.0000000 0.4472136

Parametric Bootstrap

R <- 5L # use at least 1000 in actual research
path <- getwd()
prefix <- "ssm"

We use the PBSSMFixed function from the bootStateSpace package to perform parametric bootstraping using the parameters described above. 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 state space model multiple times is computationally intensive.

library(bootStateSpace)
pb <- PBSSMFixed(
  R = R,
  path = path,
  prefix = prefix,
  n = n,
  time = time,
  mu0 = mu0,
  sigma0_l = sigma0_l,
  alpha = alpha,
  beta = beta,
  psi_l = psi_l,
  nu = nu,
  lambda = lambda,
  theta_l = theta_l,
  ncores = parallel::detectCores(),
  seed = 42
)
summary(pb)
#> Call:
#> PBSSMFixed(R = R, path = path, prefix = prefix, n = n, time = time, 
#>     mu0 = mu0, sigma0_l = sigma0_l, alpha = alpha, beta = beta, 
#>     psi_l = psi_l, nu = nu, lambda = lambda, theta_l = theta_l, 
#>     ncores = parallel::detectCores(), seed = 42)
#> 
#> Parametric bootstrap confidence intervals.
#> type = "pc"
#>             est     se R    2.5%   97.5%
#> beta_1_1    0.7 0.4170 5  0.3261  1.3107
#> beta_2_1    0.5 0.5956 5  0.1454  1.5145
#> beta_3_1   -0.1 0.1380 5 -0.1107  0.1986
#> beta_1_2    0.0 0.2327 5 -0.2269  0.3291
#> beta_2_2    0.6 0.3072 5  0.1526  0.8931
#> beta_3_2    0.4 0.1743 5  0.1203  0.5438
#> beta_1_3    0.0 0.0484 5 -0.0955  0.0258
#> beta_2_3    0.0 0.0407 5 -0.1454 -0.0515
#> beta_3_3    0.5 0.1057 5  0.4244  0.6536
#> psi_1_1     0.1 0.1220 5  0.0108  0.2875
#> psi_2_2     0.1 0.0307 5  0.0567  0.1305
#> psi_3_3     0.1 0.0202 5  0.0124  0.0638
#> theta_1_1   0.2 0.1334 5  0.0048  0.3213
#> theta_2_2   0.2 0.0552 5  0.1362  0.2610
#> theta_3_3   0.2 0.0248 5  0.2147  0.2765
#> mu0_1_1     0.0 0.2816 5 -0.2773  0.3742
#> mu0_2_1     0.0 0.2465 5 -0.3723  0.1683
#> mu0_3_1     0.0 0.4300 5 -0.0405  0.9278
#> sigma0_1_1  1.0 0.8096 5  0.1506  1.7762
#> sigma0_2_1  0.2 0.0816 5  0.2369  0.4238
#> sigma0_3_1  0.2 0.5290 5 -0.3161  1.0294
#> sigma0_2_2  1.0 0.3498 5  0.1730  0.9489
#> sigma0_3_2  0.2 0.4938 5 -0.2606  0.9392
#> sigma0_3_3  1.0 0.6481 5  0.5469  1.9850
summary(pb, type = "bc")
#> Call:
#> PBSSMFixed(R = R, path = path, prefix = prefix, n = n, time = time, 
#>     mu0 = mu0, sigma0_l = sigma0_l, alpha = alpha, beta = beta, 
#>     psi_l = psi_l, nu = nu, lambda = lambda, theta_l = theta_l, 
#>     ncores = parallel::detectCores(), seed = 42)
#> 
#> Parametric bootstrap confidence intervals.
#> type = "bc"
#>             est     se R    2.5%   97.5%
#> beta_1_1    0.7 0.4170 5  0.5577  1.3789
#> beta_2_1    0.5 0.5956 5  0.1712  1.5954
#> beta_3_1   -0.1 0.1380 5 -0.1131 -0.0559
#> beta_1_2    0.0 0.2327 5 -0.2435  0.3136
#> beta_2_2    0.6 0.3072 5  0.1017  0.7274
#> beta_3_2    0.4 0.1743 5  0.1686  0.5489
#> beta_1_3    0.0 0.0484 5 -0.0494  0.0304
#> beta_2_3    0.0 0.0407 5 -0.0504 -0.0504
#> beta_3_3    0.5 0.1057 5  0.4241  0.6424
#> psi_1_1     0.1 0.1220 5  0.0001  0.1906
#> psi_2_2     0.1 0.0307 5  0.0590  0.1329
#> psi_3_3     0.1 0.0202 5  0.0658  0.0658
#> theta_1_1   0.2 0.1334 5  0.0139  0.3309
#> theta_2_2   0.2 0.0552 5  0.1337  0.2567
#> theta_3_3   0.2 0.0248 5  0.2110  0.2110
#> mu0_1_1     0.0 0.2816 5 -0.3137  0.1883
#> mu0_2_1     0.0 0.2465 5 -0.4163  0.1375
#> mu0_3_1     0.0 0.4300 5 -0.0510  0.4663
#> sigma0_1_1  1.0 0.8096 5  0.1522  1.7854
#> sigma0_2_1  0.2 0.0816 5  0.2347  0.2347
#> sigma0_3_1  0.2 0.5290 5 -0.3577  0.8921
#> sigma0_2_2  1.0 0.3498 5  0.9555  0.9555
#> sigma0_3_2  0.2 0.4938 5 -0.2970  0.8925
#> sigma0_3_3  1.0 0.6481 5  0.5782  2.0099

References

Ou, L., Hunter, M. D., & Chow, S.-M. (2019). What’s for dynr: A package for linear and nonlinear dynamic modeling in R. The R Journal, 11(1), 91. https://doi.org/10.32614/rj-2019-012
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. https://doi.org/10.1037/met0000779
R Core Team. (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/