A mixture of shifted lognormal and uniform distribution suitable for modelling reaction times.
Source:R/shifted_lognormal_uniform.R
shifted_lognormal_uniform.RdA contaminated response time distribution. The mixture can be described as
$$y_i =
\begin{cases}
u_i & \mathrm{if} \quad z_i = 0 \\
p_i s_i + r_i & \mathrm{if} \quad z_i = 1
\end{cases}
\\
u_i \sim Uniform(0, \alpha) \\
\log(r_i) \sim Normal(\mu_i, \sigma) \\
P(z_i = 0) = \theta \\
0 < p_i < 1
$$
Here \(\mu, \sigma, \theta\) (mix) and \(p\) (shiftprop)
are estimated, whereas \(s_i\) (max_shift)
and \(\alpha\) (max_uniform) are given as data via vreal().
Usage
shifted_lognormal_uniform(
link = "identity",
link_sigma = "log",
link_mix = "logit",
link_shiftprop = "logit"
)Source
Idea by Nathaniel Haines (https://twitter.com/Nate__Haines), code by Martin Modrák. Some background, discussion and examples at http://www.martinmodrak.cz/2021/04/01/using-brms-to-model-reaction-times-contaminated-with-errors/
Details
Note that you cannot build this distribution with the built-in support
for mixtures in brms,
because the uniform component is effectively a zero-parameter distribution
which cannot be expressed in brms.
Examples
library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.22.13). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#>
#> Attaching package: ‘brms’
#> The following object is masked from ‘package:stats’:
#>
#> ar
set.seed(31546522)
# Bounds of the data
max_shift <- 0.3
shift <- runif(1) * max_shift
max_uniform <- 10
mix <- 0.1
# Generate parameters
N <- 100
Intercept <- 0.3
beta <- 0.5
X <- rnorm(N)
mu <- rep(Intercept, N) + beta * X
sigma <- 0.5
rt <- rshifted_lognormal_uniform(N, meanlog = mu, sdlog = sigma, mix = mix,
shift = shift, max_uniform = max_uniform)
dd <- data.frame(rt = rt, x = X,
max_shift = max_shift, max_uniform = max_uniform)
fam <- shifted_lognormal_uniform()
fit_mix <- brm(rt | vreal(max_shift, max_uniform) ~ x, data = dd, family = fam,
stanvars = fam$stanvars,
prior = c(prior(beta(1, 5), class = "mix")))
#> Compiling Stan program...
#> Error in .fun(model_code = .x1): Boost not found; call install.packages('BH')
plot(fit_mix)
#> Error: object 'fit_mix' not found