Q-diffusion IRT model interpretation in brms

I have been working with a simulated Q-diffusion model in brms, like the one described in van der Maas et al. (2011; APA PsycNet) and in Bürkner (2021; doi: 10.18637/jss.v100.i05). I’ve been referencing those resources as well as Henrik Singmann’s blog posts on diffusion models in brms (http://singmann.org/wiener-model-analysis-with-brms-part-i/). I have a few questions I’m hoping I could get some input on. Apologies if some of these are obvious. Code to simulate data and run the Q-diffusion model is below.

  1. I understand the interpretation of each q-diffusion parameter, but I’m having difficulty understanding the scales of drift rate, boundary separation, and non-decision time. E.g., in Figure 16 of the Bürkner IRT paper, are the y-axes for each graph log units since the link for each parameter in the model is log? If these y-axes are log units and they are transformed back into their raw units, what are those raw units for drift, bs, and ndt? When I do pp_check with my model, I’m also unclear what scale that distribution is on.

  2. I’d like to obtain the item characteristic curves (ICCs) for each item, which is pretty straight forward for regular 1PL or 2PL IRT models in brms (Make ICC plots for your brms IRT models | A. Solomon Kurz). van der Maas et al. (2011) show ICCs for the Q-diffusion in Figure 3 where ability is strictly positive. Following Henrik Sigmann’s blog, I’ve obtained the predicted response times for upper- and lower-boundaries and calculated the probability of a correct response as a response time > 0 (a correct/upper-boundary response), but I’m not sure if I’m on the right track with this. Alternatively, would I use the same approach that Solomon Kurz uses in the linked blog above (i.e., getting posterior samples for each item, creating a theta variable, calculating probability from the model intercept, item, and theta using inv_logit_scaled)?

  3. I would like to be able to estimate the Q-diffusion model in brms with open-ended responses (like van der Maas et al. do in WinBUGS – code for that from one of the author’s websites is attached). Would this be possible in brms, or would this require more custom specification in Stan?
    Syntax Q-diffusion model.pdf (101.1 KB)
    ).

Any input on these questions is appreciated!

Simulation and model fitting code

library(diffIRT)
library(brms)
library(tidyverse)
library(job)

#simulate data using diffIRT pacakge
#create dataset of 50 individuals and 30 items
simdif <- simdiff(N=50,nit=30, model = "Q")

#reaction times
d_rt <- simdif[["rt"]] %>% as_tibble() %>% rowid_to_column("id")

#items
d_ir <- simdif[["x"]] %>% as_tibble() %>% rowid_to_column("id")  

#rename columns 
colnames(d_rt) <- c("id","1":"30")
colnames(d_ir) <- c("id","1":"30")

#convert to long format
d_rtl <- d_rt %>% gather(item, time, 2:31)
d_irl <- d_ir %>% gather(item, resp, 2:31)

#merge reaction time and item datasets
d <- merge(d_rtl, d_irl, by = c("id", "item")) 
d$item <- as.numeric(d$item)

#create model--major won't run if I estimate  random effects for ndt, so just average ndt is estimated
mod <- bf(
  time|dec(resp) ~ (1|p|id) + (1|i|item),
  bs ~ (1|p|id) + (1|i|item),
  ndt ~ 1,
  bias = 0.5)

#run q-diffusion model
qdiff2  <- brm(mod, 
          data = d,
          family = brmsfamily("wiener", "log", link_bs = "log", link_ndt = "log"),
          chains = 3, 
          cores = 8, 
          backend = "cmdstanr", 
          control =
            list(c(adapt_delta = .99, max_treedepth = 15)))

#get posterior predictions
pred <- predict(qdiff2, 
                       summary = FALSE, 
                       negative_rt = TRUE, 
                       ndraws = 10)
postpreds <- as_tibble(cbind(d, as_tibble(t(pred)))) %>% 
  gather(postsamp, rt, V1:V10) %>% 
  group_by(id, item) %>% 
  mutate(acc = if_else(rt > 0, 1, 0)) %>% ungroup()

probs <- postpreds %>% group_by(id, item) %>% summarise_at(c("acc", "rt"), mean)```


* Operating System: MacOS Big Sur, 2.4 GHz 8-Core Intel Core i9
* brms Version: 2.17.0