Section 13 Latent Curve Analysis (or Latent Growth Models)

13.1 Design Issues/Considerations

• Longitudinal research involve situations in which we make repeated measurements on same individual (or groups or families or…)
• Three general versions:
  1. Completely Randomized Design (CRD)
     • Goal is to eliminate ordering effects while increasing power/control
     • Experimental protocols/features
     • Consider the example of message effects
  2. Growth Curve Models
     • Change in a trait or variable measured at different (usually regularly-spaced) times. Examples include:
       • Changes in attitude
       • Performance (e.g., GPA)
       • Relationship development
       • Perceptions of competence/status across group meetings/tasks
     • In general, there is reason to believe that time is related to some sort of underlying developmental mechanism – otherwise why bother?
  3. Treatment X Time Designs
     • Development/time an important feature yet intervention is involved. Examples include:
       • Linguistic stimuli (see Baayen et al.)
       • Establish baseline before intervention/modification
       • Or one of my favorite examlples
• Note that analysis depends somewhat on whether the design includes fixed and/or varying occasions
  • Fixed: all individuals provide measurements at the same set of occasions, usually regularly spaced (e.g., every year)
  • Varying: different set of measures taken at different points in time for different individuals (e.g., birthdays at different parts of the year)
• We will cover the last two design/models, though the concepts apply in some degree to the first model/design

13.2 What We’ll Cover

1. General MLM applied to longitudinal data
   • We’ll spend the most time on growth curves (with covariates)
   • We will look at a treatment X time model
2. Testing different random covariance structures
   • We have already seen some of these, for example unstructured and compound symmetry, but there are more options to consider with longitudinal data with more than two measurement points
   • We have already done similar types of model tests, but keeping track of the number of parameters is the important issue
3. Piecemeal approach to nonlinear data
   • Nonlinear data can occur “naturally” or it could be the result of an intervention
   • The point of the piecemeal approach is to divide the data and analyze the different linear trends separately
4. LCM approaches to nested data

13.3 Curran et al: 12 Questions About Growth Curve Modeling

1. What is Growth Curve Modeling?
   • Methods that allow for the estimation of inter-individual variability in intra-individual patterns of change over time
   • Estimate between-person differences in within-person change
2. How Do Growth Models Differ From More Traditional Longitudinal Models?
   • More flexible, including
     • partially missing data
     • unequally spaced time points
     • non-normal distributions
     • complex trajectories
     • time-varying covariates
   • More powerful statistical tests
3. How Are Growth Models Fit to Data?
   • Multilevel modeling (allows nesting)
   • SEM uses latent variables and allows for more comprehensive and customized hypothesis testing
   • Mplus and lavvan do allow for latent variable modeling in nested designs
4. What Are the Data Requirements to Use Growth Models?
   • Large sample size
     • Depends on design factors (number of measurement times, number of items, intervention or observational)
   • Need at least three measurement times (otherwise you don’t have a pattern or trajectory)
   • ML estimation for normally distributed data (other options for non-normal data)
5. Can Growth Models be Estimated with Partially Missing Data?
   • ML gives more weight to cases with more data points than those with fewer
   • Imputation methods (assumes MCAR or CAR)
6. What Are the Different Shapes of Growth Curves That Can Be Modeled?
   • Often start with “random intercept-only” model
     • Stable level of RM over time and individuals vary randomly around that level over time
   • Then more complex models, including linear, quadratic, cubic, and so on
7. How Is the Adequacy of Fit for Growth Models Assessed?
   • In SEM it’s MOSTLY the usual case of \(\chi^{2}\) and other related model fit indices
     • Remember, the fit is compared to baseline models (either best or worst fitting model); absolute and relative fit
   • In MLM model comparison is used, with significant \(\chi^{2}\) statistic indicating better fit
8. How Can Predictors Be Incorporated Into the Growth Model?
   • Can be added once trajectory has been established
   • Models with predictors are called “conditional”
   • Time-invariate covariates (TICs) and time-varying covariates (TVCs)
9. Can Growth in Two Constructs Be Measured Over Time
   • Yup: multivariate growth models
   • Curves might be different for each LV–pretty cool!
10. Can Growth Models be Simultaneously Estimated Within Two or More Groups?
    • One way is to use group as a TIC
      • Shifts conditional means of intercept and slope higher or lower
    • Another way is invariance tests of the types described previously
11. What If There Is a Potentially Important Grouping Variable That Was Not Directly Observed?
    • Can evaluate because theory might latent groups in the data (e.g., developmental)
    • Mixture models come to the rescue here
12. Where Do I Go From Here?
    • Lots of great sources
    • I started with Singer (1998) and Singer and Willet (2003)

13.4 Why Not Repeated Measures?

• RM analyzes only mean differences, not linear trends (at least not in the same way)
  • Time is treated as a categorical variable
  • In LCM it is continuous (though can be evaluated categorically as needed)
• RM treats differences among individuals (e.g., lowest level) as error variance
  • In growth models, the residual is within individuals (i.e., across measurement times)
• RM within an MLM context allows for only a single response, not anything resembling a latent variable approach
  • Growth models within an SEM can include measurement models
• Finally, RM makes assumptions about the covariance matrix
  • Univariate RM assumes compound symmetry
  • Multivariate RM assumes unstructured matrix

13.5 Preliminaries

• A latent curve model examines change over time
  • Also called growth models (though not every variable of interest “grows” over time)
• An extension of SEM, the latent parts are as follows:
  • The intercept is measured by item means at each measurement time
  • The slope is measured by the change in means at each time period
• The loadings for the manifest variables on the latent intercept are usually fixed to 1, indicating that the LV intercept has equal influence on the indicators at all time periods
  • Need good reason to assume otherwise
• The loadings on the slope generally indicate some increment of change, and often have some relationship to the actual time between measurements
  • The generic formula is to have T1 = 0, T2 =1, T3 = 2, and so on.
  • Remember that 0 is an important value in regression!
  • The “loadings” may change if the intervals between measurements are not equal or have some other important characteristic
• LCM accounts for both mean of the MVs as well as changes in them.
  • In general, MV intercepts are constrained to 0
    • MV means modeled entirely by LV means
    • This is the lavaan default! (And probably true of other programs, too.)
  • Residual variances and covariances, as well as latent variances and covariances are modeled, as usual
• LCM graph for 4 measurement times

Figure 13.1: Basic LCM, 4 Measurement Times

13.6 General Questions When Creating an LCM

1. What is the average trajectory for all respondents?
   • What is the initial level (mean) of the variable?
   • Does the variable change over time (i.e., is the slope “flat”)?
   • What is the form of the change?
2. Is the average trajectory sufficient, or should one allow for variance in trajectories?
   • Do the latent intercept and slope require a variance term or can they be constrained to 0?
   • In MLM analyses these are often referred to as “random effects”
3. If distinct trajectories (i.e., random effects) are needed, what other variables explain them?
   • Time-invariant covariates (i.e., value does not change across time, e.g., gender)
   • Time-varying covariates (i.e., value does change over time, e.g., age)

13.7 Data Description and Setup

• Here we will work with the GPA data described in Hox
• The data are measured over six semesters (i.e., fixed occasions). Other variables include:
  • Student gender, high school GPA, and hours per week working part time
• Notice that some variables have different values each semester (i.e., vary within person) and others don’t (vary between persons) – this will matter to us in the analyses that follow
• The first 12 rows (i.e., first two children–see the student column) of the data file in long format follow
  • Long format useful for MLM regression (see Table 13.1)

read.table("https://commresearch.arizona.edu/classes/comm640/data/Hox/rawdata/gpaflat.dat", 
           header = T) %>% 
  janitor::clean_names() %>% 
  {. ->> gpa0} %>% 
  slice(1:12) %>% 
  gt(caption="Hox GPA Data, Long Format") 

Table 13.1: Hox GPA Data, Long Format

