Parameter-transformations

To allow for more flexible model estimation, lessSEM offers parameter transformations. This is a feature inspired by OpenMx (Neale et al., 2016; Pritikin et al., 2015; Hunter 2018). Parameter transformations can be very powerful and simplify implementing some specific forms of regularization.

Motivation

In longitudinal SEM, it is important to investigate if parameters stay the same over time (e.g., measurement invariance of loadings). This can be difficult to decide and may require setting up many different models manually. Here, regularization techniques can be very handy. For instance, in the seminal political democracy example, the model is typically set up as follows (see ?lavaan::sem):

modelSyntax <- ' 
  # latent variable definitions
     ind60 =~ x1 + x2 + x3
     dem60 =~ y1 + a*y2 + b*y3 + c*y4
     dem65 =~ y5 + a*y6 + b*y7 + c*y8

  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60

  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'

Note that the loadings a, b, and c are assumed to stay the same over time. That is, measurement invariance is assumed! Relaxing this assumption, we could define the model as follows:

modelSyntax <- ' 
  # latent variable definitions
     ind60 =~ x1 + x2 + x3
     dem60 =~ y1 + a1*y2 + b1*y3 + c1*y4
     dem65 =~ y5 + a2*y6 + b2*y7 + c2*y8

  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60

  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'

Here, each loading is estimated separately. This results in a more complex model. How do we know which model to use? There are many procedures to answer this question (e.g., using modification indexes, setting up separate models by hand, etc.). In the following, we will show how regularization could be used (see e.g., Belzak & Bauer, 2020; Jacobucci & Grimm, 2018).

Using Regularization

First, note that measurement invariance can be rephrased as \(a_1-a_2 = 0\), \(b_1-b_2 = 0\), and \(c_1-c_2 = 0\). Thus, regularizing the differences between these parameters may allow for testing measurement invariance (e.g., Belzak & Bauer, 2020; Liang et al., 2018; Muthen & Asparouhov, 2013). In fact, this is used in Bayesian SEM to test approximate measurement invariance (Liang et al., 2018; Muthen & Asparouhov, 2013). Similar procedures have also been developed by Huang (2018) for multi-group differences in parameter estimates and by Fisher et al. (2022) in vector autoregressive models. Furthermore, Jacobucci & Grimm (2018) proposed regularizing differences in latent change score models to test equivalence of autoproportion parameters over time using a two-step procedure. To this end, there they implemented the diff_lasso in regsem (Jacobucci et al., 2019). Such a diff_lasso is not available in lessSEM. Instead, lessSEM provides a flexible workaround: parameter transformations. To make this work, we have to re-define our parameters.

Redefine:

\[ \begin{align} a_2 &= a_1 + \Delta a_2\\ b_2 &= b_1 + \Delta b_2\\ c_2 &= c_1 + \Delta c_2 \end{align} \] By regularizing \(\Delta a_2\), \(\Delta b_2\), and \(\Delta c_2\) towards zero, we can enforce measurement invariance over time.

Setting up the Model

We first start with the most flexible model which we want to test:

library(lavaan)
modelSyntax <- ' 
  # latent variable definitions
     ind60 =~ x1 + x2 + x3
     dem60 =~ y1 + a1*y2 + b1*y3 + c1*y4
     dem65 =~ y5 + a2*y6 + b2*y7 + c2*y8

  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60

  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'

lavaanFit <- sem(model = modelSyntax,
                 data = PoliticalDemocracy)

Note that the model defined above estimates all parameters time-point specific. That is, no measurement invariance is assumed.

Now, we want to redefine the parameters as outlined above:

\[ \begin{align} a_2 &= a_1 + \Delta a_2\\ b_2 &= b_1 + \Delta b_2\\ c_2 &= c_1 + \Delta c_2 \end{align} \]

In lessSEM such redefinitions are called transformations and can be passed to the penalty functions (e.g., to lasso) using the modifyModel command.

First, we have to create a definition of our transformations:

transformations <- "
// IMPORTANT: Our transformations always have to start with the follwing line:
parameters: a1, a2, b1, b2, c1, c2, delta_a2, delta_b2, delta_c2

// In the line above, we defined the names of the parameters which we
// want to use in our transformations. EACH AND EVERY PARAMETER USED IN
// THE FOLLOWING MUST BE STATED ABOVE. The line must always start with
// the keyword 'parameters' followed by a colon. The parameters must be
// separated by commata.
// Comments can be added by using double backslash as shown here.

// Now we can state our transformations:

a2 = a1 + delta_a2; // Note: Each declaration must end with a semi-colon!
b2 = b1 + delta_b2;
c2 = c1 + delta_c2;
"

Next, we have to pass the transformations variable to the penalty function:

lassoFit <- lasso(lavaanModel = lavaanFit, 
                  regularized = c("delta_a2", "delta_b2", "delta_c2"),# we want to regularize 
                  # the differences between the parameters
                  nLambdas = 100,
                  # Our model modification must make use of the modifyModel - function:
                  modifyModel = modifyModel(transformations = transformations)
)

Let’s have a look at the parameter estimates:

coef(lassoFit)@estimates[seq(1,100,10),c("a1", "b1", "c1", "delta_a2", "delta_b2", "delta_c2")]
#>             a1       b1       c1     delta_a2    delta_b2 delta_c2
#>  [1,] 1.210993 1.167936 1.234011  0.000000000 0.000000000        0
#>  [2,] 1.211114 1.166390 1.234046  0.000000000 0.002702789        0
#>  [3,] 1.212584 1.150753 1.234650  0.000000000 0.030080792        0
#>  [4,] 1.214347 1.135698 1.235563  0.000000000 0.057049529        0
#>  [5,] 1.216419 1.121154 1.236776  0.000000000 0.083716111        0
#>  [6,] 1.218762 1.107040 1.238259  0.000000000 0.110187427        0
#>  [7,] 1.221414 1.093407 1.240046  0.000000000 0.136498978        0
#>  [8,] 1.226629 1.080063 1.242051 -0.003960392 0.162516933        0
#>  [9,] 1.246435 1.067714 1.244329 -0.032144145 0.186302860        0
#> [10,] 1.266524 1.055803 1.247136 -0.060258715 0.210141745        0

Note that the differences between the parameters get smaller with larger \(\lambda\) values. We can also plot the differences:

plot(lassoFit)
plot of chunk unnamed-chunk-9
plot of chunk unnamed-chunk-9

To check if measurement invariance can be assumed, we can select the best model using information criteria:

