Growth Curve Modeling

About the data

The data are from a sleep deprivation study described in Belenky et al. (2003) [Journal of Sleep Research, 12, 1-12]. Eighteen participants were studied over a 10 day period. On day 0 the subjects had their normal amount of sleep. Starting that night they were restricted to 3 hours of sleep per night. Each day the participants were given a series of reaction time tests and the average was computed. The response variable is the average reaction time in ms. There are 10 columns of reaction time with column labels, Reaction.0, Reaction.1, …, Reaction.9. The other variables are SubNum (subject number), female (1 = female, 0 = male) and gpa. Two participants, 105 and 116, dropped out after day 6 and their missing data are coded with -999.

We will use the lavaan package in R to do growth curve modeling. We’ll compare our findings to what we would obtain in using lme4. The data are in correct format for lavaan but not for lme4. We’ll need to reshape the data into long format for lme4 and for plotting.

Reading in the data

sleep.wide <- read.table("/Users/cdesjard/Google Drive/teaching-courses/ICD_summer2016/data/sleeplab.txt", na.strings = "-999", header = T)
head(sleep.wide)

  SubNum female      gpa Reaction.0 Reaction.1 Reaction.2 Reaction.3 Reaction.4 Reaction.5 Reaction.6 Reaction.7
1    101      0 3.090871   249.5600   258.7047   250.8006   321.4398   356.8519   414.6901   382.2038   290.1486
2    102      1 2.572551   222.7339   205.2658   202.9778   204.7070   207.7161   215.9618   213.6303   217.7272
3    103      1 2.015587   199.0539   194.3322   234.3200   232.8416   229.3074   220.4579   235.4208   255.7511
4    104      1 2.524083   321.5426   300.4002   283.8565   285.1330   285.7973   297.5855   280.2396   318.2613
5    105      0 3.077358   287.6079   285.0000   301.8206   320.1153   316.2773   293.3187   290.0750         NA
6    106      0 3.419003   234.8606   242.8118   272.9613   309.7688   317.4629   309.9976   454.1619   346.8311
  Reaction.8 Reaction.9
1   430.5853   466.3535
2   224.2957   237.3142
3   261.0125   247.5153
4   305.3495   354.0487
5         NA         NA
6   330.3003   253.8644

We’ll use the reshape() function in the reshape package to convert from wide to long format. Sidenote, there’s lot of packages to reshape but I have always had the best success with the reshape() function.

head(sleep.long)

      SubNum female      gpa time reaction
101.1    101      0 3.090871    0 249.5600
101.2    101      0 3.090871    1 258.7047
101.3    101      0 3.090871    2 250.8006
101.4    101      0 3.090871    3 321.4398
101.5    101      0 3.090871    4 356.8519
101.6    101      0 3.090871    5 414.6901

Ordinarily, we might use ggplot2 to create facet plots to view individual trajectories but because we’re trying to stick with the base graphics, we’ll just use plot() and write a loop. First, let’s plot the mean growth over time.

mean.react <- tapply(sleep.long$reaction, sleep.long$time, FUN = mean, na.rm = T)
react.data <- data.frame(react = mean.react, time = 0:9)
plot(mean.react ~ time, react.data, type = "b", xlab = "Time", ylab = "Reaction Time (Mean)")

Again, we could pretty up this plot if we wanted to for publication but this works fine for us. In general, it appears that we a linear growth curve could probably do a reasonably decent job.

Now let’s do a trellis plot to view individual trajectories.

We should also look at a plot which juxtaposes these trajectories.

It looks like there are both variations in where people start and their growth trajectories over time. If we consider a model where we only allow for variations in starting places, this is known as a random intercept model. If we also allow for variations in growth over time this is called random slopes.

Fitting the growth model in lavaan

The syntax is very similar in to fitting a CFA or SEM. As was mentioned in lecture, we can consider the random intercepts and random slopes, collectively referred to as a random effects, as continuous latent variables with a mean of 0 and some non-zero variance.

library("lavaan")

This is lavaan 0.5-20
lavaan is BETA software! Please report any bugs.

library("semPlot")
growth.model <- '
  i =~ 1*Reaction.0 +  1*Reaction.1 +  1*Reaction.2 +  1*Reaction.3 +  1*Reaction.4 +  1*Reaction.5 +  1*Reaction.6 +  1*Reaction.7 +  1*Reaction.8 +  1*Reaction.9
  s =~ 0*Reaction.0 +  1*Reaction.1 +  2*Reaction.2 +  3*Reaction.3 +  4*Reaction.4 +  5*Reaction.5 +  6*Reaction.6 +  7*Reaction.7 +  8*Reaction.8 +  9*Reaction.9
  i ~~ s
