Nonlinear latent variable models

ML-inference in non-linear SEMs is complex. Computational intensive methods based on numerical integration are needed and results are sensitive to distributional assumptions.

In a recent paper: A two-stage estimation procedure for non-linear structural equation models by Klaus Kähler Holst & Esben Budtz-Jørgensen (https://doi.org/10.1093/biostatistics/kxy082), we consider two-stage estimators as a computationally simple alternative to MLE. Here both steps are based on linear models: first we predict the non-linear terms and then these are related to latent outcomes in the second step.

$$$$

Introduction: the measurement error problem

Measurement error in covariates is a common source of bias in regression analyses. In a simple linear regression setting, let $Y$ denote the response variable, $\xi$ the true exposure, for which we only observe replicated proxy measurements $X_1,\cdots ,X_q$, and $Z$ an additional confounder:

$$\begin{array}{rl}Y& =\mu +\beta \xi +\gamma Z+\delta \\ \xi & =\alpha +\kappa Z+\zeta \\ X_j& ={\nu }_j+\xi +{ϵ}_j,j=1,\dots ,q,\end{array}$$

where the measurement error terms $\delta ,\zeta ,{ϵ}_j$ are assumed to be independent of the true exposure $\xi$. Standard techniques, e.g. linear regression (MLE) on $X_j$ and $Z$ fails to provide consistent estimates. To see this observe that

$$\begin{array}{r}{\hat{\beta }}_{MLE}=\frac{\mathbb{C}\mathrm{ov}(Y,X_j\mid Z)}{\mathbb{V}\mathrm{ar}(X_j\mid Z)}=\frac{\beta \mathbb{V}\mathrm{ar}(\xi )}{\mathbb{V}\mathrm{ar}(\xi )+\mathbb{V}\mathrm{ar}({ϵ}_j)}<\beta \end{array}$$

In the following, generalizations of the above problem to general latent variable models including non-linear effects of the exposure variable are considered.

Statistical model

For subjects $i=1,\dots ,n$, we have latent response variable: ${\eta }_i=({\eta }_{i1},\dots ,{\eta }_{ip}{)}^t$ non-linearly related to latent predictor: ${\xi }_i=({\xi }_{i1},\dots ,{\xi }_{iq}{)}^t$ after adjustment for covariates $Z_i=(Z_{i1},\dots ,Z_{ir})$

$$\begin{array}{}\text{(1)}& \begin{array}{lcl}{\eta }_i& =& \alpha +B\phi ({\xi }_i)+\mathrm{\Gamma }{Z}_i+{\zeta }_i,\end{array}\end{array}$$

where $\phi :{\mathbb{R}}^p\to {\mathbb{R}}^l$ is a known measurable function such that $\phi ({\xi }_i)=({\phi }_1({\xi }_i),\dots ,{\phi }_l({\xi }_i){)}^t$ has finite variance (see Figure 1). The target parameter of interest parametrizing $B\in {\mathbb{R}}^{p\times l}$ describes the non-linear relation between ${\xi }_i$ and ${\eta }_i$. Note, that $\phi$ may also depend on some of the covariates thereby allowing the introduction of interaction terms.

The latent predictors (${\xi }_i$) are related to each other and the covariates in a linear structural equation:

$$\begin{array}{}\text{(2)}& \begin{array}{lcl}{\xi }_i& =& \tilde{\alpha }+\tilde{B}{\xi }_i+\tilde{\mathrm{\Gamma }}{Z}_i+{\tilde{\zeta }}_i.\end{array}\end{array}$$

where diagonal elements in $\tilde{B}$ are zero.

The observed variables $X_i=(X_{i1},\dots ,X_{ih}{)}^t$ and $Y_i=(Y_{i1},\dots ,Y_{im}{)}^t$ are linked to the latent variable in the two measurement models

$$\begin{array}{}\text{(3)}& Y_i=\nu +\mathrm{\Lambda }{\eta }_i+K{Z}_i+{\epsilon }_i\text{(4)}& X_i=\tilde{\nu }+\tilde{\mathrm{\Lambda }}{\xi }_i+\tilde{K}{Z}_i+{\tilde{\epsilon }}_i,\end{array}$$

where the error terms ${\epsilon }_i$ and ${\tilde{\epsilon }}_i$ are assumed to be independent with mean zero and covariance matrices of $\mathrm{\Omega }$ and $\tilde{\mathrm{\Omega }}$, respectively. The Øparameters* are collected into $\theta =({\theta }_1^t,{\theta }_2^t{)}^t$, where ${\theta }_1=(\tilde{\alpha },\tilde{B},\tilde{\gamma },\tilde{\nu },\tilde{\mathrm{\Lambda }},\tilde{K},\tilde{\mathrm{\Omega }},\tilde{\mathrm{\Psi }})$ are the parameters of the linear SEM describing the conditional distribution of $X_i$ given $Z_i$. The rest of the parameters are collected into ${\theta }_2$.

Two-stage estimator (2SSEM)

1. The linear SEM given by equations ($\text{2}$) and is ($\text{4}$) is fitted (e.g., MLE) to $(X_i,Z_i),i=1,\dots ,n$ to obtain consistent estimate of ${\theta }_1$.
2. The parameter ${\theta }_2$ is estimated via a linear SEM with measurement model given by equation ($\text{3}$) and structural model given by equation ($\text{1}$), where the latent predictor, $\phi ({\xi }_i)$, is replaced by the Empirical Bayes plug-in estimate ${\phi }_n^{\ast }({\xi }_i)={\mathbb{E}}_{{\hat{\theta }}_1}[\phi ({\xi }_i)\mid X_i,Z_i]$.

Prediction of non-linear latent terms

For important classes of non-linear functions ($\phi$) closed form expressions can be derived for ${\phi }^{\ast }({\xi }_i)=\mathbb{E}[\phi ({\xi }_i)|X_i,Z_i]$. Under Gaussian assumptions the conditional distribution of ${\xi }_i$ given $X_i$ and $Z_i$ is normal with mean and variance

$$\begin{array}{rl}{m}_{x,z}& =\mathbb{E}({\xi }_i|X_i,Z_i)=\tilde{\alpha }+\tilde{\mathrm{\Gamma }}{Z}_i+{\mathrm{\Sigma }}_{X\xi }{\mathrm{\Sigma }}_X^{-1}(X_i-{\mu }_X)\\ v& =\mathbb{V}\mathrm{ar}({\xi }_i|X_i,Z_i)=\tilde{\mathrm{\Psi }}-{\mathrm{\Sigma }}_{X\xi }{\mathrm{\Sigma }}_X^{-1}{\mathrm{\Sigma }}_{X\xi }^t,\end{array}$$

where ${\mu }_X=\tilde{\nu }+\tilde{\mathrm{\Lambda }}(I-\tilde{B}{)}^{-1}\alpha +\tilde{\mathrm{\Lambda }}(I-\tilde{B}{)}^{-1}\tilde{\mathrm{\Gamma }}{Z}_i+\tilde{K}{Z}_i,{\mathrm{\Sigma }}_X=\tilde{\mathrm{\Lambda }}(I-\tilde{B}{)}^{-1}\tilde{\mathrm{\Psi }}(I-\tilde{B}{)}^{-1t}{\tilde{\mathrm{\Lambda }}}^t+\tilde{\mathrm{\Omega }}$ and ${\mathrm{\Sigma }}_{X\xi }=\tilde{\mathrm{\Lambda }}(I-\tilde{B}{)}^{-1}\tilde{\mathrm{\Psi }}(I-\tilde{B}{)}^{-1t}$. Note, that the conditional variance does not depend on $X$, $Z$.

Polynomials

$$\begin{array}{r}\begin{array}{lcl}{\eta }_i& =& \alpha +\sum _{m=1}^k{\beta }_m{\xi }_i^m+{\zeta }_i,\end{array}\end{array}$$

Here $\phi ({\xi }_i)={\xi }_i^m,(m\in \mathbb{N})$ and conditional means are given by

$$\begin{array}{r}\begin{array}{l}\mathbb{E}({\xi }_i^m|X_i,Z_i)=\sum _{k=0}^{[m/2]},m_{x,z}^{m-2k},v^k,\frac{m!}{2^{k}k!(m-2k)!}\end{array}\end{array}$$

Exponentials

For the exponential function $\phi ({\xi }_i)=\exp ({\xi }_i)$ an expression can be obtained as $\exp ({\xi }_i)$ will follow a logarithmic normal distribution where the mean is

$$\begin{array}{r}\mathbb{E}[\exp ({\xi }_i)|X_i,Z_i]=\exp (0.5v+m_{x,y})\end{array}$$

The conditional mean of functions on the form $\phi ({\xi }_i)=\exp ({\xi }_i{)}^m$ is straightforward to calculate as this variable again follows a logarithmic normal distribution.

Interactions

Product-interaction model

$$\begin{array}{rl}{\eta }_i& =\alpha +{\beta }_1{\xi }_{1i}+{\beta }_2{\xi }_{2i}+{\beta }_{12}{\xi }_{1i}{\xi }_{2i}+{\zeta }_i,\end{array}$$

now $\mathbb{E}({\xi }_{1i}{\xi }_{2i}|X_i,Z_i)=\mathbb{C}\mathrm{ov}({\xi }_{1i},{\xi }_{2i}|X_i,Z_i)+\mathbb{E}({\xi }_{1i}|X_i,Z_i)\mathbb{E}({\xi }_{2i}|X_i,Z_i)$, where terms on the right-hand side are directly available from the bivariate normal distribution of ${\xi }_{1i},{\xi }_{2i}$ given $X_i,Z_i$. Regression calibration leads to the correct mean expect that the intercept will be biased as this method will not include the covariance term above.

Splines

A natural cubic spline with $k$ knots $t_1<t_2<\dots <t_k$ is given by

$$\begin{array}{rl}{\eta }_i& =\alpha +{\beta }_0{\xi }_i+\sum _{j=1}^{k-2}{\beta }_j{f}_j({\xi }_i)+{\zeta }_i,\end{array}$$

with $f_j({\xi }_i)=g_j({\xi }_i)-\frac{t_k-t_j}{t_k-t_{k-1}}{g}_{k-1}({\xi }_i)+\frac{t_{k-1}-t_j}{t_k-t_{k-1}}{g}_k({\xi }_i),,j=1,\dots ,k-2$ and $g_j({\xi }_i)=({\xi }_i-t_j{)}^3{1}_{{\xi }_i>t_j},,j=1,\dots ,k$. Thus, predictions $\mathbb{E}[f_j({\xi }_i)|X_i,Z_i]$ are linear functions of $\mathbb{E}[g_j({\xi }_i)|X_i,Z_i]$, and

$$\begin{array}{rl}\mathbb{E}[g_j({\xi }_i)|X_i,Z_i]=& \frac{s}{\sqrt{2\pi }}[(2s^2+(m_{x,z}-t_j{)}^2)\exp (-[\frac{(m_{x,z}-t_j)}{s\sqrt{2}}{]}^2)]+\\ & (m_{x,z}-t_j)[(m_{x,z}-t_j{)}^2+3s^2]p_{x,z,j},\end{array}$$

where $s=\sqrt{v}$ and $p_{x,z,j}=P({\xi }_i>t_j|X_i,Z_i)=1-\mathrm{\Phi }(\frac{t_j-m_{x,z}}{s})$.

Relaxing distributional assumptions

Let $G_i\sim \text{multinom}(\pi )$ be class indicator $G_i\in \{1,\dots ,K\}$, and ${\xi }_i$ the $q$-dim. latent predictor from the mixture

$$\begin{array}{r}{\xi }_i=\sum _{i=1}^{K}I(G_i=k){\xi }_{ki}\end{array}$$

with ${\tilde{\zeta }}_{ki}\sim N(0,{\tilde{\mathrm{\Psi }}}_k)$ where ${\xi }_{ki}={\tilde{\alpha }}_k+\tilde{B}{\xi }_{ki}+\tilde{\mathrm{\Gamma }}{Z}_i+{\tilde{\zeta }}_{ki}$. Assuming $G_i$ to be independent of $({\tilde{\zeta }}_{1i},\dots ,{\tilde{\zeta }}_{Ki})$ and $({\epsilon }_i,{\tilde{\epsilon }}_i)$,

$$\begin{array}{r}\begin{array}{lcl}\mathbb{E}[\phi ({\xi }_i)|X_i,Z_i]& =& \mathbb{E}\{\mathbb{E}[\phi ({\xi }_i)|X_i,Z_i,G_i]|X_i,Z_i\}\\ & =& \sum _{k=1}^{K}P(G_i=k|X_i,Z_i),\mathbb{E}[\phi ({\xi }_{ki})|X_i,Z_i,G_i=k]\\ & =& \sum _{k=1}^{K}P(G_i=k|X_i,Z_i),\mathbb{E}[\phi ({\xi }_{ki})|X_i,Z_i].\end{array}\end{array}$$

Results can be extended also to the case where $\pi$ depends on covariates. Under regularity conditions (Redner 1984) the MLE of the stage 1 mixture model is regular and asymptotic linear (RAL).

Asymptotics

Theorem. Under correctly specified model 2SSEM will yield consistent estimation of all parameters $\theta$ except for the residual covariance, $\mathrm{\Psi }$, of the latent variables in step 2.

Proof: (intuitive version - linear models and Berkson error)

Let $Y=\alpha +\beta X+ϵ$ with Berkson error $X=W+U$ and $\mathbb{C}\mathrm{ov}(W,U)=0$. Let $Y=\alpha +\beta W+\widetilde{ϵ}$. Berkson error does not lead to bias in $\beta$, but residual variance is too high Follows from noting that we have Berkson error in the second step: ${\eta }_i=\alpha +B\phi ({\xi }_i)+\mathrm{\Gamma }{Z}_i+{\zeta }_i$. Iterated conditional means: $\mathbb{C}\mathrm{ov}\{\mathbb{E}\phi (\xi )|X,Z],\mathbb{E}[\phi (\xi )|X,Z]-\phi ({\xi }_i)\}=0$.