coef(lassoFit, criterion = "BIC")
#>                                                                                                                              
#>   Tuning         ||--||  Estimates                                                                                           
#>  ------- ------- ||--|| ---------- ---------- ---------- ---------- ---------- ----------- ----------- ----------- ----------
#>   lambda   alpha ||--||  ind60=~x2  ind60=~x3         a1         b1         c1 dem60~ind60 dem65~ind60 dem65~dem60     y1~~y5
#>  ======= ======= ||--|| ========== ========== ========== ========== ========== =========== =========== =========== ==========
#>   0.2216  1.0000 ||--||     2.1825     1.8189     1.2110     1.1679     1.2340      1.4534      0.5935      0.8659     0.5552
#>                                                                                                                          
#>                                                                                                                          
#>  ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
#>      y2~~y4     y3~~y7     y4~~y8     y6~~y8     x1~~x1     x2~~x2     x3~~x3     y1~~y1     y2~~y2     y3~~y3     y4~~y4
#>  ========== ========== ========== ========== ========== ========== ========== ========== ========== ========== ==========
#>      1.5947     0.7807     0.6537     1.5350     0.0820     0.1177     0.4675     1.7929     7.3843     5.0175     3.4074
#>                                                                                                                            
#>                                                                                                                      ||--||
#>  ---------- ---------- ---------- ---------- ------------ ------------ ------------ ---------- ---------- ---------- ||--||
#>      y5~~y5     y6~~y6     y7~~y7     y8~~y8 ind60~~ind60 dem60~~dem60 dem65~~dem65   delta_a2   delta_b2   delta_c2 ||--||
#>  ========== ========== ========== ========== ============ ============ ============ ========== ========== ========== ||--||
#>      2.2857     4.8977     3.5510     3.4511       0.4480       3.9408       0.2034          .          .          . ||--||
#>                                  
#>   Transform                      
#>  ---------- ---------- ----------
#>          a2         b2         c2
#>  ========== ========== ==========
#>      1.2110     1.1679     1.2340

Note that all differences have been zeroed – that is, a model with full measurement invariance did fit best.

We can also access the transformed parameters:

head(lassoFit@transformations)
#>      lambda alpha       a2       b2       c2
#> 1 0.2410976     1 1.210993 1.167936 1.234011
#> 2 0.2386623     1 1.210996 1.167933 1.234011
#> 3 0.2362270     1 1.210997 1.167935 1.234010
#> 4 0.2337916     1 1.210992 1.167933 1.234006
#> 5 0.2313563     1 1.210991 1.167933 1.234001
#> 6 0.2289210     1 1.210987 1.167932 1.234000

Limitation: Above, we did not take into account that the variables may have different scales; a thorough use of the method should scale the data first.

Some Guidelines

When using transformations, please make sure to give your parameters names which are compatible with standard naming conventions in R. The default names of lavaan (e.g., f=~y1 for loadings) are not supported. That is, all parameters used in the transformations should be given names in the lavaan syntax.

In the example above, we used the following syntax:

modelSyntax <- ' 
  # latent variable definitions
     ind60 =~ x1 + x2 + x3
     dem60 =~ y1 + a1*y2 + b1*y3 + c1*y4
     dem65 =~ y5 + a2*y6 + b2*y7 + c2*y8

  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60

  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'

Importantly, all parameters used in the transformation (a1, b1, c1, a2, b2, c2) were labeled in the lavaan syntax. As a counter example:

modelSyntax <- ' 
  # latent variable definitions
     ind60 =~ x1 + x2 + x3
     dem60 =~ y1 + y2 + y3 + y4
     dem65 =~ y5 + y6 + y7 + y8

  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60

  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'

The syntax above specifies the same model, but will use the lavaan-specific naming convention for the parameters. a1, for example, will be named dem60=~y2. These names are not compatible with the current implementation of transformations used in lessSEM.

Further Examples

Another example where the transformations could be useful is when detecting non-stationarity in autoregressive and cross-lagged parameters (e.g., Liang et al., 2018). In the following, we will demonstrate this with an autoregressive model. The model is defined as:

\[ \begin{align} \eta_t &= a_t\eta_{t-1} + \zeta_t\\ \begin{bmatrix} y_{1,t}\\ y_{2,t}\\ y_{3,t}\\ \end{bmatrix} &= \begin{bmatrix} l_1\\ l_2\\ l_3\\ \end{bmatrix} \eta_t + \pmb\varepsilon \end{align} \]

It is often assumed that the autoregressive effect \(a_t\) is constant over time; that is, the same autoregressive effect is used for all time points. This is a strong assumption and we may want to test it. One way to do so is by using the same procedure outlined above, where we define \(a_t = a_1 + \Delta a_t\) (see Jacobucci & Grimm, 2018 for a similar procedure in latent change score models). In this case, each autoregressive effect is composed of the first autoregressive effect (\(a_1\)) and the difference between the parameters (\(\Delta a_t\)). By regularizing \(\Delta a_t\), we can enforce stationarity.

A drawback of the approach outlined above is that the first autoregressive effect is treated differently from the rest: After all, why should \(a_1\) serve as baseline and not \(a_2\) or \(a_5\)? We will take a slightly different approach that is basically identical to the fused lasso proposed by Tibshirani et al. (2005). Let’s define the autoregressive effect as

\[a_t = a_{t-1} + \Delta a_t\]

Note that \(\Delta a_t\) is no longer the difference with respect to the initial autoregressive effect \(a_1\) but the difference with respect to the directly preceding time point. When regularizing \(\Delta a_t\), we can now detect sudden changes in the parameter – e.g., due to an intervention. This can also be thought of as a regime switching model, where the underlying model changes over time (see also Ou et al., 2019 for regime switching models). With our regularization procedure, we want to detect if and when the process changes.

We won’t go into the details of how to set up the model here, but you can find them in the source of this file (e.g., in GitHub). We simulated a data set with 200 individuals measured at 10 time points. The autoregressive effect \(a_t\) changes at \(t=4\) from \(.6\) to \(.2\).

The data looks as follows:

