Interaction effects between two endogenous (i.e., dependent)
variables work as you would expect for the product indicator methods
(`"dblcent", "rca", "ca", "uca"`

). For the lms- and qml
approach however, it is not as straight forward.

The lms- and qml approach can (by default) handle interaction effects between endogenous and exogenous (i.e., independent) variables, but not interaction effects between two endogenous variables. When there is an interaction effect between two endogenous variables, the equations cannot easily be written in ‘reduced’ form – meaning that normal estimation procedures won’t work.

This being said, there is a work-around for these limitations for both the lms- and qml-approach. In essence, the model can be split into two parts, one linear and one non-linear. Basically, you can replace the covariance matrix used in the estimation of the non-linear model, with the model-implied covariance matrix from a linear model. Thus you can treat an endogenous variable as if it were exogenous – given that it can be expressed in a linear model.

Let’s consider the the theory of planned behaviour (TPB) where we wish to estimate the quadratic effect of INT on BEH (INT:INT). With the following model:

```
tpb <- '
# Outer Model (Based on Hagger et al., 2007)
ATT =~ att1 + att2 + att3 + att4 + att5
SN =~ sn1 + sn2
PBC =~ pbc1 + pbc2 + pbc3
INT =~ int1 + int2 + int3
BEH =~ b1 + b2
# Inner Model (Based on Steinmetz et al., 2011)
INT ~ ATT + SN + PBC
BEH ~ INT + PBC
BEH ~ INT:INT
'
```

Since INT is an endogenous variable, its quadratic term (i.e., an interaction effect with itself) would include two endogenous variables. Thus we would ordinarily not be able to estimate this model using the lms- or qml-approach. However, we can split the model into two parts, one linear and one non-linear. While INT is an endogenous variable, it can be expressed in a linear model – since it is not affected by any interaction terms:

We could then remove this part from the original model, giving us:

```
tpb_nonlinear <- '
# Outer Model (Based on Hagger et al., 2007)
ATT =~ att1 + att2 + att3 + att4 + att5
SN =~ sn1 + sn2
PBC =~ pbc1 + pbc2 + pbc3
INT =~ int1 + int2 + int3
BEH =~ b1 + b2
# Inner Model (Based on Steinmetz et al., 2011)
BEH ~ INT + PBC
BEH ~ INT:INT
'
```

We could now just estimate the non-linear model, since INT now is an
exogenous variable. This would however not incorporate the structural
model for INT. To address this, we can make modsem replace the
covariance matrix (phi) of (INT, PBC, ATT, SN) with the model-implied
covariance matrix from the linear model, whilst estimating both models
simultaneously. To acheive this, we can use the `cov.syntax`

argument in `modsem`

: