# Non-linear latent variable models and error-in-variable models

#### Klaus Kähler Holst

#### 2025-10-29

    library('lava')

We consider the measurement models given by

\[X_{j} = \eta_{1} + \epsilon_{j}^{x}, \quad j=1,2,3\] \[Y_{j} = \eta_{2} + \epsilon_{j}^{y}, \quad j=1,2,3\] and with a structural model given by \[\eta_{2} = f(\eta_{1}) + Z + \zeta_{2}\label{ex:eta2}\] \[\eta_{1} = Z + \zeta_{1}\label{ex:eta1}\] with iid measurement errors \(\epsilon_{j}^{x},\epsilon_{j}^{y},\zeta_{1},\zeta_{2}\sim\mathcal{N}(0,1), j=1,2,3.\) and standard normal distributed covariate \(Z\). To simulate from this model we use the following syntax:

    f <- function(x) cos(1.25*x) + x - 0.25*x^2
    m <- lvm(x1+x2+x3 ~ eta1, y1+y2+y3 ~ eta2, latent=~eta1+eta2)
    regression(m) <- eta1+eta2 ~ z
    functional(m, eta2~eta1) <- f

    d <- sim(m, n=200, seed=42) # Default is all parameters are 1

    ## plot(m, plot.engine="visNetwork")
    plot(m)

We refer to (Holst and Budtz-Jørgensen 2013) for details on the syntax for model specification.

# Estimation

To estimate the parameters using the two-stage estimator described in (Holst and Budtz-Jørgensen 2020), the first step is now to specify the measurement models

    m1 <- lvm(x1+x2+x3 ~ eta1, eta1 ~ z, latent=~eta1)
    m2 <- lvm(y1+y2+y3 ~ eta2, eta2 ~ z, latent=~eta2)

Next, we specify a quadratic relationship between the two latent variables

    nonlinear(m2, type="quadratic") <- eta2 ~ eta1

and the model can then be estimated using the two-stage estimator

    e1 <- twostage(m1, m2, data=d)
    e1
    #>                     Estimate Std. Error Z-value  P-value
    #> Measurements:
    #>    y2~eta2             0.977      0.035  28.309   <1e-12
    #>    y3~eta2             1.045      0.035  29.982   <1e-12
    #> Regressions:
    #>    eta2~z              0.885      0.208   4.260    2e-05
    #>    eta2~eta1_1         1.141      0.174   6.552  5.7e-11
    #>    eta2~eta1_2        -0.451      0.072  -6.292  3.1e-10
    #> Intercepts:
    #>    y2                 -0.122      0.109  -1.117     0.26
    #>    y3                 -0.099      0.105  -0.937     0.35
    #>    eta2                0.678      0.174   3.906  9.4e-05
    #> Residual Variances:
    #>    y1                  1.307      0.177   7.368
    #>    y2                  1.111      0.145   7.671
    #>    y3                  0.810      0.132   6.132
    #>    eta2                2.085      0.290   7.193

We see a clear statistically significant effect of the second order term (`eta2~eta1_2`). For comparison we can also estimate the full MLE of the linear model:

    e0 <- estimate(regression(m1%++%m2, eta2~eta1), d)
    estimate(e0,keep="^eta2~[a-z]",regex=TRUE) ## Extract coef. matching reg.ex.
    #>           Estimate Std.Err    2.5% 97.5%   P-value
    #> eta2~eta1   1.4140  0.2261 0.97083 1.857 4.014e-10
    #> eta2~z      0.6374  0.2778 0.09291 1.182 2.177e-02

Next, we calculate predictions from the quadratic model using the estimated parameter coefficients \[ \mathbb{E}_{\widehat{\theta}_{2}}(\eta_{2} \mid \eta_{1}, Z=0), \]

    newd <- expand.grid(eta1=seq(-4, 4, by=0.1), z=0)
    pred1 <- predict(e1, newdata=newd, x=TRUE)
    head(pred1)
    #>           y1      y2      y3    eta2
    #> [1,] -11.094 -10.959 -11.690 -11.094
    #> [2,] -10.624 -10.500 -11.199 -10.624
    #> [3,] -10.163 -10.049 -10.717 -10.163
    #> [4,]  -9.711  -9.608 -10.245  -9.711
    #> [5,]  -9.268  -9.175  -9.782  -9.268
    #> [6,]  -8.834  -8.751  -9.329  -8.834

