Age-related changes in the dynamics of fear-related regulation in early childhood
Overview
This script provides a brief overview of the data preparation, model fitting, and plotting for the main analyses presented in Morales et al. (2017). Individuals using this code are encouraged to read the published paper. For context, this is the abstract of the study:
Self-regulation is a dynamic process wherein executive processes (EP) delay, minimize or desist prepotent responses (PR) that arise in situations that threaten well-being. It is generally assumed that, over the course of early childhood, children expand and more effectively deploy their repertoire of EP-related strategies to regulate PR. However, longitudinal tests of these assumptions are scarce in part because self-regulation has been mostly studied as a static construct. This study engages dynamic systems modeling to examine developmental changes in self-regulation between ages 2 and 5 years. Second-by- second time-series data derived from behavioral observations of 112 children (63 boys) faced with novel laboratory-based situations designed to elicit wariness, hesitation, and fear were modeled using differential equation models designed to capture age-related changes in the intrinsic dynamics and bidirectional coupling of PR (fear/wariness) and EP (strategy use). Results revealed that dynamic models allow for the conceptualization and measurement of fear regulation as intrinsic processes as well as direct and indirect coupling between PR and EP. Several patterns of age-related changes were in line with developmental theory suggesting that PR weakened and was regulated more quickly and efficiently by EP at age 5 than at age 2. However, most findings were in the intrinsic dynamics and moderating influences between PR and EP rather than direct influences. The findings illustrate the precision with which specific aspects of self-regulation can be articulated using dynamic systems models, and how such models can be used to describe the development of self-regulation in nuanced and theoretically meaningful ways.
Please feel free to reach out (smoralespam@gmail.com), if you have any questions.
Citation:
Morales S, Ram N, Buss KA, Cole PM, Helm JL, Chow S‐M. Age‐related changes in the dynamics of fear‐related regulation in early childhood. Dev Sci. 2017; e12633. https://doi.org/10.1111/desc.12633
Outline
In this script, we will cover the set-up, creating epochs (horizontal aggregation), calculating the derivatives, fitting a coupled ordinary differential equations (ODE) model, and plotting the results.
Set-up: Loading packages and reading in the data
# Installing packages
list.of.packages <- c("nlme","psych", "zoo", "reshape", "car","ggplot2","grid","deSolve")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
# Loading packages
lapply(list.of.packages, require, character.only = TRUE)
# Reading in the data
df1 <- read.csv("Clown_age2_composites_4.16.2017.csv")
df2 <- read.csv("Lion_age5_composites_4.16.2017.csv")Descriptives of the raw second-by-second data.
describe(df1[,c("PR", "EP")])## vars n mean sd median trimmed mad min max range skew kurtosis se
## PR 1 31744 0.35 1.05 0 0.05 0.00 0 9 9 3.99 18.94 0.01
## EP 2 31744 1.14 0.99 1 1.04 1.48 0 4 4 0.37 -0.92 0.01describe(df2[,c("PR", "EP")])## vars n mean sd median trimmed mad min max range skew kurtosis se
## PR 1 21769 0.10 0.45 0 0.00 0 0 6 6 5.04 27.1 0
## EP 2 21769 1.04 0.64 1 1.02 0 0 4 4 0.38 0.7 0Creating 5-second epochs (horizontal aggregation) and calculating the derivatives
# For Age 2
df_data_list = split(df1[,c('id', 'EP','PR')],, f = factor(df1$id))
ma_EP_PR_2 = list()
for(x in 1:length(df_data_list)){
curr_data = df_data_list[[x]][,c('id', 'PR', 'EP')]
curr_id = unique(curr_data$id)
curr_PR = curr_data$PR
curr_EP = curr_data$EP
curr_remainder_PR = length(curr_PR) %% 5
curr_remainder_EP = length(curr_EP) %% 5
if(curr_remainder_PR > 0){
curr_PR_2 = c(curr_PR, rep(NA, curr_remainder_PR))
}
if(curr_remainder_PR == 0){
curr_PR_2 = curr_PR
}
if(curr_remainder_EP > 0){
curr_EP_2 = c(curr_EP, rep(NA, curr_remainder_EP))
}
if(curr_remainder_EP == 0){
curr_EP_2 = curr_EP
}
windowed_EP = rollapply(curr_EP_2, width = 5, by = 5, FUN = mean, align = 'left')
windowed_PR = rollapply(curr_PR_2, width = 5, by = 5, FUN = mean, align = 'left')
ma_EP_PR_2[[x]] = cbind(curr_id, windowed_EP, windowed_PR)
}
#-----------------------------------------------------------#
# Calculate derivatives for each of these #
GLLA <- function(x, embed, tau, deltaT, order=2) {
L <- rep(1,embed)
for(i in 1:order) {
L <- cbind(L,(((c(1:embed)-mean(1:embed))*tau*deltaT)^i)/factorial(i))
}
W <- solve(t(L)%*%L)%*%t(L)
EMat <- Embed(x,embed,1)
Estimates <- EMat%*%t(W)
return(Estimates)
}
Embed <- function(x,E,tau) {
len <- length(x)
out <- x[1:(len-(E*tau)+tau)]
for(i in 2:E) { out <- cbind(out,x[(1+((i-1)*tau)):(len-(E*tau)+(i*tau))]) }
return(out)
}
deriv_df=NULL
for(x in 1:length(ma_EP_PR_2)){
curr_data = ma_EP_PR_2[[x]]
curr_id = unique(curr_data[,'curr_id'])
curr_EP_derivs = GLLA(curr_data[,'windowed_EP'], embed = 4, tau = 1, deltaT = 1, order = 2)
curr_PR_derivs = GLLA(curr_data[,'windowed_PR'], embed = 4, tau = 1, deltaT = 1, order = 2)
curr_EP_derivs = data.frame(curr_EP_derivs)
curr_PR_derivs = data.frame(curr_PR_derivs)
names(curr_EP_derivs) = c('EP_0th', 'EP_1st', 'EP_2nd')
names(curr_PR_derivs) = c('PR_0th', 'PR_1st', 'PR_2nd')
deriv_df[[x]] = cbind('id' = curr_id, 'time' = 1:nrow(curr_EP_derivs), curr_EP_derivs, curr_PR_derivs)
}
derivs_df1_out = data.frame(do.call(rbind, deriv_df))
#*******Also, writing the windowed or epoched data
windowed_df_2 = data.frame(do.call(rbind, ma_EP_PR_2)) # For age 2
##################################################################
# For age 5
df_data_list = split(df2[,c('id', 'EP','PR')],, f = factor(df2$id))
ma_EP_PR_5 = list()
for(x in 1:length(df_data_list)){
curr_data = df_data_list[[x]][,c('id', 'PR', 'EP')]
curr_id = unique(curr_data$id)
curr_PR = curr_data$PR
curr_EP = curr_data$EP
curr_remainder_PR = length(curr_PR) %% 5
curr_remainder_EP = length(curr_EP) %% 5
if(curr_remainder_PR > 0){
curr_PR_2 = c(curr_PR, rep(NA, curr_remainder_PR))
}
if(curr_remainder_PR == 0){
curr_PR_2 = curr_PR
}
if(curr_remainder_EP > 0){
curr_EP_2 = c(curr_EP, rep(NA, curr_remainder_EP))
}
if(curr_remainder_EP == 0){
curr_EP_2 = curr_EP
}
windowed_EP = rollapply(curr_EP_2, width = 5, by = 5, FUN = mean, align = 'left')
windowed_PR = rollapply(curr_PR_2, width = 5, by = 5, FUN = mean, align = 'left')
ma_EP_PR_5[[x]] = cbind(curr_id, windowed_EP, windowed_PR)
}
#-----------------------------------------------------------#
# Calculate derivatives for each of these #
GLLA <- function(x, embed, tau, deltaT, order=2) {
L <- rep(1,embed)
for(i in 1:order) {
L <- cbind(L,(((c(1:embed)-mean(1:embed))*tau*deltaT)^i)/factorial(i))
}
W <- solve(t(L)%*%L)%*%t(L)
EMat <- Embed(x,embed,1)
Estimates <- EMat%*%t(W)
return(Estimates)
}
Embed <- function(x,E,tau) {
len <- length(x)
out <- x[1:(len-(E*tau)+tau)]
for(i in 2:E) { out <- cbind(out,x[(1+((i-1)*tau)):(len-(E*tau)+(i*tau))]) }
return(out)
}
deriv_df=NULL
for(x in 1:length(ma_EP_PR_5)){
curr_data = ma_EP_PR_5[[x]]
curr_id = unique(curr_data[,'curr_id'])
curr_EP_derivs = GLLA(curr_data[,'windowed_EP'], embed = 4, tau = 1, deltaT = 1, order = 2)
curr_PR_derivs = GLLA(curr_data[,'windowed_PR'], embed = 4, tau = 1, deltaT = 1, order = 2)
curr_EP_derivs = data.frame(curr_EP_derivs)
curr_PR_derivs = data.frame(curr_PR_derivs)
names(curr_EP_derivs) = c('EP_0th', 'EP_1st', 'EP_2nd')
names(curr_PR_derivs) = c('PR_0th', 'PR_1st', 'PR_2nd')
deriv_df[[x]] = cbind('id' = curr_id, 'time' = 1:nrow(curr_EP_derivs), curr_EP_derivs, curr_PR_derivs)
}
derivs_df2_out = data.frame(do.call(rbind, deriv_df))
#*******Also, writing the windowed or epoched data
windowed_df_5 = data.frame(do.call(rbind, ma_EP_PR_5)) # for age 5
# Descriptives of the windowed data
describe(windowed_df_2[-1])## vars n mean sd median trimmed mad min max range skew kurtosis se
## windowed_EP 1 6306 1.14 0.95 1 1.06 1.48 0 4 4 0.35 -0.91 0.01
## windowed_PR 2 6306 0.35 1.03 0 0.06 0.00 0 9 9 4.04 19.47 0.01describe(windowed_df_5[-1])## vars n mean sd median trimmed mad min max range skew kurtosis se
## windowed_EP 1 4325 1.04 0.53 1 1.04 0.3 0 3.4 3.4 0.13 0.9 0.01
## windowed_PR 2 4325 0.10 0.42 0 0.00 0.0 0 4.0 4.0 4.96 26.0 0.01Getting the derivatives into an analyzable format
# Creating a an age indicator
derivs_df1_out$age_ind <- 0
derivs_df2_out$age_ind <- 1
# Putting together both datasets
df <- rbind(derivs_df1_out,derivs_df2_out)
# Getting some descriptives
describeBy(df[-1], group=df$age_ind) ##
## Descriptive statistics by group
## group: 0
## vars n mean sd median trimmed mad min max range skew kurtosis se
## time 1 6020 28.50 17.23 28.00 27.78 20.76 1.00 87.00 86.00 0.34 -0.58 0.22
## EP_0th 2 5973 1.19 0.93 1.04 1.13 1.35 -0.22 4.00 4.22 0.29 -0.86 0.01
## EP_1st 3 5973 0.01 0.28 0.00 0.01 0.18 -1.20 1.16 2.36 0.17 1.52 0.00
## EP_2nd 4 5973 -0.01 0.39 0.00 0.00 0.30 -1.80 1.80 3.60 -0.03 2.22 0.00
## PR_0th 5 5973 0.35 1.04 0.00 0.07 0.00 -0.22 8.69 8.91 4.02 19.38 0.01
## PR_1st 6 5973 0.00 0.21 0.00 0.00 0.00 -2.54 2.40 4.94 0.45 21.60 0.00
## PR_2nd 7 5973 0.00 0.30 0.00 0.00 0.00 -3.00 2.30 5.30 -0.80 17.53 0.00
## age_ind 8 6020 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 NaN NaN 0.00
## -------------------------------------------------------------------------------------------------
## group: 1
## vars n mean sd median trimmed mad min max range skew kurtosis se
## time 1 4125 28.90 17.10 28.00 28.40 20.76 1.00 81.00 80.00 0.22 -0.82 0.27
## EP_0th 2 4100 1.07 0.48 1.03 1.07 0.33 -0.14 3.32 3.46 0.10 1.13 0.01
## EP_1st 3 4100 0.00 0.20 0.00 0.00 0.18 -1.14 0.72 1.86 -0.20 0.86 0.00
## EP_2nd 4 4100 0.00 0.33 0.00 0.00 0.30 -1.40 1.30 2.70 0.04 0.61 0.01
## PR_0th 5 4100 0.10 0.41 0.00 0.00 0.00 -0.19 3.56 3.75 4.87 25.15 0.01
## PR_1st 6 4100 0.00 0.14 0.00 0.00 0.00 -1.32 1.18 2.50 -0.67 24.72 0.00
## PR_2nd 7 4100 -0.01 0.22 0.00 0.00 0.00 -2.90 2.00 4.90 -2.78 39.98 0.00
## age_ind 8 4125 1.00 0.00 1.00 1.00 0.00 1.00 1.00 0.00 NaN NaN 0.00# Melting the data - to create double-entry file
dflong <- melt(df, id=c("id","time", "age_ind", "EP_1st", "PR_1st", "EP_0th", "PR_0th"))
# reordering by id and time
dflong <- dflong[order(dflong$id, dflong$age_ind, dflong$time, dflong$variable),]
#making binary indicators to use for turning coefficents on-off
dflong$PR <- ifelse(dflong$variable=="PR_2nd", 1, 0)
dflong$EP <- ifelse(dflong$variable=="EP_2nd", 1, 0)
head(dflong[dflong$id==28,], 20) # This is the same ID as shown in Fig 1 of the manuscript## id time age_ind EP_1st PR_1st EP_0th PR_0th variable value PR EP
## 1243 28 1 0 0.0 0.00 0.0000 0.0000 EP_2nd 0.0 0 1
## 11388 28 1 0 0.0 0.00 0.0000 0.0000 PR_2nd 0.0 1 0
## 1244 28 2 0 0.0 0.18 0.0000 -0.0375 EP_2nd 0.0 0 1
## 11389 28 2 0 0.0 0.18 0.0000 -0.0375 PR_2nd 0.3 1 0
## 1245 28 3 0 0.0 0.36 0.0000 0.2750 EP_2nd 0.0 0 1
## 11390 28 3 0 0.0 0.36 0.0000 0.2750 PR_2nd 0.2 1 0
## 1246 28 4 0 0.0 0.34 0.0000 0.8375 EP_2nd 0.0 0 1
## 11391 28 4 0 0.0 0.34 0.0000 0.8375 PR_2nd -0.3 1 0
## 1247 28 5 0 0.0 0.12 0.0000 1.0250 EP_2nd 0.0 0 1
## 11392 28 5 0 0.0 0.12 0.0000 1.0250 PR_2nd -0.2 1 0
## 1248 28 6 0 0.0 -0.06 0.0000 1.0125 EP_2nd 0.0 0 1
## 11393 28 6 0 0.0 -0.06 0.0000 1.0125 PR_2nd -0.1 1 0
## 1249 28 7 0 0.0 -0.32 0.0000 0.9500 EP_2nd 0.0 0 1
## 11394 28 7 0 0.0 -0.32 0.0000 0.9500 PR_2nd -0.4 1 0
## 1250 28 8 0 0.0 -0.38 0.0000 0.3875 EP_2nd 0.0 0 1
## 11395 28 8 0 0.0 -0.38 0.0000 0.3875 PR_2nd 0.1 1 0
## 1251 28 9 0 0.3 -0.24 -0.0625 -0.0500 EP_2nd 0.5 0 1
## 11396 28 9 0 0.3 -0.24 -0.0625 -0.0500 PR_2nd 0.4 1 0
## 1252 28 10 0 0.4 0.18 0.5000 -0.0375 EP_2nd 0.0 0 1
## 11397 28 10 0 0.4 0.18 0.5000 -0.0375 PR_2nd 0.3 1 0Running the model
model.bi2age = nlme(value ~ - (PR*etaPR*(PR_0th-(baselinePR+(age_ind*baselinePR5))))
- (PR*etaPRage*(PR_0th-(baselinePR+(age_ind*baselinePR5)))*age_ind)
- (PR*zetaPR*PR_1st)
- (PR*zetaPRage*PR_1st*age_ind)
- (PR*betaPR*(EP_0th-(baselineEP+(age_ind*baselineEP5))))
- (PR*betaPRage*(EP_0th-(baselineEP+(age_ind*baselineEP5)))*age_ind)
- (PR*deltaPR*EP_1st)
- (PR*deltaPRage*EP_1st*age_ind)
- (PR*piPR*(PR_0th-(baselinePR+(age_ind*baselinePR5)))*(EP_0th-(baselineEP+(age_ind*baselineEP5))))
- (PR*piPRage*(PR_0th-(baselinePR+(age_ind*baselinePR5)))*(EP_0th-(baselineEP+(age_ind*baselineEP5)))*age_ind)
- (PR*xiPR*PR_1st*(EP_0th-(baselineEP+(age_ind*baselineEP5))))
- (PR*xiPRage*PR_1st*(EP_0th-(baselineEP+(age_ind*baselineEP5)))*age_ind )
- (EP*etaEP*(EP_0th-(baselineEP+(age_ind*baselineEP5))))
- (EP*etaEPage*(EP_0th-(baselineEP+(age_ind*baselineEP5)))*age_ind)
- (EP*zetaEP*EP_1st)
- (EP*zetaEPage*EP_1st*age_ind)
- (EP*betaEP*(PR_0th-(baselinePR+(age_ind*baselinePR5))))
- (EP*betaEPage*(PR_0th-(baselinePR+(age_ind*baselinePR5)))*age_ind)
- (EP*deltaEP*PR_1st)
- (EP*deltaEPage*PR_1st*age_ind )
- (EP*piEP*(EP_0th-(baselineEP+(age_ind*baselineEP5)))*(PR_0th-(baselinePR+(age_ind*baselinePR5))))
- (EP*piEPage*(EP_0th-(baselineEP+(age_ind*baselineEP5)))*(PR_0th-(baselinePR+(age_ind*baselinePR5)))*age_ind)
- (EP*xiEP*EP_1st*(PR_0th-(baselinePR+(age_ind*baselinePR5))))
- (EP*xiEPage*EP_1st*(PR_0th-(baselinePR+(age_ind*baselinePR5)))*age_ind),
data=dflong,
fixed=baselinePR + baselinePR5 +
etaPR + etaPRage + zetaPR + zetaPRage +
betaPR + betaPRage + deltaPR + deltaPRage +
piPR + piPRage + xiPR + xiPRage +
baselineEP + baselineEP5 +
etaEP + etaEPage + zetaEP + zetaEPage +
betaEP + betaEPage + deltaEP + deltaEPage +
piEP + piEPage + xiEP + xiEPage ~ 1,
random = baselinePR + baselineEP + baselinePR5 + baselineEP5
~ 1|id,
weights=varIdent(form=~1| PR), corr=corCompSymm(form=~1|id/time),
start = c(etaPR=.25,zetaPR=0.10, baselinePR = .70, baselinePR5 = -.30,
etaPRage=.05, zetaPRage=.1,
betaPR=.01, betaPRage=.1, deltaPR=-.05, deltaPRage=.1,
piPR=-.10, piPRage=-.10, xiPR=.01, xiPRage=.01,
etaEP=.35,zetaEP=0.01, baselineEP = 1.65, baselineEP5 = -.20,
etaEPage=.1, zetaEPage=.01,
betaEP=.10, betaEPage=.1, deltaEP=.05, deltaEPage=.1,
piEP=.01, piEPage=.01, xiEP=.01, xiEPage=.01),
control = list(maxiter=300), na.action = na.omit)
summary(model.bi2age)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: value ~ -(PR * etaPR * (PR_0th - (baselinePR + (age_ind * baselinePR5)))) - (PR * etaPRage * (PR_0th - (baselinePR + (age_ind * baselinePR5))) * age_ind) - (PR * zetaPR * PR_1st) - (PR * zetaPRage * PR_1st * age_ind) - (PR * betaPR * (EP_0th - (baselineEP + (age_ind * baselineEP5)))) - (PR * betaPRage * (EP_0th - (baselineEP + (age_ind * baselineEP5))) * age_ind) - (PR * deltaPR * EP_1st) - (PR * deltaPRage * EP_1st * age_ind) - (PR * piPR * (PR_0th - (baselinePR + (age_ind * baselinePR5))) * (EP_0th - (baselineEP + (age_ind * baselineEP5)))) - (PR * piPRage * (PR_0th - (baselinePR + (age_ind * baselinePR5))) * (EP_0th - (baselineEP + (age_ind * baselineEP5))) * age_ind) - (PR * xiPR * PR_1st * (EP_0th - (baselineEP + (age_ind * baselineEP5)))) - (PR * xiPRage * PR_1st * (EP_0th - (baselineEP + (age_ind * baselineEP5))) * age_ind) - (EP * etaEP * (EP_0th - (baselineEP + (age_ind * baselineEP5)))) - (EP * etaEPage * (EP_0th - (baselineEP + (age_ind * baselineEP5))) * age_ind) - (EP * zetaEP * EP_1st) - (EP * zetaEPage * EP_1st * age_ind) - (EP * betaEP * (PR_0th - (baselinePR + (age_ind * baselinePR5)))) - (EP * betaEPage * (PR_0th - (baselinePR + (age_ind * baselinePR5))) * age_ind) - (EP * deltaEP * PR_1st) - (EP * deltaEPage * PR_1st * age_ind) - (EP * piEP * (EP_0th - (baselineEP + (age_ind * baselineEP5))) * (PR_0th - (baselinePR + (age_ind * baselinePR5)))) - (EP * piEPage * (EP_0th - (baselineEP + (age_ind * baselineEP5))) * (PR_0th - (baselinePR + (age_ind * baselinePR5))) * age_ind) - (EP * xiEP * EP_1st * (PR_0th - (baselinePR + (age_ind * baselinePR5)))) - (EP * xiEPage * EP_1st * (PR_0th - (baselinePR + (age_ind * baselinePR5))) * age_ind)
## Data: dflong
## AIC BIC logLik
## 4012.106 4336.447 -1965.053
##
## Random effects:
## Formula: list(baselinePR ~ 1, baselineEP ~ 1, baselinePR5 ~ 1, baselineEP5 ~ 1)
## Level: id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## baselinePR 0.9942618 bslnPR bslnEP bslPR5
## baselineEP 0.5845750 -0.453
## baselinePR5 1.0050746 -0.967 0.446
## baselineEP5 0.6591746 0.301 -0.917 -0.339
## Residual 0.3117288
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | id/time
## Parameter estimate(s):
## Rho
## -0.03696949
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | PR
## Parameter estimates:
## 0 1
## 1.0000000 0.6940374
## Fixed effects: baselinePR + baselinePR5 + etaPR + etaPRage + zetaPR + zetaPRage + betaPR + betaPRage + deltaPR + deltaPRage + piPR + piPRage + xiPR + xiPRage + baselineEP + baselineEP5 + etaEP + etaEPage + zetaEP + zetaEPage + betaEP + betaEPage + deltaEP + deltaEPage + piEP + piEPage + xiEP + xiEPage ~ 1
## Value Std.Error DF t-value p-value
## baselinePR 0.4461862 0.09510647 20007 4.69144 0.0000
## baselinePR5 -0.3167300 0.09835930 20007 -3.22013 0.0013
## etaPR 0.2616361 0.00802823 20007 32.58952 0.0000
## etaPRage 0.0474020 0.01405966 20007 3.37149 0.0007
## zetaPR 0.1256797 0.01597025 20007 7.86962 0.0000
## zetaPRage 0.0649534 0.03135192 20007 2.07175 0.0383
## betaPR 0.0070875 0.00472762 20007 1.49917 0.1338
## betaPRage 0.0143810 0.00997241 20007 1.44208 0.1493
## deltaPR -0.0011223 0.01009696 20007 -0.11115 0.9115
## deltaPRage -0.0026946 0.01984130 20007 -0.13581 0.8920
## piPR -0.1959870 0.00772634 20007 -25.36609 0.0000
## piPRage -0.1129750 0.01894540 20007 -5.96319 0.0000
## xiPR 0.1224790 0.02225455 20007 5.50355 0.0000
## xiPRage -0.0963441 0.05805772 20007 -1.65945 0.0970
## baselineEP 1.1409302 0.05810806 20007 19.63463 0.0000
## baselineEP5 -0.0975591 0.06799328 20007 -1.43483 0.1513
## etaEP 0.2500435 0.00561669 20007 44.51794 0.0000
## etaEPage 0.2101846 0.01275648 20007 16.47669 0.0000
## zetaEP -0.0121149 0.01510152 20007 -0.80223 0.4224
## zetaEPage 0.0069440 0.02886452 20007 0.24057 0.8099
## betaEP 0.0475614 0.00762409 20007 6.23831 0.0000
## betaEPage 0.0601552 0.01504901 20007 3.99729 0.0001
## deltaEP 0.0113538 0.01894593 20007 0.59928 0.5490
## deltaEPage -0.0332415 0.03938798 20007 -0.84395 0.3987
## piEP -0.0350854 0.01023890 20007 -3.42667 0.0006
## piEPage 0.1105234 0.02507409 20007 4.40787 0.0000
## xiEP -0.0634076 0.03135948 20007 -2.02196 0.0432
## xiEPage 0.1140571 0.07176547 20007 1.58930 0.1120
## Correlation:
## bslnPR bslPR5 etaPR etaPRg zetaPR zetPRg betaPR betPRg deltPR dltPRg piPR piPRag xiPR xiPRag bslnEP bslEP5
## baselinePR5 -0.946
## etaPR -0.004 0.005
## etaPRage 0.004 0.003 -0.574
## zetaPR 0.000 -0.001 -0.175 0.100
## zetaPRage 0.000 -0.002 0.089 -0.081 -0.510
## betaPR -0.058 0.057 0.201 -0.115 -0.016 0.008
## betaPRage 0.029 -0.068 -0.097 0.170 0.009 -0.013 -0.479
## deltaPR 0.004 -0.004 -0.003 0.001 0.067 -0.034 0.005 -0.002
## deltaPRage -0.002 0.002 0.001 -0.002 -0.034 0.044 -0.003 0.000 -0.509
## piPR -0.029 0.028 0.151 -0.088 -0.047 0.024 0.286 -0.136 0.013 -0.007
## piPRage 0.012 -0.031 -0.063 0.118 0.019 -0.037 -0.117 0.202 -0.005 -0.003 -0.409
## xiPR 0.002 -0.002 -0.051 0.029 0.510 -0.259 -0.021 0.010 -0.010 0.005 -0.037 0.015
## xiPRage -0.001 0.002 0.019 -0.039 -0.196 0.462 0.007 0.007 0.004 -0.006 0.015 -0.002 -0.383
## baselineEP -0.437 0.420 -0.132 0.077 0.043 -0.022 0.026 -0.010 -0.001 0.001 0.005 -0.002 0.001 0.000
## baselineEP5 0.282 -0.333 0.114 -0.113 -0.037 0.022 -0.021 0.028 0.001 -0.001 -0.005 0.007 -0.001 0.000 -0.882
## etaEP -0.008 0.008 0.023 -0.011 -0.008 0.004 0.012 -0.008 -0.004 0.003 -0.019 0.006 0.011 -0.005 -0.004 0.004
## etaEPage 0.003 0.005 -0.008 -0.039 0.004 -0.001 -0.006 -0.052 0.003 0.000 0.008 -0.058 -0.005 0.003 0.001 -0.003
## zetaEP -0.006 0.006 -0.005 0.002 -0.002 0.001 0.029 -0.014 -0.039 0.020 0.012 -0.005 -0.005 0.001 0.012 -0.010
## zetaEPage 0.003 -0.001 0.002 -0.001 0.001 -0.002 -0.015 0.000 0.020 -0.039 -0.006 -0.002 0.003 0.000 -0.006 0.004
## betaEP 0.003 -0.002 0.181 -0.105 -0.058 0.030 -0.010 0.005 0.000 0.000 0.037 -0.017 -0.006 0.001 -0.055 0.048
## betaEPage -0.002 -0.001 -0.092 0.017 0.028 -0.013 0.008 -0.018 0.000 0.000 -0.018 0.022 0.002 -0.004 0.028 -0.020
## deltaEP 0.000 0.000 -0.018 0.009 -0.021 0.009 -0.007 0.004 -0.002 0.001 -0.005 0.002 0.005 -0.003 0.003 -0.002
## deltaEPage 0.000 0.001 0.010 0.005 0.009 -0.031 0.004 -0.001 0.001 -0.002 0.003 -0.003 -0.001 0.002 -0.001 -0.003
## piEP 0.001 -0.001 0.248 -0.142 -0.068 0.035 -0.001 0.000 -0.001 0.001 0.032 -0.015 -0.008 0.004 -0.066 0.057
## piEPage -0.001 0.002 -0.102 0.183 0.028 -0.017 -0.001 -0.026 0.000 -0.001 -0.014 0.009 0.003 -0.003 0.027 -0.051
## xiEP 0.000 0.000 0.009 -0.006 -0.010 0.005 0.005 -0.003 -0.002 0.001 0.009 -0.005 -0.014 0.005 0.003 -0.002
## xiEPage 0.000 0.000 -0.003 0.004 0.004 -0.002 -0.001 0.002 0.000 -0.003 -0.003 0.006 0.006 -0.002 -0.001 0.002
## etaEP etaEPg zetaEP zetEPg betaEP betEPg deltEP dltEPg piEP piEPag xiEP
## baselinePR5
## etaPR
## etaPRage
## zetaPR
## zetaPRage
## betaPR
## betaPRage
## deltaPR
## deltaPRage
## piPR
## piPRage
## xiPR
## xiPRage
## baselineEP
## baselineEP5
## etaEP
## etaEPage -0.443
## zetaEP 0.042 -0.018
## zetaEPage -0.021 0.013 -0.523
## betaEP 0.203 -0.088 -0.011 0.006
## betaEPage -0.100 0.210 0.005 -0.002 -0.508
## deltaEP 0.005 -0.001 0.076 -0.040 0.031 -0.017
## deltaEPage -0.002 -0.021 -0.036 0.052 -0.014 0.002 -0.482
## piEP 0.130 -0.057 0.025 -0.013 0.109 -0.056 -0.024 0.013
## piEPage -0.055 0.082 -0.011 0.002 -0.046 0.013 0.011 -0.038 -0.410
## xiEP 0.028 -0.013 0.269 -0.141 -0.048 0.023 -0.014 0.007 0.030 -0.013
## xiEPage -0.011 0.002 -0.118 0.094 0.022 0.010 0.006 0.033 -0.012 -0.021 -0.436
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -12.42212701 -0.32766484 -0.02244346 0.25416591 10.71588597
##
## Number of Observations: 20146
## Number of Groups: 112Plotting the model
# To plot the model, we need to extract the estimated parameters and then create a function
#...specifying the model. Finally, we can use the function and the parameters to create the plot.
# Creating an object with the effects of the model
t<-data.frame(summary(model.bi2age)$tTable)
base1<-t$Value[rownames(t)=="baselinePR"] # Saving the coefficient of baselinePR into a variable
base1.5<-t$Value[rownames(t)=="baselinePR5"] # Saving the coefficient of baselinePR5 into a variable
base2 <- t$Value[rownames(t)=="baselineEP"] # Saving the coefficient of baselineEP into a variable
base2.5 <- t$Value[rownames(t)=="baselineEP5"] # Saving the coefficient of baselineEP5 into a variable
# Saving the other parameters from the model
parameters = "a=0.261636114803884,b=age*0.0474019516318653,c=0.125679694055614,d=age*0.0649533605533045,e=0.00708750464547767,f=age*0.014380981530435,g=-0.00112231261221321,h=age*-0.00269456673039096,i=-0.195986988258638,j=age*-0.112975036154197,k=0.122478952303695,l=age*-0.0963441349977665,m=0.250043479046963,n=age*0.210184590451423,o=-0.0121149055677911,p=age*0.0069439754499835,q=0.0475613784731889,r=age*0.0601552129036021,s=0.0113538257604348,t=age*-0.0332414580225045,u=-0.035085377014092,v=age*0.110523421403763,w=-0.0634075623159452,x=age*0.114057126173259"
# Creating a function specifying the model
predictTraj <- function(t, prevState, parms) {
EP <- prevState[1] # EP[t] or yhat
PR <- prevState[2] # PR[t] or xhat
y1 <- prevState[3] # dEP[t] or dxhat
y2 <- prevState[4] # dPR[t] or dyhat
with(as.list(parms), {
dEP <- y1
dPR <- y2
dy1 <- - m*(EP-(base2+(base2.5*age))) - n*(EP-(base2+(base2.5*age))) - o*y1 - p*y1 - q*(PR-(base1+(base1.5*age))) - r*(PR-(base1+(base1.5*age))) - s*y2 - t*y2 - u*(EP-(base2+(base2.5*age)))*(PR-(base1+(base1.5*age))) - v*(EP-(base2+(base2.5*age)))*(PR-(base1+(base1.5*age))) - w*y1*(PR-(base1+(base1.5*age))) - x*y1*(PR-(base1+(base1.5*age)))
dy2 <- - a*(PR-(base1+(base1.5*age))) - b*(PR-(base1+(base1.5*age))) - c*y2 - d*y2 - e*(EP-(base2+(base2.5*age))) - f*(EP-(base2+(base2.5*age))) - g*y1 - h*y1 - i*(PR-(base1+(base1.5*age)))*(EP-(base2+(base2.5*age))) - j*(PR-(base1+(base1.5*age)))*(EP-(base2+(base2.5*age))) - k*y2*(EP-(base2+(base2.5*age))) - l*y2*(EP-(base2+(base2.5*age)))
res<-c(dEP,dPR,dy1,dy2)
list(res)
}
)
}
#----------------------------------------------------------#
# Age 2
times <- seq(0, 55, by=.05) # Creating a time variable
age=0 # Specifying the age indicator
parms= eval( parse(text=paste("c(", parameters, ")" ) ) ) # Inputting the parameters
## Start values for steady state
xstart <- c(EP=0, PR=1,y1=-.01, y2=-.01) # Make sure to match the order of the name above
## This will create the data
out1 <- as.data.frame(lsoda(xstart, times, predictTraj, parms, maxsteps=10000))
## This will plot the predicted data
plot(ts(out1))
#----------------------------------------------------------#
# Age 5
times <- seq(0, 55, by=.05)
age=1
parms= eval( parse(text=paste("c(", parameters, ")" ) ) )
## Start values for steady state
xstart <- c(EP=0, PR=1,y1=-.01, y2=-.01) # Make sure to match the order of the name above
## This will create the data
out2 <- as.data.frame(lsoda(xstart, times, predictTraj, parms,maxsteps=10000))
## This will plot the predicted data
plot(ts(out2))
###
# Creating plots with shaded under 0
#...coding a variable valence so that we can use it for the different colors
out1$valencePR[out1$PR>=0] <- "posPR"
out1$valencePR[out1$PR<0] <- "negPR"
out1$valenceEP[out1$EP>=0] <- "posEP"
out1$valenceEP[out1$EP<0] <- "negEP"
out2$valencePR[out2$PR>=0] <- "posPR"
out2$valencePR[out2$PR<0] <- "negPR"
out2$valenceEP[out2$EP>=0] <- "posEP"
out2$valenceEP[out2$EP<0] <- "negEP"
# Saving my options in a list to reproduce them later
myopts <- list(geom_area(aes(fill=valencePR), alpha = .4, position='stack'),
geom_line(colour = 'red', size = 1.5),
geom_area(aes(x = time, y = EP, fill=valenceEP), alpha = .4, position = 'stack'),
geom_line(aes(x = time, y = EP),colour = 'blue', size = 2),
theme_bw(),
xlab('Time'),
ylab('Predicted PR and EP Levels'),
xlim(0, 55),
ylim(-1.5, 3),
theme(axis.text=element_text(size=15), axis.title=element_text(size=15), panel.grid.major=element_blank(),
panel.grid.minor=element_blank()),
scale_x_continuous(breaks = seq(0,50,by=10), expand = c(0,0)),
scale_fill_manual(values=c("posPR"="red","negPR"="black","posEP"="blue","negEP"="dark blue"), guide=FALSE),
geom_hline(yintercept=0))
# For age 2
p.age2 <- ggplot(out1, aes(x = time, y = PR)) + myopts +
geom_text(data = NULL, x = 4, y = 2.7, label = "Age 2", size=7) +
geom_hline(yintercept=base1, colour="red") + geom_hline(yintercept=base2, colour="blue")## Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.# For age 5
p.age5 <- ggplot(out2, aes(x = time, y = PR)) + myopts +
geom_text(data = NULL, x = 4, y = 2.7, label = "Age 5", size=7) +
geom_hline(yintercept=base1 + base1.5, colour="red") + geom_hline(yintercept=base2 + base2.5, colour="blue")## Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.# Finally, this plot generates Figure 3 of the manuscript.
multiplot(p.age2,p.age5) # Function from http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/