Linear Growth Model with Time-Invariant Covariates - Multilevel & SEM Implementation in R

Overview

This tutorial illustrates fitting of linear growth models with time-invariant covariates in the multilevel and SEM frameworks in R.

Example data and code are drawn from Chapter 5 of Grimm, Ram, and Estabrook (2017). Specifically, using the NLSY-CYA Dataset we examine how change in children’s mathematics achievement across grade is related to low birth weight and early (kindergarten to first grade) antisocial behaviors. Please see the book chapter for additional interpretations and insights about the analyses.

Preliminaries - Loading libraries used in this script.

library(psych)  #for basic functions
library(plyr)   #for data management
library(ggplot2)  #for plotting
library(nlme) #for mixed effects models
library(lme4) #for mixed effects models
library(lavaan) #for SEM 
library(semPlot) #for making SEM diagrams

Preliminaries - Data Preparation and Description

For our examples, we use the mathematics achievement scores from the NLSY-CYA Long Data.

Load the repeated measures data

#set filepath for data file
filepath <- "https://raw.githubusercontent.com/LRI-2/Data/main/GrowthModeling/nlsy_math_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_long <- dat  

#Add names the columns of the data set
names(nlsy_math_long) = c('id'     , 'female', 'lb_wght', 
                          'anti_k1', 'math'  , 'grade'  ,
                          'occ'    , 'age'   , 'men'    ,
                          'spring' , 'anti')


#view the first few observations in the data set 
head(nlsy_math_long, 10)
##      id female lb_wght anti_k1 math grade occ age men spring anti
## 1   201      1       0       0   38     3   2 111   0      1    0
## 2   201      1       0       0   55     5   3 135   1      1    0
## 3   303      1       0       1   26     2   2 121   0      1    2
## 4   303      1       0       1   33     5   3 145   0      1    2
## 5  2702      0       0       0   56     2   2 100  NA      1    0
## 6  2702      0       0       0   58     4   3 125  NA      1    2
## 7  2702      0       0       0   80     8   4 173  NA      1    2
## 8  4303      1       0       0   41     3   2 115   0      0    1
## 9  4303      1       0       0   58     4   3 135   0      1    2
## 10 5002      0       0       4   46     4   2 117  NA      1    4

Our specific interest is intraindividual change in the repeated measures of math change across grade, and how the interindividual differences in those trajectories (intercept and linear slope) are related to lb_wght and anti_k1.

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. Here, we want to examine the data to see how the repeated measures of math are structured with respect to age.

Longitudinal Plot of Math across Grade at Testing

#intraindividual change trajetories
ggplot(data=nlsy_math_long,                    #data set
       aes(x = grade, y = math, group = id)) + #setting variables
  geom_point(size=.5) + #adding points to plot
  geom_line(aes(linetype=factor(lb_wght), alpha= anti_k1)) +  #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")

Note that we are able to get all variables of interest into one plot. Cooll!

Multilevel: Linear Growth Model with Time-Invariant Covariates

To capture the systematic growth in the mathematics scores over grade, we use the grade variable, so that the rate of growth indicates (linear) school-related change in mathematics scores.

Using the notation from the lecture, the model equations are …

\[y_{it} = b_{1i}1 + b_{2i}grade_{it} + u_{it}\]

same as we used previously. However, now we also include person-level (Level-2) predictors for the person-specific intercepts captured by \(b_{1i}\) and linear slopes captured by \(b_{2i}\). Specifically, \[b_{1i} = \beta_{01} + \beta_{11}LowBirthWeight_{i} + \beta_{21}AntiSocial_{i} + d_{1i}\] and \[b_{2i} = \beta_{02} + \beta_{12}LowBirthWeight_{i} + \beta_{22}AntiSocial_{i} + d_{2i}\]

Substitution of one equation into the other we get … \[y_{it} = (\beta_{01} + \beta_{11}LowBirthWeight_{i} + \beta_{21}AntiSocial_{i} + d_{1i}) + (\beta_{02} + \beta_{12}LowBirthWeight_{i} + \beta_{22}AntiSocial_{i} + d_{2i})grade_{it} + u_{it}\]