student	occas	gpa	job	sex	highgpa	admitted
1	1	2.3	2	2	2.8	1
1	2	2.1	1	2	2.8	1
1	3	3.0	1	2	2.8	1
1	4	3.0	4	2	2.8	1
1	5	3.0	3	2	2.8	1
1	6	3.3	2	2	2.8	1
2	1	2.2	2	1	2.5	0
2	2	2.5	3	1	2.5	0
2	3	2.6	2	1	2.5	0
2	4	2.6	0	1	2.5	0
2	5	3.0	2	1	2.5	0
2	6	2.8	2	1	2.5	0

• We’ll need to construct a wide file, too
  • Wide file useful for SEM–see Table 13.2
  • Notice I include jobs in the wide file
    • It is a time-varying covariate

gpa0 %>% 
  janitor::clean_names() %>% 
  pivot_wider(id_cols = c(student, sex, 
                          highgpa, admitted), 
              names_from=occas, 
              values_from = c(job,gpa),
              names_sep = "") %>% 
  {. ->> gpa.w} %>% #save new file and pipe
  slice(1:10) %>% 
  gt(caption = "Hox GPA Data, Wide Format") 

Table 13.2: Hox GPA Data, Wide Format

student	sex	highgpa	admitted	job1	job2	job3	job4	job5	job6	gpa1	gpa2	gpa3	gpa4	gpa5	gpa6
1	2	2.8	1	2	1	1	4	3	2	2.3	2.1	3.0	3.0	3.0	3.3
2	1	2.5	0	2	3	2	0	2	2	2.2	2.5	2.6	2.6	3.0	2.8
3	2	2.5	1	2	2	2	3	2	1	2.4	2.9	3.0	2.8	3.3	3.4
4	1	3.8	0	3	2	2	1	2	2	2.5	2.7	2.4	2.7	2.9	2.7
5	1	3.1	1	1	2	2	1	2	2	2.8	2.8	2.8	3.0	2.9	3.1
6	2	2.9	0	2	3	4	1	3	3	2.5	2.4	2.4	2.3	2.7	2.8
7	1	2.3	1	3	1	3	1	3	3	2.4	2.4	2.8	2.6	3.0	3.0
8	2	3.9	1	2	1	3	2	2	1	2.8	2.8	3.1	3.3	3.3	3.4
9	1	2.0	9	1	1	3	3	3	1	2.8	2.7	2.7	3.1	3.1	3.5
10	1	2.8	1	1	1	1	3	1	1	2.8	2.8	3.0	2.7	3.0	3.0

13.7.1 Descriptive Stats

• The time-based correlations can be revealing
• Provides a sense of whether and how measurements are related across time
  • For example, gpa more highly correlated across the last 3 measurement times than in the first 3

gpa.w %>% 
  select(gpa1:gpa6) %>% 
  mlmCorrs::corstars(title = "Correlations for GPA Data")

Correlations for GPA Data
Variable	Mean	SD	1	2	3	4	5	6
1. Gpa1	2.59	0.31	--
2. Gpa2	2.72	0.34	.34**	--
3. Gpa3	2.81	0.35	.24**	.43**	--
4. Gpa4	2.92	0.36	.19**	.47**	.66**	--
5. Gpa5	3.02	0.36	.22**	.51**	.71**	.84**	--
6. Gpa6	3.13	0.38	.22**	.43**	.67**	.83**	.88**	--
Note. *p < .05. **p < .01.

13.7.2 Let’s Plot the Data

• We can plot growth curves for each child in the data set.
• Raw data in long format, which is good
• See http://www.ats.ucla.edu/stat/r/faq/longitudinal.htm
• Look at the intercepts and trajectories for each
  • It is plausible that a fixed intercept works for all?
  • Does a single slope describe the trajectory for all children?

#each institution is identified by student number, and each has 4 lines per data (assuming no missing values)
ggplot(data = gpa0, aes(x = as.factor(occas), y = gpa, group=student)) +
  geom_line() +
  labs(x="Time", y="GPA") +
  ggtitle("GPA Plots")

Figure 13.2: Plot of GPA Trajectories

• Here’s a variation on the theme–each child has it’s own plot and slope
• Note I produce plots for the first 12 children in order to save space

#here's the panel or lattice plots.  Pretty cool!
#note that this routine plots regression line and observed values
#First 12 students to save space
gpa0 %>% 
  slice(1:72) %>% 
ggplot(aes(x = occas, y = gpa)) + 
  geom_point() +
  stat_smooth(method = "lm", se = FALSE) +
  labs(x="Time", y="GPA") +
  facet_wrap(~student) +
  ggtitle("GPA Trends by Student")

Figure 13.3: Lattice Plots for GPA Data

13.8 Analysis Strategy: Model Building

• Note that this is a model building process
  • Here we look for improvement in model fit because we will add parameters to the model
    • Recall from previous chapters that adding a parameter always improves fit
    • The question is does adding a parameter improve fit enough?
  • We want \(\chi^{2}\) to be significant
  • Compare with invariance tests, where we want \(\chi^{2}\) not to be significant
    • Invariance tests decrease the number or estimated parameters
    • Reducing the number of parameters always worsens model fit

13.8.1 Baseline Model

• Mean Latent Intercept and Constrained Residual Variances
  • Only the mean across all time periods is modeled (and residual variance)
  • If this were the best fitting model, it would indicate that all children have more or less the same GPA that does not change over time
• lavaan’s growth routine has an interesting set of defaults
  • Assumes means and intercepts
  • Observed intercepts fixed to 0 (we care only about the latent intercept)
  • Latent means freely estimated
• Note I include slope estimates, all of which are fixed to zero
  • We will need this for model comparisons to follow and it’s a simple matter to copy the code and remove constraints as needed
• Degrees of freedom based number of observed values - estimated parameters
  • Observed valus: 15 covariances + 6 variances + 6 means/intercepts = 27
  • Estimated parameters: 1 intercept + 1 residual variance = 2
  • df = 27-2 = 25

gpa.model1 <- '

# intercept
i =~ 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4 + 1*gpa5 + 1*gpa6
#intercept variance fixed to 0
i~~0*i

# slope
#slope set to zero for model comparison test
s =~ 0*gpa1 + 0*gpa2 + 0*gpa3 + 0*gpa4 + 0*gpa5 + 0*gpa6

#slope intercept set to 0
s ~ 0*1

#set slope variance to 0
s ~~ 0*s

#set slope/intercept covariance to 0
s~~0*i

# residual variances
gpa1~~res*gpa1
gpa2~~res*gpa2
gpa3~~res*gpa3
gpa4~~res*gpa4
gpa5~~res*gpa5
gpa6~~res*gpa6
'

• Here’s a graph of the model with the weights on the edges (arrows)

Figure 13.4: Baseline LCM

• Let’s look at the output using the summary options

gpa.fit1 <- growth(gpa.model1, data=gpa.w)
summary(gpa.fit1)

lavaan 0.6-11 ended normally after 16 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         7
  Number of equality constraints                     5
                                                      
  Number of observations                           200
                                                      
Model Test User Model:
                                                      
  Test statistic                              1108.453
  Degrees of freedom                                25
  P-value (Chi-square)                           0.000

Parameter Estimates:

  Standard errors                             Standard
  Information                                 Expected
  Information saturated (h1) model          Structured

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)
  i =~                                                
    gpa1              1.000                           
    gpa2              1.000                           
    gpa3              1.000                           
    gpa4              1.000                           
    gpa5              1.000                           
    gpa6              1.000                           
  s =~                                                
    gpa1              0.000                           
    gpa2              0.000                           
    gpa3              0.000                           
    gpa4              0.000                           
    gpa5              0.000                           
    gpa6              0.000                           

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)
  i ~~                                                
    s                 0.000                           

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
    s                 0.000                           
   .gpa1              0.000                           
   .gpa2              0.000                           
   .gpa3              0.000                           
   .gpa4              0.000                           
   .gpa5              0.000                           
   .gpa6              0.000                           
    i                 2.865    0.011  252.610    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
    i                 0.000                           
    s                 0.000                           
   .gpa1     (res)    0.154    0.006   24.495    0.000
   .gpa2     (res)    0.154    0.006   24.495    0.000
   .gpa3     (res)    0.154    0.006   24.495    0.000
   .gpa4     (res)    0.154    0.006   24.495    0.000
   .gpa5     (res)    0.154    0.006   24.495    0.000
   .gpa6     (res)    0.154    0.006   24.495    0.000

