Due to NCP Day, the course on 19.12.2018 is reschedule to 09.01.2019.
RTs are left skewed with long tail on the right side. (Question: why?)
raw = readRDS('data/rt_raw.rds')
head(raw)
Let’s first visualize the RTs with historgram.
raw %>% filter(!outlier & !error) %>%
group_by(dimension, BlkType) %>%
ggplot(aes(x=rt)) + geom_histogram() + facet_grid(BlkType~dimension)
This kind of distribution can be predicted by the drift-diffusion model (DDM).
Reaction time can be viewed as the first passage time of the accumulated decision variable, which can be described by the drift-diffusion model, one of the most successful models of choice RT in cognitive psychology. DDM has 4 key parameters:
To analyze the RT distribution with diffusion model, we need R package RWiener
, which provides R functions for the Wiener diffusion model.
#install.packages("RWiener")
library(RWiener)
Let’s get familiar with the RWiener
package.
We use rwiener
to generate 100 samples with boundary separation of 2, non-decision time 0.3 second, initial unbias (0.5) and drift rate of 0.5. The returned data is a table with responses (resp
) either ‘upper’ or ‘lower’ (2AFC) and RTs (q
). Then, we use dwiener
to visualize the density function.
# randomly generate the RT data
dat = rwiener(n = 100, alpha = 2, tau = 0.3, beta = .5, delta = .5)
# visualize the density function
curve(dwiener(x, 2, .3, .5, .5, rep("upper", length(x))),
xlim=c(0,3), main="Density of upper responses",
ylab="density", xlab="quantile")
Or we can easily visualize from the data by using
wiener_plot
.
wiener_plot(dat)
In the above section, we specify parameters and simulate RT data. In reality, we need to estimate the parameters of the DDM. This can be done by using R’s optim function, in which we need to specify the initial guess parameter (\(\alpha, \tau, \beta, \delta\)), wiener_deviance
and data table.
# estimate parameters
est = optim(c(1,0.1,0.1,1), wiener_deviance, dat = dat)
# show parameters (true values: 2, 0.3, 0.5, 0.5)
est$par
## [1] 1.8581722 0.3040059 0.4767229 0.4639543
# Note the new version use wdm function.
Note that, the upper and lower are two possible responses during a drift-diffusion process. Here we use ‘color’ for the upper bound and the ‘orientation’ for the lower bound. This also means that the rate is not purely information intaking, rather the ratio of the information intaking speed. Alternatively, one can use ‘error’ and ‘correct’ for the lower and upper bound. Doing this you will get absolute information intaking speed, but this requires a certain amount of error trials. Otherwise, it is hard to find the lower boundary.
Let’s first estimate one condition from one participant.
# select subject data
raw %>% filter(!outlier & !error & BlkType == '3:1' & sub =='aleg') %>%
mutate(resp = ifelse(dimension=='Color', 'upper','lower'),
q = rt) %>% # add upper, lower
select(q,resp) -> sub1
# find parameter
optim(c(2,.1,.5,0.3), wiener_deviance, dat = sub1) -> est
# print estimated parameter
est$par
## [1] 1.0517691 0.2637321 0.4653969 1.3890350
Let’s estimate parameters for all conditions.
# we need a function for estimating parameter
estPar <- function(df){
est = optim(c(2,0.1,0.5,0.3), wiener_deviance, dat = data.frame(df))
data.frame(a = est$par[1], tau =est$par[2], beta = est$par[3], delta = est$par[4] ) # return parameter
}
# subsecting data, estimate parameter, then combine all parameters together
raw %>% filter(!outlier & !error & sub =='aleg' ) %>%
mutate(resp = ifelse(dimension == 'Color', 'upper','lower'), q = rt) %>%
select(q, resp, BlkType, sub) %>%
group_by(sub, BlkType) %>% nest() %>%
mutate(par = data %>% map(estPar)) %>%
select(sub, BlkType, par) %>%
unnest()