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.

Load required libraries and read in the data

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
hght_long <- read.csv(file=url(filepath.long),header=TRUE,na.string='.')
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

#contraints to define factor loadings
  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

#contraints to define factor loadings
  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