The Cortisol Data

Growth Models: Simple to Complex

Overview

This script shares The Cortisol Data, an N = 34, T = 9 time points data set we have used to illustrate a variety of growth modeling and mixture modeling methods. We describe the data and then walk through R-based implementations of the models covered in:

         Ram, N., & Grimm, K. (2007). Using simple and complex growth models to
         articulate developmental change: Matching theory to method.
         International Journal of Behavioral Development, 31(4), 303-316.
         https://doi.org/10.1177/0165025407077751
        .

The data are shared with intent that others may find them useful for learning about growth modeling or developing new methods for analysis of change.

New publications based on these data require citation and acknowledgement of the full set of papers that have used the data, inlcuding the above paper, and …

         Ram, N., & Grimm, K. (2009). Growth mixture modeling: A method for
         identifying differences in longitudinal change among unobserved groups.
         International Journal of Behavioral Development, 33(6), 565-576.
         https://doi.org/10.1177/0165025409343765

         Ram, N., Grimm, K., Gatzke-Kopp, L. & Molenaar, P.C.M. (2011).
         Longitudinal mixture models and the identification of archetypes:
         Action-adventure, mystery, science fiction, fantasy, or romance? In B.
         Laursen, T. Little, & N. Card (Eds.) Handbook of Developmental Research
         Methods (pp. 481-500). New York: Guilford.

         Grimm, K.J., Steele, J.S., Ram, N., & Nesselroade, J.R. (2013).
         Exploratory latent growth models in the structural equation modeling
         framework. Structural Equation Modeling, 20(4), 568-591.
         https://doi.org/10.1080/10705511.2013.824775
        .

Nilam & Kevin

Outline

This script covers …

Loading the Cortisol Data
Reshaping the Cortisol Data
Fitting the series of growth models from Ram & Grimm (2007):
- Linear
- Quadratic
- Latent Basis
- Exponential
- Multiphase

Preliminaries

Loading libraries used in this script.

library(psych)           # data descriptives
library(nlme)            # Linear & Nonlinear Mixed Effects Models
library(lme4)            # Linear Mixed-Effects Models using 'Eigen' and S4
library(lmerTest)        # Model p-values
library(tidyverse)       # ggplot, stringr, dplyr, readr, forcats, tidyr, purr

Loading the public data

#set filepath for data file
filepath <- "https://raw.githubusercontent.com/The-Change-Lab/collaborations/main/TheCortisolData/TheCortisolData.csv"

#read in the .csv file using the url() function
cortisol_wide <- read.csv(file=url(filepath), header=TRUE)

Looking at the top few rows of the wide data.

head(cortisol_wide, 6)

##   id cort_0 cort_1 cort_2 cort_3 cort_4 cort_5 cort_6 cort_7 cort_8
## 1  1    4.2    4.1    9.7   14.0   19.0   18.0   20.0   23.0   24.0
## 2  2    5.5    5.6   14.0   16.0   19.0   17.0   18.0   20.0   19.0
## 3  3    4.0    3.8    7.5   12.0   14.0   13.0    9.1    8.2    7.9
## 4  4    6.1    5.6   14.0   20.0   26.0   23.0   26.0   25.0   26.0
## 5  5    4.6    4.4    7.2   12.3   15.8   16.1   17.0   17.8   19.1
## 6  6    6.8    9.5   14.2   19.6   19.0   13.9   13.4   12.5   11.7

Reshaping the data

Two main data schema are used to accommodate repeated measures data - “Wide Format” and “Long Format”. Different functions work with different kinds of data input. We already have the wide format data. We make a set of long format data.

Reshape from wide to long

#reshaping wide to long
cortisol_long <- reshape(data=cortisol_wide, 
                         timevar=c("time"), 
                         idvar="id",
                         varying=c("cort_0","cort_1","cort_2","cort_3",
                                   "cort_4","cort_5","cort_6","cort_7","cort_8"),
                         direction="long", sep="_")

#sorting for easy viewing
#order by id and time
cortisol_long <- cortisol_long[order(cortisol_long$id,cortisol_long$time), ]