where \[u_{it} \tilde{} N(0,\sigma^{2}_{u})\], and \[\left[\begin{array} {r} d_{1i} \\ d_{2i} \end{array}\right] \tilde{} N(\left[\begin{array} {r} 0 \\ 0 \end{array}\right] ,\left[\begin{array} {r} \sigma^{2}_{d_1} & \sigma_{d_1d_2} \\ \sigma_{d_1d_2} & \sigma^{2}_{d_2} \end{array}\right]\].

This is then the model that gets coded into the multilevel modeling functions. This can be done in a variety of ways. As usual, we concentrate on use of the nlme:nlme() function.

Linear Growth Model with Time-Invariant Covariates using nlme() function

This is then the model that gets coded into the multilevel modeling functions. This can be done in a variety of ways. We concentrate on use of the nlme:nlme() function, with an explicit coding style.

Note that in fitting the linear growth model, we recenter the time variable, so that the intercept is the expected score at Grade 2 (the beginning of the data).

#fitting linear growth model with tic and assigning it to an object
lg_math_tic_nlme <- nlme(math ~ (beta_01 + beta_11*lb_wght + beta_21*anti_k1 + d_1i) +
                                (beta_02 + beta_12*lb_wght + beta_22*anti_k1 +d_2i)*(grade-2),
                          data=nlsy_math_long,
                          fixed=beta_01+beta_11+beta_21+beta_02+beta_12+beta_22~1,
                          random=d_1i+d_2i~1,
                          group=~id,
                          start=c(beta_01=30, beta_11=0, beta_21=0, 
                                  beta_02=4, beta_12=0, beta_22=0),
                          na.action=na.exclude)

#obtaining summary of the model using the object we just created
summary(lg_math_tic_nlme)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: math ~ (beta_01 + beta_11 * lb_wght + beta_21 * anti_k1 + d_1i) +      (beta_02 + beta_12 * lb_wght + beta_22 * anti_k1 + d_2i) *          (grade - 2) 
##   Data: nlsy_math_long 
##        AIC     BIC   logLik
##   15943.14 16000.2 -7961.57
## 
## Random effects:
##  Formula: list(d_1i ~ 1, d_2i ~ 1)
##  Level: id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev    Corr  
## d_1i     7.9413099 d_1i  
## d_2i     0.8445068 -0.012
## Residual 6.0213575       
## 
## Fixed effects:  beta_01 + beta_11 + beta_21 + beta_02 + beta_12 + beta_22 ~ 1 
##            Value Std.Error   DF  t-value p-value
## beta_01 36.28987 0.4969765 1284 73.02130  0.0000
## beta_11 -2.71621 1.2953408 1284 -2.09691  0.0362
## beta_21 -0.55088 0.2327789 1284 -2.36655  0.0181
## beta_02  4.31516 0.1207520 1284 35.73572  0.0000
## beta_12  0.62464 0.3335739 1284  1.87256  0.0614
## beta_22 -0.01929 0.0589361 1284 -0.32731  0.7435
##  Correlation: 
##         bet_01 bet_11 bet_21 bet_02 bet_12
## beta_11 -0.194                            
## beta_21 -0.671 -0.025                     
## beta_02 -0.529  0.096  0.358              
## beta_12  0.091 -0.532  0.026 -0.168       
## beta_22  0.343  0.026 -0.529 -0.660 -0.055
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.080148313 -0.525141719 -0.008659962  0.530786044  2.534567886 
## 
## Number of Observations: 2221
## Number of Groups: 932

We get information about how the model fits the data, as well as some information about the relations of the trajectories with the time-invariant covariates. Specifically, we see from the model that there is evidence that differneces in low birth weight status and antisocial behavior are related to individual differences in intercept, but not in linear slope. More detailed interpetation is given in the accompanying video.

Plotting the predicted trajectories.

#obtaining predicted scores for individuals
nlsy_math_long$pred_lg <- predict(lg_math_tic_nlme)

#obtaining predicted scores for prototype
nlsy_math_long$proto_lg <- predict(lg_math_tic_nlme, level=0)

#intraindividual change trajetories
#plotting predicted trajectory
ggplot(data = nlsy_math_long, aes(x = grade, y = pred_lg, group = id)) +
  geom_point(color="gray", size=.5) + 
  geom_line(aes(linetype=factor(lb_wght), alpha= anti_k1)) +  #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 = "Predicted PIAT Mathematics")

Linear Growth Model with Time-Invariant Covariates using lmer() function (in the lme4 library)

The model can also be fit in other multilevel package, here lme4.

#making the recentered time variable
nlsy_math_long$grade_c2 <- nlsy_math_long$grade-2

#fitting linear growth model and assigning it to an object
lg_math_tic_lmer <- lmer(math ~ 1 + lb_wght + anti_k1 + 
                                grade_c2 + I(grade_c2*lb_wght) + I(grade_c2*anti_k1)  + 
                           (1 + grade_c2 | id), 
                     data = nlsy_math_long, 
                     REML = TRUE,
                     na.action = na.exclude)

summary(lg_math_tic_lmer)
## Linear mixed model fit by REML ['lmerMod']
## Formula: math ~ 1 + lb_wght + anti_k1 + grade_c2 + I(grade_c2 * lb_wght) +  
##     I(grade_c2 * anti_k1) + (1 + grade_c2 | id)
##    Data: nlsy_math_long
## 
## REML criterion at convergence: 15930.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.07031 -0.52511 -0.00901  0.53039  2.53205 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  id       (Intercept) 63.4359  7.9647        
##           grade_c2     0.7344  0.8569   -0.02
##  Residual             36.2650  6.0220        
## Number of obs: 2221, groups:  id, 932
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)           36.28899    0.49725  72.979
## lb_wght               -2.71513    1.29602  -2.095
## anti_k1               -0.55066    0.23290  -2.364
## grade_c2               4.31611    0.12089  35.703
## I(grade_c2 * lb_wght)  0.62408    0.33389   1.869
## I(grade_c2 * anti_k1) -0.01945    0.05899  -0.330
## 
## Correlation of Fixed Effects:
##             (Intr) lb_wgh ant_k1 grd_c2 I(_2*l
## lb_wght     -0.194                            
## anti_k1     -0.671 -0.025                     
## grade_c2    -0.529  0.097  0.358              
## I(grd_2*l_)  0.091 -0.532  0.025 -0.168       
## I(gr_2*a_1)  0.343  0.026 -0.529 -0.660 -0.055