head(data)
#>           y1_t1       y2_t1       y3_t1      y1_t2      y2_t2      y3_t2      y1_t3       y2_t3      y3_t3       y1_t4
#> [1,]  0.2627294 -0.07745288 -0.06014729  1.1052896 -0.5125258  0.1352253  0.6775861  0.48382200  0.5163607  0.06936298
#> [2,] -0.5795123  0.27989039 -1.76295329 -1.1646535  0.0607009 -0.7770722 -0.4571539 -0.30575453  0.7552239  0.06618071
#> [3,]  0.5746526  0.67203557  1.30772224  0.5620735  1.3873269  1.2415508  0.8255529 -0.09010965  1.5877450  0.75756487
#> [4,]  1.6376973  1.30880439  1.35159277  0.4681794  1.5370870  0.6961462  0.1872212 -1.57377418  1.1081811 -0.50832183
#> [5,]  0.4995286 -0.23722728  0.61244148  0.3469646  0.3052115  0.2639443  1.3768225  0.02394048  1.1573912  1.30350064
#> [6,] -0.9766373  0.94700260 -0.59421930 -0.3197840  0.8221090 -0.6808501 -0.3434095 -0.90348981 -1.1760506 -1.72146444
#>           y2_t4      y3_t4      y1_t5      y2_t5      y3_t5      y1_t6      y2_t6       y3_t6      y1_t7      y2_t7      y3_t7
#> [1,]  1.6842420  0.5604284 -1.6171372 -1.5100658 -2.1212448 -0.9014921 -0.9732405 -0.89209997 -1.4697023 -0.9088437 -2.1031570
#> [2,] -1.5803092  0.1672394  0.1045129  0.3483139 -0.3669688 -0.6347357 -1.1095700  0.01821975  0.3845160  1.0140172  1.8585956
#> [3,]  1.3579232 -0.1379793 -1.1968405 -0.3623344 -0.8732056 -0.8290071 -0.3555817 -0.98268524 -1.8270782 -0.2335533 -1.6441484
#> [4,]  1.1968847 -1.6473363  0.4423361 -1.0372594 -0.8153181 -1.6003121 -2.1169045 -0.46541710  1.3226249  1.7154104  1.8113749
#> [5,]  0.3049751 -0.0146624  1.3351806  0.1229075  0.8274846  0.4508487  0.1328104  0.65883759 -0.3915346 -0.5549143 -0.2476716
#> [6,] -1.2828716 -0.4379788  1.0261730  0.1686185  0.6083559 -0.7427341  0.4514472 -1.15435404  1.7453347  1.2488910  0.6649188
#>           y1_t8       y2_t8      y3_t8       y1_t9      y2_t9      y3_t9      y1_t10      y2_t10     y3_t10
#> [1,]  0.9915954  0.26021125  1.2988304  1.27416201  1.2010039  1.1426722  0.04668192 -0.24123204 -0.4421784
#> [2,] -1.1280623  0.62646559 -0.8323908 -0.41574193 -0.6890046 -1.4714503 -0.01788118  0.02058262  0.3191417
#> [3,] -0.2693417  0.09757043  0.5752644  0.37098052 -0.3486185 -0.4273464 -0.87844796 -1.36443913 -1.0240934
#> [4,]  0.3459617  0.29614603  0.2463911  0.55566216  1.2444869  0.1788493  0.81351163 -0.39685221  3.3330365
#> [5,] -1.3659303  1.14721685 -0.1968243 -1.58477183 -1.8546910  0.7432728  0.89776939  1.05624257  0.3635211
#> [6,] -0.6131110 -0.27533871  0.1895920  0.02486185 -0.5151950 -0.2018260  1.13097234  0.74034559  2.5292256

The lavaan model is defined as follows:

#> eta2 ~ a1*eta1
#> eta3 ~ a2*eta2
#> eta4 ~ a3*eta3
#> eta5 ~ a4*eta4
#> eta6 ~ a5*eta5
#> eta7 ~ a6*eta6
#> eta8 ~ a7*eta7
#> eta9 ~ a8*eta8
#> eta10 ~ a9*eta9
#> 
#> 
#> eta1 ~~ eta1
#> eta2 ~~ v*eta2
#> eta3 ~~ v*eta3
#> eta4 ~~ v*eta4
#> eta5 ~~ v*eta5
#> eta6 ~~ v*eta6
#> eta7 ~~ v*eta7
#> eta8 ~~ v*eta8
#> eta9 ~~ v*eta9
#> eta10 ~~ v*eta10
#> 
#> eta1 =~ 1*y1_t1 + l2*y2_t1 + l3*y3_t1
#> y1_t1 ~~ mvar1*y1_t1
#> y2_t1 ~~ mvar2*y2_t1
#> y3_t1 ~~ mvar3*y3_t1
#> eta2 =~ 1*y1_t2 + l2*y2_t2 + l3*y3_t2
#> y1_t2 ~~ mvar1*y1_t2
#> y2_t2 ~~ mvar2*y2_t2
#> y3_t2 ~~ mvar3*y3_t2
#> eta3 =~ 1*y1_t3 + l2*y2_t3 + l3*y3_t3
#> y1_t3 ~~ mvar1*y1_t3
#> y2_t3 ~~ mvar2*y2_t3
#> y3_t3 ~~ mvar3*y3_t3
#> eta4 =~ 1*y1_t4 + l2*y2_t4 + l3*y3_t4
#> y1_t4 ~~ mvar1*y1_t4
#> y2_t4 ~~ mvar2*y2_t4
#> y3_t4 ~~ mvar3*y3_t4
#> eta5 =~ 1*y1_t5 + l2*y2_t5 + l3*y3_t5
#> y1_t5 ~~ mvar1*y1_t5
#> y2_t5 ~~ mvar2*y2_t5
#> y3_t5 ~~ mvar3*y3_t5
#> eta6 =~ 1*y1_t6 + l2*y2_t6 + l3*y3_t6
#> y1_t6 ~~ mvar1*y1_t6
#> y2_t6 ~~ mvar2*y2_t6
#> y3_t6 ~~ mvar3*y3_t6
#> eta7 =~ 1*y1_t7 + l2*y2_t7 + l3*y3_t7
#> y1_t7 ~~ mvar1*y1_t7
#> y2_t7 ~~ mvar2*y2_t7
#> y3_t7 ~~ mvar3*y3_t7
#> eta8 =~ 1*y1_t8 + l2*y2_t8 + l3*y3_t8
#> y1_t8 ~~ mvar1*y1_t8
#> y2_t8 ~~ mvar2*y2_t8
#> y3_t8 ~~ mvar3*y3_t8
#> eta9 =~ 1*y1_t9 + l2*y2_t9 + l3*y3_t9
#> y1_t9 ~~ mvar1*y1_t9
#> y2_t9 ~~ mvar2*y2_t9
#> y3_t9 ~~ mvar3*y3_t9
#> eta10 =~ 1*y1_t10 + l2*y2_t10 + l3*y3_t10
#> y1_t10 ~~ mvar1*y1_t10
#> y2_t10 ~~ mvar2*y2_t10
#> y3_t10 ~~ mvar3*y3_t10

We fit the model using lavaan:

lavaanFit <- sem(model = lavaanSyntax, 
                 data = data,
                 orthogonal.y = TRUE, 
                 orthogonal.x = TRUE,
                 missing = "ml")
