-
-
Save jvcasillas/628232a0c759e97b9911ab3eeff26241 to your computer and use it in GitHub Desktop.
library(tidyverse) | |
# DDM simulation function | |
sim_ddm <- function(drift_rate, boundary_separation, bias, ndt, n_sims, seed = NULL) { | |
set.seed(seed) | |
sim_df = NULL | |
for (sim in 1:n_sims){ | |
step = 1 | |
value = bias | |
n = 1 | |
while (abs(value[n]) < boundary_separation) { | |
n = n + 1 | |
value[n] = (value[n - 1] + rnorm(1, 0, abs(drift_rate))) | |
step[n] = n | |
} | |
dd_df <- tibble(sim_n = sim, step, value) %>% | |
mutate( | |
sim_n = as.factor(sim_n), | |
value = case_when( | |
value > boundary_separation ~ boundary_separation, | |
value < -boundary_separation ~ -boundary_separation, | |
TRUE ~ .$value) | |
) | |
sim_df = rbind(sim_df, dd_df) | |
} | |
return(sim_df) | |
} | |
# Get α, β, δ, and τ from DDM | |
# Plot 20 random simulations after fitting model | |
sim_ddm( | |
drift_rate = 0.8, | |
boundary_separation = 1, | |
bias = 0, | |
ndt = 0.2, | |
n_sims = 20, | |
seed = 12345) %>% | |
ggplot(., aes(x = step, y = value, color = sim_n)) + | |
geom_line(show.legend = F) | |
# You can increase the number of sims and plot the average trajectories too |
Thank you! Have you tired to use brms to run your ddm?
Yep. I looked into a few other options (there is an interesting one in python but I couldnt get it to work: http://ski.clps.brown.edu/hddm_docs/index.html), but brms is pretty much home for me now. How do you usually fit yours?
I haven't done much DDM. I am planning to collect some data online next semester that I want to look at with DDM so I am exploring options.
@jvcasillas Can you share the code you used to make the figure as well? I think it would be a good way to depict predictions one has about the different parameters.
Ah the plot you are referring to is just a schematic to explain the model (it's in the method section of the paper Im working on). I use it when I explain what each parameter represents. Here is the code:
# DDM explanation -------------------------------------------------------------
# Generate skewed RT-like data
my_skewed_x <- fGarch::rsnorm(100, mean = 2.5, sd = 1.5, xi = 3.5)
# Correct response dist
inset_top <- tibble(x = my_skewed_x) %>%
ggplot(., aes(x = x)) +
geom_density(fill = viridis::viridis_pal(option = "B", begin = 0.3)(1)) +
theme_void()
# Incorrect response dist
inset_bottom <- tibble(x = my_skewed_x) %>%
ggplot(., aes(x = x)) +
geom_density(fill = viridis::viridis_pal(option = "B", begin = 0.75)(1)) +
scale_y_reverse() +
theme_void()
threshold_labs <- tibble(
step = 80, value = c(1.4, -1.4),
labs = c("Correct response", "Incorrect response")
)
# Main plot
main_plot <- sim_ddm(q_type = "yn", eq = "low", lt = "low", seed = 20180119,
drift_rate = 0.11, boundary_separation = 1.25, bias = 0, ndt = 0, n_sims = 2) %>%
ggplot(., aes(x = step, y = value, color = sim_n)) +
geom_vline(xintercept = 1) +
geom_hline(yintercept = c(-1.25, 1.25), lty = 3) +
geom_line(show.legend = F, size = 1.25) +
geom_segment(inherit.aes = F, aes(x = 1, xend = 98, y = 0, yend = 0),
lty = 2, size = 0.25) +
coord_cartesian(ylim = c(-2, 2), xlim = c(-5, 100)) +
scale_x_continuous(breaks = c(1, seq(10, 100, 10)),
labels = c("0", seq(10, 100, 10))) +
scale_y_continuous(breaks = NULL) +
labs(y = expression(alpha), x = "Time step") +
annotate("text", x = -4.5, y = 0.15, size = 5, label = expression(tau)) +
geom_segment(inherit.aes = F, aes(x = -5, y = 0, xend = 0, yend = 0),
arrow = arrow(length = unit(0.2, "cm")), size = 0.5) +
geom_segment(inherit.aes = F, aes(x = -5, y = 0, xend = -9, yend = 0),
arrow = arrow(length = unit(0.2, "cm"))) +
annotate("text", x = 101, y = 0, size = 5, label = expression(beta)) +
geom_segment(inherit.aes = F, aes(x = 101, y = 0.2, xend = 101, yend = 0.5),
arrow = arrow(length = unit(0.2, "cm")), size = 0.5) +
geom_segment(inherit.aes = F, aes(x = 101, y = -0.2, xend = 101, yend = -0.5),
arrow = arrow(length = unit(0.2, "cm"))) +
annotate("text", x = 32, y = 0.25, size = 5, label = expression(delta)) +
geom_segment(inherit.aes = F, aes(x = 30.5, y = 0.35, xend = 25, yend = 0.7),
arrow = arrow(length = unit(0.2, "cm"))) +
geom_text(inherit.aes = F, data = threshold_labs, hjust = 0,
aes(x = step, y = value, label = labs), family = "Times") +
scale_color_viridis_d(option = "B", begin = 0.3, end = 0.75) +
theme_minimal(base_family = 'Times', base_size = 13) +
theme(panel.grid.major = element_line(size = 0.2),
panel.grid.minor = element_blank(),
axis.title.y = element_text(size = rel(.9), hjust = 0.5))
# Combine plots
ddm_explanation <- main_plot + annotation_custom(
grob = ggplotGrob(inset_top),
ymin = 1.25, ymax = 2, xmin = -2.5, xmax = 75) +
annotation_custom(
grob = ggplotGrob(inset_bottom),
ymin = -1.25, ymax = -2, xmin = -2.5, xmax = 75)
# -----------------------------------------------------------------------------
This is how Im using it:
@jgeller112 👍