Distribution functions for accumulator models in R Henrik Singmann Matthew Gretton Scott Brown Andrew Heathcote Speed versus accuracy in Observed Behavior Lexical Decision Task Evidence Accumulation Models ID: 504402
Download Presentation The PPT/PDF document "rtdists" is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
Slide1
rtdists: Distribution functions for accumulator models in R
Henrik SingmannMatthew GrettonScott BrownAndrew HeathcoteSlide2
Speed versus accuracy in Observed Behavior
Lexical Decision TaskEvidence Accumulation Models
Latent Parameters:quality/rate of evidence accumulation
response caution (position of threshold)(non-decision time, position of start point, …)
Exp
. 1;
Wagenmakers
, Ratcliff, Gomez, &
McKoon
,
2008Slide3
RT distributions
Types of Errors
Ratcliff
/Drift Difussion Model
Fast Errors
Slow Errors
v:
for
identifiability
:
s
within
-trials
= 1 [
or
0.1]
t
0Slide4
The Math: ugly and slow(ish)
Density of simple diffusion (an infinite sum, converges slowly for large t)
Another approximation (better for large t, Navarro & Fuss, 2009)
Between-trial variability's require numerical integration (very slow!), the best code for this has been made available by the Voss brothers.Initial fast-dm, Voss & Voss (2007),
fast method for CDF, Voss & Voss (2008), solve PDE:initial conditions
Voss
, Voss &
Lerche
(2015) released f
ast
-dm-30
likelihoods, which does likelihoods.and boundary conditions
f
fSlide5
LBA
Linear Ballistic AccumulatorThe “Race Equation
”
Brown
& Heathcote (2008)
"
simplest
complete model of choice
response time"
The density function for choice
i
reaching threshold at time
t
.
The probability that choice
j
has not reached threshold by time
t
.
+
nondecision
time t
0
and
variability
s
t0
ASlide6
RTdists
R package providing distribution functions for Diffusion model and LBA:density functions (PDF): d… [
n1PDF for LBA]
cumulative distribution function (CDF): p… [n1CDF
for LBA]random derivatives (RNG): r…already available from CRAN:
install.packages
("
rtdists
")
Diffusion model uses C code by Voss brothers.
LBA model for four different drift rate distributions (Normal, Gamma,
Frechet
, log-normal)
Contains no fitting interface, but allows any type of fitting: maximum likelihood, Bayesian, Chi-Square
random derivatives allow
assessment
of
model
fit via
simulationSlide7
> require(
rtdists)
> s1
<- rlba_norm
(1000, A = 0.5, b = 1, t0 = 0.4, mean_v
= c(2, 1),
sd_v
= 1)
>
prop.table
(table(s1$response
))
# 1 2
# 0.714
0.286
> head(s1)
#
rt
response
#
1
1.1367528 2
#
2
0.8015599 2
#
3
0.8084981 2
#
4
0.7680607 2
# 5
0.6380967 1
#
6
0.6702449 1Slide8
> s1
<- rlba_norm
(1000, A = 0.5, b = 1, t0 = 0.4, mean_v
= c(2, 1), sd_v = 1)
> prop.table(table(s1$response))
# 1 2
# 0.714
0.286
ll_lba
<- function(pars,
rt
, response) {
ds <- split(
rt
, f = response)
d1 <- n1PDF(ds[[1]], A = pars["A"], b = pars["A"]+pars["b"], t0 = pars["t0"],
mean_v
= pars[c("v1", "v2")],
sd_v
= 1, silent=TRUE)
d2 <- n1PDF(ds[[2]], A = pars["A"], b = pars["A"]+pars["b"], t0 = pars["t0"],
mean_v
= pars[c("v2", "v1")],
sd_v
= 1, silent=TRUE)
d <- c(d1, d2)
if (any(d == 0)) return(1e6)
else return(-sum(log(d)))
}
start
<- c(
runif
(3, 0.3, 3), runif(2, 0, 0.2))
names
(
start
) <- c("A", "v1", "v2", "b", "t0
")
nlminb
(
start
,
ll_lba
,
lower
= 0,
rt
=s1$rt,
response
=s1$response
)
# $par
# A v1 v2 b t0
# 0.4336521 1.9278475 0.9856471 0.5177903 0.3962212
#
# $
objective
# [1] 137.3075
#
# $
convergence
# [1] 0
#
# $
iterations
# [1] 41
#
# $
evaluations
#
function
gradient
# 57 242
# # $message# [1] "relative convergence (4)"Slide9
> require(
rtdists)
> s1
<- rlba_norm
(1000, A = 0.5, b = 1, t0 = 0.4, mean_v
= c(2, 1),
sd_v
= 1)
> s2
<-
rlba_norm
(1000, A = 0.5, b = 0.6, t0 = 0.4,
mean_v
= c(2, 1),
sd_v
= c(0.5, 1
))
>
prop.table
(table(s2$response
))
# 1 2
# 0.723 0.277Slide10
> s2
<- rlba_norm
(1000, A = 0.5, b = 0.6, t0 = 0.4,
mean_v
= c(2, 1), sd_v = c(0.5, 1))
>
prop.table
(table(s2$response
))
# 1 2
# 0.723 0.277
ll_lba
<-
function
(pars,
rt
,
response
) {
ds
<-
split
(
rt
, f =
response
)
d1 <- n1PDF(
ds
[[1]], A = pars["A"], b = pars["A"]+pars["b"], t0 = pars["t0"],
mean_v = pars[c("v1", "v2")],
sd_v = c(pars["sv"], 1),
silent
=TRUE)
d2 <- n1PDF(
ds
[[2]], A = pars["A"], b = pars["A"]+pars["b"], t0 = pars["t0"],
mean_v
= pars[c("v2", "v1")],
sd_v
= c(1, pars["
sv
"]),
silent
=TRUE)
d <- c(d1, d2)
if
(
any
(d == 0))
return
(1e6)
else
return
(-
sum
(log(d)))
}
start
<- c(
runif
(3, 0.3, 3),
runif
(3, 0, 0.2))
names
(
start
) <- c("A", "v1", "v2", "b", "t0", "
sv
")
nlminb(start, ll_lba, lower = 0, rt=s2$rt, response=s2$response)# $par# A v1 v2 b t0 sv # 0.5832459 2.2496939 1.2850416 0.1522494 0.3838269 0.4766951 # # $
objective
# [1] -601.8689
#
# $
convergence# [1] 0# # $iterations# [1] 51# # $evaluations# function gradient # 87 341 # # $message# [1] "relative convergence (4)"Slide11
data(
speed_acc)
# Exp. 1; Wagenmakers, Ratcliff, Gomez, & McKoon
, 2008
# remove excluded trials:speed_acc <-
speed_acc
[!
speed_acc$censor
,]
speed_acc$corr
<- with(
speed_acc
, 1-as.numeric(
stim_cat
== response))
p11 <-
speed_acc
[
speed_acc$id
== 11 &
speed_acc$condition
== "accuracy" &
speed_acc$stim_cat
== "
nonword
",]
prop.table
(table(p11$corr
))
# 0 1 # 0.95833333 0.04166667
ll_lba <- function(pars,
rt, response) { ds <- split(rt
, f = response) d1 <- n1PDF(ds[[1]], A = pars["A"], b = pars["A"]+pars["b"], t0 = pars["t0"], mean_v = pars[c("v1", "v2")],
sd_v
= c(pars["
sv
"], 1), silent=TRUE)
d2 <- n1PDF(ds[[2]], A = pars["A"], b = pars["A"]+pars["b"], t0 = pars["t0"],
mean_v
= pars[c("v2", "v1")],
sd_v
= c(1, pars["
sv
"]), silent=TRUE)
d <- c(d1, d2)
if (any(d == 0)) return(1e6)
else return(-sum(log(d)))
}
start <- c(
runif
(3, 0.5, 3),
runif
(2, 0, 0.2),
rnorm
(1))
names(start) <- c("A", "v1", "v2", "b", "t0", "
sv
")
p11_norm <-
nlminb
(start,
ll_lba
, lower = c(0, 0, -
Inf
, 0, 0, 0),
rt
=p11$rt, response=p11$corr)
p11_norm
# $par
# A v1 v2 b t0
sv
# 0.1182946 1.0449834 -2.7409705 0.4513529 0.1243451 0.2609940
#
# $objective
# [1] -211.4202# # $convergence# [1] 0# # $iterations# [1] 65# # $evaluations# function gradient # 109 474 # # $message
# [1] "relative convergence (4)"Slide12
q <- c(0.1, 0.3, 0.5, 0.7, 0.9)
p11_q_c <- quantile(p11[p11$corr == 0, "
rt
"], probs = q)
# 10% 30% 50% 70% 90%
# 0.4900
0.5557 0.6060 0.6773 0.8231
p11_q_e <- quantile(p11[p11$corr == 1, "
rt
"],
probs
= q)
#
10% 30% 50% 70% 90%
# 0.4908
0.5391 0.5905 0.6413 1.0653
max_pred
<- n1CDF(
Inf
, A = p11_norm$par["A"], b = p11_norm$par["A"]+p11_norm$par["b"],
t0
= p11_norm$par["t0"],
mean_v
= c(p11_norm$par["v1"], p11_norm$par["v2"]),
sd_v
= c(p11_norm$par["
sv
"], 1))
# [
1] 0.9581341
inv_cdf
<- function(x, A, b, t0, mean_v
,
sd_v
, value) {
abs(value - n1CDF(x, A=A, b = b, t0 = t0,
mean_v
=
mean_v
,
sd_v
=
sd_v
, silent=TRUE))
}
pred_correct
<-
sapply
(q*
max_pred
, function(y) optimize(
inv_cdf
, c(0, 3), A = p11_norm$par["A"],
b
= p11_norm$par["A"]+p11_norm$par["b"], t0 = p11_norm$par["t0"],
mean_v
= c(p11_norm$par["v1"], p11_norm$par["v2"]),
sd_v
= c(p11_norm$par["
sv
"], 1),
value
= y)[[1
]])
# [
1] 0.4871693 0.5510531 0.6081732 0.6809982 0.8301352
pred_error
<-
sapply
(q*(1-max_pred), function(y) optimize(
inv_cdf
, c(0, 3), A = p11_norm$par["A"],
b = p11_norm$par["A"]+p11_norm$par["b"], t0 = p11_norm$par["t0"], mean_v = c(p11_norm$par["v2"], p11_norm$par["v1"]), sd_v = c(1, p11_norm$par["sv"]), value = y)[[1]])# [1] 0.4684436 0.5529506 0.6273629 0.7233882 0.9314874Slide13
q <- c(0.1, 0.3, 0.5, 0.7, 0.9)
plot(
pred_correct
, q*max_pred, type = "b",
ylim=c(0, 1), xlim = c(0.4, 1.3),
ylab
= "Cumulative Probability",
xlab
= "Response Time (sec)")
points(p11_q_c, q*
prop.table
(table(p11$corr))[1],
pch
= 2)
lines(
pred_error
, q*(1-max_pred), type = "b")
points(p11_q_e, q*
prop.table
(table(p11$corr))[2],
pch
= 2)
legend("right", legend = c("data", "LBA predictions"),
pch
= c(2, 1),
lty
= c(0, 1))Slide14
speed_acc$boundary
<- factor(speed_acc$corr
, labels = c("upper", "lower"))
ll_diffusion
<- function(pars, rt, boundary)
{
densities <-
tryCatch
(abs(
drd
(
rt
, boundary=boundary,
a=pars[1
], v=pars[2], t0=pars[3],
z=0.5,
sz
=pars[4], st0=pars[5],
sv
=pars[6
])
), error = function(e) 0)
if (any(densities == 0)) return(1e6)
return(-sum(log(densities)))
}
start <- c(
runif
(2, 0.5, 3), 0.1,
runif
(3, 0, 0.5))names(start) <- c("a", "v", "t0", "sz", "st0", "
sv")p11_fit <- nlminb
(start,
ll_diffusion
, lower = 0,
rt
=p11$rt, boundary=
as.character
(p11$boundary))
#
$par
# a v t0
sz
st0
sv
# 1.3203928 3.2746681 0.3385756 0.3504362 0.2019925 1.0567170
#
# $objective
# [1] -207.5513
#
# $convergence
# [1] 0
#
# $iterations
# [1] 38
#
# $evaluations
# function gradient
# 68 272
#
# $message
# [1] "relative convergence (4)"Slide15
max_pred
<-
do.call(
prd, args = c(t =
Inf, as.list(p11_fit$par), boundary = "upper"))
# [1] 0.9389375
i_prd
<- function(x,
args
, value, boundary) {
abs(value -
do.call
(
prd
,
args
= c(t = x,
args
, boundary = boundary)))
}
pred_correct
<-
sapply
(q*
max_pred
, function(y) optimize(
i_prd
, c(0, 3),
args
=
as.list(p11_fit$par), value = y, boundary = "upper")[[1
]])# [1] 0.4963202 0.5736658 0.6361402 0.7147534 0.8815223pred_error
<- sapply(q*(1-max_pred), function(y) optimize(i_prd, c(0, 3),
args
=
as.list
(p11_fit$par), value = y, boundary = "lower")[[1
]])
# [1] 0.4480115 0.5222731 0.5824374 0.6665744
0.8829667
plot(
pred_correct
, q*
max_pred
, type = "b",
ylim
=c(0, 1),
xlim
= c(0.4, 1.3),
ylab
= "Cumulative Probability",
xlab
= "Response Time (sec)")
points(p11_q_c, q*
prop.table
(table(p11$correct))[1],
pch
= 2)
lines(
pred_error
, q*(1-max_pred), type = "b")
points(p11_q_e, q*
prop.table
(table(p11$correct))[2],
pch
= 2)
legend("right", legend = c("data", "Diffusion predictions"),
pch
= c(2, 1),
lty
= c(0, 1
))Slide16
nlminb
(start,
ll_lba
, lower = c(0, 0, -
Inf
, 0, 0, 0),
rt
=p11$rt, response=p11$corr
)
# $par
# A v1 v2 b t0
sv
# 0.1182946 1.0449834 -2.7409705 0.4513529 0.1243451 0.2609940
# $objective
# [1]
-
211.4202
nlminb
(
start_frechet
,
ll_lba_frechet
, lower = 0,
rt
=p11$rt, response=p11$corr)
#
$par
# A v1 v2 b t0
sv
# 1.781754 2.279085 1.959484 1.007131 0.353327 5.536361
# $objective
# [1]
-174.0797
nlminb
(
start_gamma
,
ll_lba_gamma
, lower = 0,
rt
=p11$rt, response=p11$corr)
#
$par
# A v1 v2 b t0
sv
# 0.07309753 7.85791431 0.62676205 0.96593679 0.24999751 0.36519063
# $objective
# [1]
-
208.1622
nlminb
(start
,
ll_lba_lnorm
, lower = 0,
rt
=p11$rt, response=p11$corr
)
# $par
#
A v1 v2 b t0
sv
# 0.001157895 2.056824618 0.140372535 2.360040196 0.315294542 0.453494110
# $objective
# [1]
-204.272
Terry
,
Marley
,
Barnwal
,
Wagenmakers
,
Heathcote, & Brown (in press)Slide17
Summary
rtdists intended to provide reference implementation of two most common accumulator models in R: Diffusion model and LBAFunctions are fast and fully vectorized (allowing trialwise parameters)
Allows implementation of all major estimation methods.