• Prettier estimates in Table 13.3
• I didn’t print the loadings as they are fixed and not really needed
• Note the poor model fit
• This gives us a comparison baseline
• I left the loadings in for your inspection but will exclude them in models that follow
• I don’t reveal code for subsequent table to save space, but it is mostly the same to what you see below

semTable::semTable(gpa.fit1,
                   paramSets = c("residualvariances",
                                           "latentvariances",
                                           "latentmeans",
                                           "latentcovariances"),
                   type="html", columns = c("estsestars"),
                   print.results = F) %>%
  rvest::minimal_html() %>% 
  rvest::html_table(., fill=T, header = F) %>% 
  as.data.frame() %>% 
  slice(2:n()) %>% #gets rid of X? variable names
  gt::gt(caption = ("Growth Models for GPA Data**")) %>% 
  fmt_missing(columns = 1:2, missing_text="") %>% 
  cols_label(X1=md("**Term**"), X2=md("**Fixed Intercept Model**")) %>% 
  tab_options(table_body.hlines.width   = 0,
                table.border.bottom.width = 0,
)

Table 13.3: Growth Models for GPA Data**

Term	Fixed Intercept Model
	Estimate(Std.Err.)
	Residual Variances
gpa1	0.15(0.01)***
gpa2	0.15(0.01)***
gpa3	0.15(0.01)***
gpa4	0.15(0.01)***
gpa5	0.15(0.01)***
gpa6	0.15(0.01)***
	Latent Intercepts
s	0.00+
i	2.87(0.01)***
	Latent Variances
i	0.00+
s	0.00+
	Latent Covariances
i w/s	0.00+
	Fit Indices
χ2	1108.45(25)***
CFI	0.00
TLI	0.18
RMSEA	0.47
+Fixed parameter	+Fixed parameter
*p<0.05, **p<0.01, ***p<0.001	*p<0.05, **p<0.01, ***p<0.001

• Where do the estimates come from?
  • Remember that long file? Look at the mean and variance for gpa (across all 6 times)
  • Look familiar–intercept value and residual variance for this first model

mean(gpa0$gpa)

[1] 2.86

var(gpa0$gpa)

[1] 0.154

• Notice the residual variance is 0.154 (which is the variance of the MVs, all constrained to be equal)

13.8.2 Unconditional Model

• Mean Latent Intercept That Is Allowed to Vary, and Constrained Residual Variances
• This model is akin to the unconditional model in MLM parlance
  • Unconditional means that the estimates are not conditioned on other factors in the model
  • Grand mean collapsed across time, though between-person variance is modeled
• If the intercept is allowed to vary and the model fit is significantly improved
  • The average intercept is not sufficient
  • Some children have greater and smaller “starting points” (T1) than average
  • Notice the difference in the code–the constraint on the intercept variance as been removed
• Importantly, we can ask, if variances of intercepts is significant, what might “cause” the variation
  • We’ll address this issue in subsequent models
• One last thing: notice the new residual estimate is quite a bit smaller
  • The intercept now has a variance
  • The intercept variance + residual = total variance in baseline model
• df now = 24, as one additional parameter, intercept variance, is estimated

gpa.model2 <- '
# intercept
i =~ 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4 + 1*gpa5 + 1*gpa6

#intercept variance NOT fixed to 0
#i~~0*i

# slope
#slope set to zero for model comparison test
s =~ 0*gpa1 + 0*gpa2 + 0*gpa3 + 0*gpa4 + 0*gpa5 + 0*gpa6

#slope intercept set to 0
s ~ 0*1

#set slope variance to 0
s ~~ 0*s

#set slope/intercept covariance to 0
s~~0*i

# residual variances
gpa1~~res*gpa1
gpa2~~res*gpa2
gpa3~~res*gpa3
gpa4~~res*gpa4
gpa5~~res*gpa5
gpa6~~res*gpa6
'

Figure 13.5: Unconditional LCM

Table 13.4: Growth Models for GPA Data

Term	Fixed Intercept Model	Random Intercept Model
	Estimate(Std.Err.)	Estimate(Std.Err.)
	Residual Variances	Residual Variances
gpa1	0.15(0.01)***	0.10(0.00)***
gpa2	0.15(0.01)***	0.10(0.00)***
gpa3	0.15(0.01)***	0.10(0.00)***
gpa4	0.15(0.01)***	0.10(0.00)***
gpa5	0.15(0.01)***	0.10(0.00)***
gpa6	0.15(0.01)***	0.10(0.00)***
	Latent Intercepts	Latent Intercepts
s	0.00+	0.00+
i	2.87(0.01)***	2.87(0.02)***
	Latent Variances	Latent Variances
i	0.00+	0.06(0.01)***
s	0.00+	0.00+
	Latent Covariances	Latent Covariances
i w/s	0.00+	0.00+
	Fit Indices	Fit Indices
χ2	1108.45(25)***	858.63(24)***
CFI	0.00	0.00
TLI	0.18	0.35
RMSEA	0.47	0.42
+Fixed parameter	+Fixed parameter	+Fixed parameter
*p<0.05, **p<0.01, ***p<0.001	*p<0.05, **p<0.01, ***p<0.001	*p<0.05, **p<0.01, ***p<0.001

• Compare models
• Allowing the intercept to vary across participants substantially improves model fit

#test the models
anova(gpa.fit1, gpa.fit2)

Chi-Squared Difference Test

         Df  AIC  BIC Chisq Chisq diff Df diff          Pr(>Chisq)    
gpa.fit2 24  919  929   859                                           
gpa.fit1 25 1167 1174  1108        250       1 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Notice that the model \(\chi^{2}\) estimates are significantly different, which is GOOD in this case.
• Residual variance = 0.098, which is much lower than the estimate in model 1: 0.154

13.8.3 Fixed Slope Model

• Mean Latent Intercept That Is Allowed to Vary, Estimated Slope (Not Allowed to Vary), and Constrained Residual Variances
• We’ve not estimated the slope until now, and doing so influences the estimates
  • Notice the intercept value has changed
  • Mean at T1 is 2.594
• Slope variance is NOT estimated (which means the slope is assumed to be sufficient for all students)
• Slope and intercept covariance set to zero because we have no slope variance
  • Is akin to the MLM model with random intercept variance but no random slope variance
• How do we interpret the slope?
• Why did the intercept value change?
  • Compare to Table ??

gpa.model3 <- '
# intercept
i =~ 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4 + 1*gpa5 + 1*gpa6

#intercept variance NOT fixed to 0
i~~i

# slope
#slope set to zero for model comparison test
s =~ 0*gpa1 + 1*gpa2 + 2*gpa3 + 3*gpa4 + 4*gpa5 + 5*gpa6

#slope intercept NOT set to 0
s ~ 1

#set slope variance to 0
s ~~ 0*s

#set slope/intercept covariance to 0
s~~0*i

# residual variances
gpa1~~res*gpa1
gpa2~~res*gpa2
gpa3~~res*gpa3
gpa4~~res*gpa4
gpa5~~res*gpa5
gpa6~~res*gpa6
'

Figure 13.6: Fixed Slope LCM

Table 13.5: Growth Models for GPA Data

Term	Fixed Intercept Model	Random Intercept Model	Fixed Slope Model
	Estimate(Std.Err.)	Estimate(Std.Err.)	Estimate(Std.Err.)
	Residual Variances	Residual Variances	Residual Variances
