R Code to simulate and describe a reaction-time decision process

bdd.R

library(ggplot2)
library(plyr)

theme_set(theme_bw())
theme_update(panel.border=element_blank())

# A reaction time experiment works as follows: On each trial, an
# observer is asked to view a stimulus and categorize it into one of
# multiple values. For instance, it may be a noisy motion stimulus,
# and the observer is asked to classify the motion as "leftward" or
# "rightward."
#
# The diffusion model for reaction time experiment data supposes that
# observers maintain an internal "decision variable," which evolves
# from the start of the trial. The stimulus adds a drift to this
# internal variable; when it drifts past either a lower or upper bound
# the observer makes a response at that time. The dicision variable
# also accumulates noise, so the evolution of the decision variable
# follows a Wiener process.
#
# (In prosaic terms, we are supposing that observers form decisions by
# implementing a sequential probability ratio test or more generally a
# Kalman filter, then we are trying to perform inference _about_ the
# observer's inferential process.)

# This function simulates a single reaction time trial, producing a named vector
# of the reaction time and the binary response.
sim_single_trial <- function
(
 strength=0,           #strength of signal
 bias=0,               #constant added to drift (parameter)
 startingpoint=0,      #another bias parameter
 noise=1,              #how much variance accumulates per unit time (>0)
 bound=1,              #how high the bounds on the accumulator (>0)
 timestep=0.02,        #how small the timestep
 ndt = 0.2,            #residual decision/reaction time (>0)
 fuzz = 0.05,          #response time fuzz (to hide timesteps) (>0)
 nsteps=250)
{
  decision_variable <- cumsum(rnorm(n=nsteps,
                             mean=(bias + strength)*timestep,
                             sd=sqrt(noise) * sqrt(timestep))) + startingpoint
  escape_step <- which(abs(decision_variable) > bound)[1]
  c(rt=ndt + escape_step * timestep + rnorm(1, sd=fuzz),
    response = decision_variable[escape_step] > 0)
}

#here is a randomly simulated trial:
sim_single_trial(strength=0.3)

#"strength" is the strength of the signal in the stimulus; a typical
#experiment will use a mixture of strengths, positive and negative,
#Here are the strengths we will use:
strengths <- as.vector((2^(10:15) / 10000) %o% c(-1, 1))
strengths

#Simulate 1000 trials at each strength.
rt.data <- rdply(1000, cbind(strength=strengths, ldply(strengths, sim_single_trial)))
rt.data$.n <- NULL

#We don't have patience for simulating extremely long trials, so we
#have limited the number of time steps. This resulted in some trials
#where no decision was recorded (in which case `rt` and `response` are
#`NA`.) In practice, response times are limited, so we drop extremely
#long reaction times (and hopefully model the data as right truncated.)
rt.data <- subset(rt.data, !(is.na(rt) | rt > 4))

#Here are ten randomly selected rows of 'rt.data':
rt.data[sample(nrow(rt.data), 10),]

#Function to fold the data so that "correct" responses (where observer
#responds in same direction as stimulus) are considered positive:
fold <- function(data)
    mutate(data,
           response = xor(response, strength<0),
           strength = abs(strength))

#Function to compute bernoulli means and CIs for plotting
bernoulli_mean_ci <- function(x) {
  conf = prop.test(sum(as.logical(x)), length(x))$conf.int
  data.frame(y = mean(x), ymin = conf[1], ymax = conf[2])
}

# As stimulus strength increases, observers are more likely to respond
# correctly. Qualitatively their responses are fit well by a
# logit
psi.plot <- (
  ggplot(fold(rt.data))
  + aes(x=strength, y=as.numeric(response))
  + stat_summary(fun.data=bernoulli_mean_ci, geom="pointrange")
  + labs(x="Coherence", y="Proportion correct")
  + geom_hline(y=0.5) + geom_vline(x=0)
  + geom_smooth(method="glm", formula=y~x-1,
                family=binomial(link=logit)))
print(psi.plot)

# The model also predicts a distribution of response times.
rt.plot <- (
  ggplot(fold(rt.data)) + aes(x=rt)
  + stat_density(position="identity", fill="NA", geom="line")
  + labs(x="Response time", y="Density", title="Response time distribution")
  + theme(axis.text.y = element_blank())
  + geom_hline(y=0)
  )
print(rt.plot)

#There are subtleties to how response times vary.  One hallmark of
#typical reaction time data is that (_when aggregating trials with
#varying stimulus strengths_) errors tend to have longer response times
#than correct responses:
slow.error.plot <- (
  rt.plot
  + aes(color=factor(response))
  + labs(color="Response\ncorrect"))
print(slow.error.plot)

#Another is to look at response time stratified by coherence. With
#stronger stimuli, response times are shorter.
rt.strength.plot <- (
  ggplot(fold(rt.data))
  + aes(rt, color=factor(strength))
  + stat_density(geom="line", position="identity")
  + scale_color_discrete("Strength", h=c(270,0), direction=1)
  + labs(x="Response time", y="Density",
         title="Distribution of response times by stimulus strength")
 )
print(rt.strength.plot)

#But it turns out that "errors are slow" is false once you also
#condition on stimulus strength. Errors seem slow because observers
#encounter a mixture of stimulus strengths; at any particular stimulus
#strength RT is independent of error.)
(ggplot(fold(rt.data)) + aes(x=factor(strength), y=rt, color=response)
 + stat_summary(fun.data=mean_cl_boot, geom="pointrange",
                position=position_dodge(width=0.25))
 + labs(title="RTs conditional on correct answer and stimulus strength",
        x="Stimulus strength",
        y="Response time",
        color="Response\ncorrect"))