---
title: "Reanalysis of Ratcliff and Rouder (1998)"
author: "Henrik Singmann"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Reanalysis of Ratcliff and Rouder (1998)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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 the data and 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 which 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.

```{r, fig.height=4, fig.width=7}
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

data(rr98)
rr98 <- rr98[!rr98$outlier,]  #remove outliers

# aggregate data for first plot:
agg_rr98 <- rr98  %>% group_by(id, instruction, strength) %>% 
  summarise(prop = mean(response == "dark"), mean_rt = mean(rt), median_rt = mean(rt)) %>% 
  ungroup()

xyplot(prop ~ strength|id, agg_rr98, 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. 

```{r, fig.height=6, fig.width=7}

quantiles <- c(0.1, 0.3, 0.5, 0.7, 0.9)
## aggregate data for quantile plot
quantiles_rr98 <- rr98  %>% group_by(id, instruction, strength) %>% 
  do(as.data.frame(t(quantile(.$rt, probs = quantiles)))) %>%
  ungroup() %>%
  gather(quantile, rt,  -id, - instruction, - strength)
quantiles_rr98$quantile <- factor(quantiles_rr98$quantile, levels = c("90%", "70%", "50%", "30%", "10%"))

xyplot(rt ~ strength|id + instruction, quantiles_rr98, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed")

xyplot(rt ~ strength|id + instruction, quantiles_rr98, 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. 

However, as this plot is getting a little bit messy, let us next bin the strength levels to get a somewhat clearer overview and less noise. For this, we will make up seven separate strength bins with approximately equal response behavior. This is also what was done originally by Ratcliff and Rouder (1998). We first look again and the response proportions and see more clearly the difference between the strength conditions at the outer bins.


```{r, fig.height=4, fig.width=7}

bins <- c(-0.5, 5.5, 10.5, 13.5, 16.5, 19.5, 25.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 first plot:
agg_rr98_bin <- rr98  %>% group_by(id, instruction, strength_bin) %>% 
  summarise(prop = mean(response == "dark"), mean_rt = mean(rt), median_rt = mean(rt)) %>% 
  ungroup()

xyplot(prop ~ strength_bin|id, agg_rr98_bin, group = instruction, type = "b", auto.key = list(lines = TRUE), ylab = "Proportion of 'dark' responses")


```

Now we also look again at the RT quantiles and see the symmetrical inverse u-shaped increase in RTs for the middle bins described above more clearly.

```{r, fig.height=6, fig.width=7}

## 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")


```

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).

```{r, fig.height=6, fig.width=7}

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.


## Full Diffusion Model

In the first step of the analysis we will follow Ratcliff and Rouder (1998) and basically analyze the data with the full diffusion model. The analysis will be performed in two steps. First, we will fit a separate model to each participant and instruction condition. After that, we will fit a reduced joint-model to both conditions of each participant in which only the boundary separation `a` is allowed to differ between the two conditions and all other parameters are shared. Like Ratcliff and Rouder we will fit the data to the strength bins instead of the full strength manipulation (although the results do not differ between both approaches).

### Separate Fits for Instructions

The combined model follows Ratcliff and Rouder (1998) and has 20 parameters per participant:

- 2 $\times$ 7 drift rates (i.e., one per strength bin) = 14 drift rates $v$
- 2 boundary separations $a$
- 2 non-decision times $t_0$
- 2 drift rate variabilities $s_v$

Also following Ratcliff and Rouder (1998) we fix the start point $z$ to 0.5 (i.e., $a/2$) and fix the start point variability $s_z$ to 0.1. The two boundaries are the two response options "dark" and "light". To estimate the model we diverge from Ratcliff and Rouder 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:

```{r}
# objective function for diffusion with one a. loops over drift to assign drift rates to strength
objective_diffusion_separate <- function(pars, rt, boundary, drift, ...) {
  base_par <- 3  # 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=0.5, 
                 sz=0.1, sv=pars["sv"], #st0=pars["st0"], 
                 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. 

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

```{r}

# function that creates random start values
get_start <- function(base_par) {
  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), 
    z = runif(1, 0.4, 0.6), 
    sz = runif(1, 0, 0.5),
    sv = runif(1, 0, 0.5),
    st0 = runif(1, 0, 0.5),
    d = rnorm(1, 0, 0.05)
  )
  start2 <- sort(rnorm(7), decreasing = FALSE)
  names(start2) <- paste0("v_", 1:7)
  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) {
  
  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, 1:7), 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, 1:7), instruction = data$instruction,
                lower = c(rep(0, length(base_pars)), -Inf,
                          rep(-Inf, length(start_function(base_pars))-length(base_pars))))
  
  fit
}

```

```{r, echo=FALSE}
load("rr98_full-diffusion_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 the `dplyr` `do` function.


```{r, eval = FALSE}

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", "sv"))) %>% ungroup()

```

The following table gives the parameters and the negative summed log-likelihoods obtained from this fit:

```{r}
pars_separate <- as.data.frame(fits_separate %>% group_by(id, instruction) %>% do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>% ungroup())
pars_separate$ll <- (fits_separate %>% group_by(id, instruction) %>% do(ll = .$diffusion[[1]][["objective"]]) %>%  summarize(ll2 = mean(ll[[1]])) %>% as.data.frame())[[1]]

knitr::kable(pars_separate)

```

We can see from these values that there is a large effect of instruction on $a$. However, instruction also has effects on other parameters:

- $t_0$ is consistently larger in the accuracy compared to the speed condition.
- $s_v$ is estimated at 0 for all participants in the speed condition, but 0.46 for one and around 0.8 for the other two participants.
- The differences in $v$ are not consistent across participants, but also do not seem completely negligible.

### Joint Fits

Despite this apparent differences in more than just $a$ we now fit the joint model in which only $a$ is allowed to vary between instruction conditions and all other parameters are shared. The joint (i.e., reduced) model has 11 parameters per participant:

- 7 drift rates (i.e., one per strength bin)
- 2 boundary separations $a$
- 1 non-decision time $t_0$
- 1 drift rate variability $s_v$
- We again fix the start point $z$ to 0.5 (i.e., $a/2$) and fix the start point variability $s_z$ to 0.1.

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. 

```{r}
objective_diffusion_joint <- function(pars, rt, boundary, drift, instruction) {
  base_par <- 4
  densities <- vector("numeric", length(rt))
  as <- c(pars["a_1"], pars["a_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=pars["t0"], z=0.5, 
                   sz=0.1, sv=pars["sv"], #st0=pars["st0"], 
                   v=pars[base_par+i]), 
        error = function(e) 0)  
    }  
  }
  if (any(densities == 0)) return(1e6)
  return(-sum(log(densities)))
}


