Misclassification or measurement error for binomial model with brms (Statistical Rethinking 2023)

Interfaces brms
1

I am working through the 2023 version of the Statistical Rethinking course, and I have gotten stuck embarrassingly early in my attempt to use brms.

In lecture “bonus” material, McElreath asks how we could model the proportion of water on the globe but with a twist. The twist is that for x proportion of the time, the outsome is misclassified. He codes this as

sim_globe <- function(p = 0.7, N = 9, x = 0.1) {
  true_sample <- sample(c("W", "L"), size = N, prob = c(p, 1 - p), replace = TRUE)
  obs_sample <- ifelse(runif(N) < x,
                       ifelse(true_sample == "W", "L", "W"), # error
                       true_sample)
  return(obs_sample)
}

sim_globe()
#> [1] "W" "W" "L" "L" "L" "W" "W" "W" "W"

My question is: how could in model this with brms?

Here is a brms model without incorporating measurement error taken from Kurz’s brms translation of Statistical Rethinking:

library(brms)
b2.1 <-
  brm(
    w | trials(n) ~ 0 + Intercept,
    data = data. Frame(w = 24, n = 36),
    family = binomial(link = "identity"),
    prior =  prior(beta(1, 1), class = b, lb = 0, ub = 1)
  )

My initial thought was to do something like this, but here’s what I get:

b2.2 <-
  brm(
    w | mi() + trials(n) ~ 0 + Intercept,
    data = data.frame(w = 24, n = 36),
    family = binomial(link = "identity"),
    prior =  prior(beta(1, 1), class = b, lb = 0, ub = 1),
    seed = 2
  )
#> Error: Argument 'mi' is not supported for family 'binomial(identity)'.

Is there a way to properly specify measurement error with a binomial model with brms?

Session info:

  • R version 4.2.2 (2022-10-31)
  • Platform: aarch64-apple-darwin20 (64-bit)
  • Running under: macOS Ventura 13.0.1
  • brms Version: 2.18
2 Likes
2

If you know the misclassification rates, then I think the most general way to approach a problem like this is marginalize over the possible true data values according to the classification probabilities. We could achieve this in brms, if desired, via a custom family that takes the known misclassification probability as a vreal argument (note that this also allows us to pass observation-specific misclassification probabilities if these are known to vary).

I think (but I haven’t verified for myself) that we could estimate the same model in brms without writing any custom Stan code by splitting every observation in two, one with the same response and one with the opposite response, and giving a weight of 1 - x to the observation with the same response and a weight of x to the observation with the opposite response. This again assumes that the misclassification probabilities x are known.

It’s an interesting and worthwhile question!

2 Likes
3

If anyone has the code for this, I’d love to see it.

4

In the very most general case, where every pixel comes with its own misclassification probability, we would set up the linear predictor and logistic inverse link, but in place of binary 1/0 data the data would be the probability that the response is 1, and the likelihood would be

target += log_sum_exp(
  log(y[i]) + bernoulli_lpmf(1 | lin_pred[i]),
  log1m(y[i]) + bernoulli_lpmf(0 | lin_pred[i])
)

where lin_pred is the inverse logit of the linear predictor.

1 Like
5

Hi all, my previous post here is incorrect. We need to marginalize over two possibilities: that the true state is zero and that the true state is one. However, to achieve this we need likelihoods not based on y[i] the a priori probability that the true state is one, but rather on z[i], the likelihood of observing the data given that the true state is zero, and o[i], the likelihood of observing the data given that the true state is one.

To see that these are different, consider the case where y contains a thousand zeros and a hundred 0.5’s. In expectation, we think there are 1050 true zeros and 50 true ones, and of the 1050 true zeros 50 of them are misclassified. Thus, the likelihood of observing a 1 when the true state is zero is not 0.5, but rather (something like) 50/1050.

6

Does the Rethinking package let you generate the Stan code? If so, can you post it?