To match the scaling of time used in some of the papers we add an additional time variable that runs from 0 to 1.

cortisol_long$timescaled <- (cortisol_long$time - 0)/8

Looking at the top few rows of the long data.

head(cortisol_long, 18)

##     id time cort timescaled
## 1.0  1    0  4.2      0.000
## 1.1  1    1  4.1      0.125
## 1.2  1    2  9.7      0.250
## 1.3  1    3 14.0      0.375
## 1.4  1    4 19.0      0.500
## 1.5  1    5 18.0      0.625
## 1.6  1    6 20.0      0.750
## 1.7  1    7 23.0      0.875
## 1.8  1    8 24.0      1.000
## 2.0  2    0  5.5      0.000
## 2.1  2    1  5.6      0.125
## 2.2  2    2 14.0      0.250
## 2.3  2    3 16.0      0.375
## 2.4  2    4 19.0      0.500
## 2.5  2    5 17.0      0.625
## 2.6  2    6 18.0      0.750
## 2.7  2    7 20.0      0.875
## 2.8  2    8 19.0      1.000

Fitting the series of growth models from Ram & Grimm (2007)

Linear Growth Model

#linear model
cort_linear <- nlme(cort ~ g0 + g1*timescaled,
                    fixed = g0 + g1 ~ 1,
                    random = g0 + g1 ~ 1,
                    group = ~id,
                    start = c(g0=8.0, g1=10.9),
                    data = cortisol_long,
                    na.action = "na.exclude")

## Warning in nlme.formula(cort ~ g0 + g1 * timescaled, fixed = g0 + g1 ~ 1, :
## Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase
## 'msMaxIter'!

summary(cort_linear)

## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: cort ~ g0 + g1 * timescaled 
##   Data: cortisol_long 
##        AIC      BIC    logLik
##   1837.972 1860.313 -912.9859
## 
## Random effects:
##  Formula: list(g0 ~ 1, g1 ~ 1)
##  Level: id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev    Corr 
## g0       0.2031221 g0   
## g1       4.4910644 0.999
## Residual 4.3824247      
## 
## Fixed effects:  g0 + g1 ~ 1 
##        Value Std.Error  DF   t-value p-value
## g0  8.000588 0.4647811 271 17.213670       0
## g1 10.881176 1.0970972 271  9.918152       0
##  Correlation: 
##    g0    
## g1 -0.542
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.8427253 -0.7510109 -0.1820519  0.5969253  2.9508200 
## 
## Number of Observations: 306
## Number of Groups: 34

VarCorr(cort_linear)

## id = pdLogChol(list(g0 ~ 1,g1 ~ 1)) 
##          Variance    StdDev    Corr 
## g0        0.04125857 0.2031221 g0   
## g1       20.16965946 4.4910644 0.999
## Residual 19.20564640 4.3824247

Plotting the predicted trajectories

#obtaining predicted scores for individuals
cortisol_long$pred_linear <- predict(cort_linear)

#obtaining predicted scores for prototype
cortisol_long$proto_linear <- predict(cort_linear, level=0)

#plotting predicted trajectories
#intraindividual change trajetories
cortisol_long %>%
  ggplot( aes(x = time, y = pred_linear, group = id)) +
  geom_line(color="black") +
  geom_line(aes(x = time, y = proto_linear), 
            color="red", linewidth = 1.5) + 
  xlab("Time") + 
  ylab("Cortisol") + ylim(0,30) +
  scale_x_continuous(breaks=seq(0,8, by = 1)) +
  theme_minimal()

Note that the model has some convergence issues, and the solution has hit a parameter boundary. In this model, and some of the models that follow, the correlation between the random effects in intercept and slope is questionable.

Quadratic Growth Model

#creating the quadratic time variable
cortisol_long$timesq <- cortisol_long$timescaled^2

cort_quad <- lmer(cort ~ 1 + timescaled + timesq + 
                    (1 + timescaled + timesq | id),
                    data = cortisol_long,
                    na.action = "na.exclude")