```


```{r, include=FALSE, eval=FALSE}

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", "sv"))) %>% ungroup()

```


The following table gives the parameters and the negative summed log-likelihoods obtained from this fit:

```{r}
pars_joint <- as.data.frame(fits_joint %>% group_by(id) %>% do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>% ungroup())

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)

```


```{r obtain_fits_not_run, eval = FALSE, include = FALSE}

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", "sv"))) %>% ungroup()

objective_diffusion_separate_alt <- function(pars, rt, boundary, drift, ...) {
  base_par <- 7  # 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=pars["z"], 
                 sz=pars["sz"], sv=pars["sv"], st0=pars["st0"], 
                 d = pars["d"],
                 v=pars[base_par+i]), 
      error = function(e) 0)  
  }
  if (any(densities == 0)) return(1e6)
  return(-sum(log(densities)))
}


fits_separate_alt <- 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_alt, 
                            base_pars = c("a", "t0", "sv", "sz", "z", "st0", "d"))) %>% ungroup()

fits_separate_alt_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_alt, 
                            base_pars = c("a", "t0", "sv", "sz", "z"))) %>% ungroup()
pars_separate <- as.data.frame(fits_separate_alt %>% group_by(id, instruction) %>% do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>% ungroup())

pars_separate_b <- as.data.frame(fits_separate_alt_b %>% group_by(id, instruction) %>% do(as.data.frame(t(.$diffusion[[1]][["par"]]))) %>% ungroup())


all.equal(pars_separate, pars_separate_b, tolerance = 0.001)

pars_separate$ll <- (fits_separate %>% group_by(id, instruction) %>% do(ll = .$diffusion[[1]][["objective"]]) %>%  summarize(ll2 = mean(ll[[1]])) %>% as.data.frame())[[1]]