coef(lavaanFit)
#>         a1         a2         a3         a4         a5         a6         a7         a8         a9 eta1~~eta1          v 
#>      0.484      0.535      0.531      0.135      0.124      0.040      0.173      0.296      0.185      1.108      0.804 
#>          v          v          v          v          v          v          v          v         l2         l3      mvar1 
#>      0.804      0.804      0.804      0.804      0.804      0.804      0.804      0.804      0.567      0.672      0.049 
#>      mvar2      mvar3         l2         l3      mvar1      mvar2      mvar3         l2         l3      mvar1      mvar2 
#>      0.744      0.652      0.567      0.672      0.049      0.744      0.652      0.567      0.672      0.049      0.744 
#>      mvar3         l2         l3      mvar1      mvar2      mvar3         l2         l3      mvar1      mvar2      mvar3 
#>      0.652      0.567      0.672      0.049      0.744      0.652      0.567      0.672      0.049      0.744      0.652 
#>         l2         l3      mvar1      mvar2      mvar3         l2         l3      mvar1      mvar2      mvar3         l2 
#>      0.567      0.672      0.049      0.744      0.652      0.567      0.672      0.049      0.744      0.652      0.567 
#>         l3      mvar1      mvar2      mvar3         l2         l3      mvar1      mvar2      mvar3         l2         l3 
#>      0.672      0.049      0.744      0.652      0.567      0.672      0.049      0.744      0.652      0.567      0.672 
#>      mvar1      mvar2      mvar3    y1_t1~1    y2_t1~1    y3_t1~1    y1_t2~1    y2_t2~1    y3_t2~1    y1_t3~1    y2_t3~1 
#>      0.049      0.744      0.652     -0.155     -0.088     -0.175     -0.157     -0.049     -0.085     -0.150     -0.030 
#>    y3_t3~1    y1_t4~1    y2_t4~1    y3_t4~1    y1_t5~1    y2_t5~1    y3_t5~1    y1_t6~1    y2_t6~1    y3_t6~1    y1_t7~1 
#>     -0.157     -0.095     -0.022     -0.080      0.037     -0.023      0.058      0.050      0.065     -0.002     -0.047 
#>    y2_t7~1    y3_t7~1    y1_t8~1    y2_t8~1    y3_t8~1    y1_t9~1    y2_t9~1    y3_t9~1   y1_t10~1   y2_t10~1   y3_t10~1 
#>      0.068     -0.038     -0.022     -0.028     -0.048     -0.077     -0.058     -0.041      0.013     -0.012     -0.004

Note that no constraints on autoregressive effects are implemented – each effect (a1-a9) is estimated separately.

We now define transformations as follows:

#> parameters: a1, a2, a3, a4, a5, a6, a7, a8, a9, delta2, delta3, delta4, delta5, delta6, delta7, delta8, delta9
#> 
#> a2 = a1 + delta2;
#> a3 = a2 + delta3;
#> a4 = a3 + delta4;
#> a5 = a4 + delta5;
#> a6 = a5 + delta6;
#> a7 = a6 + delta7;
#> a8 = a7 + delta8;
#> a9 = a8 + delta9;

Finally, we can fit our model:

lassoFit <- lasso(lavaanModel = lavaanFit, 
                  regularized = paste0("delta", 2:9),# we want to regularize 
                  # the differences between the parameters
                  nLambdas = 100,
                  # glmnet is considerably faster here:
                  method = "glmnet",
                  control = controlGlmnet(),
                  # Our model modification must make use of the modifyModel - function:
                  modifyModel = modifyModel(transformations = transformations)
)

Extracting the best fitting model:

coef(lassoFit, criterion = "BIC")
#>                                                                                                                           
#>   Tuning         ||--||  Estimates                                                                                        
#>  ------- ------- ||--|| ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
#>   lambda   alpha ||--||         a1 eta1~~eta1          v         l2         l3      mvar1      mvar2      mvar3    y1_t1~1
#>  ======= ======= ||--|| ========== ========== ========== ========== ========== ========== ========== ========== ==========
#>   0.2790  1.0000 ||--||     0.4665     1.1055     0.8051     0.5696     0.6754     0.0533     0.7428     0.6500    -0.1545
#>                                                                                                                          
#>                                                                                                                          
#>  ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
#>     y2_t1~1    y3_t1~1    y1_t2~1    y2_t2~1    y3_t2~1    y1_t3~1    y2_t3~1    y3_t3~1    y1_t4~1    y2_t4~1    y3_t4~1
#>  ========== ========== ========== ========== ========== ========== ========== ========== ========== ========== ==========
#>     -0.0878    -0.1745    -0.1565    -0.0492    -0.0849    -0.1499    -0.0302    -0.1572    -0.0952    -0.0224    -0.0798
#>                                                                                                                          
#>                                                                                                                          
#>  ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
#>     y1_t5~1    y2_t5~1    y3_t5~1    y1_t6~1    y2_t6~1    y3_t6~1    y1_t7~1    y2_t7~1    y3_t7~1    y1_t8~1    y2_t8~1
#>  ========== ========== ========== ========== ========== ========== ========== ========== ========== ========== ==========
#>      0.0374    -0.0231     0.0576     0.0502     0.0651    -0.0020    -0.0466     0.0678    -0.0383    -0.0221    -0.0279
#>                                                                                                                          
#>                                                                                                                          
#>  ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
#>     y3_t8~1    y1_t9~1    y2_t9~1    y3_t9~1   y1_t10~1   y2_t10~1   y3_t10~1     delta2     delta3     delta4     delta5
#>  ========== ========== ========== ========== ========== ========== ========== ========== ========== ========== ==========
#>     -0.0476    -0.0775    -0.0579    -0.0407     0.0131    -0.0118    -0.0043          .          .    -0.2807          .
#>                                                                                                                      
#>                                              ||--||  Transform                                                       
#>  ---------- ---------- ---------- ---------- ||--|| ---------- ---------- ---------- ---------- ---------- ----------
#>      delta6     delta7     delta8     delta9 ||--||         a2         a3         a4         a5         a6         a7
#>  ========== ========== ========== ========== ||--|| ========== ========== ========== ========== ========== ==========
#>           .          .          .          . ||--||     0.4665     0.4665     0.1859     0.1859     0.1859     0.1859
#>                       
#>                       
#>  ---------- ----------
#>          a8         a9
#>  ========== ==========
#>      0.1859     0.1859

The true autoregressive effects are given by

#> [1] 0.6 0.6 0.6 0.2 0.2 0.2 0.2 0.2 0.2

while the estimates are

#> [1] 0.4665209 0.4665209 0.4665209 0.1858669 0.1858669 0.1858669 0.1858669 0.1858669 0.1858669

The result is not perfect, but lessSEM correctly identified a change in the autoregressive parameter.

Looking under the hood

The transformations used above are implemented using RcppArmadillo. This allows for a lot more complicated transformations than those outlined before. In general, lessSEM will take your transformations and try to translate them to C++. Let’s assume that the model is given by our first example:

modelSyntax <- ' 
  # latent variable definitions
     ind60 =~ x1 + x2 + x3
     dem60 =~ y1 + a1*y2 + b1*y3 + c1*y4
     dem65 =~ y5 + a2*y6 + b2*y7 + c2*y8

  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60

  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'

lavaanFit <- sem(model = modelSyntax, 
                 data = PoliticalDemocracy)

We defined the transformations to be:

transformations <- "
parameters: a1, a2, b1, b2, c1, c2, delta_a2, delta_b2, delta_c2

a2 = a1 + delta_a2;
b2 = b1 + delta_b2;
c2 = c1 + delta_c2;
"

When this transformation is passed to lessSEM, lessSEM will first try to figure out which parameters are already in the model and which ones are new. In our case a1, a2, b1, b2, c1, and c2 are already known, while delta_a2, delta_b2, and delta_c2 are new. lessSEM will now add the new parameters to the internal parameter vector. Next, lessSEM will scan the names of the parameters which appear on the left hand side of an equation (a2, b2, and c2) in our case. This will tell lessSEM which of your parameters are functions of other parameters (i.e., transformations). Knowing that a2 is a function of other parameters will tell lessSEM, that a2 should no longer be estimated. Instead, the parameters which make up a2 are estimated: a1 and delta_a2.

