## Overview

This tutorial walks through the fitting of growth models with nonlinearity in random coefficients in several different frameworks (e.g., multilevel modeling framework, structural equation modeling framework), and demonstrates these models using different R packages (knowing how to fit the models in different packages can be helpful when trying to fit more complex models as each packages as its own advantages and limitations). In this tutorial, we will be using a sample data set that includes repeated measures of individuals’ height during infancy and early childhood.

The code and example provided in this tutorial are from Chapter 12 of Grimm, Ram, and Estabrook (2016), with a few additions in code and commentary. Please refer to the chapter for further interpretations and insights about the analyses.

## Outline

This tutorial provides line-by-line code to examine growth models with nonlineary in parameters using the:

1. Jenss-Bayley Growth Model,
• Direct Optimization -“nlme” library -“lavaan” library
• First Order Approximation -“nlme” library
2. Bilinear Spline Growth Model with Variation in the Knot Point
• “nlme” library
• “lavaan” library

## Step 0: Read in the data and call needed libraries.

library(nlme)
library(lavaan)
## This is lavaan 0.5-20
## lavaan is BETA software! Please report any bugs.
#set filepath for long data file
filepath.long<-"http://quant-d.ssri.psu.edu/sites/qdev/files/bgs_height_long_0.csv"
#set filepath for wide data file
filepath.wide<-"http://quant-d.ssri.psu.edu/sites/qdev/files/bgs_height_wide.csv"

#read in the .csv files using the url() function
head(hght_long)
##   id age hght
## 1  1   1 53.5
## 2  1   3 60.4
## 3  1   6 67.2
## 4  1   9 70.9
## 5  1  12 76.2
## 6  1  15 81.1
hght_wide <- read.csv(file=url(filepath.wide),header=TRUE,na.string='.')
head(hght_wide)
##   id hght01 hght03 hght06 hght09 hght12 hght15 hght18 hght24 hght36
## 1  1   53.5   60.4   67.2   70.9   76.2   81.1     NA   89.5   95.9
## 2  2   54.5   58.8   66.0   69.1   74.5   78.2   81.2   86.7     NA
## 3  3   49.5   56.0   61.9   64.1   67.9   69.6   73.8     NA   85.1
## 4  4   56.9   64.6   70.0   72.5   76.1   80.4   82.3   87.9   94.3
## 5  5   56.7   63.1     NA   73.1   77.3   81.8   83.5   88.9   99.1
## 6  6   56.1   61.9   67.9   72.3   77.3   79.9   82.3   88.3   97.5

## Step 1: Jenss-Bayley Growth Model.

### Direct Optimization.

#### Fitting the model using the nlme library.

hght.jb.nlme <- nlme(hght~b_1i + b_2i * (age/12) + b_3i * (exp(b_4i*(age/12))-1),
data = hght_long,
fixed = b_1i + b_2i + b_3i + b_4i ~ 1,
random = b_1i + b_2i + b_3i + b_4i ~ 1,
groups = ~ id,
start = c(50, 10, -18, -2),
na.action = na.exclude)
summary(hght.jb.nlme)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: hght ~ b_1i + b_2i * (age/12) + b_3i * (exp(b_4i * (age/12)) -      1)
##  Data: hght_long
##        AIC      BIC   logLik
##   2005.178 2070.649 -987.589
##
## Random effects:
##  Formula: list(b_1i ~ 1, b_2i ~ 1, b_3i ~ 1, b_4i ~ 1)
##  Level: id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev    Corr
## b_1i     2.6774685 b_1i   b_2i   b_3i
## b_2i     0.8717526  0.319
## b_3i     3.4099947  0.491  0.449
## b_4i     0.2746913 -0.120 -0.275 -0.642
## Residual 0.8041717
##
## Fixed effects: b_1i + b_2i + b_3i + b_4i ~ 1
##          Value Std.Error  DF   t-value p-value
## b_1i  50.99063 0.3422353 495 148.99288       0
## b_2i   9.31946 0.1607672 495  57.96865       0
## b_3i -17.84870 0.4814401 495 -37.07357       0
## b_4i  -2.13255 0.0769100 495 -27.72791       0
##  Correlation:
##      b_1i   b_2i   b_3i
## b_2i  0.022
## b_3i  0.376  0.607
## b_4i  0.263 -0.642 -0.536
##
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max
## -3.482492379 -0.518214899  0.001140808  0.517854147  2.745161193
##
## Number of Observations: 581
## Number of Groups: 83