summary(cort_quad)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: cort ~ 1 + timescaled + timesq + (1 + timescaled + timesq | id)
##    Data: cortisol_long
## 
## REML criterion at convergence: 1649.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.24859 -0.66952  0.03488  0.52129  2.55586 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  id       (Intercept)  1.308   1.143               
##           timescaled  28.449   5.334     0.63      
##           timesq      26.582   5.156    -0.97 -0.44
##  Residual              9.353   3.058               
## Number of obs: 306, groups:  id, 34
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   3.6009     0.4692  51.9152   7.674 4.22e-10 ***
## timescaled   41.0508     2.1882  62.2810  18.760  < 2e-16 ***
## timesq      -30.1696     2.1071  51.2034 -14.318  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) tmscld
## timescaled -0.552       
## timesq      0.369 -0.870

cort_quad <- nlme(cort ~ g0 + g1*timescaled + g2*timesq,
                    fixed = g0 + g1 + g2 ~ 1,
                    random = g0 + g1 + g2 ~ 1,
                    group = ~id,
                    start = c(g0=3.6, g1=41.1, g2=-30.2),
                    data = cortisol_long,
                    na.action = "na.exclude")

## Warning in nlme.formula(cort ~ g0 + g1 * timescaled + g2 * timesq, fixed = g0 +
## : Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase
## 'msMaxIter'!

summary(cort_quad)

## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: cort ~ g0 + g1 * timescaled + g2 * timesq 
##   Data: cortisol_long 
##        AIC     BIC    logLik
##   1675.114 1712.35 -827.5571
## 
## Random effects:
##  Formula: list(g0 ~ 1, g1 ~ 1, g2 ~ 1)
##  Level: id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev   Corr         
## g0       1.116417 g0     g1    
## g1       5.246087  0.627       
## g2       5.033705 -0.964 -0.436
## Residual 3.051628              
## 
## Fixed effects:  g0 + g1 + g2 ~ 1 
##        Value Std.Error  DF    t-value p-value
## g0   3.60086 0.4687746 270   7.681422       0
## g1  41.05077 2.1888421 270  18.754561       0
## g2 -30.16959 2.1050257 270 -14.332174       0
##  Correlation: 
##    g0     g1    
## g1 -0.560       
## g2  0.381 -0.873
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.24520649 -0.67204216  0.03085329  0.53148561  2.57505404 
## 
## Number of Observations: 306
## Number of Groups: 34

VarCorr(cort_quad)

## id = pdLogChol(list(g0 ~ 1,g1 ~ 1,g2 ~ 1)) 
##          Variance  StdDev   Corr         
## g0        1.246386 1.116417 g0     g1    
## g1       27.521432 5.246087  0.627       
## g2       25.338182 5.033705 -0.964 -0.436
## Residual  9.312435 3.051628

Plotting the predicted trajectories

#obtaining predicted scores for individuals
cortisol_long$pred_quad <- predict(cort_quad)

#obtaining predicted scores for prototype
cortisol_long$proto_quad <- predict(cort_quad, level=0)

#plotting predicted trajectories
#intraindividual change trajetories
cortisol_long %>%
  ggplot(aes(x = time, y = pred_quad, group = id)) +
  geom_line(color="black") +
  geom_line(aes(x = time, y = proto_quad), 
            color = "red", linewidth = 1.5) + 
  xlab("Time") + 
  ylab("Cortisol") + ylim(0,30) +
  scale_x_continuous(breaks=seq(0,8, by = 1)) +
  theme_minimal()

Latent Basis Growth Model

#creating time-dummy variables
cortisol_long$time0 <- ifelse(cortisol_long$time == 0, 1, 0)
cortisol_long$time1 <- ifelse(cortisol_long$time == 1, 1, 0)
cortisol_long$time2 <- ifelse(cortisol_long$time == 2, 1, 0)
cortisol_long$time3 <- ifelse(cortisol_long$time == 3, 1, 0)
cortisol_long$time4 <- ifelse(cortisol_long$time == 4, 1, 0)
cortisol_long$time5 <- ifelse(cortisol_long$time == 5, 1, 0)
cortisol_long$time6 <- ifelse(cortisol_long$time == 6, 1, 0)
cortisol_long$time7 <- ifelse(cortisol_long$time == 7, 1, 0)
cortisol_long$time8 <- ifelse(cortisol_long$time == 8, 1, 0)

