Skip to contents

Data Generation Using the SimSSMOUFixed Function from the simStateSpace Package

n
#> [1] 100
time
#> [1] 1000
delta_t
#> [1] 0.1
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
mu
#> [1] 0 0 0
phi
#>        [,1]   [,2]   [,3]
#> [1,] -0.357  0.000  0.000
#> [2,]  0.771 -0.511  0.000
#> [3,] -0.450  0.729 -0.693
sigma
#>             [,1]       [,2]        [,3]
#> [1,]  0.24455556 0.02201587 -0.05004762
#> [2,]  0.02201587 0.07067800  0.01539456
#> [3,] -0.05004762 0.01539456  0.07553061
sigma_l # sigma_l <- t(chol(sigma))
#>             [,1]      [,2]     [,3]
#> [1,]  0.49452559 0.0000000 0.000000
#> [2,]  0.04451917 0.2620993 0.000000
#> [3,] -0.10120330 0.0759256 0.243975
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
library(simStateSpace)
sim <- SimSSMOUFixed(
  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,
  type = 0
)
data <- as.data.frame(sim)
colnames(data) <- c("id", "time", "x", "m", "y")
head(data)
#>   id time          x          m        y
#> 1  1  0.0 -0.3504435 0.41877429 2.611996
#> 2  1  0.1 -0.5920330 1.07433208 1.669272
#> 3  1  0.2 -0.7619855 1.21483834 2.369837
#> 4  1  0.3 -1.6964652 0.21209722 2.128531
#> 5  1  0.4 -1.2282686 0.09950326 1.891140
#> 6  1  0.5  0.1433985 0.66784226 2.036033

Model Fitting

We use the dynr package to fit the continuous-time vector autoregressive model.

Prepare the Data

dynr_data <- dynr.data(
  dataframe = data,
  id = "id",
  time = "time",
  observed = c("x", "m", "y")
)

Prepare the Initial Condition

dynr_initial <- prep.initial(
  values.inistate = c(
    0, 0, 0
  ),
  params.inistate = c(
    "mu0_1_1", "mu0_2_1", "mu0_3_1"
  ),
  values.inicov = matrix(
    data = c(
      1.0, 0.2, 0.2,
      0.2, 1.0, 0.2,
      0.2, 0.2, 1.0
    ),
    nrow = 3
  ),
  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 the Measurement Model

dynr_measurement <- prep.measurement(
  values.load = diag(3),
  params.load = matrix(
    data = "fixed",
    nrow = 3,
    ncol = 3
  ),
  state.names = c("eta_x", "eta_m", "eta_y"),
  obs.names = c("x", "m", "y")
)

Prepare the Dynamic Model

dynr_dynamics <- prep.formulaDynamics(
  formula = list(
    eta_x ~ phi_1_1 * eta_x + phi_1_2 * eta_m + phi_1_3 * eta_y,
    eta_m ~ phi_2_1 * eta_x + phi_2_2 * eta_m + phi_2_3 * eta_y,
    eta_y ~ phi_3_1 * eta_x + phi_3_2 * eta_m + phi_3_3 * eta_y
  ),
  startval = c(
    phi_1_1 = -0.2, phi_2_1 = 0.0, phi_3_1 = 0.0,
    phi_1_2 = 0.0, phi_2_2 = -0.2, phi_3_2 = 0.0,
    phi_1_3 = 0.0, phi_2_3 = 0.0, phi_3_3 = -0.2
  ),
  isContinuousTime = TRUE
)

Prepare the Noise Matrices

dynr_noise <- prep.noise(
  values.latent = matrix(
    data = 0.2 * diag(3),
    nrow = 3
  ),
  params.latent = matrix(
    data = 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
  ),
  values.observed = 0.2 * diag(3),
  params.observed = matrix(
    data = c(
      "theta_1_1", "fixed", "fixed",
      "fixed", "theta_2_2", "fixed",
      "fixed", "fixed", "theta_3_3"
    ),
    nrow = 3
  )
)

Prepare the Model

In this example, we increase the maximum evaluations in the optimization process and constrain the lower and upper bounds of the parameters for the drift matrix.

dynr_model <- dynr.model(
  data = dynr_data,
  initial = dynr_initial,
  measurement = dynr_measurement,
  dynamics = dynr_dynamics,
  noise = dynr_noise,
  outfile = "fit-ct-var-dynr.c"
)
dynr_model@options$maxeval <- 100000
lb <- ub <- rep(NA, times = length(dynr_model$xstart))
names(ub) <- names(lb) <- names(dynr_model$xstart)
lb[
  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"
  )
] <- -1.5
ub[
  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"
  )
] <- 1.5
dynr_model$lb <- lb
dynr_model$ub <- ub

