Trying to understand a latent curve model in terms of mixed effects regression I'm trying to understand exactly what the following model is trying to represent:

(taken from Beaujean's Latent Variable Modeling Using R book)
The text indicates that this is a random intercept/random slope model, but it is being estimated through nothing but a covariance matrix and a list of means for $X_1 ... X_4$ so I have a hard time visualizing the implied structure. Here's my intuition:
Let the intercept be $I$, the slope $S$ and error $E$::
$$X_j = I_j + (j-1)\times S_j + E_{X1}$$
$$I_j \sim N(\mu_{I_j}, b)\qquad S_j \sim N(\mu_{S_j}, c)\qquad E_{XJ} \sim N(0, \sigma_{E_{XJ}}^2)$$
Ignoring the covariance $a$, would the above specification correctly describe the path model above? 
 A: The algorithms for fitting SEM models are different. SEM models use information from the first two moments, means and variances. Multilevel models use all of the data. To set this up as a multilevel model, the data need to be in long form where each individual repeats across several rows - 4 in OP's example.
Once that is done, then the model is:
\begin{equation}
y = \beta_0 + \beta_1 \times Time + \mathbf{Z}g + \epsilon_t
\end{equation}
where:


*

*$Time = \{0,1,2,3\}$ representing the fixed loadings

*$\mathbf{Z}$ is the model matrix for the random intercept and random slope on the $Time$ variable

*$g \sim \mathcal{N}_2(\mathbf{0}, \mathbf{\Sigma})$. These are the random intercept and slope. $\mathbf{\Sigma}$ is a variance-covariance matrix with $b$ and $c$ on the diagonal representing the variances of the random intercept and slope respectively, and $a$ on the off-diagonal representing their covariance.

*$\epsilon_t$ is heteroskedastic and varies by time.


This is a model you can easily fit with nlme. Take the example on this UCLA page. The corresponding syntax to fit the model would be:
dat <- read.table(
  "https://stats.idre.ucla.edu/wp-content/uploads/2016/03/ex6.1.dat")
head(dat)
#          V1       V2       V3       V4
# 1  0.036879 1.473688 1.683264 1.949616
# 2 -2.692610 1.658458 2.212524 4.025853
# 3  2.753869 5.125832 5.201779 6.776908
names(dat) <- paste0("t", 0:3)
dat$ID <- 1:nrow(dat)
dat <- tidyr::gather(dat, time, y, t0:t3) # wide to long form
dat$time <- as.integer(gsub("t", "", dat$time))
library(nlme)
fit <- lme(
  y ~ time, dat, ~ 1 + time | ID, weights = varIdent(form = ~ 1 | time),
  method = "ML")
logLik(fit)
# 'log Lik.' -3016.386 (df=9)

This is the same log likelihood on the UCLA page.
