Model Without Interaction
It can be instructive to start by considering a model in which we fix
\(\\lambda_{4} = 0\). This type of
model would be estimated with the following code:
lambda <- matrix(c(1, NA, NA), ncol = 1)
mod0 <- galamm(
formula = y ~ type + x:response + (0 + loading | id),
data = latent_covariates,
load.var = "type",
lambda = lambda,
factor = "loading"
)
In the data generating simulations, the true values were \(\lambda_{1}=1\), \(\lambda_{2} = 1.3\) and \(\lambda_{3} = -0.3\). The former two are
very well recovered, but the latter is too positive, which is likely due
to us omitting the interaction \(\lambda_{4}\), whose true value was
0.2.
summary(mod0)
#> GALAMM fit by maximum marginal likelihood.
#> Formula: y ~ type + x:response + (0 + loading | id)
#> Data: latent_covariates
#>
#> AIC BIC logLik deviance df.resid
#> 197.5 232.7 -90.8 181.5 592
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.92130 -0.52315 -0.01058 0.52532 3.01104
#>
#> Lambda:
#> loading SE
#> lambda1 1.0000 .
#> lambda2 1.3033 0.013687
#> lambda3 -0.2063 0.008515
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> id loading 0.98112 0.9905
#> Residual 0.01357 0.1165
#> Number of obs: 600, groups: id, 200
#>
#> Fixed effects:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.010588 0.07052 -0.15014 8.807e-01
#> typemeasurement2 -0.002173 0.02423 -0.08968 9.285e-01
#> typeresponse 0.029353 0.08678 0.33826 7.352e-01
#> x:response 0.470444 0.03093 15.20943 3.062e-52
Linear Interaction Between Observed and Latent Covariates
The measurement model can be equivalently written as
\[
\begin{pmatrix}
y_{1} \\
y_{2} \\
y_{3}
\end{pmatrix}
=
\boldsymbol{\beta}_{0} +
\begin{pmatrix}
1 \\
\lambda_{2} \\
\lambda_{3} + \lambda_{4} x
\end{pmatrix}
\eta
+
\begin{pmatrix}
0 \\
0 \\
x \beta
\end{pmatrix}
+ \boldsymbol{\epsilon}.
\]
This way of writing shows more explicitly which factor loadings are
connected with which observation. In order to fit this model with
galamm, we must provide formulas for the terms in the
loading matrix
\[
\begin{pmatrix}
1 \\
\lambda_{2} \\
\lambda_{3} + \lambda_{4} x
\end{pmatrix}.
\]
We specify the factor interactions with a list, one for each row of
lambda:
factor_interactions <- list(~ 1, ~ 1, ~ x)
This specifies that for the first two rows, there are no covariates,
but for the third row, we want a linear regression with \(x\) as covariate. Next, we specify the
loading matrix without the interaction parameter, i.e.,
we reuse the lambda object that was specified for
mod0 above. This lets us fit the model as follows:
mod <- galamm(
formula = y ~ type + x:response + (0 + loading | id),
data = latent_covariates,
load.var = "type",
lambda = lambda,
factor = "loading",
factor_interactions = factor_interactions
)
A model comparison shows overwhelming evidence in favor of this
model, which is not surprising since this is how the data were
simulated.
anova(mod, mod0)
#> Data: latent_covariates
#> Models:
#> mod0: y ~ type + x:response + (0 + loading | id)
#> mod: y ~ type + x:response + (0 + loading | id)
#> npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
#> mod0 8 197.52 232.70 -90.761 120.31
#> mod 9 138.31 177.88 -60.155 120.31 61.21 1 5.129e-15 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The summary also shows that the bias in \(\lambda_{3}\) has basically disappeared, as
it is up to -0.318 from -0.195, with the true value being -0.3. The
interaction is estimated at 0.233, which is also very close to the true
value 0.2. It should of course be noted here that the noise level in
this simulated dataset was set unrealistically low, to let us confirm
that the implementation itself is correct.
summary(mod)
#> GALAMM fit by maximum marginal likelihood.
#> Formula: y ~ type + x:response + (0 + loading | id)
#> Data: latent_covariates
#>
#> AIC BIC logLik deviance df.resid
#> 138.3 177.9 -60.2 120.3 591
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.2033 -0.5251 -0.0273 0.5146 3.5029
#>
#> Lambda:
#> loading SE
#> lambda1 1.0000 .
#> lambda2 1.3034 0.01268
#> lambda3 -0.3183 0.01604
#> lambda4_x 0.2331 0.02873
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> id loading 0.98175 0.9908
#> Residual 0.01164 0.1079
#> Number of obs: 600, groups: id, 200
#>
#> Fixed effects:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.010589 0.07048 -0.15024 8.806e-01
#> typemeasurement2 -0.002173 0.02384 -0.09116 9.274e-01
#> typeresponse 0.034005 0.09417 0.36109 7.180e-01
#> x:response 0.462507 0.03300 14.01556 1.252e-44
Interaction Between Latent Covariate and a Quadratic Term
We can also try to add interactions between the \(x^{2}\) and \(\eta\). We first update the formula in
factor_interactions:
factor_interactions <- list(~ 1, ~ 1, ~ x + I(x^2))
Then we fit the model as before:
mod2 <- galamm(
formula = y ~ type + x:response + (0 + loading | id),
data = latent_covariates,
load.var = "type",
lambda = lambda,
factor = "loading",
factor_interactions = factor_interactions
)
As can be seen, the coefficient for this squared interaction is not
significantly different from zero.
summary(mod2)
#> GALAMM fit by maximum marginal likelihood.
#> Formula: y ~ type + x:response + (0 + loading | id)
#> Data: latent_covariates
#>
#> AIC BIC logLik deviance df.resid
#> 140.3 184.2 -60.1 120.3 590
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.2012 -0.5212 -0.0228 0.5173 3.4968
#>
#> Lambda:
#> loading SE
#> lambda1 1.00000 .
#> lambda2 1.30340 0.01267
#> lambda3 -0.31418 0.02446
#> lambda4_x 0.20906 0.11252
#> lambda5_I(x^2) 0.02453 0.11097
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> id loading 0.98175 0.9908
#> Residual 0.01164 0.1079
#> Number of obs: 600, groups: id, 200
#>
#> Fixed effects:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.009919 0.07054 -0.14062 8.882e-01
#> typemeasurement2 -0.001969 0.02385 -0.08256 9.342e-01
#> typeresponse 0.033289 0.09426 0.35316 7.240e-01
#> x:response 0.462228 0.03305 13.98397 1.953e-44
