Growth Models with Nonlinearity in Random Coefficients - Multilevel & SEM Implementation in R
Overview
This tutorial walks through the fitting of growth models with nonlinearity in random coefficients in several different frameworks (e.g., multilevel modeling framework, structural equation modeling framework), and demonstrates these models using different R packages (knowing how to fit the models in different packages can be helpful when trying to fit more complex models as each packages as its own advantages and limitations). Our examples make use of data that are repeated measures of individuals’ height during infancy and early childhood.
The code and example provided in this tutorial are from Chapter 12 of Grimm, Ram, and Estabrook (2016), with a few additions in code and commentary. Please refer to the chapter for further interpretations and insights about the analyses.
Practical Prelimininaries
Load required libraries.
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
We read in the long data file
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 :166Data set must be in wide format for SEM implementations
Converting to wide data
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.5Plot the longitudinal data
Plotting the observed individual trajectories using the long 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)")
Jenss-Bayley Growth Model
Direct Optimization in MLM framework
Fitting the model using the nlme library.
hght.jb.nlme <- nlme(hght~b_1i + b_2i * (age/12) + b_3i * (exp(b_4i*(age/12))-1),
data = hght_long,
fixed = b_1i + b_2i + b_3i + b_4i ~ 1,
random = b_1i + b_2i + b_3i + b_4i ~ 1,
groups = ~ id,
start = c(50, 10, -18, -2),
#control=nlmeControl(msMaxIter=100, msVerbose =TRUE, opt="nlminb"),
na.action = na.exclude)
summary(hght.jb.nlme)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: hght ~ b_1i + b_2i * (age/12) + b_3i * (exp(b_4i * (age/12)) - 1)
## Data: hght_long
## AIC BIC logLik
## 2005.178 2070.649 -987.589
##
## Random effects:
## Formula: list(b_1i ~ 1, b_2i ~ 1, b_3i ~ 1, b_4i ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b_1i 2.6774717 b_1i b_2i b_3i
## b_2i 0.8717594 0.319
## b_3i 3.4099982 0.491 0.449
## b_4i 0.2746997 -0.120 -0.275 -0.642
## Residual 0.8041709
##
## Fixed effects: b_1i + b_2i + b_3i + b_4i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 50.99063 0.3422356 495 148.99274 0
## b_2i 9.31946 0.1607677 495 57.96847 0
## b_3i -17.84870 0.4814404 495 -37.07355 0
## b_4i -2.13255 0.0769104 495 -27.72776 0
## Correlation:
## b_1i b_2i b_3i
## b_2i 0.022
## b_3i 0.376 0.607
## b_4i 0.263 -0.642 -0.536
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.482497934 -0.518215058 0.001139802 0.517855522 2.745168922
##
## Number of Observations: 581
## Number of Groups: 83Notice that there was a warning. It might be useful to try a few different options in the (currently commented out) nlmeControl() function in order to see how stable the resutls are with different optimizer and iterations. The msVerbose = TRUE option is useful for tracking how the optimization is progressing.
Plotting the predicted trajectories. Note that this is not totally accurate - as it does not quite get the trajectories smooth enough. Rather than just using the model object predicted scores, we should generate fine-grained data based on the parameter estimates so that we can see the full continuous curvature of each individual’s predicted trajectory more precisely.
#obtaining predicted scores for individuals
hght_long$pred_jb <- predict(hght.jb.nlme, newdata=hght_long)
#obtaining predicted scores for prototype
hght_long$proto_jb <- predict(hght.jb.nlme, newdata=hght_long, level=0)
#plotting predicted trajectories
#intraindividual change trajetories
ggplot(data = hght_long, aes(x = age, y = pred_jb, group = id)) +
#geom_point(color="black") +
geom_line(color="black") +
geom_line(aes(x = age, y = proto_jb), color="red",size=2) +
theme_classic(base_size = 18) +
theme(legend.position = "none") +
labs(title = "Predicted Individual Height Trajectories", y = "Height (cm)", x = "Age (months)")
First Order Approximation in the MLM framework
To ease the computational difficulty of the estimation, it may be useful to estimate the model using a first-order approximation. Details on how the model is rewritten as a first-order approximation are given in the chapter.
Fitting the first-order model expansion using the nlme
library.
hght.jb.nlme.fo <- nlme(hght ~ (b_1 + d_1i) + (b_2 + d_2i) * (age/12) +
(b_3 + d_3i) * (exp(b_4 * (age/12))-1) +
d_4i * (b_3 * (age/12) * exp(b_4 * (age/12))),
data = hght_long,
fixed = b_1 + b_2 + b_3 + b_4 ~ 1,
random = d_1i + d_2i + d_3i + d_4i ~ 1,
groups =~ id,
start = c(50, 10, -18, -2),
na.action = na.exclude)
summary(hght.jb.nlme.fo)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: hght ~ (b_1 + d_1i) + (b_2 + d_2i) * (age/12) + (b_3 + d_3i) * (exp(b_4 * (age/12)) - 1) + d_4i * (b_3 * (age/12) * exp(b_4 * (age/12)))
## Data: hght_long
## AIC BIC logLik
## 2005.212 2070.684 -987.6062
##
## Random effects:
## Formula: list(d_1i ~ 1, d_2i ~ 1, d_3i ~ 1, d_4i ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## d_1i 2.7167333 d_1i d_2i d_3i
## d_2i 0.8928622 0.267
## d_3i 3.3840056 0.482 0.433
## d_4i 0.2775733 0.027 -0.299 -0.525
## Residual 0.8027898
##
## Fixed effects: b_1 + b_2 + b_3 + b_4 ~ 1
## Value Std.Error DF t-value p-value
## b_1 50.96847 0.3478285 495 146.53333 0
## b_2 9.33386 0.1602783 495 58.23531 0
## b_3 -17.83237 0.4667396 495 -38.20624 0
## b_4 -2.11292 0.0752808 495 -28.06722 0
## Correlation:
## b_1 b_2 b_3
## b_2 -0.006
## b_3 0.385 0.578
## b_4 0.326 -0.634 -0.462
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.43002233 -0.51726774 -0.02358415 0.52150224 2.78477911
##
## Number of Observations: 581
## Number of Groups: 83These results are slightly different but similar to those obtained in the direct fitting.
Janss-Bayley model in teh SEM framework using
lavaan
Specifying the model
hght.jb.model <- '
#latent variable definitions
#defining the intercept
eta_1 =~ 1*hght01 +
1*hght03 +
1*hght06 +
1*hght09 +
1*hght12 +
1*hght15 +
1*hght18 +
1*hght24 +
1*hght36
#defining the slope
eta_2 =~ .0833*hght01 +
.25*hght03 +
.5*hght06 +
.75*hght09 +
1*hght12 +
1.25*hght15 +
1.5*hght18 +
2*hght24 +
3*hght36
#defining eta_3
eta_3 =~ start(-.15)*L13*hght01 +
start(-.40)*L23*hght03 +
start(-.64)*L33*hght06 +
start(-.78)*L43*hght09 +
start(-.87)*L53*hght12 +
start(-.92)*L63*hght15 +
start(-.95)*L73*hght18 +
start(-.98)*L83*hght24 +
start(-.99)*L93*hght36
#defining eta_4
eta_4 =~ start(-1.25)*L14*hght01 +
start(-2.65)*L24*hght03 +
start(-3.16)*L34*hght06 +
start(-2.82)*L44*hght09 +
start(-2.24)*L54*hght12 +
start(-1.66)*L64*hght15 +
start(-1.19)*L74*hght18 +
start(-0.56)*L84*hght24 +
start(-0.11)*L94*hght36
#latent means
eta_1 ~ start(51)*1
eta_2 ~ start(9.25)*1
eta_3 ~ start(-17.8)*alpha3*1
eta_4 ~ 0 # start(0.07)*1
#introducing phantom factor, setting all parameters to 0 and mean to alpha4
phantom =~ 0*hght01
phantom ~~ 0*phantom
phantom ~~ 0*eta_1
phantom ~~ 0*eta_2
phantom ~~ 0*eta_3
phantom ~~ 0*eta_4
phantom ~ start(-2)*alpha4*1
#contraints to define factor loadings
alpha4 < -2 #add this constraint in so that it finds the viable solution
L13 == exp(alpha4*0.0833)-1
L23 == exp(alpha4*0.25)-1
L33 == exp(alpha4*0.50)-1
L43 == exp(alpha4*0.75)-1
L53 == exp(alpha4*1.00)-1
L63 == exp(alpha4*1.25)-1
L73 == exp(alpha4*1.50)-1
L83 == exp(alpha4*2.00)-1
L93 == exp(alpha4*3.00)-1
L14 == alpha3 * 0.0833 * exp(alpha4*0.0833)
L24 == alpha3 * 0.25 * exp(alpha4*0.25)
L34 == alpha3 * 0.50 * exp(alpha4*0.50)
L44 == alpha3 * 0.75 * exp(alpha4*0.75)
L54 == alpha3 * 1.00 * exp(alpha4*1.00)
L64 == alpha3 * 1.25 * exp(alpha4*1.25)
L74 == alpha3 * 1.50 * exp(alpha4*1.50)
L84 == alpha3 * 2.00 * exp(alpha4*2.00)
L94 == alpha3 * 3.00 * exp(alpha4*3.00)
#factor variances
eta_1 ~~ start(7.4)*eta_1
eta_2 ~~ start(0.81)*eta_2
eta_3 ~~ start(11.63)*eta_3
eta_4 ~~ start(0.08)*eta_4
#factor covariances
eta_1 ~~ start(0.64)*eta_2
eta_1 ~~ start(4.45)*eta_3
eta_1 ~~ start(0.03)*eta_4
eta_2 ~~ start(1.35)*eta_3
eta_2 ~~ start(-0.08)*eta_4
eta_3 ~~ start(-0.49)*eta_4
#manifest variances
hght01 ~~ start(.64)*theta*hght01
hght03 ~~ theta*hght03
hght06 ~~ theta*hght06
hght09 ~~ theta*hght09
hght12 ~~ theta*hght12
hght15 ~~ theta*hght15
hght18 ~~ theta*hght18
hght24 ~~ theta*hght24
hght36 ~~ theta*hght36
#manifest means (fixed to zero)
hght01 ~ 0*1
hght03 ~ 0*1
hght06 ~ 0*1
hght09 ~ 0*1
hght12 ~ 0*1
hght15 ~ 0*1
hght18 ~ 0*1
hght24 ~ 0*1
hght36 ~ 0*1'Fitting the model object
hght.jb.lavaan <- lavaan(hght.jb.model,
data = hght_wide,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE,
mimic="mplus",
control=list(iter.max=500),
verbose=FALSE)
summary(hght.jb.lavaan, fit.measures=TRUE)## lavaan 0.6-9 ended normally after 1865 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 41
## Number of inequality constraints 1
##
## Number of observations 83
## Number of missing patterns 31
##
## Model Test User Model:
##
## Test statistic 103.452
## Degrees of freedom 39
## 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.930
## Tucker-Lewis Index (TLI) 0.935
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -987.467
## Loglikelihood unrestricted model (H1) -935.741
##
## Akaike (AIC) 2004.935
## Bayesian (BIC) 2041.217
## Sample-size adjusted Bayesian (BIC) 1993.904
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.141
## 90 Percent confidence interval - lower 0.108
## 90 Percent confidence interval - upper 0.174
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.306
##
## 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.159 0.005 -28.939 0.000
## hght03 (L23) -0.405 0.012 -34.751 0.000
## hght06 (L33) -0.646 0.014 -46.575 0.000
## hght09 (L43) -0.789 0.012 -63.764 0.000
## hght12 (L53) -0.875 0.010 -89.056 0.000
## hght15 (L63) -0.925 0.007 -126.680 0.000
## hght18 (L73) -0.956 0.005 -183.199 0.000
## hght24 (L83) -0.984 0.002 -399.685 0.000
## hght36 (L93) -0.998 0.000 -2154.943 0.000
## eta_4 =~
## hght01 (L14) -1.250 0.038 -32.517 0.000
## hght03 (L24) -2.653 0.107 -24.720 0.000
## hght06 (L34) -3.158 0.182 -17.378 0.000
## hght09 (L44) -2.819 0.214 -13.184 0.000
## hght12 (L54) -2.236 0.212 -10.563 0.000
## hght15 (L64) -1.663 0.189 -8.792 0.000
## hght18 (L74) -1.188 0.158 -7.520 0.000
## hght24 (L84) -0.561 0.096 -5.825 0.000
## hght36 (L94) -0.105 0.026 -4.008 0.000
## phantom =~
## hght01 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## eta_1 ~~
## phantom 0.000
## eta_2 ~~
## phantom 0.000
## eta_3 ~~
## phantom 0.000
## eta_4 ~~
## phantom 0.000
## eta_1 ~~
## eta_2 0.639 0.528 1.212 0.226
## eta_3 4.448 1.572 2.829 0.005
## eta_4 0.030 0.262 0.113 0.910
## eta_2 ~~
## eta_3 1.346 0.849 1.586 0.113
## eta_4 -0.078 0.143 -0.542 0.588
## eta_3 ~~
## eta_4 -0.488 0.383 -1.274 0.203
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## eta_1 51.076 0.349 146.464 0.000
## eta_2 9.296 0.168 55.402 0.000
## eta_3 (alp3) -17.836 0.481 -37.051 0.000
## eta_4 0.000
## phantom (alp4) -2.076 0.078 -26.508 0.000
## .hght01 0.000
## .hght03 0.000
## .hght06 0.000
## .hght09 0.000
## .hght12 0.000
## .hght15 0.000
## .hght18 0.000
## .hght24 0.000
## .hght36 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## phantom 0.000
## eta_1 7.410 1.516 4.889 0.000
## eta_2 0.809 0.340 2.380 0.017
## eta_3 11.626 3.006 3.868 0.000
## eta_4 0.076 0.079 0.959 0.338
## .hght01 (thet) 0.644 0.053 12.175 0.000
## .hght03 (thet) 0.644 0.053 12.175 0.000
## .hght06 (thet) 0.644 0.053 12.175 0.000
## .hght09 (thet) 0.644 0.053 12.175 0.000
## .hght12 (thet) 0.644 0.053 12.175 0.000
## .hght15 (thet) 0.644 0.053 12.175 0.000
## .hght18 (thet) 0.644 0.053 12.175 0.000
## .hght24 (thet) 0.644 0.053 12.175 0.000
## .hght36 (thet) 0.644 0.053 12.175 0.000
##
## Constraints:
## |Slack|
## -2 - (alpha4) 0.076
## L13 - (exp(alpha4*0.0833)-1) 0.000
## L23 - (exp(alpha4*0.25)-1) 0.000
## L33 - (exp(alpha4*0.50)-1) 0.000
## L43 - (exp(alpha4*0.75)-1) 0.000
## L53 - (exp(alpha4*1.00)-1) 0.000
## L63 - (exp(alpha4*1.25)-1) 0.000
## L73 - (exp(alpha4*1.50)-1) 0.000
## L83 - (exp(alpha4*2.00)-1) 0.000
## L93 - (exp(alpha4*3.00)-1) 0.000
## L14 - (alpha3*0.0833*exp(alpha4*0.0833)) 0.000
## L24 - (alpha3*0.25*exp(alpha4*0.25)) 0.000
## L34 - (alpha3*0.50*exp(alpha4*0.50)) 0.000
## L44 - (alpha3*0.75*exp(alpha4*0.75)) 0.000
## L54 - (alpha3*1.00*exp(alpha4*1.00)) 0.000
## L64 - (alpha3*1.25*exp(alpha4*1.25)) 0.000
## L74 - (alpha3*1.50*exp(alpha4*1.50)) 0.000
## L84 - (alpha3*2.00*exp(alpha4*2.00)) 0.000
## L94 - (alpha3*3.00*exp(alpha4*3.00)) 0.000This matches Output 12.7 in the book! Yay!
We can try to make a diagram to check if we specified and estimated the model as intended, but the semPaths package does not know what to do witht eh phantom variables and constraints.
#diagram of fitted model
#semPaths(hght.jb.lavaan,what = "path", whatLabels = "par")Bilinear Spline Growth Model with Variation in the Knot Point.
Bilinear Spline Growth Model with Variation in the Knot Point in MLM
framework with nlme
Specifying and fitting the full spline model
hght.spline.nlme <- nlme(hght ~ b_1i + b_2i * (pmin(0, age - b_4i)) +
b_3i * (pmax(0, age - b_4i)),
data = hght_long,
fixed = b_1i + b_2i + b_3i + b_4i ~ 1,
random = b_1i + b_2i + b_3i + b_4i ~ 1,
groups =~ id,
start = c(60, 5, 2, 8),
na.action = na.omit)## Warning in nlme.formula(hght ~ b_1i + b_2i * (pmin(0, age - b_4i)) + b_3i
## * : Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase
## 'msMaxIter'!## Warning in nlme.formula(hght ~ b_1i + b_2i * (pmin(0, age - b_4i)) + b_3i
## * : Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase
## 'msMaxIter'!summary(hght.spline.nlme) ## Nonlinear mixed-effects model fit by maximum likelihood
## Model: hght ~ b_1i + b_2i * (pmin(0, age - b_4i)) + b_3i * (pmax(0, age - b_4i))
## Data: hght_long
## AIC BIC logLik
## 2293.05 2358.521 -1131.525
##
## Random effects:
## Formula: list(b_1i ~ 1, b_2i ~ 1, b_3i ~ 1, b_4i ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b_1i 2.70815098 b_1i b_2i b_3i
## b_2i 0.28583819 0.108
## b_3i 0.05569232 0.475 0.445
## b_4i 0.29285714 0.794 -0.504 0.233
## Residual 1.21256519
##
## Fixed effects: b_1i + b_2i + b_3i + b_4i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 71.92526 0.3602998 495 199.62612 0
## b_2i 2.51075 0.0540562 495 46.44704 0
## b_3i 0.90033 0.0098081 495 91.79507 0
## b_4i 7.72069 0.1427210 495 54.09638 0
## Correlation:
## b_1i b_2i b_3i
## b_2i -0.153
## b_3i 0.027 0.168
## b_4i 0.617 -0.685 -0.190
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.638686277 -0.530024577 0.009102385 0.573875184 2.623835752
##
## Number of Observations: 581
## Number of Groups: 83Plotting the predicted trajectories. Note that this is not totally accurate - as it does not quite get the splines constructed properly. Rather than just using the model object predicted scores, we should generate fine-grained data based on the parameter estimates so that we can see the differences in knot-point and the slopes of each individual’s predicted trajectory more precisely.
#obtaining predicted scores for individuals
hght_long$pred_spline <- predict(hght.spline.nlme, newdata=hght_long)
#obtaining predicted scores for prototype
hght_long$proto_spline <- predict(hght.spline.nlme, newdata=hght_long, level=0)
#plotting predicted trajectories
#intraindividual change trajetories
ggplot(data = hght_long, aes(x = age, y = pred_spline, group = id)) +
#geom_point(color="black") +
geom_line(color="black") +
geom_line(aes(x = age, y = proto_spline), color="red",size=2) +
theme_classic(base_size = 18) +
theme(legend.position = "none") +
labs(title = "Precicted Individual Height Trajectories", y = "Height (cm)", x = "Age (years)")
###Fitting the spline growth model with estimated knot point in teh
SEM framework using lavaan.
Note: this model results in estimation issues (e.g., variance of the knot point was estimated to be negative), thus the model results should not be interpreted. Specifying the model
hght.spline.model <- '
#latent variable definitions
#defining the intercept
eta_1 =~ 1*hght01 +
1*hght03 +
1*hght06 +
1*hght09 +
1*hght12 +
1*hght15 +
1*hght18 +
1*hght24 +
1*hght36
#defining the slope
eta_2 =~ 1*L12*hght01 +
3*L22*hght03 +
6*L32*hght06 +
9*L42*hght09 +
12*L52*hght12 +
15*L62*hght15 +
18*L72*hght18 +
24*L82*hght24 +
36*L92*hght36
#defining eta_3
eta_3 =~ NA*L13*hght01 +
L23*hght03 +
L33*hght06 +
L43*hght09 +
L53*hght12 +
L63*hght15 +
L73*hght18 +
L83*hght24 +
L93*hght36
#defining eta_4
eta_4 =~ start(-0.90)*L14*hght01 +
start(-0.90)*L24*hght03 +
start(-0.90)*L34*hght06 +
start(-2.57)*L44*hght09 +
start(-2.57)*L54*hght12 +
start(-2.57)*L64*hght15 +
start(-2.57)*L74*hght18 +
start(-2.57)*L84*hght24 +
start(-2.57)*L94*hght36
#latent means
eta_1 ~ start(71)*1
eta_2 ~ start(1.7)*alpha2*1
eta_3 ~ start(-0.8)*alpha3*1
eta_4 ~ start(7.7)*alpha4*1
#factor variances
#fixing variances to zero b/c getting negative estimates
eta_1 ~~ start(0.01)*eta_1
eta_2 ~~ 0*eta_2
eta_3 ~~ start(0.01)*eta_3
eta_4 ~~ start(0.08)*eta_4
#factor covariances
eta_1 ~~ start(.08)*eta_2
eta_1 ~~ 0*eta_3
eta_1 ~~ start(.03)*eta_4
eta_2 ~~ 0*eta_3
eta_2 ~~ start(.40)*eta_4
eta_3 ~~ 0*eta_4
#introducing phantom factors, setting all parameters to 0 and mean to alpha4
phantom2 =~ 0*hght01
phantom2 ~~ 0*phantom2
phantom2 ~~ 0*eta_1
phantom2 ~~ 0*eta_2
phantom2 ~~ 0*eta_3
phantom2 ~~ 0*eta_4
phantom2 ~ start(.29)*alpha2*1
phantom3 =~ 0*hght01
phantom3 ~~ 0*phantom3
phantom3 ~~ 0*phantom2
phantom3 ~~ 0*eta_1
phantom3 ~~ 0*eta_2
phantom3 ~~ 0*eta_3
phantom3 ~~ 0*eta_4
phantom3 ~ start(.29)*alpha3*1
phantom4 =~ 0*hght01
phantom4 ~~ 0*phantom4
phantom4 ~~ 0*phantom3
phantom4 ~~ 0*phantom2
phantom4 ~~ 0*eta_1
phantom4 ~~ 0*eta_2
phantom4 ~~ 0*eta_3
phantom4 ~~ 0*eta_4
phantom4 ~ start(.29)*alpha4*1
#contraints to define factor loadings
L12 == 1 - alpha4
L22 == 3 - alpha4
L32 == 6 - alpha4
L42 == 9 - alpha4
L52 == 12 - alpha4
L62 == 15 - alpha4
L72 == 18 - alpha4
L82 == 24 - alpha4
L92 == 36 - alpha4
L13 == sqrt(( 1 - alpha4)^2)
L23 == sqrt(( 3 - alpha4)^2)
L33 == sqrt(( 6 - alpha4)^2)
L43 == sqrt(( 9 - alpha4)^2)
L53 == sqrt((12 - alpha4)^2)
L63 == sqrt((15 - alpha4)^2)
L73 == sqrt((18 - alpha4)^2)
L83 == sqrt((24 - alpha4)^2)
L93 == sqrt((36 - alpha4)^2)
L14 == (-alpha2 * sqrt(( 1 - alpha4)^2) + alpha3 * 1 - alpha3 * alpha4) / sqrt(( 1 - alpha4)^2)
L24 == (-alpha2 * sqrt(( 3 - alpha4)^2) + alpha3 * 3 - alpha3 * alpha4) / sqrt(( 3 - alpha4)^2)
L34 == (-alpha2 * sqrt(( 6 - alpha4)^2) + alpha3 * 6 - alpha3 * alpha4) / sqrt(( 6 - alpha4)^2)
L44 == (-alpha2 * sqrt(( 9 - alpha4)^2) + alpha3 * 9 - alpha3 * alpha4) / sqrt(( 9 - alpha4)^2)
L54 == (-alpha2 * sqrt((12 - alpha4)^2) + alpha3 * 12 - alpha3 * alpha4) / sqrt((12 - alpha4)^2)
L64 == (-alpha2 * sqrt((15 - alpha4)^2) + alpha3 * 15 - alpha3 * alpha4) / sqrt((15 - alpha4)^2)
L74 == (-alpha2 * sqrt((18 - alpha4)^2) + alpha3 * 18 - alpha3 * alpha4) / sqrt((18 - alpha4)^2)
L84 == (-alpha2 * sqrt((24 - alpha4)^2) + alpha3 * 24 - alpha3 * alpha4) / sqrt((24 - alpha4)^2)
L94 == (-alpha2 * sqrt((36 - alpha4)^2) + alpha3 * 36 - alpha3 * alpha4) / sqrt((36 - alpha4)^2)
#manifest variances
hght01 ~~ start(1.2)*theta*hght01
hght03 ~~ theta*hght03
hght06 ~~ theta*hght06
hght09 ~~ theta*hght09
hght12 ~~ theta*hght12
hght15 ~~ theta*hght15
hght18 ~~ theta*hght18
hght24 ~~ theta*hght24
hght36 ~~ theta*hght36
#manifest means (fixed to zero)
hght01 ~ 0*1
hght03 ~ 0*1
hght06 ~ 0*1
hght09 ~ 0*1
hght12 ~ 0*1
hght15 ~ 0*1
hght18 ~ 0*1
hght24 ~ 0*1
hght36 ~ 0*1'Fitting the model
hght.spline.lavaan <- lavaan(hght.spline.model,
data = hght_wide,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE,
mimic="mplus",
control=list(iter.max=500),
verbose=FALSE)## Warning in lav_start_check_cov(lavpartable = lavpartable, start = START): lavaan WARNING: non-zero covariance element set to zero, due to fixed-to-zero variances
## variables involved are: eta_1 eta_2## Warning in lav_start_check_cov(lavpartable = lavpartable, start = START): lavaan WARNING: starting values imply a correlation larger than 1;
## variables involved are: eta_1 eta_4## Warning in lav_start_check_cov(lavpartable = lavpartable, start = START): lavaan WARNING: non-zero covariance element set to zero, due to fixed-to-zero variances
## variables involved are: eta_2 eta_4## Warning in lavaan(hght.spline.model, data = hght_wide, meanstructure = TRUE, : lavaan WARNING:
## the optimizer warns that a solution has NOT been found!summary(hght.spline.lavaan, fit.measures=TRUE)## lavaan 0.6-9 did NOT end normally after 444 iterations
## ** WARNING ** Estimates below are most likely unreliable
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 49
##
## Number of observations 83
## Number of missing patterns 31
##
## Model Test User Model:
##
## Test statistic NA
## Degrees of freedom NA## Warning in .local(object, ...): lavaan WARNING: fit measures not available if model did not converge##
## 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 (L12) -16.959 NA
## hght03 (L22) -15.024 NA
## hght06 (L32) -12.049 NA
## hght09 (L42) -9.011 NA
## hght12 (L52) -6.031 NA
## hght15 (L62) -3.001 NA
## hght18 (L72) 0.069 NA
## hght24 (L82) 5.961 NA
## hght36 (L92) 18.024 NA
## eta_3 =~
## hght01 (L13) 17.016 NA
## hght03 (L23) 14.993 NA
## hght06 (L33) 11.984 NA
## hght09 (L43) 9.013 NA
## hght12 (L53) 5.963 NA
## hght15 (L63) 2.963 NA
## hght18 (L73) 0.008 NA
## hght24 (L83) 5.968 NA
## hght36 (L93) 18.021 NA
## eta_4 =~
## hght01 (L14) 1.621 NA
## hght03 (L24) 1.522 NA
## hght06 (L34) 1.397 NA
## hght09 (L44) 1.316 NA
## hght12 (L54) 1.229 NA
## hght15 (L64) 1.172 NA
## hght18 (L74) 1.084 NA
## hght24 (L84) 0.691 NA
## hght36 (L94) 0.119 NA
## phantom2 =~
## hght01 0.000
## phantom3 =~
## hght01 0.000
## phantom4 =~
## hght01 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## eta_1 ~~
## eta_2 22.517 NA
## eta_3 0.000
## eta_4 30.717 NA
## eta_2 ~~
## eta_3 0.000
## eta_4 -14.058 NA
## eta_3 ~~
## eta_4 0.000
## eta_1 ~~
## phantom2 0.000
## eta_2 ~~
## phantom2 0.000
## eta_3 ~~
## phantom2 0.000
## eta_4 ~~
## phantom2 0.000
## phantom2 ~~
## phantom3 0.000
## eta_1 ~~
## phantom3 0.000
## eta_2 ~~
## phantom3 0.000
## eta_3 ~~
## phantom3 0.000
## eta_4 ~~
## phantom3 0.000
## phantom3 ~~
## phantom4 0.000
## phantom2 ~~
## phantom4 0.000
## eta_1 ~~
## phantom4 0.000
## eta_2 ~~
## phantom4 0.000
## eta_3 ~~
## phantom4 0.000
## eta_4 ~~
## phantom4 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## eta_1 53.909 NA
## eta_2 (alp2) -0.833 NA
## eta_3 (alp3) -0.495 NA
## eta_4 (alp4) 18.000 NA
## phantm2 (alp2) -0.869 NA
## phantm3 (alp3) -0.508 NA
## phantm4 (alp4) 17.998 NA
## .hght01 0.000
## .hght03 0.000
## .hght06 0.000
## .hght09 0.000
## .hght12 0.000
## .hght15 0.000
## .hght18 0.000
## .hght24 0.000
## .hght36 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## eta_1 -14.780 NA
## eta_2 0.000
## eta_3 0.185 NA
## eta_4 10.817 NA
## phantm2 0.000
## phantm3 0.000
## phantm4 0.000
## .hght01 (thet) 7.413 NA
## .hght03 (thet) 7.478 NA
## .hght06 (thet) 7.388 NA
## .hght09 (thet) 7.391 NA
## .hght12 (thet) 7.383 NA
## .hght15 (thet) 7.415 NA
## .hght18 (thet) 7.409 NA
## .hght24 (thet) 7.382 NA
## .hght36 (thet) 7.419 NA
##
## Constraints:
## |Slack|
## L12 - (1-alpha4) 0.041
## L22 - (3-alpha4) 0.024
## L32 - (6-alpha4) 0.049
## L42 - (9-alpha4) 0.011
## L52 - (12-alpha4) 0.031
## L62 - (15-alpha4) 0.001
## L72 - (18-alpha4) 0.069
## L82 - (24-alpha4) 0.039
## L92 - (36-alpha4) 0.024
## L13 - (sqrt((1-alpha4)^2)) 0.016
## L23 - (sqrt((3-alpha4)^2)) 0.007
## L33 - (sqrt((6-alpha4)^2)) 0.016
## L43 - (sqrt((9-alpha4)^2)) 0.013
## L53 - (sqrt((12-alpha4)^2)) 0.037
## L63 - (sqrt((15-alpha4)^2)) 0.037
## L73 - (sqrt((18-alpha4)^2)) 0.008
## L83 - (sqrt((24-alpha4)^2)) 0.032
## L93 - (sqrt((36-alpha4)^2)) 0.021
## L14-((-lph2*((1-4)^2)+3*1-3*4)/((1-4)^2)) 0.293
## L24-((-lph2*((3-4)^2)+3*3-3*4)/((3-4)^2)) 0.193
## L34-((-lph2*((6-4)^2)+3*6-3*4)/((6-4)^2)) 0.069
## L44-((-lph2*((9-4)^2)+3*9-3*4)/((9-4)^2)) 0.012
## L54-((-2*((12-4)^2)+3*12-3*4)/((12-4)^2)) 0.099
## L64-((-2*((15-4)^2)+3*15-3*4)/((15-4)^2)) 0.156
## L74-((-2*((18-4)^2)+3*18-3*4)/((18-4)^2)) 0.244
## L84-((-2*((24-4)^2)+3*24-3*4)/((24-4)^2)) 0.353
## L94-((-2*((36-4)^2)+3*36-3*4)/((36-4)^2)) 0.219Something is wrong in the estimation of the model. After some hours of trying to sort out what the issue are, it was time to post the code. So it is getting posted with the caveat that additional trouble-shooting needs to be done.
Conclusion
This code has demonstrated some of the possibiites for moving into more complex models with interesting functional forms and nonlinearity in how the interindividual differences manifest in the trajectories of change.
WOWSERS!