# Using least squares to estimate variance of latent variable

I am having trouble understanding why I can't use least squares to solve an overdetermined system of linear equations using $$\bf{x} = (\bf{A}'\bf{A})^{-1} \bf{A}'\bf{b}$$. The same model estimated using maximimum likelihood in structural equation modeling works fine.

I want to understand analytically how a latent individual effects variable ($$\alpha$$) is being estimated in a "simple" three-wave autoregressive model with unobserved heterogeneity. I want to use least squares rather than maximum likelihood to make it more accessible to those unfamiliar with structural equation modeling.

I simulated a sample covariance matrix using lavaan based on the model:

\begin{align} x_{1} & = \alpha + \delta_{1} & \\ x_{t} & = \alpha + \varphi x_{t-1} + \delta_{t} & \ \text{for} \ t = 2,3 \end{align} My lavaan code for generating the observed covariance matrix is:

library( lavaan)
popx <- '
x1 ~ 1*a
x2 ~ 1*a + 0.2*x1
x3 ~ 1*a + 0.2*x2
'
simdfx <- simulateData( popx, sample.nobs = 200, empirical = TRUE)
cov( simdfx)


The empirical = TRUE option means I don't have to worry about sampling error.

I am running an Autoregressive Latent Trajectory (ALT) model with no slope factor and predetermined first observation $$x_{1}$$, i.e. unconstrained $$Cov(x_{1},\alpha)$$. The model is estimated in lavaan with the following code:

m1 <- '
# Individual effects
a =~ 1*x2 + 1*x3
# Regressions, AR parameter fixed over time
x2 ~ p*x1
x3 ~ p*x2
# Covariance
a ~~ x1
# Fixed error variance for sake of simplicity
x2 ~~ u*x2
x3 ~~ u*x3
'
summary( m1.fit <- sem( m1, simdfx))
inspect( m1.fit, 'est')


I then use the package Ryacas (computer algebra software for R) to double-check the symbolic model-implied covariance matrix $$\bf{\Sigma}(\bf{\theta})$$. I equate the six unique model-implied equations to the observed covariance matrix $$\bf{S}$$ and rewrite in $$\bf{A}\bf{x} = \bf{b}$$ format to try and solve for $$Var(\alpha)$$. There are four unknowns: $$Var(\alpha)$$, $$Cov(\alpha,x_{1})$$, $$Var(\delta_{1})$$, $$Var(\delta_{2})$$. My $$m \times n$$ matrix $$\bf{A}$$ is overdetermined with six rows, and five columns. I substituted the true value of $$\varphi = 0.2$$ in $$\bf{x}$$ because I thought that would make things a little simpler, plus right now I just want to see what is going on in the estimation of $$Var(\alpha)$$.

My system of equations are:

\begin{alignat}{2} & x_{1}^{2} & = x_{1}x_{1} & = 2.01 \\ & 0.2x_{1}^2 + x_{1}\alpha & = x_{1}x_{2} & = 1.41\\ & 0.04x_{1}^2 + 1.2x_{1}\alpha & = x_{1}x_{3} & = 1.29 \\ & 0.04x_{1}^2 + 0.4x_{1}\alpha + \alpha^2 + \delta_{2}^2 & = x_{2}x_{2} & = 2.49 \\ & 0.008x_{1}^2 + 0.28x_{1}\alpha + 1.2\alpha^2 + 0.28\delta_{2}^2 & = x_{2}x_{3} & = 1.70 \\ & 0.0016x_{1}^2 + 0.096x_{1}\alpha + 1.44\alpha^2 + 0.04\delta_{2}^2 + \delta_{3}^2 & = x_{3}x_{3} & = 2.59 \end{alignat}

which I rewrote as

\begin{align} \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0.2 & 1 & 0 & 0 & 0 \\ 0.04 & 1.2 & 0 & 0 & 0 \\ 0.04 & 0.4 & 1 & 1 & 0 \\ 0.008 & 0.28 & 1.2 & 0.28 & 0 \\ 0.0016 & 0.096 & 1.44 & 0.04 & 1 \end{bmatrix} \begin{bmatrix} x_{1}^2 \\ x_{1}\alpha \\ \alpha^2 \\ \delta_{2}^2 \\ \delta_{3}^2 \end{bmatrix} = \begin{bmatrix} 2.01 \\ 1.41 \\ 1.29 \\ 2.49 \\ 1.70 \\ 2.59 \end{bmatrix}. \end{align}

Now I am trying to solve for $$\bf{x}$$ with $$\bf{x} = (\bf{A}'\bf{A})^{-1}\bf{A}'\bf{b}$$ and keep coming up with a vector $$\bf{x}' = (\frac{\infty}{\infty}, \frac{\infty}{\infty}, \frac{\infty}{\infty}, \frac{\infty}{\infty}, \frac{\infty}{\infty})$$.

It seems like the elements of the inverse of $$\bf{A}'\bf{A}$$ are extremely small which means I am basically dividing $$\bf{A}'\bf{b}$$ by zero. On the other hand, $$\bf{A}'\bf{A}$$ looks to be invertable because (at least) the determinant doesn't seem to be zero. I tried dropping the fifth column for $$Var(\delta_{3})$$ since the error variances are constrained to be equal in the SEM model but this didn't seem to help. I'm thinking now that there might be an issue with multicolinearity because due to the autoregression, I suppose everything could be rewritten as a function of $$Var(\alpha)$$.

However, my SEM model works fine and gives me correct estimates of $$\widehat{Var}(x_{1}) = 2.0$$, $$\widehat{Var}(\alpha) = 1.0$$, $$\widehat{Cov}(x_{1},\alpha) = 1.0$$, $$\widehat{Var}(\delta_{2}) = \widehat{Var}(\delta_{3}) = 1.0$$.

I used the same procedure to come up with a least squares estimate for $$Var(\alpha)$$ in an intercept-only model which worked out fine. I would appreciate if anyone could point me in the right direction and let me know what I am obviously overlooking.

• Sorry, this always happens. I spend ages trying to figure something out, then as soon as I ask I end up figuring it out myself. Or at least I think I know what was wrong: the variance of x1 is known, so I first subtracted it from b and eliminated that column along with the first row (because I don't think it provides any new information). Then I tried the same thing with the new matrix equation and now I am getting reasonable least squares estimates: Cov(x1,a) = 1.00, Var(a) = 0.91, Var(d) = 1.11. – hendogg87 Aug 15 at 9:43