Week 7 Model evaluation
We extended the model based on experimental manipulation last week. To find out the best that balance goodness-of-fit and model complexity, we will conduct
1. model comparison, and
2. Ex post simulation
Agian, we need three elements: Data, Models, Deviation functions. You might have noticed that these materials are needed since week 4. That is because the model comparison relies on negative log-likelihood and/or corrected negative log-likelihood (e.g., AIC or BIC), which is computed by the combination of fitting the model to data and deviation function.
Now, suppose we include two models (DDM with varied drifts rate and DDM with varied thresholds) in our model space. To simplify the code, we only fit the data at group-level.
# Week 7: Model evaluation
# First and Last name
# date- To start the new project, clear working environment with
rm(list=ls()).
# clear working environment
rm(list=ls())- Load the package
# Open "rtdists" package
library(rtdists)- Load data
# Load data from the "rtdists" package
data(rr98)
# Double-click rr98 in your Environment panel (right-top window)
# In the dataset, there are 3 participants, who completed ~ 4000 trials in each condition (speed/accuracy).- Preprocess data to make sure the data and data structure are suitable for follow-up model fitting.
# Preprocess data
rr98 <- rr98[!rr98$outlier,] # remove outliers, as in original paper
rr98[,'RT'] <- rr98$rt
rr98[,'accuracy'] = 0 # add a new column labeled as accuracy, where incorrect responses are coded as 0
rr98[rr98$correct, 'accuracy'] = 1 # and correct responses are coded as 1
rr98$choice = rr98$response
rr98$response = as.character(rr98$response)
rr98[rr98$correct == 1, "response"] = "upper"
rr98[rr98$correct == 0, "response"] = "lower"
head(rr98)7.1 Set up models
7.1.1 DDM1
# Extend the model by splitting drift rate (v) based on instruction manipulations.
# parameters: v_speed, v_accuracy, a, t0
# Important: order of parameter matters!! always keep tracking the order of parameter by writing them down as a note.
DDM1 <- function(pars, dat)
{
log_lik <- sum(log(ddiffusion(rt=dat[dat["accuracy"] == 1 & dat["instruction"] == "speed", "rt"],
response="upper", v=pars[1], a=pars[3], t0=pars[4])))
log_lik <- log_lik +
sum(log(ddiffusion(rt=dat[dat["accuracy"] == 1 & dat["instruction"] == "accuracy", "rt"],
response="upper", v=pars[2], a=pars[3], t0=pars[4])))
log_lik <- log_lik +
sum(log(ddiffusion(rt=dat[dat["accuracy"] == 0 & dat["instruction"] == "speed", "rt"],
response="lower", v=pars[1], a=pars[3], t0=pars[4])))
log_lik <- log_lik +
sum(log(ddiffusion(rt=dat[dat["accuracy"] == 0 & dat["instruction"] == "accuracy", "rt"],
response="lower", v=pars[2], a=pars[3], t0=pars[4])))
if (any((-2*log_lik == 0))) return(1e6)
return(-2*log_lik)
}7.1.2 DDM2
# Extend the model by splitting threshold based on instruction manipulations
# parameters: v, a_speed, a_accuracy, t0
# Note: order of parameter matters!! always keep tracking the order of parameter.
DDM2 <- function(pars, dat)
{
log_lik <- sum(log(ddiffusion(rt=dat[dat["accuracy"] == 1 & dat["instruction"] == "speed", "rt"],
response="upper", v=pars[1], a=pars[2], t0=pars[4])))
log_lik <- log_lik +
sum(log(ddiffusion(rt=dat[dat["accuracy"] == 1 & dat["instruction"] == "accuracy", "rt"],
response="upper", v=pars[1], a=pars[3], t0=pars[4])))
log_lik <- log_lik +
sum(log(ddiffusion(rt=dat[dat["accuracy"] == 0 & dat["instruction"] == "speed", "rt"],
response="lower", v=pars[1], a=pars[2], t0=pars[4])))
log_lik <- log_lik +
sum(log(ddiffusion(rt=dat[dat["accuracy"] == 0 & dat["instruction"] == "accuracy", "rt"],
response="lower", v=pars[1], a=pars[3], t0=pars[4])))
if (any((-2*log_lik == 0))) return(1e6)
return(-2*log_lik)
}7.2 Model fitting
7.2.1 Fitting DDM1
# 1. set up lower and upper boundary for each parameter
# Also, check if the same range is used for the same parameter
lb <-c(-5,-5,0.5,0.1) # Lower Boundary for each parameter. Order is important. Here, drift rate, threshold and non-decision time
ub <-c(5,5,5,0.5) # Upper Boundary for each parameter. Order is important. Here, drift rate, threshold and non-decision time
# 2. Create optimization function with an additional input: number of repetition.
# (The code is adapted from Mikhail Spektor's).
results.DDM1<-setNames(data.frame(matrix(ncol = 6, nrow = 0)), c("nRepetition","vspeed","vaccuracy", "a","t0 (ms)", "nll"))
run <- 1
optimization_runs <-5
# Run the parameter estimation for several times (this will take few seconds/minutes to complete)
while (run <= optimization_runs){
starting_values <- c(runif(n=1, min=-5, max=5), # selected from the same range for the same parameter
runif(n=1, min=-5, max=5), # selected from the same range for the same parameter
runif(n=1, min=0.5, max=5),
runif(n=1, min=0, max=0.1))
# parameter estimation procedure using nlminb
fitDDM1 <- nlminb(start=starting_values, objective=DDM1,
lower = lb,
upper= ub,
dat=rr98)
deviance <- fitDDM1[["objective"]]
params <- fitDDM1[["par"]]
# save results for each iterate
results.DDM1[run,] <- c(run,
round(params[1],4),
round(params[2],4),
round(params[3],4),
round(params[4],4),
deviance)
run <- run + 1}
# display results in the Console window
results.DDM17.2.2 Fitting DDM2
# 1. set up lower and upper boundary for each parameter
# Also, check if the same range is used for the same parameter
lb <-c(-5,0.5,0.5,0.1) # Lower Boundary for each parameter. Order is important. Here, drift rate, threshold and non-decision time
ub <-c(5,5,5,0.5) # Upper Boundary for each parameter. Order is important. Here, drift rate, threshold and non-decision time
# 2. Create optimization function with an additional input: number of repetition.
# (The code is adapted from Mikhail Spektor's).
results.DDM2<-setNames(data.frame(matrix(ncol = 6, nrow = 0)), c("nRepetition","v","a_speed","a_accuracy", "t0 (ms)", "nll"))
run <- 1
optimization_runs <-5
# Run the parameter estimation for several times (this will take few seconds/minutes to complete)
while (run <= optimization_runs){
starting_values <- c(runif(n=1, min=-5, max=5),
runif(n=1, min=0.5, max=5), # selected from the same range for the same parameter
runif(n=1, min=0.5, max=5), # selected from the same range for the same parameter
runif(n=1, min=0, max=0.1))
# parameter estimation procedure using nlminb
fitDDM2 <- nlminb(start=starting_values, objective=DDM2,
lower = lb,
upper= ub,
dat=rr98)
deviance <- fitDDM2[["objective"]]
params <- fitDDM2[["par"]]
# save results for each iterate
results.DDM2[run,] <- c(run,
round(params[1],4),
round(params[2],4),
round(params[3],4),
round(params[4],4),
deviance)
run <- run + 1}
# display results in the Console window
results.DDM27.3 Information criteria
7.3.1 Compute AIC and BIC
# predetermined the number of parameters
nParam.DDM1 <-4;
nParam.DDM2 <-4;
# compute AIC and BIC based on nLL, number of parameter and number of data points.
AIC.DDM1 <- results.DDM1[1,"nll"]+2*nParam.DDM1 # the nll is already -2*negative log-likelihood
AIC.DDM2 <- results.DDM2[1,"nll"]+2*nParam.DDM2
BIC.DDM1 <- results.DDM1[1,"nll"]+nParam.DDM1*log(length(rr98[,1]))
BIC.DDM2 <- results.DDM2[1,"nll"]+nParam.DDM2*log(length(rr98[,1])) 7.3.2 Comparison
With AIC and BIC for each model, we can make comparison between models.
# smaller is better
if (AIC.DDM1 > AIC.DDM2){
print("DDM2 is better")
}else {
print("DDM1 is better")}
if (BIC.DDM1 > BIC.DDM2){
print("DDM2 is better")
}else {
print("DDM1 is better")}7.4 Ex Post Simulation
After the model comparison using AIC and/or BIC, we only recognize the relatively good model with relative balance of goodness-of-fit and complexity). To further verify if the model captures the important behavioral patterns, we use the estimated parameter values to run the simulation and ideally, reproduce the patterns of interest.
# Extract the data we need for simulating data using DDM1 (with varied drift rate)
# if drift rate is divided by instruction manipulation rather than strength and the order of manipulations, then the only information influencing model simulation is number of trials for each manipulation.
nTrial.speed <- sum(rr98['instruction'] == 'speed')
nTrial.accuracy <- sum(rr98['instruction'] == 'accuracy')
nTrial.speed
nTrial.accuracy7.4.1 DDM1: Simulate data
First, we simulate data with estimated parameters using rdiffusion, which requires number of data and parameters.
# Simulate data with existing function
DDM1simu.speed <- rdiffusion(nTrial.speed, v = results.DDM1[1,2], a = results.DDM1[1,4], t0 = results.DDM1[1,5])
DDM1simu.accuracy <- rdiffusion(nTrial.accuracy, v = results.DDM1[1,3], a = results.DDM1[1,4], t0 = results.DDM1[1,5])Second, summarize the accuracy rate.
#summarize accuracy rate for simulated data and real data
sum(DDM1simu.speed["response"] == "upper")/length(DDM1simu.speed[,1])
sum(DDM1simu.accuracy["response"] == "upper")/length(DDM1simu.speed[,1])
sum(rr98["instruction"] == "speed" & rr98["response"] == "upper")/sum(rr98["instruction"] == "speed")
sum(rr98["instruction"] == "accuracy" & rr98["response"] == "upper")/sum(rr98["instruction"] == "accuracy")Third, we can then visualize the simulated and real data via either histogram or ttest or any function you want to compare the real and simulated data. Agian, ideally the simulated data from the model should be able to reproduce important behavioral patterns/effects.
# Visualize the simulated and real data (adapt from Jannnick's code)
#Plot RTs distributions
par(mfrow=c(2,1)) #show multiple plots in 1 Graphic
###RTs in speed condition
## real data
ID.speed = rr98["instruction"] == "speed"
hist(rr98[ID.speed,"rt"]*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM1: varied drift. Speed condition (real data)"), col = c("gray"))
abline(v=mean(rr98[ID.speed,"rt"]*1000),lwd=5, lty=2) #Add abline to visualize the mean of RT
## simulated data
hist(DDM1simu.speed[,"rt"]*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM1: varied drift. Speed condition (simulated data)"), col = c("blue"))
#Add abline to visualize the mean of RT
abline(v=mean(DDM1simu.speed[,"rt"]*1000),lwd=5, lty=2)#Add abline to visualize the mean of RT
###RTs in accuracy condition
## real data
par(mfrow=c(2,1)) #show multiple plots in 1 Graphic
ID.accuracy = rr98["instruction"] == "accuracy"
hist(rr98[ID.accuracy,"rt"]*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM1: varied drift. Accuracy condition (real data)"), col = c("gray"))
abline(v=mean(rr98[ID.accuracy,"rt"]*1000),lwd=5, lty=2)#Add abline to visualize the mean of RT
## simulated data
hist(DDM1simu.accuracy[,"rt"]*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM1: varied drift. Accuracy condition (simulated data)"), col = c("green"))
abline(v=mean(DDM1simu.accuracy[,"rt"]*1000),lwd=5,lty=2)#Add abline to visualize the mean of RT
###RTs in both conditions
## real data
par(mfrow=c(2,1)) #show multiple plots in 1 Graphic
hist(rr98[,"rt"]*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM1: varied drift. Real data regardless of condition"), col = c("gray"))
abline(v=mean(rr98[,"rt"]*1000),lwd=5, lty=2)#Add abline to visualize the mean of RT
## simulated data
hist(c(DDM1simu.speed[,"rt"], DDM1simu.accuracy[,"rt"])*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM1: varied drift. Simulated data regardless of condition"), col = c("yellow"))
abline(v=mean(c(DDM1simu.speed[,"rt"], DDM1simu.accuracy[,"rt"])*1000),lwd=5,lty=2)#Add abline to visualize the mean of RT 7.4.2 DDM2: Simulate data
Here, we repeat the same procedure to test DDM2.
# Simulate data with existing function
DDM2simu.speed <- rdiffusion(nTrial.speed, v = results.DDM2[1,2], a = results.DDM2[1,3], t0 = results.DDM2[1,5])
DDM2simu.accuracy <- rdiffusion(nTrial.accuracy, v = results.DDM2[1,2], a = results.DDM2[1,4], t0 = results.DDM2[1,5])#summarize accuracy rate for simulated data and real data
sum(DDM2simu.speed["response"] == "upper")/length(DDM2simu.speed[,1])
sum(DDM2simu.accuracy["response"] == "upper")/length(DDM2simu.speed[,1])
sum(rr98["instruction"] == "speed" & rr98["response"] == "upper")/sum(rr98["instruction"] == "speed")
sum(rr98["instruction"] == "accuracy" & rr98["response"] == "upper")/sum(rr98["instruction"] == "accuracy")# Visualize the simulated and real data (adapt from Jannnick's code)
#Plot RTs distributions
par(mfrow=c(2,1)) #show multiple plots in 1 Graphic
###RTs in speed condition
## real data
ID.speed = rr98["instruction"] == "speed"
hist(rr98[ID.speed,"rt"]*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM2: varied threshold. Speed condition (real data)"), col = c("gray"))
#Add abline to visualize the mean of RT
abline(v=mean(rr98[ID.speed,"rt"]*1000),lwd=5, lty=2)
## simulated data
hist(DDM2simu.speed[,"rt"]*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM2: varied threshold. Speed condition (simulated data)"), col = c("blue"))
#Add abline to visualize the mean of RT
abline(v=mean(DDM2simu.speed[,"rt"]*1000),lwd=5, lty=2)
###RTs in accuracy condition
## real data
par(mfrow=c(2,1)) #show multiple plots in 1 Graphic
ID.accuracy = rr98["instruction"] == "accuracy"
hist(rr98[ID.accuracy,"rt"]*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM2: varied threshold. Accuracy condition (real data)"), col = c("gray"))
abline(v=mean(rr98[ID.accuracy,"rt"]*1000),lwd=5, lty=2)#Add abline to visualize the mean of RT
## simulated data
hist(DDM2simu.accuracy[,"rt"]*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM2: varied threshold. Accuracy condition (simulated data)"), col = c("green"))
abline(v=mean(DDM2simu.accuracy[,"rt"]*1000),lwd=5,lty=2)#Add abline to visualize the mean of RT
###RTs in both conditions
## real data
par(mfrow=c(2,1)) #show multiple plots in 1 Graphic
hist(rr98[,"rt"]*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM2. Real data regardless of condition"), col = c("gray"))
abline(v=mean(rr98[,"rt"]*1000),lwd=5, lty=2)#Add abline to visualize the mean of RT
## simulated data
hist(c(DDM2simu.speed[,"rt"], DDM2simu.accuracy[,"rt"])*1000,xlab="RT (ms)", xlim=c(200,1000), breaks="FD", main=paste("DDM2. Simulated data regardless of condition"), col = c("yellow"))
abline(v=mean(c(DDM2simu.speed[,"rt"], DDM2simu.accuracy[,"rt"])*1000),lwd=5,lty=2)#Add abline to visualize the mean of RT According to the visual inspection, the difference between real RT distribution and simulated RT distribution are similar if we didn’t split the RT distribution based on the instruction manipulation. This is the case for both DDM1 and DDM2. When we split the data into speed condition and accuracy condition, the difference is more obvious if data is simulated by DDM1, which is not worse model in the model space.
7.5 Exercise
7.5.1 Task
Include original model (with only three parameters: v, a, t0) and model with varied t0 in your model space. Is the DDM2 (varied threshold) is still the winning model? Can the models reproduce the behavioral patterns (i.e., rt and accuracy rate)?