SEM: Linear Growth Model with Time-Invariant Covariates

Preliminaries - Data Preparation

For the SEM implementation, we use the mathematics achievement scores and time-invariant covariates from the NLSY-CYA Wide Data.

Load the repeated measures data

#set filepath for data file
filepath <- "https://raw.githubusercontent.com/LRI-2/Data/main/GrowthModeling/nlsy_math_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_wide <- dat  

# Give the variable names
names(nlsy_math_wide)=c('id'     ,'female' ,'lb_wght','anti_k1',
                        'math2'  ,'math3'  ,'math4'  ,'math5'  ,'math6'  ,'math7'  ,'math8'  ,
                        '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')


#view the first few observations (and columns) in the data set 
head(nlsy_math_wide[ ,1:11], 10)
##      id female lb_wght anti_k1 math2 math3 math4 math5 math6 math7 math8
## 1   201      1       0       0    NA    38    NA    55    NA    NA    NA
## 2   303      1       0       1    26    NA    NA    33    NA    NA    NA
## 3  2702      0       0       0    56    NA    58    NA    NA    NA    80
## 4  4303      1       0       0    NA    41    58    NA    NA    NA    NA
## 5  5002      0       0       4    NA    NA    46    NA    54    NA    66
## 6  5005      1       0       0    35    NA    50    NA    60    NA    59
## 7  5701      0       0       2    NA    62    61    NA    NA    NA    NA
## 8  6102      0       0       0    NA    NA    55    67    NA    81    NA
## 9  6801      1       0       0    NA    54    NA    62    NA    66    NA
## 10 6802      0       0       0    NA    55    NA    66    NA    68    NA

Our specific interest is in change across the the repeated measures of math, math2 to math8 and how those changes are related to individual differences in lb_wght and anti_k1.

Linear Growth Model with Time-Invariant Covariates in lavaan

Specifying the SEM model