Assume that we have i.i.d. observations $(Y_i,X_i,Z_i),i=1,\dots ,n$, and we also restrict attention only to the consistent parameter estimates i.e., ${\theta }_2$ does not contain any of the parameters belonging to $\mathrm{\Psi }$.

We will assume that the stage 1 model estimator is obtained as the solution to the following score equation:

$$\begin{array}{r}{\mathcal{U}}_1({\theta }_1;X,Z)=\sum _{i=1}^n{U}_1({\theta }_1;X_i,Z_i)=0,\end{array}$$

and that

1. The estimator of the stage 1 model is consistent, linear, regular, and asymptotic normal.
2. $\mathcal{U}$ is twice continuous differentiable in a neighbourhood around the true (limiting) parameters $({\theta }_{01}^T,{\theta }_{02}^T{)}^T$. Further, $$n^{-1}\sum _{i=1}^n\mathrm{\nabla }{U}_2(Y_i,X_i,Z_i;{\theta }_1,{\theta }_2)$$ converges uniformly to $\mathbb{E}[\mathrm{\nabla }{U}_2(Y_i,X_i,Z_i;{\theta }_1,{\theta }_2)]$ in a neighbourhood around $({\theta }_{01}^T,{\theta }_{02}^T{)}^T$,
3. and when evaluated here $-\mathbb{E}(\mathrm{\nabla }{U}_2)$ is positive definite.