Fit the Model

fit <- dynr.cook(
  dynr_model,
  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.3518392 0.7442816 -0.4586796 0.01731083 
#> -0.4888207 0.7267998 -0.02381434 -0.009810166 -0.6883342 -1.418073 0.09609769 
#> -0.2088286 -2.681139 0.2898059 -2.881283 -1.615149 -1.611839 -1.603597 
#> 0.006324029 -0.04252988 0.1300433 0.1400118 0.3596041 0.1964667 0.07065217 
#> 0.1435497 -0.1097269 
#> 
#> Transformed fitted parameters:  -0.3518392 0.7442816 -0.4586796 0.01731083 
#> -0.4888207 0.7267998 -0.02381434 -0.009810166 -0.6883342 0.2421803 0.02327296 
#> -0.05057416 0.07072156 0.01498732 0.07237598 0.1988611 0.1995203 0.2011716 
#> 0.006324029 -0.04252988 0.1300433 1.150287 0.413648 0.2259932 1.221957 
#> 0.2353267 0.962594 
#> 
#> Doing end processing
#> Successful trial
#> Total Time: 42.69519 
#> Backend Time: 42.68898
summary(fit)
#> Coefficients:
#>             Estimate Std. Error t value  ci.lower  ci.upper Pr(>|t|)    
#> phi_1_1    -0.351839   0.036416  -9.662 -0.423213 -0.280465   <2e-16 ***
#> phi_2_1     0.744282   0.021777  34.177  0.701599  0.786964   <2e-16 ***
#> phi_3_1    -0.458680   0.023534 -19.490 -0.504806 -0.412554   <2e-16 ***
#> phi_1_2     0.017311   0.031705   0.546 -0.044829  0.079451   0.2925    
#> phi_2_2    -0.488821   0.019277 -25.358 -0.526602 -0.451039   <2e-16 ***
#> phi_3_2     0.726800   0.020871  34.824  0.685894  0.767706   <2e-16 ***
#> phi_1_3    -0.023814   0.024025  -0.991 -0.070903  0.023275   0.1608    
#> phi_2_3    -0.009810   0.014718  -0.667 -0.038657  0.019036   0.2525    
#> phi_3_3    -0.688334   0.016040 -42.913 -0.719773 -0.656896   <2e-16 ***
#> sigma_1_1   0.242180   0.006794  35.646  0.228864  0.255496   <2e-16 ***
#> sigma_2_1   0.023273   0.002545   9.146  0.018285  0.028261   <2e-16 ***
#> sigma_3_1  -0.050574   0.002749 -18.395 -0.055963 -0.045186   <2e-16 ***
#> sigma_2_2   0.070722   0.001907  37.093  0.066985  0.074458   <2e-16 ***
#> sigma_3_2   0.014987   0.001381  10.854  0.012281  0.017694   <2e-16 ***
#> sigma_3_3   0.072376   0.002099  34.475  0.068261  0.076491   <2e-16 ***
#> theta_1_1   0.198861   0.001170 169.909  0.196567  0.201155   <2e-16 ***
#> theta_2_2   0.199520   0.001000 199.500  0.197560  0.201480   <2e-16 ***
#> theta_3_3   0.201172   0.001016 198.052  0.199181  0.203162   <2e-16 ***
#> mu0_1_1     0.006324   0.111110   0.057 -0.211447  0.224095   0.4773    
#> mu0_2_1    -0.042530   0.114320  -0.372 -0.266593  0.181533   0.3549    
#> mu0_3_1     0.130043   0.102109   1.274 -0.070086  0.330172   0.1014    
#> sigma0_1_1  1.150287   0.168811   6.814  0.819425  1.481150   <2e-16 ***
#> sigma0_2_1  0.413648   0.133495   3.099  0.152003  0.675293   0.0010 ***
#> sigma0_3_1  0.225993   0.123478   1.830 -0.016019  0.468006   0.0336 *  
#> sigma0_2_2  1.221957   0.182233   6.705  0.864787  1.579128   <2e-16 ***
#> sigma0_3_2  0.235327   0.117629   2.001  0.004779  0.465875   0.0227 *  
#> sigma0_3_3  0.962594   0.142152   6.772  0.683981  1.241207   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log-likelihood value at convergence = 429365.49
#> AIC = 429419.49
#> BIC = 429676.34
coefs <- coef(fit)
vcovs <- vcov(fit)

Extract Matrices from the Fitted Model to use in cTMed

phi_names <- 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"
)
sigma_names <- 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"
)
sigma_vech_names <- c(
  "sigma_1_1", "sigma_2_1", "sigma_3_1",
  "sigma_2_2", "sigma_3_2",
  "sigma_3_3"
)
theta_names <- 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",
  "sigma_1_1", "sigma_2_1", "sigma_3_1",
  "sigma_2_2", "sigma_3_2",
  "sigma_3_3"
)
phi <- matrix(
  data = coefs[phi_names],
  nrow = 3,
  ncol = 3
)
sigma <- matrix(
  data = coefs[sigma_names],
  nrow = 3,
  ncol = 3
)
theta <- coefs[theta_names]
vcov_phi_vec <- vcovs[phi_names, phi_names]
vcov_sigma_vech <- vcovs[sigma_vech_names, sigma_vech_names]
vcov_theta <- vcovs[theta_names, theta_names]

