Differential Equation Model
APA ATI Session V
Overview
The differential equation model was developed to examine dynamic process-related research questions.
Outline
A. Differentiate Data
B. Plot Derivatives
C. Fit Damped Oscillator Model
D. Impulse Plot
Preliminaries
Loading Libraries
Loading libraries used in this script.
# Check to see if necessary packages are installed, and install if not
packages <- c("rstudioapi", "psych", "ggplot2", "data.table", "plyr",
"reshape", "nlme", "lme4")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
install.packages(setdiff(packages, rownames(installed.packages())))
}
# Load packages
library(rstudioapi)
library(psych)
library(deSolve) #for the differential equation solvers
library(ggplot2) #for plotting
library(nlme) #for fitting mixed effets models
library(lme4) # for fitting mixed effects models
library(zoo) #for making time-series embedded data
library(reshape2) #for reshaping data
library(plyr) # for person-mean centeringLoading Data
Our example makes use of the day-level (EMA-type) AMIB data files wherein a subset of n = 30 individuals provided up to 21 days of daily assessments.
Loading day-level file (T = 8) and subsetting to variables of interest.
#set filepath for data file
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/AMIBshare_phase2_daily_2019_0501.csv"
#read in the .csv file using the url() function
AMIB_daily <- read.csv(file=url(filepath),header=TRUE)Add day of week variable. For IDs in the 100s, 200s, and 400s, day=0 was a Thursday and for IDs in the 300s, day=0 was a Monday. Recode day so that a Thursday will be a zero for the ids in the 300s. This will allow us to examine for weekly cycles
# recode day
AMIB_daily$daynew <- ifelse(AMIB_daily$id >=300 & AMIB_daily$id < 400, AMIB_daily$day-3, AMIB_daily$day)
#check new data
head(AMIB_daily[which(AMIB_daily$id == 203),c("id","day","daynew")])## id day daynew
## 1 203 0 0
## 2 203 1 1
## 3 203 2 2
## 4 203 3 3
## 5 203 4 4
## 6 203 5 5head(AMIB_daily[which(AMIB_daily$id == 301),c("id","day","daynew")])## id day daynew
## 228 301 0 -3
## 229 301 1 -2
## 230 301 2 -1
## 231 301 3 0
## 232 301 4 1
## 233 301 5 2In addition to the daynew variable, we will make use of affbalance, which was calculated as posaff-negaff.
AMIB_daily <- AMIB_daily[,c("id","daynew","affbalance")]We now person-center the affbalance scores to remove any intraindividual differences in equilibrium.
#calculate intraindividual means
AMIB_imeans <- ddply(AMIB_daily, "id", summarize,
affbalance_trait = mean(affbalance, na.rm=TRUE))
#merging person-level data into daily data
AMIB_daily <- merge(AMIB_daily,AMIB_imeans,by="id")
#calculating state variables
AMIB_daily$affbalance_state <- AMIB_daily$affbalance - AMIB_daily$affbalance_traitExamine the detrended data:
cols <- c("Raw"="#FC6170","Detrended"="#696CFF")
#pdf("Raw&DetrendedAffBalanceDataViz.pdf",width=9,height=5.75)
ggplot(data = AMIB_daily[which(AMIB_daily$id <= 226),], aes(x = daynew)) +
geom_line(aes(y=affbalance, colour="Raw"), size=1, linetype="solid") +
geom_line(aes(y=affbalance_state, colour="Detrended"), size=1, linetype="solid") +
scale_colour_manual(name="Data Type",
values=cols,
guide = guide_legend(override.aes = list(
linetype = c("solid","solid"),
size=1))) +
xlab("Day") +
ylab("Affect Balance") +
facet_wrap(~id) +
theme_classic() +
theme(axis.title=element_text(size=18),
axis.text=element_text(size=14, colour="black"),
plot.title=element_text(size=18, hjust=.5,face="bold"),
strip.text=element_text(size=12),
legend.text=element_text(size=14),
legend.title=element_text(size=14)) +
ggtitle("Plot of Affect Balance Time Series\nfor 9 Individuals")
#dev.off()Differentiate the Data
Differential equation models can be fit directly - after differentiating the data. This means that we will obtain estimates of derivatives for each observation.
There are several methods that can be used to estimate derivatives - for this lab section, we use the generalized local linear approximation (GLLA) approach. Further materials and code can be found on Steve Boker’s site (http://people.virginia.edu/~smb3u/GLLAfunctions.R)
For this lab section, we use slightly adapted (short) versions of Boker’s GLLA functions. The more general functions can be adapted for other settings (along with other functions for doing the time-delay embedding, e.g., embedd() in the tseriesChaos package)
GLLA
The GLLA function specific to this example is:
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)
#replacement for Embed for embed = 5
xzoo <- zoo(x)
xzooembed <- cbind(lag(xzoo,k=-2), lag(xzoo,k=-1), xzoo, lag(xzoo,k=1), lag(xzoo,k=2))
EMat <- as.matrix(xzooembed)
#Continuation
Estimates <- EMat%*%t(W)
return(Estimates)
}We then loop our “observed” data through the GLLA function to obtain estimates of the “observed” derivatives, that is, to obtain the “differentiated data”.
# This is the for-loop that will go through the data and calculate derivatives.
# it goes through all of the elements of 'curr_data', and uses GLLA
idlist <- unique(AMIB_daily$id)
#create open object
derivdf <- data.frame()
for(i in 1:length(idlist)){
#i <- 1 #for testing
curr_data = AMIB_daily[which(AMIB_daily$id == idlist[i]),]
derivs_x = GLLA(x = curr_data[,'affbalance_state'], embed = 5, tau = 1, deltaT = 1, order = 2)
colnames(derivs_x) = c('x_est', 'dx_est', 'd2x_est')
derivs_plus = cbind(curr_data, derivs_x)
derivdf<- rbind(derivdf,derivs_plus)
}Check out the new dataframe with the derivative estimates:
round(head(derivdf[,c("id","daynew","affbalance_state","x_est","dx_est","d2x_est")]),2)## id daynew affbalance_state x_est dx_est d2x_est
## 1 203 0 1.44 NA NA NA
## 2 203 1 1.24 NA NA NA
## 3 203 2 0.54 0.85 -0.68 -0.43
## 4 203 3 0.44 0.13 -1.33 -0.87
## 5 203 4 -1.56 -2.12 -0.64 1.09
## 6 203 5 -4.36 -2.80 0.09 1.70Summary statistics:
describe(derivdf[,c("affbalance_state","x_est","dx_est","d2x_est")])## vars n mean sd median trimmed mad min max range
## affbalance_state 1 627 0.00 1.44 0.15 0.07 1.29 -5.68 4.40 10.08
## x_est 2 503 -0.04 1.15 0.06 -0.02 1.07 -5.03 3.49 8.52
## dx_est 3 503 0.00 0.54 -0.03 -0.01 0.49 -1.66 1.93 3.59
## d2x_est 4 503 0.02 0.76 0.04 0.01 0.66 -2.60 2.83 5.43
## skew kurtosis se
## affbalance_state -0.53 0.61 0.06
## x_est -0.29 0.85 0.05
## dx_est 0.27 0.55 0.02
## d2x_est 0.07 0.84 0.03- Notice that the mean for
affbalance_state(= 0) andx_est(= -0.04) is approximately the same.
Derivative Plots
Lets plot the derivative estimates of \(x_{it}\), \(dx_{it}\), and \(d2x_{it}\) for a subset of individuals.
#pdf("DerivEstimatesDataViz.pdf",width=9, height=5.75)
cols <- c("x"="#FC6170", "dx"="#26BFBF", "d2x"="#FF8A47")
ggplot(data = derivdf[which(derivdf$id <= 205),], aes(x = daynew)) +
geom_line(aes(y=x_est, colour="x"), linetype = "solid", size=1) +
geom_line(aes(y=dx_est, colour="dx"), linetype = "solid", size=1) +
geom_line(aes(y=d2x_est, colour="d2x"), linetype = "solid", size=1) +
scale_colour_manual(name="Derivative",
values=cols,
guide = guide_legend(override.aes = list(
linetype = c("solid","solid", "solid"),
size=1))) +
xlab("Day") +
ylab("GLLA Derivative Estimates") +
facet_wrap(~id,nrow=5) +
theme_classic() +
theme(axis.title=element_text(size=16),
axis.text=element_text(size=10, colour="black"),
plot.title=element_text(size=16, hjust=.5,face="bold"),
legend.title=element_text(size=14),
legend.text=element_text(size=12),
strip.text=element_text(size=12)) +
ggtitle("Plot of x, dx, and d2x")
#dev.off()Fit Damped Linear Oscillator (DLO)
dlo_modfit <- nlme(d2x_est ~ eta_x*x_est + zeta_x*dx_est,
fixed = eta_x + zeta_x ~ 1,
random = eta_x ~ 1|id,
start = c(eta_x = -0.80, # frequency of x
zeta_x = -0.01), # damping of x
data=derivdf,
na.action = na.exclude,
verbose=TRUE)##
## **Iteration 1
## LME step: Loglik: -314.5108, nlminb iterations: 1
## reStruct parameters:
## id
## 1.130059
## PNLS step: RSS = 96.47977
## fixed effects: -0.5438015 -0.0253445
## iterations: 3
## Convergence crit. (must all become <= tolerance = 1e-05):
## fixed reStruct
## 6.054371e-01 2.168754e-06
##
## **Iteration 2
## LME step: Loglik: -314.5108, nlminb iterations: 1
## reStruct parameters:
## id
## 1.130062
## PNLS step: RSS = 96.47977
## fixed effects: -0.5438015 -0.0253445
## iterations: 1
## Convergence crit. (must all become <= tolerance = 1e-05):
## fixed reStruct
## 0.000000e+00 1.734997e-13summary(dlo_modfit)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: d2x_est ~ eta_x * x_est + zeta_x * dx_est
## Data: derivdf
## AIC BIC logLik
## 637.0216 653.9039 -314.5108
##
## Random effects:
## Formula: eta_x ~ 1 | id
## eta_x Residual
## StdDev: 0.1414668 0.4379597
##
## Fixed effects: eta_x + zeta_x ~ 1
## Value Std.Error DF t-value p-value
## eta_x -0.5438015 0.03293993 472 -16.508886 0.000
## zeta_x -0.0253445 0.03643417 472 -0.695625 0.487
## Correlation:
## eta_x
## zeta_x 0.004
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.033033909 -0.608529547 0.003441901 0.602124636 3.200024742
##
## Number of Observations: 503
## Number of Groups: 30intervals(dlo_modfit)## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## eta_x -0.60839970 -0.5438015 -0.47920324
## zeta_x -0.09679526 -0.0253445 0.04610626
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: id
## lower est. upper
## sd(eta_x) 0.09463504 0.1414668 0.2114742
##
## Within-group standard error:
## lower est. upper
## 0.4109590 0.4379597 0.4667344- Note, could only estimate the model with a random effect for \(\eta\) (not for both \(\eta\) and \(\zeta\))
Interpret Results from the Fitted DLO Model
- The \(\eta\) parameter (= -0.54, p = 0) indicates that the frequency by which affect balance fluctuates around its equilibrium is significantly different from zero.
- The \(\zeta\) parameter (= -0.03, p = 0.49), is non-significant, indicating that there is not damping in the system (i.e., individuals’ affect balance does not tend to move toward a stable equilibrium)
- There were also noteworthy individual differences in the frequency by which affect balance fluctuated around its equilibrium (\(\sigma^{2}\) = 0.02, 95%CI= [0.01, 0.04]).
Plot trajectories for select individuals
Save model implied scores for d2x_est
# remove missing data in d2x to then align with model predicted values
derivdf_sub <- na.omit(derivdf, cols = c("d2x_est"))
# get model predicted scores for d2x
derivdf_sub$pred_dlo <- predict(dlo_modfit)Get \(\eta\) value for each person:
dlo_re <- data.frame(dlo_modfit$coefficients$random)pdf("LowerEtaPeople.pdf",width=9, height=5.75)
ggplot(data = derivdf_sub[which(derivdf_sub$id %in% c(321, 316, 244)),], aes(x = daynew)) +
geom_line(aes(y=pred_dlo, colour="pred_dlo"), linetype = "solid", size=1) +
scale_x_continuous(breaks=c(0,7,14,21), limits=c(0,21)) +
xlab("Day") +
ylab("Model Implied d2x_est") +
facet_wrap(~id,nrow=3,
labeller = labeller(id = c("244" = "eta = -0.69",
"316" = "eta = -0.71",
"321" = "eta = -0.75"))) +
theme_classic() +
theme(axis.title=element_text(size=16),
axis.text.x=element_text(size=14, colour="black"),
axis.text.y=element_text(size=10, colour="black"),
plot.title=element_text(size=16, hjust=.5,face="bold"),
legend.position="none",
strip.text=element_text(size=14)) +
ggtitle(paste0("Lower Eta People"))## Warning: Removed 1 rows containing missing values (geom_path).dev.off()## quartz_off_screen
## 2pdf("HigherEtaPeople.pdf",width=9, height=5.75)
ggplot(data = derivdf_sub[which(derivdf_sub$id %in% c(439, 226, 403)),], aes(x = daynew)) +
geom_line(aes(y=pred_dlo, colour="pred_dlo"), linetype = "solid", size=1) +
scale_x_continuous(breaks=c(0,7,14,21), limits=c(0,21)) +
xlab("Day") +
ylab("Model Implied d2x_est") +
facet_wrap(~id,nrow=3,
labeller = labeller(id = c("226" = "eta = -0.36",
"403" = "eta = -0.36",
"439" = "eta = -0.36"))) +
theme_classic() +
theme(axis.title=element_text(size=16),
axis.text.y=element_text(size=10, colour="black"),
axis.text.x=element_text(size=14, colour="black"),
plot.title=element_text(size=16, hjust=.5,face="bold"),
legend.position="none",
strip.text=element_text(size=14)) +
ggtitle(paste0("Higher Eta People"))
dev.off()## quartz_off_screen
## 2Impulse Plot
Get equations from the nlme model:
dlo_modfit$call$model## d2x_est ~ eta_x * x_est + zeta_x * dx_estWrite out the model as a function:
dlo_mod <- function(times, y , parms) {
with(as.list(c(y, parms)), {
derivs = c(
dx = dx,
d2x = eta_x*x + zeta_x*dx
)
return(list(derivs))
})
}Extract parameter estimates from the model object:
parms_est <- summary(dlo_modfit)$tTable[,"Value"]Set initial values:
initvalues_dlo1 <- c(x = -5, dx = 0)
initvalues_dlo2 <- c(x = -0, dx = 0)
initvalues_dlo3 <- c(x = 5, dx = 0)Set times to integrate over:
times <- seq(from=-3, to=21, by = 1)Solve to get trajectory:
out1 <- ode(y = initvalues_dlo1, times = times, func = dlo_mod, parms = parms_est)
out2 <- ode(y = initvalues_dlo2, times = times, func = dlo_mod, parms = parms_est)
out3 <- ode(y = initvalues_dlo3, times = times, func = dlo_mod, parms = parms_est)
# convert to data.frame
outdata1 <- as.data.frame(out1)
outdata2 <- as.data.frame(out2)
outdata3 <- as.data.frame(out3)ggplot(data = outdata1, aes(x = time)) +
geom_line(aes(y=x), color="#122140", size=2) +
xlab("Day") +
ylab("Affect Balance") +
theme_classic() +
theme(axis.title=element_text(size=16),
axis.text=element_text(size=14, colour="black"),
plot.title=element_text(size=16, hjust=.5,face="bold"),
legend.text=element_text(size=14)) +
ggtitle("Damped Linear Oscillator Model\nInitial Value = -5")
ggplot(data = outdata2, aes(x = time)) +
geom_line(aes(y=x), color="#254280", size=2) +
xlab("Day") +
ylab("Affect Balance") +
theme_classic() +
theme(axis.title=element_text(size=16),
axis.text=element_text(size=14, colour="black"),
plot.title=element_text(size=16, hjust=.5,face="bold"),
legend.text=element_text(size=14)) +
ggtitle("Damped Linear Oscillator Model\nInitial Value = 0")
ggplot(data = outdata3, aes(x = time)) +
geom_line(aes(y=x), color="#4376E6", size=2) +
xlab("Day") +
ylab("Affect Balance") +
theme_classic() +
theme(axis.title=element_text(size=16),
axis.text=element_text(size=14, colour="black"),
plot.title=element_text(size=16, hjust=.5,face="bold"),
legend.text=element_text(size=14)) +
ggtitle("Damped Linear Oscillator Model\nInitial Value = 5")