gpa1	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***
gpa2	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***
gpa3	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***
gpa4	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***
gpa5	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***
gpa6	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***
	Latent Intercepts	Latent Intercepts	Latent Intercepts
s	0.00+	0.00+	0.11(0.00)***
i	2.87(0.01)***	2.87(0.02)***	2.60(0.02)***
	Latent Variances	Latent Variances	Latent Variances
i	0.00+	0.06(0.01)***	0.06(0.01)***
s	0.00+	0.00+	0.00+
	Latent Covariances	Latent Covariances	Latent Covariances
i w/s	0.00+	0.00+	0.00+
	Fit Indices	Fit Indices	Fit Indices
χ2	1108.45(25)***	858.63(24)***	338.82(23)***
CFI	0.00	0.00	0.60
TLI	0.18	0.35	0.74
RMSEA	0.47	0.42	0.26
+Fixed parameter	+Fixed parameter	+Fixed parameter	+Fixed parameter
*p<0.05, **p<0.01, ***p<0.001	*p<0.05, **p<0.01, ***p<0.001	*p<0.05, **p<0.01, ***p<0.001	*p<0.05, **p<0.01, ***p<0.001

• Compare models
• Adding a fixed slope improves model fit

anova(gpa.fit3, gpa.fit2)

Chi-Squared Difference Test

         Df AIC BIC Chisq Chisq diff Df diff          Pr(>Chisq)    
gpa.fit3 23 402 415   339                                           
gpa.fit2 24 919 929   859        520       1 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

13.8.4 Random Slope Model

• Mean Latent Intercept That Is Allowed to Vary, Mean Latent Slope That Is Estimated and Allowed to Vary, and Constrained Residual Variances
• Slope is estimated
• Variance of slope is estimated (akin to the random slope model in MLM)
• Covariance between slope and intercept estimated
  • If significant, the starting point (score at T1) has some bearing on the rate of change

gpa.model4 <- '
# intercept
i =~ 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4 + 1*gpa5 + 1*gpa6

#intercept variance NOT fixed to 0
i~~i

# slope
#slope set to zero for model comparison test
s =~ 0*gpa1 + 1*gpa2 + 2*gpa3 + 3*gpa4 + 4*gpa5 + 5*gpa6

#slope intercept NOT set to 0
s ~ 1

#set slope variance NOT set to 0
s ~~ s

#NOT set slope/intercept covariance to 0
s~~i

# residual variances
gpa1~~res*gpa1
gpa2~~res*gpa2
gpa3~~res*gpa3
gpa4~~res*gpa4
gpa5~~res*gpa5
gpa6~~res*gpa6
'

Figure 13.7: Random Slope LCM

	Fixed Intercept	Random Intercept	Fixed Slope	Random Slope
	Estimate(Std.Err.)	Estimate(Std.Err.)	Estimate(Std.Err.)	Estimate(Std.Err.)
	Residual Variances
gpa1	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***
gpa2	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***
gpa3	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***
gpa4	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***
gpa5	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***
gpa6	0.15(0.01)***	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***
	Latent Intercepts
s	0.00+	0.00+	0.11(0.00)***	0.11(0.01)***
i	2.87(0.01)***	2.87(0.02)***	2.60(0.02)***	2.60(0.02)***
	Latent Variances
i	0.00+	0.06(0.01)***	0.06(0.01)***	0.04(0.01)***
s	0.00+	0.00+	0.00+	0.00(0.00)***
	Latent Covariances
i w/s	0.00+	0.00+	0.00+	-0.00(0.00)
	Fit Indices
χ2	1108.45(25)***	858.63(24)***	338.82(23)***	191.41(21)***
CFI	0.00	0.00	0.60	0.79
TLI	0.18	0.35	0.74	0.85
RMSEA	0.47	0.42	0.26	0.20
+Fixed parameter
*p<0.05, **p<0.01, ***p<0.001

• Model comparison indicates that adding a random slope improves model fit
• The interpretation is that the slope varies across children

# model fit improved
anova(gpa.fit3, gpa.fit4)

Chi-Squared Difference Test

         Df AIC BIC Chisq Chisq diff Df diff          Pr(>Chisq)    
gpa.fit4 21 258 278   191                                           
gpa.fit3 23 402 415   339        147       2 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

13.8.5 Unconstrained Model

• Same as previous model except that item residual variances estimated rather than fixed to be equal

gpa.model5 <- '
# intercept
# intercept
i =~ 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4 + 1*gpa5 + 1*gpa6

#intercept variance NOT fixed to 0
i~~i

# slope
#slope set to zero for model comparison test
s =~ 0*gpa1 + 1*gpa2 + 2*gpa3 + 3*gpa4 + 4*gpa5 + 5*gpa6

#slope intercept NOT set to 0
s ~ 1

#set slope variance NOT set to 0
s ~~ s

#NOT set slope/intercept covariance to 0
s~~i

# residual variances allowed to freely vary
gpa1~~gpa1
gpa2~~gpa2
gpa3~~gpa3
gpa4~~gpa4
gpa5~~gpa5
gpa6~~gpa6
'

Figure 13.8: Unconstrained LCM

Table 13.6: Growth Models for GPA Data

Term	Random Intercept Model	Fixed Slope Model	Random Slope Model	Residual Variances
	Estimate(Std.Err.)	Estimate(Std.Err.)	Estimate(Std.Err.)	Estimate(Std.Err.)
	Residual Variances	Residual Variances	Residual Variances	Residual Variances
gpa1	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***	0.08(0.01)***
gpa2	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***	0.07(0.01)***
gpa3	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***	0.05(0.01)***
gpa4	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***	0.03(0.00)***
gpa5	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***	0.01(0.00)***
gpa6	0.10(0.00)***	0.06(0.00)***	0.04(0.00)***	0.02(0.00)***
	Latent Intercepts	Latent Intercepts	Latent Intercepts	Latent Intercepts
s	0.00+	0.11(0.00)***	0.11(0.01)***	0.11(0.01)***
i	2.87(0.02)***	2.60(0.02)***	2.60(0.02)***	2.60(0.02)***
	Latent Variances	Latent Variances	Latent Variances	Latent Variances
s	0.00+	0.00+	0.00(0.00)***	0.00(0.00)***
i	0.06(0.01)***	0.06(0.01)***	0.04(0.01)***	0.03(0.01)***
	Latent Covariances	Latent Covariances	Latent Covariances	Latent Covariances
i w/s	0.00+	0.00+	-0.00(0.00)	0.00(0.00)
	Fit Indices	Fit Indices	Fit Indices	Fit Indices
χ2	858.63(24)***	338.82(23)***	191.41(21)***	43.95(16)***
CFI	0.00	0.60	0.79	0.96
TLI	0.35	0.74	0.85	0.97
RMSEA	0.42	0.26	0.20	0.09
+Fixed parameter	+Fixed parameter	+Fixed parameter	+Fixed parameter	+Fixed parameter
*p<0.05, **p<0.01, ***p<0.001	*p<0.05, **p<0.01, ***p<0.001	*p<0.05, **p<0.01, ***p<0.001	*p<0.05, **p<0.01, ***p<0.001	*p<0.05, **p<0.01, ***p<0.001

• Link to MLM Table ??
• One last model fit
• Improvement, hence the random slope model seems the best fitting model

anova(gpa.fit5, gpa.fit4)

Chi-Squared Difference Test

         Df AIC BIC Chisq Chisq diff Df diff          Pr(>Chisq)    
gpa.fit5 16 121 157    44                                           
gpa.fit4 21 258 278   191        148       5 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

13.8.6 Add Gender as a Predictor

• Gender is a time invariant predictor (maale = 0)
  • Dummy coding interpretations same as always
    • LV intercept and slope estimates are for group coded 0
    • Regression estimates indicate difference of group coded 1 from the reference group
• Notice that both intercept and slope are regressed on gender
• Builds on the “random slope” model from above, which was the best fitting model
• Note that the structural model is different, hence:
  • Model comparisons with previous models not warranted because covariance matrices not the same