To see this in action, we can create the C++ function without compilation:

transformationFunction <- lessSEM:::.compileTransformations(syntax = transformations, 
                                                            parameterLabels = names(getLavaanParameters(lavaanFit)),
                                                            compile = FALSE)

First, let’s have a look at the extended parameter vector:

transformationFunction$parameters
#> [1] "a1"       "a2"       "b1"       "b2"       "c1"       "c2"       "delta_a2" "delta_b2" "delta_c2"

Note that delta_a2, delta_b2, and delta_c2 have been added. Some of these are transformations of other parameters:

transformationFunction$isTransformation
#> [1] "a2" "b2" "c2"

These will not be estimated but computed based on the other model parameters.

Finally, the C++ function has been returned:

cat(transformationFunction$armaFunction)
#> 
#>   // [[Rcpp::depends(RcppArmadillo)]]
#>   #include <RcppArmadillo.h>
#>   // [[Rcpp::export]]
#>   Rcpp::NumericVector transformationFunction(Rcpp::NumericVector& parameterValues, Rcpp::List transformationList)
#>   {
#>   using namespace Rcpp;
#>   using namespace arma;
#>   
#>   // extract required parameters from parameterValues
#>   
#> double a1 = parameterValues["a1"];
#> double a2 = parameterValues["a2"];
#> double b1 = parameterValues["b1"];
#> double b2 = parameterValues["b2"];
#> double c1 = parameterValues["c1"];
#> double c2 = parameterValues["c2"];
#> double delta_a2 = parameterValues["delta_a2"];
#> double delta_b2 = parameterValues["delta_b2"];
#> double delta_c2 = parameterValues["delta_c2"];
#> 
#> 
#> // add user defined functions
#> 
#> 
#> a2 = a1 + delta_a2;
#> b2 = b1 + delta_b2;
#> c2 = c1 + delta_c2;
#> 
#> 
#> 
#> // update parameters
#> parameterValues["a1"] = a1;
#> parameterValues["a2"] = a2;
#> parameterValues["b1"] = b1;
#> parameterValues["b2"] = b2;
#> parameterValues["c1"] = c1;
#> parameterValues["c2"] = c2;
#> parameterValues["delta_a2"] = delta_a2;
#> parameterValues["delta_b2"] = delta_b2;
#> parameterValues["delta_c2"] = delta_c2;
#> 
#>   
#>   return(parameterValues);
#>   }
#> 
#>   
#>   
#>   // Dirk Eddelbuettel at
#>   // https://gallery.rcpp.org/articles/passing-cpp-function-pointers/
#> typedef Rcpp::NumericVector (*transformationFunctionPtr)(Rcpp::NumericVector&, //parameters
#> Rcpp::List // transformationList
#> );
#> 
#> typedef Rcpp::XPtr<transformationFunctionPtr> transformationFunctionPtr_t;
#> 
#> // [[Rcpp::export]]
#> transformationFunctionPtr_t getPtr() {
#>         return(transformationFunctionPtr_t(new transformationFunctionPtr(&transformationFunction)));
#> }

Most importantly, note that the first step here is to extract the required parameters from the parameter vector (e.g., double a1 = parameterValues["a1"];). Next, all of these parameters are directly available for use in your transformations. This is why we can simply write a2 = a1 + delta_a2;. Finally, the transformed parameters are returned. To pass our function to lessSEM, we also create a pointer to our function, but this is beyond the scope here.

At this point you may be wondering where all the more complicated transformations are that we promised above. Importantly, you can use any of the functions implemented in Rcpp or RcppArmadillo which can be applied to variables of type double within your transformations without any in-depth knowledge of C++. For instance, RcppArmadillo comes with an exponential-function (exp), a pow and a log function. Making use of this, we can implement a univariate continuous time SEM (e.g., Voelkle & Oud, 2012). Far superior versions of this model are implemented in ctsem and dynr) and Arnold et al. (in submission) recently derived close form solutions for the gradients of these models which should outperform lessSEM considerably in terms of runtime.

We will use the same model from above, but remove the change in the autoregressive effect. The code to simulate the data set can be found in the source of this file (e.g., on GitHub).

The initial model is the same as before, however the autoregressive effect is constrained to equality over time and so are the manifest means. We also added an initial mean for latent variable \(\eta\) and changed the names of some variables to make using ctsem with this data easier:

