In rtdists/rtdists: Response Time Distributions
This vignette provides the R scripts for a reanalysis of Experiment 1 of Ratcliff and Rouder (1998). In contrast to the original analysis, which used RT bins, we will employ trial-wise maximum likelihood estimation.
Description of the Experiment
In the experiment, three participants were asked to decide whether the overall brightness of pixel arrays displayed on a computer monitor was "high" or "low". To this end, the number of white versus black pixels (i.e., the brightness strength) was manipulated in 33 levels from 0% white pixels (level 0) to 100% white pixels (level 32). In addition, instruction manipulated speed and accuracy between blocks. In total, each participant contributed around XX trials per instruction condition.
The experiment contained another manipulation, the distribution (or brightness source) from which the pixel array was drawn. One distribution mean was on the "high" brightness side and one distribution mean was on the "low" brightness side. However, as the distributions were unbounded and overlapping, the same strength level could come from either distribution. Participant also received feedback whether or not they had picked the correct distribution (e.g., for the middle strength level 16 probability of belonging to either source was 50%). We do not further consider this manipulation in the following, which seems to be in line with the analysis of Ratcliff and Rouder (1998).
Descriptive data
As a first step, we load necessary packages, the data, and remove a few outliers. Then, we bin the strength levels to remove some of the noise in the extreme strength levels. For this, we construct five separate strength bins with approximately equal response behavior and comparable numbers of trials. This is similar to what was done originally by Ratcliff and Rouder (1998). The next table shows the number of trials per participant, bin, and response.
require(rtdists)
require(dplyr) # for data manipulations and looping
require(tidyr) # for data manipulations
require(lattice) # for plotting and corresponding themes
require(latticeExtra)
lattice.options(default.theme = standard.theme(color = FALSE))
lattice.options(default.args = list(as.table = TRUE))
options(digits = 3) # only three decimal digits
require(binom) # for binomial confidence intervals
data(rr98)
rr98 <- rr98[!rr98$outlier,] #remove outliers
#bins <- c(-0.5, 5.5, 10.5, 13.5, 16.5, 19.5, 25.5, 32.5) # seven bins like RR98
bins <- c(-0.5, 10.5, 13.5, 16.5, 19.5, 32.5)
rr98$strength_bin <- cut(rr98$strength, breaks = bins, include.lowest = TRUE)
levels(rr98$strength_bin) <- as.character(1:7)
# aggregate data for response probability plot:
agg_rr98_bin <- rr98 %>% group_by(id, instruction, strength_bin) %>%
summarise(n = n(),
dark = sum(response == "dark"),
light = sum(response == "light")) %>%
ungroup() %>%
mutate(prop = binom.confint(dark, n, methods = "agresti-coull")[,"mean"],
lower = binom.confint(dark, n, methods = "agresti-coull")$lower,
upper = binom.confint(dark, n, methods = "agresti-coull")$upper)
knitr::kable(
rr98 %>% group_by(id, instruction, strength_bin, response) %>%
summarise(n = n()) %>%
spread(strength_bin, n)
)
We then plot the probability with which each response (i.e., "dark" or "light") is given as a function of strength and instruction condition. This clearly shows that there is a massive effect of strength on the probability with which each response is given while at the same time the instruction only seems to have a minor effect and more on the extremes than in the middle.
xyplot(prop ~ strength_bin|id, agg_rr98_bin, group = instruction, type = "b",
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses")
Next, we want to get an overview of the response time distributions. For this we look at the response times of the five quantiles (i.e., 0.1, 0.3, 0.5/median, 0.7, 0.9) across the strength manipulations. This time, we also separate the plots by condition as the speed condition resulted in, as expected, vastly shorter response times. These two plots reveal considerable differences between the two instruction conditions.
quantiles <- c(0.1, 0.3, 0.5, 0.7, 0.9)
## aggregate data for quantile plot
quantiles_rr98_bin <- rr98 %>% group_by(id, instruction, strength_bin) %>%
do(as.data.frame(t(quantile(.$rt, probs = quantiles)))) %>%
ungroup() %>%
gather(quantile, rt, -id, -instruction, -strength_bin)
quantiles_rr98_bin$quantile <- factor(quantiles_rr98_bin$quantile,
levels = c("90%", "70%", "50%", "30%", "10%"))
xyplot(rt ~ strength_bin|id + instruction, quantiles_rr98_bin, group = quantile, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed")
xyplot(rt ~ strength_bin|id + instruction, quantiles_rr98_bin, group = quantile, type = "b",
auto.key = FALSE, ylab = "RT (in seconds)", subset = instruction == "accuracy")
In the speed conditions, response times were, as expected, generally fast and there seemed to be hardly any effect of strength. Only for one participant, nh, we can see a small increase in RTs for the higher quantiles for strength values near the middle. In contrast, in the accuracy condition strength has a considerable effect on response times for all participants. Again, this increase was especially strong for the slower responses (i.e., the higher quantiles). For those we see a strong inverse u-shaped effect, symmetrically around the middle -- where the probability for each response is 50% -- with very high response times for strength values near the middle.
With this clear pattern we now take a look at the RT distributions separately for both responses to see if they are simply mirror images of each other or not. For this, we overlay the two RT quantile plots for all trials in which the responses was "dark" in black (there are more "dark" pixels for the bins on the left side of the plot) with the same plot in which the responses was "light" in grey (there are more "light" pixels for the bins on the right side of the plot).
agg2_rr98_response <- rr98 %>% group_by(id, instruction, strength_bin, response) %>%
do(as.data.frame(t(quantile(.$rt, probs = c(0.1, 0.3, 0.5, 0.7, 0.9))))) %>%
ungroup() %>%
gather(quantile, rt, -id, - instruction, - strength_bin, -response)
agg2_rr98_response$quantile <- factor(agg2_rr98_response$quantile,
levels = c("90%", "70%", "50%", "30%", "10%"))
p1 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "speed" & response == "dark", layout = c(3,1))
p2 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "speed" & response == "light", col = "grey")
p1 + as.layer(p2)
p1 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "accuracy" & response == "dark", layout = c(3,1))
p2 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "accuracy" & response == "light", col = "grey")
p1 + as.layer(p2)
These two plots reveal an interesting pattern. In the speed condition (upper plot), we particularly see fast "errors" (i.e., responses to "dark" when there are more light pixels or the other way round). When "dark" is the more likely response (i.e. on the left side) the "light" responses in grey are faster and this is especially true for the lower quantiles. The opposite pattern seems to hold on the opposite side where "dark" responses in black are faster than "light" responses in grey. At intermediate bins the difference seems to be rather at the higher quantiles. This is particularly noticeable for participant kr for which for there seem to be slow "light"-"errors" just to the left to the middle bin and slow "right"-"errors" just to the right of the middle bin.
For the accuracy condition in the lower plot the pattern is noticeably different. First of all, there are only very few or no "error" responses in the extreme bins. Consequently, there does not seem to be any evidence for fast errors at the extremes (and also not at intermediate strength levels). However, we here more clearly see the slow errors at the intermediate bins. When "dark" is somewhat more probably (i.e., to the left of the middle) "light" responses are noticeably slower than "dark" responses. The same holds for "dark" responses if "light" is more probable. Importantly, this shows that the symmetrical inverse u-shaped increase for the middle bins described above is actually a consequence of a mixture of slow "errors", two asymmetric increases for the two different responses.
Wiener Model Analysis
We will analyze the data with the 4 parameter diffusion model also known as the Wiener model which has no trial-by-trial variabilities of the parameters. For this, we will fit a separate model to each participant and instruction condition. Like Ratcliff and Rouder (1998) we will fit the data to the strength bins instead of the full strength manipulation. Our model has 16 parameters per participant:
- 2 $\times$ 5 drift rates (i.e., one per strength bin) = 10 drift rates $v$
- 2 boundary separations $a$ (i.e., one per instruction condition)
- 2 non-decision times $t_0$
- 2 starting points $z$
The two boundaries are the two response options "dark" and "light". To estimate the model we diverge from Ratcliff and Rouder (1998) and employ trial wise maximum likelihood estimation (i.e., no binning of responses).
To do so, we simply need to have a wrapper function which returns us the summed log-likelihood of the data (i.e., RTs and corresponding responses) given a set of parameters. The following function for which we simply loop across drift rates will do so for the first model:
# objective function for diffusion with 1 a. loops over drift to assign drift rates to strength
objective_diffusion_separate <- function(pars, rt, boundary, drift, ...) {
non_v_pars <- grep("^v", names(pars), invert = TRUE, value = TRUE)
base_par <- length(non_v_pars) # number of non-drift parameters
densities <- vector("numeric", length(rt))
for (i in seq_along(levels(drift))) {
densities[drift == levels(drift)[i]] <- tryCatch(
ddiffusion(rt[drift == levels(drift)[i]], boundary=boundary[drift == levels(drift)[i]],
a=pars["a"], t0=pars["t0"],
z=if ("z" %in% non_v_pars) pars["z"] else 0.5,
v=pars[base_par+i]),
error = function(e) 0)
}
if (any(densities == 0)) return(1e6)
return(-sum(log(densities)))
}
As one can see at this function we also wrap the call to ddiffusion (i.e., the PDF of the diffusion model) into a tryCatch statement. This prevents failures during optimization as ddifussion can fail for some impossible parameters values. Note also that the function is written in such a way that we could easily fix certain parameters without the necessity to change it (using if-then on the parameters names passed via pars).
Additionally, we also need a function that generates a set of random starting values. And, as any random set of starting values may be impossible, another wrapper function that generates starting values until a set of valid starting values is found and then passes those to the optimization routine. As optimization routine we will be using nlminb. These functions are given next and are already specified in a way that it will be usable for the second
# function that creates random start values, also
get_start <- function(base_par, n_drift = 5) {
start1 <- c(
a = runif(1, 0.5, 3),
a_1 = runif(1, 0.5, 3),
a_2 = runif(1, 0.5, 3),
t0 = runif(1, 0, 0.5),
t0_1 = runif(1, 0, 0.5),
t0_2 = runif(1, 0, 0.5),
z = runif(1, 0.4, 0.6),
z = runif(1, 0.4, 0.6),
z = runif(1, 0.4, 0.6)
)
start2 <- sort(rnorm(n_drift), decreasing = FALSE)
names(start2) <- paste0("v_", seq_len(n_drift))
c(start1[base_par], start2)
}
# function that tries different random start values until it works:
ensure_fit <- function(data, start_function, objective_function, base_pars, n_drift = 5) {
start_ll <- 1e+06
#browser()
while(start_ll == 1e+06) {
start <- start_function(base_pars)
start_ll <- objective_function(start,
rt = data$rt, boundary = data$response_num,
drift = factor(data$strength_bin, seq_len(n_drift)),
instruction = data$instruction)
}
cat("\nstart fitting.\n") # just for information to see if it is stuck
fit <- nlminb(start, objective_function,
rt = data$rt, boundary = data$response_num,
drift = factor(data$strength_bin, seq_len(n_drift)),
instruction = data$instruction,
lower = c(rep(0, length(base_pars)), -Inf,
rep(-Inf, length(start_function(base_pars))-length(base_pars))))
fit
}
load("rr98_wiener_fits.rda")
With these functions in place, we now simply need to loop over participants and items to obtain the fit. We do this (as above) using convenient dplyr syntax and dplyr::do.
fits_separate <- rr98 %>%
group_by(id, instruction) %>% # we loop across both, id and instruction
do(diffusion = ensure_fit(data = ., start_function = get_start,
objective_function = objective_diffusion_separate,
base_pars = c("a", "t0", "z"))) %>% ungroup()
The following table gives the parameters and the negative summed log-likelihoods obtained from this fit:
pars_separate <- fits_separate %>% group_by(id, instruction) %>%
do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>% ungroup() %>%
as.data.frame()
pars_separate$ll <- (fits_separate %>% group_by(id, instruction) %>%
do(ll = .$diffusion[[1]][["objective"]]) %>%
summarize(ll2 = mean(ll[[1]])) %>% as.data.frame())[[1]]
if (!("z" %in% colnames(pars_separate))) pars_separate$z <- 0.5
knitr::kable(pars_separate)
We can see from these values that there is a large effect of instruction on $a$ as can be expected. However, instruction also has effects on most of the other parameters:
- $t_0$ is consistently larger in the accuracy compared to the speed condition, although this effect is small.
- The differences in $v$ are not consistent across participants, but also do not seem completely negligible.
- There does not seem to be a systematic effect on $z$, which hovers around 0.5.
fits_separate <- rr98 %>%
group_by(id, instruction) %>% # we loop across both, id and instruction
do(diffusion = ensure_fit(data = ., start_function = get_start,
objective_function = objective_diffusion_separate,
base_pars = c("a", "t0", "z"))) %>% ungroup()
fits_separate_b <- rr98 %>%
group_by(id, instruction) %>% # we loop across both, id and instruction
do(diffusion = ensure_fit(data = ., start_function = get_start,
objective_function = objective_diffusion_separate,
base_pars = c("a", "t0", "z"))) %>% ungroup()
pars_separate_b <- fits_separate_b %>% group_by(id, instruction) %>%
do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>% ungroup() %>%
as.data.frame()
pars_separate_b$ll <- (fits_separate_b %>% group_by(id, instruction) %>%
do(ll = .$diffusion[[1]][["objective"]]) %>%
summarize(ll2 = mean(ll[[1]])) %>% as.data.frame())[[1]]
all.equal(pars_separate, pars_separate_b, tolerance = 0.001)
fits_joint <- rr98 %>%
group_by(id) %>% # we loop only across id
do(diffusion = ensure_fit(data = ., start_function = get_start,
objective_function = objective_diffusion_joint,
base_pars = c("a_1", "a_2", "t0_1", "t0_2", "z"))) %>% ungroup()
fits_joint_b <- rr98 %>%
group_by(id) %>% # we loop only across id
do(diffusion = ensure_fit(data = ., start_function = get_start,
objective_function = objective_diffusion_joint,
base_pars = c("a_1", "a_2", "t0_1", "t0_2", "z"))) %>% ungroup()
pars_joint_b <- fits_joint_b %>% group_by(id) %>%
do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>%
ungroup() %>% as.data.frame()
pars_joint_b$ll <- (fits_joint_b %>% group_by(id) %>%
do(ll = .$diffusion[[1]][["objective"]]) %>%
summarize(ll2 = mean(ll[[1]])) %>% as.data.frame())[[1]]
all.equal(pars_joint, pars_joint_b, tolerance = 0.0001)
save(fits_separate, fits_separate_b, fits_joint, fits_joint_b, file = "rr98_wiener_fits.rda")
Graphical Model Fit
Predicted Response Rates
To evaluate the fits graphically we first compare the actual response rates for the two responses with the predicted responses rates. The grey lines and points show the observed data and the error bars are binomial confidence intervals. The black lines and points show the predicted response rates.
# get predicted response proportions
pars_separate_l <- pars_separate %>% gather("strength_bin", "v", starts_with("v"))
pars_separate_l$strength_bin <- factor(substr(pars_separate_l$strength_bin, 3,3),
levels = as.character(seq_len(length(bins)-1)))
#pars_separate_l <- inner_join(pars_separate_l, agg_rr98_bin)
pars_separate_l <- pars_separate_l %>% group_by(id, instruction, strength_bin) %>%
mutate(resp_prop = pdiffusion(rt=20, boundary="lower", a=a, v=v, t0=t0, z=z))
p1 <- xyplot(prop ~ strength_bin|id + instruction, agg_rr98_bin, type = "b", auto.key =
list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey")
p2 <- segplot(strength_bin ~ upper+lower|id + instruction, agg_rr98_bin,
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "grey", horizontal = FALSE, segments.fun = panel.arrows,
draw.bands = FALSE, angle = 90, length = 0.05, ends = "both")
p3 <- xyplot(resp_prop ~ strength_bin|id + instruction, pars_separate_l, type = "b",
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "black")
p2 + as.layer(p1) + as.layer(p3)
This figure show that overall the model can predict the actual response rates quite accurately. However, especially for the medium strengths and specifically in the accuracy condition the model mispredicts the rate of "dark" responses. In the speed conditions it slightly under-predicts and in the accuracy over-predicts the rate of "dark" responses.
Predicted Median RTs
Next we compare the central tendency of the RTs with the prediction. For this we evaluate the CDF at the quantiles of the predicted response proportions. Again, the data is shown in grey (and error bars show the standard errors of the median) and the predicted RTs in black. We first show the predictions for the speed condition, separated by response.
# get predicted quantiles (uses predicted response proportions)
separate_pred_dark <- pars_separate_l %>% do(as.data.frame(t(
qdiffusion(quantiles*.$resp_prop, boundary="lower",
a=.$a, v=.$v, t0=.$t0, z=.$z)))) %>%
ungroup() %>% gather("quantiles", "dark", V1:V5)
separate_pred_light <- pars_separate_l %>% do(as.data.frame(t(
qdiffusion(quantiles*(1-.$resp_prop), boundary="upper",
a=.$a, v=.$v, t0=.$t0, z=.$z)))) %>%
ungroup() %>% gather("quantiles", "light", V1:V5)
#separate_pred_light %>% filter(is.na(light))
separate_pred <- inner_join(separate_pred_dark, separate_pred_light)
separate_pred$quantiles <- factor(separate_pred$quantiles,
levels = c("V5", "V4", "V3", "V2", "V1"),
labels = c("90%", "70%", "50%", "30%", "10%"))
separate_pred <- separate_pred %>% gather("response", "rt", dark, light)
# get SE for observed quantiles
agg2_rr98_response_se <- rr98 %>% group_by(id, instruction, strength_bin, response) %>%
summarise(se_median = sqrt(pi/2)*(sd(rt)/sqrt(n()))) %>%
ungroup()
# calculate error bars for quantiles.
agg2_rr98_response <- left_join(agg2_rr98_response, agg2_rr98_response_se)
agg2_rr98_response <- agg2_rr98_response %>%
mutate(lower = rt-se_median, upper = rt+se_median)
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "speed" & quantile == "50%",
layout = c(3,2), col = "grey")
p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response,
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "grey", horizontal = FALSE, segments.fun = panel.arrows,
draw.bands = FALSE, angle = 90, length = 0.05, ends = "both",
subset = instruction == "speed" & quantile == "50%", layout = c(3,2))
p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "speed" & quantiles == "50%",
scales = list(y = list(limits = c(0.25, 0.5))))
p2 + as.layer(p1) + as.layer(p1e)
While the model seems to be overall able to describe the general pattern which is rather flat for two of the three participants, it has more problems for participant nh. Here we see a strong increase in median RT for the middle strength bin with misfits in excess of 50 ms.
Next shows the same plot for the accuracy condition. Here the misfit is even larger, there is almost no data point that is accurately described. Additionally, even the pattern seems to be missed for nh and light responses (albeit only for the bins with a few responses) and jf and light responses. The misfits here are a lot larger and often way above 250 ms.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "accuracy" & quantile == "50%",
layout = c(3,2), col = "grey")
p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response,
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "grey", horizontal = FALSE, segments.fun = panel.arrows,
draw.bands = FALSE, angle = 90, length = 0.05, ends = "both",
subset = instruction == "accuracy" & quantile == "50%", layout = c(3,2))
p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "accuracy" & quantiles == "50%",
scales = list(y = list(limits = c(0.2, 1.5))))
p2 + as.layer(p1) + as.layer(p1e)
All quantiles
Next, we investigate the full RT distribution by comparing observed and predicted quantiles. The observed quantiles are again displayed in grey and the predictions in black. What will be apparent easily is the fact that the model cannot make differential predictions for the two responses. Their pattern needs to be a mirror image of each other. Consequently, the asymmetric increase in RTs for the high quantiles cannot be described at all.
The first plot shows the speed condition separated by response.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, group = quantile, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "speed", layout = c(3,2), col = "grey")
p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response,
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "grey", horizontal = FALSE, segments.fun = panel.arrows,
draw.bands = FALSE, angle = 90, length = 0.05, ends = "both",
subset = instruction == "speed")
p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, group = quantiles, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "speed", scales = list(y = list(limits = c(0.2, 0.9))))
p2 + as.layer(p1) + as.layer(p1e)
The next plot shows the accuracy condition separated by response.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, group = quantile, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "accuracy", layout = c(3,2), col = "grey")
p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response,
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "grey", horizontal = FALSE, segments.fun = panel.arrows,
draw.bands = FALSE, angle = 90, length = 0.05, ends = "both",
subset = instruction == "accuracy")
p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, group = quantiles, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "accuracy", scales = list(y = list(limits = c(0.1, 3.0))))
p2 + as.layer(p1) + as.layer(p1e)
Besides revealing the problems in accounting for the slow errors, these two plots show that the somewhat meager performance for the central tendency, the median, extends to the other quantiles. However, while there clearly are problem with the Wiener model in accounting for the data, the predictions are still near the actual observations and the predicted patterns in line with what is observed.
Restricted Model
As shown above, the parameters seem to suggest that there are differences between the two instructions in the $a$ and $t_0$ parameters. Consequently, we fit a joint (i.e., reduced) model in which the other parameters (i.e., $v$ and $z$) are shared across instruction conditions and only $a$ and $t_0$ differ. This model has 10 parameters per participant:
- 5 drift rates (i.e., one per strength bin)
- 2 boundary separations $a$
- 2 non-decision time $t_0$
- 1 start point $z$.
For the joint model we need a new objective function that also loops across instruction. Then we use again a dplyr call to loop across participants to get the fit.
objective_diffusion_joint <- function(pars, rt, boundary, drift, instruction) {
non_v_pars <- grep("^v", names(pars), invert = TRUE, value = TRUE)
base_par <- length(non_v_pars) # number of non-drift parameters
densities <- vector("numeric", length(rt))
as <- c(pars["a_1"], pars["a_2"])
ts <- c(pars["t0_1"], pars["t0_2"])
for (j in seq_along(levels(instruction))) {
for (i in seq_along(levels(drift))) {
densities[drift == levels(drift)[i] & instruction == levels(instruction)[j]] <- tryCatch(
ddiffusion(rt[drift == levels(drift)[i] & instruction == levels(instruction)[j]],
boundary=boundary[drift==levels(drift)[i]&instruction==levels(instruction)[j]],
a=as[j], t0=ts[j], z=pars["z"],
v=pars[base_par+i]),
error = function(e) 0)
}
}
if (any(densities == 0)) return(1e6)
return(-sum(log(densities)))
}
fits_joint <- rr98 %>%
group_by(id) %>% # we loop only across id
do(diffusion = ensure_fit(data = ., start_function = get_start,
objective_function = objective_diffusion_joint,
base_pars = c("a_1", "a_2", "t0_1", "t0_2", "z"))) %>% ungroup()
The following table gives the parameters and the negative summed log-likelihoods obtained from this fit:
pars_joint <- fits_joint %>% group_by(id) %>%
do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>%
ungroup() %>% as.data.frame()
pars_joint$ll <- (fits_joint %>% group_by(id) %>%
do(ll = .$diffusion[[1]][["objective"]]) %>%
summarize(ll2 = mean(ll[[1]])) %>% as.data.frame())[[1]]
knitr::kable(pars_joint)
Compare Model Fits
To compare the fits between the two model version we simply perform a likelihood ratio test between the separate and the joint/restricted model. For this we first need to sum the two log-likelihoods of the separate model. Then we simply take their difference and can get a p-value for the likelihood ratio tests with 6 degrees of freedom (i.e., their difference in the number of parameters). As we can see the test is highly significant. However, given the large number of trials (between 7581 and 8532) this might perhaps be not too surprising. Consequently, we also perform a graphical fit assessment for the restricted model.
ll_tab <- left_join(pars_joint[,c("id", "ll")],
pars_separate %>% group_by(id) %>% summarise(ll_sep = sum(ll))) %>%
mutate(ll_diff_2 = 2*(ll-ll_sep),
p = round(pchisq(ll_diff_2, df = 5, lower.tail = FALSE), 4))
#rr98 %>% group_by(id) %>% summarise(n())
knitr::kable(ll_tab)
Graphical Model Fit
Predicted Response Rates
First we again assess the predicted response proportions. As before, the grey lines and points show the observed data and the error bars are binomial confidence intervals and the black lines and points show the predicted response rates of the restricted model. In addition, the dark grey dashed lines with squares show the unrestricted (i.e., separate) model fitted before.
# get predicted response proportions
pars_joint_l <- pars_joint %>% gather("strength_bin", "v", starts_with("v")) %>%
gather("instruction_tmp", "pars", a_1, a_2, t0_1, t0_2) %>%
separate(instruction_tmp, c("para", "instruction")) %>%
spread(para, pars)
pars_joint_l$instruction <- factor(pars_joint_l$instruction,
levels = c("1", "2"), labels = c("speed", "accuracy"))
pars_joint_l$strength_bin <- factor(substr(pars_joint_l$strength_bin, 3,3),
levels = as.character(seq_len(length(bins)-1)))
#pars_separate_l <- inner_join(pars_separate_l, agg_rr98_bin)
pars_joint_l <- pars_joint_l %>% group_by(id, instruction, strength_bin) %>%
mutate(resp_prop = pdiffusion(rt=20, boundary="lower", a=a, v=v, t0=t0, z=z))
p1 <- xyplot(prop ~ strength_bin|id + instruction, agg_rr98_bin, type = "b", auto.key =
list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey")
p1e <- segplot(strength_bin ~ upper+lower|id + instruction, agg_rr98_bin,
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "grey", horizontal = FALSE, segments.fun = panel.arrows,
draw.bands = FALSE, angle = 90, length = 0.05, ends = "both")
p2 <- xyplot(resp_prop ~ strength_bin|id + instruction, pars_separate_l, type = "b",
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "darkgrey", lty = 3, pch = 0)
p3 <- xyplot(resp_prop ~ strength_bin|id + instruction, pars_joint_l, type = "b",
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "black")
p2 + as.layer(p3) + as.layer(p1) + as.layer(p1e)
We see that the joint model makes predictions that are overall very similar to the ones from the separate fits, especially for the accuracy condition. In the speed condition the predictions are somewhat different, but not necessarily a lot worse.
Predicted Median RTs
Next we compare the central tendency of the RTs with the prediction. Again, the data is shown in grey (and error bars show the standard errors of the median) and the predicted RTs in black and the fit of the unrestricted model in dark grey with squares and a dashed line. We first show the predictions for the speed condition, separated by response.
# get predicted quantiles (uses predicted response proportions)
joint_pred_dark <- pars_joint_l %>% do(as.data.frame(t(
qdiffusion(quantiles*.$resp_prop, boundary="lower",
a=.$a, v=.$v, t0=.$t0, z=.$z)))) %>%
ungroup() %>% gather("quantiles", "dark", V1:V5)
joint_pred_light <- pars_joint_l %>% do(as.data.frame(t(
qdiffusion(quantiles*(1-.$resp_prop), boundary="upper",
a=.$a, v=.$v, t0=.$t0, z=.$z)))) %>%
ungroup() %>% gather("quantiles", "light", V1:V5)
#joint_pred_light %>% filter(is.na(light))
joint_pred <- inner_join(joint_pred_dark, joint_pred_light)
joint_pred$quantiles <- factor(joint_pred$quantiles,
levels = c("V5", "V4", "V3", "V2", "V1"),
labels = c("90%", "70%", "50%", "30%", "10%"))
joint_pred <- joint_pred %>% gather("response", "rt", dark, light)
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "speed" & quantile == "50%",
layout = c(3,2), col = "grey")
p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response,
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "grey", horizontal = FALSE, segments.fun = panel.arrows,
draw.bands = FALSE, angle = 90, length = 0.05, ends = "both",
subset = instruction == "speed" & quantile == "50%", layout = c(3,2))
p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "speed" & quantiles == "50%",
scales = list(y = list(limits = c(0.25, 0.5))),
col = "darkgrey", lty = 3, pch = 0)
p3 <- xyplot(rt ~ strength_bin|id + response, joint_pred, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "speed" & quantiles == "50%",
scales = list(y = list(limits = c(0.25, 0.5))))
p2 + as.layer(p3) + as.layer(p1) + as.layer(p1e)
Now we show the predictions for the accuracy condition, separated by response.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "accuracy" & quantile == "50%",
layout = c(3,2), col = "grey")
p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response,
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "grey", horizontal = FALSE, segments.fun = panel.arrows,
draw.bands = FALSE, angle = 90, length = 0.05, ends = "both",
subset = instruction == "accuracy" & quantile == "50%", layout = c(3,2))
p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "accuracy" & quantiles == "50%",
scales = list(y = list(limits = c(0.2, 1.5))),
col = "darkgrey", lty = 3, pch = 0)
p3 <- xyplot(rt ~ strength_bin|id + response, joint_pred, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "accuracy" & quantiles == "50%",
scales = list(y = list(limits = c(0.25, 0.5))))
p2 + as.layer(p3) + as.layer(p1) + as.layer(p1e)
We can see that the unrestricted model provides a somewhat better fit, especially in the speed condition, but given the considerable misfit of both model this seems rather negligible.
All quantiles
Finally, we also take a look at the full RT distribution by comparing observed and predicted quantiles. However, as the plots are already busy with one model, we now only show the prediction of the restricted model and the observed data. The observed quantiles are again displayed in grey and the predictions in black.
The first plot shows the speed condition separated by response.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, group = quantile, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "speed", layout = c(3,2), col = "grey")
p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response,
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "grey", horizontal = FALSE, segments.fun = panel.arrows,
draw.bands = FALSE, angle = 90, length = 0.05, ends = "both",
subset = instruction == "speed")
p2 <- xyplot(rt ~ strength_bin|id + response, joint_pred, group = quantiles, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "speed", scales = list(y = list(limits = c(0.2, 0.9))))
p2 + as.layer(p1) + as.layer(p1e)
The next plot shows the accuracy condition separated by response.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, group = quantile, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "accuracy", layout = c(3,2), col = "grey")
p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response,
auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses",
col = "grey", horizontal = FALSE, segments.fun = panel.arrows,
draw.bands = FALSE, angle = 90, length = 0.05, ends = "both",
subset = instruction == "accuracy")
p2 <- xyplot(rt ~ strength_bin|id + response, joint_pred, group = quantiles, type = "b",
auto.key = list(lines = TRUE), ylab = "RT (in seconds)",
subset = instruction == "accuracy", scales = list(y = list(limits = c(0.1, 3.0))))
p2 + as.layer(p1) + as.layer(p1e)
We can again see quite some misfit, but it does not appear to be dramatically different from the one for the restricted model.
References
- Ratcliff, R., & Rouder, J. N. (1998). Modeling Response Times for Two-Choice Decisions. Psychological Science, 9(5), 347--356. http://doi.org/10.1111/1467-9280.00067
rtdists/rtdists documentation built on Jan. 6, 2022, 2:31 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
This vignette provides the R scripts for a reanalysis of Experiment 1 of Ratcliff and Rouder (1998). In contrast to the original analysis, which used RT bins, we will employ trial-wise maximum likelihood estimation.
Description of the Experiment
In the experiment, three participants were asked to decide whether the overall brightness of pixel arrays displayed on a computer monitor was "high" or "low". To this end, the number of white versus black pixels (i.e., the brightness strength) was manipulated in 33 levels from 0% white pixels (level 0) to 100% white pixels (level 32). In addition, instruction manipulated speed and accuracy between blocks. In total, each participant contributed around XX trials per instruction condition.
The experiment contained another manipulation, the distribution (or brightness source) from which the pixel array was drawn. One distribution mean was on the "high" brightness side and one distribution mean was on the "low" brightness side. However, as the distributions were unbounded and overlapping, the same strength level could come from either distribution. Participant also received feedback whether or not they had picked the correct distribution (e.g., for the middle strength level 16 probability of belonging to either source was 50%). We do not further consider this manipulation in the following, which seems to be in line with the analysis of Ratcliff and Rouder (1998).
Descriptive data
As a first step, we load necessary packages, the data, and remove a few outliers. Then, we bin the strength levels to remove some of the noise in the extreme strength levels. For this, we construct five separate strength bins with approximately equal response behavior and comparable numbers of trials. This is similar to what was done originally by Ratcliff and Rouder (1998). The next table shows the number of trials per participant, bin, and response.
require(rtdists) require(dplyr) # for data manipulations and looping require(tidyr) # for data manipulations require(lattice) # for plotting and corresponding themes require(latticeExtra) lattice.options(default.theme = standard.theme(color = FALSE)) lattice.options(default.args = list(as.table = TRUE)) options(digits = 3) # only three decimal digits require(binom) # for binomial confidence intervals data(rr98) rr98 <- rr98[!rr98$outlier,] #remove outliers #bins <- c(-0.5, 5.5, 10.5, 13.5, 16.5, 19.5, 25.5, 32.5) # seven bins like RR98 bins <- c(-0.5, 10.5, 13.5, 16.5, 19.5, 32.5) rr98$strength_bin <- cut(rr98$strength, breaks = bins, include.lowest = TRUE) levels(rr98$strength_bin) <- as.character(1:7) # aggregate data for response probability plot: agg_rr98_bin <- rr98 %>% group_by(id, instruction, strength_bin) %>% summarise(n = n(), dark = sum(response == "dark"), light = sum(response == "light")) %>% ungroup() %>% mutate(prop = binom.confint(dark, n, methods = "agresti-coull")[,"mean"], lower = binom.confint(dark, n, methods = "agresti-coull")$lower, upper = binom.confint(dark, n, methods = "agresti-coull")$upper) knitr::kable( rr98 %>% group_by(id, instruction, strength_bin, response) %>% summarise(n = n()) %>% spread(strength_bin, n) )
We then plot the probability with which each response (i.e., "dark" or "light") is given as a function of strength and instruction condition. This clearly shows that there is a massive effect of strength on the probability with which each response is given while at the same time the instruction only seems to have a minor effect and more on the extremes than in the middle.
xyplot(prop ~ strength_bin|id, agg_rr98_bin, group = instruction, type = "b", auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses")
Next, we want to get an overview of the response time distributions. For this we look at the response times of the five quantiles (i.e., 0.1, 0.3, 0.5/median, 0.7, 0.9) across the strength manipulations. This time, we also separate the plots by condition as the speed condition resulted in, as expected, vastly shorter response times. These two plots reveal considerable differences between the two instruction conditions.
quantiles <- c(0.1, 0.3, 0.5, 0.7, 0.9) ## aggregate data for quantile plot quantiles_rr98_bin <- rr98 %>% group_by(id, instruction, strength_bin) %>% do(as.data.frame(t(quantile(.$rt, probs = quantiles)))) %>% ungroup() %>% gather(quantile, rt, -id, -instruction, -strength_bin) quantiles_rr98_bin$quantile <- factor(quantiles_rr98_bin$quantile, levels = c("90%", "70%", "50%", "30%", "10%")) xyplot(rt ~ strength_bin|id + instruction, quantiles_rr98_bin, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed") xyplot(rt ~ strength_bin|id + instruction, quantiles_rr98_bin, group = quantile, type = "b", auto.key = FALSE, ylab = "RT (in seconds)", subset = instruction == "accuracy")
In the speed conditions, response times were, as expected, generally fast and there seemed to be hardly any effect of strength. Only for one participant, nh, we can see a small increase in RTs for the higher quantiles for strength values near the middle. In contrast, in the accuracy condition strength has a considerable effect on response times for all participants. Again, this increase was especially strong for the slower responses (i.e., the higher quantiles). For those we see a strong inverse u-shaped effect, symmetrically around the middle -- where the probability for each response is 50% -- with very high response times for strength values near the middle.
With this clear pattern we now take a look at the RT distributions separately for both responses to see if they are simply mirror images of each other or not. For this, we overlay the two RT quantile plots for all trials in which the responses was "dark" in black (there are more "dark" pixels for the bins on the left side of the plot) with the same plot in which the responses was "light" in grey (there are more "light" pixels for the bins on the right side of the plot).
agg2_rr98_response <- rr98 %>% group_by(id, instruction, strength_bin, response) %>% do(as.data.frame(t(quantile(.$rt, probs = c(0.1, 0.3, 0.5, 0.7, 0.9))))) %>% ungroup() %>% gather(quantile, rt, -id, - instruction, - strength_bin, -response) agg2_rr98_response$quantile <- factor(agg2_rr98_response$quantile, levels = c("90%", "70%", "50%", "30%", "10%")) p1 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed" & response == "dark", layout = c(3,1)) p2 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed" & response == "light", col = "grey") p1 + as.layer(p2) p1 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "accuracy" & response == "dark", layout = c(3,1)) p2 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "accuracy" & response == "light", col = "grey") p1 + as.layer(p2)
These two plots reveal an interesting pattern. In the speed condition (upper plot), we particularly see fast "errors" (i.e., responses to "dark" when there are more light pixels or the other way round). When "dark" is the more likely response (i.e. on the left side) the "light" responses in grey are faster and this is especially true for the lower quantiles. The opposite pattern seems to hold on the opposite side where "dark" responses in black are faster than "light" responses in grey. At intermediate bins the difference seems to be rather at the higher quantiles. This is particularly noticeable for participant kr for which for there seem to be slow "light"-"errors" just to the left to the middle bin and slow "right"-"errors" just to the right of the middle bin.
For the accuracy condition in the lower plot the pattern is noticeably different. First of all, there are only very few or no "error" responses in the extreme bins. Consequently, there does not seem to be any evidence for fast errors at the extremes (and also not at intermediate strength levels). However, we here more clearly see the slow errors at the intermediate bins. When "dark" is somewhat more probably (i.e., to the left of the middle) "light" responses are noticeably slower than "dark" responses. The same holds for "dark" responses if "light" is more probable. Importantly, this shows that the symmetrical inverse u-shaped increase for the middle bins described above is actually a consequence of a mixture of slow "errors", two asymmetric increases for the two different responses.
Wiener Model Analysis
We will analyze the data with the 4 parameter diffusion model also known as the Wiener model which has no trial-by-trial variabilities of the parameters. For this, we will fit a separate model to each participant and instruction condition. Like Ratcliff and Rouder (1998) we will fit the data to the strength bins instead of the full strength manipulation. Our model has 16 parameters per participant:
- 2 $\times$ 5 drift rates (i.e., one per strength bin) = 10 drift rates $v$
- 2 boundary separations $a$ (i.e., one per instruction condition)
- 2 non-decision times $t_0$
- 2 starting points $z$
The two boundaries are the two response options "dark" and "light". To estimate the model we diverge from Ratcliff and Rouder (1998) and employ trial wise maximum likelihood estimation (i.e., no binning of responses).
To do so, we simply need to have a wrapper function which returns us the summed log-likelihood of the data (i.e., RTs and corresponding responses) given a set of parameters. The following function for which we simply loop across drift rates will do so for the first model:
# objective function for diffusion with 1 a. loops over drift to assign drift rates to strength objective_diffusion_separate <- function(pars, rt, boundary, drift, ...) { non_v_pars <- grep("^v", names(pars), invert = TRUE, value = TRUE) base_par <- length(non_v_pars) # number of non-drift parameters densities <- vector("numeric", length(rt)) for (i in seq_along(levels(drift))) { densities[drift == levels(drift)[i]] <- tryCatch( ddiffusion(rt[drift == levels(drift)[i]], boundary=boundary[drift == levels(drift)[i]], a=pars["a"], t0=pars["t0"], z=if ("z" %in% non_v_pars) pars["z"] else 0.5, v=pars[base_par+i]), error = function(e) 0) } if (any(densities == 0)) return(1e6) return(-sum(log(densities))) }
As one can see at this function we also wrap the call to ddiffusion (i.e., the PDF of the diffusion model) into a tryCatch statement. This prevents failures during optimization as ddifussion can fail for some impossible parameters values. Note also that the function is written in such a way that we could easily fix certain parameters without the necessity to change it (using if-then on the parameters names passed via pars).
Additionally, we also need a function that generates a set of random starting values. And, as any random set of starting values may be impossible, another wrapper function that generates starting values until a set of valid starting values is found and then passes those to the optimization routine. As optimization routine we will be using nlminb. These functions are given next and are already specified in a way that it will be usable for the second
# function that creates random start values, also get_start <- function(base_par, n_drift = 5) { start1 <- c( a = runif(1, 0.5, 3), a_1 = runif(1, 0.5, 3), a_2 = runif(1, 0.5, 3), t0 = runif(1, 0, 0.5), t0_1 = runif(1, 0, 0.5), t0_2 = runif(1, 0, 0.5), z = runif(1, 0.4, 0.6), z = runif(1, 0.4, 0.6), z = runif(1, 0.4, 0.6) ) start2 <- sort(rnorm(n_drift), decreasing = FALSE) names(start2) <- paste0("v_", seq_len(n_drift)) c(start1[base_par], start2) } # function that tries different random start values until it works: ensure_fit <- function(data, start_function, objective_function, base_pars, n_drift = 5) { start_ll <- 1e+06 #browser() while(start_ll == 1e+06) { start <- start_function(base_pars) start_ll <- objective_function(start, rt = data$rt, boundary = data$response_num, drift = factor(data$strength_bin, seq_len(n_drift)), instruction = data$instruction) } cat("\nstart fitting.\n") # just for information to see if it is stuck fit <- nlminb(start, objective_function, rt = data$rt, boundary = data$response_num, drift = factor(data$strength_bin, seq_len(n_drift)), instruction = data$instruction, lower = c(rep(0, length(base_pars)), -Inf, rep(-Inf, length(start_function(base_pars))-length(base_pars)))) fit }
load("rr98_wiener_fits.rda")
With these functions in place, we now simply need to loop over participants and items to obtain the fit. We do this (as above) using convenient dplyr syntax and dplyr::do.
fits_separate <- rr98 %>% group_by(id, instruction) %>% # we loop across both, id and instruction do(diffusion = ensure_fit(data = ., start_function = get_start, objective_function = objective_diffusion_separate, base_pars = c("a", "t0", "z"))) %>% ungroup()
The following table gives the parameters and the negative summed log-likelihoods obtained from this fit:
pars_separate <- fits_separate %>% group_by(id, instruction) %>% do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>% ungroup() %>% as.data.frame() pars_separate$ll <- (fits_separate %>% group_by(id, instruction) %>% do(ll = .$diffusion[[1]][["objective"]]) %>% summarize(ll2 = mean(ll[[1]])) %>% as.data.frame())[[1]] if (!("z" %in% colnames(pars_separate))) pars_separate$z <- 0.5 knitr::kable(pars_separate)
We can see from these values that there is a large effect of instruction on $a$ as can be expected. However, instruction also has effects on most of the other parameters:
- $t_0$ is consistently larger in the accuracy compared to the speed condition, although this effect is small.
- The differences in $v$ are not consistent across participants, but also do not seem completely negligible.
- There does not seem to be a systematic effect on $z$, which hovers around 0.5.
fits_separate <- rr98 %>% group_by(id, instruction) %>% # we loop across both, id and instruction do(diffusion = ensure_fit(data = ., start_function = get_start, objective_function = objective_diffusion_separate, base_pars = c("a", "t0", "z"))) %>% ungroup() fits_separate_b <- rr98 %>% group_by(id, instruction) %>% # we loop across both, id and instruction do(diffusion = ensure_fit(data = ., start_function = get_start, objective_function = objective_diffusion_separate, base_pars = c("a", "t0", "z"))) %>% ungroup() pars_separate_b <- fits_separate_b %>% group_by(id, instruction) %>% do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>% ungroup() %>% as.data.frame() pars_separate_b$ll <- (fits_separate_b %>% group_by(id, instruction) %>% do(ll = .$diffusion[[1]][["objective"]]) %>% summarize(ll2 = mean(ll[[1]])) %>% as.data.frame())[[1]] all.equal(pars_separate, pars_separate_b, tolerance = 0.001) fits_joint <- rr98 %>% group_by(id) %>% # we loop only across id do(diffusion = ensure_fit(data = ., start_function = get_start, objective_function = objective_diffusion_joint, base_pars = c("a_1", "a_2", "t0_1", "t0_2", "z"))) %>% ungroup() fits_joint_b <- rr98 %>% group_by(id) %>% # we loop only across id do(diffusion = ensure_fit(data = ., start_function = get_start, objective_function = objective_diffusion_joint, base_pars = c("a_1", "a_2", "t0_1", "t0_2", "z"))) %>% ungroup() pars_joint_b <- fits_joint_b %>% group_by(id) %>% do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>% ungroup() %>% as.data.frame() pars_joint_b$ll <- (fits_joint_b %>% group_by(id) %>% do(ll = .$diffusion[[1]][["objective"]]) %>% summarize(ll2 = mean(ll[[1]])) %>% as.data.frame())[[1]] all.equal(pars_joint, pars_joint_b, tolerance = 0.0001) save(fits_separate, fits_separate_b, fits_joint, fits_joint_b, file = "rr98_wiener_fits.rda")
Graphical Model Fit
Predicted Response Rates
To evaluate the fits graphically we first compare the actual response rates for the two responses with the predicted responses rates. The grey lines and points show the observed data and the error bars are binomial confidence intervals. The black lines and points show the predicted response rates.
# get predicted response proportions pars_separate_l <- pars_separate %>% gather("strength_bin", "v", starts_with("v")) pars_separate_l$strength_bin <- factor(substr(pars_separate_l$strength_bin, 3,3), levels = as.character(seq_len(length(bins)-1))) #pars_separate_l <- inner_join(pars_separate_l, agg_rr98_bin) pars_separate_l <- pars_separate_l %>% group_by(id, instruction, strength_bin) %>% mutate(resp_prop = pdiffusion(rt=20, boundary="lower", a=a, v=v, t0=t0, z=z)) p1 <- xyplot(prop ~ strength_bin|id + instruction, agg_rr98_bin, type = "b", auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey") p2 <- segplot(strength_bin ~ upper+lower|id + instruction, agg_rr98_bin, auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey", horizontal = FALSE, segments.fun = panel.arrows, draw.bands = FALSE, angle = 90, length = 0.05, ends = "both") p3 <- xyplot(resp_prop ~ strength_bin|id + instruction, pars_separate_l, type = "b", auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "black") p2 + as.layer(p1) + as.layer(p3)
This figure show that overall the model can predict the actual response rates quite accurately. However, especially for the medium strengths and specifically in the accuracy condition the model mispredicts the rate of "dark" responses. In the speed conditions it slightly under-predicts and in the accuracy over-predicts the rate of "dark" responses.
Predicted Median RTs
Next we compare the central tendency of the RTs with the prediction. For this we evaluate the CDF at the quantiles of the predicted response proportions. Again, the data is shown in grey (and error bars show the standard errors of the median) and the predicted RTs in black. We first show the predictions for the speed condition, separated by response.
# get predicted quantiles (uses predicted response proportions) separate_pred_dark <- pars_separate_l %>% do(as.data.frame(t( qdiffusion(quantiles*.$resp_prop, boundary="lower", a=.$a, v=.$v, t0=.$t0, z=.$z)))) %>% ungroup() %>% gather("quantiles", "dark", V1:V5) separate_pred_light <- pars_separate_l %>% do(as.data.frame(t( qdiffusion(quantiles*(1-.$resp_prop), boundary="upper", a=.$a, v=.$v, t0=.$t0, z=.$z)))) %>% ungroup() %>% gather("quantiles", "light", V1:V5) #separate_pred_light %>% filter(is.na(light)) separate_pred <- inner_join(separate_pred_dark, separate_pred_light) separate_pred$quantiles <- factor(separate_pred$quantiles, levels = c("V5", "V4", "V3", "V2", "V1"), labels = c("90%", "70%", "50%", "30%", "10%")) separate_pred <- separate_pred %>% gather("response", "rt", dark, light) # get SE for observed quantiles agg2_rr98_response_se <- rr98 %>% group_by(id, instruction, strength_bin, response) %>% summarise(se_median = sqrt(pi/2)*(sd(rt)/sqrt(n()))) %>% ungroup() # calculate error bars for quantiles. agg2_rr98_response <- left_join(agg2_rr98_response, agg2_rr98_response_se) agg2_rr98_response <- agg2_rr98_response %>% mutate(lower = rt-se_median, upper = rt+se_median) p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed" & quantile == "50%", layout = c(3,2), col = "grey") p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response, auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey", horizontal = FALSE, segments.fun = panel.arrows, draw.bands = FALSE, angle = 90, length = 0.05, ends = "both", subset = instruction == "speed" & quantile == "50%", layout = c(3,2)) p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed" & quantiles == "50%", scales = list(y = list(limits = c(0.25, 0.5)))) p2 + as.layer(p1) + as.layer(p1e)
While the model seems to be overall able to describe the general pattern which is rather flat for two of the three participants, it has more problems for participant nh. Here we see a strong increase in median RT for the middle strength bin with misfits in excess of 50 ms.
Next shows the same plot for the accuracy condition. Here the misfit is even larger, there is almost no data point that is accurately described. Additionally, even the pattern seems to be missed for nh and light responses (albeit only for the bins with a few responses) and jf and light responses. The misfits here are a lot larger and often way above 250 ms.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "accuracy" & quantile == "50%", layout = c(3,2), col = "grey") p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response, auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey", horizontal = FALSE, segments.fun = panel.arrows, draw.bands = FALSE, angle = 90, length = 0.05, ends = "both", subset = instruction == "accuracy" & quantile == "50%", layout = c(3,2)) p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "accuracy" & quantiles == "50%", scales = list(y = list(limits = c(0.2, 1.5)))) p2 + as.layer(p1) + as.layer(p1e)
All quantiles
Next, we investigate the full RT distribution by comparing observed and predicted quantiles. The observed quantiles are again displayed in grey and the predictions in black. What will be apparent easily is the fact that the model cannot make differential predictions for the two responses. Their pattern needs to be a mirror image of each other. Consequently, the asymmetric increase in RTs for the high quantiles cannot be described at all.
The first plot shows the speed condition separated by response.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed", layout = c(3,2), col = "grey") p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response, auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey", horizontal = FALSE, segments.fun = panel.arrows, draw.bands = FALSE, angle = 90, length = 0.05, ends = "both", subset = instruction == "speed") p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, group = quantiles, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed", scales = list(y = list(limits = c(0.2, 0.9)))) p2 + as.layer(p1) + as.layer(p1e)
The next plot shows the accuracy condition separated by response.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "accuracy", layout = c(3,2), col = "grey") p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response, auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey", horizontal = FALSE, segments.fun = panel.arrows, draw.bands = FALSE, angle = 90, length = 0.05, ends = "both", subset = instruction == "accuracy") p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, group = quantiles, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "accuracy", scales = list(y = list(limits = c(0.1, 3.0)))) p2 + as.layer(p1) + as.layer(p1e)
Besides revealing the problems in accounting for the slow errors, these two plots show that the somewhat meager performance for the central tendency, the median, extends to the other quantiles. However, while there clearly are problem with the Wiener model in accounting for the data, the predictions are still near the actual observations and the predicted patterns in line with what is observed.
Restricted Model
As shown above, the parameters seem to suggest that there are differences between the two instructions in the $a$ and $t_0$ parameters. Consequently, we fit a joint (i.e., reduced) model in which the other parameters (i.e., $v$ and $z$) are shared across instruction conditions and only $a$ and $t_0$ differ. This model has 10 parameters per participant:
- 5 drift rates (i.e., one per strength bin)
- 2 boundary separations $a$
- 2 non-decision time $t_0$
- 1 start point $z$.
For the joint model we need a new objective function that also loops across instruction. Then we use again a dplyr call to loop across participants to get the fit.
objective_diffusion_joint <- function(pars, rt, boundary, drift, instruction) { non_v_pars <- grep("^v", names(pars), invert = TRUE, value = TRUE) base_par <- length(non_v_pars) # number of non-drift parameters densities <- vector("numeric", length(rt)) as <- c(pars["a_1"], pars["a_2"]) ts <- c(pars["t0_1"], pars["t0_2"]) for (j in seq_along(levels(instruction))) { for (i in seq_along(levels(drift))) { densities[drift == levels(drift)[i] & instruction == levels(instruction)[j]] <- tryCatch( ddiffusion(rt[drift == levels(drift)[i] & instruction == levels(instruction)[j]], boundary=boundary[drift==levels(drift)[i]&instruction==levels(instruction)[j]], a=as[j], t0=ts[j], z=pars["z"], v=pars[base_par+i]), error = function(e) 0) } } if (any(densities == 0)) return(1e6) return(-sum(log(densities))) }
fits_joint <- rr98 %>% group_by(id) %>% # we loop only across id do(diffusion = ensure_fit(data = ., start_function = get_start, objective_function = objective_diffusion_joint, base_pars = c("a_1", "a_2", "t0_1", "t0_2", "z"))) %>% ungroup()
The following table gives the parameters and the negative summed log-likelihoods obtained from this fit:
pars_joint <- fits_joint %>% group_by(id) %>% do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>% ungroup() %>% as.data.frame() pars_joint$ll <- (fits_joint %>% group_by(id) %>% do(ll = .$diffusion[[1]][["objective"]]) %>% summarize(ll2 = mean(ll[[1]])) %>% as.data.frame())[[1]] knitr::kable(pars_joint)
Compare Model Fits
To compare the fits between the two model version we simply perform a likelihood ratio test between the separate and the joint/restricted model. For this we first need to sum the two log-likelihoods of the separate model. Then we simply take their difference and can get a p-value for the likelihood ratio tests with 6 degrees of freedom (i.e., their difference in the number of parameters). As we can see the test is highly significant. However, given the large number of trials (between 7581 and 8532) this might perhaps be not too surprising. Consequently, we also perform a graphical fit assessment for the restricted model.
ll_tab <- left_join(pars_joint[,c("id", "ll")], pars_separate %>% group_by(id) %>% summarise(ll_sep = sum(ll))) %>% mutate(ll_diff_2 = 2*(ll-ll_sep), p = round(pchisq(ll_diff_2, df = 5, lower.tail = FALSE), 4)) #rr98 %>% group_by(id) %>% summarise(n()) knitr::kable(ll_tab)
Graphical Model Fit
Predicted Response Rates
First we again assess the predicted response proportions. As before, the grey lines and points show the observed data and the error bars are binomial confidence intervals and the black lines and points show the predicted response rates of the restricted model. In addition, the dark grey dashed lines with squares show the unrestricted (i.e., separate) model fitted before.
# get predicted response proportions pars_joint_l <- pars_joint %>% gather("strength_bin", "v", starts_with("v")) %>% gather("instruction_tmp", "pars", a_1, a_2, t0_1, t0_2) %>% separate(instruction_tmp, c("para", "instruction")) %>% spread(para, pars) pars_joint_l$instruction <- factor(pars_joint_l$instruction, levels = c("1", "2"), labels = c("speed", "accuracy")) pars_joint_l$strength_bin <- factor(substr(pars_joint_l$strength_bin, 3,3), levels = as.character(seq_len(length(bins)-1))) #pars_separate_l <- inner_join(pars_separate_l, agg_rr98_bin) pars_joint_l <- pars_joint_l %>% group_by(id, instruction, strength_bin) %>% mutate(resp_prop = pdiffusion(rt=20, boundary="lower", a=a, v=v, t0=t0, z=z)) p1 <- xyplot(prop ~ strength_bin|id + instruction, agg_rr98_bin, type = "b", auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey") p1e <- segplot(strength_bin ~ upper+lower|id + instruction, agg_rr98_bin, auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey", horizontal = FALSE, segments.fun = panel.arrows, draw.bands = FALSE, angle = 90, length = 0.05, ends = "both") p2 <- xyplot(resp_prop ~ strength_bin|id + instruction, pars_separate_l, type = "b", auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "darkgrey", lty = 3, pch = 0) p3 <- xyplot(resp_prop ~ strength_bin|id + instruction, pars_joint_l, type = "b", auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "black") p2 + as.layer(p3) + as.layer(p1) + as.layer(p1e)
We see that the joint model makes predictions that are overall very similar to the ones from the separate fits, especially for the accuracy condition. In the speed condition the predictions are somewhat different, but not necessarily a lot worse.
Predicted Median RTs
Next we compare the central tendency of the RTs with the prediction. Again, the data is shown in grey (and error bars show the standard errors of the median) and the predicted RTs in black and the fit of the unrestricted model in dark grey with squares and a dashed line. We first show the predictions for the speed condition, separated by response.
# get predicted quantiles (uses predicted response proportions) joint_pred_dark <- pars_joint_l %>% do(as.data.frame(t( qdiffusion(quantiles*.$resp_prop, boundary="lower", a=.$a, v=.$v, t0=.$t0, z=.$z)))) %>% ungroup() %>% gather("quantiles", "dark", V1:V5) joint_pred_light <- pars_joint_l %>% do(as.data.frame(t( qdiffusion(quantiles*(1-.$resp_prop), boundary="upper", a=.$a, v=.$v, t0=.$t0, z=.$z)))) %>% ungroup() %>% gather("quantiles", "light", V1:V5) #joint_pred_light %>% filter(is.na(light)) joint_pred <- inner_join(joint_pred_dark, joint_pred_light) joint_pred$quantiles <- factor(joint_pred$quantiles, levels = c("V5", "V4", "V3", "V2", "V1"), labels = c("90%", "70%", "50%", "30%", "10%")) joint_pred <- joint_pred %>% gather("response", "rt", dark, light) p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed" & quantile == "50%", layout = c(3,2), col = "grey") p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response, auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey", horizontal = FALSE, segments.fun = panel.arrows, draw.bands = FALSE, angle = 90, length = 0.05, ends = "both", subset = instruction == "speed" & quantile == "50%", layout = c(3,2)) p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed" & quantiles == "50%", scales = list(y = list(limits = c(0.25, 0.5))), col = "darkgrey", lty = 3, pch = 0) p3 <- xyplot(rt ~ strength_bin|id + response, joint_pred, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed" & quantiles == "50%", scales = list(y = list(limits = c(0.25, 0.5)))) p2 + as.layer(p3) + as.layer(p1) + as.layer(p1e)
Now we show the predictions for the accuracy condition, separated by response.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "accuracy" & quantile == "50%", layout = c(3,2), col = "grey") p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response, auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey", horizontal = FALSE, segments.fun = panel.arrows, draw.bands = FALSE, angle = 90, length = 0.05, ends = "both", subset = instruction == "accuracy" & quantile == "50%", layout = c(3,2)) p2 <- xyplot(rt ~ strength_bin|id + response, separate_pred, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "accuracy" & quantiles == "50%", scales = list(y = list(limits = c(0.2, 1.5))), col = "darkgrey", lty = 3, pch = 0) p3 <- xyplot(rt ~ strength_bin|id + response, joint_pred, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "accuracy" & quantiles == "50%", scales = list(y = list(limits = c(0.25, 0.5)))) p2 + as.layer(p3) + as.layer(p1) + as.layer(p1e)
We can see that the unrestricted model provides a somewhat better fit, especially in the speed condition, but given the considerable misfit of both model this seems rather negligible.
All quantiles
Finally, we also take a look at the full RT distribution by comparing observed and predicted quantiles. However, as the plots are already busy with one model, we now only show the prediction of the restricted model and the observed data. The observed quantiles are again displayed in grey and the predictions in black.
The first plot shows the speed condition separated by response.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed", layout = c(3,2), col = "grey") p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response, auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey", horizontal = FALSE, segments.fun = panel.arrows, draw.bands = FALSE, angle = 90, length = 0.05, ends = "both", subset = instruction == "speed") p2 <- xyplot(rt ~ strength_bin|id + response, joint_pred, group = quantiles, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed", scales = list(y = list(limits = c(0.2, 0.9)))) p2 + as.layer(p1) + as.layer(p1e)
The next plot shows the accuracy condition separated by response.
p1 <- xyplot(rt ~ strength_bin|id+response, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "accuracy", layout = c(3,2), col = "grey") p1e <- segplot(strength_bin ~ upper+lower|id+response, agg2_rr98_response, auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses", col = "grey", horizontal = FALSE, segments.fun = panel.arrows, draw.bands = FALSE, angle = 90, length = 0.05, ends = "both", subset = instruction == "accuracy") p2 <- xyplot(rt ~ strength_bin|id + response, joint_pred, group = quantiles, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "accuracy", scales = list(y = list(limits = c(0.1, 3.0)))) p2 + as.layer(p1) + as.layer(p1e)
We can again see quite some misfit, but it does not appear to be dramatically different from the one for the restricted model.
References
- Ratcliff, R., & Rouder, J. N. (1998). Modeling Response Times for Two-Choice Decisions. Psychological Science, 9(5), 347--356. http://doi.org/10.1111/1467-9280.00067
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.