Conditional_effects plot for monotonic predictors with logit link

[Famondir]

Hi,
I fitted a model with a model with monotonic predictor and came across two problems.

binomial_base_model_fs_mono2 <- brm(
  data = success_moduls %>% mutate(SemesterOrd = ordered(Semester)),
  n_Prüfungsstatus | trials(n_Prüfungen) ~ 1 + mo(SemesterOrd),
  family=binomial(link='logit'),
  cores = 4,
  sample_prior = TRUE,
  file = "models/binomial_base_model_fs_mono2.rds"
)

1. lineplot vs catergorial plot

When variable Semester is numeric the conditional_effects plot shows a linegraph what was not expected:

Queston 1: Is there a easy way to get a categorial plot with errorbars like this?:

I had to refit the model and change the variable type to ordered to get this.

I tried to use the categorical, ordinal and re_formula options from conditional_effect but had no luck.

2. Calculation of conditional effects

When I tried to do a plot of this with tidybayes stat_halfeye style I came across a second problem.

I found this post that links to the paper where the calculation of the monotonic effects is described. When I follow the description I end up with a different plot than the conditional_effects function gave (D = 5):

The intercept / first category seems to be correct. But afterwards it does not work out well.

When I set D to 1 instead I get this:

The big step between category 5 and 6 is not present in both plots. I wonder if the D parameter has to be variable or if the terms have to be added differently with a logit link function.

Queston 2: Can someone explain me the difference or what I have to change in my calculation for the conditional effects?

I used this transformation:

D <- 5

binomial_base_model_fs_mono2 %>%
  spread_draws(b_Intercept, bsp_moSemesterOrd,
               simo_moSemesterOrd1[semester]) %>%
  pivot_wider(names_from = semester, values_from = simo_moSemesterOrd1) %>% 
  transmute(
    .draw = .draw,
    `1` = inv_logit_scaled(b_Intercept),
    `2` = inv_logit_scaled(b_Intercept+bsp_moSemesterOrd*D*(`1`)),
    `3` = inv_logit_scaled(b_Intercept+bsp_moSemesterOrd*D*(`1`+`2`)),
    `4` = inv_logit_scaled(b_Intercept+bsp_moSemesterOrd*D*(`1`+`2`+`3`)),
    `5` = inv_logit_scaled(b_Intercept+bsp_moSemesterOrd*D*(`1`+`2`+`3`+`4`)),
    `6` = inv_logit_scaled(b_Intercept+bsp_moSemesterOrd*D*(`1`+`2`+`3`+`4`+`5`))
  ) %>% pivot_longer(-.draw) %>% 
  # mutate(semester = as.character(semester+1)) %>% 
  ggplot(aes(y = name, x = value)) +
  tidybayes::stat_halfeye() +
  xlab("Success probability") +
  ylab("Semester") +
  coord_flip()

The model summary is:

 Family: binomial 
  Links: mu = logit 
