Week 5 Parameter Estimation2

We will run the model fitting again and improve the optimization method a bit. Briefly say, we will repeat the parameter estimation with a number of different starting values and check if the estimates are always the same.
This optimization procedure

# Week 5: Parameter estimation
# First and Last name
# date
  • To start the new project, clear working environment with rm(list=ls()).
# clear working environment
rm(list=ls())
  • Check if you have installed packages
# Open "rtdists" package
library(rtdists)

5.1 Maximum Likelihood (ML)

Toss a coin 20 times (N=20) and observe 10 Heads. What is the likelihood of observing the data? To compute it, we can compute likelihood for each hypothetical θ. Here we compute likelihood for θ = .5 and θ = .75.

Likelihood.theta50<- .5^20 
Likelihood.theta75<- .75^10 *.25^10
Likelihood.Ratio<-Likelihood.theta50/Likelihood.theta75
Likelihood.Ratio
## [1] 17.75773

We can compute the Likelihood by our created likelihood function (see week2 for details).

# Create a LF function given binomial model
# inputs: theta, x and N
# output: Likelihood
LF <- function(LF.theta,x,N){ dbinom(x,N,LF.theta) }
LF50 <- LF(0.5, 10, 20)
LF75 <- LF(0.75, 10, 20)
LF50
## [1] 0.1761971
LF75
## [1] 0.009922275
LFratio <-LF50/LF75
LFratio
## [1] 17.75773

From the examples above, you might notice that the likelihood values are meaningful when you make comparison between them (e.g., likelihood ratio).

5.2 ML in DDM

# 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
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)
# Create a deviance function to quantify the difference between predicted and actual data given model and parameter values
# That is, we compute deviance as log likelihood: how likely the set of parameter values is true given data and model?
# Note1: The log likelihood is computed separately for correct and incorrect trials.
# Note2: The ddiffusion model is the model, which requires rt, response and DDM-related parameters
diffusion_deviance2 <- function(pars, dat)
{
  log_lik <- sum(log(ddiffusion(rt=dat[dat["accuracy"] == 1, "rt"],
                                response="upper",  v=pars[1], a=pars[2], t0=pars[3], z=pars[2]/2, s= 1)))
  log_lik <- log_lik +
    sum(log(ddiffusion(rt=dat[dat["accuracy"] == 0, "rt"],
                       response="lower",  v=pars[1], a=pars[2], t0=pars[3], z=pars[2]/2, s= 1)))
  if (any((-2*log_lik == 0))) return(1e6)
  return(-2*log_lik)
}

Last week, we estimated parameters for each subject as followed:

# Initialization
n<-0

# Pre-allocation (i.e., create an empty data frame and name each column)
results<-setNames(data.frame(matrix(ncol = 7, nrow = 0)), c("name", "Mean RT (ms)","Mean Accuracy","v", "a","t0 (ms)", "nll"))

# Main for loop to estimate parameters for each subject
for (pp in unique(rr98$id)) {
  print(pp) # display the current subject in the Console window, so that you know the current progress
  
  data_pp = rr98[rr98$id == pp,]
  
  fit <- nlminb(start=starting_values, 
                objective=diffusion_deviance2,
                lower=lb,# the lower boundary for each parameter. 
                upper=ub, # the upper boundary for each parameter.
                dat=data_pp)
  deviance <- fit[["objective"]]
  params <- fit[["par"]]
  
  # store result for the first to the last subject
  n <- n+1 #update the current number of subject
  results[n,] <- c(pp,                      #subject's name
                   round(agg.id[n,2]*1000), #Mean RT (ms)
                   agg.id[n,3],             #Mean Accuracy
                   round(params[1],2),      #drift rate
                   round(params[2],2),      #threshold 
                   round(params[3]*1000),   #non-decision time (ms)
                   deviance)                #negative log likelihood
}
results

This week, we will repeat the optimization procedure more than once and check if the results are consistent.

# 1. set up lower and upper boundary for each parameter
lb <-c(-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,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<-setNames(data.frame(matrix(ncol = 5, nrow = 0)), c("nRepetition","v", "a","t0 (ms)", "nll"))
run <- 1
optimization_runs <-10

# 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),
                       runif(n=1, min=0,  max=0.1))
  fit <- nlminb(start=starting_values, objective=diffusion_deviance2,
                lower = lb, 
                upper= ub, 
                dat=rr98) 
  deviance <- fit[["objective"]]
  params <- fit[["par"]]
  
  # save results (i.e., run#, a, v, t0, Discrepancy) for each iterate
  results[run,] <- c(run,
                     round(params[1],2),
                     round(params[2],2),
                     round(params[3]*1000), 
                     deviance)
  run <- run + 1} 

# display results in the Console window
results