#### Fitting the model using the lavaan library.

jb.hght.lavaan <- '
#latent variable definitions

#defining the intercept
eta_1 =~ 1*hght01 +
1*hght03 +
1*hght06 +
1*hght09 +
1*hght12 +
1*hght15 +
1*hght18 +
1*hght24 +
1*hght36

#defining the slope
eta_2 =~ .0833*hght01 +
.25*hght03 +
.5*hght06 +
.75*hght09 +
1*hght12 +
1.25*hght15 +
1.5*hght18 +
2*hght24 +
3*hght36

#defining eta_3
eta_3 =~ start(-.15)*L13*hght01 +
start(-.40)*L23*hght03 +
start(-.64)*L33*hght06 +
start(-.78)*L43*hght09 +
start(-.87)*L53*hght12 +
start(-.92)*L63*hght15 +
start(-.95)*L73*hght18 +
start(-.98)*L83*hght24 +
start(-.99)*L93*hght36

#defining eta_4
eta_4 =~ start(-1.25)*L14*hght01 +
start(-2.65)*L24*hght03 +
start(-3.16)*L34*hght06 +
start(-2.82)*L44*hght09 +
start(-2.24)*L54*hght12 +
start(-1.66)*L64*hght15 +
start(-1.19)*L74*hght18 +
start(-0.56)*L84*hght24 +
start(-0.11)*L94*hght36

#latent means
eta_1 ~ start(51)*1
eta_2 ~ start(9.25)*1
eta_3 ~ start(-17.8)*alpha3*1
eta_4 ~ 0 # start(0.07)*1

#introducing phantom factor, setting all parameters to 0 and mean to alpha4
phantom =~ 0*hght01
phantom ~~ 0*phantom
phantom ~~ 0*eta_1
phantom ~~ 0*eta_2
phantom ~~ 0*eta_3
phantom ~~ 0*eta_4
phantom ~ start(-2)*alpha4*1

alpha4 < -2   #add this constraint in so that it finds the viable solution
L13 == exp(alpha4*0.0833)-1
L23 == exp(alpha4*0.25)-1
L33 == exp(alpha4*0.50)-1
L43 == exp(alpha4*0.75)-1
L53 == exp(alpha4*1.00)-1
L63 == exp(alpha4*1.25)-1
L73 == exp(alpha4*1.50)-1
L83 == exp(alpha4*2.00)-1
L93 == exp(alpha4*3.00)-1

L14 == alpha3 * 0.0833 * exp(alpha4*0.0833)-1
L24 == alpha3 * 0.25 * exp(alpha4*0.25)-1
L34 == alpha3 * 0.50 * exp(alpha4*0.50)-1
L44 == alpha3 * 0.75 * exp(alpha4*0.75)-1
L54 == alpha3 * 1.00 * exp(alpha4*1.00)-1
L64 == alpha3 * 1.25 * exp(alpha4*1.25)-1
L74 == alpha3 * 1.50 * exp(alpha4*1.50)-1
L84 == alpha3 * 2.00 * exp(alpha4*2.00)-1
L94 == alpha3 * 3.00 * exp(alpha4*3.00)-1

#factor variances
eta_1 ~~ start(7.4)*eta_1
eta_2 ~~ start(0.81)*eta_2
eta_3 ~~ start(11.63)*eta_3
eta_4 ~~ start(0.08)*eta_4

#factor covariances
eta_1 ~~ start(0.64)*eta_2
eta_1 ~~ start(4.45)*eta_3
eta_1 ~~ start(0.03)*eta_4

eta_2 ~~ start(1.35)*eta_3
eta_2 ~~ start(-0.08)*eta_4

eta_3 ~~ start(-0.49)*eta_4

#manifest variances
hght01 ~~ start(.64)*theta*hght01
hght03 ~~ theta*hght03
hght06 ~~ theta*hght06
hght09 ~~ theta*hght09
hght12 ~~ theta*hght12
hght15 ~~ theta*hght15
hght18 ~~ theta*hght18
hght24 ~~ theta*hght24
hght36 ~~ theta*hght36

