Using brms to model reaction times contaminated with errors

· 2021/04/01 · 23 minute read

Nathaniel Haines made a neat tweet showing off his model of reaction times that handles possible contamination with both implausibly short reaction times (e.g., if people make an anticipatory response that is not actually based on processing the stimulus of interest) or implausibly large reaction times (e.g., if their attention drifts away from the task, but they snap back to it after having “zoned out” for a few seconds). Response times that arise from such processes are not actually what we aim to measure in most cognitive tasks — we are instead interested in how people process and respond to a particular stimulus. Therefore, by explicitly modeling these “contamination” response times, we can get better estimates of the decision-making parameters that we actually care about. Such a model often makes more sense than just throwing away a part of the data.

Several people asked, if you can do that in brms. This started a vortex of productive procrastrination on my side - it sure should be easy to do this, right? And while Nathaniel didn’t have a brms code ready, I assure you that, yes, it is possible in brms, it is not completely straightforward, but I’ll show you the path code.

Nathaniel was kind enough to provide a bit of feedback on the post (I have no experience with reaction-time data or cogsci in general), but I should repeat that the clarity of the idea is his while all errors are mine. The overall idea of using mixtures of “real” and “contaminating” distributions to model responses is however not new - see e.g. Ratcliff & Tuerlinckx 2002.

In this model we will take a shifted lognormal representing the actual decision process and a uniform distribution modelling the contamination. For this to make sense, we need to have some upper limit on the possible contamination times, representing the maximum times we could have observed. In most cases, the limit should be larger than the maximum time we observed, although this is not strictly mathematically necessary. We then assume that each trial has a small probability of being contaminated.

Here is how generating a single data point for such a model could look in R code:

shift <- 0.1 # Shortest reaction time possible if not contaminated
mu <- log(0.5)
sigma <- 0.6
mix <- 0.06 # Probability of contamination
upper <- 5 # Maximum time of contamination

if(runif(1) < mix) {
  # Contaminated
  y <- runif(1, 0, upper)
} else {
  # Non-contaminated
  y <- shift + rlnorm(1, mu, sigma)
}

The same could be expressed in math as:

\[ y_i = \begin{cases} u_i & \mathrm{if} \quad z_i = 0 \\ 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 \]

Where \(\theta\) corresponds to mix, \(\alpha\) to upper and \(s_i\) to shift.

Technically, the non-contaminated signal is allowed to take values larger than upper. In practice we would however usually want upper to be large enough that larger values do not really occur.

There is one important detail in how brms does handle the shifted lognormal: brms does treat shift as unknown and estimates it, but does not allow the shift parameter to be larger than any actually observed y. We will therefore mimic this behaviour, but since we also have the contamination process, shift can in principle be larger than some y. This can potentially introduce problems for the sampler as the posterior density is not smooth when shift crosses some of the observed y values (the lognormal component is added/removed, resulting in a sharp change).

It however turns out that if shift crossing some y is rare enough, the sampling works just fine. To ensure this rarity we introduce max_shift as the upper bound for shift. In most cases, this will be the same for the whole dataset. Instead of shift, the model would then work with shiftprop = shift / max_shift - a value between 0 and 1 that is easier to work with mathematically.

Of the model parameters, we take max_shift and upper as known (but possibly differ between observations) while mu, sigma, mix and shiftprop are to be estimated and can depend on predictors. However, shiftprop is a bit complicated here and the model will make most sense if observations that have different max_shift are also allowed to have different shiftprop by putting a suitable predictor on shiftprop. Different shiftprop with the same max_shift is however definitely not an issue. So while you need to be careful with varying max_shift, varying shiftprop is OK, just note the implied logit scale. For a review on why varying shift might be important see e.g. Dully, McGovern & O’Connell 2018.

For some use cases, one could also want to set the lower bound of the contamination distribution. To keep things simple we don’t do that here, but basically the same result can then be achieved by adding/subtracting a suitable number to the response (y) and bounds (max_shift, upper)

Some experimental designs also involve a limit on the maximum time the response could have taken. In such contexts, it might make sense to treat the values as right-censored. brms supports censoring for most families, so we want our implementation to be compatible with it.