To obtain a potential better fit we next proceed with a natural cubic spline

    kn <- seq(-3,3,length.out=5)
    nonlinear(m2, type="spline", knots=kn) <- eta2 ~ eta1
    e2 <- twostage(m1, m2, data=d)
    e2
    #>                     Estimate Std. Error Z-value  P-value
    #> Measurements:
    #>    y2~eta2             0.978      0.035  28.313   <1e-12
    #>    y3~eta2             1.045      0.035  29.971   <1e-12
    #> Regressions:
    #>    eta2~z              0.867      0.203   4.278  1.9e-05
    #>    eta2~eta1_1         2.862      0.673   4.255  2.1e-05
    #>    eta2~eta1_2         0.003      0.101   0.034     0.97
    #>    eta2~eta1_3        -0.263      0.294  -0.894     0.37
    #>    eta2~eta1_4         0.508      0.352   1.443     0.15
    #> Intercepts:
    #>    y2                 -0.122      0.109  -1.116     0.26
    #>    y3                 -0.099      0.105  -0.936     0.35
    #>    eta2                1.838      1.664   1.104     0.27
    #> Residual Variances:
    #>    y1                  1.313      0.178   7.396
    #>    y2                  1.104      0.145   7.639
    #>    y3                  0.811      0.132   6.153
    #>    eta2                1.994      0.269   7.402

Confidence limits can be obtained via the Delta method using the `estimate` method:

    p <- cbind(eta1=newd$eta1,
          estimate(e2,f=function(p) predict(e2,p=p,newdata=newd))$coefmat)
    head(p)
    #>    eta1 Estimate Std.Err   2.5%  97.5%   P-value
    #> p1 -4.0   -9.611  1.2647 -12.09 -7.132 2.978e-14
    #> p2 -3.9   -9.325  1.2051 -11.69 -6.963 1.012e-14
    #> p3 -3.8   -9.039  1.1464 -11.29 -6.792 3.152e-15
    #> p4 -3.7   -8.752  1.0886 -10.89 -6.619 8.959e-16
    #> p5 -3.6   -8.466  1.0319 -10.49 -6.444 2.320e-16
    #> p6 -3.5   -8.180  0.9766 -10.09 -6.266 5.494e-17

The fitted function can be obtained with the following code:

    plot(I(eta2-z) ~ eta1, data=d, col=Col("black",0.5), pch=16,
         xlab=expression(eta[1]), ylab=expression(eta[2]), xlim=c(-4,4))
    lines(Estimate ~ eta1, data=as.data.frame(p), col="darkblue", lwd=5)
    confband(p[,1], lower=p[,4], upper=p[,5], polygon=TRUE,
         border=NA, col=Col("darkblue",0.2))

# Cross-validation

A more formal comparison of the different models can be obtained by cross-validation. Here we specify linear, quadratic and cubic spline models with 4 and 9 degrees of freedom.

    m2a <- nonlinear(m2, type="linear", eta2~eta1)
    m2b <- nonlinear(m2, type="quadratic", eta2~eta1)
    kn1 <- seq(-3,3,length.out=5)
    kn2 <- seq(-3,3,length.out=8)
    m2c <- nonlinear(m2, type="spline", knots=kn1, eta2~eta1)
    m2d <- nonlinear(m2, type="spline", knots=kn2, eta2~eta1)

To assess the model fit average RMSE is estimated with 5-fold cross-validation repeated two times

    ## Scale models in stage 2 to allow for a fair RMSE comparison
    d0 <- d
    for (i in endogenous(m2))
        d0[,i] <- scale(d0[,i],center=TRUE,scale=TRUE)
    ## Repeated 5-fold cross-validation:
    ff <- lapply(list(linear=m2a,quadratic=m2b,spline4=m2c,spline6=m2d),
            function(m) function(data,...) twostage(m1,m,data=data,stderr=FALSE,control=list(start=coef(e0),contrain=TRUE)))
    fit.cv <- lava:::cv(ff,data=d,K=5,rep=2,mc.cores=parallel::detectCores(),seed=1)

    summary(fit.cv)
    #>       Length Class  Mode
    #> cv    40     -none- numeric
    #> call   7     -none- call
    #> names  4     -none- character
    #> rep    1     -none- numeric
    #> folds  1     -none- numeric