#manifest means (fixed to zero)
hght01 ~ 0*1
hght03 ~ 0*1
hght06 ~ 0*1
hght09 ~ 0*1
hght12 ~ 0*1
hght15 ~ 0*1
hght18 ~ 0*1
hght24 ~ 0*1
hght36 ~ 0*1'

lavaan.fit <- lavaan(jb.hght.lavaan,
data = hght_wide,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE,
mimic="mplus",
control=list(iter.max=500),
verbose=FALSE)
summary(lavaan.fit, fit.measures=TRUE)
## lavaan (0.5-20) converged normally after 1829 iterations
##
##   Number of observations                            83
##
##   Number of missing patterns                        31
##
##   Estimator                                         ML
##   Minimum Function Test Statistic              103.452
##   Degrees of freedom                                39
##   P-value (Chi-square)                           0.000
##
## Model test baseline model:
##
##   Minimum Function Test Statistic              955.844
##   Degrees of freedom                                36
##   P-value                                        0.000
##
## User model versus baseline model:
##
##   Comparative Fit Index (CFI)                    0.930
##   Tucker-Lewis Index (TLI)                       0.935
##
## Loglikelihood and Information Criteria:
##
##   Loglikelihood user model (H0)               -987.467
##   Loglikelihood unrestricted model (H1)       -935.741
##
##   Number of free parameters                         15
##   Akaike (AIC)                                2004.935
##   Bayesian (BIC)                              2041.217
##   Sample-size adjusted Bayesian (BIC)         1993.904
##
## Root Mean Square Error of Approximation:
##
##   RMSEA                                          0.141
##   90 Percent Confidence Interval          0.108  0.174
##   P-value RMSEA <= 0.05                          0.000
##
## Standardized Root Mean Square Residual:
##
##   SRMR                                           0.306
##
## Parameter Estimates:
##
##   Information                                 Observed
##   Standard Errors                             Standard
##
## Latent Variables:
##                    Estimate  Std.Err  Z-value   P(>|z|)
##   eta_1 =~
##     hght01            1.000
##     hght03            1.000
##     hght06            1.000
##     hght09            1.000
##     hght12            1.000
##     hght15            1.000
##     hght18            1.000
##     hght24            1.000
##     hght36            1.000
##   eta_2 =~
##     hght01            0.083
##     hght03            0.250
##     hght06            0.500
##     hght09            0.750
##     hght12            1.000
##     hght15            1.250
##     hght18            1.500
##     hght24            2.000
##     hght36            3.000
##   eta_3 =~
##     hght01   (L13)   -0.159    0.005   -28.938    0.000
##     hght03   (L23)   -0.405    0.012   -34.751    0.000
##     hght06   (L33)   -0.646    0.014   -46.575    0.000
##     hght09   (L43)   -0.789    0.012   -63.764    0.000
##     hght12   (L53)   -0.875    0.010   -89.055    0.000
##     hght15   (L63)   -0.925    0.007  -126.678    0.000
##     hght18   (L73)   -0.956    0.005  -183.196    0.000
##     hght24   (L83)   -0.984    0.002  -399.676    0.000
##     hght36   (L93)   -0.998    0.000 -2154.873    0.000
##   eta_4 =~
##     hght01   (L14)   -2.250    0.038   -58.536    0.000
##     hght03   (L24)   -3.653    0.107   -34.037    0.000
##     hght06   (L34)   -4.158    0.182   -22.881    0.000
##     hght09   (L44)   -3.819    0.214   -17.862    0.000
##     hght12   (L54)   -3.236    0.212   -15.287    0.000
##     hght15   (L64)   -2.663    0.189   -14.077    0.000
##     hght18   (L74)   -2.188    0.158   -13.852    0.000
##     hght24   (L84)   -1.561    0.096   -16.214    0.000
##     hght36   (L94)   -1.105    0.026   -42.016    0.000
##   phantom =~
##     hght01            0.000
##
## Covariances:
##                    Estimate  Std.Err  Z-value   P(>|z|)
##   eta_1 ~~
##     phantom           0.000
##   eta_2 ~~
##     phantom           0.000
##   eta_3 ~~
##     phantom           0.000
##   eta_4 ~~
##     phantom           0.000
##   eta_1 ~~
##     eta_2             0.561    0.585     0.960    0.337
##     eta_3             3.960    1.641     2.413    0.016
##     eta_4             0.105    0.304     0.347    0.729
##   eta_2 ~~
##     eta_3             1.346    0.849     1.586    0.113
##     eta_4            -0.078    0.143    -0.542    0.588
##   eta_3 ~~
##     eta_4            -0.488    0.383    -1.274    0.203
##
## Intercepts:
##                    Estimate  Std.Err  Z-value   P(>|z|)
##     eta_1            51.076    0.349   146.465    0.000
##     eta_2             9.296    0.168    55.402    0.000
##     eta_3   (alp3)  -17.836    0.481   -37.051    0.000
##     eta_4             0.000
##     phantom (alp4)   -2.076    0.078   -26.508    0.000
##     hght01            0.000
##     hght03            0.000
##     hght06            0.000
##     hght09            0.000
##     hght12            0.000
##     hght15            0.000
##     hght18            0.000
##     hght24            0.000
##     hght36            0.000
##
## Variances:
##                    Estimate  Std.Err  Z-value   P(>|z|)
##     phantom           0.000
##     eta_1             7.545    1.828     4.126    0.000
##     eta_2             0.809    0.340     2.380    0.017
##     eta_3            11.626    3.006     3.868    0.000
##     eta_4             0.076    0.079     0.959    0.338
##     hght01  (thet)    0.644    0.053    12.175    0.000
##     hght03  (thet)    0.644    0.053    12.175    0.000
##     hght06  (thet)    0.644    0.053    12.175    0.000
##     hght09  (thet)    0.644    0.053    12.175    0.000
##     hght12  (thet)    0.644    0.053    12.175    0.000
##     hght15  (thet)    0.644    0.053    12.175    0.000
##     hght18  (thet)    0.644    0.053    12.175    0.000
##     hght24  (thet)    0.644    0.053    12.175    0.000
##     hght36  (thet)    0.644    0.053    12.175    0.000
##
## Constraints:
##                                                |Slack|
##     -2 - (alpha4)                                0.076
##     L13 - (exp(alpha4*0.0833)-1)                 0.000
##     L23 - (exp(alpha4*0.25)-1)                   0.000
##     L33 - (exp(alpha4*0.50)-1)                   0.000
##     L43 - (exp(alpha4*0.75)-1)                   0.000
##     L53 - (exp(alpha4*1.00)-1)                   0.000
##     L63 - (exp(alpha4*1.25)-1)                   0.000
##     L73 - (exp(alpha4*1.50)-1)                   0.000
##     L83 - (exp(alpha4*2.00)-1)                   0.000
##     L93 - (exp(alpha4*3.00)-1)                   0.000
##     L14-(alpha3*0.0833*exp(alpha4*0.0833)-1)     0.000
##     L24 - (alpha3*0.25*exp(alpha4*0.25)-1)       0.000
##     L34 - (alpha3*0.50*exp(alpha4*0.50)-1)       0.000
##     L44 - (alpha3*0.75*exp(alpha4*0.75)-1)       0.000
##     L54 - (alpha3*1.00*exp(alpha4*1.00)-1)       0.000
##     L64 - (alpha3*1.25*exp(alpha4*1.25)-1)       0.000
##     L74 - (alpha3*1.50*exp(alpha4*1.50)-1)       0.000
##     L84 - (alpha3*2.00*exp(alpha4*2.00)-1)       0.000
##     L94 - (alpha3*3.00*exp(alpha4*3.00)-1)       0.000