#writing out linear growth model with tic in full SEM way 
lg_math_tic_lavaan_model <- '
    #latent variable definitions
            #intercept
              eta1 =~ 1*math2+
                      1*math3+
                      1*math4+
                      1*math5+
                      1*math6+
                      1*math7+
                      1*math8
            #linear slope
              eta2 =~ 0*math2+
                      1*math3+
                      2*math4+
                      3*math5+
                      4*math6+
                      5*math7+
                      6*math8

          #factor variances
            eta1 ~~ eta1
            eta2 ~~ eta2

          #factor covariance
            eta1 ~~ eta2

          #manifest variances (set equal by naming theta)
            math2 ~~ theta*math2
            math3 ~~ theta*math3
            math4 ~~ theta*math4
            math5 ~~ theta*math5
            math6 ~~ theta*math6
            math7 ~~ theta*math7
            math8 ~~ theta*math8

          #latent means (freely estimated)
            eta1 ~ 1
            eta2 ~ 1

          #manifest means (fixed to zero)
            math2 ~ 0*1
            math3 ~ 0*1
            math4 ~ 0*1
            math5 ~ 0*1
            math6 ~ 0*1
            math7 ~ 0*1
            math8 ~ 0*1

        #Time invariant covaraite
          #regression of time-invariant covariate on intercept and slope factors
            eta1~lb_wght+anti_k1
            eta2~lb_wght+anti_k1

        #variance of TIV covariates
            lb_wght ~~ lb_wght
            anti_k1 ~~ anti_k1

        #covariance of TIV covaraites
            lb_wght ~~ anti_k1

        #means of TIV covariates (freely estimated)
            lb_wght ~ 1
            anti_k1 ~ 1
' #end of model definition

And fit the model to the data … using maximum likelihood estimator and full information maximum likelihood to handle missing values.

#estimating the model using sem() function
lg_math_tic_lavaan_fit <- sem(lg_math_tic_lavaan_model, 
                          data = nlsy_math_wide,
                          meanstructure = TRUE,
                          estimator = "ML",
                          missing = "fiml")
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan
## WARNING: some observed variances are (at least) a factor 1000 times larger than
## others; use varTable(fit) to investigate
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan
## WARNING: due to missing values, some pairwise combinations have less than 10%
## coverage
## Warning in lav_mvnorm_missing_h1_estimate_moments(Y = X[[g]], wt = WT[[g]], : lavaan WARNING:
##     Maximum number of iterations reached when computing the sample
##     moments using EM; use the em.h1.iter.max= argument to increase the
##     number of iterations

We can look more specifically at the numerical summaries of the model.

#obtaining summary of the model using the object we just created
summary(lg_math_tic_lavaan_fit, fit.measures=TRUE)
## lavaan 0.6-9 ended normally after 113 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        21
##   Number of equality constraints                     6
##                                                       
##   Number of observations                           933
##   Number of missing patterns                        61
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                               220.221
##   Degrees of freedom                                39
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                               892.616
##   Degrees of freedom                                36
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.788
##   Tucker-Lewis Index (TLI)                       0.805
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -9785.085
##   Loglikelihood unrestricted model (H1)      -9674.975
##                                                       
##   Akaike (AIC)                               19600.171
##   Bayesian (BIC)                             19672.747
##   Sample-size adjusted Bayesian (BIC)        19625.108
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.071
##   90 Percent confidence interval - lower         0.062
##   90 Percent confidence interval - upper         0.080
##   P-value RMSEA <= 0.05                          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|)
##   eta1 =~                                             
##     math2             1.000                           
##     math3             1.000                           
##     math4             1.000                           
##     math5             1.000                           
##     math6             1.000                           
##     math7             1.000                           
##     math8             1.000                           
##   eta2 =~                                             
##     math2             0.000                           
##     math3             1.000                           
##     math4             2.000                           
##     math5             3.000                           
##     math6             4.000                           
##     math7             5.000                           
##     math8             6.000                           
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   eta1 ~                                              
##     lb_wght          -2.716    1.294   -2.099    0.036
##     anti_k1          -0.551    0.232   -2.369    0.018
##   eta2 ~                                              
##     lb_wght           0.625    0.333    1.873    0.061
##     anti_k1          -0.019    0.059   -0.327    0.743
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##  .eta1 ~~                                             
##    .eta2             -0.078    1.145   -0.068    0.945
##   lb_wght ~~                                          
##     anti_k1           0.007    0.014    0.548    0.584
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .eta1             36.290    0.497   73.052    0.000
##    .eta2              4.315    0.122   35.420    0.000
##    .math2             0.000                           
##    .math3             0.000                           
##    .math4             0.000                           
##    .math5             0.000                           
##    .math6             0.000                           
##    .math7             0.000                           
##    .math8             0.000                           
##     lb_wght           0.080    0.009    9.031    0.000
##     anti_k1           1.454    0.050   29.216    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .eta1             63.064    5.609   11.242    0.000
##    .eta2              0.713    0.326    2.185    0.029
##    .math2   (thet)   36.257    1.868   19.406    0.000
##    .math3   (thet)   36.257    1.868   19.406    0.000
##    .math4   (thet)   36.257    1.868   19.406    0.000
##    .math5   (thet)   36.257    1.868   19.406    0.000
##    .math6   (thet)   36.257    1.868   19.406    0.000
##    .math7   (thet)   36.257    1.868   19.406    0.000
##    .math8   (thet)   36.257    1.868   19.406    0.000
##     lb_wght           0.074    0.003   21.599    0.000
##     anti_k1           2.312    0.107   21.599    0.000
#Other summaries
#parameterEstimates(lg_math_tic_lavaan_fit)
#inspect(lg_math_tic_lavaan_fit, what="est")