Here the RMSE is in favour of the splines model with 4 degrees of freedom:

    fit <- lapply(list(m2a,m2b,m2c,m2d),
             function(x) {
             e <- twostage(m1,x,data=d)
             pr <- cbind(eta1=newd$eta1,predict(e,newdata=newd$eta1,x=TRUE))
             return(list(estimate=e,predict=as.data.frame(pr)))
             })

    plot(I(eta2-z) ~ eta1, data=d, col=Col("black",0.5), pch=16,
         xlab=expression(eta[1]), ylab=expression(eta[2]), xlim=c(-4,4))
    col <- c("orange","darkred","darkgreen","darkblue")
    lty <- c(3,4,1,5)
    for (i in seq_along(fit)) {
        with(fit[[i]]$pr, lines(eta2 ~ eta1, col=col[i], lwd=4, lty=lty[i]))
    }
    legend("bottomright",
          c("linear","quadratic","spline(df=4)","spline(df=6)"),
          col=col, lty=lty, lwd=3)

For convenience, the function `twostageCV` can be used to do the cross-validation (also for choosing the mixture distribution via the ``nmix`` argument, see the section below). For example,

    selmod <- twostageCV(m1, m2, data=d, df=2:4, nmix=1:2,
            nfolds=2, rep=1, mc.cores=parallel::detectCores())

applies cross-validation (here just 2 folds for simplicity) to select the best splines with degrees of freedom varying from from 1-3 (the linear model is automatically included)

    selmod
    #>        Length Class               Mode
    #> model1 11     summary.lvm.mixture list
    #> AIC1    2     -none-              numeric
    #> cv      4     -none-              numeric
    #> knots   4     -none-              list
    #> model2 11     summary.lvmfit      list

# Specification of general functional forms

Next, we show how to specify a general functional relation of multiple different latent or exogenous variables. This is achieved via the `predict.fun` argument. To illustrate this we include interactions between the latent variable \(\eta_{1}\) and a dichotomized version of the covariate \(z\)

    d$g <- (d$z<0)*1 ## Group variable
    mm1 <- regression(m1, ~g)  # Add grouping variable as exogenous variable (effect specified via 'predict.fun')
    mm2 <- regression(m2, eta2~ u1+u2+u1:g+u2:g+z)
    pred <- function(mu,var,data,...) {
        cbind("u1"=mu[,1],"u2"=mu[,1]^2+var[1],
          "u1:g"=mu[,1]*data[,"g"],"u2:g"=(mu[,1]^2+var[1])*data[,"g"])
    }
    ee1 <- twostage(mm1, model2=mm2, data=d, predict.fun=pred)
    estimate(ee1,keep="eta2~u",regex=TRUE)
    #>           Estimate Std.Err    2.5%   97.5%  P-value
    #> eta2~u1     0.9891  0.3020  0.3971  1.5810 0.001057
    #> eta2~u2    -0.3962  0.1443 -0.6791 -0.1133 0.006047
    #> eta2~u1:g   0.4487  0.4620 -0.4568  1.3543 0.331409
    #> eta2~u2:g   0.0441  0.2166 -0.3804  0.4686 0.838667

A formal test show no statistically significant effect of this interaction

    summary(estimate(ee1,keep="(:g)", regex=TRUE))
    #> Call: estimate.default(x = ee1, keep = "(:g)", regex = TRUE)
    #> ────────────────────────────────────────────────────────────────────────────────
    #>           Estimate Std.Err    2.5%  97.5% P-value
    #> eta2~u1:g   0.4487  0.4620 -0.4568 1.3543  0.3314
    #> eta2~u2:g   0.0441  0.2166 -0.3804 0.4686  0.8387
    #>
    #>  Null Hypothesis:
    #>   [eta2~u1:g] = 0
    #>   [eta2~u2:g] = 0
    #>
    #> chisq = 0.9441, df = 2, p-value = 0.6237

# Mixture models