#recode sex
gpa.w$sex0 <- gpa.w$sex-1

gpa.model6 <- '
# intercept
# intercept
i =~ 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4 + 1*gpa5 + 1*gpa6

#intercept variance NOT fixed to 0
i~~i

# slope
#slope set to zero for model comparison test
s =~ 0*gpa1 + 1*gpa2 + 2*gpa3 + 3*gpa4 + 4*gpa5 + 5*gpa6

#slope intercept NOT set to 0
s ~ 1

#set slope variance NOT set to 0
s ~~ s

#NOT set slope/intercept covariance to 0
s~~i

#regress slope and intercept on sex0
i + s ~ sex0

# residual variances allowed to freely vary
gpa1~~gpa1
gpa2~~gpa2
gpa3~~gpa3
gpa4~~gpa4
gpa5~~gpa5
gpa6~~gpa6
'

Figure 13.9: Model with Gender as Predictor

• How do we interpret the effect of gender on the model intercepts and slopes?

Table 13.7: Intercept-Only Growth Model With Intercept Variance

Gender Differences Estimate(Std.Err.) Regression Slopes i sex0 0.09(0.04)* s sex0 0.03(0.01)* Residual Variances gpa1 0.08(0.01)*** gpa2 0.07(0.01)*** gpa3 0.05(0.01)*** gpa4 0.03(0.00)*** gpa5 0.01(0.00)*** gpa6 0.02(0.00)*** sex0 0.25+ Latent Intercepts s 0.09(0.01)*** i 2.55(0.03)*** Latent Variances i 0.03(0.01)*** s 0.00(0.00)*** Latent Covariances i w/s 0.00(0.00)  Fit Indices χ2 45.48(20)** CFI 0.97 TLI 0.97 RMSEA 0.08 +Fixed parameter p<0.05, p<0.01, p<0.001

13.8.7 Include Hours Worked as a Predictor

• Hours worked (“jobs”) is a time-varying covariate
• The dataset does not have a time-variant covariate but imagine school or athletic budget

gpa.model7 <- '

# intercept
i =~ 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4 + 1*gpa5 + 1*gpa6

#intercept variance NOT fixed to 0
i~~i

# slope
#slope set to zero for model comparison test
s =~ 0*gpa1 + 1*gpa2 + 2*gpa3 + 3*gpa4 + 4*gpa5 + 5*gpa6

#slope intercept NOT set to 0
s ~ 1

#set slope variance NOT set to 0
s ~~ s

#NOT set slope/intercept covariance to 0
s~~i

#regress slope and intercept on sex0
i + s ~ sex0

#regress jobs on gpa
gpa1 ~ job1
gpa2 ~ job2
gpa3 ~ job3
gpa4 ~ job4
gpa5 ~ job5
gpa6 ~ job6


# residual variances allowed to freely vary
gpa1~~gpa1
gpa2~~gpa2
gpa3~~gpa3
gpa4~~gpa4
gpa5~~gpa5
gpa6~~gpa6
'

Figure 7.1: Include Hours Worked as Predictor

Table 13.8: Intercept-Only Growth Model With Intercept Variance

Hours Worked Estimate(Std.Err.) Regression Slopes i sex0 0.09(0.04)* s sex0 0.03(0.01)* gpa1 job1 -0.01(0.02)  gpa2 job2 -0.01(0.02)  gpa3 job3 -0.01(0.01)  gpa4 job4 -0.01(0.01)  gpa5 job5 -0.01(0.01)  gpa6 job6 -0.01(0.01)  Residual Variances gpa1 0.08(0.01)*** gpa2 0.07(0.01)*** gpa3 0.05(0.01)*** gpa4 0.03(0.00)*** gpa5 0.01(0.00)*** gpa6 0.02(0.00)*** sex0 0.25+ job1 0.47+ job2 0.53+ job3 0.95+ job4 0.91+ job5 1.18+ job6 1.15+ Latent Intercepts s 0.09(0.01)*** i 2.57(0.04)*** Latent Variances i 0.03(0.01)*** s 0.00(0.00)*** Latent Covariances i w/s 0.00(0.00)  Fit Indices χ2 160.81(50)*** CFI 0.88 TLI 0.86 RMSEA 0.11 +Fixed parameter p<0.05, p<0.01, p<0.001

13.9 Multilevel Regression Estimates for LCM

• Need to use long file for MLM estimates
• Notice the progression in models
  • mlm2 below equals gpa.fit4 above
  • mlm3 equals the model (6a) above
• MLM identifies the variances for the intercepts and slopes as random
  • Random is an important concept, as most inferential models include at least one random component
  • MLM models often include several random components, which makes this an attractive option for more complex models
• Let’s model without nesting
• lme4 won’t work on non-nested models
• tab_model won’t print variances because they aren’t random
• Notice how the output lines up nicely with baseline model

head(gpa0)

  student occas gpa job sex highgpa admitted
1       1     1 2.3   2   2     2.8        1
2       1     2 2.1   1   2     2.8        1
3       1     3 3.0   1   2     2.8        1
4       1     4 3.0   4   2     2.8        1
5       1     5 3.0   3   2     2.8        1
6       1     6 3.3   2   2     2.8        1

mlm.lm <- lm(gpa ~ 1, data=gpa0)
summary(mlm.lm)

Call:
lm(formula = gpa ~ 1, data = gpa0)

Residuals:
   Min     1Q Median     3Q    Max 
-1.165 -0.265 -0.065  0.235  1.135 

Coefficients:
            Estimate Std. Error t value            Pr(>|t|)    
(Intercept)   2.8650     0.0113     252 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.393 on 1199 degrees of freedom

• Notice the unconditional MLM decomposes the total variance into the residual and between-individual
• Compare these estimates with those in Table 13.6

#occas0 is the zero-based time variable 
#unconditional model with intercept variance estimated
#data are nested within student (which is community)

gpa0$occas0 <- gpa0$occas-1
gpa0$sex0 <- gpa0$sex-1

mlm0 <- lmer(gpa ~ (1|student), data=gpa0, REML=FALSE) 

#include time fixed effect
mlm1 <- lmer(gpa ~ occas0 + (1|student) , data=gpa0, REML=FALSE)

#include time as a random effect (and produce unstructured matrix)
mlm2 <- lmer(gpa ~ occas0 + (1+occas0|student) , data=gpa0, REML=FALSE) 

#include sex as a predictor 
mlm3 <- lmer(gpa ~ occas0 + sex0 + (1+occas0|student) , data=gpa0, REML=FALSE) 

#awesome package for producing side-by-side models in one table


sjPlot::tab_model(list(mlm0, mlm1, mlm2, mlm3),
                  show.dev = TRUE, show.ci = F,
                   digits = 3, p.style = "stars",
                  dv.labels = c("Unconditional","Fixed Slope", 
                                "Random Slope", "Gender"))$knitr %>% 
  html()

Registered S3 method overwritten by 'parameters':
  method                         from      
  format.parameters_distribution datawizard

	Unconditional	Fixed Slope	Random Slope	Gender
Predictors	Estimates	Estimates	Estimates	Estimates
(Intercept)	2.865 ***	2.599 ***	2.599 ***	2.541 ***
occas0		0.106 ***	0.106 ***	0.106 ***
sex0				0.110 **
Random Effects
σ2	0.10	0.06	0.04	0.04
τ00	0.06 student	0.06 student	0.04 student	0.04 student
τ11			0.00 student.occas0	0.00 student.occas0
ρ01			-0.10 student	-0.15 student
ICC	0.37	0.52	0.65	0.64
N	200 student	200 student	200 student	200 student
Observations	1200	1200	1200	1200
Marginal R2 / Conditional R2	0.000 / 0.368	0.214 / 0.624	0.214 / 0.725	0.236 / 0.722
Deviance	913.456	393.649	246.234	236.068
* p<0.05   ** p<0.01   *** p<0.001

13.10 A Word About Covariance Structures