cat(lavaanSyntax)
#> eta1 ~ a*eta0
#> eta2 ~ a*eta1
#> eta3 ~ a*eta2
#> eta4 ~ a*eta3
#> eta5 ~ a*eta4
#> eta6 ~ a*eta5
#> eta7 ~ a*eta6
#> eta8 ~ a*eta7
#> eta9 ~ a*eta8
#> 
#> 
#> eta0 ~~ eta0
#> eta1 ~~ v*eta1
#> eta2 ~~ v*eta2
#> eta3 ~~ v*eta3
#> eta4 ~~ v*eta4
#> eta5 ~~ v*eta5
#> eta6 ~~ v*eta6
#> eta7 ~~ v*eta7
#> eta8 ~~ v*eta8
#> eta9 ~~ v*eta9
#> 
#> 
#> eta0~1
#> 
#> eta0 =~ 1*y1_T0 + l2*y2_T0 + l3*y3_T0
#> y1_T0 ~~ mvar1*y1_T0
#> y2_T0 ~~ mvar2*y2_T0
#> y3_T0 ~~ mvar3*y3_T0
#> y1_T0 ~ mMean1*1
#> y2_T0 ~ mMean2*1
#> y3_T0 ~ mMean3*1
#> eta1 =~ 1*y1_T1 + l2*y2_T1 + l3*y3_T1
#> y1_T1 ~~ mvar1*y1_T1
#> y2_T1 ~~ mvar2*y2_T1
#> y3_T1 ~~ mvar3*y3_T1
#> y1_T1 ~ mMean1*1
#> y2_T1 ~ mMean2*1
#> y3_T1 ~ mMean3*1
#> eta2 =~ 1*y1_T2 + l2*y2_T2 + l3*y3_T2
#> y1_T2 ~~ mvar1*y1_T2
#> y2_T2 ~~ mvar2*y2_T2
#> y3_T2 ~~ mvar3*y3_T2
#> y1_T2 ~ mMean1*1
#> y2_T2 ~ mMean2*1
#> y3_T2 ~ mMean3*1
#> eta3 =~ 1*y1_T3 + l2*y2_T3 + l3*y3_T3
#> y1_T3 ~~ mvar1*y1_T3
#> y2_T3 ~~ mvar2*y2_T3
#> y3_T3 ~~ mvar3*y3_T3
#> y1_T3 ~ mMean1*1
#> y2_T3 ~ mMean2*1
#> y3_T3 ~ mMean3*1
#> eta4 =~ 1*y1_T4 + l2*y2_T4 + l3*y3_T4
#> y1_T4 ~~ mvar1*y1_T4
#> y2_T4 ~~ mvar2*y2_T4
#> y3_T4 ~~ mvar3*y3_T4
#> y1_T4 ~ mMean1*1
#> y2_T4 ~ mMean2*1
#> y3_T4 ~ mMean3*1
#> eta5 =~ 1*y1_T5 + l2*y2_T5 + l3*y3_T5
#> y1_T5 ~~ mvar1*y1_T5
#> y2_T5 ~~ mvar2*y2_T5
#> y3_T5 ~~ mvar3*y3_T5
#> y1_T5 ~ mMean1*1
#> y2_T5 ~ mMean2*1
#> y3_T5 ~ mMean3*1
#> eta6 =~ 1*y1_T6 + l2*y2_T6 + l3*y3_T6
#> y1_T6 ~~ mvar1*y1_T6
#> y2_T6 ~~ mvar2*y2_T6
#> y3_T6 ~~ mvar3*y3_T6
#> y1_T6 ~ mMean1*1
#> y2_T6 ~ mMean2*1
#> y3_T6 ~ mMean3*1
#> eta7 =~ 1*y1_T7 + l2*y2_T7 + l3*y3_T7
#> y1_T7 ~~ mvar1*y1_T7
#> y2_T7 ~~ mvar2*y2_T7
#> y3_T7 ~~ mvar3*y3_T7
#> y1_T7 ~ mMean1*1
#> y2_T7 ~ mMean2*1
#> y3_T7 ~ mMean3*1
#> eta8 =~ 1*y1_T8 + l2*y2_T8 + l3*y3_T8
#> y1_T8 ~~ mvar1*y1_T8
#> y2_T8 ~~ mvar2*y2_T8
#> y3_T8 ~~ mvar3*y3_T8
#> y1_T8 ~ mMean1*1
#> y2_T8 ~ mMean2*1
#> y3_T8 ~ mMean3*1
#> eta9 =~ 1*y1_T9 + l2*y2_T9 + l3*y3_T9
#> y1_T9 ~~ mvar1*y1_T9
#> y2_T9 ~~ mvar2*y2_T9
#> y3_T9 ~~ mvar3*y3_T9
#> y1_T9 ~ mMean1*1
#> y2_T9 ~ mMean2*1
#> y3_T9 ~ mMean3*1
lavaanFit <- sem(model = lavaanSyntax, 
                 data = data,
                 orthogonal.y = TRUE, 
                 orthogonal.x = TRUE,
                 missing = "ml")
getLavaanParameters(lavaanFit)
#>            a   eta0~~eta0            v       eta0~1           l2           l3        mvar1        mvar2        mvar3 
#>  0.592770935  0.873870399  0.387065602 -0.091482670  0.578004929  0.745298413  0.127999058  0.750360430  0.577566347 
#>       mMean1       mMean2       mMean3 
#> -0.026015718  0.005858067 -0.027032832

Now, we define the transformations for the latent variables to turn our model in a continuous time SEM:

transformations <- "
parameters: a, ctA, v, ctV
// NOTE: We can define starting values for our parameters. This
// is implemented with the 'start:' keyword:
start: ctA = -.1, ctV = .1

// We changed the starting values for the ct parameters
// because the auto-effect ctA should be negative.

a = exp(ctA);
v = log((1.0/(2.0*ctA))*(exp(2.0*ctA)-1)*pow(ctV,2.0)); // we take
// the log because lessSEM internally takes the exponential of
// any variance parameter (v in our case) to avoid negative variances.
"
lessSEMFit <- bfgs(lavaanModel = lavaanFit, 
                   # Our model modification must make use of the modifyModel - function:
                   modifyModel = modifyModel(transformations = transformations)
)

Let’s have a look at the parameter estimates:

coef(lessSEMFit)
#>                                                                                                                           
#>   Tuning         ||--||  Estimates                                                                                        
#>  ------- ------- ||--|| ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
#>   lambda   alpha ||--|| eta0~~eta0     eta0~1         l2         l3      mvar1      mvar2      mvar3     mMean1     mMean2
#>  ======= ======= ||--|| ========== ========== ========== ========== ========== ========== ========== ========== ==========
#>   0.0000  0.0000 ||--||     0.8736    -0.0902     0.5780     0.7453     0.1280     0.7504     0.5775    -0.0267     0.0054
#>                                                               
#>                                   ||--||  Transform           
#>  ---------- ---------- ---------- ||--|| ---------- ----------
#>      mMean3        ctA        ctV ||--||          a          v
#>  ========== ========== ========== ||--|| ========== ==========
#>     -0.0276    -0.5228     0.7899 ||--||     0.5928     0.3870

For comparison, we will run the same model with ctsem:

library(ctsemOMX)
dataCt <- cbind(data,
                data.frame("dT1" = rep(1,nrow(data)),
                           "dT2" = rep(1,nrow(data)),
                           "dT3" = rep(1,nrow(data)),
                           "dT4" = rep(1,nrow(data)),
                           "dT5" = rep(1,nrow(data)),
                           "dT6" = rep(1,nrow(data)),
                           "dT7" = rep(1,nrow(data)),
                           "dT8" = rep(1,nrow(data)),
                           "dT9" = rep(1,nrow(data))))
cModel <- ctModel(type = "omx", 
                  n.manifest = 3, 
                  n.latent = 1, 
                  Tpoints = 10,
                  manifestNames = c("y1","y2", "y3"), 
                  latentNames = "eta",
                  LAMBDA = matrix(c(1,
                                    "l2",
                                    "l3"),3,1,TRUE), 
                  DRIFT = matrix("a",1,1)
)

cFit <- ctFit(dat = dataCt, ctmodelobj = cModel)
ctSummary <- summary(cFit)

The parameter ctA in our model corresponds to the DRIFT parameter in the ctsem summary and the parameter ctV corresponds to the root of the DIFFUSION parameter in the ctsem summary:

coef(lessSEMFit)@estimates[,c("ctA", "ctV")]
#>        ctA        ctV 
#> -0.5228340  0.7899115

# drift value from ctsem:
ctSummary$DRIFT
#>            eta
#> eta -0.5229469
# sqrt(diffusion) value from ctsem:
sqrt(ctSummary$DIFFUSION)
#>           eta
#> eta 0.7900236

Making use of C++

In the example above, we used a univariate ctsem. Because of this, all our functions were fairly simple in that we only needed the log, exp, and pow functions for single variables. However, because lessSEM creates a C++ function, we can build much more powerful transformations if we are familiar with RcppArmadillo. In the following, we will therefore extend the example to a multivariate continuous time SEM. We will use the AnomAuth data set from ctsem. The data is included in the ctsemOMX package:

