Diffusion model initialization issue using brms

Modeling
, ,
1

Dear all,
I am trying to estimate the drift rate in a task-switching paradigm.

In my experiment, which is a two-by-two design, I am exploring two stimulus types(“congruent” vs “incongruent”) and two task types (“repeat trial” vs “switch trial”). And I used the Wiener diffusion model together with brms package to fit the parameters. Generally I followed the tutorial by Dr. Singmann here.

Unfortunately, I got the error message upon the initialization phase.
I ran one chain, and it says,

Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: Exception: wiener_lpdf: Nondecision time is -0.000770665, but must be positive finite! (in ‘string’, line 29, column 7 to column 55) (in ‘string’, line 29, column 7 to column 55)
Chain 1:
Chain 1: Initialization between (-2, 2) failed after 100 attempts.
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done”

To address this, I tried

  1. filter the fast RT trials, set a threshold of 0.1s and remove those trials below.
  2. set a small non-decision time, I tried 0.001, 0.05 and a uniform(0.001, 0.05) as random initial value for ndt.
  • Operating system: x86_64-w64-mingw32
  • R version: 4.3.1 (2023-06-16 ucrt)
  • brms version: 2.20.4

Below, I have included the relevant portion of my code for your reference:

library(brms)

formula <- bf(
  rt | dec(acc) ~ 0 + stimtype:tasktype + (0 + stimtype:tasktype | p | subj_num),
  bs ~ 0 + stimtype:tasktype  + (0 + stimtype:tasktype | p | subj_num),
  ndt ~ 0 + stimtype:tasktype + (0 + stimtype:tasktype | p | subj_num),
  bias = 0.5
)

prior <- c(
  prior("cauchy(0, 5)", class = "b"),
  set_prior("gamma(3, .5)", class = "b", dpar = "bs",lb=0),
  set_prior("gamma(1, .5)" , class = "b", dpar = "ndt",lb=0)
)

make_stancode(formula,
              family = wiener(link_bs = "identity",
                              link_ndt = "identity"),
              data = data_pre,
              prior = prior)

pre_dat <- make_standata(
  formula,
  family = wiener(
    link_bs = "identity",
    link_ndt = "identity"
  ),
  data = data_pre,
  prior = prior
)
initfun <- function() {
  list(
    b = rnorm(pre_dat$K),
    b_bs = runif(pre_dat$K_bs, 1, 2),
    b_ndt = runif(pre_dat$K_ndt, 0.001, 0.05),
    # b_ndt = 0.05,
    sd_1 = runif(pre_dat$M_1, 0.5, 1),
    z_1 = matrix(rnorm(pre_dat$M_1*pre_dat$N_1, 0, 0.01),
                 pre_dat$M_1, pre_dat$N_1),
    L_1 = diag(pre_dat$M_1)
  )
}


fit_wiener_pre <- brm(
  formula,
  data = data_pre,
  family = wiener(
    link_bs = "identity",
    link_ndt = "identity"
  ),
  prior = prior,
  init = initfun,
  iter = 1000, warmup = 500,
  chains = 1, cores = 4,
  control = list(max_treedepth = 15,adapt_delta = 0.9)
)

Could someone help me out what went wrong with the initialization? Or how could I debug the model?

Help is much appreciated!

Best, Shuai

2

With these link functions bs and ndt parameters are not constrained to be positive. Based on the exception, specifically, there is a problem with ndt. If you don’t like the default link functions

  link_bs = "log",
  link_ndt = "log",
)

Maybe “softplus” would be better than “identity”?

1 Like
3

Hey, avehtari, thanks for your fast response! I’d like to clarify something regarding the ‘ndt’ parameter in my model. Even though I’ve set a Gamma prior( which inherently constrains the parameter to be positive) and initialized it with values drawn from a uniform dist between 0.001 and 0.5, will there be any circumstances the ‘ndt’ parameter might not remain positive? Would you please explain a bit why this might be the case?

And btw, I got the error message stating, “Chain 1: Initialization between (-2, 2) failed after 100 attempts.” Does this suggest that I might not succeed passing initifun to brm()? The message seems to imply that Stan is defaulting to random initialization values between -2 and 2, instead of using the initfun I passed to it. Could this indicate a problem with the integrating of my initfun into the brm function?

Thank you for your insights on these matters.

4

As data_pre was no defined in your code I can’t check, but maybe the problem is the following (quoting from the brm doc)

If specifying initial values using a list or a function then currently the parameter names must correspond to the names used in the generated Stan code (not the names used in R).

5

Just looking quickly at your code, could it be due to the group-level effects in the model? That is, you’ve set a prior for ndt that covers the overall mean of each condition, but you haven’t specified any priors for the distribution of subjects around that mean. If you use get_prior() this will show you all the possible priors you can set, and you’ll probably see a bunch of ndt-related group-level priors which are not covered by those you’ve set, and so are using whatever the default is. This is probably leading to the negative ndt values you’re observing.

If that’s the case, then you could try to select priors for the group-level effects which ensure the ndt values never go below 0. However, it’s probably easier to just use a link function on ndt which prevents it from going below 0, as avehtari suggested.

3 Likes
6

Thank you, @ps2 and @avehtari .

Now I figured out what the problem is.
I added individual level random effects on top of the fixed effects in the formula:

formula <- bf(
  rt | dec(acc) ~ 0 + stimtype:tasktype + (0 + stimtype:tasktype | p | subj_num),
  bs ~ 0 + stimtype:tasktype  + (0 + stimtype:tasktype | p | subj_num),
  ndt ~ 0 + stimtype:tasktype + (0 + stimtype:tasktype | p | subj_num),
  bias = 0.5
)

Therefore, even I intentionally choose the prior and starting values to be positive, those random effects could drag the values below 0, and caused the ndt error.

Thank you again, guys, for your precious time and valuable help!
Wish you a nice day!

1 Like