• We’ve seen several different random covariance structures to this point (e.g., compound symmetry)
• For growth models, one might be interested in the covariances of the measurement times
  • For example, one might assume that the association among all measurement points is approximately equal
  • As another example, one might assume that adjacent measurement times have higher associations than do distal ones
• The trick is to determine the best fitting covariance structure for the residual
  • Parsimony is the ideal, which is the most useful feature of the compound symmetry matrix
  • The best fitting model will always be the unstructured matrix
• The question is what is the pattern of variance and covariance over time?

13.10.1 Some Covariance Structures

• Compound symmetry is the simplest residual structure discussed here
• Note that variance at each time are equal and all covariances are constrained to be equal

\[ \begin{bmatrix} \sigma^{2}_{e}+\tau^{2}_{U_{0}}& \tau^{2}_{U_{0}} &\tau^{2}_{U_{0}} &\tau^{2}_{U_{0}} \\ \tau^{2}_{U_{0}}&\sigma^{2}_{e}+\tau^{2}_{U_{0}} &\tau^{2}_{U_{0}} &\tau^{2}_{U_{0}} \\ \tau^{2}_{U_{0}}& \tau^{2}_{U_{0}} & \sigma^{2}_{e}+\tau^{2}_{U_{0}} & \tau^{2}_{U_{0}}\\ \tau^{2}_{U_{0}}& \tau^{2}_{U_{0}} &\tau^{2}_{U_{0}} & \sigma^{2}_{e}+\tau^{2}_{U_{0}} \end{bmatrix} \]

• Unstructured is the most complex structure in the sense there are no constraints

• Always the best fitting random variance structure (akin to an SEM model in which all covariances are modeled) \[ \begin{bmatrix} \sigma^{2}_{11}& \sigma_{12} &\sigma_{13} &\sigma_{14} \\ \sigma_{21}& \sigma^{2}_{22} &\sigma_{23} &\sigma_{24} \\ \sigma_{31}& \sigma_{32} & \sigma^{2}_{33} & \sigma_{34}\\ \sigma_{41}& \sigma_{42} &\sigma_{43} & \sigma^{2}_{44} \end{bmatrix} \]

• 1st Order Auto-Regressive

  • 2 parameters:
    • 1 constant total variance \(\sigma^{2}_{T}\)
    • 1 “AR1” total auto-correlation \(r_{T}\) across occasions
    • \(r^{1}_{T}\) is lag-1 correlation, \(r^{2}_{T}\) is lag-2 correlation, \(r^{3}_{T}\)is lag-3 correlation, and so on
    • The superscripts are exponents, so the correlation at T3 is raised to the second power and so on

\[ \begin{bmatrix} \sigma^{2}_{T} & r^{1}_{T}\sigma^{2}_{T} & r^{2}_{T}\sigma^{2}_{T} & r^{3}_{T}\sigma^{2}_{T} \\ r^{1}_{T}\sigma^{2}_{T} & \sigma^{2}_{T} & r^{1}_{T}\sigma^{2}_{T} & r^{2}_{T}\sigma^{2}_{T} \\ r^{2}_{T}\sigma^{2}_{T} & r^{1}_{T}\sigma^{2}_{T} & \sigma^{2}_{T} & r^{1}_{T}\sigma^{2}_{T} \\ r^{3}_{T}\sigma^{2}_{T} & r^{2}_{T}\sigma^{2}_{T} & r^{1}_{T}\sigma^{2}_{T} & \sigma^{2}_{T} \end{bmatrix} \]

13.10.2 Modeling Different Covarance Structures

• Easiest way is with multilevel regression using nlme

13.10.2.1 Unstructured Matrix

• Note that I use corSymm here to get the unstructured matrix
• I extract the covariance matrix and print:
  • Standardized covariances (correlations)
  • Unstandardized covariances by multiplying the matrix by the residual variance (*un$sigma^2)
    • Note the diagonal is the same value for each time, a consequnces of multiplying the entire matrix by a scalar

un <- gls(gpa ~ as.factor(occas),
                correlation=corSymm(form=~1|student),
                method="REML", data=gpa0)

un.r <- corMatrix(un$modelStruct$corStruct)[[5]]
un.v <- corMatrix(un$modelStruct$corStruct)[[5]]*un$sigma^2

colnames(un.r) <- rownames(un.r) <-  c("GPA1","GPA2","GPA3","GPA4","GPA5","GPA6")
colnames(un.v) <- rownames(un.v) <- c("GPA1","GPA2","GPA3","GPA4","GPA5","GPA6")

rbind(un.r, un.v) %>% 
  as.data.frame %>% 
gt(rownames_to_stub = T, caption = "Unstructured Matrix") %>% 
  tab_row_group(group="Standardized",
                rows = 1:6) %>% 
   tab_row_group(group="Unstandardized",
                rows = 7:12)

Table 13.9: Unstructured Matrix

	GPA1	GPA2	GPA3	GPA4	GPA5	GPA6
Unstandardized
GPA1.1	0.1163	0.0425	0.0279	0.0213	0.0251	0.0246
GPA2.1	0.0425	0.1163	0.0489	0.0532	0.0580	0.0483
GPA3.1	0.0279	0.0489	0.1163	0.0743	0.0804	0.0743
GPA4.1	0.0213	0.0532	0.0743	0.1163	0.0956	0.0941
GPA5.1	0.0251	0.0580	0.0804	0.0956	0.1163	0.1002
GPA6.1	0.0246	0.0483	0.0743	0.0941	0.1002	0.1163
Standardized
GPA1	1.0000	0.3656	0.2399	0.1834	0.2158	0.2119
GPA2	0.3656	1.0000	0.4203	0.4572	0.4987	0.4152
GPA3	0.2399	0.4203	1.0000	0.6387	0.6911	0.6393
GPA4	0.1834	0.4572	0.6387	1.0000	0.8217	0.8089
GPA5	0.2158	0.4987	0.6911	0.8217	1.0000	0.8619
GPA6	0.2119	0.4152	0.6393	0.8089	0.8619	1.0000

13.10.2.2 Compound Symmetry Matrix

• Note the use of the corCompSymm option here

cs <- gls(gpa ~ as.factor(occas),
                correlation=corCompSymm(form=~1|student),
                method="REML", data=gpa0)

cs.r <- corMatrix(cs$modelStruct$corStruct)[[5]]
cs.v <- corMatrix(cs$modelStruct$corStruct)[[5]]*cs$sigma^2

colnames(cs.r) <- rownames(cs.r) <-  c("GPA1","GPA2","GPA3","GPA4","GPA5","GPA6")

colnames(cs.v) <- rownames(cs.v) <- c("GPA1","GPA2","GPA3","GPA4","GPA5","GPA6")

rbind(cs.r, cs.v) %>% 
  as.data.frame %>% 
gt(rownames_to_stub = T, caption = "Compound Symmetry Matrix") %>% 
  tab_row_group(group="Standardized",
                rows = 1:6) %>% 
   tab_row_group(group="Unstandardized",
                rows = 7:12)

Table 13.10: Compound Symmetry Matrix

	GPA1	GPA2	GPA3	GPA4	GPA5	GPA6
Unstandardized
GPA1.1	0.1220	0.0637	0.0637	0.0637	0.0637	0.0637
GPA2.1	0.0637	0.1220	0.0637	0.0637	0.0637	0.0637
GPA3.1	0.0637	0.0637	0.1220	0.0637	0.0637	0.0637
GPA4.1	0.0637	0.0637	0.0637	0.1220	0.0637	0.0637
GPA5.1	0.0637	0.0637	0.0637	0.0637	0.1220	0.0637
GPA6.1	0.0637	0.0637	0.0637	0.0637	0.0637	0.1220
Standardized
GPA1	1.0000	0.5221	0.5221	0.5221	0.5221	0.5221
GPA2	0.5221	1.0000	0.5221	0.5221	0.5221	0.5221
GPA3	0.5221	0.5221	1.0000	0.5221	0.5221	0.5221
GPA4	0.5221	0.5221	0.5221	1.0000	0.5221	0.5221
GPA5	0.5221	0.5221	0.5221	0.5221	1.0000	0.5221
GPA6	0.5221	0.5221	0.5221	0.5221	0.5221	1.0000