### First Order Approximation.

#### Fitting the model using the nlme library.

hght.jb.nlme.fo <- nlme(hght ~ (b_1 + d_1i) + (b_2 + d_2i) * (age/12) +
(b_3 + d_3i) * (exp(b_4 * (age/12))-1) +
d_4i * (b_3 * (age/12) * exp(b_4 * (age/12))),
data = hght_long,
fixed = b_1 + b_2 + b_3 + b_4 ~ 1,
random = d_1i + d_2i + d_3i + d_4i ~ 1,
groups =~ id,
start = c(50, 10, -18, -2),
na.action = na.exclude)
summary(hght.jb.nlme.fo)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: hght ~ (b_1 + d_1i) + (b_2 + d_2i) * (age/12) + (b_3 + d_3i) *      (exp(b_4 * (age/12)) - 1) + d_4i * (b_3 * (age/12) * exp(b_4 *      (age/12)))
##  Data: hght_long
##        AIC      BIC    logLik
##   2005.213 2070.684 -987.6063
##
## Random effects:
##  Formula: list(d_1i ~ 1, d_2i ~ 1, d_3i ~ 1, d_4i ~ 1)
##  Level: id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev    Corr
## d_1i     2.7167272 d_1i   d_2i   d_3i
## d_2i     0.8928535  0.267
## d_3i     3.3840109  0.482  0.433
## d_4i     0.2775674  0.027 -0.299 -0.525
## Residual 0.8027907
##
## Fixed effects: b_1 + b_2 + b_3 + b_4 ~ 1
##         Value Std.Error  DF   t-value p-value
## b_1  50.96848 0.3478280 495 146.53357       0
## b_2   9.33387 0.1602774 495  58.23568       0
## b_3 -17.83233 0.4667397 495 -38.20615       0
## b_4  -2.11293 0.0752805 495 -28.06740       0
##  Correlation:
##     b_1    b_2    b_3
## b_2 -0.006
## b_3  0.385  0.578
## b_4  0.326 -0.634 -0.462
##
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max
## -3.43001311 -0.51726171 -0.02358487  0.52150038  2.78477882
##
## Number of Observations: 581
## Number of Groups: 83