Estimated Drift Matrix with Corresponding Sampling Covariance Matrix

phi
#>            [,1]        [,2]         [,3]
#> [1,] -0.3518392  0.01731083 -0.023814339
#> [2,]  0.7442816 -0.48882067 -0.009810166
#> [3,] -0.4586796  0.72679980 -0.688334177
vcov_phi_vec
#>               phi_1_1       phi_2_1       phi_3_1       phi_1_2       phi_2_2
#> phi_1_1  1.326121e-03  9.158296e-05 -2.258193e-04 -1.108000e-03 -8.829144e-05
#> phi_2_1  9.158296e-05  4.742430e-04 -3.596064e-06 -6.566903e-05 -4.021299e-04
#> phi_3_1 -2.258193e-04 -3.596064e-06  5.538551e-04  1.845433e-04  1.077730e-05
#> phi_1_2 -1.108000e-03 -6.566903e-05  1.845433e-04  1.005190e-03  7.060966e-05
#> phi_2_2 -8.829144e-05 -4.021299e-04  1.077730e-05  7.060966e-05  3.715859e-04
#> phi_3_2  1.994853e-04 -3.470399e-06 -4.716662e-04 -1.780693e-04 -3.561584e-06
#> phi_1_3  7.414387e-04  3.913166e-05 -1.151977e-04 -6.965449e-04 -4.344659e-05
#> phi_2_3  6.974956e-05  2.704005e-04 -1.452237e-05 -5.931610e-05 -2.595864e-04
#> phi_3_3 -1.424724e-04  7.176532e-06  3.208390e-04  1.318148e-04 -2.963173e-06
#>               phi_3_2       phi_1_3       phi_2_3       phi_3_3
#> phi_1_1  1.994853e-04  7.414387e-04  6.974956e-05 -1.424724e-04
#> phi_2_1 -3.470399e-06  3.913166e-05  2.704005e-04  7.176532e-06
#> phi_3_1 -4.716662e-04 -1.151977e-04 -1.452237e-05  3.208390e-04
#> phi_1_2 -1.780693e-04 -6.965449e-04 -5.931610e-05  1.318148e-04
#> phi_2_2 -3.561584e-06 -4.344659e-05 -2.595864e-04 -2.963173e-06
#> phi_3_2  4.355884e-04  1.163975e-04  1.028072e-05 -3.075120e-04
#> phi_1_3  1.163975e-04  5.772234e-04  4.373490e-05 -1.045143e-04
#> phi_2_3  1.028072e-05  4.373490e-05  2.166149e-04 -2.958830e-06
#> phi_3_3 -3.075120e-04 -1.045143e-04 -2.958830e-06  2.572924e-04