This implies the following i.i.d. decomposition of the two-stage estimator

$$\begin{array}{rl}\sqrt{n}({\hat{\theta }}_2({\hat{\theta }}_1)-{\theta }_2)& =n^{-\frac{1}{2}}\sum _{i=1}^n{\xi }_2(Y_i,X_i,Z_i;{\theta }_2)\\ & +n^{-\frac{1}{2}}\mathbb{E}[-{\mathrm{\nabla }}_{{\theta }_2}{U}_2(Y,Z,X;{\theta }_2,{\theta }_1){]}^{-1}\\ & \times \mathbb{E}[{\mathrm{\nabla }}_{{\theta }_1}U(Y,Z,X;{\theta }_2,{\theta }_1)]\sum _{i=1}^n{\xi }_1(X_i,Z_i;{\theta }_1)+o_p(1)\\ & =n^{-\frac{1}{2}}\sum _{i=1}^n{\xi }_3(Y_i,X_i,Z_i;{\theta }_2,{\theta }_1)+o_p(1).\end{array}$$

where ${\xi }_1$ and ${\xi }_2$ are the influence functions from the stage 1 and stage 2 models.

Implementation

Installation in R:

  library(lava)
  library(magrittr)

Simulation

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

      # Default: all parameter values are 1. Here we change the covariate effect on eta1
      d <- sim(m, n=200, seed=1, p=c('eta1~z'=-1))