13.10.2.3 Autoregressive

• Finally, I use corAR1 to get the autoregressive matrix

ar1 <- gls(gpa ~ as.factor(occas),
                correlation=corAR1(form=~1|student),
                method="REML", data=gpa0)

ar1.r <- corMatrix(ar1$modelStruct$corStruct)[[5]]
ar1.v <- corMatrix(ar1$modelStruct$corStruct)[[5]]*ar1$sigma^2

colnames(ar1.r) <- rownames(ar1.r) <-  c("GPA1","GPA2","GPA3","GPA4","GPA5","GPA6")

colnames(ar1.v) <- rownames(ar1.v) <- c("GPA1","GPA2","GPA3","GPA4","GPA5","GPA6")

rbind(ar1.r, ar1.v) %>% 
  as.data.frame %>% 
gt(rownames_to_stub = T, caption = "Autogressive Matrix") %>% 
  tab_row_group(group="Standardized",
                rows = 1:6) %>% 
   tab_row_group(group="Unstandardized",
                rows = 7:12)

Table 13.11: Autogressive Matrix

	GPA1	GPA2	GPA3	GPA4	GPA5	GPA6
Unstandardized
GPA1.1	0.1215	0.0784	0.0506	0.0326	0.0211	0.0136
GPA2.1	0.0784	0.1215	0.0784	0.0506	0.0326	0.0211
GPA3.1	0.0506	0.0784	0.1215	0.0784	0.0506	0.0326
GPA4.1	0.0326	0.0506	0.0784	0.1215	0.0784	0.0506
GPA5.1	0.0211	0.0326	0.0506	0.0784	0.1215	0.0784
GPA6.1	0.0136	0.0211	0.0326	0.0506	0.0784	0.1215
Standardized
GPA1	1.0000	0.6452	0.4163	0.2686	0.1733	0.1118
GPA2	0.6452	1.0000	0.6452	0.4163	0.2686	0.1733
GPA3	0.4163	0.6452	1.0000	0.6452	0.4163	0.2686
GPA4	0.2686	0.4163	0.6452	1.0000	0.6452	0.4163
GPA5	0.1733	0.2686	0.4163	0.6452	1.0000	0.6452
GPA6	0.1118	0.1733	0.2686	0.4163	0.6452	1.0000

13.10.2.4 Which Covariance Structure Fits Best?

• Use the anova function to compare models
• CS and AR1 have the same number of parameters, so no test
• Can compare each to UN
• Remember, UN is best fitting model by default
  • UN has 22 parameters
    • 15 covariance terms: (n*(n-1))/2 = 6*5/2 = 15
    • 6 coefficients (one for each time/occasion)
    • 1 residual variance (the matrix diagonal)
  • CS has 8 parameters
    • 1 residual variance
    • 1 covariance (the same value in each off-diagonal cell)
    • 6 coefficients (one for each time/occasion)
  • AR has 8 parameters
    • 1 residual variance
    • 1 correlation (rho, which is squared for each successive time interval)
    • 6 coefficients (one for each time/occasion)
  • We compare UN to simpler models, we want \(\chi^2\) not to be significant
  • Recall that models with fewer parameters always reduce model fit–we want that reduction not to be significant

anova(un, cs) 

   Model df AIC BIC logLik   Test L.Ratio p-value
un     1 22 145 257  -50.5                       
cs     2  8 445 486 -214.4 1 vs 2     328  <.0001

anova(un, ar1)

    Model df AIC BIC logLik   Test L.Ratio p-value
un      1 22 145 257  -50.5                       
ar1     2  8 384 425 -184.0 1 vs 2     267  <.0001

13.11 Non-linear Trends: Quadratic Growth Model

• What if we expect or want to test if the growth rate is NOT linear?
• In regression or MLM models, we’d use the square of the time variable
• That is done here by squaring each of the time “loadings” from the original slope (now called s1) and using them in s2 to assess a potential quadratic trend

# quadratic growth--second latent slope's loadings are the square of the first latent slope's loadings
model.quad <- '
# intercept
i =~ 1*Time1 + 1*Time2 + 1*Time3 + 1*Time4

# linear slope
s =~ 0*Time1 + 1*Time2 + 2*Time3 + 3*Time4

# slope 2 (Quadratic)
s2 =~ 0*Time1 + 1*Time2 + 4*Time3 + 9*Time4
'

Figure 13.10: Quadratic LCM

13.11.1 Let’s See an Example

• Mplus data (example 6.9)
• Data in this example are wide–convert to long for graph

pw <- read.table("https://www.statmodel.com/usersguide/chap6/ex6.9.dat")

names(pw) <- c("Time1","Time2","Time3","Time4")

#convert to long
pw %>% 
  mutate(student = 1:n()) %>% #generate student
  pivot_longer(-student, names_to = "time", values_to = "Y") -> pw.long

pw.long %>% 
ggplot(aes(x = time, y = Y, group=student)) +
  geom_line()

• Here’s the model that estimates non-linear trend

# piecewise linear slopes
model.mp <- '
# intercept
i =~ 1*Time1 + 1*Time2 + 1*Time3 + 1*Time4

# linear slope
s =~ 0*Time1 + 1*Time2 + 2*Time3 + 3*Time4

# slope 2 (Quadratic)
s2 =~ 0*Time1 + 1*Time2 + 4*Time3 + 9*Time4

#set covariances to 0
#covariance assumes linearity
i ~~ 0*s
i ~~ 0*s2
s ~~ 0*s2
'

model.mp.out <- growth(model.mp, data=pw)

Table 13.12: Quadratic Growth Model

Quadratic Trend Estimate(Std.Err.) Factor Loadings i Time1 1.00+ Time2 1.00+ Time3 1.00+ Time4 1.00+ s Time1 0.00+ Time2 1.00+ Time3 2.00+ Time4 3.00+ s2 Time1 0.00+ Time2 1.00+ Time3 4.00+ Time4 9.00+ Intercepts Time1 0.00+ Time2 0.00+ Time3 0.00+ Time4 0.00+ Residual Variances Time1 0.83(0.10)*** Time2 0.93(0.09)*** Time3 0.85(0.14)*** Time4 1.60(0.55)** Latent Intercepts i 0.52(0.06)*** s 1.03(0.08)*** s2 0.51(0.03)*** Latent Variances i 1.11(0.11)*** s 0.64(0.12)*** s2 0.20(0.02)*** Latent Covariances i w/s 0.00+ i w/s2 0.00+ s w/s2 0.00+ Fit Indices χ2 5.00(4) CFI 1.00 TLI 1.00 RMSEA 0.02 +Fixed parameter p<0.05, p<0.01, p<0.001

13.12 Piecewise Model

• Interesting question: if slope is nonlinear, how is change analyzed?
  • Piecewise approach “isolates” nonlinear pattern via use of coefficients
• Notice s1 and s2 below
  • s1 says the DV decreases at a steady rate through T4 then flattens
  • s2 says no change through the first four time periods, but consistent increase after
• If interested, see Bonito et al. (2014) for an example of how a piecemeal approach, in an MLM context, might be done
  • Hoffman (2014) is an excellent source, as is her website

13.12.1 Piecewise Example

model.pw <- '
# intercept
i =~ 1*Time1 + 1*Time2 + 1*Time3 + 1*Time4

# linear slope from t1 to t2, flat after
s =~ 0*Time1 + 1*Time2 + 1*Time3 + 1*Time4

# linear slope from t2 to t4, flat to start
s2 =~ 0*Time1 + 0*Time2 + 1*Time3 + 2*Time4

