This example is designed to give you a better feel for the importance of evidence accumulaton to a the threshold in a DDM and its effects on decision-making. It is written using R Markdown. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

You can download the code file used to create this page at http://decisionneurolab.com/resources/Intro_to_DDM/exercises/threshold.Rmd.

Author: Cendri Hutcherson

Date last modified: May 1, 2019

I’m assuming that you’ve already read the introductory material describing the basic notion of sequential accumulation to bound. Here, we’ll play around with a function that computes a choice and an RT using these two principles.

Example 2.1. Accumulation of evidence to bound

Let’s start by writing a simple version of this function in R.

ddm = function(driftrate = .08, threshold = .15, noise = .1, precision = .001){
  # function ddm:
  # driftrate: average strength of the evidence
  # threshold: the choice-determining threshold
  # noise: standard deviation of the evidence, by convention set to .1
  # precision: the size of the time-step for accumulating evidence (default = .001s)
  
  # 1. start at time t = 1
  t = 1
  
  # 2. start with evidence at 0
  accumulatedEvidence = 0 
  
  # 3. use a while loop that will break when the evidence passes +/- the threshold
  while(abs(accumulatedEvidence[t]) < threshold){
    t = t+1
    
    # 4. take a noisy sample, scaled for the temporal precision with which we are estimating
    noisySample = driftrate*precision + rnorm(1, mean = 0, sd = noise*sqrt(precision))
    
    # 5. add it to the previous evidence
    accumulatedEvidence[t] = accumulatedEvidence[t - 1] + noisySample
  }
  
  # 6. return a data structure with the choice, the RT (in s), and the observed accumulated evidence
  choiceInfo = list(choice = sign(accumulatedEvidence[t]), rt = t*precision, drift = accumulatedEvidence)
  
  }

Above is a very simple function to compute a DDM. It takes as input 4 arguments: the desired drift rate, the threshold, the noise, and the temporal step size used to calculate the total evidence in a single sample.

Let’s take this baby out for a stroll and see what it does!

aChoice = ddm() # compute information about a single choice, with the default options

# what was the choice?
aChoice$choice
## [1] 1
# how long did it take?
aChoice$rt
## [1] 1.99
# what did the accumulation of evidence look like?
plot(aChoice$drift, type = 'l', col = 'blue', xlab = 'Time (in ms)', ylim = c(-.15, .15), ylab = "Accumulated Evidence")

# draw in the thresholds and the zero line
segments(c(0,0),c(.15, -.15), c(1000*aChoice$rt+50,1000*aChoice$rt+50), c(.15,-.15), lwd = 2) 
segments(0,0,1000*aChoice$rt+50,0, lty = 2)

Well, that looks like some accumulation of evidence right there, doesn’t it? You can see that the accumulated evidence should cross one of the thresholds right at the point where the rt indicated.

Example 2.2. DDMs are like fingerprints. No two are ever the same.

What would happen if we ran another instantiation of the DDM, with exactly the same parameters?

anotherChoice = ddm()

# what was the choice?
anotherChoice$choice
## [1] -1
# how long did it take?
anotherChoice$rt
## [1] 1.54
# what did the accumulation of evidence look like?
plot(anotherChoice$drift, type = 'l', col = 'blue', xlab = 'Time (in ms)', ylim = c(-.15, .15), 
     ylab = "Accumulated Evidence")

# draw in the thresholds and the zero line
segments(c(0,0),c(.15, -.15), c(1000*anotherChoice$rt+50,1000*anotherChoice$rt+50), c(.15,-.15), lwd = 2) 
segments(0,0,1000*anotherChoice$rt+50,0, lty = 2)

As you can see, due to the randomness of the noise, no two accumulation trajectories ever look precisely the same.

Example 2.3. Constructing choice distributions

Given that no two DDM realizations ever look exactly the same, one might ask “Are there some regularities that we can use to describe the resulting choices and RTs?” Well, let’s take a peek. To do this, we’ll create a function to simulate and plot the frequency of choices at different RTs. Then, we’ll simulate with the default parameters, and graph the resulting choices and RTs in a histogram.

simAndPlot = function(driftrate = .08, threshold = .15, nSim = 1000){
  # initialize vectors to hold the simulated choices and RTs
  Choices = RTs = vector(mode = 'numeric', length = nSim)
  
  # run the simulations and save the data
  for(sim in 1:nSim){
    temp = ddm(driftrate = driftrate, threshold = threshold)
    Choices[sim] = temp$choice
    RTs[sim] = temp$rt
  }
  
  # Plot the data!
  # 1. Get the frequencies of YES and NO responses in specific RT bins, determined by the maximum RT
  maxRT = ceiling(max(RTs))
  rtBins = seq(0,maxRT, length = 40)
  freqYes =   hist(RTs[Choices == 1], breaks = rtBins, plot = FALSE)$counts
  freqNo =   hist(RTs[Choices == -1], breaks = rtBins, plot = FALSE)$counts
  
  
  # 2. plot the frequencies
  par(mfrow = c(2,1), mar = c(0,4,5,1) + .2)
  yMax = ceiling(max(c(freqYes, freqNo)/50))*50
  hist(RTs[Choices == 1], breaks = rtBins, ylim = c(0,yMax), main = '', xaxt = 'n'
       , col = 'red')
  par(mar = c(4,4,1,1) + .2)
  hist(RTs[Choices == -1], breaks = rtBins, ylim = c(yMax,0), main = '', xlab = 'RT'
       , col = 'blue')
  
  # 3. Add descriptive text to the plot
  text(maxRT * .5, yMax * .4, paste('% YES: ', format(mean(Choices == 1), digits = 3)), adj = 0)
  text(maxRT * .5, yMax * .6, paste('Ave. YES RT: ', format(mean(RTs[Choices == 1]), digits = 3)), adj = 0)
  text(maxRT * .5, yMax * .8, paste('Ave. NO RT: ', format(mean(RTs[Choices == -1]), digits = 3)), adj = 0)
}

# Now let's use our function!
simAndPlot(nSim = 1000)

What you can see is that, given a drift rate of +.08 and threshold of .15, YES responses are much more frequent than NO responses (~90% vs. 10%), and the average RT tends to be about 1.5 seconds or so.

Example 2.4. Another distribution

Let’s examine what happens if we double the drift rate (to .16). First, think about what you predict will happen to the distribution of the choices and RTs. Ok, let’s test your intuition.

simAndPlot(driftrate = .16)

Was it what you predicted?

Example 2.5. Negative drift rates.

Let’s see what happens if the evidence is on average negative rather than positive.

simAndPlot(driftrate = -.08, threshold = .15)

Example 2.6. Changing the threshold.

Let’s say that instead of changing the drift rate, we changed the treshold. We’ll try a couple of examples here. Note how changing the threshold changes the likelihood of a yes response, and the associated RTs.

simAndPlot(driftrate = .08, threshold = .05)

simAndPlot(driftrate = .08, threshold = .1)

simAndPlot(driftrate = .08, threshold = .2)

Exercises/Test Yourself

  1. What is the acceptance rate and average RT when the evidence is .2 and the threshold is .2?
  2. What is the acceptance rate and average RT when the evidence is .1 and the threshold is .1? Compare this to the case in Exercise 1.
  3. Play around with different instantiations of the threshold and drift rate, using 5000 simulations instead of 1000. For each instantiation, compare the average RT for YES and NO. How do they tend to compare? Are the RTs for YES faster, slower, or the same? Does it depend on the drift rate or the threshold?