head(d)

          x1         x2         x3       eta1         y1         y2         y3
1 -1.3117148 -0.2758592  0.3891800 -0.6852610  0.6565382  2.8784121  0.1864112
2  1.6733480  3.1785780  3.3853595  1.4897047 -0.9867733  1.9512415  2.7624733
3  0.5661292  2.9883463  0.7987605  1.4017578  0.5694039 -1.2966555 -2.2827075
4  1.7946719 -0.1315167 -0.1914767  0.1993911  0.7221921  0.9447854 -1.3720646
5  0.3702222 -2.2445211 -0.3755076  0.0407144 -0.3144152  0.3546089  0.9828617
6 -2.8786355  0.4394945 -2.4338245 -2.0581671 -4.1310534 -5.6157543 -3.0456611
        eta2           z
1  1.7434470  0.34419403
2  0.8393097  0.01271984
3 -0.4258779 -0.87345013
4  0.7340538  0.34280028
5  0.2852132 -0.17738775
6 -3.9531055  0.92143325

Specification of stage 1

    m1 <- lvm(x1+x2+x3 ~ eta1, eta1 ~ z, latent = ~ eta1)

Specification of stage 2:

      m2 <- lvm() %>%
	  regression(y1+y2+y3 ~ eta2) %>%
	  regression(eta2 ~ z) %>%
	  latent(~ eta2) %>%
	  nonlinear(type="quadratic", eta2 ~ eta1)

Estimation:

(mm <- twostage(m1, m2, data=d))

                    Estimate Std. Error  Z-value   P-value
Measurements:
   y2~eta2           0.98628    0.04612 21.38536    <1e-12
   y3~eta2           0.97614    0.05203 18.76166    <1e-12
Regressions:
   eta2~z            1.08890    0.17735  6.13990 8.257e-10
   eta2~eta1_1       1.13932    0.15548  7.32770    <1e-12
   eta2~eta1_2      -0.38404    0.05770 -6.65582 2.817e-11
Intercepts:
   y2               -0.09581    0.11350 -0.84413    0.3986
   y3                0.01440    0.10849  0.13273    0.8944
   eta2              0.50414    0.17550  2.87264  0.004071
Residual Variances:
   y1                1.27777    0.18980  6.73206
   y2                1.02924    0.13986  7.35895
   y3                0.82589    0.14089  5.86181
   eta2              1.94918    0.26911  7.24305

  pf <- function(p) p["eta2"]+p["eta2~eta1_1"]*u + p["eta2~eta1_2"]*u^2
  plot(mm,f=pf,data=data.frame(u=seq(-2,2,length.out=100)),lwd=2)