'

The path diagram is:

The model is fitted and summary information is produced. Note, this time we are using the growth() function.

summary(fitted.model)

lavaan (0.5-20) converged normally after  89 iterations

  Number of observations                            16

  Estimator                                         ML
  Minimum Function Test Statistic              155.467
  Degrees of freedom                                59
  P-value (Chi-square)                           0.000

Parameter Estimates:

  Information                                 Expected
  Standard Errors                             Standard

Latent Variables:
                   Estimate  Std.Err  Z-value  P(>|z|)
  i =~                                                
    Reaction.0        1.000                           
    Reaction.1        1.000                           
    Reaction.2        1.000                           
    Reaction.3        1.000                           
    Reaction.4        1.000                           
    Reaction.5        1.000                           
    Reaction.6        1.000                           
    Reaction.7        1.000                           
    Reaction.8        1.000                           
    Reaction.9        1.000                           
  s =~                                                
    Reaction.0        0.000                           
    Reaction.1        1.000                           
    Reaction.2        2.000                           
    Reaction.3        3.000                           
    Reaction.4        4.000                           
    Reaction.5        5.000                           
    Reaction.6        6.000                           
    Reaction.7        7.000                           
    Reaction.8        8.000                           
    Reaction.9        9.000                           

Covariances:
                   Estimate  Std.Err  Z-value  P(>|z|)
  i ~~                                                
    s                39.232   42.643    0.920    0.358

Intercepts:
                   Estimate  Std.Err  Z-value  P(>|z|)
    Reaction.0        0.000                           
    Reaction.1        0.000                           
    Reaction.2        0.000                           
    Reaction.3        0.000                           
    Reaction.4        0.000                           
    Reaction.5        0.000                           
    Reaction.6        0.000                           
    Reaction.7        0.000                           
    Reaction.8        0.000                           
    Reaction.9        0.000                           
    i               251.819    6.671   37.748    0.000
    s                10.318    1.589    6.494    0.000

Variances:
                   Estimate  Std.Err  Z-value  P(>|z|)
    Reactin.0 (Re)  666.144   83.268    8.000    0.000
    Reactin.1 (Re)  666.144   83.268    8.000    0.000
    Reactin.2 (Re)  666.144   83.268    8.000    0.000
    Reactin.3 (Re)  666.144   83.268    8.000    0.000
    Reactin.4 (Re)  666.144   83.268    8.000    0.000
    Reactin.5 (Re)  666.144   83.268    8.000    0.000
    Reactin.6 (Re)  666.144   83.268    8.000    0.000
    Reactin.7 (Re)  666.144   83.268    8.000    0.000
    Reactin.8 (Re)  666.144   83.268    8.000    0.000
    Reactin.9 (Re)  666.144   83.268    8.000    0.000
    i               481.933  253.388    1.902    0.057
    s                32.311   14.314    2.257    0.024

The path diagram for the fitted model is:

A similar, though not identical model, can be fit in nlme::lme with the following syntax.

summary(growth.lme)

Linear mixed-effects model fit by REML
 Data: sleep.long 
     AIC      BIC  logLik
  1561.5 1579.875 -774.75

Random effects:
 Formula: ~1 + time | SubNum
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 23.008054 (Intr)
time         5.916318 0.29  
Residual    25.809856       

Fixed effects: reaction ~ time 
                Value Std.Error  DF  t-value p-value
(Intercept) 251.81899  6.889733 143 36.54989       0
time         10.31806  1.640833 143  6.28831       0
 Correlation: 
     (Intr)
time 0.017 

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-3.903444378 -0.456212342  0.002653454  0.436323779  5.133051327 

Number of Observations: 160
Number of Groups: 16 

Note that the estimates are similar but not identical and that our estimates from lavaan are actually included in teh confidence intervals.

If we want a fitted plot fitted, we need to extract the factor scores and do a little math. First, let’s combine the data and just look at the correlation between our predicted scores and our observed scores.

# View the data
head(predict.data)

   reaction.pred SubNum time reaction.orig
1       246.1592    101    0      249.5600
17      269.2409    101    1      258.7047
33      292.3226    101    2      250.8006
49      315.4043    101    3      321.4398
65      338.4860    101    4      356.8519
81      361.5676    101    5      414.6901

# Correlation
cor(predict.data$reaction.pred, predict.data$reaction.orig)

[1] 0.9103764

Let’s create a new trellis plot to see how our model looks.

The model fits well for some people but not well for everyone. Finally, our overall mean growth curve.

Overall, the model seems to fit reasonably well.