Growth Modeling Chapter 10: Growth Models with Nonlinearity in Time
Overview
This tutorial walks through a few helpful initial steps before conducting nonlinear growth curve analyses (or any analyses for that matter).
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; however, the chpater should be referred to for further interpretations and insights about the analyses.
Outline
This tutorial provides line-by-line code to
- plot raw data
- fit a multilevel quadratic growth model with
nlme - fit a multilevel spline model with
nlme - fit a quadratic growth model with
lavaan - fit a spline model with
lavaan
Preliminaries: Read in the data and call needed libraries.
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/bgs_height_long.csv"
hght_long <- read.csv(file=url(filepath),header=TRUE)
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.1Quick descriptives of our dataset
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 :166Loading necessary libraries
library(nlme)
library(lavaan)
library(ggplot2)Plot 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 (cm)", x = "Age (months)")
Multilevel Modeling: Fit a Quadratic Growth in 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.exclude omits observations with incomplete values for age using FIML.
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)
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.8722010 b_1i b_2i
## b_2i 0.9542761 0.698
## b_3i 0.4289492 -0.425 -0.944
## Residual 1.5717470
##
## Fixed effects: b_1i + b_2i + b_3i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 83.08276 0.3445299 496 241.14815 0
## b_2i 13.80616 0.1380442 496 100.01263 0
## b_3i -3.43111 0.0976758 496 -35.12752 0
## Correlation:
## b_1i b_2i
## b_2i 0.614
## b_3i -0.419 -0.405
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.63071485 -0.62260082 0.06748581 0.64194540 2.59829087
##
## Number of Observations: 581
## Number of Groups: 83Quadratic 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.8885105 0.62
## Residual 1.6303875
##
## 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.9850630 -0.5640788 0.1010256 0.6286861 2.8576172
##
## Number of Observations: 581
## Number of Groups: 83The results are the same as those given in Chapter 10 (p. 213).
Multilevel Modeling: Spline Model in nlme
The first line of code within the nlme() function resembles formula 10.3, where the outcome variable, hght, is specified as,
\(y_{ti} = b_{1i} + b_{2i}*(min(t - k_1,0)) + b_{3i}*(max(t - k_1,0)) + u_{ti}\)
where \(k_{1}\) is the knot point, \(b_{2i}\) is the random pre knot point slope for an individual, and \(b_{3i}\) is the random post knot point slope for an individual. Here, the knot point is set to 9 months, so we subtract age by 9. na.action = na.exclude omits observations with incomplete values for age using FIML.
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.4747762 b_1i b_2i
## b_2i 0.2088762 0.451
## b_3i 0.0466204 0.727 0.862
## Residual 1.3461151
##
## Fixed effects: b_1i + b_2i + b_3i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 73.35059 0.29455024 496 249.02575 0
## b_2i 2.22669 0.03329293 496 66.88180 0
## b_3i 0.88742 0.00949459 496 93.46634 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.71701453 -0.56032355 0.06312546 0.63763290 2.54232230
##
## Number of Observations: 581
## Number of Groups: 83The results are the same as those given in Chapter 10 (p. 217).
Here’s the same spline model. We created new variables pre and post which indicate the knot point to use in the model specification.
hght_long$pre = pmin(0,hght_long$age-9)
hght_long$post = pmax(0,hght_long$age-9)
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.4747762 b_1i b_2i
## b_2i 0.2088762 0.451
## b_3i 0.0466204 0.727 0.862
## Residual 1.3461151
##
## Fixed effects: b_1i + b_2i + b_3i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 73.35059 0.29455024 496 249.02575 0
## b_2i 2.22669 0.03329293 496 66.88180 0
## b_3i 0.88742 0.00949459 496 93.46634 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.71701453 -0.56032355 0.06312546 0.63763290 2.54232230
##
## Number of Observations: 581
## Number of Groups: 83These results are identical to the previous model.
SEM: Quadratic Growth in lavaan
Data set must be in wide format for SEM implementation, so we will use the reshape() function.
hght_long <- hght_long[,c("id", "age", "hght")] # excluding the pre and post variables we had created
hght_wide <- reshape(hght_long, #data set
v.names='hght', #repeated measures variable
idvar='id', #id variable
timevar='age', #time metric/occasion variable
direction='wide') #direction of re-structuring
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
## 10 2 54.5 58.8 66.0 69.1 74.5 78.2 81.2 86.7 NA
## 19 3 49.5 56.0 61.9 64.1 67.9 69.6 73.8 NA 85.1
## 28 4 56.9 64.6 70.0 72.5 76.1 80.4 82.3 87.9 94.3
## 37 5 56.7 63.1 NA 73.1 77.3 81.8 83.5 88.9 99.1
## 46 6 56.1 61.9 67.9 72.3 77.3 79.9 82.3 88.3 97.5Now, each column aside from the id column corresponds to one age’s height values.
Age is centered at 18 months and scaled by year. We also provide starting values for the three latent factors to aid in model estimation.
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*65
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
'quad_fit <- sem(Quad_model, data = hght_wide, meanstructure = T, mimic="mplus")## Warning in lav_object_post_check(lavobject): lavaan WARNING: covariance matrix of latent variables
## is not positive definite;
## use inspect(fit,"cov.lv") to investigate.summary(quad_fit, fit.measures = T, standardized = T)## lavaan (0.5-22) converged normally after 106 iterations
##
## Number of observations 83
##
## Number of missing patterns 31
##
## Estimator ML
## Minimum Function Test Statistic 574.655
## Degrees of freedom 44
## 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.423
## Tucker-Lewis Index (TLI) 0.528
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1223.069
## Loglikelihood unrestricted model (H1) -935.741
##
## Number of free parameters 10
## 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 0.354 0.409
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.976
##
## Parameter Estimates:
##
## Information Observed
## Standard Errors Standard
##
## 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
## (Eieee:cnff") 83.075 0.343 241.977 0.000 28.913 28.913
## 13.760 0.137 100.344 0.000 14.386 14.386
## -3.389 0.089 -38.066 0.000 -10.002 -10.002
## . 0.000 0.000 0.000
## . 0.000 0.000 0.000
## . 0.000 0.000 0.000
## . 0.000 0.000 0.000
## . 0.000 0.000 0.000
## . 0.000 0.000 0.000
## . 0.000 0.000 0.000
## . 0.000 0.000 0.000
## . 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 from lavaan that the covariance of latent variables is not positive definite, so we will constrain the variance of \(eta_3\) similiar to the variance constraint implemented in the quadratic multilevel model.
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
'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.5-22) converged normally after 45 iterations
##
## Number of observations 83
##
## Number of missing patterns 31
##
## Estimator ML
## Minimum Function Test Statistic 593.684
## Degrees of freedom 47
## 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.406
## Tucker-Lewis Index (TLI) 0.545
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1232.583
## Loglikelihood unrestricted model (H1) -935.741
##
## Number of free parameters 7
## 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 0.348 0.402
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.910
##
## Parameter Estimates:
##
## Information Observed
## Standard Errors Standard
##
## 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 NaNThe results are the same as those given in Chapter 10 (p. 225-7).
SEM: Spline Model fit in lavaan
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
'Running the spline model
spline_fit <- sem(Spline_model, data = hght_wide, mimic="mplus")
summary(spline_fit)## lavaan (0.5-22) converged normally after 55 iterations
##
## Number of observations 83
##
## Number of missing patterns 31
##
## Estimator ML
## Minimum Function Test Statistic 467.700
## Degrees of freedom 44
## P-value (Chi-square) 0.000
##
## 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 -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).
Overall, this tutorial provides two nonlinear methods (quadratic growth models and spline models) and two frameworks (multilevel modeling and structural equation modeling) with which to characterize longitudinal measures.