Lastly, we demonstrate how the distributional assumptions of stage 1 model can be relaxed by letting the conditional distribution of the latent variable given covariates follow a Gaussian mixture distribution. The following code explictly defines the parameter constraints of the model by setting the intercept of the first indicator variable, \(x_{1}\), to zero and the factor loading parameter of the same variable to one.

    m1 <- baptize(m1)  ## Label all parameters
    intercept(m1, ~x1+eta1) <- list(0,NA) ## Set intercept of x1 to zero. Remove the label of eta1
    regression(m1,x1~eta1) <- 1 ## Factor loading fixed to 1

The mixture model may then be estimated using the `mixture` method (note, this requires the `mets` package to be installed), where the Parameter names shared across the different mixture components given in the `list` will be constrained to be identical in the mixture model. Thus, only the intercept of \(\eta_{1}\) is allowed to vary between the mixtures.

    set.seed(1)
    em0 <- mixture(m1, k=2, data=d)

To decrease the risk of using a local maximizer of the likelihood we can rerun the estimation with different random starting values

    em0 <- NULL
    ll <- c()
    for (i in 1:5) {
        set.seed(i)
        em <- mixture(m1, k=2, data=d, control=list(trace=0))
        ll <- c(ll,logLik(em))
        if (is.null(em0) || logLik(em0)<tail(ll,1))
        em0 <- em
    }

    summary(em0)
    #> Cluster 1 (n=162, Prior=0.776):
    #> --------------------------------------------------
    #>                     Estimate Std. Error Z value Pr(>|z|)
    #> Measurements:
    #>    x1~eta1           1.000
    #>    x2~eta1           0.996    0.079     12.541    <1e-12
    #>    x3~eta1           1.063    0.084     12.605    <1e-12
    #> Regressions:
    #>    eta1~z            1.067    0.085     12.510    <1e-12
    #> Intercepts:
    #>    x1                0.000
    #>    x2                0.038    0.099      0.389  0.7
    #>    x3               -0.025    0.103     -0.247  0.81
    #>    eta1              0.209    0.132      1.590  0.11
    #> Residual Variances:
    #>    x1                0.985    0.133      7.400
    #>    x2                0.972    0.132      7.387
    #>    x3                1.013    0.143      7.088
    #>    eta1              0.290    0.111      2.610
    #>
    #> Cluster 2 (n=38, Prior=0.224):
    #> --------------------------------------------------
    #>                     Estimate Std. Error Z value Pr(>|z|)
    #> Measurements:
    #>    x1~eta1           1.000
    #>    x2~eta1           0.996    0.079     12.541    <1e-12
    #>    x3~eta1           1.063    0.084     12.605    <1e-12
    #> Regressions:
    #>    eta1~z            1.067    0.085     12.510    <1e-12
    #> Intercepts:
    #>    x1                0.000
    #>    x2                0.038    0.099      0.389  0.7
    #>    x3               -0.025    0.103     -0.247  0.81
    #>    eta1             -1.443    0.259     -5.578  2.4e-08
    #> Residual Variances:
    #>    x1                0.985    0.133      7.400
    #>    x2                0.972    0.132      7.387
    #>    x3                1.013    0.143      7.088
    #>    eta1              0.290    0.111      2.610
    #> --------------------------------------------------
    #> AIC= 1959
    #> ||score||^2= 8.818e-06

Measured by AIC there is a slight improvement in the model fit using the mixture model

    e0 <- estimate(m1,data=d)
    AIC(e0,em0)
    #>     df  AIC
    #> e0  10 1962
    #> em0 12 1959

The spline model may then be estimated as before with the `two-stage` method

    em2 <- twostage(em0,m2,data=d)
    em2
    #>                     Estimate Std. Error Z-value  P-value
    #> Measurements:
    #>    y2~eta2             0.978      0.035  28.237   <1e-12
    #>    y3~eta2             1.045      0.035  30.040   <1e-12
    #> Regressions:
    #>    eta2~z              1.029      0.223   4.607  4.1e-06
    #>    eta2~eta1_1         2.804      0.655   4.278  1.9e-05
    #>    eta2~eta1_2        -0.022      0.100  -0.225     0.82
    #>    eta2~eta1_3        -0.173      0.289  -0.599     0.55
    #>    eta2~eta1_4         0.387      0.340   1.138     0.26
    #> Intercepts:
    #>    y2                 -0.122      0.109  -1.114     0.27
    #>    y3                 -0.099      0.105  -0.936     0.35
    #>    eta2                2.124      1.666   1.275      0.2
    #> Residual Variances:
    #>    y1                  1.319      0.177   7.470
    #>    y2                  1.097      0.145   7.564
    #>    y3                  0.813      0.133   6.135
    #>    eta2                1.996      0.283   7.055