Process Noise Covariance Matrix with Corresponding Sampling Covariance Matrix

sigma
#>             [,1]       [,2]        [,3]
#> [1,]  0.24218026 0.02327296 -0.05057416
#> [2,]  0.02327296 0.07072156  0.01498732
#> [3,] -0.05057416 0.01498732  0.07237598
vcov_sigma_vech
#>               sigma_1_1     sigma_2_1     sigma_3_1     sigma_2_2     sigma_3_2
#> sigma_1_1  4.616001e-05 -3.795706e-07 -5.585381e-06 -1.692509e-07  3.054065e-07
#> sigma_2_1 -3.795706e-07  6.475588e-06  1.002384e-07 -4.387889e-07 -9.708718e-07
#> sigma_3_1 -5.585381e-06  1.002384e-07  7.558692e-06  9.058849e-08 -3.163361e-07
#> sigma_2_2 -1.692509e-07 -4.387889e-07  9.058849e-08  3.635137e-06  3.523268e-08
#> sigma_3_2  3.054065e-07 -9.708718e-07 -3.163361e-07  3.523268e-08  1.906780e-06
#> sigma_3_3  6.193436e-07  3.360567e-08 -1.947016e-06 -4.509608e-08  4.773520e-08
#>               sigma_3_3
#> sigma_1_1  6.193436e-07
#> sigma_2_1  3.360567e-08
#> sigma_3_1 -1.947016e-06
#> sigma_2_2 -4.509608e-08
#> sigma_3_2  4.773520e-08
#> sigma_3_3  4.407331e-06

Estimated Drift Matrix and Process Noise Covariance Matrix with Corresponding Sampling Covariance Matrix