fits_joint <- rr98 %>% 
  group_by(id) %>% # we loop only across instruction
  do(diffusion = ensure_fit(data = ., start_function = get_start, 
                            objective_function = objective_diffusion_joint, 
                            base_pars = c("a_1", "a_2", "t0", "sv"))) %>% ungroup()


save(fits_separate, fits_joint, fits_separate_alt, file = "rr98_full-diffusion_fits.rda")

# save for fast-dm

tmp_out <- rr98 %>% filter(id == "jf" & instruction == "speed") %>% mutate(resp = response_num-1) %>% select(strength_bin, resp, rt) 
write.table(tmp_out, file = "jf_speed.dat", sep = "\t", row.names = FALSE, col.names = FALSE, quote =FALSE)  

rr98 %>% group_by(id, instruction, strength_bin, response) %>% summarise(n = n()) %>% spread(strength_bin, n)

## single tries 
dput(pars_separate[1,])

pdiffusion(rt = 30, boundary = "upper", a = 0.579802547413743, t0 = 0.195986367994755, sv = 0.897142805898331, sz = 0.596868027475543, z = 0.430283354091839, st0 = 0.144057563742441, d = 0.00190627271452611, v = -4.92625351825116)

# JF speed, v_1
(xp <- pdiffusion(rt = 20, boundary = "lower", a = 0.579802547413743, t0 = 0.195986367994755, sv = 0.897142805898331, sz = 0.596868027475543, z = 0.430283354091839, st0 = 0.144057563742441, d = 0.00190627271452611, v = -4.92625351825116))

qdiffusion(xp*quantiles, "lower", a = 0.579802547413743, t0 = 0.195986367994755, sv = 0.897142805898331, sz = 0.596868027475543, z = 0.430283354091839, st0 = 0.144057563742441, d = 0.00190627271452611, v = -4.92625351825116)

(xp <- pdiffusion(rt = 20, boundary = "lower", a = 0.667287, t0 = 0.254885, sv = 0.448566, sz = 0.416655, z = 0.500000, st0 = 0, d = 0, v = -3.388013))

qdiffusion(xp*quantiles, "lower",  a = 0.667287, t0 = 0.254885, sv = 0.448566, sz = 0.416655, z = 0.500000, st0 = 0, d = 0, v = -3.388013)

(xp <- pdiffusion(rt = 20, boundary = "lower", a = 0.667287, t0 = 0.254885, sv = 0.448566, sz = 0.416655, z = 0.500000, st0 = 0, d = 0, v = 3.005193))

qdiffusion(xp*quantiles, "lower",  a = 0.667287, t0 = 0.254885, sv = 0.448566, sz = 0.416655, z = 0.500000, st0 = 0, d = 0, v = 3.005193)


head(rr98)
rr98 %>% filter(id == "jf" & instruction=="speed") %>% group_by(strength_bin, response) %>% 
  summarise(n = n(), med_rt = median(rt), q1 = quantile(rt, 0.1), q5 = quantile(rt, 0.9))


```

### Compare Model Fits

To compare the fits between the two diffusion version we simply perform a likelihood ratio test between the separate and the joint model. For this we first need to sum the two log-likelihoods of the joint model. Then we simply take their difference and can get the p-value for the likelihood ratio tests with 9 degrees of freedom (i.e., their difference in the number of parameters). As we can see the test is hgihly significant. Even if we consider the large sample size and obtain the ciritcal chi-square value from a compromise power test with a small effect size of $w=0.1$, $df=9$, and $\alpha = \beta$ the critical value for the largest $n$ (for participant `nh`, $n=8532$) is 35.06 which is way below our value of 238 for this participant.

```{r, message=FALSE}

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 = afex::round_ps(pchisq(ll_diff_2, df = 9, lower.tail = FALSE)))
#rr98 %>% group_by(id) %>% summarise(n())
knitr::kable(ll_tab)

```

### Graphical Model Fit

TTo evaluate the fits graphically we first compare the actual response rates for the two responses with the predicted responses rates.
he following plots show the quantiles of the data in grey and the predicted qunatiles from the separate fits in black. As one can see, the models have some problems adequately describing the asymmetric data patterns.

#### Separate Fits

```{r, eval = FALSE, include = FALSE}