## Step 2: Bilinear Spline Growth Model with Variation in the Knot Point.

### Fitting the model using the nlme library.

hght.spline.nlme <- nlme(hght ~ b_1i + b_2i * (pmin(0, age - b_4i)) +
b_3i * (pmax(0, age - b_4i)),
data = hght_long,
fixed = b_1i + b_2i + b_3i + b_4i ~ 1,
random = b_1i + b_2i + b_3i + b_4i ~ 1,
groups =~ id,
start = c(60, 5, 2, 8),
na.action = na.omit)
summary(hght.spline.nlme) 
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: hght ~ b_1i + b_2i * (pmin(0, age - b_4i)) + b_3i * (pmax(0,      age - b_4i))
##  Data: hght_long
##        AIC      BIC    logLik
##   2293.053 2358.524 -1131.526
##
## Random effects:
##  Formula: list(b_1i ~ 1, b_2i ~ 1, b_3i ~ 1, b_4i ~ 1)
##  Level: id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev     Corr
## b_1i     2.70810506 b_1i   b_2i   b_3i
## b_2i     0.28603962  0.107
## b_3i     0.05567495  0.475  0.444
## b_4i     0.29364264  0.792 -0.506  0.235
## Residual 1.21258074
##
## Fixed effects: b_1i + b_2i + b_3i + b_4i ~ 1
##         Value Std.Error  DF   t-value p-value
## b_1i 71.92542 0.3602927 495 199.63051       0
## b_2i  2.51072 0.0540711 495  46.43362       0
## b_3i  0.90033 0.0098067 495  91.80777       0
## b_4i  7.72086 0.1427460 495  54.08808       0
##  Correlation:
##      b_1i   b_2i   b_3i
## b_2i -0.153
## b_3i  0.028  0.168
## b_4i  0.617 -0.685 -0.190
##
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max
## -2.639119082 -0.529769751  0.009075796  0.574138864  2.623652491
##
## Number of Observations: 581
## Number of Groups: 83

### Fitting the model using the lavaan library.

Note: this model results in estimation issues (e.g., variance of the knot point was estimated to be negative), thus the model results should not be interpreted.

spline.hght.lavaan <- '
#latent variable definitions

#defining the intercept
eta_1 =~ 1*hght01 +
1*hght03 +
1*hght06 +
1*hght09 +
1*hght12 +
1*hght15 +
1*hght18 +
1*hght24 +
1*hght36

#defining the slope
eta_2 =~ 1*L12*hght01 +
3*L22*hght03 +
6*L32*hght06 +
9*L42*hght09 +
12*L52*hght12 +
15*L62*hght15 +
18*L72*hght18 +
24*L82*hght24 +
36*L92*hght36

#defining eta_3
eta_3 =~ NA*L13*hght01 +
L23*hght03 +
L33*hght06 +
L43*hght09 +
L53*hght12 +
L63*hght15 +
L73*hght18 +
L83*hght24 +
L93*hght36