theta
#>      phi_1_1      phi_2_1      phi_3_1      phi_1_2      phi_2_2      phi_3_2 
#> -0.351839184  0.744281649 -0.458679613  0.017310833 -0.488820674  0.726799800 
#>      phi_1_3      phi_2_3      phi_3_3    sigma_1_1    sigma_2_1    sigma_3_1 
#> -0.023814339 -0.009810166 -0.688334177  0.242180259  0.023272963 -0.050574160 
#>    sigma_2_2    sigma_3_2    sigma_3_3 
#>  0.070721559  0.014987323  0.072375983
vcov_theta
#>                 phi_1_1       phi_2_1       phi_3_1       phi_1_2       phi_2_2
#> phi_1_1    1.326121e-03  9.158296e-05 -2.258193e-04 -1.108000e-03 -8.829144e-05
#> phi_2_1    9.158296e-05  4.742430e-04 -3.596064e-06 -6.566903e-05 -4.021299e-04
#> phi_3_1   -2.258193e-04 -3.596064e-06  5.538551e-04  1.845433e-04  1.077730e-05
#> phi_1_2   -1.108000e-03 -6.566903e-05  1.845433e-04  1.005190e-03  7.060966e-05
#> phi_2_2   -8.829144e-05 -4.021299e-04  1.077730e-05  7.060966e-05  3.715859e-04
#> phi_3_2    1.994853e-04 -3.470399e-06 -4.716662e-04 -1.780693e-04 -3.561584e-06
#> phi_1_3    7.414387e-04  3.913166e-05 -1.151977e-04 -6.965449e-04 -4.344659e-05
#> phi_2_3    6.974956e-05  2.704005e-04 -1.452237e-05 -5.931610e-05 -2.595864e-04
#> phi_3_3   -1.424724e-04  7.176532e-06  3.208390e-04  1.318148e-04 -2.963173e-06
#> sigma_1_1 -1.982932e-04 -2.184773e-05  3.456919e-05  1.543789e-04  1.849510e-05
#> sigma_2_1  1.502203e-05 -3.437997e-05 -6.620570e-06 -1.591184e-05  2.574007e-05
#> sigma_3_1  1.794495e-05  4.674418e-06 -4.279970e-05 -1.117672e-05 -3.814434e-06
#> sigma_2_2  4.361454e-07  1.424515e-05  1.010140e-06 -1.579752e-07 -1.451605e-05
#> sigma_3_2 -3.741407e-06  1.607279e-06  8.491226e-06  3.585007e-06  4.201294e-07
#> sigma_3_3 -1.043117e-06 -8.607286e-07  3.442663e-06  7.089937e-08  4.858637e-07
#>                 phi_3_2       phi_1_3       phi_2_3       phi_3_3     sigma_1_1
#> phi_1_1    1.994853e-04  7.414387e-04  6.974956e-05 -1.424724e-04 -1.982932e-04
#> phi_2_1   -3.470399e-06  3.913166e-05  2.704005e-04  7.176532e-06 -2.184773e-05
#> phi_3_1   -4.716662e-04 -1.151977e-04 -1.452237e-05  3.208390e-04  3.456919e-05
#> phi_1_2   -1.780693e-04 -6.965449e-04 -5.931610e-05  1.318148e-04  1.543789e-04
#> phi_2_2   -3.561584e-06 -4.344659e-05 -2.595864e-04 -2.963173e-06  1.849510e-05
#> phi_3_2    4.355884e-04  1.163975e-04  1.028072e-05 -3.075120e-04 -2.844219e-05
#> phi_1_3    1.163975e-04  5.772234e-04  4.373490e-05 -1.045143e-04 -9.445101e-05
#> phi_2_3    1.028072e-05  4.373490e-05  2.166149e-04 -2.958830e-06 -1.266847e-05
#> phi_3_3   -3.075120e-04 -1.045143e-04 -2.958830e-06  2.572924e-04  1.859927e-05
#> sigma_1_1 -2.844219e-05 -9.445101e-05 -1.266847e-05  1.859927e-05  4.616001e-05
#> sigma_2_1  6.267263e-06  1.038566e-05 -1.506230e-05 -4.238014e-06 -3.795706e-07
#> sigma_3_1  3.275583e-05  1.446600e-06  2.154429e-06 -1.897455e-05 -5.585381e-06
#> sigma_2_2 -1.392458e-06  1.948765e-08  9.194261e-06  1.027499e-06 -1.692509e-07
#> sigma_3_2 -8.374628e-06 -2.326402e-06 -2.643065e-06  5.086543e-06  3.054065e-07
#> sigma_3_3  6.621115e-07  1.340588e-06  1.970612e-07 -5.788873e-06  6.193436e-07
#>               sigma_2_1     sigma_3_1     sigma_2_2     sigma_3_2     sigma_3_3
#> phi_1_1    1.502203e-05  1.794495e-05  4.361454e-07 -3.741407e-06 -1.043117e-06
#> phi_2_1   -3.437997e-05  4.674418e-06  1.424515e-05  1.607279e-06 -8.607286e-07
#> phi_3_1   -6.620570e-06 -4.279970e-05  1.010140e-06  8.491226e-06  3.442663e-06
#> phi_1_2   -1.591184e-05 -1.117672e-05 -1.579752e-07  3.585007e-06  7.089937e-08
#> phi_2_2    2.574007e-05 -3.814434e-06 -1.451605e-05  4.201294e-07  4.858637e-07
#> phi_3_2    6.267263e-06  3.275583e-05 -1.392458e-06 -8.374628e-06  6.621115e-07
#> phi_1_3    1.038566e-05  1.446600e-06  1.948765e-08 -2.326402e-06  1.340588e-06
#> phi_2_3   -1.506230e-05  2.154429e-06  9.194261e-06 -2.643065e-06  1.970612e-07
#> phi_3_3   -4.238014e-06 -1.897455e-05  1.027499e-06  5.086543e-06 -5.788873e-06
#> sigma_1_1 -3.795706e-07 -5.585381e-06 -1.692509e-07  3.054065e-07  6.193436e-07
#> sigma_2_1  6.475588e-06  1.002384e-07 -4.387889e-07 -9.708718e-07  3.360567e-08
#> sigma_3_1  1.002384e-07  7.558692e-06  9.058849e-08 -3.163361e-07 -1.947016e-06
#> sigma_2_2 -4.387889e-07  9.058849e-08  3.635137e-06  3.523268e-08 -4.509608e-08
#> sigma_3_2 -9.708718e-07 -3.163361e-07  3.523268e-08  1.906780e-06  4.773520e-08
#> sigma_3_3  3.360567e-08 -1.947016e-06 -4.509608e-08  4.773520e-08  4.407331e-06

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/