Advice on distribution for modelling reaction times? (and setting priors)


I am trying to model reaction times with a relatively simply formula:

rt ~ 0 + condition + n:condition + (1 | p_id)

where condition is a three level factor, and n is an integer (0 - 36).

A standard lmm works well enough, does a good of fitting the means, etc But really, we want a model that won’t predict values <0, and can match the skew of the distribution.

Using lognormal() helps quite a lot, and I’m almost tempted to call it a day and get on with submitting my paper! But, when the posterior predictions are plotted, you can see that model is predicting too may short rts . We next tried using a shifted_lognormal() to account for this, but we didn’t get very far.

I’m currently trying to work with family = inverse_gaussian(), and this is working more or less. Although it it often fails to find good starting values.

We’re having quite a bit of trouble in specifying sensible priors. We’d really really like to have priors for each of the three distributions that result in broadly similar prior predictions (we know that most our of data is going to be between 0.2 and 2 seconds!). Whenever I try to run the brms model with sample_prior = "only", we get lots of initalising warnings etc.

I’ve even tried running the model with the default priors, then using those posteriors to inform a new set of priors to test, but I can’t get this working either!

Does anybody have any experience working with any of these distributions? It would be great to get some advice! I’m happy to switch back to shifted_lognormal, stick with inverse_gaussian, or move to something else.

Thanks in advance,


Have you taken a look at a shifted Weibull distribution? Jeff Rouder did some work back in the early oughts with it and I find it a really nice descriptive distribution for RTs.

For reading times, I always found that the shifted lognormal worked well. Other people use the reciprocal normal, 1/rt ~ normal, where you are actually modelling speed rather than rt.

But this is a nice paper about the Weibull :
And if you want to explore other possible distributions, check

BTW, I added the tag cognitive-science to your post. And after you submit a preprint, you can add your paper to this website I’m maintaining:

1 Like

Do you have any advice on how to set up the shifted lognormal? If you have any example code you’re willing to share, that would be great!

Or, I will sort out a small example of what we’re trying to do and post here :)

I have code, but not with brms.

This is really old code with Stan:

There must be better ways to fit it now with Stan.

Brms hast the shifted_lognormal family() for that purpose in case you have overlooked it.


hi Paul,

first of all, thanks as always for the amazing work you are doing!

We were trying to get shifted_lognormal to work, but we weren’t getting anywhere very quickly.

Am I right in thinking that I probably want something like:

bf( y ~ mx+c + (1|pid), ndt ~ 1)

at least to start with. Once we get that working, we can try something more sophisticated for ndt!

Thanks, that’s still helpful… at the very least, it suggests some good starting points for setting priors!


still not much luck here.

If I :

remove the random intercepts
use the default priors
only use 10% of the data (to speed things up)

the model appears to run fine, although it does throw a load of initialization errors:

SAMPLING FOR MODEL ‘f4f04fdddfc097a21f0b3152cccaf0e9’ NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: lognormal_lpdf: Random variable[1] is -0.228194, but must be >= 0! (in ‘model2c9438c04639_f4f04fdddfc097a21f0b3152cccaf0e9’ at line 35)

But the resulting posterior fit looks really good, Rhats are good, etc

Should I try and fix these initialization problems before trying to add the random effects and my own priors back in?

m <- brm(
rt ~ 0 + d_feature + d_feature:log(N_T+1) ,#+ (1|p_id),
ndt ~ 1),
family = shifted_lognormal(),
data = d,
chains = 2,
iter = 3000,
control =list(adapt_delta = 0.9)