5  Longitudinal Models with lavaan

5.1 Introduction

In longitudinal models, time is usually considered as a continuous variable.

It is common to examine - the fixed effects: what is the nature of the average trend (linear, quadratic)

• the random effects: individual differences

In addition, these individual differences can be related to

• time-invariant covariates (age, sex, …)

• time-varying covariates (measured at each time point)

SEM can be closely related to ‘mixed models’ (linear mixed models, generalised mixed models) but requires balanced data. The ability of SEM is the addition of indirect paths and latent variables.

5.2 Grow Curve Modelling (GCM)

5.2.1 Back to the example

The teacher’s relation was measured at six times.

The data are available here

library(readxl)
library(semTools)

Warning: le package 'semTools' a été compilé avec la version R 4.3.3

Le chargement a nécessité le package : lavaan

This is lavaan 0.6-17
lavaan is FREE software! Please report any bugs.

###############################################################################

This is semTools 0.5-6

All users of R (or SEM) are invited to submit functions or ideas for functions.

###############################################################################

library(lmerTest)

Le chargement a nécessité le package : lme4

Le chargement a nécessité le package : Matrix

Attachement du package : 'lmerTest'

L'objet suivant est masqué depuis 'package:lme4':

    lmer

L'objet suivant est masqué depuis 'package:stats':

    step

library(dplyr)

Attachement du package : 'dplyr'

Les objets suivants sont masqués depuis 'package:stats':

    filter, lag

Les objets suivants sont masqués depuis 'package:base':

    intersect, setdiff, setequal, union

df<-read_excel("data_long.xlsx")

Long format data for adjustment with lmer models.

df_long <- reshape(data = df, 
                  idvar="id",
                  v.names="TR",
                  varying  = paste("TR",1:6,sep=""),
                  timevar = "time",
                  times = 1:6,
                  direction = "long")

Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
Warning in value[[jj]][ri] <- if (is.factor(xij)) as.vector(xij) else xij: le
nombre d'objets à remplacer n'est pas multiple de la taille du remplacement

boxplot(TR~time,data=df_long)

fit1<-lmer(TR~1+(1|time),
           data=df_long,
           REML=F)

boundary (singular) fit: see help('isSingular')

summary(fit1)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: TR ~ 1 + (1 | time)
   Data: df_long

     AIC      BIC   logLik deviance df.resid 
  1108.9   1119.0   -551.4   1102.9      216 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1789 -0.4637  0.1360  0.6708  1.8020 

Random effects:
 Groups   Name        Variance Std.Dev.
 time     (Intercept) 0.000    0.000   
 Residual             9.007    3.001   
Number of obs: 219, groups:  time, 6

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)   5.5917     0.2028 219.0000   27.57   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

On lavaan:

model<-'
i =~ 1*TR1+1*TR2+1*TR3+1*TR4+1*TR5+1*TR6

# zero intercepts for observed variables

TR1+TR2+TR3+TR4+TR5+TR6 ~0*1

# fixed effect on intercept
i~1

# random effect on intercept
i~~i

# force same variance
TR1 ~~ v*TR1
TR2 ~~ v*TR2
TR3 ~~ v*TR3
TR4 ~~ v*TR4
TR5 ~~ v*TR5
TR6 ~~ v*TR6
'

library(lavaan)
fit2<-lavaan(model,data=df)
summary(fit2)

lavaan 0.6.17 ended normally after 14 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         8
  Number of equality constraints                     5

  Number of observations                           219

Model Test User Model:
                                                      
  Test statistic                               220.165
  Degrees of freedom                                24
  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 =~                                                
    TR1               1.000                           
    TR2               1.000                           
    TR3               1.000                           
    TR4               1.000                           
    TR5               1.000                           
    TR6               1.000                           

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
   .TR1               0.000                           
   .TR2               0.000                           
   .TR3               0.000                           
   .TR4               0.000                           
   .TR5               0.000                           
   .TR6               0.000                           
    i                 6.371    0.166   38.341    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
    i                 5.322    0.579    9.197    0.000
   .TR1        (v)    4.347    0.186   23.399    0.000
   .TR2        (v)    4.347    0.186   23.399    0.000
   .TR3        (v)    4.347    0.186   23.399    0.000
   .TR4        (v)    4.347    0.186   23.399    0.000
   .TR5        (v)    4.347    0.186   23.399    0.000
   .TR6        (v)    4.347    0.186   23.399    0.000

5.3 Random intercept and fixed slope

On lmer

fit1<-lmer(TR~1+time+(1|time),
           data=df_long,
           REML=F)

boundary (singular) fit: see help('isSingular')

summary(fit1)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: TR ~ 1 + time + (1 | time)
   Data: df_long

     AIC      BIC   logLik deviance df.resid 
  1108.5   1122.0   -550.2   1100.5      215 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1873 -0.4518  0.1512  0.6872  1.8263 

Random effects:
 Groups   Name        Variance Std.Dev.
 time     (Intercept) 0.00     0.000   
 Residual             8.91     2.985   
Number of obs: 219, groups:  time, 6

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)   4.9201     0.4787 219.0000  10.278   <2e-16 ***
time          0.6285     0.4063 219.0000   1.547    0.123    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
     (Intr)
time -0.907
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

On lavaan:

model<-'
i =~ 1*TR1+1*TR2+1*TR3+1*TR4+1*TR5+1*TR6
s =~ 1*TR1+2*TR2+3*TR3+4*TR4+5*TR5+6*TR6

#fixed effect on intercept and slope
i~1
s~1

# random intercept 
i~~i

# fixed slope
s~~ 0*s # no variance
s~~ 0*i # no covariance

# force same variance
TR1 ~~ v*TR1
TR2 ~~ v*TR2
TR3 ~~ v*TR3
TR4 ~~ v*TR4
TR5 ~~ v*TR5
TR6 ~~ v*TR6
'

fit<-lavaan(model,data=df)
summary(fit)

lavaan 0.6.17 ended normally after 17 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         9
  Number of equality constraints                     5

  Number of observations                           219

Model Test User Model:
                                                      
  Test statistic                               146.622
  Degrees of freedom                                23
  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 =~                                                
    TR1               1.000                           
    TR2               1.000                           
    TR3               1.000                           
    TR4               1.000                           
    TR5               1.000                           
    TR6               1.000                           
  s =~                                                
    TR1               1.000                           
    TR2               2.000                           
    TR3               3.000                           
    TR4               4.000                           
    TR5               5.000                           
    TR6               6.000                           

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

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
    i                 5.377    0.202   26.683    0.000
    s                 0.284    0.033    8.722    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
    i                 5.369    0.579    9.280    0.000
    s                 0.000                           
   .TR1        (v)    4.065    0.174   23.399    0.000
   .TR2        (v)    4.065    0.174   23.399    0.000
   .TR3        (v)    4.065    0.174   23.399    0.000
   .TR4        (v)    4.065    0.174   23.399    0.000
   .TR5        (v)    4.065    0.174   23.399    0.000
   .TR6        (v)    4.065    0.174   23.399    0.000