Our goal is that at the end of the post we will be able to write models like

brm(bf(y | vreal(max_shift, upper) + cens(censoring) ~ 1 + condition + (1 | subject_id),
       sigma ~ condition,
       mix ~ (1 | subject_id),
       family = RTmixture), ...)

And let brms handle all the rest. The final result, packaged in a single file you can just load into your project is at https://github.com/martinmodrak/blog/blob/master/content/post/RTmixture.R However, be advised that the code was only quite shallowly tested, so look both ways before crossing and test if you can recover parameters from simulated data before trusting me too much.

You may know that, brms has good support for mixtures, so why not just write family = mixture(uniform, shifted_lognormal)? It turns out brms has as one of its core assumptions that every family has at least one parameter to be estimated - our uniform distribution for the contamination parameter however does not have that and thus cannot be used with brms directly. So instead we’ll have to implement a full blown custom family.

The necessary background for implementing custom families in brms can be found in the vignette on custom distributions. Here, we will explain only the more weird stuff.

Setting up

Let’s set up and get our hands dirty.

library(cmdstanr)
library(brms)
library(tidyverse)
library(knitr)
library(patchwork)
library(bayesplot)

source("RTmixture.R")

ggplot2::theme_set(cowplot::theme_cowplot())
options(mc.cores = parallel::detectCores(), brms.backend = "cmdstanr")

cache_dir <- "_RTmixture_cache"
if(!dir.exists(cache_dir)) {
  dir.create(cache_dir)
}

First, we’ll generate some fake data to test the model against. Below is just a more concise and optimized version of the random generation scheme I showed earlier.

rRTmixture <- function(n, meanlog, sdlog, mix, shift, upper) { 
  ifelse(runif(n) < mix,  
         runif(n, 0, upper),  
         shift + rlnorm(n, meanlog = meanlog, sdlog = sdlog)) 
} 
 

Then let us generate some data

set.seed(31546522)
# Bounds of the data
max_shift <- 0.3
shift <- runif(1) * max_shift
upper <- 10
mix <- 0.1

N <- 100
Intercept <- 0.3
beta <- 0.5
X <- rnorm(N)
mu <- rep(Intercept, N) + beta * X
sigma <- 0.5

rt <- rRTmixture(N, meanlog = mu, sdlog = sigma, mix = mix, shift = shift, upper = upper)

dd <- data.frame(rt = rt, x = X, max_shift = max_shift, upper = upper)

Looking nice!

ggplot(dd, aes(x = rt)) + geom_density()

Core of the family

Now we need the Stan implementation of the family. That is probably the most technical part. Stan user’s guide has some background on mixture models in Stan. We’ll note that times before shift can only come from the uniform component and times after upper can only come from the lognormal component. For others we mix both a lognormal and the uniform via log_mix.

With the Stan code ready, we then define the parameters of the distribution in a way that brms understands.