# pars_separate_l <- pars_separate %>% gather("strength_bin", "v", v_1:v_7) %>% group_by(id, instruction, strength_bin) %>%
#   mutate(resp_prop = pdiffusion(rt=50, boundary="lower", a=a, v=v, t0=t0, sz = sz, z=z, sv=sv, d=d, st0=st0)) 
# separate_pred_dark <- pars_separate_l %>% do(as.data.frame(t(qdiffusion(quantiles*.$resp_prop, boundary="lower", a=.$a, v=.$v, t0=.$t0, sz = .$sz, z = .$z, sv=.$sv, d=.$d, st0=.$st0)))) %>% 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, sz = .$sz, z = .$z, sv=.$sv, d=.$d, st0=.$st0)))) %>% ungroup() %>% gather("quantiles", "light", V1:V5)


pars_separate_l <- pars_separate %>% gather("strength_bin", "v", v_1:v_7) %>% group_by(id, instruction, strength_bin) %>%
  mutate(resp_prop = pdiffusion(rt=20, boundary="lower", a=a, v=v, t0=t0, sz = 0.1, sv=sv)) 
separate_pred_dark <- pars_separate_l %>% do(as.data.frame(t(qdiffusion(quantiles*.$resp_prop, boundary="lower", a=.$a, v=.$v, t0=.$t0, sz = 0.1, sv=.$sv)))) %>% 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, sz = 0.1, sv=.$sv)))) %>% ungroup() %>% gather("quantiles", "light", V1:V5)
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$strength_bin <- factor(substr(separate_pred$strength_bin, 3,3), levels = as.character(1:7))

pars_joint_l <- pars_joint %>% gather("strength_bin", "v", v_1:v_7) %>% gather("instruction", "a", a_1, a_2)
pars_joint_l$instruction <- factor(pars_joint_l$instruction, levels = c("a_1", "a_2"), labels = c("speed", "accuracy"))
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, sz = 0.1, sv=sv)) 
joint_pred_dark <- pars_joint_l %>% do(as.data.frame(t(qdiffusion(quantiles*.$resp_prop, boundary="lower", a=.$a, v=.$v, t0=.$t0, sz = 0.1, sv=.$sv)))) %>% 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, sz = 0.1, sv=.$sv)))) %>% ungroup() %>% gather("quantiles", "light", V1:V5)
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$strength_bin <- factor(substr(joint_pred$strength_bin, 3,3), levels = as.character(1:7))




```


```{r, fig.height=5, fig.width=7}

p1 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = response == "dark" & instruction == "speed", layout = c(3,1), col = "grey")
p2 <- xyplot(dark ~ strength_bin|id + instruction, separate_pred, group = quantiles, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed")
p2 + as.layer(p1)


```

Response = dark, speed condition:

```{r, fig.height=5, fig.width=7}

p1 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = response == "dark" & instruction == "speed", layout = c(3,1), col = "grey")
p2 <- xyplot(dark ~ strength_bin|id + instruction, separate_pred, group = quantiles, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed")
p2 + as.layer(p1)


```

Response = dark, accuracy condition:

```{r, fig.height=5, fig.width=7}

p1 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = response == "dark" & instruction == "accuracy", layout = c(3,1), col = "grey")
p2 <- xyplot(dark ~ strength_bin|id + instruction, 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, 2.4))))
p2 + as.layer(p1)


```

Response = light, speed condition:

```{r, fig.height=5, fig.width=7}

p1 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = response == "light" & instruction == "speed", layout = c(3,1), col = "grey")
p2 <- xyplot(light ~ strength_bin|id + instruction, separate_pred, group = quantiles, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = instruction == "speed")
p2 + as.layer(p1)


```

Response = light, accuracy condition:

```{r, fig.height=5, fig.width=7}

p1 <- xyplot(rt ~ strength_bin|id, agg2_rr98_response, group = quantile, type = "b", auto.key = list(lines = TRUE), ylab = "RT (in seconds)", subset = response == "light" & instruction == "accuracy", layout = c(3,1), col = "grey", ylim = c(0.2, 2))
p2 <- xyplot(light ~ strength_bin|id + instruction, 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.2, 3))))
p2 + as.layer(p1)


```

## 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

