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")