Growth Models with Nonlinearity in Parameters - Multilevel & SEM Implementation in R
Overview
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.
Preliminaries
Loading libraries used in this script.
library(tidyverse) #for data management
library(psych) #for basic functions
library(ggplot2) #for plotting
library(nlme) #for mixed effects models
library(lavaan) #for SEM
library(semPlot) #for making SEM diagramsRead in and prep the long and wide data.
filepath <- "https://raw.githubusercontent.com/LRI-2/Data/main/GrowthModeling/bgs_height_long.dat"
#read in the text data file using the url() function
hght_long <- read.table(file=url(filepath),
na.strings = ".") #indicates the missing data designator
#Add names the columns of the data set
names(hght_long) <- c('id', 'age', 'hght')
#view the first few observations in the data set
head(hght_long, 10)## 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
## 7 1 18 NA
## 8 1 24 89.5
## 9 1 36 95.9
## 10 2 1 54.5#quick summary of the data
summary(hght_long)## id age hght
## Min. : 1 Min. : 1.00 Min. : 49.50
## 1st Qu.:21 1st Qu.: 6.00 1st Qu.: 65.00
## Median :42 Median :12.00 Median : 74.20
## Mean :42 Mean :13.78 Mean : 74.34
## 3rd Qu.:63 3rd Qu.:18.00 3rd Qu.: 83.20
## Max. :83 Max. :36.00 Max. :106.00
## NA's :166#creating a wide version of the data file for the SEM framework
hght_wide <- hght_long %>%
spread(age, hght)
names(hght_wide) <- c("id","hght01","hght03", "hght06", "hght09","hght12", "hght15", "hght18", "hght24", "hght36")
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.5Jenss-Bayley Growth Model in MLM using nlme
Fitting the model using the nlme 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.0567735 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.0731175 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.62957185 -0.51594951 -0.02464761 0.52653168 2.87228852
##
## Number of Observations: 581
## Number of Groups: 83Jenss-Bayley Growth Model in SEM using lavaan
Specifying the model.
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, 2006)
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 ~~ start(.67)*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'Fitting the model.
hght.jb.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(hght.jb.fit, fit.measures=TRUE)## lavaan 0.6-9 ended normally after 1255 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 28
## Number of inequality constraints 1
##
## Number of observations 83
## Number of missing patterns 31
##
## Model Test User Model:
##
## Test statistic 106.612
## Degrees of freedom 43
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## 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
##
## Akaike (AIC) 2000.094
## Bayesian (BIC) 2026.701
## Sample-size adjusted Bayesian (BIC) 1992.005
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.134
## 90 Percent confidence interval - lower 0.102
## 90 Percent confidence interval - upper 0.166
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.309
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## 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.760 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.792 0.000
## hght36 (L93) -0.998 0.000 -2212.578 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.217 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.000We can try to make a diagram to check if we specified and estimated the model as intended.
#diagram of fitted model
semPaths(hght.jb.fit, what = "path", whatLabels = "par")However, it does not work, probably becasue the function does not know how to deal with the phantom variables and constraints.
Latent Basis Growth Model in MLM using nlme
Fitting the model using the nlme library.
#creating a set of dummy variables that can be used to "turn on" and "turn off" each possible age value in each row of data
for(i in 1:length(hght_long$age)){
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}
}
#fitting the model in nlme
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.58052646 -0.53575511 0.02044749 0.52029421 4.35506534
##
## Number of Observations: 581
## Number of Groups: 83Latent Basis Growth Model in SEM using lavaan
Specifying the model
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'Fitting the model
lb.hght.fit <- lavaan(lb.hght.lavaan,
data = hght_wide,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE)
summary(lb.hght.fit, fit.measures=TRUE)## lavaan 0.6-9 ended normally after 143 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 21
## Number of equality constraints 8
##
## Number of observations 83
## Number of missing patterns 31
##
## Model Test User Model:
##
## Test statistic 124.722
## Degrees of freedom 41
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## 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
##
## 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 - lower 0.126
## 90 Percent confidence interval - upper 0.189
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.054
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## 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.000We can try to make a diagram to check if we specified and estimated the model as intended.
#diagram of fitted model
semPaths(lb.hght.fit, what = "path", whatLabels = "par")
Bilinear Spline Growth Model with estimated knot point in MLM using nlme
Fitting the model using the nlme 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.exclude)
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.43555766 b_1i b_2i
## b_2i 0.26036626 0.418
## b_3i 0.05674639 0.548 0.550
## Residual 1.21623023
##
## 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.64758397 -0.55058442 0.01787727 0.58480211 2.59728429
##
## Number of Observations: 581
## Number of Groups: 83Bilinear Spline Growth Model with estimated knot point in SEM using lavaan
Note: this code contains a few bugs.
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'Fitting the model.
spline.hght.fit <- 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(spline.hght.fit, fit.measures=TRUE)## lavaan 0.6-9 ended normally after 845 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 37
## Number of inequality constraints 2
## Row rank of the constraints matrix 27
##
## Number of observations 83
## Number of missing patterns 31
##
## Model Test User Model:
##
## Test statistic 393.205
## Degrees of freedom 43
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## 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
##
## 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 - lower 0.285
## 90 Percent confidence interval - upper 0.342
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.494
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## 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.243 0.000
## hght18 (L73) 10.483 0.133 78.791 0.000
## hght24 (L83) 16.483 0.133 123.886 0.000
## hght36 (L93) 28.483 0.133 214.077 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.323 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.434 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.000We can try to make a diagram to check if we specified and estimated the model as intended.
#diagram of fitted model
semPaths(spline.hght.fit, what = "path", whatLabels = "par")However, again this does not work, probably because the function does not know how to work with the phantom variables.
Conclusion
We have demonstrated how to fit growth models that have nonlinearity in their parameters. This greatly expands what it possible. However, be cautious in interpretaion regarding the homoggeneity of shape or knot point across persons.
All that change is so cool!