#defining eta_4
eta_4 =~ start(-0.90)*L14*hght01 +
start(-0.90)*L24*hght03 +
start(-0.90)*L34*hght06 +
start(-2.57)*L44*hght09 +
start(-2.57)*L54*hght12 +
start(-2.57)*L64*hght15 +
start(-2.57)*L74*hght18 +
start(-2.57)*L84*hght24 +
start(-2.57)*L94*hght36

#latent means
eta_1 ~ start(71)*1
eta_2 ~ start(2.5)*alpha2*1
eta_3 ~ start(.09)*alpha3*1
eta_4 ~ start(7.7)*alpha4*1

#factor variances
#fixing variances to zero b/c getting negative estimates
eta_1 ~~ start(7.3)*eta_1
eta_2 ~~ 0*eta_2
eta_3 ~~ start(0)*eta_3
eta_4 ~~ start(0.08)*eta_4

#factor covariances
eta_1 ~~ start(.08)*eta_2
eta_1 ~~ 0*eta_3
eta_1 ~~ start(.03)*eta_4

eta_2 ~~ 0*eta_3
eta_2 ~~ start(.40)*eta_4

eta_3 ~~ 0*eta_4

#introducing phantom factors,  setting all parameters to 0 and mean to alpha4
phantom2 =~ 0*hght01
phantom2 ~~ 0*phantom2
phantom2 ~~ 0*eta_1
phantom2 ~~ 0*eta_2
phantom2 ~~ 0*eta_3
phantom2 ~~ 0*eta_4
phantom2 ~ start(.29)*alpha2*1

phantom3 =~ 0*hght01
phantom3 ~~ 0*phantom3
phantom3 ~~ 0*phantom2
phantom3 ~~ 0*eta_1
phantom3 ~~ 0*eta_2
phantom3 ~~ 0*eta_3
phantom3 ~~ 0*eta_4
phantom3 ~ start(.29)*alpha3*1

phantom4 =~ 0*hght01
phantom4 ~~ 0*phantom4
phantom4 ~~ 0*phantom3
phantom4 ~~ 0*phantom2
phantom4 ~~ 0*eta_1
phantom4 ~~ 0*eta_2
phantom4 ~~ 0*eta_3
phantom4 ~~ 0*eta_4
phantom4 ~ start(.29)*alpha4*1

L12 ==  1 - alpha4
L22 ==  3 - alpha4
L32 ==  6 - alpha4
L42 ==  9 - alpha4
L52 == 12 - alpha4
L62 == 15 - alpha4
L72 == 18 - alpha4
L82 == 24 - alpha4
L92 == 36 - alpha4

L13 == sqrt(( 1 - alpha4)^2)
L23 == sqrt(( 3 - alpha4)^2)
L33 == sqrt(( 6 - alpha4)^2)
L43 == sqrt(( 9 - alpha4)^2)
L53 == sqrt((12 - alpha4)^2)
L63 == sqrt((15 - alpha4)^2)
L73 == sqrt((18 - alpha4)^2)
L83 == sqrt((24 - alpha4)^2)
L93 == sqrt((36 - alpha4)^2)

L14 == (-alpha2 * sqrt(( 1 - alpha4)^2) + alpha3 * 1 - alpha3 * alpha4) / sqrt(( 1 - alpha4)^2)
L24 == (-alpha2 * sqrt(( 3 - alpha4)^2) + alpha3 * 3 - alpha3 * alpha4) / sqrt(( 3 - alpha4)^2)
L34 == (-alpha2 * sqrt(( 6 - alpha4)^2) + alpha3 * 6 - alpha3 * alpha4) / sqrt(( 6 - alpha4)^2)
L44 == (-alpha2 * sqrt(( 9 - alpha4)^2) + alpha3 * 9 - alpha3 * alpha4) / sqrt(( 9 - alpha4)^2)
L54 == (-alpha2 * sqrt((12 - alpha4)^2) + alpha3 * 12 - alpha3 * alpha4) / sqrt((12 - alpha4)^2)
L64 == (-alpha2 * sqrt((15 - alpha4)^2) + alpha3 * 15 - alpha3 * alpha4) / sqrt((15 - alpha4)^2)
L74 == (-alpha2 * sqrt((18 - alpha4)^2) + alpha3 * 18 - alpha3 * alpha4) / sqrt((18 - alpha4)^2)
L84 == (-alpha2 * sqrt((24 - alpha4)^2) + alpha3 * 24 - alpha3 * alpha4) / sqrt((24 - alpha4)^2)
L94 == (-alpha2 * sqrt((36 - alpha4)^2) + alpha3 * 36 - alpha3 * alpha4) / sqrt((36 - alpha4)^2)