#set covariances to 0
#covariance assumes linearity
i ~~ 0*s
i ~~ 0*s2
s ~~ 0*s2
'

model.pw.out <- growth(model.pw, data=pw)

Table 13.13: Piecewise Model

Piecewise Estimate(Std.Err.) Factor Loadings i Time1 1.00+ Time2 1.00+ Time3 1.00+ Time4 1.00+ s Time1 0.00+ Time2 1.00+ Time3 1.00+ Time4 1.00+ s2 Time1 0.00+ Time2 0.00+ Time3 1.00+ Time4 2.00+ Intercepts Time1 0.00+ Time2 0.00+ Time3 0.00+ Time4 0.00+ Residual Variances Time1 0.76(0.10)*** Time2 0.18(0.12)  Time3 0.74(0.14)*** Time4 4.42(0.57)*** Latent Intercepts i 0.52(0.06)*** s 1.54(0.07)*** s2 2.77(0.09)*** Latent Variances i 1.18(0.12)*** s 1.55(0.18)*** s2 3.50(0.26)*** Latent Covariances i w/s 0.00+ i w/s2 0.00+ s w/s2 0.00+ Fit Indices χ2 155.81(4)*** CFI 0.88 TLI 0.81 RMSEA 0.28 +Fixed parameter p<0.05, p<0.01, p<0.001

13.13 LCM with Nested Data

• We can add a level and nest time within individuals, who are then nested within groups

13.13.1 Graph the Model

• Conceptual graph showing nesting
• Three measurement times
• Two 3-person groups (you’d obviously have more in a study but done here for convenience)
  • L1 is time, which provides the residual variance
  • L2 is person, which provides between-person variance
  • L3 is group, which provides between-groups variance

Figure 13.11: Multilevel Latent Growth Model–Design Version

13.13.2 An Example

• An example from Heck and Thomas (2015)
• Math scores for grades 6, 7, and 8 nested within students, who are nested within schools
• Covariates
  • L3 variables
    • School quality, cses (student SES composition)
  • L2 variables
    • female (females coded 1), lowses (dummy coded, 1 = free lunch program)

13.13.3 Model Graph

• Notice here that the graph looks much like multilevel SEM graphs we’ve seen
• Many modeling choices–here’s I’ve provided only the model with identical structures at the individual and group levels
• Used red edges (arrows) for the intercepts for readability

Figure 13.12: Multilevel Latent Growth Model for Heck and Thomas Data

13.13.4 Read the Data

• Data are in wide format, with each student having one record
• This data setup is appropriate for SEM–long file needed if multivariate regression used

ht7.2 <- read.table("/Volumes/GoogleDrive/My Drive/SAS/Classes/comm640/data/HT2015/ch7datasyntaxes/ch7 ex2.txt")

names(ht7.2) <- c("schcode","math1","math2","math3","female","lowses",
                  "slep","sped","middle","notrans","cses","quality")

head(ht7.2) %>% 
  gt(caption = "Heck and Thomas Math Achievement Data")

Table 13.14: Heck and Thomas Math Achievement Data

schcode	math1	math2	math3	female	lowses	slep	sped	middle	notrans	cses	quality
104	592	631	639	0	1	1	0	1	0	0.97	-1.14
104	570	652	659	0	1	1	0	1	0	0.97	-1.14
104	553	631	645	1	1	1	0	1	0	0.97	-1.14
104	627	695	723	0	1	1	0	1	0	0.97	-1.14
104	580	619	640	0	1	1	0	1	0	0.97	-1.14
104	588	685	657	1	1	1	0	1	0	0.97	-1.14

13.13.5 Descriptive Statistics

13.13.5.1 Individual Level

• Let’s have a look at the multilevel correlation model

ht7.2 %>% 
  select(schcode, math1:math3) %>%
  mlmCorrs::lgm(group="schcode", title = "LGM Matrix for Math Achievement Data")

LGM Matrix for Math Achievement Data
Variable	Mean	SD	1	2	3
1. Math1	598.36	41.56	.06**	.73**	.74**
2. Math2	660.00	37.91	.79**	.09**	.83**
3. Math3	683.49	39.51	.76**	.87**	.09**
Note. *p < .05. **p < .01.

• Plot the data
• Need long file

ht7.2 %>% 
  rowid_to_column("ID") %>% # no subject ID in file--easy to do
  select(math1:math3, ID) %>% 
  pivot_longer(cols=math1:math3, values_to = "math", names_to = "time") %>% 
  ggplot(aes(x=time, y=math, group=ID)) + 
  geom_line()

13.13.5.2 School Level

• Graph by school
• Interesting trajectory at the school level, and looks a bit different compared to that for students

ht7.2 %>% 
  group_by(schcode) %>% 
  summarize(m.math1 = mean(math1), m.math2=mean(math2), m.math3=mean(math3)) %>%
  pivot_longer(cols=m.math1:m.math3, values_to = "math", names_to = "time") %>% 
  ggplot(aes(x=time, y=math, group=schcode)) + 
  geom_line()

13.13.6 Quick Analysis

• These models sometimes have difficulty converging
• This is shown for illustrative purposes–more work is needed to find a well-behaved model

13.13.6.1 Model Without Nesting

• Note the slope had negative variance
• Set slope variance and slope/intercept covariance to 0

ht.model0 <- '
iw =~ 1*math1 + 1*math2 + 1*math3 
sw =~ 0*math1 + 1*math2 +  2*math3 
sw~~0*sw
iw~~0*sw
'

ht.model0.out <- growth(ht.model0, data=ht7.2)

summary(ht.model0.out)

lavaan 0.6-11 ended normally after 153 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         6
                                                      
  Number of observations                          7029
                                                      
Model Test User Model:
                                                      
  Test statistic                              5414.872
  Degrees of freedom                                 3
  P-value (Chi-square)                           0.000

Parameter Estimates:

  Standard errors                             Standard
  Information                                 Expected
  Information saturated (h1) model          Structured

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)
  iw =~                                               
    math1             1.000                           
    math2             1.000                           
    math3             1.000                           
  sw =~                                               
    math1             0.000                           
    math2             1.000                           
    math3             2.000                           

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)
  iw ~~                                               
    sw                0.000                           

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
   .math1             0.000                           
   .math2             0.000                           
   .math3             0.000                           
    iw              607.523    0.499 1216.695    0.000
    sw               39.508    0.165  239.077    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
    sw                0.000                           
   .math1           626.232   13.008   48.141    0.000
   .math2           443.120   10.390   42.647    0.000
   .math3           208.498    7.754   26.888    0.000
    iw             1276.676   23.589   54.121    0.000

13.13.6.2 What the Nested Version Would Look like

• Heck and Thomas had trouble with this model in Mplus
• I couldn’t get it to produce useable/interpretable results
• Some options below

ht.model <- '

level: 1
#used this for cross level constraints
#iw =~ 1*math1 + a*math1 + 1*math2 + b*math2 + 1*math3 + c*math3
#sw =~ 0*math1 + d*math1 + 1*math2 + e*math2 + 2*math3 + f*math3

#freely estimated latent variables
iw =~ 1*math1 + 1*math2 + 1*math3 
sw =~ 0*math1 + 1*math2 +  2*math3 

sw~~0*sw
sw~~0*iw

level: 2

#tried L1 model with L2 covariances--didnt work
#math1 ~~ math2 + math3
#math2 ~~ math3

#cross level constraints
#ib =~ 1*math1 + a*math1 + 1*math2 + b*math2 + 1*math3 + c*math3
#sb =~ 0*math1 + d*math1 + 1*math2 + e*math2 + 2*math3 + f*math3

#freely estimated
ib =~ 1*math1 + 1*math2 + 1*math3 
sb =~ 0*math1 + 1*math2 + 2*math3 

#playing with residual variances
#math1~~0*math1
#math2~~0*math2
#math3~~0*math3
ib~~0*ib
ib~~0*sb
'

ht.model.out <- growth(ht.model, data=ht7.2, cluster = "schcode")