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

Klaus Kähler Holst

2024-03-05

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 (K. K. 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 (Klaus Kähler 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.3.2 (2023-10-31)
#> Platform: aarch64-apple-darwin22.6.0 (64-bit)
#> Running under: macOS Sonoma 14.3.1
#> 
#> Matrix products: default
#> BLAS:   /Users/kkzh/.asdf/installs/R/4.3.2/lib/R/lib/libRblas.dylib 
#> LAPACK: /Users/kkzh/.asdf/installs/R/4.3.2/lib/R/lib/libRlapack.dylib;  LAPACK version 3.11.0
#> 
#> 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.5-7 lava_1.8.0    
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyr_1.3.0          sass_0.4.7           utf8_1.2.4          
#>  [4] future_1.33.1        generics_0.1.3       futile.options_1.0.1
#>  [7] lattice_0.22-5       pracma_2.4.4         listenv_0.9.1       
#> [10] digest_0.6.34        magrittr_2.0.3       evaluate_0.23       
#> [13] grid_4.3.2           mvtnorm_1.2-4        fastmap_1.1.1       
#> [16] jsonlite_1.8.8       Matrix_1.6-5         backports_1.4.1     
#> [19] graph_1.80.0         formatR_1.14         targeted_0.5        
#> [22] purrr_1.0.2          fansi_1.0.6          Rgraphviz_2.46.0    
#> [25] codetools_0.2-19     numDeriv_2016.8-1.1  jquerylib_0.1.4     
#> [28] cli_3.6.2            rlang_1.1.3          futile.logger_1.4.3 
#> [31] mets_1.3.4           parallelly_1.37.1    future.apply_1.11.1 
#> [34] splines_4.3.2        geepack_1.3.9        cachem_1.0.8        
#> [37] yaml_2.3.7           tools_4.3.2          parallel_4.3.2      
#> [40] nloptr_2.0.3         dplyr_1.1.3          optimx_2023-10.21   
#> [43] lambda.r_1.2.4       globals_0.16.2       BiocGenerics_0.48.0 
#> [46] broom_1.0.5          vctrs_0.6.5          R6_2.5.1            
#> [49] stats4_4.3.2         lifecycle_1.0.4      MASS_7.3-60         
#> [52] pkgconfig_2.0.3      timereg_2.0.5        progressr_0.14.0    
#> [55] bslib_0.5.1          pillar_1.9.0         data.table_1.15.2   
#> [58] glue_1.7.0           Rcpp_1.0.12          tidyselect_1.2.0    
#> [61] xfun_0.41            tibble_3.2.1         highr_0.10          
#> [64] knitr_1.45           htmltools_0.5.6.1    rmarkdown_2.25      
#> [67] compiler_4.3.2

Bibliography

Holst, K. K., and E. Budtz-Jørgensen. 2013. “Linear Latent Variable Models: The Lava-Package.” Computational Statistics 28 (4): 1385–1452. 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.