Parameter-transformations.Rmd

---
title: "Parameter-transformations"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Parameter-transformations}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---



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`):


```r
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:


```r
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:


```r
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:


```r
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:




```r
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:


```r
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:


```r
plot(lassoFit)
```

![plot of chunk unnamed-chunk-9](lessSEMFigures/lessSEM-unnamed-chunk-9-1.png)

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


```r
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:


```r
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:


```r
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:


```r
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](https://github.com/jhorzek/lessSEM/blob/main/vignettes/Parameter-transformations.Rmd)).
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:

```r
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:


```r
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:




```r
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:


```r
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](https://github.com/RcppCore/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:


```r
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:


```r
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:


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

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

```r
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:

```r
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:

```r
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](https://github.com/cdriveraus/ctsem) and [dynr](https://github.com/mhunter1/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](https://github.com/jhorzek/lessSEM/blob/main/vignettes/Parameter-transformations.Rmd)).



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:


```r
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
```


```r
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:


```r
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.
"
```




```r
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:

```r
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**:



```r
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:


```r
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:

```r
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:


```r
# 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:

```r
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.


```r
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**




```r
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**:



```r
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:

```r
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:

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

The diffusions are given by:

lessSEM:

```r
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:

```r
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):




```r
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)
)
```


```r
plot(lassoFit)
```

![plot of chunk unnamed-chunk-55](lessSEMFigures/lessSEM-unnamed-chunk-55-1.png)

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](https://en.wikipedia.org/wiki/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

* Arnold, M., Cancér, P. F., Estrada, E., & Voelkle, M. C. (in submission). Score-Guided Recursive Partitioning of Continuous-Time Structural Equation Models.
* Belzak, W. C. M., & Bauer, D. J. (2020). Improving the assessment of measurement invariance: Using regularization to select anchor items and identify differential item functioning. Psychological Methods, 25(6), 673–690. https://doi.org/10.1037/met0000253
* Driver, C. C., Oud, J. H. L., & Voelkle, M. C. (2017). Continuous time structural equation modelling with R package ctsem. Journal of Statistical Software, 77(5), 1–36. https://doi.org/10.18637/jss.v077.i05
* Fisher, Z. F., Kim, Y., Fredrickson, B. L., & Pipiras, V. (2022). Penalized Estimation and Forecasting of Multiple Subject Intensive Longitudinal Data. Psychometrika, 87(2), 1–29. https://doi.org/10.1007/s11336-021-09825-7
* Huang, P.-H. (2018). A penalized likelihood method for multi-group structural equation modelling. British Journal of Mathematical and Statistical Psychology, 71(3), 499–522. https://doi.org/10.1111/bmsp.12130
* Hunter, M. D. (2018). State space modeling in an open source, modular, structural equation modeling
  environment. Structural Equation Modeling, 25(2), 307-324. doi: 10.1080/10705511.2017.1369354
* Jacobucci, R., & Grimm, K. J. (2018). Regularized Estimation of Multivariate Latent Change Score Models. In E. Ferrer, S. M. Boker, & K. J. Grimm (Eds.), Longitudinal Multivariate Psychology (1st ed., pp. 109–125). Routledge. https://doi.org/10.4324/9781315160542-6
* Jacobucci, R., Grimm, K. J., Brandmaier, A. M., Serang, S., Kievit, R. A., & Scharf, F. (2019). regsem: Regularized structural equation modeling. https://CRAN.R-project.org/package=regsem
* Liang, X., Yang, Y., & Huang, J. (2018). Evaluation of Structural Relationships in Autoregressive Cross-Lagged Models Under Longitudinal Approximate Invariance:A Bayesian Analysis. Structural Equation Modeling: A Multidisciplinary Journal, 25(4), 558–572. https://doi.org/10.1080/10705511.2017.1410706
* Muthen, B., & Asparouhov, T. (2013). BSEM Measurement Invariance Analysis. Mplus Web Notes: No. 17.
* Neale, M. C., Hunter, M. D., Pritikin, J. N., Zahery, M., Brick, T. R., Kirkpatrick, R. M., Estabrook, R., Bates, T. C., Maes, H. H., & Boker, S. M. (2016). OpenMx 2.0: Extended structural equation and statistical modeling. Psychometrika, 81(2), 535–549. https://doi.org/10.1007/s11336-014-9435-8
* Ou, L., Hunter, M., D., & Chow, S.-M. (2019). Whats for dynr: A package for linear and nonlinear dynamic modeling in r. The R Journal, 11(1), 91–111. https://doi.org/10.32614/RJ-2019-012
* Orzek, J. H., & Voelkle, M. C. (under review). Regularized continuous time structural equation models: A network perspective.
* Pritikin, J. N., Hunter, M. D., & Boker, S. M. (2015). Modular open-source software for Item Factor
  Analysis. Educational and Psychological Measurement, 75(3), 458-474
* Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., & Knight, K. (2005). Sparsity and smoothness via the fused lasso. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(1), 91–108. https://doi.org/10.1111/j.1467-9868.2005.00490.x
* Voelkle, M. C., Oud, J. H. L., Davidov, E., & Schmidt, P. (2012). An sem approach to continuous time modeling of panel data: Relating authoritarianism and anomia. Psychological Methods, 17(2), 176–192. https://doi.org/10.1037/a0027543