Skip to contents

Model

The measurement model is given by 𝐲i,t=π›ˆi,t\begin{equation} \mathbf{y}_{i, t} = \boldsymbol{\eta}_{i, t} \end{equation} where 𝐲i,t\mathbf{y}_{i, t} represents a vector of observed variables and π›ˆi,t\boldsymbol{\eta}_{i, t} a vector of latent variables for individual ii and time tt. Since the observed and latent variables are equal, we only generate data from the dynamic structure.

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} .

Data Generation

Notation

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

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=(0.19607840.11832320.02985390.11832320.34377110.13818550.02985390.13818550.2663828).\begin{equation} \boldsymbol{\Sigma}_{\boldsymbol{\eta} \mid 0} = \left( \begin{array}{ccc} 0.1960784 & 0.1183232 & 0.0298539 \\ 0.1183232 & 0.3437711 & 0.1381855 \\ 0.0298539 & 0.1381855 & 0.2663828 \\ \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] 100
time
#> [1] 500
mu0
#> [1] 0 0 0
sigma0
#>            [,1]      [,2]       [,3]
#> [1,] 0.19607843 0.1183232 0.02985385
#> [2,] 0.11832319 0.3437711 0.13818551
#> [3,] 0.02985385 0.1381855 0.26638284
sigma0_l # sigma0_l <- t(chol(sigma0))
#>            [,1]      [,2]     [,3]
#> [1,] 0.44280744 0.0000000 0.000000
#> [2,] 0.26721139 0.5218900 0.000000
#> [3,] 0.06741949 0.2302597 0.456966
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

Visualizing the Dynamics Without Process Noise (n = 5 with Different Initial Condition)

Using the SimSSMVARFixed Function from the simStateSpace Package to Simulate Data

library(simStateSpace)
sim <- SimSSMVARFixed(
  n = n,
  time = time,
  mu0 = mu0,
  sigma0_l = sigma0_l,
  alpha = alpha,
  beta = beta,
  psi_l = psi_l
)
data <- as.data.frame(sim)
head(data)
#>   id time          y1          y2          y3
#> 1  1    0 -0.81728749  0.01319023  0.57975035
#> 2  1    1 -0.62264147 -0.94877859  0.54422379
#> 3  1    2 -0.06731464 -0.69530903 -0.09626706
#> 4  1    3 -0.09174342 -0.44804477 -0.27983813
#> 5  1    4 -0.04104160 -0.72769549 -0.19238403
#> 6  1    5  0.18907246 -0.81602072 -0.73831208
summary(data)
#>        id              time             y1                  y2           
#>  Min.   :  1.00   Min.   :  0.0   Min.   :-1.831372   Min.   :-2.644556  
#>  1st Qu.: 25.75   1st Qu.:124.8   1st Qu.:-0.297127   1st Qu.:-0.395276  
#>  Median : 50.50   Median :249.5   Median :-0.003762   Median :-0.003310  
#>  Mean   : 50.50   Mean   :249.5   Mean   :-0.002119   Mean   :-0.002361  
#>  3rd Qu.: 75.25   3rd Qu.:374.2   3rd Qu.: 0.290142   3rd Qu.: 0.390052  
#>  Max.   :100.00   Max.   :499.0   Max.   : 1.975304   Max.   : 2.318903  
#>        y3            
#>  Min.   :-2.1828250  
#>  1st Qu.:-0.3458073  
#>  Median :-0.0015880  
#>  Mean   :-0.0004972  
#>  3rd Qu.: 0.3443516  
#>  Max.   : 2.1134391
plot(sim)

Model Fitting

Prepare Data

dynr_data <- dynr::dynr.data(
  data = data,
  id = "id",
  time = "time",
  observed = c("y1", "y2", "y3")
)

Prepare Initial Condition

dynr_initial <- dynr::prep.initial(
  values.inistate = mu0,
  params.inistate = c("mu0_1_1", "mu0_2_1", "mu0_3_1"),
  values.inicov = sigma0,
  params.inicov = matrix(
    data = 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
  )
)

Prepare Measurement Model

dynr_measurement <- dynr::prep.measurement(
  values.load = diag(3),
  params.load = matrix(data = "fixed", nrow = 3, ncol = 3),
  state.names = c("eta_1", "eta_2", "eta_3"),
  obs.names = c("y1", "y2", "y3")
)