In this example the results are very similar to the Gaussian model:

    plot(I(eta2-z) ~ eta1, data=d, col=Col("black",0.5), pch=16,
         xlab=expression(eta[1]), ylab=expression(eta[2]))

    lines(Estimate ~ eta1, data=as.data.frame(p), col="darkblue", lwd=5)
    confband(p[,1], lower=p[,4], upper=p[,5], polygon=TRUE,
         border=NA, col=Col("darkblue",0.2))

    pm <- cbind(eta1=newd$eta1,
            estimate(em2, f=function(p) predict(e2,p=p,newdata=newd))$coefmat)
    lines(Estimate ~ eta1, data=as.data.frame(pm), col="darkred", lwd=5)
    confband(pm[,1], lower=pm[,4], upper=pm[,5], polygon=TRUE,
         border=NA, col=Col("darkred",0.2))
    legend("bottomright", c("Gaussian","Mixture"),
           col=c("darkblue","darkred"), lwd=2, bty="n")

# SessionInfo

    sessionInfo()
    #> R version 4.5.1 (2025-06-13)
    #> Platform: aarch64-apple-darwin25.0.0
    #> Running under: macOS Tahoe 26.0.1
    #>
    #> Matrix products: default
    #> BLAS:   /Users/kkzh/.asdf/installs/r/4.5.1/lib/R/lib/libRblas.dylib
    #> LAPACK: /Users/kkzh/.asdf/installs/r/4.5.1/lib/R/lib/libRlapack.dylib;  LAPACK version 3.12.1
    #>
    #> locale:
    #> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
    #>
    #> time zone: Europe/Copenhagen
    #> tzcode source: internal
    #>
    #> attached base packages:
    #> [1] stats     graphics  grDevices utils     datasets  methods   base
    #>
    #> other attached packages:
    #> [1] survival_3.8-3 lava_1.8.2
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] tidyr_1.3.1         sass_0.4.10         future_1.67.0
    #>  [4] generics_0.1.4      lattice_0.22-7      listenv_0.9.1
    #>  [7] digest_0.6.37       magrittr_2.0.4      evaluate_1.0.5
    #> [10] grid_4.5.1          mvtnorm_1.3-3       fastmap_1.2.0
    #> [13] jsonlite_2.0.0      Matrix_1.7-4        backports_1.5.0
    #> [16] graph_1.86.0        purrr_1.1.0         Rgraphviz_2.52.0
    #> [19] codetools_0.2-20    numDeriv_2016.8-1.1 jquerylib_0.1.4
    #> [22] cli_3.6.5           rlang_1.1.6         mets_1.3.8
    #> [25] parallelly_1.45.1   future.apply_1.20.0 splines_4.5.1
    #> [28] geepack_1.3.13      cachem_1.1.0        yaml_2.3.10
    #> [31] tools_4.5.1         parallel_4.5.1      dplyr_1.1.4
    #> [34] globals_0.18.0      BiocGenerics_0.54.1 broom_1.0.10
    #> [37] vctrs_0.6.5         R6_2.6.1            stats4_4.5.1
    #> [40] lifecycle_1.0.4     MASS_7.3-65         pkgconfig_2.0.3
    #> [43] timereg_2.0.7       progressr_0.17.0    bslib_0.9.0
    #> [46] pillar_1.11.1       glue_1.8.0          Rcpp_1.1.0
    #> [49] xfun_0.53           tibble_3.3.0        tidyselect_1.2.1
    #> [52] knitr_1.50          htmltools_0.5.8.1   rmarkdown_2.30
    #> [55] compiler_4.5.1

# Bibliography

Holst, K. K., and E. Budtz-Jørgensen. 2013. “Linear Latent Variable Models: The Lava-Package.” Computational Statistics 28 (4): 1385–452. https://doi.org/10.1007/s00180-012-0344-y.
Holst, Klaus Kähler, and Esben Budtz-Jørgensen. 2020. “A Two-Stage Estimation Procedure for Non-Linear Structural Equation Models.” Biostatistics (in press). https://doi.org/10.1093/biostatistics/kxy082.