data("AnomAuth")
head(AnomAuth)
#>   Y1_T0 Y2_T0 Y1_T1 Y2_T1 Y1_T2 Y2_T2 Y1_T3 Y2_T3 Y1_T4 Y2_T4 dT1 dT2 dT3 dT4
#> 1  2.67  3.50  3.33   3.5    NA    NA    NA    NA    NA    NA   1   1   2   2
#> 2  3.33  3.25    NA    NA    NA    NA    NA    NA    NA    NA   1   1   2   2
#> 3  3.33  2.75  3.33   3.0  3.33   2.5  2.33     3  2.33     3   1   1   2   2
#> 4  3.33  3.25    NA    NA    NA    NA    NA    NA    NA    NA   1   1   2   2
#> 5  4.00  4.00    NA    NA    NA    NA    NA    NA    NA    NA   1   1   2   2
#> 6  3.67  4.00    NA    NA    NA    NA  4.00     4  4.00     4   1   1   2   2

The five measurement occasions are unequally spaced. In a discrete time model, we could take care of this by implementing a model with different autoregressive and cross-lagged effects for the different time intervals:

# initial time point
lavaanSyntax <- 
  "eta1_T0 =~ 1 * Y1_T0
eta2_T0 =~ 1 * Y2_T0
Y1_T0 ~~ 0*Y1_T0
Y2_T0 ~~ 0*Y2_T0\n"

# variances
lavaanSyntax <- c(lavaanSyntax,
                  "eta1_T0 ~~ v0_11 * eta1_T0 + v0_12 * eta2_T0\neta2_T0 ~~ v0_22 * eta2_T0\n"
)

# means
lavaanSyntax <- c(lavaanSyntax,
                  "eta1_T0 ~ 1\neta2_T0 ~ 1\nY1_T0~mMean1*1\nY2_T0~mMean2*1\n" 
)

for(tp in c(0,1,2,3)){
  if(tp < 2) {
    a <- "a1"
    v <- "v1"
  }else{
    a <-"a2"
    v <- "v2"
  }
  
  # autoregressive and cross-lagged
  lavaanSyntax <- c(lavaanSyntax,
                    paste0(
                      "eta1_T", tp+1, " ~ ", a, "_11 * eta1_T", tp, " + ", a, "_12 * eta2_T", tp,"\n",
                      "eta2_T", tp+1, " ~ ", a, "_21 * eta1_T", tp, " + ", a, "_22 * eta2_T", tp, "\n"
                    )
  )
  
  # variances
  lavaanSyntax <- c(lavaanSyntax,
                    paste0(
                      "eta1_T", tp+1, " ~~ ", v, "_11 * eta1_T", tp+1, " + ", v, "_12 * eta2_T", tp+1,"\n",
                      "eta2_T", tp+1, " ~~ ", v, "_22 * eta2_T", tp+1, "\n"
                    )
  )
  
  # loadings
  lavaanSyntax <- c(lavaanSyntax,
                    paste0(
                      "eta1_T", tp+1, " =~ 1 * Y1_T", tp+1,"\n",
                      "eta2_T", tp+1, " =~ 1 * Y2_T", tp+1,"\n"
                    )
  )
  
  # manifest variances
  lavaanSyntax <- c(lavaanSyntax,
                    paste0(
                      "Y1_T", tp+1, " ~~ 0* Y1_T", tp+1, "\n",
                      "Y2_T", tp+1, " ~~ 0* Y2_T", tp+1, "\n"
                    )
  )
  
  # manifest means
  lavaanSyntax <- c(lavaanSyntax,
                    paste0(
                      "Y1_T", tp+1, " ~ mMean1 * 1\n",
                      "Y2_T", tp+1, " ~ mMean2 * 1\n"
                    )
  )
}
lavaanSyntax <- paste0(lavaanSyntax, collapse = "")
cat(lavaanSyntax)
#> eta1_T0 =~ 1 * Y1_T0
#> eta2_T0 =~ 1 * Y2_T0
#> Y1_T0 ~~ 0*Y1_T0
#> Y2_T0 ~~ 0*Y2_T0
#> eta1_T0 ~~ v0_11 * eta1_T0 + v0_12 * eta2_T0
#> eta2_T0 ~~ v0_22 * eta2_T0
#> eta1_T0 ~ 1
#> eta2_T0 ~ 1
#> Y1_T0~mMean1*1
#> Y2_T0~mMean2*1
#> eta1_T1 ~ a1_11 * eta1_T0 + a1_12 * eta2_T0
#> eta2_T1 ~ a1_21 * eta1_T0 + a1_22 * eta2_T0
#> eta1_T1 ~~ v1_11 * eta1_T1 + v1_12 * eta2_T1
#> eta2_T1 ~~ v1_22 * eta2_T1
#> eta1_T1 =~ 1 * Y1_T1
#> eta2_T1 =~ 1 * Y2_T1
#> Y1_T1 ~~ 0* Y1_T1
#> Y2_T1 ~~ 0* Y2_T1
#> Y1_T1 ~ mMean1 * 1
#> Y2_T1 ~ mMean2 * 1
#> eta1_T2 ~ a1_11 * eta1_T1 + a1_12 * eta2_T1
#> eta2_T2 ~ a1_21 * eta1_T1 + a1_22 * eta2_T1
#> eta1_T2 ~~ v1_11 * eta1_T2 + v1_12 * eta2_T2
#> eta2_T2 ~~ v1_22 * eta2_T2
#> eta1_T2 =~ 1 * Y1_T2
#> eta2_T2 =~ 1 * Y2_T2
#> Y1_T2 ~~ 0* Y1_T2
#> Y2_T2 ~~ 0* Y2_T2
#> Y1_T2 ~ mMean1 * 1
#> Y2_T2 ~ mMean2 * 1
#> eta1_T3 ~ a2_11 * eta1_T2 + a2_12 * eta2_T2
#> eta2_T3 ~ a2_21 * eta1_T2 + a2_22 * eta2_T2
#> eta1_T3 ~~ v2_11 * eta1_T3 + v2_12 * eta2_T3
#> eta2_T3 ~~ v2_22 * eta2_T3
#> eta1_T3 =~ 1 * Y1_T3
#> eta2_T3 =~ 1 * Y2_T3
#> Y1_T3 ~~ 0* Y1_T3
#> Y2_T3 ~~ 0* Y2_T3
#> Y1_T3 ~ mMean1 * 1
#> Y2_T3 ~ mMean2 * 1
#> eta1_T4 ~ a2_11 * eta1_T3 + a2_12 * eta2_T3
#> eta2_T4 ~ a2_21 * eta1_T3 + a2_22 * eta2_T3
#> eta1_T4 ~~ v2_11 * eta1_T4 + v2_12 * eta2_T4
#> eta2_T4 ~~ v2_22 * eta2_T4
#> eta1_T4 =~ 1 * Y1_T4
#> eta2_T4 =~ 1 * Y2_T4
#> Y1_T4 ~~ 0* Y1_T4
#> Y2_T4 ~~ 0* Y2_T4
#> Y1_T4 ~ mMean1 * 1
#> Y2_T4 ~ mMean2 * 1