#latent basis model
cort_latentb <- nlme(cort ~ time0*(g0 + g1*0) +
                            time1*(g0 + g1*A_1) +
                            time2*(g0 + g1*A_2) +
                            time3*(g0 + g1*A_3) +
                            time4*(g0 + g1*A_4) +
                            time5*(g0 + g1*A_5) +
                            time6*(g0 + g1*A_6) +
                            time7*(g0 + g1*A_7) +
                            time8*(g0 + g1*1),
                            fixed = g0 + g1 + A_1 + A_2 + A_3 + A_4 + A_5 + A_6 + A_7 ~ 1,
                            random = g0 + g1 ~ 1,
                         groups =~ id,
                         start = c(g0=5.5, g1=11.5, 
                                   A_1=.3, A_2=.4, A_3=.5, A_4=.6, 
                                   A_5=.7, A_6=.8, A_7=.9),
                         data = cortisol_long,
                         na.action = na.exclude) 

summary(cort_latentb)

## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: cort ~ time0 * (g0 + g1 * 0) + time1 * (g0 + g1 * A_1) + time2 *      (g0 + g1 * A_2) + time3 * (g0 + g1 * A_3) + time4 * (g0 +      g1 * A_4) + time5 * (g0 + g1 * A_5) + time6 * (g0 + g1 *      A_6) + time7 * (g0 + g1 * A_7) + time8 * (g0 + g1 * 1) 
##   Data: cortisol_long 
##        AIC      BIC    logLik
##   1569.164 1617.571 -771.5821
## 
## Random effects:
##  Formula: list(g0 ~ 1, g1 ~ 1)
##  Level: id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev   Corr 
## g0       1.698578 g0   
## g1       2.735103 0.001
## Residual 2.529224      
## 
## Fixed effects:  g0 + g1 + A_1 + A_2 + A_3 + A_4 + A_5 + A_6 + A_7 ~ 1 
##         Value Std.Error  DF   t-value p-value
## g0   5.342011 0.5207451 264 10.258399  0.0000
## g1  10.447225 0.7660285 264 13.638168  0.0000
## A_1 -0.041829 0.0582607 264 -0.717961  0.4734
## A_2  0.522586 0.0500425 264 10.442847  0.0000
## A_3  1.118113 0.0611387 264 18.288132  0.0000
## A_4  1.336652 0.0689221 264 19.393645  0.0000
## A_5  1.081029 0.0599665 264 18.027218  0.0000
## A_6  0.951968 0.0562971 264 16.909722  0.0000
## A_7  1.008691 0.0578269 264 17.443266  0.0000
##  Correlation: 
##     g0     g1     A_1    A_2    A_3    A_4    A_5    A_6   
## g1  -0.455                                                 
## A_1 -0.581  0.404                                          
## A_2 -0.314 -0.023  0.262                                   
## A_3  0.056 -0.446 -0.066  0.287                            
## A_4  0.151 -0.534 -0.150  0.261  0.589                     
## A_5  0.037 -0.427 -0.050  0.291  0.540  0.577              
## A_6 -0.034 -0.355  0.013  0.306  0.501  0.526  0.495       
## A_7 -0.002 -0.388 -0.015  0.300  0.519  0.550  0.512  0.481
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.4004317 -0.6084323 -0.1077279  0.5806294  2.8252173 
## 
## Number of Observations: 306
## Number of Groups: 34

VarCorr(cort_latentb)

## id = pdLogChol(list(g0 ~ 1,g1 ~ 1)) 
##          Variance StdDev   Corr 
## g0       2.885168 1.698578 g0   
## g1       7.480787 2.735103 0.001
## Residual 6.396974 2.529224

Plotting the predicted trajectories