Prepare Dynamic Process

dynr_dynamics <- dynr::prep.formulaDynamics(
  formula = list(
    eta_1 ~ alpha_1_1 * 1 + beta_1_1 * eta_1 + beta_1_2 * eta_2 + beta_1_3 * eta_3,
    eta_2 ~ alpha_2_1 * 1 + beta_2_1 * eta_1 + beta_2_2 * eta_2 + beta_2_3 * eta_3,
    eta_3 ~ alpha_3_1 * 1 + beta_3_1 * eta_1 + beta_3_2 * eta_2 + beta_3_3 * eta_3
  ),
  startval = c(
    alpha_1_1 = alpha[1], alpha_2_1 = alpha[2], alpha_3_1 = alpha[3],
    beta_1_1 = beta[1, 1], beta_1_2 = beta[1, 2], beta_1_3 = beta[1, 3],
    beta_2_1 = beta[2, 1], beta_2_2 = beta[2, 2], beta_2_3 = beta[2, 3],
    beta_3_1 = beta[3, 1], beta_3_2 = beta[3, 2], beta_3_3 = beta[3, 3]
  ),
  isContinuousTime = FALSE
)

Prepare Process Noise

dynr_noise <- dynr::prep.noise(
  values.latent = psi,
  params.latent = matrix(
    data = c(
      "psi_1_1", "psi_2_1", "psi_3_1",
      "psi_2_1", "psi_2_2", "psi_3_2",
      "psi_3_1", "psi_3_2", "psi_3_3"
    ),
    nrow = 3
  ),
  values.observed = matrix(data = 0, nrow = 3, ncol = 3),
  params.observed = matrix(data = "fixed", nrow = 3, ncol = 3)
)

Prepare the Model

model <- dynr::dynr.model(
  data = dynr_data,
  initial = dynr_initial,
  measurement = dynr_measurement,
  dynamics = dynr_dynamics,
  noise = dynr_noise,
  outfile = "var.c"
)

Fit the Model

results <- dynr::dynr.cook(
  model,
  debug_flag = TRUE,
  verbose = FALSE
)
#> [1] "Get ready!!!!"
#> using C compiler: β€˜gcc (Ubuntu 13.3.0-6ubuntu2~24.04) 13.3.0’
#> Optimization function called.
#> Starting Hessian calculation ...
#> Finished Hessian calculation.
#> Original exit flag:  3 
#> Modified exit flag:  3 
#> Optimization terminated successfully: ftol_rel or ftol_abs was reached. 
#> Original fitted parameters:  -0.0006607231 3.211261e-05 0.0003389175 0.6935267 
#> 0.0005250874 -0.0005643153 0.4984296 0.6011778 -0.002500733 -0.1056564 
#> 0.4056584 0.4971785 -2.308965 -0.007690335 0.004142653 -2.309886 -0.001043527 
#> -2.300854 -0.002408683 0.001315326 0.02286415 -1.674994 0.3533811 -0.060354 
#> -1.299978 0.4453573 -1.556499 
#> 
#> Transformed fitted parameters:  -0.0006607231 3.211261e-05 0.0003389175 
#> 0.6935267 0.0005250874 -0.0005643153 0.4984296 0.6011778 -0.002500733 
#> -0.1056564 0.4056584 0.4971785 0.09936408 -0.000764143 0.0004116308 0.09927844 
#> -0.0001067592 0.1001751 -0.002408683 0.001315326 0.02286415 0.1873094 
#> 0.06619159 -0.01130487 0.2959286 0.1173817 0.2656113 
#> 
#> Doing end processing
#> Successful trial
#> Total Time: 39.83028 
#> Backend Time: 39.8202

Summary

