Find the predictor value at which posterior_epred(model) = x (in a multilevel logistic regression)

Interfaces brms
, , ,
7

I’ve been following this thread with interest because it’s pretty similar to a problem I’ve been working on with a binomial model. I don’t know the clever / programmatic way to do this (perhaps something with the model matrix?) but it’s pretty much just algebra + patience to do the dumb method below. It uses spread_draws as @matti suggests.

Using @matti’s example fit:

library(tidybayes)

# note order of b_ coef's is same as in model matrix
coef_draws <- spread_draws(fit, `b_.*`, 
             # `sd_.*`, `cor_.*`, # for varying effects etc.
             regex = T, ndraws = 4000)

# combine with a grid of interest:
expand_grid(
  test = factor(c("a", "b")),
  trial = -5:5,
  coef_draws
) %>% 
  mutate(
    # dummy for test b
    x_b = if_else(test == "b", 1, 0),
    # the value of stim to yield logit(0.5)
    x_stim_for_0.5 = (qlogis(0.5) - b_Intercept - b_testb * x_b - b_trial * trial - 
      `b_testb:trial` * trial * x_b) / (
        b_stim + (`b_testb:stim` * x_b) + (`b_stim:trial` * trial) + 
          (`b_testb:stim:trial` * x_b * trial)
      )
  ) %>% 
  # plot
  ggplot(aes(trial, x_stim_for_0.5, fill = test, color = test)) + 
  ggdist::stat_dist_halfeye(alpha = 0.3)

image
image679×501 7.43 KB

Essentially arrange the linear model to solve for the stim predictor. I.e.:
{logit(0.5) - all the stuff not associated with stim} / {all the stuff with stim but with stim factored out}

Note this works just as well if you varying intercepts/slopes – that spread_draws call can also grab the relevant sd_ and cor_ parameters. Though if you do have cor_ terms, that means you have to simulate the varying effects from the multivariate normal distribution… which I also did in a dumb manner :^)

I’d be keen to know if there’s a more programmatic method for this (e.g. if one wanted to know the same thing but for trial without the manual re-arrangement of model terms).