#manifest variances
hght01 ~~ start(1.2)*theta*hght01
hght03 ~~ theta*hght03
hght06 ~~ theta*hght06
hght09 ~~ theta*hght09
hght12 ~~ theta*hght12
hght15 ~~ theta*hght15
hght18 ~~ theta*hght18
hght24 ~~ theta*hght24
hght36 ~~ theta*hght36

#manifest means (fixed to zero)
hght01 ~ 0*1
hght03 ~ 0*1
hght06 ~ 0*1
hght09 ~ 0*1
hght12 ~ 0*1
hght15 ~ 0*1
hght18 ~ 0*1
hght24 ~ 0*1
hght36 ~ 0*1'

lavaan.fit1 <- lavaan(spline.hght.lavaan,
data = hght_wide,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE)
## Warning in lavaan(spline.hght.lavaan, data = hght_wide, meanstructure =
## TRUE, : lavaan WARNING: model has NOT converged!
summary(lavaan.fit1, fit.measures=TRUE)
## ** WARNING ** lavaan (0.5-20) did NOT converge after 668 iterations
## ** WARNING ** Estimates below are most likely unreliable
##
##   Number of observations                            83
##
##   Number of missing patterns                        31
##
##   Estimator                                         ML
##   Minimum Function Test Statistic                   NA
##   Degrees of freedom                                NA
##   P-value                                           NA
## Warning in .local(object, ...): lavaan WARNING: fit measures not available if model did not converge
##
## Parameter Estimates:
##
##   Information                                 Observed
##   Standard Errors                             Standard
##
## Latent Variables:
##                    Estimate  Std.Err  Z-value  P(>|z|)
##   eta_1 =~
##     hght01            1.000
##     hght03            1.000
##     hght06            1.000
##     hght09            1.000
##     hght12            1.000
##     hght15            1.000
##     hght18            1.000
##     hght24            1.000
##     hght36            1.000
##   eta_2 =~
##     hght01   (L12)  -13.998       NA
##     hght03   (L22)  -12.000       NA
##     hght06   (L32)   -9.000       NA
##     hght09   (L42)   -6.000       NA
##     hght12   (L52)   -3.000       NA
##     hght15   (L62)    0.000       NA
##     hght18   (L72)    2.999       NA
##     hght24   (L82)    9.000       NA
##     hght36   (L92)   21.001       NA
##   eta_3 =~
##     hght01   (L13)   14.000       NA
##     hght03   (L23)   12.000       NA
##     hght06   (L33)    9.000       NA
##     hght09   (L43)    6.000       NA
##     hght12   (L53)    3.000       NA
##     hght15   (L63)    0.000       NA
##     hght18   (L73)    3.000       NA
##     hght24   (L83)    9.000       NA
##     hght36   (L93)   21.000       NA
##   eta_4 =~
##     hght01   (L14)   -1.197       NA
##     hght03   (L24)   -1.186       NA
##     hght06   (L34)   -1.183       NA
##     hght09   (L44)   -1.191       NA
##     hght12   (L54)   -1.199       NA
##     hght15   (L64)   -1.200       NA
##     hght18   (L74)   -1.163       NA
##     hght24   (L84)   -1.152       NA
##     hght36   (L94)   -1.134       NA
##   phantom2 =~
##     hght01            0.000
##   phantom3 =~
##     hght01            0.000
##   phantom4 =~
##     hght01            0.000
##
## Covariances:
##                    Estimate  Std.Err  Z-value  P(>|z|)
##   eta_1 ~~
##     eta_2             0.529       NA
##     eta_3             0.000
##     eta_4             2.512       NA
##   eta_2 ~~
##     eta_3             0.000
##     eta_4             0.418       NA
##   eta_3 ~~
##     eta_4             0.000
##   eta_1 ~~
##     phantom2          0.000
##   eta_2 ~~
##     phantom2          0.000
##   eta_3 ~~
##     phantom2          0.000
##   eta_4 ~~
##     phantom2          0.000
##   phantom2 ~~
##     phantom3          0.000
##   eta_1 ~~
##     phantom3          0.000
##   eta_2 ~~
##     phantom3          0.000
##   eta_3 ~~
##     phantom3          0.000
##   eta_4 ~~
##     phantom3          0.000
##   phantom3 ~~
##     phantom4          0.000
##   phantom2 ~~
##     phantom4          0.000
##   eta_1 ~~
##     phantom4          0.000
##   eta_2 ~~
##     phantom4          0.000
##   eta_3 ~~
##     phantom4          0.000
##   eta_4 ~~
##     phantom4          0.000
##
## Intercepts:
##                    Estimate  Std.Err  Z-value  P(>|z|)
##     eta_1            98.002       NA
##     eta_2   (alp2)    1.170       NA
##     eta_3   (alp3)    0.026       NA
##     eta_4   (alp4)   15.000       NA
##     phantm2 (alp2)    1.170       NA
##     phantm3 (alp3)    0.026       NA
##     phantm4 (alp4)   15.000       NA
##     hght01            0.000
##     hght03            0.000
##     hght06            0.000
##     hght09            0.000
##     hght12            0.000
##     hght15            0.000
##     hght18            0.000
##     hght24            0.000
##     hght36            0.000
##
## Variances:
##                    Estimate  Std.Err  Z-value  P(>|z|)
##     eta_1             9.417       NA
##     eta_2             0.000
##     eta_3             0.215       NA
##     eta_4             1.145       NA
##     phantm2           0.000
##     phantm3           0.000
##     phantm4           0.000
##     hght01  (thet)    1.531       NA
##     hght03  (thet)    1.530       NA
##     hght06  (thet)    1.525       NA
##     hght09  (thet)    1.533       NA
##     hght12  (thet)    1.533       NA
##     hght15  (thet)    1.533       NA
##     hght18  (thet)    1.530       NA
##     hght24  (thet)    1.528       NA
##     hght36  (thet)    1.531       NA
##
## Constraints:
##                                                |Slack|
##     L12 - (1-alpha4)                             0.002
##     L22 - (3-alpha4)                             0.000
##     L32 - (6-alpha4)                             0.000
##     L42 - (9-alpha4)                             0.000
##     L52 - (12-alpha4)                            0.000
##     L62 - (15-alpha4)                            0.000
##     L72 - (18-alpha4)                            0.001
##     L82 - (24-alpha4)                            0.000
##     L92 - (36-alpha4)                            0.001
##     L13 - (sqrt((1-alpha4)^2))                   0.000
##     L23 - (sqrt((3-alpha4)^2))                   0.000
##     L33 - (sqrt((6-alpha4)^2))                   0.000
##     L43 - (sqrt((9-alpha4)^2))                   0.000
##     L53 - (sqrt((12-alpha4)^2))                  0.000
##     L63 - (sqrt((15-alpha4)^2))                  0.000
##     L73 - (sqrt((18-alpha4)^2))                  0.000
##     L83 - (sqrt((24-alpha4)^2))                  0.000
##     L93 - (sqrt((36-alpha4)^2))                  0.000
##     L14-((-lph2*((1-4)^2)+3*1-3*4)/((1-4)^2))    0.001
##     L24-((-lph2*((3-4)^2)+3*3-3*4)/((3-4)^2))    0.009
##     L34-((-lph2*((6-4)^2)+3*6-3*4)/((6-4)^2))    0.013
##     L44-((-lph2*((9-4)^2)+3*9-3*4)/((9-4)^2))    0.005
##     L54-((-2*((12-4)^2)+3*12-3*4)/((12-4)^2))    0.003
##     L64-((-2*((15-4)^2)+3*15-3*4)/((15-4)^2))    0.004
##     L74-((-2*((18-4)^2)+3*18-3*4)/((18-4)^2))    0.019
##     L84-((-2*((24-4)^2)+3*24-3*4)/((24-4)^2))    0.008
##     L94-((-2*((36-4)^2)+3*36-3*4)/((36-4)^2))    0.010