This tutorial walks through the fitting of growth models with nonlinearity in parameters 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 11 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.
This tutorial provides line-by-line code to examine growth models with nonlineary in parameters using the:
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.1hght_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.5nlme library.hght.jb.nlme <- nlme(hght~b_1i + b_2i * (age/12) + b_3i * (exp(gamma*(age/12))-1),
data = hght_long,
fixed = b_1i + b_2i + b_3i + gamma ~ 1,
random = b_1i + b_2i + b_3i ~ 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(gamma * (age/12)) - 1)
## Data: hght_long
## AIC BIC logLik
## 2000.239 2048.252 -989.1197
##
## Random effects:
## Formula: list(b_1i ~ 1, b_2i ~ 1, b_3i ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b_1i 2.7199531 b_1i b_2i
## b_2i 0.8590517 0.260
## b_3i 3.0567734 0.583 0.248
## Residual 0.8193595
##
## Fixed effects: b_1i + b_2i + b_3i + gamma ~ 1
## Value Std.Error DF t-value p-value
## b_1i 51.04753 0.3501220 495 145.79927 0
## b_2i 9.31413 0.1649162 495 56.47799 0
## b_3i -17.81609 0.4537913 495 -39.26054 0
## gamma -2.08950 0.0731174 495 -28.57738 0
## Correlation:
## b_1i b_2i b_3i
## b_2i -0.034
## b_3i 0.412 0.543
## gamma 0.350 -0.655 -0.392
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.62957197 -0.51594954 -0.02464767 0.52653169 2.87228859
##
## Number of Observations: 581
## Number of Groups: 83lavaan 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 the "vertical distance between the actual intercept and the intercept of the linear asymptote" (p.252, Grimm, Ram, & Estabrook, 2016)
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
#introducing phantom factor, setting all parameters to 0 and mean to gamma
phantom =~ 0*hght01
phantom ~~ 0*phantom
phantom ~~ 0*eta_1
phantom ~~ 0*eta_2
phantom ~~ 0*eta_3
phantom ~ start(-2)*gamma*1
#contraints to define factor loadings
gamma < -2 #this constraint is added so that it finds the viable solution
L13 == exp(gamma*0.0833)-1
L23 == exp(gamma*0.25)-1
L33 == exp(gamma*0.50)-1
L43 == exp(gamma*0.75)-1
L53 == exp(gamma*1.00)-1
L63 == exp(gamma*1.25)-1
L73 == exp(gamma*1.50)-1
L83 == exp(gamma*2.00)-1
L93 == exp(gamma*3.00)-1
#factor variances
eta_1 ~~ start(7.3)*eta_1
eta_2 ~~ start(0.7)*eta_2
eta_3 ~~ start(9.4)*eta_3
#factor covariances
eta_1 ~~ start(0.61)*eta_2
eta_2 ~~ start(0.67)*eta_3
eta_3 ~~ start(4.85)*eta_1
#latent means
eta_1 ~ start(51)*1
eta_2 ~ start(9.2)*1
eta_3 ~ start(-17.8)*1
#manifest variances
hght01 ~~ start(.67)*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(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.fit1, fit.measures=TRUE)## lavaan (0.5-20) converged normally after 1218 iterations
##
## Number of observations 83
##
## Number of missing patterns 31
##
## Estimator ML
## Minimum Function Test Statistic 106.612
## Degrees of freedom 43
## 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.931
## Tucker-Lewis Index (TLI) 0.942
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -989.047
## Loglikelihood unrestricted model (H1) -935.741
##
## Number of free parameters 11
## Akaike (AIC) 2000.094
## Bayesian (BIC) 2026.702
## Sample-size adjusted Bayesian (BIC) 1992.005
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.134
## 90 Percent Confidence Interval 0.102 0.166
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.309
##
## 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.158 0.005 -30.722 0.000
## hght03 (L23) -0.403 0.011 -36.846 0.000
## hght06 (L33) -0.644 0.013 -49.279 0.000
## hght09 (L43) -0.787 0.012 -67.307 0.000
## hght12 (L53) -0.873 0.009 -93.759 0.000
## hght15 (L63) -0.924 0.007 -132.998 0.000
## hght18 (L73) -0.955 0.005 -191.771 0.000
## hght24 (L83) -0.984 0.002 -415.790 0.000
## hght36 (L93) -0.998 0.000 -2212.563 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_1 ~~
## eta_2 0.615 0.375 1.642 0.101
## eta_2 ~~
## eta_3 0.671 0.489 1.371 0.170
## eta_1 ~~
## eta_3 4.859 1.420 3.422 0.001
##
## Intercepts:
## Estimate Std.Err Z-value P(>|z|)
## phantom (gamm) -2.063 0.073 -28.158 0.000
## eta_1 51.093 0.349 146.216 0.000
## eta_2 9.274 0.167 55.411 0.000
## eta_3 -17.881 0.458 -39.001 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.389 1.352 5.466 0.000
## eta_2 0.740 0.181 4.082 0.000
## eta_3 9.423 2.145 4.393 0.000
## hght01 (thet) 0.671 0.050 13.411 0.000
## hght03 (thet) 0.671 0.050 13.411 0.000
## hght06 (thet) 0.671 0.050 13.411 0.000
## hght09 (thet) 0.671 0.050 13.411 0.000
## hght12 (thet) 0.671 0.050 13.411 0.000
## hght15 (thet) 0.671 0.050 13.411 0.000
## hght18 (thet) 0.671 0.050 13.411 0.000
## hght24 (thet) 0.671 0.050 13.411 0.000
## hght36 (thet) 0.671 0.050 13.411 0.000
##
## Constraints:
## |Slack|
## -2 - (gamma) 0.063
## L13 - (exp(gamma*0.0833)-1) 0.000
## L23 - (exp(gamma*0.25)-1) 0.000
## L33 - (exp(gamma*0.50)-1) 0.000
## L43 - (exp(gamma*0.75)-1) 0.000
## L53 - (exp(gamma*1.00)-1) 0.000
## L63 - (exp(gamma*1.25)-1) 0.000
## L73 - (exp(gamma*1.50)-1) 0.000
## L83 - (exp(gamma*2.00)-1) 0.000
## L93 - (exp(gamma*3.00)-1) 0.000nlme library.for(i in 1:747){
if (hght_long$age[i] == 1) {hght_long$age1[i]=1}
if (hght_long$age[i] != 1) {hght_long$age1[i]=0}
if (hght_long$age[i] == 3) {hght_long$age3[i]=1}
if (hght_long$age[i] != 3) {hght_long$age3[i]=0}
if (hght_long$age[i] == 6) {hght_long$age6[i]=1}
if (hght_long$age[i] != 6) {hght_long$age6[i]=0}
if (hght_long$age[i] == 9) {hght_long$age9[i]=1}
if (hght_long$age[i] != 9) {hght_long$age9[i]=0}
if (hght_long$age[i] == 12) {hght_long$age12[i]=1}
if (hght_long$age[i] != 12) {hght_long$age12[i]=0}
if (hght_long$age[i] == 15) {hght_long$age15[i]=1}
if (hght_long$age[i] != 15) {hght_long$age15[i]=0}
if (hght_long$age[i] == 18) {hght_long$age18[i]=1}
if (hght_long$age[i] != 18) {hght_long$age18[i]=0}
if (hght_long$age[i] == 24) {hght_long$age24[i]=1}
if (hght_long$age[i] != 24) {hght_long$age24[i]=0}
if (hght_long$age[i] == 36) {hght_long$age36[i]=1}
if (hght_long$age[i] != 36) {hght_long$age36[i]=0}
}
hght.latent.nlme <- nlme(hght ~ age1 *(b_1i + b_2i*0) +
age3 *(b_1i + b_2i*A_2) +
age6 *(b_1i + b_2i*A_3) +
age9 *(b_1i + b_2i*A_4) +
age12*(b_1i + b_2i*A_5) +
age15*(b_1i + b_2i*A_6) +
age18*(b_1i + b_2i*A_7) +
age24*(b_1i + b_2i*A_8) +
age36*(b_1i + b_2i*1),
data = hght_long,
fixed = b_1i + b_2i + A_2 + A_3 + A_4 + A_5 + A_6 + A_7 + A_8 ~ 1,
random = b_1i + b_2i ~ 1,
groups =~ id,
start = c(60, 20, .3, .4, .5, .6, .7, .8, .9),
na.action = na.exclude)
summary(hght.latent.nlme)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: hght ~ age1 * (b_1i + b_2i * 0) + age3 * (b_1i + b_2i * A_2) + age6 * (b_1i + b_2i * A_3) + age9 * (b_1i + b_2i * A_4) + age12 * (b_1i + b_2i * A_5) + age15 * (b_1i + b_2i * A_6) + age18 * (b_1i + b_2i * A_7) + age24 * (b_1i + b_2i * A_8) + age36 * (b_1i + b_2i * 1)
## Data: hght_long
## AIC BIC logLik
## 2022.224 2078.966 -998.112
##
## Random effects:
## Formula: list(b_1i ~ 1, b_2i ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b_1i 2.2669367 b_1i
## b_2i 3.1982959 -0.128
## Residual 0.9031948
##
## Fixed effects: b_1i + b_2i + A_2 + A_3 + A_4 + A_5 + A_6 + A_7 + A_8 ~ 1
## Value Std.Error DF t-value p-value
## b_1i 54.54834 0.2797289 490 195.00427 0
## b_2i 42.08689 0.4113588 490 102.31187 0
## A_2 0.14678 0.0035631 490 41.19609 0
## A_3 0.30874 0.0034520 490 89.43782 0
## A_4 0.40956 0.0034124 490 120.02170 0
## A_5 0.50299 0.0033746 490 149.05122 0
## A_6 0.58723 0.0035781 490 164.11636 0
## A_7 0.65704 0.0035792 490 183.57248 0
## A_8 0.78358 0.0036713 490 213.43595 0
## Correlation:
## b_1i b_2i A_2 A_3 A_4 A_5 A_6 A_7
## b_2i -0.239
## A_2 -0.290 0.148
## A_3 -0.244 0.061 0.429
## A_4 -0.215 0.021 0.386 0.398
## A_5 -0.183 -0.005 0.338 0.358 0.365
## A_6 -0.140 -0.047 0.282 0.318 0.338 0.346
## A_7 -0.119 -0.076 0.248 0.293 0.321 0.345 0.341
## A_8 -0.076 -0.122 0.189 0.253 0.291 0.329 0.340 0.365
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.58052647 -0.53575511 0.02044748 0.52029420 4.35506540
##
## Number of Observations: 581
## Number of Groups: 83lavaan library.lb.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 =~ 0*hght01 +
hght03 +
hght06 +
hght09 +
hght12 +
hght15 +
hght18 +
hght24 +
1*hght36
#factor variances
eta_1 ~~ start(60)*eta_1
eta_2 ~~ start(4)*eta_2
#factor covariances
eta_1 ~~ eta_2
#latent means
eta_1 ~ 1
eta_2 ~ 1
#manifest variances (set equal by naming theta)
hght01 ~~ 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.fit2 <- lavaan(lb.hght.lavaan,
data = hght_wide,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE)
summary(lavaan.fit2, fit.measures=TRUE)## lavaan (0.5-20) converged normally after 158 iterations
##
## Number of observations 83
##
## Number of missing patterns 31
##
## Estimator ML
## Minimum Function Test Statistic 124.722
## Degrees of freedom 41
## 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.909
## Tucker-Lewis Index (TLI) 0.920
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -998.102
## Loglikelihood unrestricted model (H1) -935.741
##
## Number of free parameters 13
## Akaike (AIC) 2022.204
## Bayesian (BIC) 2053.649
## Sample-size adjusted Bayesian (BIC) 2012.644
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.157
## 90 Percent Confidence Interval 0.126 0.189
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.054
##
## 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.000
## hght03 0.147 0.004 41.572 0.000
## hght06 0.309 0.003 90.148 0.000
## hght09 0.410 0.003 120.745 0.000
## hght12 0.503 0.003 150.125 0.000
## hght15 0.587 0.004 165.354 0.000
## hght18 0.657 0.004 184.921 0.000
## hght24 0.783 0.004 215.121 0.000
## hght36 1.000
##
## Covariances:
## Estimate Std.Err Z-value P(>|z|)
## eta_1 ~~
## eta_2 -0.928 0.973 -0.954 0.340
##
## Intercepts:
## Estimate Std.Err Z-value P(>|z|)
## eta_1 54.543 0.278 196.354 0.000
## eta_2 42.100 0.409 102.997 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|)
## eta_1 5.139 0.889 5.783 0.000
## eta_2 10.231 1.981 5.165 0.000
## hght01 (thet) 0.816 0.056 14.476 0.000
## hght03 (thet) 0.816 0.056 14.476 0.000
## hght06 (thet) 0.816 0.056 14.476 0.000
## hght09 (thet) 0.816 0.056 14.476 0.000
## hght12 (thet) 0.816 0.056 14.476 0.000
## hght15 (thet) 0.816 0.056 14.476 0.000
## hght18 (thet) 0.816 0.056 14.476 0.000
## hght24 (thet) 0.816 0.056 14.476 0.000
## hght36 (thet) 0.816 0.056 14.476 0.000nlme library.hght.spline.nlme <- nlme(hght ~ b_1i + b_2i*(pmin(0,age-gamma)) + b_3i*(pmax(0,age-gamma)),
data = hght_long,
fixed = b_1i + b_2i + b_3i + gamma ~ 1,
random = b_1i + b_2i + b_3i ~ 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 - gamma)) + b_3i * (pmax(0, age - gamma))
## Data: hght_long
## AIC BIC logLik
## 2287.346 2335.359 -1132.673
##
## Random effects:
## Formula: list(b_1i ~ 1, b_2i ~ 1, b_3i ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b_1i 2.43555767 b_1i b_2i
## b_2i 0.26036627 0.418
## b_3i 0.05674639 0.548 0.550
## Residual 1.21623022
##
## Fixed effects: b_1i + b_2i + b_3i + gamma ~ 1
## Value Std.Error DF t-value p-value
## b_1i 71.82716 0.3353841 495 214.16390 0
## b_2i 2.53762 0.0525520 495 48.28782 0
## b_3i 0.90117 0.0098949 495 91.07447 0
## gamma 7.62649 0.1372391 495 55.57085 0
## Correlation:
## b_1i b_2i b_3i
## b_2i -0.040
## b_3i 0.041 0.199
## gamma 0.514 -0.650 -0.234
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.64758402 -0.55058441 0.01787727 0.58480212 2.59728431
##
## Number of Observations: 581
## Number of Groups: 83lavaan library.est.knot.spline.hght.lavaan <- '
#latent variable definitions
#defining the intercept
eta_1s =~ 1*hght01 +
1*hght03 +
1*hght06 +
1*hght09 +
1*hght12 +
1*hght15 +
1*hght18 +
1*hght24 +
1*hght36
#defining the slope
eta_2s =~ 1*hght01 +
3*hght03 +
6*hght06 +
9*hght09 +
12*hght12 +
15*hght15 +
18*hght18 +
24*hght24 +
36*hght36
#defining eta_3
eta_3s =~ NA*L13*hght01 +
L23*hght03 +
L33*hght06 +
L43*hght09 +
L53*hght12 +
L63*hght15 +
L73*hght18 +
L83*hght24 +
L93*hght36
#latent means
eta_1s ~ start(58)*alpha_1s*1
eta_2s ~ start(1.7)*alpha_2s*1
eta_3s ~ start(-0.8)*alpha_3s*1
#factor variances
eta_1s ~~ start(4.6)*psi_11s*eta_1s
eta_2s ~~ start(0.03)*psi_22s*eta_2s
eta_3s ~~ start(0.02)*psi_33s*eta_3s
#factor covariances
#fixing variances to zero b/c getting negative estimates
eta_1s ~~ start(0)*psi_21s*eta_2s
eta_1s ~~ start(0)*psi_31s*eta_3s
eta_2s ~~ start(-0.01)*psi_32s*eta_3s
#introducing phantom factors, setting all parameters to 0 and mean to gamma
gamma > 7 #added to get convergence
gamma < 8 #added to get convergence
phantom1 =~ 0*hght01
phantom1 ~~ 0*phantom1
phantom1 ~~ 0*eta_1s
phantom1 ~~ 0*eta_2s
phantom1 ~~ 0*eta_3s
phantom1 ~ start(8)*gamma*1
phantom2 =~ 0*hght01
phantom2 ~~ 0*phantom2
phantom2 ~~ 0*phantom1
phantom2 ~~ 0*eta_1s
phantom2 ~~ 0*eta_2s
phantom2 ~~ 0*eta_3s
phantom2 ~ start(71.7)*alpha_1*1
phantom3 =~ 0*hght01
phantom3 ~~ 0*phantom3
phantom3 ~~ 0*phantom2
phantom3 ~~ 0*phantom1
phantom3 ~~ 0*eta_1s
phantom3 ~~ 0*eta_2s
phantom3 ~~ 0*eta_3s
phantom3 ~ start(2.5)*alpha_2*1
phantom4 =~ 0*hght01
phantom4 ~~ 0*phantom4
phantom4 ~~ 0*phantom3
phantom4 ~~ 0*phantom2
phantom4 ~~ 0*phantom1
phantom4 ~~ 0*eta_1s
phantom4 ~~ 0*eta_2s
phantom4 ~~ 0*eta_3s
phantom4 ~ start(0.9)*alpha_3*1
phantom5 =~ 0*hght01
phantom5 ~~ 0*phantom5
phantom5 ~~ 0*phantom4
phantom5 ~~ 0*phantom3
phantom5 ~~ 0*phantom2
phantom5 ~~ 0*phantom1
phantom5 ~~ 0*eta_1s
phantom5 ~~ 0*eta_2s
phantom5 ~~ 0*eta_3s
phantom5 ~ start(5.9)*psi_11*1
phantom6 =~ 0*hght01
phantom6 ~~ 0*phantom6
phantom6 ~~ 0*phantom5
phantom6 ~~ 0*phantom4
phantom6 ~~ 0*phantom3
phantom6 ~~ 0*phantom2
phantom6 ~~ 0*phantom1
phantom6 ~~ 0*eta_1s
phantom6 ~~ 0*eta_2s
phantom6 ~~ 0*eta_3s
phantom6 ~ start(0.05)*psi_22*1
phantom7 =~ 0*hght01
phantom7 ~~ 0*phantom7
phantom7 ~~ 0*phantom6
phantom7 ~~ 0*phantom5
phantom7 ~~ 0*phantom4
phantom7 ~~ 0*phantom3
phantom7 ~~ 0*phantom2
phantom7 ~~ 0*phantom1
phantom7 ~~ 0*eta_1s
phantom7 ~~ 0*eta_2s
phantom7 ~~ 0*eta_3s
phantom7 ~ start(0.001)*psi_33*1
phantom8 =~ 0*hght01
phantom8 ~~ 0*phantom8
phantom8 ~~ 0*phantom7
phantom8 ~~ 0*phantom6
phantom8 ~~ 0*phantom5
phantom8 ~~ 0*phantom4
phantom8 ~~ 0*phantom3
phantom8 ~~ 0*phantom2
phantom8 ~~ 0*phantom1
phantom8 ~~ 0*eta_1s
phantom8 ~~ 0*eta_2s
phantom8 ~~ 0*eta_3s
phantom8 ~ start(0.26)*psi_21*1
phantom9 =~ 0*hght01
phantom9 ~~ 0*phantom9
phantom9 ~~ 0*phantom8
phantom9 ~~ 0*phantom7
phantom9 ~~ 0*phantom6
phantom9 ~~ 0*phantom5
phantom9 ~~ 0*phantom4
phantom9 ~~ 0*phantom3
phantom9 ~~ 0*phantom2
phantom9 ~~ 0*phantom1
phantom9 ~~ 0*eta_1s
phantom9 ~~ 0*eta_2s
phantom9 ~~ 0*eta_3s
phantom9 ~ start(0.07)*psi_31*1
phantom10 =~ 0*hght01
phantom10 ~~ 0*phantom10
phantom10 ~~ 0*phantom9
phantom10 ~~ 0*phantom8
phantom10 ~~ 0*phantom7
phantom10 ~~ 0*phantom6
phantom10 ~~ 0*phantom5
phantom10 ~~ 0*phantom4
phantom10 ~~ 0*phantom3
phantom10 ~~ 0*phantom2
phantom10 ~~ 0*phantom1
phantom10 ~~ 0*eta_1s
phantom10 ~~ 0*eta_2s
phantom10 ~~ 0*eta_3s
phantom10 ~ start(0.008)*psi_32*1
#contraints to define factor loadings
L13 == sqrt(( 1 - gamma)^2)
L23 == sqrt(( 3 - gamma)^2)
L33 == sqrt(( 6 - gamma)^2)
L43 == sqrt(( 9 - gamma)^2)
L53 == sqrt((12 - gamma)^2)
L63 == sqrt((15 - gamma)^2)
L73 == sqrt((18 - gamma)^2)
L83 == sqrt((24 - gamma)^2)
L93 == sqrt((36 - gamma)^2)
alpha_1 == alpha_1s + gamma*alpha_2s
alpha_2 == alpha_2s - alpha_3s
alpha_3 == alpha_2s + alpha_3s
psi_11 == (psi_11s + gamma*gamma*psi_22s) + 2*psi_21s*gamma
psi_22 == (psi_22s + psi_33s) - 2*psi_32s
psi_33 == (psi_22s + psi_33s) + 2*psi_32s
psi_21 == psi_21s - psi_31s + gamma*(psi_22s-psi_32s)
psi_31 == psi_21s + psi_31s + gamma*(psi_22s+psi_32s)
psi_32 == psi_22s - psi_33s
#manifest variances
hght01 ~~ start(1.47)*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(est.knot.spline.hght.lavaan,
data = hght_wide,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE,
mimic="mplus",
control=list(iter.max=1000),
verbose=FALSE)
summary(lavaan.fit1, fit.measures=TRUE)## lavaan (0.5-20) converged normally after 830 iterations
##
## Number of observations 83
##
## Number of missing patterns 31
##
## Estimator ML
## Minimum Function Test Statistic 393.205
## Degrees of freedom 43
## 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.619
## Tucker-Lewis Index (TLI) 0.681
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1132.344
## Loglikelihood unrestricted model (H1) -935.741
##
## Number of free parameters 11
## Akaike (AIC) 2286.687
## Bayesian (BIC) 2313.294
## Sample-size adjusted Bayesian (BIC) 2278.598
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.313
## 90 Percent Confidence Interval 0.285 0.342
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.494
##
## Parameter Estimates:
##
## Information Observed
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err Z-value P(>|z|)
## eta_1s =~
## 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_2s =~
## hght01 1.000
## hght03 3.000
## hght06 6.000
## hght09 9.000
## hght12 12.000
## hght15 15.000
## hght18 18.000
## hght24 24.000
## hght36 36.000
## eta_3s =~
## hght01 (L13) 6.517 0.133 48.979 0.000
## hght03 (L23) 4.517 0.133 33.947 0.000
## hght06 (L33) 1.517 0.133 11.400 0.000
## hght09 (L43) 1.483 0.133 11.148 0.000
## hght12 (L53) 4.483 0.133 33.696 0.000
## hght15 (L63) 7.483 0.133 56.244 0.000
## hght18 (L73) 10.483 0.133 78.792 0.000
## hght24 (L83) 16.483 0.133 123.887 0.000
## hght36 (L93) 28.483 0.133 214.079 0.000
## phantom1 =~
## hght01 0.000
## phantom2 =~
## hght01 0.000
## phantom3 =~
## hght01 0.000
## phantom4 =~
## hght01 0.000
## phantom5 =~
## hght01 0.000
## phantom6 =~
## hght01 0.000
## phantom7 =~
## hght01 0.000
## phantom8 =~
## hght01 0.000
## phantom9 =~
## hght01 0.000
## phantom10 =~
## hght01 0.000
##
## Covariances:
## Estimate Std.Err Z-value P(>|z|)
## eta_1s ~~
## eta_2s (p_21) 0.003 0.048 0.057 0.955
## eta_3s (p_31) 0.030 0.046 0.642 0.521
## eta_2s ~~
## eta_3s (p_32) -0.017 0.005 -3.244 0.001
## eta_1s ~~
## phantm1 0.000
## eta_2s ~~
## phantm1 0.000
## eta_3s ~~
## phantm1 0.000
## phantom1 ~~
## phantm2 0.000
## eta_1s ~~
## phantm2 0.000
## eta_2s ~~
## phantm2 0.000
## eta_3s ~~
## phantm2 0.000
## phantom2 ~~
## phantm3 0.000
## phantom1 ~~
## phantm3 0.000
## eta_1s ~~
## phantm3 0.000
## eta_2s ~~
## phantm3 0.000
## eta_3s ~~
## phantm3 0.000
## phantom3 ~~
## phantm4 0.000
## phantom2 ~~
## phantm4 0.000
## phantom1 ~~
## phantm4 0.000
## eta_1s ~~
## phantm4 0.000
## eta_2s ~~
## phantm4 0.000
## eta_3s ~~
## phantm4 0.000
## phantom4 ~~
## phantm5 0.000
## phantom3 ~~
## phantm5 0.000
## phantom2 ~~
## phantm5 0.000
## phantom1 ~~
## phantm5 0.000
## eta_1s ~~
## phantm5 0.000
## eta_2s ~~
## phantm5 0.000
## eta_3s ~~
## phantm5 0.000
## phantom5 ~~
## phantm6 0.000
## phantom4 ~~
## phantm6 0.000
## phantom3 ~~
## phantm6 0.000
## phantom2 ~~
## phantm6 0.000
## phantom1 ~~
## phantm6 0.000
## eta_1s ~~
## phantm6 0.000
## eta_2s ~~
## phantm6 0.000
## eta_3s ~~
## phantm6 0.000
## phantom6 ~~
## phantm7 0.000
## phantom5 ~~
## phantm7 0.000
## phantom4 ~~
## phantm7 0.000
## phantom3 ~~
## phantm7 0.000
## phantom2 ~~
## phantm7 0.000
## phantom1 ~~
## phantm7 0.000
## eta_1s ~~
## phantm7 0.000
## eta_2s ~~
## phantm7 0.000
## eta_3s ~~
## phantm7 0.000
## phantom7 ~~
## phantm8 0.000
## phantom6 ~~
## phantm8 0.000
## phantom5 ~~
## phantm8 0.000
## phantom4 ~~
## phantm8 0.000
## phantom3 ~~
## phantm8 0.000
## phantom2 ~~
## phantm8 0.000
## phantom1 ~~
## phantm8 0.000
## eta_1s ~~
## phantm8 0.000
## eta_2s ~~
## phantm8 0.000
## eta_3s ~~
## phantm8 0.000
## phantom8 ~~
## phantm9 0.000
## phantom7 ~~
## phantm9 0.000
## phantom6 ~~
## phantm9 0.000
## phantom5 ~~
## phantm9 0.000
## phantom4 ~~
## phantm9 0.000
## phantom3 ~~
## phantm9 0.000
## phantom2 ~~
## phantm9 0.000
## phantom1 ~~
## phantm9 0.000
## eta_1s ~~
## phantm9 0.000
## eta_2s ~~
## phantm9 0.000
## eta_3s ~~
## phantm9 0.000
## phantom9 ~~
## phntm10 0.000
## phantom8 ~~
## phntm10 0.000
## phantom7 ~~
## phntm10 0.000
## phantom6 ~~
## phntm10 0.000
## phantom5 ~~
## phntm10 0.000
## phantom4 ~~
## phntm10 0.000
## phantom3 ~~
## phntm10 0.000
## phantom2 ~~
## phntm10 0.000
## phantom1 ~~
## phntm10 0.000
## eta_1s ~~
## phntm10 0.000
## eta_2s ~~
## phntm10 0.000
## eta_3s ~~
## phntm10 0.000
##
## Intercepts:
## Estimate Std.Err Z-value P(>|z|)
## et_1 (alph_1s) 58.653 0.264 222.322 0.000
## et_2 (alph_2s) 1.734 0.028 62.401 0.000
## et_3 (alph_3s) -0.831 0.026 -32.163 0.000
## phn1 (gamm) 7.517 0.133 56.495 0.000
## phn2 (alpha_1) 71.688 0.333 215.282 0.000
## phn3 (alpha_2) 2.565 0.053 48.639 0.000
## phn4 (alpha_3) 0.903 0.010 90.984 0.000
## phn5 (p_11) 5.906 1.045 5.650 0.000
## phn6 (p_22) 0.070 0.020 3.403 0.001
## phn7 (p_33) 0.003 0.001 2.915 0.004
## phn8 (p_21) 0.266 0.112 2.365 0.018
## phn9 (p_31) 0.076 0.025 3.037 0.002
## ph10 (p_32) 0.008 0.003 2.433 0.015
## hg01 0.000
## hg03 0.000
## hg06 0.000
## hg09 0.000
## hg12 0.000
## hg15 0.000
## hg18 0.000
## hg24 0.000
## hg36 0.000
##
## Variances:
## Estimate Std.Err Z-value P(>|z|)
## eta_1s (p_11) 4.603 0.820 5.615 0.000
## eta_2s (p_22) 0.022 0.006 4.056 0.000
## eta_3s (p_33) 0.014 0.005 2.667 0.008
## phantm1 0.000
## phantm2 0.000
## phantm3 0.000
## phantm4 0.000
## phantm5 0.000
## phantm6 0.000
## phantm7 0.000
## phantm8 0.000
## phantm9 0.000
## phntm10 0.000
## hght01 (thet) 1.476 0.110 13.379 0.000
## hght03 (thet) 1.476 0.110 13.379 0.000
## hght06 (thet) 1.476 0.110 13.379 0.000
## hght09 (thet) 1.476 0.110 13.379 0.000
## hght12 (thet) 1.476 0.110 13.379 0.000
## hght15 (thet) 1.476 0.110 13.379 0.000
## hght18 (thet) 1.476 0.110 13.379 0.000
## hght24 (thet) 1.476 0.110 13.379 0.000
## hght36 (thet) 1.476 0.110 13.379 0.000
##
## Constraints:
## |Slack|
## gamma - (7) 0.517
## 8 - (gamma) 0.483
## L13 - (sqrt((1-gamma)^2)) 0.000
## L23 - (sqrt((3-gamma)^2)) 0.000
## L33 - (sqrt((6-gamma)^2)) 0.000
## L43 - (sqrt((9-gamma)^2)) 0.000
## L53 - (sqrt((12-gamma)^2)) 0.000
## L63 - (sqrt((15-gamma)^2)) 0.000
## L73 - (sqrt((18-gamma)^2)) 0.000
## L83 - (sqrt((24-gamma)^2)) 0.000
## L93 - (sqrt((36-gamma)^2)) 0.000
## alpha_1 - (alpha_1s+gamma*alpha_2s) 0.000
## alpha_2 - (alpha_2s-alpha_3s) 0.000
## alpha_3 - (alpha_2s+alpha_3s) 0.000
## ps_11-((ps_11s+gmm*gmm*ps_22s)+2*ps_21s*) 0.000
## psi_22 - ((psi_22s+psi_33s)-2*psi_32s) 0.000
## psi_33 - ((psi_22s+psi_33s)+2*psi_32s) 0.000
## ps_21-(ps_21s-ps_31s+gmm*(ps_22s-ps_32s)) 0.000
## ps_31-(ps_21s+ps_31s+gmm*(ps_22s+ps_32s)) 0.000
## psi_32 - (psi_22s-psi_33s) 0.000