Formula: n_Prüfungsstatus | trials(n_Prüfungen) ~ 1 + mo(SemesterOrd) 
   Data: success_moduls %>% mutate(SemesterOrd = ordered(Se (Number of observations: 4501) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept        -0.13      0.01    -0.16    -0.11 1.00     3516     2899
moSemesterOrd     0.31      0.01     0.29     0.33 1.00     3743     2546

Monotonic Simplex Parameters:
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
moSemesterOrd1[1]     0.14      0.01     0.12     0.17 1.00     3433     2681
moSemesterOrd1[2]     0.12      0.01     0.09     0.15 1.00     3431     2421
moSemesterOrd1[3]     0.09      0.01     0.06     0.12 1.00     3677     2481
moSemesterOrd1[4]     0.17      0.03     0.10     0.23 1.00     4290     2386
moSemesterOrd1[5]     0.48      0.03     0.41     0.55 1.00     4142     2636

Draws were sampled using sampling(NUTS). 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).

[zacho]

g’day @Famondir ,
For your first question, it doesn’t seem like it. I couldn’t find any arguments for it, I just confirmed what you found in that conditional_effects just responds to how the monotonic predictor is formatted. If it’s integer or numeric, then you get the line + ribbon. I don’t know for sure, but I think this might be to respect the cases where the predictor really is continuous – the use of mo() is akin to monotonic regression splines.

For your 2nd question:
First, D isn’t variable, it’s defined as per the number of categories minus one. So D=5 should be correct for your example.
I think you were on the right path. The key error I see in the code is that, for each category’s posterior prediction, you add up in the last set of parentheses the previous categories’ response scale predictions; rather, you want the simplex parameters in there. (The \zeta_i in that excellent Bürkner and Charpentier 2020 paper you found; look at eqn 3.) These are the moSemesterOrd1[i] bits in the model (they sum up to 1, hence the highest category gets the full effect of moSemesterOrd).

So, for Semesters 2 and 3 in your model, it should be:

`2` = inv_logit_scaled(b_Intercept+bsp_moSemesterOrd*D*(moSemesterOrd1[1]))
`3` = inv_logit_scaled(b_Intercept+bsp_moSemesterOrd*D*(moSemesterOrd1[1] + moSemesterOrd1[2]))

Here’s a reproducible example with a toy dataset where the conditional_effects output is (more or less) done the manual way:

reprex

library(brms)
library(tidyverse)
library(tidybayes)
library(patchwork) # to look at plots side-by-side

# example data from ?brms::mo

# generate some data
income_options <- c("below_20", "20_to_40", "40_to_100", "greater_100")
income <- factor(sample(income_options, 100, TRUE),
                 levels = income_options, ordered = TRUE)
mean_ls <- c(30, 60, 70, 75)
ls <- mean_ls[income] + rnorm(100, sd = 7)
dat <- data.frame(income, ls)

# plot raw data:
ggplot(dat, aes(income, ls))+geom_point()

# fit a simple monotonic model
get_prior(ls ~ mo(income), data = dat)
fit1 <- brm(ls ~ mo(income), data = dat, chains = 4, cores = 4, backend = 'cmdstanr')
summary(fit1)
p_ce <- plot(conditional_effects(fit1))

# now try manually with spread_draws:
get_variables(fit1) # what's available

p_manual_ce <- fit1 %>% 
  spread_draws(b_Intercept, bsp_moincome,
               simo_moincome1[i]) %>% 
  # simplex parameters are grouped by variable `i`, pivot into same row
  pivot_wider(names_from = `i`, names_prefix = 'simo_moincome1_', values_from = simo_moincome1) %>% 
  mutate(
    # D is equal to number of categories minus 1
    D = 3,
    # for each level of income, get posterior mu
    mu_cat_1 = b_Intercept,
    mu_cat_2 = b_Intercept + (bsp_moincome * D * (simo_moincome1_1)),
    mu_cat_3 = b_Intercept + (bsp_moincome * D * (simo_moincome1_1 + simo_moincome1_2)),
    mu_cat_4 = b_Intercept + (bsp_moincome * D * (simo_moincome1_1 + simo_moincome1_2 + simo_moincome1_3))
  ) %>% 
  # at this point or after pivoting again, you could apply the inverse link function on mu
  pivot_longer(starts_with('mu_cat'), names_to = 'mu_cat', values_to = 'post_mu') %>% 
  ggplot(aes(y = post_mu, x = mu_cat))+
  stat_halfeye()

(p_ce$income + scale_y_continuous(limits = c(20, 81))) + 
  (p_manual_ce + scale_y_continuous(limits = c(20, 81)))

# looks equivalent

[Famondir]

Thank you very much, zacho. Great that you even gave a short reprex.

I didn’t recognized that transmute (as mutate does as well) does the steps one by one so 1 becomes altered for the second step because I reuse the column names. I changed it and now it works as expected.

I realy appretiate your response!