For a simple model

y_ij | u_i ~ N(u_i, s^2)

we usually assume u_i ~ N(0, tau^2) when obtaining the posterior distribution for u_i. However, if a bimodal (instead of unimodal) distribution is more appropriate for u_i, is there any way to handle the situation?

u_i ~ N(0, tau^2)

u_i

Ideally, you would condition on whatever is the source of the bimodality. Failing that, a mixture of two normals. See the log_mix function.

log_mix

Thanks, Ben! Is there a specification available with stan_lmer() in rstanarm? Or I have to code it directly using Stan?

stan_lmer()

rstanarm

Stan

Nothing like that is in the rstanarm package