#obtaining predicted scores for individuals
cortisol_long$pred_latentb <- predict(cort_latentb)

#obtaining predicted scores for prototype
cortisol_long$proto_latentb <- predict(cort_latentb, level=0)

#plotting predicted trajectories
#intraindividual change trajetories
cortisol_long %>%
  ggplot(aes(x = time, y = pred_latentb, group = id)) +
  geom_line(color="black") +
  geom_line(aes(x = time, y = proto_latentb), 
            color = "red", linewidth = 1.5) + 
  xlab("Time") + 
  ylab("Cortisol") + ylim(0,30) +
  scale_x_continuous(breaks=seq(0,8, by = 1)) +
  theme_minimal()

Exponential Growth Model

#exponential model
cort_expo <- nlme(cort ~ g0 + g1*(exp(-1*alpha*timescaled)),
                    fixed = g0 + g1 + alpha ~ 1,
                    random = g0 + g1  ~ 1,
                    group = ~id,
                    start = c(g0=17, g1=-14, alpha=0.5),
                    data = cortisol_long,
                    na.action = "na.exclude")

## Warning in nlme.formula(cort ~ g0 + g1 * (exp(-1 * alpha * timescaled)), :
## Iteration 4, LME step: nlminb() did not converge (code = 1). Do increase
## 'msMaxIter'!

summary(cort_expo)

## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: cort ~ g0 + g1 * (exp(-1 * alpha * timescaled)) 
##   Data: cortisol_long 
##        AIC      BIC    logLik
##   1730.236 1756.302 -858.1182
## 
## Random effects:
##  Formula: list(g0 ~ 1, g1 ~ 1)
##  Level: id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev   Corr  
## g0       3.672342 g0    
## g1       3.988619 -0.999
## Residual 3.594282       
## 
## Fixed effects:  g0 + g1 + alpha ~ 1 
##            Value Std.Error  DF    t-value p-value
## g0     17.203250 0.7689413 270  22.372645       0
## g1    -13.735218 0.9554844 270 -14.375136       0
## alpha   4.108994 0.4851947 270   8.468752       0
##  Correlation: 
##       g0     g1    
## g1    -0.788       
## alpha -0.442  0.077
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.06657120 -0.69580376 -0.03326631  0.60565773  2.73037715 
## 
## Number of Observations: 306
## Number of Groups: 34

VarCorr(cort_expo)

## id = pdLogChol(list(g0 ~ 1,g1 ~ 1)) 
##          Variance StdDev   Corr  
## g0       13.48610 3.672342 g0    
## g1       15.90908 3.988619 -0.999
## Residual 12.91886 3.594282

Plotting the predicted trajectories

#obtaining predicted scores for individuals
cortisol_long$pred_expo <- predict(cort_expo)

#obtaining predicted scores for prototype
cortisol_long$proto_expo <- predict(cort_expo, level=0)

#plotting predicted trajectories
#intraindividual change trajetories
cortisol_long %>%
  ggplot(aes(x = time, y = pred_expo, group = id)) +
  geom_line(color="black") +
  geom_line(aes(x = time, y = proto_expo), 
            color="red", linewidth = 1.5) + 
  xlab("Time") + 
  ylab("Cortisol") + ylim(0,30) +
  scale_x_continuous(breaks=seq(0,8, by = 1)) +
  theme_minimal()

Multiphase Growth Model

#multiphase model
cort_multi <- nlme(cort ~ time0*(g0 + g1*0 + g2*0) +
                          time1*(g0 + g1*0 + g2*0) +
                          time2*(g0 + g1*A_2 + g2*0) +
                          time3*(g0 + g1*A_3 + g2*0) +
                          time4*(g0 + g1*1 + g2*0) +
                          time5*(g0 + g1*1 + g2*A_5) +
                          time6*(g0 + g1*1 + g2*A_6) +
                          time7*(g0 + g1*1 + g2*A_7) +
                          time8*(g0 + g1*1 + g2*1),
                   fixed = g0 + g1 + g2 + A_2 + A_3 + A_5 + A_6 + A_7 ~ 1,
                   random = g0 + g1 + g2 ~ 1,
                   groups =~ id,
                   start = c(g0=15, g1=10, g2=-4, 
                             A_2=.4, A_3=.5, A_5=.7, 
                             A_6=.8, A_7=.9),
                   data = cortisol_long,
                   na.action = na.exclude) 