We get information about how the model fits the data, as well as some information about the data itself. Interpetation is covered in more detail in the accompanying video.

We can use the fitted model object to make a diagram (and also check if we specified and estimated the model as intended).

#diagram of fitted model
semPaths(lg_math_tic_lavaan_fit,what = "path", whatLabels = "par")

Testing the significance of the Time-Invariant Covariates in lavaan

To determine whether the addition of low birth weight and kindergarten antisocial behaviors is useful for model fit, we can fit the linear growth model with regression effects set = zero.

Specifying the SEM model with predictor effects at 0.

#writing out linear growth model with tic in full SEM way 
lg_math_ticZERO_lavaan_model <- '
    #latent variable definitions
            #intercept
              eta1 =~ 1*math2+
                      1*math3+
                      1*math4+
                      1*math5+
                      1*math6+
                      1*math7+
                      1*math8
            #linear slope
              eta2 =~ 0*math2+
                      1*math3+
                      2*math4+
                      3*math5+
                      4*math6+
                      5*math7+
                      6*math8

          #factor variances
            eta1 ~~ eta1
            eta2 ~~ eta2

          #factor covariance
            eta1 ~~ eta2

          #manifest variances (set equal by naming theta)
            math2 ~~ theta*math2
            math3 ~~ theta*math3
            math4 ~~ theta*math4
            math5 ~~ theta*math5
            math6 ~~ theta*math6
            math7 ~~ theta*math7
            math8 ~~ theta*math8

          #latent means (freely estimated)
            eta1 ~ 1
            eta2 ~ 1

          #manifest means (fixed to zero)
            math2 ~ 0*1
            math3 ~ 0*1
            math4 ~ 0*1
            math5 ~ 0*1
            math6 ~ 0*1
            math7 ~ 0*1
            math8 ~ 0*1

        #Time invariant covaraite
          #regression of time-invariant covariate on intercept and slope factors
          #FIXED to 0
            eta1~0*lb_wght+0*anti_k1
            eta2~0*lb_wght+0*anti_k1

        #variance of TIV covariates
            lb_wght ~~ lb_wght
            anti_k1 ~~ anti_k1

        #covariance of TIV covaraites
            lb_wght ~~ anti_k1

        #means of TIV covariates (freely estimated)
            lb_wght ~ 1
            anti_k1 ~ 1
' #end of model definition

And fit the model to the data … using maximum likelihood estimator and full information maximum likelihood to handle missing values.

#estimating the model using sem() function
lg_math_ticZERO_lavaan_fit <- sem(lg_math_ticZERO_lavaan_model, 
                          data = nlsy_math_wide,
                          meanstructure = TRUE,
                          estimator = "ML",
                          missing = "fiml")
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan
## WARNING: some observed variances are (at least) a factor 1000 times larger than
## others; use varTable(fit) to investigate
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan
## WARNING: due to missing values, some pairwise combinations have less than 10%
## coverage
## Warning in lav_mvnorm_missing_h1_estimate_moments(Y = X[[g]], wt = WT[[g]], : lavaan WARNING:
##     Maximum number of iterations reached when computing the sample
##     moments using EM; use the em.h1.iter.max= argument to increase the
##     number of iterations

