Growth Models with Nonlinearity in Time - Multilevel & SEM Implementation in R
Overview
This code walks through fitting of growth models where we have nonlinearity in time.
The code and example provided in this tutorial are from Chapter 10 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 the 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 :166Plot the longitudinal data
ggplot(hght_long, aes(x = age, y = hght, color = as.factor(id), group = id)) +
geom_point() +
geom_line() +
theme_classic(base_size = 18) +
theme(legend.position = "none") +
labs(title = "Individual Height Trajectories", y = "Height (inches)", x = "Age (years)")
Quadratic Growth in MLM using nlme
The first line of code within the nlme() function
resembles formula 10.1, where the outcome variable, hght,
is specified as,
\(y_{ti} = b_{1i} + b_{2i}*((t - k_1)/ k_2) + b_{3i}*((t - k_1)/ k_2)^2 + u_{ti}\)
where \(b_{3i}\) is the quadratic
component. Age is centered at 18 months and scaled to years (divided by
12). na.action = na.omit omits observations with incomplete
values for hght or age.
hght.quad.nlme <- nlme(hght~b_1i+b_2i*((age-18)/12)+b_3i*((age-18)/12)^2,
data=hght_long,
fixed=b_1i+b_2i+b_3i~1,
random=b_1i+b_2i+b_3i~1,
groups=~id,
start=c(30, 10, -3),
na.action=na.exclude)## Warning in nlme.formula(hght ~ b_1i + b_2i * ((age - 18)/12) + b_3i * ((age
## - : Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase
## 'msMaxIter'!## Warning in nlme.formula(hght ~ b_1i + b_2i * ((age - 18)/12) + b_3i * ((age
## - : Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase
## 'msMaxIter'!summary(hght.quad.nlme) ## Nonlinear mixed-effects model fit by maximum likelihood
## Model: hght ~ b_1i + b_2i * ((age - 18)/12) + b_3i * ((age - 18)/12)^2
## Data: hght_long
## AIC BIC logLik
## 2505.736 2549.383 -1242.868
##
## 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.8721557 b_1i b_2i
## b_2i 0.9542989 0.698
## b_3i 0.4289106 -0.425 -0.944
## Residual 1.5717467
##
## Fixed effects: b_1i + b_2i + b_3i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 83.08275 0.3445249 496 241.15162 0
## b_2i 13.80616 0.1380466 496 100.01090 0
## b_3i -3.43110 0.0976734 496 -35.12831 0
## Correlation:
## b_1i b_2i
## b_2i 0.614
## b_3i -0.418 -0.405
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.63073749 -0.62262089 0.06750143 0.64190425 2.59829973
##
## Number of Observations: 581
## Number of Groups: 83Note that there are some warnings that emerge in the fitting process. As well, the correlation between two of the random effects is -0.944. That is very close to -1.00. Given these issues, we may consider a model where the random effect structure has been changed.
Quadratic Growth Model without variation in \(b_{3i}\)
hght.quad.nlme <- nlme(hght~b_1i+b_2i*((age-18)/12)+b_3i*((age-18)/12)^2,
data=hght_long,
fixed=b_1i+b_2i+b_3i~1,
random=b_1i+b_2i~1,
groups=~id,
start=c(30, 10, -3),
na.action=na.exclude)
summary(hght.quad.nlme)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: hght ~ b_1i + b_2i * ((age - 18)/12) + b_3i * ((age - 18)/12)^2
## Data: hght_long
## AIC BIC logLik
## 2515.275 2545.828 -1250.637
##
## Random effects:
## Formula: list(b_1i ~ 1, b_2i ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b_1i 2.7040353 b_1i
## b_2i 0.8885104 0.62
## Residual 1.6303876
##
## Fixed effects: b_1i + b_2i + b_3i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 83.01636 0.3248686 496 255.5383 0
## b_2i 13.78432 0.1344658 496 102.5117 0
## b_3i -3.42943 0.0870154 496 -39.4118 0
## Correlation:
## b_1i b_2i
## b_2i 0.532
## b_3i -0.249 0.010
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.9850629 -0.5640788 0.1010256 0.6286861 2.8576172
##
## Number of Observations: 581
## Number of Groups: 83Quadratic Growth in SEM using lavaan
Data set must be in wide format for SEM implementation
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 <- 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.5Setting up the model.
Quad_model <- '
#intercept factor loadings
eta_1 =~ 1*hght01 +
1*hght03 +
1*hght06 +
1*hght09 +
1*hght12 +
1*hght15 +
1*hght18 +
1*hght24 +
1*hght36
#linear change factor loadings
eta_2 =~ -1.4167*hght01 +
-1.25 *hght03 +
-1 *hght06 +
-0.75 *hght09 +
-0.5 *hght12 +
-0.25 *hght15 +
0 *hght18 +
0.5 *hght24 +
1.5 *hght36
#quadratic change factor loadings
eta_3 =~ 2.1169*hght01 +
1.5625*hght03 +
1 *hght06 +
0.5625*hght09 +
0.25 *hght12 +
0.0625 *hght15 +
0 *hght18 +
0.25 *hght24 +
2.25 *hght36
#manifest variances to be same over time
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
#latent variances
eta_1 ~~ eta_1
eta_2 ~~ eta_2
eta_3 ~~ eta_3
#latent covariances
eta_1 ~~ eta_2
eta_1 ~~ eta_3
eta_2 ~~ eta_3
#latent means
eta_1 ~ start(65)*1
eta_2 ~ start(2)*1
eta_3 ~ start(0.1)*1
#manifest means fixed to 0
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.
quad_fit <- sem(Quad_model, data = hght_wide, meanstructure = T, mimic="mplus")## Warning in lav_object_post_check(object): lavaan WARNING: covariance matrix of latent variables
## is not positive definite;
## use lavInspect(fit, "cov.lv") to investigate.summary(quad_fit, fit.measures = T, standardized = T)## lavaan 0.6-9 ended normally after 60 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 18
## Number of equality constraints 8
##
## Number of observations 83
## Number of missing patterns 31
##
## Model Test User Model:
##
## Test statistic 574.655
## Degrees of freedom 44
## 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.423
## Tucker-Lewis Index (TLI) 0.528
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1223.069
## Loglikelihood unrestricted model (H1) -935.741
##
## Akaike (AIC) 2466.138
## Bayesian (BIC) 2490.326
## Sample-size adjusted Bayesian (BIC) 2458.783
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.381
## 90 Percent confidence interval - lower 0.354
## 90 Percent confidence interval - upper 0.409
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.976
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## eta_1 =~
## hght01 1.000 2.873 1.026
## hght03 1.000 2.873 1.048
## hght06 1.000 2.873 1.045
## hght09 1.000 2.873 1.016
## hght12 1.000 2.873 0.973
## hght15 1.000 2.873 0.927
## hght18 1.000 2.873 0.883
## hght24 1.000 2.873 0.814
## hght36 1.000 2.873 0.761
## eta_2 =~
## hght01 -1.417 -1.355 -0.484
## hght03 -1.250 -1.196 -0.436
## hght06 -1.000 -0.956 -0.348
## hght09 -0.750 -0.717 -0.254
## hght12 -0.500 -0.478 -0.162
## hght15 -0.250 -0.239 -0.077
## hght18 0.000 0.000 0.000
## hght24 0.500 0.478 0.136
## hght36 1.500 1.435 0.380
## eta_3 =~
## hght01 2.117 0.717 0.256
## hght03 1.562 0.529 0.193
## hght06 1.000 0.339 0.123
## hght09 0.562 0.191 0.067
## hght12 0.250 0.085 0.029
## hght15 0.062 0.021 0.007
## hght18 0.000 0.000 0.000
## hght24 0.250 0.085 0.024
## hght36 2.250 0.762 0.202
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## eta_1 ~~
## eta_2 1.967 0.500 3.931 0.000 0.716 0.716
## eta_3 -0.468 0.296 -1.580 0.114 -0.481 -0.481
## eta_2 ~~
## eta_3 -0.411 0.122 -3.378 0.001 -1.268 -1.268
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## eta_1 83.075 0.343 241.978 0.000 28.913 28.913
## eta_2 13.760 0.137 100.344 0.000 14.386 14.386
## eta_3 -3.389 0.089 -38.066 0.000 -10.002 -10.002
## .hght01 0.000 0.000 0.000
## .hght03 0.000 0.000 0.000
## .hght06 0.000 0.000 0.000
## .hght09 0.000 0.000 0.000
## .hght12 0.000 0.000 0.000
## .hght15 0.000 0.000 0.000
## .hght18 0.000 0.000 0.000
## .hght24 0.000 0.000 0.000
## .hght36 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hght01 (thet) 2.325 0.176 13.203 0.000 2.325 0.296
## .hght03 (thet) 2.325 0.176 13.203 0.000 2.325 0.309
## .hght06 (thet) 2.325 0.176 13.203 0.000 2.325 0.307
## .hght09 (thet) 2.325 0.176 13.203 0.000 2.325 0.291
## .hght12 (thet) 2.325 0.176 13.203 0.000 2.325 0.267
## .hght15 (thet) 2.325 0.176 13.203 0.000 2.325 0.242
## .hght18 (thet) 2.325 0.176 13.203 0.000 2.325 0.220
## .hght24 (thet) 2.325 0.176 13.203 0.000 2.325 0.187
## .hght36 (thet) 2.325 0.176 13.203 0.000 2.325 0.163
## eta_1 8.256 1.497 5.517 0.000 1.000 1.000
## eta_2 0.915 0.242 3.783 0.000 1.000 1.000
## eta_3 0.115 0.108 1.063 0.288 1.000 1.000We get a warning Warning message: In lav_object_post_check(object) : lavaan WARNING: covariance matrix of latent variables is not positive definite; use lavInspect(fit, “cov.lv”) to investigate.
As before, we change the structure of the random effects. Adjusting the factor variance-covariance structure
Quad_model_adjusted <- '
#intercept factor loadings
eta_1 =~ 1*hght01 +
1*hght03 +
1*hght06 +
1*hght09 +
1*hght12 +
1*hght15 +
1*hght18 +
1*hght24 +
1*hght36
#linear change factor loadings
eta_2 =~ -1.4167*hght01 +
-1.25 *hght03 +
-1 *hght06 +
-0.75 *hght09 +
-0.5 *hght12 +
-0.25 *hght15 +
0 *hght18 +
0.5 *hght24 +
1.5 *hght36
#quadratic change factor loadings
eta_3 =~ 2.1169*hght01 +
1.5625*hght03 +
1 *hght06 +
0.5625*hght09 +
0.25 *hght12 +
0.0625 *hght15 +
0 *hght18 +
0.25 *hght24 +
2.25 *hght36
#manifest variances to be same over time
hght01 ~~ start(1)*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
#latent variances
eta_1 ~~ eta_1
eta_2 ~~ eta_2
eta_3 ~~ 0*eta_3 #fixing the variance to 0
#latent covariances
eta_1 ~~ eta_2
eta_1 ~~ 0*eta_3 #fixing the covariance to 0
eta_2 ~~ 0*eta_3 #fixing the covariance to 0
#latent means
eta_1 ~ start(65)*1
eta_2 ~ start(2)*1
eta_3 ~ start(0.1)*1
#manifest means fixed to 0
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 adjsted quadratic model.
quad_fit_adjusted <- sem(Quad_model_adjusted, data = hght_wide, meanstructure = T, mimic="mplus")
summary(quad_fit_adjusted, fit.measures = T, standardized = T)## lavaan 0.6-9 ended normally after 45 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 15
## Number of equality constraints 8
##
## Number of observations 83
## Number of missing patterns 31
##
## Model Test User Model:
##
## Test statistic 593.684
## Degrees of freedom 47
## 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.406
## Tucker-Lewis Index (TLI) 0.545
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1232.583
## Loglikelihood unrestricted model (H1) -935.741
##
## Akaike (AIC) 2479.166
## Bayesian (BIC) 2496.098
## Sample-size adjusted Bayesian (BIC) 2474.018
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.374
## 90 Percent confidence interval - lower 0.348
## 90 Percent confidence interval - upper 0.402
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.910
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## eta_1 =~
## hght01 1.000 2.712 1.014
## hght03 1.000 2.712 1.004
## hght06 1.000 2.712 0.984
## hght09 1.000 2.712 0.960
## hght12 1.000 2.712 0.931
## hght15 1.000 2.712 0.899
## hght18 1.000 2.712 0.866
## hght24 1.000 2.712 0.799
## hght36 1.000 2.712 0.674
## eta_2 =~
## hght01 -1.417 -1.287 -0.481
## hght03 -1.250 -1.135 -0.420
## hght06 -1.000 -0.908 -0.330
## hght09 -0.750 -0.681 -0.241
## hght12 -0.500 -0.454 -0.156
## hght15 -0.250 -0.227 -0.075
## hght18 0.000 0.000 0.000
## hght24 0.500 0.454 0.134
## hght36 1.500 1.362 0.338
## eta_3 =~
## hght01 2.117 0.000 0.000
## hght03 1.562 0.000 0.000
## hght06 1.000 0.000 0.000
## hght09 0.562 0.000 0.000
## hght12 0.250 0.000 0.000
## hght15 0.062 0.000 0.000
## hght18 0.000 0.000 0.000
## hght24 0.250 0.000 0.000
## hght36 2.250 0.000 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## eta_1 ~~
## eta_2 1.519 0.431 3.527 0.000 0.617 0.617
## eta_3 0.000 NaN NaN
## eta_2 ~~
## eta_3 0.000 NaN NaN
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## eta_1 83.008 0.323 256.810 0.000 30.606 30.606
## eta_2 13.730 0.135 102.058 0.000 15.117 15.117
## eta_3 -3.395 0.082 -41.497 0.000 -Inf -Inf
## .hght01 0.000 0.000 0.000
## .hght03 0.000 0.000 0.000
## .hght06 0.000 0.000 0.000
## .hght09 0.000 0.000 0.000
## .hght12 0.000 0.000 0.000
## .hght15 0.000 0.000 0.000
## .hght18 0.000 0.000 0.000
## .hght24 0.000 0.000 0.000
## .hght36 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hght01 (thet) 2.447 0.169 14.453 0.000 2.447 0.342
## .hght03 (thet) 2.447 0.169 14.453 0.000 2.447 0.335
## .hght06 (thet) 2.447 0.169 14.453 0.000 2.447 0.322
## .hght09 (thet) 2.447 0.169 14.453 0.000 2.447 0.306
## .hght12 (thet) 2.447 0.169 14.453 0.000 2.447 0.288
## .hght15 (thet) 2.447 0.169 14.453 0.000 2.447 0.269
## .hght18 (thet) 2.447 0.169 14.453 0.000 2.447 0.250
## .hght24 (thet) 2.447 0.169 14.453 0.000 2.447 0.212
## .hght36 (thet) 2.447 0.169 14.453 0.000 2.447 0.151
## eta_1 7.356 1.263 5.825 0.000 1.000 1.000
## eta_2 0.825 0.227 3.627 0.000 1.000 1.000
## eta_3 0.000 NaN NaNThis matches Output 10.5 in the book! Yay!
We can make a diagram to check if we specified and estimated the model as intended.
#diagram of fitted model
semPaths(quad_fit_adjusted,what = "path", whatLabels = "par")
Spline Model in MLM using nlme
Fitting the spline model with a knot point at age 9.
hght.spline.nlme <- nlme(hght~b_1i+b_2i*(pmin(0,age-9))+b_3i*(pmax(0,age-9)),
data=hght_long,
fixed=b_1i+b_2i+b_3i~1,
random=b_1i+b_2i+b_3i~1,
groups=~id,
start=c(60, 10, 6),
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 - 9)) + b_3i * (pmax(0, age - 9))
## Data: hght_long
## AIC BIC logLik
## 2359.183 2402.83 -1169.591
##
## 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.47477913 b_1i b_2i
## b_2i 0.20887748 0.451
## b_3i 0.04662129 0.727 0.862
## Residual 1.34611267
##
## Fixed effects: b_1i + b_2i + b_3i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 73.35059 0.29455052 496 249.02552 0
## b_2i 2.22669 0.03329301 496 66.88164 0
## b_3i 0.88742 0.00949465 496 93.46577 0
## Correlation:
## b_1i b_2i
## b_2i 0.481
## b_3i 0.208 0.104
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.71699930 -0.56033780 0.06312659 0.63762987 2.54231124
##
## Number of Observations: 581
## Number of Groups: 83Here is the same spline model with new coding of the time variables.
We created new variables pre and post which
indicate the knot point to use in the model specification. Not that we
have in this case determined the location of the knot point a
priori.
#prepping new time variables
hght_long$pre = pmin(0,hght_long$age-9)
hght_long$post = pmax(0,hght_long$age-9)
#looking at the new data
head(hght_long)## id age hght pre post
## 1 1 1 53.5 -8 0
## 2 1 3 60.4 -6 0
## 3 1 6 67.2 -3 0
## 4 1 9 70.9 0 0
## 5 1 12 76.2 0 3
## 6 1 15 81.1 0 6#setting up the spline model with the new time variables
hght.spline2.nlme <- nlme(hght~b_1i+b_2i*pre+b_3i*post,
data=hght_long,
fixed=b_1i+b_2i+b_3i~1,
random=b_1i+b_2i+b_3i~1,
groups=~id,
start=c(60, 10, 6),
na.action=na.exclude)
summary(hght.spline2.nlme)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: hght ~ b_1i + b_2i * pre + b_3i * post
## Data: hght_long
## AIC BIC logLik
## 2359.183 2402.83 -1169.591
##
## 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.47477913 b_1i b_2i
## b_2i 0.20887748 0.451
## b_3i 0.04662129 0.727 0.862
## Residual 1.34611267
##
## Fixed effects: b_1i + b_2i + b_3i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 73.35059 0.29455052 496 249.02552 0
## b_2i 2.22669 0.03329301 496 66.88164 0
## b_3i 0.88742 0.00949465 496 93.46577 0
## Correlation:
## b_1i b_2i
## b_2i 0.481
## b_3i 0.208 0.104
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.71699930 -0.56033780 0.06312659 0.63762987 2.54231124
##
## Number of Observations: 581
## Number of Groups: 83Spline Model in SEM using lavaan
Setting up the model with a knot point at age 9.
Spline_model <- '
#intercept factor loadings
eta_1 =~ 1*hght01 +
1*hght03 +
1*hght06 +
1*hght09 +
1*hght12 +
1*hght15 +
1*hght18 +
1*hght24 +
1*hght36
#pre knot linear change factor loadings
eta_2 =~ -8*hght01 +
-6*hght03 +
-3*hght06 +
0*hght09 +
0*hght12 +
0*hght15 +
0*hght18 +
0*hght24 +
0*hght36
#post knot linear change factor loadings
eta_3 =~ 0 *hght01 +
0 *hght03 +
0 *hght06 +
0 *hght09 +
3 *hght12 +
9 *hght18 +
6 *hght15 +
15*hght24 +
27*hght36
#manifest variances to be same over time
hght01 ~~ start(1)*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
#latent variances
eta_1 ~~ start(1)*eta_1
eta_2 ~~ eta_2
eta_3 ~~ eta_3
#latent means
eta_1 ~ start(65)*1
eta_2 ~ start(2)*1
eta_3 ~ start(0.1)*1
#manifest means fixed to 0
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 spline model.
spline_fit <- sem(Spline_model, data = hght_wide, mimic="mplus")
summary(spline_fit)## lavaan 0.6-9 ended normally after 55 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 18
## Number of equality constraints 8
##
## Number of observations 83
## Number of missing patterns 31
##
## Model Test User Model:
##
## Test statistic 467.700
## Degrees of freedom 44
## P-value (Chi-square) 0.000
##
## 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 -8.000
## hght03 -6.000
## hght06 -3.000
## hght09 0.000
## hght12 0.000
## hght15 0.000
## hght18 0.000
## hght24 0.000
## hght36 0.000
## eta_3 =~
## hght01 0.000
## hght03 0.000
## hght06 0.000
## hght09 0.000
## hght12 3.000
## hght18 9.000
## hght15 6.000
## hght24 15.000
## hght36 27.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## eta_1 ~~
## eta_2 0.233 0.096 2.435 0.015
## eta_3 0.084 0.025 3.327 0.001
## eta_2 ~~
## eta_3 0.008 0.003 3.037 0.002
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## eta_1 73.351 0.294 249.542 0.000
## eta_2 2.227 0.033 66.772 0.000
## eta_3 0.887 0.010 93.397 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|)
## .hght01 (thet) 1.812 0.136 13.369 0.000
## .hght03 (thet) 1.812 0.136 13.369 0.000
## .hght06 (thet) 1.812 0.136 13.369 0.000
## .hght09 (thet) 1.812 0.136 13.369 0.000
## .hght12 (thet) 1.812 0.136 13.369 0.000
## .hght15 (thet) 1.812 0.136 13.369 0.000
## .hght18 (thet) 1.812 0.136 13.369 0.000
## .hght24 (thet) 1.812 0.136 13.369 0.000
## .hght36 (thet) 1.812 0.136 13.369 0.000
## eta_1 6.125 1.099 5.572 0.000
## eta_2 0.044 0.014 3.107 0.002
## eta_3 0.002 0.001 2.003 0.045The results are the same as those given in Chapter 10 (p. 228-31). Yay!
We can make a diagram to check if we specified and estimated the model as intended.
#diagram of fitted model
semPaths(spline_fit, what = "path", whatLabels = "par")
Conclusion
We have specified two nonlinear growth models that are nonlinear in time: the quadratic growth model and a spline (or multiphase) growth model with known knot-point. We could fit both models in both the MLM and the SEM frameworks in R. That was very cool!
We are growing too!