summary(cort_multi)

## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: cort ~ time0 * (g0 + g1 * 0 + g2 * 0) + time1 * (g0 + g1 * 0 +      g2 * 0) + time2 * (g0 + g1 * A_2 + g2 * 0) + time3 * (g0 +      g1 * A_3 + g2 * 0) + time4 * (g0 + g1 * 1 + g2 * 0) + time5 *      (g0 + g1 * 1 + g2 * A_5) + time6 * (g0 + g1 * 1 + g2 * A_6) +      time7 * (g0 + g1 * 1 + g2 * A_7) + time8 * (g0 + g1 * 1 +      g2 * 1) 
##   Data: cortisol_long 
##        AIC      BIC    logLik
##   1412.881 1468.735 -691.4407
## 
## Random effects:
##  Formula: list(g0 ~ 1, g1 ~ 1, g2 ~ 1)
##  Level: id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev   Corr         
## g0       2.161867 g0     g1    
## g1       3.580496 -0.315       
## g2       4.241375 -0.041 -0.286
## Residual 1.535205              
## 
## Fixed effects:  g0 + g1 + g2 + A_2 + A_3 + A_5 + A_6 + A_7 ~ 1 
##         Value Std.Error  DF  t-value p-value
## g0   5.118082 0.4201167 265 12.18252       0
## g1  14.248891 0.6998856 265 20.35889       0
## g2  -3.797260 0.8089645 265 -4.69398       0
## A_2  0.410745 0.0210662 265 19.49783       0
## A_3  0.860230 0.0238885 265 36.01023       0
## A_5  0.715784 0.0598000 265 11.96963       0
## A_6  1.046663 0.0683974 265 15.30267       0
## A_7  1.042029 0.0682431 265 15.26937       0
##  Correlation: 
##     g0     g1     g2     A_2    A_3    A_5    A_6   
## g1  -0.372                                          
## g2  -0.033 -0.352                                   
## A_2 -0.158 -0.032  0.111                            
## A_3 -0.033 -0.214  0.205  0.250                     
## A_5  0.000  0.059  0.073 -0.056 -0.104              
## A_6  0.000 -0.008  0.165  0.007  0.014  0.393       
## A_7  0.000 -0.007  0.164  0.006  0.012  0.392  0.504
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.40843238 -0.50779682 -0.01780372  0.49095891  3.70249846 
## 
## Number of Observations: 306
## Number of Groups: 34

VarCorr(cort_multi)

## id = pdLogChol(list(g0 ~ 1,g1 ~ 1,g2 ~ 1)) 
##          Variance  StdDev   Corr         
## g0        4.673669 2.161867 g0     g1    
## g1       12.819948 3.580496 -0.315       
## g2       17.989264 4.241375 -0.041 -0.286
## Residual  2.356854 1.535205

Plotting the predicted trajectories

#obtaining predicted scores for individuals
cortisol_long$pred_multi <- predict(cort_multi)

#obtaining predicted scores for prototype
cortisol_long$proto_multi <- predict(cort_multi, level=0)

#plotting predicted trajectories
#intraindividual change trajetories
cortisol_long %>%
  ggplot(aes(x = time, y = pred_multi, group = id)) +
  geom_line(color="black") +
  geom_line(aes(x = time, y = proto_multi), 
            color="red", linewidth = 1.5) + 
  xlab("Time") + 
  ylab("Cortisol") + ylim(0,30) +
  scale_x_continuous(breaks=seq(0,8, by = 1)) +
  theme_minimal()

These were the 5 models covered in Ram & Grimm (2007). In sharing the data, we encourage the community to discover even more interesting models.

Thanks for playing!

Nilam & Kevin