Bivariate Growth Model - Multilevel & SEM Implementation in R
Overview
This tutorial illustrates fitting of multivariate (bivariate) linear growth models in the multilevel and SEM frameworks in R.
Example data and code are drawn from Chapter 8 of Grimm, Ram, and Estabrook (2017). Specifically, using the NLSY-CYA Dataset we examine how individual differences in change in children’s mathematics achievement across grade are related to individual differences in change in children’s hyperactivity (as rated by teachers) across grade. Please see the book chapter for additional interpretations and insights about the analyses.
Preliminaries
Loading libraries used in this script.
Data Description
For our examples, we use the mathematics achievement scores from the NLSY-CYA Long Data.
Load the repeated measures data (long format)
#set filepath for data file
filepath <- "https://raw.githubusercontent.com/The-Change-Lab/collaborations/refs/heads/main/GrowthModeling/nlsy_math_hyp_long_R.dat"
#read in the text data file using the url() function
dat <- read.table(file=url(filepath),
na.strings = ".") #indicates the missing data designator
#copy data with new name
nlsy_math_hyp_long <- dat
#Add names the columns of the data set
names(nlsy_math_hyp_long) = c('id', 'female', 'lb_wght', 'anti_k1',
'math', 'comp', 'rec', 'bpi', 'as', 'anx', 'hd',
'hyp', 'dp', 'wd',
'grade', 'occ', 'age', 'men', 'spring', 'anti')
#reducing to variables of interest
nlsy_math_hyp_long <- nlsy_math_hyp_long %>%
select(id, grade, math, hyp)
#view the first few observations in the data set
head(nlsy_math_hyp_long, 10)## id grade math hyp
## 1 201 3 38 0
## 2 201 5 55 0
## 3 303 2 26 1
## 4 303 5 33 1
## 5 2702 2 56 2
## 6 2702 4 58 3
## 7 2702 8 80 3
## 8 4303 3 41 1
## 9 4303 4 58 1
## 10 5002 4 46 3Our specific interest is intraindividual change in the repeated
measures of math and hyp across
grade.
As noted in Chapter 2 , it is important to plot the data to obtain a better understanding of the structure and form of the observed phenomenon.
Longitudinal Plot of Math across Grade at Testing
#intraindividual change trajectories
nlsy_math_hyp_long %>% #dataset
ggplot(aes(x = grade, y = math, group = id)) + #setting variables
geom_point(size=.5) + #adding points to plot
geom_line() + #adding lines to plot
theme_bw() + #changing style/background
#setting the x-axis with breaks and labels
scale_x_continuous(limits=c(2,8),
breaks = c(2,3,4,5,6,7,8),
name = "Grade at Testing") +
#setting the y-axis with limits breaks and labels
scale_y_continuous(limits=c(10,90),
breaks = c(10,30,50,70,90),
name = "PIAT Mathematics")
Longitudinal Plot of Hyperactivity across Grade at Testing
#intraindividual change trajectories
nlsy_math_hyp_long %>% #dataset
ggplot(aes(x = grade, y = hyp, group = id)) + #setting variables
geom_point(size=.5) + #adding points to plot
geom_line() + #adding lines to plot
theme_bw() + #changing style/background
#setting the x-axis with breaks and labels
scale_x_continuous(limits=c(2,8),
breaks = c(2,3,4,5,6,7,8),
name = "Grade at Testing") +
#setting the y-axis with limits breaks and labels
scale_y_continuous(limits=c(0,5),
breaks = c(0,1,2,3,4,5),
name = "Hyperactivity")
Data Preparation
For the multilevel model implementation we need to collect both
math and hyp into a single outcome variable
along with dummy variables that indicate which variable is represented.
There are lots of ways to do this. We use a simple approach here. There
are more elegant ways, too.
# creating a stacked multivariate form of the data
#make math data
math <- data.frame (id = nlsy_math_hyp_long$id,
var = nlsy_math_hyp_long$math,
grade = nlsy_math_hyp_long$grade,
d_math = 1,
d_hyp = 0,
grp = 'math')
#make hyp data
hyp <- data.frame (id = nlsy_math_hyp_long$id,
var = nlsy_math_hyp_long$hyp,
grade = nlsy_math_hyp_long$grade,
d_math = 0,
d_hyp = 1,
grp = 'hyp')
#rowbinding the data
nlsy_multi_long <- rbind(math,hyp)
#reordering for easy viewing
nlsy_multi_long <- nlsy_multi_long[order(nlsy_multi_long$id,
nlsy_multi_long$grade,
nlsy_multi_long$d_hyp), ]
#view the first few observations in the data set
head(nlsy_multi_long, 20)## id var grade d_math d_hyp grp
## 1 201 38 3 1 0 math
## 2222 201 0 3 0 1 hyp
## 2 201 55 5 1 0 math
## 2223 201 0 5 0 1 hyp
## 3 303 26 2 1 0 math
## 2224 303 1 2 0 1 hyp
## 4 303 33 5 1 0 math
## 2225 303 1 5 0 1 hyp
## 5 2702 56 2 1 0 math
## 2226 2702 2 2 0 1 hyp
## 6 2702 58 4 1 0 math
## 2227 2702 3 4 0 1 hyp
## 7 2702 80 8 1 0 math
## 2228 2702 3 8 0 1 hyp
## 8 4303 41 3 1 0 math
## 2229 4303 1 3 0 1 hyp
## 9 4303 58 4 1 0 math
## 2230 4303 1 4 0 1 hyp
## 10 5002 46 4 1 0 math
## 2231 5002 3 4 0 1 hypNote the alternating row structure of the data with both
math and hyp in a single variable
var.
For the SEM implementation we need wide data with
math2-math8 and hyp2-hym8 all in one row per
person.
Load the repeated measures data (wide format)
#set filepath for data file
filepath <- "https://raw.githubusercontent.com/The-Change-Lab/collaborations/refs/heads/main/GrowthModeling/nlsy_math_hyp_wide_R.dat"
#read in the text data file using the url() function
dat <- read.table(file=url(filepath),
na.strings = ".") #indicates the missing data designator
#copy data with new name
nlsy_math_hyp_wide <- dat
#Add names the columns of the data set
names(nlsy_math_hyp_wide)<-c('id','female','lb_wght','anti_k1',
'math2','math3','math4','math5','math6','math7','math8',
'comp2','comp3','comp4','comp5','comp6','comp7','comp8',
'rec2','rec3','rec4','rec5','rec6','rec7','rec8',
'bpi2','bpi3','bpi4','bpi5','bpi6','bpi7','bpi8',
'asl2','asl3','asl4','asl5','asl6','asl7','asl8',
'ax2','ax3','ax4','ax5','ax6','ax7','ax8',
'hds2','hds3','hds4','hds5','hds6','hds7','hds8',
'hyp2','hyp3','hyp4','hyp5','hyp6','hyp7','hyp8',
'dpn2','dpn3','dpn4','dpn5','dpn6','dpn7','dpn8',
'wdn2','wdn3','wdn4','wdn5','wdn6','wdn7','wdn8',
'age2','age3','age4','age5','age6','age7','age8',
'men2','men3','men4','men5','men6','men7','men8',
'spring2','spring3','spring4',
'spring5','spring6','spring7','spring8',
'anti2','anti3','anti4','anti5','anti6','anti7','anti8')
#reducing to variables of interest
nlsy_multi_wide <- nlsy_math_hyp_wide %>%
select(id, math2, math3, math4, math5, math6, math7, math8,
hyp2, hyp3, hyp4, hyp5, hyp6, hyp7, hyp8)
#view the first few observations in the data set
head(nlsy_multi_wide, 10)## id math2 math3 math4 math5 math6 math7 math8 hyp2 hyp3 hyp4 hyp5 hyp6 hyp7
## 1 201 NA 38 NA 55 NA NA NA NA 0 NA 0 NA NA
## 2 303 26 NA NA 33 NA NA NA 1 NA NA 1 NA NA
## 3 2702 56 NA 58 NA NA NA 80 2 NA 3 NA NA NA
## 4 4303 NA 41 58 NA NA NA NA NA 1 1 NA NA NA
## 5 5002 NA NA 46 NA 54 NA 66 NA NA 3 NA 2 NA
## 6 5005 35 NA 50 NA 60 NA 59 0 NA 3 NA 0 NA
## 7 5701 NA 62 61 NA NA NA NA NA 4 3 NA NA NA
## 8 6102 NA NA 55 67 NA 81 NA NA NA 2 0 NA 0
## 9 6801 NA 54 NA 62 NA 66 NA NA 0 NA 1 NA 1
## 10 6802 NA 55 NA 66 NA 68 NA NA 0 NA 0 NA 0
## hyp8
## 1 NA
## 2 NA
## 3 3
## 4 NA
## 5 3
## 6 1
## 7 NA
## 8 NA
## 9 NA
## 10 NAWe are now ready to rock!
Multilevel: Bivariate Growth Model using nlme()
function
To capture the systematic growth in the mathematics scores and
hyperactivity scores over grade, we use the grade variable,
so that the rate of growth indicates (linear) school-related
change in mathematics scores.
We use the multivariate (stacked) data in “extralong” format.
Then we fit the model with nlme using our expansive coding style. This is a difficult model to estimate so we add control parameters to the model (change iteration number to 200).
# bivariate growth model with residual covariance
bivar_nlme <- nlme(var~d_math*(b_1i+b_2i*(grade-2))+
d_hyp *(h_1i+h_2i*(grade-2)),
data = nlsy_multi_long,
fixed = b_1i+b_2i+h_1i+h_2i~1,
random = b_1i+b_2i+h_1i+h_2i~1,
group = ~id,
start = c(35, 4, 1, -1),
#invoking separate residual variance for each variable
weights = varIdent(c(hyp=.3), form = ~1|grp),
# invokes off-diagonal symmetric structure for level 1 residuals
corr=corSymm(form = ~1 | id/grade),
na.action = na.exclude,
control=lmeControl(maxIter = 500, msMaxIter = 200, tolerance = 1e-4,
niter = 100, msTol = 1e-5, nlmStepMax = 500,
msVerbose = FALSE,
returnObject = TRUE)) #to increase iterations
#obtaining summary of the model using the object we just created
summary(bivar_nlme)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: var ~ d_math * (b_1i + b_2i * (grade - 2)) + d_hyp * (h_1i + h_2i * (grade - 2))
## Data: nlsy_multi_long
## AIC BIC logLik
## 23614.95 23723.54 -11790.48
##
## Random effects:
## Formula: list(b_1i ~ 1, b_2i ~ 1, h_1i ~ 1, h_2i ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b_1i 8.04429094 b_1i b_2i h_1i
## b_2i 0.86686900 -0.029
## h_1i 1.24180244 -0.298 0.091
## h_2i 0.06805571 0.196 -0.681 -0.235
## Residual 6.01054587
##
## Correlation Structure: General
## Formula: ~1 | id/grade
## Parameter estimate(s):
## Correlation:
## 1
## 2 -0.002
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | grp
## Parameter estimates:
## math hyp
## 1.0000000 0.1748212
## Fixed effects: b_1i + b_2i + h_1i + h_2i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 35.25868 0.3551098 3456 99.28951 0e+00
## b_2i 4.34300 0.0873776 3456 49.70380 0e+00
## h_1i 1.90348 0.0581509 3456 32.73344 0e+00
## h_2i -0.05660 0.0142232 3456 -3.97912 1e-04
## Correlation:
## b_1i b_2i h_1i
## b_2i -0.531
## h_1i -0.166 0.037
## h_2i 0.040 -0.067 -0.608
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.32770363 -0.53353287 -0.06380594 0.55396047 3.05031262
##
## Number of Observations: 4391
## Number of Groups: 932Cool that this all works!
Note that this is a slightly more precise model than is in the book (because it has residual covariance in it).
Note that here, as in the book, the slope-slope correlation looks to be highly negative. That is different than what was obtained in the other modeling frameworks.
Multilevel: Bivariate Growth Model using lme()
function
We can also fit the same model using lme() function with
a traditional (non-equation type) coding. This version is useful in some
cases (e.g. when extending to include categorical TICs).
#creating the new time variable explicitly
nlsy_multi_long$grade2c <- nlsy_multi_long$grade - 2
# bivariate growth model with residual covariance
# using dummy variables (d_math, d_hyp) to invoke intercepts and sloes
bivar_lme <- lme(var~ -1 + d_math + d_math:grade2c + #note the -1 to turn off the default intercept
d_hyp + d_hyp:grade2c, #note use of : to create the interaction terms
data = nlsy_multi_long,
random = ~ -1 + d_math + d_math:grade2c + d_hyp + d_hyp:grade2c | id,
#invoking separate residual variance for each variable
weights = varIdent(form = ~1|grp),
# invokes off-diagonal symmetric structure for level 1 residuals
corr=corSymm(form = ~1 | id/grade),
na.action = na.exclude,
control=lmeControl(maxIter = 500, msMaxIter = 200, tolerance = 1e-4,
niter = 100, msTol = 1e-5, nlmStepMax = 500,
msVerbose = FALSE,
returnObject = TRUE)) #to increase iterations
#obtaining summary of the model using the object we just created
summary(bivar_lme)## Linear mixed-effects model fit by REML
## Data: nlsy_multi_long
## AIC BIC logLik
## 23629.59 23738.15 -11797.79
##
## Random effects:
## Formula: ~-1 + d_math + d_math:grade2c + d_hyp + d_hyp:grade2c | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## d_math 8.0523554 d_math d_hyp d_mt:2
## d_hyp 1.2432375 -0.298
## d_math:grade2c 0.8714015 -0.032 0.090
## grade2c:d_hyp 0.0695732 0.191 -0.236 -0.662
## Residual 6.0105499
##
## Correlation Structure: General
## Formula: ~1 | id/grade
## Parameter estimate(s):
## Correlation:
## 1
## 2 -0.002
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | grp
## Parameter estimates:
## math hyp
## 1.0000000 0.1748186
## Fixed effects: var ~ -1 + d_math + d_math:grade2c + d_hyp + d_hyp:grade2c
## Value Std.Error DF t-value p-value
## d_math 35.25855 0.3551777 3456 99.27016 0e+00
## d_hyp 1.90351 0.0581641 3456 32.72657 0e+00
## d_math:grade2c 4.34324 0.0874137 3456 49.68602 0e+00
## grade2c:d_hyp -0.05661 0.0142302 3456 -3.97784 1e-04
## Correlation:
## d_math d_hyp d_mt:2
## d_hyp -0.165
## d_math:grade2c -0.531 0.037
## grade2c:d_hyp 0.039 -0.608 -0.067
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.3213835 -0.5333360 -0.0638431 0.5532024 3.0501886
##
## Number of Observations: 4391
## Number of Groups: 932Yay! - the results are the same as obtained above with nlme using this coding approach! - which is nice because this version extends more easily when including categorical TICs. But, as noted above the highly negative slope-slope correlation is different than what was obtained in the other modeling frameworks.
Multilevel: Bivariate Growth Model using brm()
function
We can also fit the same model in a Bayesian Multilevel modeling
framework using brm() function with a traditional
(non-equation type) coding. Here, we actually do not need the stacked
data - we can have multiple outcomes and equations, so we use the
original data instead. This version is useful because it requires less
data management, is easily extendible (e.g. when extending to include
categorical TICs), and has the benefits fo Bayesian inference.
#creating the new time variable explicitly
nlsy_math_hyp_long$grade2c <- nlsy_math_hyp_long$grade - 2
# bivariate growth model with residual covariance
#writing model for hyperactivity
bfhyp <- bf(hyp ~ 1 + grade2c + (1 + grade2c | pp | id))
#writing model for math
bfmath <- bf(math ~ 1 + grade2c + (1 + grade2c | pp | id))
#checking priors
(myprior <- get_prior(mvbf(bfhyp + bfmath + set_rescor(rescor = TRUE)),
data = nlsy_math_hyp_long, family = gaussian())) ## Warning: Rows containing NAs were excluded from the model.## prior class coef group resp dpar nlpar lb ub
## (flat) b
## lkj(1) cor
## lkj(1) cor id
## (flat) Intercept
## lkj(1) rescor
## (flat) b hyp
## (flat) b grade2c hyp
## student_t(3, 1, 2.5) Intercept hyp
## student_t(3, 0, 2.5) sd hyp 0
## student_t(3, 0, 2.5) sd id hyp 0
## student_t(3, 0, 2.5) sd grade2c id hyp 0
## student_t(3, 0, 2.5) sd Intercept id hyp 0
## student_t(3, 0, 2.5) sigma hyp 0
## (flat) b math
## (flat) b grade2c math
## student_t(3, 46, 11.9) Intercept math
## student_t(3, 0, 11.9) sd math 0
## student_t(3, 0, 11.9) sd id math 0
## student_t(3, 0, 11.9) sd grade2c id math 0
## student_t(3, 0, 11.9) sd Intercept id math 0
## student_t(3, 0, 11.9) sigma math 0
## source
## default
## default
## (vectorized)
## default
## default
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default#Combining the two formulas into one bivariate growth model
bivar_brm <- brm(bfhyp + bfmath + set_rescor(rescor = TRUE),
data = nlsy_math_hyp_long,
control = list(adapt_delta = 0.95),
chains = 4, iter=4000, thin=1,
cores=4,
seed = 20250428,
file='brms_bivar_growthmodel')
#checking convergence
plot(bivar_brm)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.## Family: MV(gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: hyp ~ 1 + grade2c + (1 + grade2c | pp | id)
## math ~ 1 + grade2c + (1 + grade2c | pp | id)
## Data: nlsy_math_hyp_long (Number of observations: 2170)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~id (Number of levels: 930)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(hyp_Intercept) 1.22 0.05 1.13 1.33 1.00
## sd(hyp_grade2c) 0.03 0.03 0.00 0.11 1.02
## sd(math_Intercept) 8.06 0.35 7.38 8.76 1.00
## sd(math_grade2c) 0.72 0.24 0.18 1.13 1.02
## cor(hyp_Intercept,hyp_grade2c) -0.11 0.37 -0.76 0.69 1.00
## cor(hyp_Intercept,math_Intercept) -0.26 0.05 -0.37 -0.16 1.02
## cor(hyp_grade2c,math_Intercept) -0.13 0.33 -0.75 0.46 1.17
## cor(hyp_Intercept,math_grade2c) -0.01 0.17 -0.36 0.33 1.01
## cor(hyp_grade2c,math_grade2c) -0.15 0.44 -0.88 0.71 1.01
## cor(math_Intercept,math_grade2c) 0.09 0.24 -0.25 0.71 1.01
## Bulk_ESS Tail_ESS
## sd(hyp_Intercept) 1733 2398
## sd(hyp_grade2c) 240 376
## sd(math_Intercept) 1545 3916
## sd(math_grade2c) 265 485
## cor(hyp_Intercept,hyp_grade2c) 6007 4331
## cor(hyp_Intercept,math_Intercept) 295 1887
## cor(hyp_grade2c,math_Intercept) 18 79
## cor(hyp_Intercept,math_grade2c) 2106 1281
## cor(hyp_grade2c,math_grade2c) 267 615
## cor(math_Intercept,math_grade2c) 488 881
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## hyp_Intercept 1.90 0.06 1.79 2.02 1.00 2768 5099
## math_Intercept 35.29 0.36 34.58 35.99 1.00 3855 5609
## hyp_grade2c -0.06 0.01 -0.08 -0.03 1.00 10411 6537
## math_grade2c 4.34 0.09 4.17 4.52 1.00 7477 6014
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_hyp 1.06 0.02 1.02 1.10 1.00 2109 3066
## sigma_math 6.07 0.16 5.75 6.38 1.01 632 2075
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(hyp,math) -0.02 0.03 -0.08 0.04 1.00 2043 3054
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
Yay! - that works - The results for the results are the same using the multivariate Bayesian multilevel model framework as we obtain in the other frameworks (see note above) - which is nice because, although it takes much longer to fit, it requires less data management. The slope-slope correlation of interest is not discernible from zero. We do see here, though, that the distribution of the hyperactivity variable (which cannot go below zero) is not accurate. We might consider using some different distributional assumptions.
SEM: Bivariate Growth Model using lavaan
For the SEM implementation, we use the mathematics achievement and hyperactivity scores and time-invariant covariates from the Wide Data.
Specifying bivariate linear growth model
#writing out linear growth model in full SEM way
bivariate_lavaan_model <- '
# latent variable definitions
#intercept for math
eta_1 =~ 1*math2
eta_1 =~ 1*math3
eta_1 =~ 1*math4
eta_1 =~ 1*math5
eta_1 =~ 1*math6
eta_1 =~ 1*math7
eta_1 =~ 1*math8
#linear slope for math
eta_2 =~ 0*math2
eta_2 =~ 1*math3
eta_2 =~ 2*math4
eta_2 =~ 3*math5
eta_2 =~ 4*math6
eta_2 =~ 5*math7
eta_2 =~ 6*math8
#intercept for hyp
eta_3 =~ 1*hyp2
eta_3 =~ 1*hyp3
eta_3 =~ 1*hyp4
eta_3 =~ 1*hyp5
eta_3 =~ 1*hyp6
eta_3 =~ 1*hyp7
eta_3 =~ 1*hyp8
#linear slope for hyp
eta_4 =~ 0*hyp2
eta_4 =~ 1*hyp3
eta_4 =~ 2*hyp4
eta_4 =~ 3*hyp5
eta_4 =~ 4*hyp6
eta_4 =~ 5*hyp7
eta_4 =~ 6*hyp8
# factor variances
eta_1 ~~ eta_1
eta_2 ~~ eta_2
eta_3 ~~ eta_3
eta_4 ~~ eta_4
# covariances among factors
eta_1 ~~ eta_2 + eta_3 + eta_4
eta_2 ~~ eta_3 + eta_4
eta_3 ~~ eta_4
# factor means
eta_1 ~ start(35)*1
eta_2 ~ start(4)*1
eta_3 ~ start(2)*1
eta_4 ~ start(.1)*1
# manifest variances for math (made equivalent by naming theta1)
math2 ~~ theta1*math2
math3 ~~ theta1*math3
math4 ~~ theta1*math4
math5 ~~ theta1*math5
math6 ~~ theta1*math6
math7 ~~ theta1*math7
math8 ~~ theta1*math8
# manifest variances for hyp (made equivalent by naming theta2)
hyp2 ~~ theta2*hyp2
hyp3 ~~ theta2*hyp3
hyp4 ~~ theta2*hyp4
hyp5 ~~ theta2*hyp5
hyp6 ~~ theta2*hyp6
hyp7 ~~ theta2*hyp7
hyp8 ~~ theta2*hyp8
# residual covariances (made equivalent by naming theta12)
math2 ~~ theta12*hyp2
math3 ~~ theta12*hyp3
math4 ~~ theta12*hyp4
math5 ~~ theta12*hyp5
math6 ~~ theta12*hyp6
math7 ~~ theta12*hyp7
math8 ~~ theta12*hyp8
# manifest means for math (fixed at zero)
math2 ~ 0*1
math3 ~ 0*1
math4 ~ 0*1
math5 ~ 0*1
math6 ~ 0*1
math7 ~ 0*1
math8 ~ 0*1
# manifest means for hyp (fixed at zero)
hyp2 ~ 0*1
hyp3 ~ 0*1
hyp4 ~ 0*1
hyp5 ~ 0*1
hyp6 ~ 0*1
hyp7 ~ 0*1
hyp8 ~ 0*1
' #end of model definitionAnd estimating the model …
#estimating the model using sem() function
bivariate_lavaan_fit <- sem(bivariate_lavaan_model,
data = nlsy_multi_wide,
meanstructure = TRUE,
estimator = "ML",
missing = "fiml")## Warning: lavaan->lav_data_full():
## some cases are empty and will be ignored: 741.## Warning: lavaan->lav_data_full():
## due to missing values, some pairwise combinations have less than 10%
## coverage; use lavInspect(fit, "coverage") to investigate.## Warning: lavaan->lav_mvnorm_missing_h1_estimate_moments():
## Maximum number of iterations reached when computing the sample moments
## using EM; use the em.h1.iter.max= argument to increase the number of
## iterationsWe can look more specifically at the numerical summaries of the model.
#obtaining summary of the model using the object we just created
summary(bivariate_lavaan_fit, fit.measures=TRUE)## lavaan 0.6-19 ended normally after 67 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 35
## Number of equality constraints 18
##
## Used Total
## Number of observations 932 933
## Number of missing patterns 96
##
## Model Test User Model:
##
## Test statistic 318.885
## Degrees of freedom 102
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 1532.041
## Degrees of freedom 91
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.849
## Tucker-Lewis Index (TLI) 0.866
##
## Robust Comparative Fit Index (CFI) 1.000
## Robust Tucker-Lewis Index (TLI) -6.557
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -11790.476
## Loglikelihood unrestricted model (H1) -11631.033
##
## Akaike (AIC) 23614.951
## Bayesian (BIC) 23697.186
## Sample-size adjusted Bayesian (SABIC) 23643.196
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.048
## 90 Percent confidence interval - lower 0.042
## 90 Percent confidence interval - upper 0.054
## P-value H_0: RMSEA <= 0.050 0.724
## P-value H_0: RMSEA >= 0.080 0.000
##
## Robust RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: Robust RMSEA <= 0.050 1.000
## P-value H_0: Robust RMSEA >= 0.080 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.097
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## eta_1 =~
## math2 1.000
## math3 1.000
## math4 1.000
## math5 1.000
## math6 1.000
## math7 1.000
## math8 1.000
## eta_2 =~
## math2 0.000
## math3 1.000
## math4 2.000
## math5 3.000
## math6 4.000
## math7 5.000
## math8 6.000
## eta_3 =~
## hyp2 1.000
## hyp3 1.000
## hyp4 1.000
## hyp5 1.000
## hyp6 1.000
## hyp7 1.000
## hyp8 1.000
## eta_4 =~
## hyp2 0.000
## hyp3 1.000
## hyp4 2.000
## hyp5 3.000
## hyp6 4.000
## hyp7 5.000
## hyp8 6.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## eta_1 ~~
## eta_2 -0.204 1.154 -0.176 0.860
## eta_3 -2.979 0.673 -4.426 0.000
## eta_4 0.107 0.164 0.654 0.513
## eta_2 ~~
## eta_3 0.098 0.161 0.608 0.543
## eta_4 -0.040 0.038 -1.061 0.289
## eta_3 ~~
## eta_4 -0.020 0.031 -0.644 0.520
## .math2 ~~
## .hyp2 (th12) -0.011 0.233 -0.046 0.963
## .math3 ~~
## .hyp3 (th12) -0.011 0.233 -0.046 0.963
## .math4 ~~
## .hyp4 (th12) -0.011 0.233 -0.046 0.963
## .math5 ~~
## .hyp5 (th12) -0.011 0.233 -0.046 0.963
## .math6 ~~
## .hyp6 (th12) -0.011 0.233 -0.046 0.963
## .math7 ~~
## .hyp7 (th12) -0.011 0.233 -0.046 0.963
## .math8 ~~
## .hyp8 (th12) -0.011 0.233 -0.046 0.963
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## eta_1 35.259 0.356 99.178 0.000
## eta_2 4.343 0.088 49.139 0.000
## eta_3 1.903 0.058 32.702 0.000
## eta_4 -0.057 0.014 -3.950 0.000
## .math2 0.000
## .math3 0.000
## .math4 0.000
## .math5 0.000
## .math6 0.000
## .math7 0.000
## .math8 0.000
## .hyp2 0.000
## .hyp3 0.000
## .hyp4 0.000
## .hyp5 0.000
## .hyp6 0.000
## .hyp7 0.000
## .hyp8 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## eta_1 64.711 5.664 11.424 0.000
## eta_2 0.752 0.329 2.282 0.023
## eta_3 1.542 0.156 9.908 0.000
## eta_4 0.005 0.009 0.539 0.590
## .math2 (tht1) 36.126 1.863 19.390 0.000
## .math3 (tht1) 36.126 1.863 19.390 0.000
## .math4 (tht1) 36.126 1.863 19.390 0.000
## .math5 (tht1) 36.126 1.863 19.390 0.000
## .math6 (tht1) 36.126 1.863 19.390 0.000
## .math7 (tht1) 36.126 1.863 19.390 0.000
## .math8 (tht1) 36.126 1.863 19.390 0.000
## .hyp2 (tht2) 1.104 0.057 19.404 0.000
## .hyp3 (tht2) 1.104 0.057 19.404 0.000
## .hyp4 (tht2) 1.104 0.057 19.404 0.000
## .hyp5 (tht2) 1.104 0.057 19.404 0.000
## .hyp6 (tht2) 1.104 0.057 19.404 0.000
## .hyp7 (tht2) 1.104 0.057 19.404 0.000
## .hyp8 (tht2) 1.104 0.057 19.404 0.000We get information about how the model fits the data, and the model parameters.
The main parameter of interest, the slope-slope covariance, is not significantly different than zero.
#from the output above
# Estimate Std.Err z-value P(>|z|)
#eta_3 ~~
# eta_4 -0.020 0.031 -0.644 0.520We can make a diagram to check if we specified and estimated the model as intended.
Thats cool!
Conclusion
This tutorial has presented how the multivariate growth model can be used to examine interrelations in change - specifically the between-person correlation of slopes.
Multivariate change is awesome, but be cautious with the presentation and inferences!
Citations
Bürkner, P.-C. (2017). brms: An R Package for Bayesian Multilevel Models Using Stan. Journal of Statistical Software, 80, 1–28. https://doi.org/10.18637/jss.v080.i01
Epskamp, S. (2022). semPlot: Path Diagrams and Visual Analysis of Various SEM Packages’ Output (Version 1.1.6). https://CRAN.R-project.org/package=semPlot
Grimm, K. J., Ram, N., & Estabrook, R. (2017). Growth Modeling: Structural Equation and Multilevel Modeling Approaches. Guilford Publications.
Pinheiro, J. C., & Bates, D. (2000). Mixed-Effects Models in S and S-PLUS. Springer-Verlag. https://doi.org/10.1007/b98882
Pinheiro, J. C., Bates, D., & R Core Team. (2025). nlme: Linear and Nonlinear Mixed Effects Models (Version 3.1-167). https://doi.org/10.1007/b98882
R Core Team. (2024). R: A Language and Environment for Statistical Computing. Foundation for Statistical Computing. https://www.R-project.org/
Revelle, W. (2024). psych: Procedures for Psychological, Psychometric, and Personality Research. Northwestern University. https://CRAN.R-project.org/package=psych
Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48, 1–36. https://doi.org/10.18637/jss.v048.i02
Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., … Yutani, H. (2019). Welcome to the Tidyverse. Journal of Open Source Software, 4(43), 1686. https://doi.org/10.21105/joss.01686