We can look more specifically at the numerical summaries of the model.

#obtaining summary of the model using the object we just created
summary(lg_math_ticZERO_lavaan_fit, fit.measures=TRUE)
## lavaan 0.6-9 ended normally after 90 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        17
##   Number of equality constraints                     6
##                                                       
##   Number of observations                           933
##   Number of missing patterns                        61
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                               234.467
##   Degrees of freedom                                43
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                               892.616
##   Degrees of freedom                                36
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.776
##   Tucker-Lewis Index (TLI)                       0.813
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -9792.208
##   Loglikelihood unrestricted model (H1)      -9674.975
##                                                       
##   Akaike (AIC)                               19606.416
##   Bayesian (BIC)                             19659.638
##   Sample-size adjusted Bayesian (BIC)        19624.703
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.069
##   90 Percent confidence interval - lower         0.061
##   90 Percent confidence interval - upper         0.078
##   P-value RMSEA <= 0.05                          0.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.103
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Observed
##   Observed information based on                Hessian
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   eta1 =~                                             
##     math2             1.000                           
##     math3             1.000                           
##     math4             1.000                           
##     math5             1.000                           
##     math6             1.000                           
##     math7             1.000                           
##     math8             1.000                           
##   eta2 =~                                             
##     math2             0.000                           
##     math3             1.000                           
##     math4             2.000                           
##     math5             3.000                           
##     math6             4.000                           
##     math7             5.000                           
##     math8             6.000                           
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   eta1 ~                                              
##     lb_wght           0.000                           
##     anti_k1           0.000                           
##   eta2 ~                                              
##     lb_wght           0.000                           
##     anti_k1           0.000                           
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##  .eta1 ~~                                             
##    .eta2             -0.181    1.150   -0.158    0.875
##   lb_wght ~~                                          
##     anti_k1           0.007    0.014    0.548    0.584
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .eta1             35.267    0.355   99.229    0.000
##    .eta2              4.339    0.088   49.136    0.000
##    .math2             0.000                           
##    .math3             0.000                           
##    .math4             0.000                           
##    .math5             0.000                           
##    .math6             0.000                           
##    .math7             0.000                           
##    .math8             0.000                           
##     lb_wght           0.080    0.009    9.031    0.000
##     anti_k1           1.454    0.050   29.216    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .eta1             64.562    5.659   11.408    0.000
##    .eta2              0.733    0.327    2.238    0.025
##    .math2   (thet)   36.230    1.867   19.410    0.000
##    .math3   (thet)   36.230    1.867   19.410    0.000
##    .math4   (thet)   36.230    1.867   19.410    0.000
##    .math5   (thet)   36.230    1.867   19.410    0.000
##    .math6   (thet)   36.230    1.867   19.410    0.000
##    .math7   (thet)   36.230    1.867   19.410    0.000
##    .math8   (thet)   36.230    1.867   19.410    0.000
##     lb_wght           0.074    0.003   21.599    0.000
##     anti_k1           2.312    0.107   21.599    0.000
#Other summaries
#parameterEstimates(lg_math_tic_lavaan_fit)
#inspect(lg_math_tic_lavaan_fit, what="est")

And can the new constrained model as a diagram (and also check if we specified and estimated the model as intended).

#diagram of fitted model
semPaths(lg_math_ticZERO_lavaan_fit,what = "path", whatLabels = "par")

The chi-square difference test (obtained from the output or using the anova() function) indicates that the time invariant predictors are useful to the model. For a more in-depth discussion of model fit and comparison, see pages 110-111 in the book.

anova(lg_math_tic_lavaan_fit,lg_math_ticZERO_lavaan_fit)
## Chi-Squared Difference Test
## 
##                            Df   AIC   BIC  Chisq Chisq diff Df diff Pr(>Chisq)
## lg_math_tic_lavaan_fit     39 19600 19673 220.22                              
## lg_math_ticZERO_lavaan_fit 43 19606 19660 234.47     14.245       4   0.006552
##                              
## lg_math_tic_lavaan_fit       
## lg_math_ticZERO_lavaan_fit **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusion

This tutorial has presented how time invariant covariates can be added to the linear growth model in both the multilevel and SEM frameworks.

Examining interindividual differences in intraindividual change! Bring it!