summary(results)
#> Coefficients:
#>              Estimate Std. Error t value   ci.lower   ci.upper Pr(>|t|)    
#> alpha_1_1  -6.607e-04  1.411e-03  -0.468 -3.426e-03  2.105e-03   0.3198    
#> alpha_2_1   3.211e-05  1.411e-03   0.023 -2.732e-03  2.797e-03   0.4909    
#> alpha_3_1   3.389e-04  1.417e-03   0.239 -2.438e-03  3.116e-03   0.4055    
#> beta_1_1    6.935e-01  3.613e-03 191.939  6.864e-01  7.006e-01   <2e-16 ***
#> beta_1_2    5.251e-04  3.035e-03   0.173 -5.424e-03  6.474e-03   0.4313    
#> beta_1_3   -5.643e-04  3.087e-03  -0.183 -6.615e-03  5.486e-03   0.4275    
#> beta_2_1    4.984e-01  3.612e-03 138.009  4.914e-01  5.055e-01   <2e-16 ***
#> beta_2_2    6.012e-01  3.034e-03 198.146  5.952e-01  6.071e-01   <2e-16 ***
#> beta_2_3   -2.501e-03  3.086e-03  -0.810 -8.549e-03  3.547e-03   0.2089    
#> beta_3_1   -1.057e-01  3.628e-03 -29.123 -1.128e-01 -9.855e-02   <2e-16 ***
#> beta_3_2    4.057e-01  3.048e-03 133.101  3.997e-01  4.116e-01   <2e-16 ***
#> beta_3_3    4.972e-01  3.100e-03 160.389  4.911e-01  5.033e-01   <2e-16 ***
#> psi_1_1     9.936e-02  6.290e-04 157.960  9.813e-02  1.006e-01   <2e-16 ***
#> psi_2_1    -7.641e-04  4.446e-04  -1.719 -1.636e-03  1.073e-04   0.0428 *  
#> psi_3_1     4.116e-04  4.466e-04   0.922 -4.637e-04  1.287e-03   0.1783    
#> psi_2_2     9.928e-02  6.285e-04 157.956  9.805e-02  1.005e-01   <2e-16 ***
#> psi_3_2    -1.068e-04  4.464e-04  -0.239 -9.817e-04  7.682e-04   0.4055    
#> psi_3_3     1.002e-01  6.342e-04 157.962  9.893e-02  1.014e-01   <2e-16 ***
#> mu0_1_1    -2.409e-03  4.327e-02  -0.056 -8.721e-02  8.239e-02   0.4778    
#> mu0_2_1     1.315e-03  5.439e-02   0.024 -1.053e-01  1.079e-01   0.4904    
#> mu0_3_1     2.286e-02  5.153e-02   0.444 -7.813e-02  1.239e-01   0.3286    
#> sigma0_1_1  1.873e-01  2.645e-02   7.083  1.355e-01  2.391e-01   <2e-16 ***
#> sigma0_2_1  6.619e-02  2.444e-02   2.709  1.829e-02  1.141e-01   0.0034 ** 
#> sigma0_3_1 -1.130e-02  2.223e-02  -0.509 -5.487e-02  3.226e-02   0.3055    
#> sigma0_2_2  2.959e-01  4.184e-02   7.073  2.139e-01  3.779e-01   <2e-16 ***
#> sigma0_3_2  1.174e-01  3.040e-02   3.861  5.779e-02  1.770e-01   0.0001 ***
#> sigma0_3_3  2.656e-01  3.754e-02   7.076  1.920e-01  3.392e-01   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log-likelihood value at convergence = 79944.89
#> AIC = 79998.89
#> BIC = 80237.02

Parameter Estimates

alpha_hat
#> [1] -6.607231e-04  3.211261e-05  3.389175e-04
beta_hat
#>            [,1]         [,2]          [,3]
#> [1,]  0.6935267 0.0005250874 -0.0005643153
#> [2,]  0.4984296 0.6011777566 -0.0025007326
#> [3,] -0.1056564 0.4056584147  0.4971784809
psi_hat
#>               [,1]          [,2]          [,3]
#> [1,]  0.0993640778 -0.0007641430  0.0004116308
#> [2,] -0.0007641430  0.0992784439 -0.0001067592
#> [3,]  0.0004116308 -0.0001067592  0.1001751099
mu0_hat
#> [1] -0.002408683  0.001315326  0.022864146
sigma0_hat
#>             [,1]       [,2]        [,3]
#> [1,]  0.18730937 0.06619159 -0.01130487
#> [2,]  0.06619159 0.29592857  0.11738173
#> [3,] -0.01130487 0.11738173  0.26561129

References

Chow, S.-M., Ho, M. R., Hamaker, E. L., & Dolan, C. V. (2010). Equivalence and differences between structural equation modeling and state-space modeling techniques. Structural Equation Modeling: A Multidisciplinary Journal, 17(2), 303–332. https://doi.org/10.1080/10705511003661553
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/