stan_funs_base <- stanvar(block = "functions", scode = " 
  real RTmixture_lpdf(real y, real mu, real sigma, real mix,  
                      real shiftprop, real max_shift, real upper) { 
    real shift = shiftprop * max_shift; 
    if(y <= shift) { 
      // Could only be created by the contamination 
      return log(mix) + uniform_lpdf(y | 0, upper); 
    } else if(y >= upper) { 
      // Could only come from the lognormal 
      return log1m(mix) + lognormal_lpdf(y - shift | mu, sigma); 
    } else { 
      // Actually mixing 
      real lognormal_llh = lognormal_lpdf(y - shift | mu, sigma); 
      real uniform_llh = uniform_lpdf(y | 0, upper); 
      return log_mix(mix, uniform_llh, lognormal_llh); 
    } 
  } 
 
") 
 
 
 
RTmixture <- custom_family( 
  "RTmixture",  
  dpars = c("mu", "sigma", "mix", "shiftprop"), # Those will be estimated 
  links = c("identity", "log", "logit", "logit"), 
  type = "real", 
  lb = c(NA, 0, 0, 0), # bounds for the parameters  
  ub = c(NA, NA, 1, 1), 
  vars = c("vreal1[n]", "vreal2[n]") # Data for max_shift and upper (known) 
) 

And we are ready to fit! We will put a weakly informative beta(1,5) prior on the proportion of contamination - this means we a prior believe that there is a 95% chance that the contamination is lower than qbeta(0.95, 1, 5) = 0.4507197. One could definitely be justified in tightening this prior even further toward zero for many tasks. vreal is just brms’s way of annotating arbitrary additional data for the distribution. We need to pass both the family and the associated stanvars.

fit_mix <- brm(rt | vreal(max_shift, upper) ~ x, data = dd, family = RTmixture, 
               stanvars = stan_funs_base, 
               refresh = 0,
               file = paste0(cache_dir, "/mix"), file_refit = "on_change",
               prior = c(prior(beta(1, 5), class = "mix")))
fit_mix
##  Family: RTmixture 
##   Links: mu = identity; sigma = identity; mix = identity; shiftprop = identity 
## Formula: rt | vreal(max_shift, upper) ~ x 
##    Data: dd (Number of observations: 100) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.33      0.07     0.20     0.48 1.00     2096     2242
## x             0.44      0.05     0.33     0.54 1.00     2368     2460
## 
## Family Specific Parameters: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma         0.45      0.05     0.36     0.56 1.00     2379     1976
## mix           0.13      0.05     0.05     0.24 1.00     3158     2589
## shiftprop     0.82      0.18     0.32     1.00 1.00     1639     1623
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

We note that we have quite good recovery of the effect of x (simulated as 0.5) and of sigma (which was 0.5), but 100 observations are not enough to constrain the mix parameter really well (simulated as 0.1).

For comparison, we also fit the default shifted lognormal as implemented in brms.

fit_base <- brm(rt ~ x, data = dd, family = shifted_lognormal, refresh = 0,
                file = paste0(cache_dir, "/base"), file_refit = "on_change")

fit_base
##  Family: shifted_lognormal 
##   Links: mu = identity; sigma = identity; ndt = identity 
## Formula: rt ~ x 
##    Data: dd (Number of observations: 100) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.34      0.09     0.16     0.53 1.00     1959     1949
## x             0.47      0.07     0.33     0.61 1.00     2712     2495
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.69      0.07     0.56     0.82 1.00     1883     1802
## ndt       0.34      0.07     0.17     0.44 1.00     1633     1420
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

We see that the inferences for sigma are a bit biased but this is not necessarily only due to the mixture, another potentially biasing is the different handling of the shift.

Censoring + constant shift

To support censoring in brms the family has to come with log CDF (cumulative distribution function) and log CCDF (complementary CDF) implementations in Stan, which we provide below. Those match the _lpdf pretty closely.

 
stan_funs <- stan_funs_base + stanvar(block = "functions", scode = " 
  real RTmixture_lcdf(real y, real mu, real sigma, real mix,  
                      real shiftprop, real max_shift, real upper) { 
    real shift = shiftprop * max_shift; 
    if(y <= shift) { 
      return log(mix) + uniform_lcdf(y | 0, upper); 
    } else if(y >= upper) { 
      // The whole uniform part is below, so the mixture part is log(1) = 0 
      return log_mix(mix, 0, lognormal_lcdf(y - shift | mu, sigma)); 
    } else { 
      real lognormal_llh = lognormal_lcdf(y - shift | mu, sigma); 
      real uniform_llh = uniform_lcdf(y | 0, upper); 
      return log_mix(mix, uniform_llh, lognormal_llh); 
    } 
  } 
   
  real RTmixture_lccdf(real y, real mu, real sigma, real mix,  
                      real shiftprop, real max_shift, real upper) { 
 
    real shift = shiftprop * max_shift; 
    if(y <= shift) { 
      // The whole lognormal part is above, so the mixture part is log(1) = 0 
      return log_mix(mix, uniform_lccdf(y | 0, upper), 0); 
    } else if(y >= upper) { 
      return log1m(mix) + lognormal_lccdf(y - shift | mu, sigma); 
    } else { 
      real lognormal_llh = lognormal_lccdf(y - shift | mu, sigma); 
      real uniform_llh = uniform_lccdf(y | 0, upper); 
      return log_mix(mix, uniform_llh, lognormal_llh); 
    } 
 
  } 
") 
 

To test if this work, we’ll do quite aggressive censoring and treat anything larger than 1.5 as censored. In most cases it makes sense to have upper be the same as the censoring bound, so we’ll do that

set.seed(25462255)
shift <- 0.15
cens_bound <- upper <- 1.5
mix <- 0.08

N <- 110
Intercept <- 0.5
beta <- -0.3
X <- rnorm(N)
mu <- rep(Intercept, N) + beta * X
sigma <- 0.4

rt <- rRTmixture(N, meanlog = mu, sdlog = sigma, 
                 mix = mix, shift = shift, upper = upper)
censored <- rt > cens_bound
rt[censored] <- cens_bound

dd_cens <- data.frame(rt = rt, 
                      censored = if_else(censored, "right", "none"),  
                      x = X, max_shift = shift, upper = upper)

Finally, this model starts to be problematic if we try to estimate shift (well, actually shiftprop) as well. An easy way to to make shift always equal to max_shift is to set a constant prior on shiftprop, as we do below.

fit_mix_cens <- brm(rt | vreal(max_shift, upper) + cens(censored) ~ x, 
                    data = dd_cens, 
                    family = RTmixture, 
                    stanvars = stan_funs, 
                    refresh = 0,
                    file = paste0(cache_dir, "/mix_cens"), 
                    file_refit = "on_change",
                    prior = c(prior(beta(1, 5), class = "mix"),
                              prior(constant(1), class = "shiftprop")))
fit_mix_cens
##  Family: RTmixture 
##   Links: mu = identity; sigma = identity; mix = identity; shiftprop = identity 
## Formula: rt | vreal(max_shift, upper) + cens(censored) ~ x 
##    Data: dd_cens (Number of observations: 110) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.46      0.08     0.34     0.65 1.00     2022     1674
## x            -0.23      0.07    -0.38    -0.12 1.00     2093     1606
## 
## Family Specific Parameters: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma         0.36      0.07     0.26     0.53 1.00     2019     2565
## mix           0.15      0.05     0.07     0.25 1.00     2942     2595
## shiftprop     1.00      0.00     1.00     1.00 1.00     4000     4000
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

It works and the inferences are reasonably close to what we simulated with. A more thorough evaluation would require simulation-based calibration, which would be nice, but would require a bit more energy than I have now. But it seems that at least the models are not completely wrong.

If you want to model varying shift but having issues fitting, it might make sense to adjust max_shift on a per-group basis to have max_shift larger only than a small proportion of observations in this group. As noted above, if you set different max_shift per for example subject_id, you should also have shiftprop ~ subject_id or the model might not make sense.

Making predictions

We successfully fitted a few models, but there are some tweaks we need to do to make full use of the family. We might for example want to make predictions - e.g. to make posterior predictive checks - so we also need to implement prediction code. You’ll notice that we are just extracting the parameters from the prepared predictions and passing those to the generator function we defined earlier.

 
posterior_predict_RTmixture <- function(i, prep, ...) { 
  if((!is.null(prep$data$lb) && prep$data$lb[i] > 0) ||  
     (!is.null(prep$data$ub) && prep$data$ub[i] < Inf)) { 
    stop("Predictions for truncated distributions not supported") 
  }   
   
  mu <- brms:::get_dpar(prep, "mu", i = i) 
  sigma <- brms:::get_dpar(prep, "sigma", i = i) 
  mix <- brms:::get_dpar(prep, "mix", i = i) 
  shiftprop <- brms:::get_dpar(prep, "shiftprop", i = i) 
   
  max_shift <- prep$data$vreal1[i] 
  upper <- prep$data$vreal2[i] 
  shift = shiftprop * max_shift 
   
  rRTmixture(prep$nsamples, meanlog = mu, sdlog = sigma,  
             mix = mix, shift = shift, upper = upper) 
} 
 

Note that the get_dpar helper that simplifies some bookeeping is currently internal in brms, but will be exposed in upcoming release.

With that, we can do a posterior predictive check for both models. We use only single core for predictions, because on Windows, multicore is slow and will not be able to access the custom prediction functions.

pp_mix <- pp_check(fit_mix, type = "dens_overlay", nsamples = 100,  cores = 1)  +
  ggtitle("Mixture")
pp_base <- pp_check(fit_base, type = "dens_overlay", nsamples = 100,  cores = 1) +
  ggtitle("Shifted lognormal")
pp_mix / pp_base

For this dataset, the mixture is not doing that much in improving the bulk of the predictions, but it manages to avoid the very long tail the lognormal-only model needs to accomodate the larger values.

We might also look at checks of the censored model. brms does not directly support predicting censored variables (because the data passed to the model are not enough to completely determine all censoring), but we can easily do this manually:

set.seed(123566)
pred_cens <- posterior_predict(fit_mix_cens, cores = 1)
pred_cens_cens <- pred_cens
# Do the censoring
pred_cens_cens[pred_cens > cens_bound] <- cens_bound 
samples_dens <- sample(1:(dim(pred_cens)[1]), size = 50)
ppc_cens1 <- ppc_dens_overlay(dd_cens$rt, pred_cens_cens[samples_dens,])  + 
  ggtitle("Censored dataset")
ppc_cens2 <- ppc_stat(1.0 * (dd_cens$censored == "right"), 
                      1.0 * (pred_cens >= cens_bound), 
                      binwidth = 0.02) + 
  ggtitle("Proportion censored")

ppc_cens1 + ppc_cens2

The model seems to do OK.

Using loo

Similarly, we might want to do model comparison or stacking with loo, so we also implement the log_lik function.

## Needed for numerical stability 
## from http://tr.im/hH5A 
logsumexp <- function (x) { 
  y = max(x) 
  y + log(sum(exp(x - y))) 
} 
 
 
RTmixture_lpdf <- function(y, meanlog, sdlog, mix, shift, upper) { 
  unif_llh = dunif(y , min = 0, max = upper, log = TRUE) 
  lognormal_llh = dlnorm(y - shift, meanlog = meanlog, sdlog = sdlog, log = TRUE) -  
    plnorm(upper - shift, meanlog = meanlog, sdlog = sdlog, log.p = TRUE) 
   
   
  # Computing logsumexp(log(mix) + unif_llh, log1p(-mix) + lognormal_llh)     
  # but vectorized 
  llh_matrix <- array(NA_real_, dim = c(2, max(length(unif_llh), length(lognormal_llh)))) 
  llh_matrix[1,] <- log(mix) + unif_llh 
  llh_matrix[2,] <- log1p(-mix) + lognormal_llh 
  apply(llh_matrix, MARGIN = 2, FUN = logsumexp) 
} 
 
log_lik_RTmixture <- function(i, draws) { 
  mu <- brms:::get_dpar(draws, "mu", i = i) 
  sigma <- brms:::get_dpar(draws, "sigma", i = i) 
  mix <- brms:::get_dpar(draws, "mix", i = i) 
  shiftprop <- brms:::get_dpar(draws, "shiftprop", i = i) 
   
  max_shift <- draws$data$vreal1[i] 
  upper <- draws$data$vreal2[i] 
  shift = shiftprop * max_shift 
   
  y <- draws$data$Y[i] 
  RTmixture_lpdf(y, meanlog = mu, sdlog = sigma,  
                 mix = mix, shift = shift, upper = upper) 
   
} 

And now, we can compare the models:

fit_mix <- add_criterion(fit_mix, "loo", cores = 1)
fit_base <- add_criterion(fit_base, "loo", cores = 1)
loo_compare(fit_mix, fit_base)
##          elpd_diff se_diff
## fit_mix   0.0       0.0   
## fit_base -7.5       5.3

No surprise here - we simulated the data with the mixture model and indeed, this is preferred to a different model. Also, the shifted-lognormal model has one very influential observation, which turns out to be the smallest observed reaction time.

dd$rt[fit_base$criteria$loo$diagnostics$pareto_k > 0.7]
## [1] 0.4909057
min(dd$rt)
## [1] 0.4909057

This once again shows that the lognormal has problem accomodating both the high and low contamination (while it is plausible it could accomodate a small amount of just high or just low contamination quite well).

Crazy models

Since brms is great, we can now do all sorts of stuff like put predictors on the mix parameter - e.g. to get a per-subject estimate of the amount of contamination.

To do this, we’ll also put a weakly informative prior on the intercept for the mixture that assumes low contamination and we don’t expect huge variability in the amount of contamination (with wider priors the model starts to diverge as we would need much more data to constrain it well).

set.seed(35486622)
dd_subj <- dd_cens
dd_subj$subject_id <- sample(1:12, size = nrow(dd_cens), replace = TRUE)

fit_mix_all <- brm(
  bf(rt | vreal(max_shift, upper) + cens(censored) ~ x, 
     mix ~ 1 + (1 | subject_id),
     family = RTmixture),
  data = dd_subj,
  stanvars = stan_funs, 
               refresh = 0,
               file = paste0(cache_dir, "/mix_all"), file_refit = "on_change",
               prior = c(prior(normal(-3, 1), class = "Intercept", dpar = "mix"),
                         prior(normal(0,0.5), class = "sd", dpar = "mix"),
                         prior(constant(1), class = "shiftprop")))

fit_mix_all
##  Family: RTmixture 
##   Links: mu = identity; sigma = identity; mix = logit; shiftprop = identity 
## Formula: rt | vreal(max_shift, upper) + cens(censored) ~ x 
##          mix ~ 1 + (1 | subject_id)
##    Data: dd_subj (Number of observations: 110) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~subject_id (Number of levels: 12) 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(mix_Intercept)     0.38      0.29     0.02     1.07 1.00     2706     2137
## 
## Population-Level Effects: 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept         0.45      0.07     0.33     0.62 1.00     3586     2305
## mix_Intercept    -2.05      0.44    -2.96    -1.25 1.00     4206     2807
## x                -0.22      0.06    -0.37    -0.11 1.00     3714     2666
## 
## Family Specific Parameters: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma         0.37      0.07     0.26     0.54 1.00     3535     3025
## shiftprop     1.00      0.00     1.00     1.00 1.00     4000     4000
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Checking that posterior predictions work:

set.seed(1233354)
pred_cens <- posterior_predict(fit_mix_all, cores = 1)
pred_cens_cens <- pred_cens
pred_cens_cens[pred_cens > cens_bound] <- cens_bound 
samples_dens <- sample(1:(dim(pred_cens)[1]), size = 50)
ppc_dens_overlay(dd_cens$rt, pred_cens_cens[samples_dens,])

We can also do multivariate models where some of the predictors are correlated across answers:

set.seed(0245562)
# Build a dataset containing two separate predictions
dd_both <- dd
dd_both$rt2 <- dd_cens$rt[1:nrow(dd_both)]
dd_both$x2 <- dd_cens$x[1:nrow(dd_both)]
dd_both$censored2 <- dd_cens$censored[1:nrow(dd_both)]
dd_both$max_shift2 <- dd_cens$max_shift[1:nrow(dd_both)]
dd_both$upper2 <- dd_cens$upper[1:nrow(dd_both)]
dd_both$subject_id <- sample(1:12, size = nrow(dd_both), replace = TRUE)

fit_mix_multivar <- brm(
  bf(rt | vreal(max_shift, upper)  ~ x, 
     mix ~ 1 + (1 | p | subject_id),
     family = RTmixture) +
  bf(rt2 | vreal(max_shift2, upper2) + cens(censored2) ~ x2, 
   mix ~ 1 + (1 | p | subject_id),
     family = RTmixture),
  data = dd_both,
  stanvars = stan_funs, 
  refresh = 0,
  file = paste0(cache_dir, "/mix_multivar"), file_refit = "on_change",
  prior = c(prior(normal(-3, 1), class = "Intercept", dpar = "mix", resp = "rt"),
           prior(normal(0,0.5), class = "sd", dpar = "mix", resp = "rt"),
           prior(constant(1), class = "shiftprop", resp = "rt"),
           prior(normal(-3, 1), class = "Intercept", dpar = "mix", resp = "rt2"),
           prior(normal(0,0.5), class = "sd", dpar = "mix", resp = "rt2"),
           prior(constant(1), class = "shiftprop", resp = "rt2")
           ),
  adapt_delta = 0.95
  )
## Setting 'rescor' to FALSE by default for this model
fit_mix_multivar
##  Family: MV(RTmixture, RTmixture) 
##   Links: mu = identity; sigma = identity; mix = logit; shiftprop = identity
##          mu = identity; sigma = identity; mix = logit; shiftprop = identity 
## Formula: rt | vreal(max_shift, upper) ~ x 
##          mix ~ 1 + (1 | p | subject_id)
##          rt2 | vreal(max_shift2, upper2) + cens(censored2) ~ x2 
##          mix ~ 1 + (1 | p | subject_id)
##    Data: dd_both (Number of observations: 100) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~subject_id (Number of levels: 12) 
##                                         Estimate Est.Error l-95% CI u-95% CI
## sd(mix_rt_Intercept)                        0.33      0.25     0.01     0.91
## sd(mix_rt2_Intercept)                       0.37      0.27     0.02     1.01
## cor(mix_rt_Intercept,mix_rt2_Intercept)     0.00      0.58    -0.95     0.95
##                                         Rhat Bulk_ESS Tail_ESS
## sd(mix_rt_Intercept)                    1.00     2712     1787
## sd(mix_rt2_Intercept)                   1.00     2587     2046
## cor(mix_rt_Intercept,mix_rt2_Intercept) 1.00     4330     2836
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rt_Intercept          0.30      0.06     0.18     0.43 1.00     4739     3085
## mix_rt_Intercept     -2.28      0.54    -3.50    -1.35 1.00     4033     1843
## rt2_Intercept         0.43      0.08     0.31     0.60 1.00     3258     1826
## mix_rt2_Intercept    -1.91      0.44    -2.80    -1.11 1.00     4455     3016
## rt_x                  0.46      0.05     0.36     0.56 1.00     5438     2978
## rt2_x2               -0.25      0.07    -0.40    -0.13 1.00     3314     2106
## 
## Family Specific Parameters: 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_rt          0.49      0.05     0.40     0.60 1.00     4904     3014
## sigma_rt2         0.36      0.07     0.25     0.53 1.00     3451     2885
## shiftprop_rt      1.00      0.00     1.00     1.00 1.00     4000     4000
## shiftprop_rt2     1.00      0.00     1.00     1.00 1.00     4000     4000
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Testing that predictions work even for multivariate models. Note that we don’t bother with censoring for rt2 so the predictions look wrong.

pp_check(fit_mix_multivar, resp = "rt", nsamples = 30, cores = 1)

pp_check(fit_mix_multivar, resp = "rt2", nsamples = 30, cores = 1)

But here we’ll also have to face possibly the biggest problem with brms: that it becomes very easy to specify a model that is too complex to be well informed by the data we have or to even build a completely broken model that no amount of data will save. The data and a few settings for the “crazy” models shown above have actually had to be tweaked for them to work well for this post. So enjoy with moderation :-).

Again, if you want the complete code, packaged in a single file you can just load into your project, go to https://github.com/martinmodrak/blog/blob/master/content/post/RTmixture.R

If you encounter problems running the models that you can’t resolve yourself, be sure to ask questions on Stan Discourse and tag me (@martinmodrak) in the question!

Original computing environment

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Czech_Czechia.1250  LC_CTYPE=Czech_Czechia.1250   
## [3] LC_MONETARY=Czech_Czechia.1250 LC_NUMERIC=C                  
## [5] LC_TIME=Czech_Czechia.1250    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bayesplot_1.8.0 patchwork_1.1.1 knitr_1.31      forcats_0.5.1  
##  [5] stringr_1.4.0   dplyr_1.0.4     purrr_0.3.4     readr_1.4.0    
##  [9] tidyr_1.1.2     tibble_3.0.6    ggplot2_3.3.3   tidyverse_1.3.0
## [13] brms_2.15.0     Rcpp_1.0.6      cmdstanr_0.3.0 
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1         backports_1.2.1      plyr_1.8.6          
##   [4] igraph_1.2.6         splines_4.0.3        crosstalk_1.1.1     
##   [7] TH.data_1.0-10       rstantools_2.1.1     inline_0.3.17       
##  [10] digest_0.6.27        htmltools_0.5.1.1    rsconnect_0.8.16    
##  [13] fansi_0.4.2          magrittr_2.0.1       modelr_0.1.8        
##  [16] RcppParallel_5.0.3   matrixStats_0.58.0   xts_0.12.1          
##  [19] sandwich_3.0-0       prettyunits_1.1.1    colorspace_2.0-0    
##  [22] rvest_0.3.6          haven_2.3.1          xfun_0.22           
##  [25] callr_3.5.1          crayon_1.4.1         jsonlite_1.7.2      
##  [28] lme4_1.1-26          survival_3.2-7       zoo_1.8-8           
##  [31] glue_1.4.2           gtable_0.3.0         emmeans_1.5.4       
##  [34] V8_3.4.0             pkgbuild_1.2.0       rstan_2.26.1        
##  [37] abind_1.4-5          scales_1.1.1         mvtnorm_1.1-1       
##  [40] DBI_1.1.1            miniUI_0.1.1.1       xtable_1.8-4        
##  [43] stats4_4.0.3         StanHeaders_2.26.1   DT_0.17             
##  [46] htmlwidgets_1.5.3    httr_1.4.2           threejs_0.3.3       
##  [49] ellipsis_0.3.1       farver_2.0.3         pkgconfig_2.0.3     
##  [52] loo_2.4.1            sass_0.3.1           dbplyr_2.1.0        
##  [55] utf8_1.1.4           labeling_0.4.2       tidyselect_1.1.0    
##  [58] rlang_0.4.10         reshape2_1.4.4       later_1.1.0.1       
##  [61] munsell_0.5.0        cellranger_1.1.0     tools_4.0.3         
##  [64] cli_2.3.1            generics_0.1.0       broom_0.7.5         
##  [67] ggridges_0.5.3       evaluate_0.14        fastmap_1.1.0       
##  [70] yaml_2.2.1           processx_3.4.5       fs_1.5.0            
##  [73] nlme_3.1-149         mime_0.10            projpred_2.0.2      
##  [76] xml2_1.3.2           rstudioapi_0.13      compiler_4.0.3      
##  [79] shinythemes_1.2.0    curl_4.3             gamm4_0.2-6         
##  [82] reprex_1.0.0         statmod_1.4.35       bslib_0.2.4         
##  [85] stringi_1.5.3        highr_0.8            ps_1.5.0            
##  [88] blogdown_1.2.3       Brobdingnag_1.2-6    lattice_0.20-41     
##  [91] Matrix_1.2-18        nloptr_1.2.2.2       markdown_1.1        
##  [94] shinyjs_2.0.0        vctrs_0.3.6          pillar_1.5.0        
##  [97] lifecycle_1.0.0      jquerylib_0.1.3      bridgesampling_1.0-0
## [100] estimability_1.3     cowplot_1.1.1        httpuv_1.5.5        
## [103] R6_2.5.0             bookdown_0.21        promises_1.2.0.1    
## [106] gridExtra_2.3        codetools_0.2-16     boot_1.3-25         
## [109] colourpicker_1.1.0   MASS_7.3-53          gtools_3.8.2        
## [112] assertthat_0.2.1     withr_2.4.1          shinystan_2.5.0     
## [115] multcomp_1.4-16      mgcv_1.8-33          parallel_4.0.3      
## [118] hms_1.0.0            grid_4.0.3           coda_0.19-4         
## [121] minqa_1.2.4          rmarkdown_2.7        shiny_1.6.0         
## [124] lubridate_1.7.9.2    base64enc_0.1-3      dygraphs_1.1.1.6

All content is licensed under the BSD 2-clause license. Source files for the blog are available at https://github.com/martinmodrak/blog.