Setting up the model in lavaan:

lavaanFit <- sem(model = lavaanSyntax, data = AnomAuth,
                 orthogonal.y = TRUE, 
                 orthogonal.x = TRUE,
                 missing = "ml",
                 do.fit = FALSE)

To transform the parameters to those of a continuous time model, we have to define the transformations again. We will not go into the details here. See Voelkle et al. (2012) for the underlying transformations.

transformations <- "
// Define all parameters which we want to use:
parameters: a1_11, a1_12, a1_21, a1_22, a2_11, a2_12, a2_21, a2_22, 
ctA_11, ctA_12, ctA_21, ctA_22, 
v1_11, v1_12, v1_22, v2_11, v2_12, v2_22, 
ctV_11, ctV_12, ctV_22

// Define the starting values for the continuous time parameters:
start: ctA_11 = -1, ctA_12 = 0, ctA_21 = 0, ctA_22 = -1, 
ctV_11 = .1, ctV_12 = 0, ctV_22 = .1

// transformations:
arma::mat drift(2,2);
arma::mat ARCL1(2,2);
arma::mat ARCL2(2,2);
arma::mat driftHash(4,4);
drift(0,0) = ctA_11;
drift(1,0) = ctA_21;
drift(0,1) = ctA_12;
drift(1,1) = ctA_22;
ARCL1 = expmat(drift);
ARCL2 = expmat(drift*2.0);

driftHash = kron(drift, arma::eye(2,2)) + kron(arma::eye(2,2), drift);

arma::mat diffusion(2,2);
arma::mat discreteDiff1(2,2);
arma::mat discreteDiff2(2,2);
diffusion(0,0) = ctV_11;
diffusion(1,0) = ctV_12;
diffusion(0,1) = ctV_12;
diffusion(1,1) = ctV_22;
discreteDiff1 = arma::reshape(arma::inv(driftHash) * 
  (expmat(driftHash) - arma::eye(arma::size(expmat(driftHash))))*
  arma::vectorise(diffusion),2,2);
discreteDiff2 = arma::reshape(arma::inv(driftHash) * 
  (expmat(driftHash*2.0) - arma::eye(arma::size(expmat(driftHash*2.0))))*
  arma::vectorise(diffusion),2,2);

// extract parameters

a1_11 = ARCL1(0,0);
a1_12 = ARCL1(0,1);
a1_21 = ARCL1(1,0);
a1_22 = ARCL1(1,1);

a2_11 = ARCL2(0,0);
a2_12 = ARCL2(0,1);
a2_21 = ARCL2(1,0);
a2_22 = ARCL2(1,1);

v1_11 = log(discreteDiff1(0,0)); // we take the log because of the internal 
// transformation in lessSEM
v1_12 = discreteDiff1(0,1);
v1_22 = log(discreteDiff1(1,1)); // we take the log because of the internal 
// transformation in lessSEM

v2_11 = log(discreteDiff2(0,0)); // we take the log because of the internal 
// transformation in lessSEM
v2_12 = discreteDiff2(0,1);
v2_22 = log(discreteDiff2(1,1)); // we take the log because of the internal 
// transformation in lessSEM
"

Most importantly, we used matrices in our transformations. This is possible, because lessSEM uses RcppArmadillo in the background and therefore also provides users with all of the functions implemented therein.

As before, we can now fit the model with lessSEM

lessSEMFit <- bfgs(lavaanModel = lavaanFit,
                   # Our model modification must make use of the modifyModel - function:
                   modifyModel = modifyModel(transformations = transformations)
)

For comparison, we will also fit the model with ctsemOMX:

AnomAuthmodel <- ctModel(LAMBDA = matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2), 
                         Tpoints = 5, n.latent = 2, n.manifest = 2, MANIFESTVAR=diag(0, 2), TRAITVAR = NULL) 
AnomAuthfit <- ctFit(AnomAuth, AnomAuthmodel)

The following matrices are the drifts of the lessSEM model and the ctsemOMX model:

lessSEM:

matrix(coef(lessSEMFit)@estimates[,c("ctA_11", "ctA_21", "ctA_12", "ctA_22")],2,2)
#>            [,1]       [,2]
#> [1,] -0.4481642  0.2320794
#> [2,]  0.0425989 -0.1171247

ctsemOMX:

AnomAuthfit$mxobj$DRIFT$values
#>             [,1]       [,2]
#> [1,] -0.44728184  0.2324980
#> [2,]  0.04329283 -0.1174662

The diffusions are given by:

lessSEM:

matrix(coef(lessSEMFit)@estimates[,c("ctV_11", "ctV_12", "ctV_12", "ctV_22")],2,2)
#>              [,1]         [,2]
#> [1,]  0.473495271 -0.003846946
#> [2,] -0.003846946  0.154505468

ctsemOMX:

AnomAuthfit$mxobj$DIFFUSIONchol$result%*%t(AnomAuthfit$mxobj$DIFFUSIONchol$result)
#>              [,1]         [,2]
#> [1,]  0.473241884 -0.004610149
#> [2,] -0.004610149  0.154509547

Regularization could be used to enforce sparsity (Orzek & Voelkle, under review):

lassoFit <- lasso(lavaanModel = lavaanFit, 
                  regularized = "ctA_21",
                  nLambdas = 30,
                  method = "glmnet",
                  control = controlGlmnet(),
                  # Our model modification must make use of the modifyModel - function:
                  modifyModel = modifyModel(transformations = transformations)
)
plot(lassoFit)
plot of chunk unnamed-chunk-55
plot of chunk unnamed-chunk-55

Again, all these steps are easier when using dedicated packages such as ctsem and ctsemOMX for continuous time SEM and regCtsem (Orzek & Voelkle, under review) for regularized continuous time SEM.

How it is implemented

The basic idea behind the transformations is as follows: Assume that the SEM parameters are given by \(\pmb\theta\). The 2-log-likelihood of the model is given by \(f(\pmb\theta)\). When using transformations, we redefine \(\pmb\theta\) to be a function of other parameters, say \(\pmb\gamma\). That is, \(\pmb\theta =\pmb g(\pmb\gamma)\), where \(\pmb g\) is a function returning a vector. As a result, we can re-write the fitting function as \(f(\pmb\theta) = f(\pmb g(\pmb\gamma))\). Instead of optimizing \(\pmb\theta\), we now optimize \(\pmb\gamma\). Within lessSEM, the gradients of \(f(\pmb\theta)\) with respect to \(\pmb\theta\) are implemented in closed form as this results in a considerably faster run time. To get the gradients of \(f(\pmb\theta) = f(\pmb g(\pmb\gamma))\) with respect to \(\pmb\gamma\), lessSEM makes use of the chain rule. The gradients of the transformation \(\pmb g(\pmb\gamma)\) are approximated numerically. This is similar to the procedure by Arnold et al. (in submission) for continuous time SEM, where the derivative of the matrix exponential is approximated numerically, while all other elements are derived with closed form solutions.

Bibliography