The Cortisol Data - Growth Models: Simple to Complex
0.1 Overview
This script shares The Cortisol Data, an N = 34, T = 9 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, 303-316.
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, 565-576.
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, 568-591.
Thank you to John Nesselroade, Teresa Seeman, and Marilyn Albert for data collection and inspiration.
Please contact us with questions about the data and their use.
Nilam & Kevin
0.1.1 Outline
This script covers …
A. Describing The Cortisol Data
B. Fitting the series of growth models from Ram & Grimm (2007): Linear, Quadratic, Latent Basis, Exponential, Multiphase
0.1.1.1 Preliminaries
Loading libraries used in this script.
library(psych)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.4.4library(nlme)
library(lme4)## Warning: package 'lme4' was built under R version 3.4.40.1.1.2 Reading in the data
Loading the public data
#set filepath for data file
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/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 | 4.2 | 4.1 | 9.7 | 14.0 | 19.0 | 18.0 | 20.0 | 23.0 | 24.0 |
| 2 | 5.5 | 5.6 | 14.0 | 16.0 | 19.0 | 17.0 | 18.0 | 20.0 | 19.0 |
| 3 | 4.0 | 3.8 | 7.5 | 12.0 | 14.0 | 13.0 | 9.1 | 8.2 | 7.9 |
| 4 | 6.1 | 5.6 | 14.0 | 20.0 | 26.0 | 23.0 | 26.0 | 25.0 | 26.0 |
| 5 | 4.6 | 4.4 | 7.2 | 12.3 | 15.8 | 16.1 | 17.0 | 17.8 | 19.1 |
| 6 | 6.8 | 9.5 | 14.2 | 19.6 | 19.0 | 13.9 | 13.4 | 12.5 | 11.7 |
0.1.1.3 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)/8Looking 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 |
0.2 A. Describing The Cortisol Data
Basic descriptives of the 9-occasion data.
#means, sd, etc.
describe(cortisol_wide)| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| id | 1 | 34 | 17.500000 | 9.958246 | 17.50 | 17.500000 | 12.60210 | 1 | 34.0 | 33.0 | 0.0000000 | -1.3062827 | 1.7078251 |
| cort_0 | 2 | 34 | 5.247059 | 2.454462 | 5.00 | 4.903571 | 2.29803 | 2 | 12.0 | 10.0 | 1.1784853 | 0.9152094 | 0.4209368 |
| cort_1 | 3 | 34 | 5.064706 | 2.593685 | 4.50 | 4.782143 | 2.22390 | 2 | 11.3 | 9.3 | 1.0159850 | -0.0131327 | 0.4448133 |
| cort_2 | 4 | 34 | 10.847059 | 3.110124 | 10.05 | 10.789286 | 3.85476 | 6 | 16.7 | 10.7 | 0.0864489 | -1.1724863 | 0.5333818 |
| cort_3 | 5 | 34 | 17.355882 | 3.685614 | 17.35 | 17.392857 | 3.48411 | 7 | 24.9 | 17.9 | -0.2934763 | 0.3862287 | 0.6320776 |
| cort_4 | 6 | 34 | 19.479412 | 3.722314 | 19.10 | 19.492857 | 2.81694 | 12 | 27.1 | 15.1 | -0.0439142 | -0.5599728 | 0.6383716 |
| cort_5 | 7 | 34 | 16.461765 | 4.180838 | 16.95 | 16.414286 | 4.52193 | 8 | 26.0 | 18.0 | 0.0779079 | -0.5824582 | 0.7170077 |
| cort_6 | 8 | 34 | 15.029412 | 4.569366 | 14.65 | 14.900000 | 4.15128 | 6 | 26.0 | 20.0 | 0.2614218 | -0.3075741 | 0.7836399 |
| cort_7 | 9 | 34 | 15.620588 | 4.809446 | 16.15 | 15.603571 | 5.55975 | 7 | 25.0 | 18.0 | -0.0237796 | -0.9863200 | 0.8248133 |
| cort_8 | 10 | 34 | 15.864706 | 5.609455 | 16.00 | 15.814286 | 6.07866 | 6 | 26.0 | 20.0 | 0.1359885 | -1.0784551 | 0.9620137 |
#boxplot
#boxplot by day
ggplot(data=cortisol_long, aes(x=factor(time), y=cort)) +
geom_boxplot(notch = TRUE) +
stat_summary(fun.y="mean", geom="point", shape=23, size=3, fill="white") +
labs(x = "Time", y = "Cortisol")
#correlations
# Correlations
round(cor(cortisol_wide[,c("cort_0","cort_1","cort_2","cort_3",
"cort_4","cort_5","cort_6","cort_7","cort_8")],
use="complete.obs",method="spearman"),2)## cort_0 cort_1 cort_2 cort_3 cort_4 cort_5 cort_6 cort_7 cort_8
## cort_0 1.00 0.97 0.45 0.27 0.23 0.17 0.28 0.22 0.20
## cort_1 0.97 1.00 0.48 0.27 0.18 0.14 0.23 0.19 0.16
## cort_2 0.45 0.48 1.00 0.71 0.66 0.47 0.42 0.39 0.34
## cort_3 0.27 0.27 0.71 1.00 0.71 0.37 0.16 0.10 0.14
## cort_4 0.23 0.18 0.66 0.71 1.00 0.81 0.58 0.50 0.52
## cort_5 0.17 0.14 0.47 0.37 0.81 1.00 0.88 0.84 0.82
## cort_6 0.28 0.23 0.42 0.16 0.58 0.88 1.00 0.97 0.93
## cort_7 0.22 0.19 0.39 0.10 0.50 0.84 0.97 1.00 0.94
## cort_8 0.20 0.16 0.34 0.14 0.52 0.82 0.93 0.94 1.00#pairs plot from the psych library
pairs.panels(cortisol_wide[,c("cort_0","cort_1","cort_2","cort_3",
"cort_4","cort_5","cort_6","cort_7","cort_8"),])
Plot of sample-level change in distribution across time
#Density distribution by day
ggplot(data=cortisol_long, aes(x=cort)) +
geom_density(aes(group=factor(time), colour=factor(time), fill=factor(time)), alpha=0.3)
Plot of individual-level trajectories
#intraindividual change trajetories
ggplot(data = cortisol_long, aes(x = time, y = cort, group = id)) +
geom_point(color="black") +
geom_line(color="black") +
xlab("Time") +
ylab("Cortisol") + ylim(0,30) +
scale_x_continuous(breaks=seq(0,8,by=1)) 
Yay!
Now lets fit some growth models …
0.3 B. Fitting the series of growth models from Ram & Grimm (2007)
0.3.0.1 Linear Growth Model
# #linear model
# cort_linear <- lme(cort ~ 1 + time,
# random = ~ 1 + time | id,
# data = cortisol_long,
# na.action = "na.exclude")
# #does not converge
#
# cort_linear <- lmer(cort ~ 1 + time + (1 + time | id),
# data = cortisol_long,
# na.action = "na.exclude")
# #does not converge
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")
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.203126 g0
## g1 4.491069 0.999
## Residual 4.382425
##
## Fixed effects: g0 + g1 ~ 1
## Value Std.Error DF t-value p-value
## g0 8.000588 0.4647811 271 17.213668 0
## g1 10.881177 1.0970976 271 9.918148 0
## Correlation:
## g0
## g1 -0.542
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.8427249 -0.7510109 -0.1820521 0.5969255 2.9508202
##
## Number of Observations: 306
## Number of Groups: 34VarCorr(cort_linear)## id = pdLogChol(list(g0 ~ 1,g1 ~ 1))
## Variance StdDev Corr
## g0 0.04126016 0.203126 g0
## g1 20.16969646 4.491069 0.999
## Residual 19.20564559 4.382425Plotting 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
ggplot(data = cortisol_long, aes(x = time, y = pred_linear, group = id)) +
#geom_point(color="black") +
geom_line(color="black") +
geom_line(aes(x = time, y = proto_linear), color="red",size=2) +
xlab("Time") +
ylab("Cortisol") + ylim(0,30) +
scale_x_continuous(breaks=seq(0,8,by=1)) 
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.
0.3.0.2 Quadratic Growth Model
#creating the quadratic time variable
cortisol_long$timesq <- cortisol_long$timescaled^2
#quadratic model
# cort_quad <- lme(cort ~ 1 + timescaled + timesq,
# random = ~ 1 + timescaled + timesq| id,
# data = cortisol_long,
# na.action = "na.exclude")
# #does not converge
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 ['lmerMod']
## 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.24861 -0.66952 0.03489 0.52127 2.55593
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 1.307 1.143
## timescaled 28.429 5.332 0.63
## timesq 26.560 5.154 -0.97 -0.44
## Residual 9.353 3.058
## Number of obs: 306, groups: id, 34
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.6009 0.4692 7.674
## timescaled 41.0508 2.1881 18.761
## timesq -30.1696 2.1070 -14.319
##
## Correlation of Fixed Effects:
## (Intr) tmscld
## timescaled -0.552
## timesq 0.369 -0.870cort_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")
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.557
##
## Random effects:
## Formula: list(g0 ~ 1, g1 ~ 1, g2 ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## g0 1.115779 g0 g1
## g1 5.243092 0.628
## g2 5.029870 -0.967 -0.435
## Residual 3.051693
##
## Fixed effects: g0 + g1 + g2 ~ 1
## Value Std.Error DF t-value p-value
## g0 3.60086 0.4687378 270 7.682025 0
## g1 41.05077 2.1886675 270 18.756057 0
## g2 -30.16960 2.1047906 270 -14.333775 0
## Correlation:
## g0 g1
## g1 -0.560
## g2 0.380 -0.872
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.24527602 -0.67203180 0.03085987 0.53159569 2.57516257
##
## Number of Observations: 306
## Number of Groups: 34VarCorr(cort_quad)## id = pdLogChol(list(g0 ~ 1,g1 ~ 1,g2 ~ 1))
## Variance StdDev Corr
## g0 1.244964 1.115779 g0 g1
## g1 27.490017 5.243092 0.628
## g2 25.299597 5.029870 -0.967 -0.435
## Residual 9.312832 3.051693Plotting 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
ggplot(data = cortisol_long, aes(x = time, y = pred_quad, group = id)) +
#geom_point(color="black") +
geom_line(color="black") +
geom_line(aes(x = time, y = proto_quad), color="red",size=2) +
xlab("Time") +
ylab("Cortisol") + ylim(0,30) +
scale_x_continuous(breaks=seq(0,8,by=1)) 
0.3.0.3 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.027217 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: 34VarCorr(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.529224Plotting 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
ggplot(data = cortisol_long, aes(x = time, y = pred_latentb, group = id)) +
#geom_point(color="black") +
geom_line(color="black") +
geom_line(aes(x = time, y = proto_latentb), color="red",size=2) +
xlab("Time") +
ylab("Cortisol") + ylim(0,30) +
scale_x_continuous(breaks=seq(0,8,by=1)) 
0.3.0.4 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")
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.239 1756.304 -858.1195
##
## Random effects:
## Formula: list(g0 ~ 1, g1 ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## g0 3.672433 g0
## g1 3.988481 -0.999
## Residual 3.594264
##
## Fixed effects: g0 + g1 + alpha ~ 1
## Value Std.Error DF t-value p-value
## g0 17.202771 0.7689113 270 22.37289 0
## g1 -13.735113 0.9554561 270 -14.37545 0
## alpha 4.109679 0.4852488 270 8.46922 0
## Correlation:
## g0 g1
## g1 -0.788
## alpha -0.442 0.077
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.06666362 -0.69579849 -0.03313565 0.60577564 2.73039108
##
## Number of Observations: 306
## Number of Groups: 34VarCorr(cort_expo)## id = pdLogChol(list(g0 ~ 1,g1 ~ 1))
## Variance StdDev Corr
## g0 13.48676 3.672433 g0
## g1 15.90798 3.988481 -0.999
## Residual 12.91874 3.594264Plotting 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
ggplot(data = cortisol_long, aes(x = time, y = pred_expo, group = id)) +
#geom_point(color="black") +
geom_line(color="black") +
geom_line(aes(x = time, y = proto_expo), color="red",size=2) +
xlab("Time") +
ylab("Cortisol") + ylim(0,30) +
scale_x_continuous(breaks=seq(0,8,by=1)) 
0.3.0.5 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.580495 -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.50779684 -0.01780371 0.49095891 3.70249850
##
## Number of Observations: 306
## Number of Groups: 34VarCorr(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.580495 -0.315
## g2 17.989264 4.241375 -0.041 -0.286
## Residual 2.356854 1.535205Plotting 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
ggplot(data = cortisol_long, aes(x = time, y = pred_multi, group = id)) +
#geom_point(color="black") +
geom_line(color="black") +
geom_line(aes(x = time, y = proto_multi), color="red",size=2) +
xlab("Time") +
ylab("Cortisol") + ylim(0,30) +
scale_x_continuous(breaks=seq(0,8,by=1)) 
These were the 5 models covered in Ram & Grimm (2007). We will eventually add R-scripts for the models forwarded in the other papers that have used these data. In sharing the data, we encourage the community to discover even more interesting models.
Thanks for playing!
Nilam & Kevin