Binomial model with brms

[daSilva5]

Hi,

For odd.ratio contrasts avg_comparison is not returning compatible estimates with emmeans.

Code for reproducible example was borrowed from https://www.andrewheiss.com/blog/2023/05/15/fancy-bayes-diffs-props/#counts-and-trials-as-formula-outcomes

Load the data and run the model:

#install.packages("likert")
library(likert)

data("pisaitems", package = "likert")


reading_wide <- pisaitems %>% 

  mutate(id = 1:n()) %>%

  dplyr::select(id, country = CNT, `Comic books` = ST25Q02, Newspapers = ST25Q05) %>% 
 
  mutate(across(c(`Comic books`, Newspapers), 
                ~fct_collapse(
                  .,
                  "Rarely" = c("Never or almost never"),
                  "Sometimes" = c("A few times a year", "About once a month"),
                  "Often" = c("Several times a month", "Several times a week")))) %>% 

  mutate(across(c(`Comic books`, Newspapers), 
                ~fct_relevel(., c("Rarely", "Sometimes", "Often")))) %>% 

  mutate(country = fct_drop(country))




reading <- reading_wide %>% 
  pivot_longer(-c(id, country),
               names_to = "book_type",
               values_to = "frequency") %>% 
  drop_na(frequency)






reading_counts <- reading %>% 
  group_by(country, book_type, frequency) %>% 
  reframe(n = n()) %>% 
  group_by(country, book_type) %>% 
  mutate(total = sum(n),
         prop = n / sum(n),
non_n=total-n) %>% 
  ungroup()



often_comics_only <- reading_counts %>% 
  filter(book_type == "Comic books", frequency == "Often")

# GLM binomial
often_comics_glm_logit<-glm(cbind(n,non_n)~country,family=binomial,data=often_comics_only)


#BRMS binomial

often_comics_model_logit <- brm(
  bf(n | trials(total) ~ country),
  data = often_comics_only,
  family = binomial(link = "logit"),
  chains = 4, iter = 2000,
  refresh = 0
)

GLM contrasts match

avg_comparisons(
  often_comics_glm_logit,
  variables = list(country = "revpairwise"),
  comparison = "lnoravg",
  transform = "exp")

                               Contrast Estimate Pr(>|z|)     S 2.5 % 97.5 %
 ln(odds(Canada) / odds(Mexico))           0.407   <0.001   Inf  0.39  0.426
 ln(odds(Canada) / odds(United States))    1.408   <0.001  37.0  1.28  1.553
 ln(odds(Mexico) / odds(United States))    3.457   <0.001 496.9  3.15  3.794

Term: country
Type: response

pairs(emmeans(often_comics_glm_logit,~country,type="response"))

contrast               odds.ratio      SE  df null z.ratio p.value
 Canada / Mexico             0.407 0.00915 Inf    1 -40.003 <0.0001
 Canada / United States      1.408 0.07040 Inf    1   6.848 <0.0001
 Mexico / United States      3.457 0.16400 Inf    1  26.113 <0.0001

P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log odds ratio scale

Emmeans GLM match results with BRMS model, but avg_comparison does not.

pairs(emmeans(often_comics_model_logit,~country,type="response")) #brms model

 contrast               odds.ratio lower.HPD upper.HPD
 Canada / Mexico             0.407      0.39     0.425
 Canada / United States      1.408      1.28     1.543
 Mexico / United States      3.458      3.16     3.782

Point estimate displayed: median 
Results are back-transformed from the log odds ratio scale 
HPD interval probability: 0.95 

Contrasts for odd.ratios using avg_comparison for BRMS:

avg_comparisons(
  often_comics_model_logit, #brms model
  variables = list(country = "revpairwise"),
  comparison = "lnoravg",
  transform = "exp")

                               Contrast Estimate 2.5 % 97.5 %
 ln(odds(Canada) / odds(Mexico))               1     1      1
 ln(odds(Canada) / odds(United States))        1     1      1
 ln(odds(Mexico) / odds(United States))        1     1      1

Term: country
Type: response

I really appreciate the help with this.

[vincentarelbundock]

I do not know what this emmeans command does, but the marginaleffects command that you are executing uses this summarizing function:

lnoravg = function(hi, lo) {
  log((mean(hi) / (1 - mean(hi))) / (mean(lo) / (1 - mean(lo))))
}

You can check that the results are correct by making a simpler comparisons() call to get counterfactual predictions, and then computing the output manually using the function above:

cmp = comparisons(
  often_comics_model_logit,
  variables = list(country = "revpairwise")
)

hi = cmp$predicted_hi
lo = cmp$predicted_lo

lnoravg = function(hi, lo) {
  log((mean(hi) / (1 - mean(hi))) / (mean(lo) / (1 - mean(lo))))
}
exp(lnoravg(hi[1:3], lo[1:3]))

    [1] 1.00017

Which explains why you seeing (rounded) values of 1 in your output.

[daSilva5]

cmp = comparisons(
often_comics_model_logit,
variables = list(country = "revpairwise")
)

hi = cmp$predicted_hi
lo = cmp$predicted_lo

lnoravg = function(hi, lo) {
log((mean(hi) / (1 - mean(hi))) / (mean(lo) / (1 - mean(lo))))
}
exp(lnoravg(hi[1:3], lo[1:3]))

[1] 1.00017

Hi,
Thank you for your response. The value is still off though. If you run this same code using "often_comics_glm_logit", you will see that the predicted values are not correct, or the expected values on "response". I have seem a post in the brms forum saying that for this model "|trials(n)", to get the estimates in the response need to use linpred with transform=TRUE. epred is not doing the expected transformation.

hi and lo using epred is just the count "n".

often_comics_only%>%add_linpred_draws(often_comics_model_logit, transform=T)%>%mean_qi(.linpred)


often_comics_only%>%add_epred_draws(often_comics_model_logit, transform=T)%>%mean_qi(.epred)

cmp_glm = comparisons(
  often_comics_glm_logit,
  variables = list(country = "revpairwise")
)

hi_glm=cmp_glm$predicted_hi
lo_glm=cmp_glm$predicted_lo

lnoravg = function(hi, lo) {
  log((mean(hi) / (1 - mean(hi))) / (mean(lo) / (1 - mean(lo))))
}

exp(lnoravg(hi_glm[1:3], lo_glm[1:3]))

Kind regards,
Isis

[vincentarelbundock]

Have you tried type="link"?

[daSilva5]

Have you tried type="link"?

Hi,

type="link", is identical for both approaches, using emmeans, or avg_comparison.
Using your code above with type= "link", then doing the inverse manually make the results compatible with emmeans type="response". If not, adding the plogis in the comparison function can be an alternative to use with avg_comparison.

avg_comparison and other functions that return type response could have a way of detecting this type of model, or an option for trials=TRUE, that then change the type of predict to .linpred( transform=TRUE).

Could you add an example for brms logit in the documentation maybe?

I do not think this is a problem for the bernoulli (link = "logit"), only family binomial with |trials(). Thus I think it would be good to show the differences in approach in the documentation.

cmp = comparisons(
  often_comics_model_logit,
  variables = list(country = "revpairwise"),
  type="link"
)

hi_plogis = plogis(cmp$predicted_hi)

lo_plogis = plogis(cmp$predicted_lo)

lnoravg = function(hi, lo) {
  log((mean(hi) / (1 - mean(hi))) / (mean(lo) / (1 - mean(lo))))
}

exp(lnoravg(hi_plogis[1:3], lo_plogis[1:3]))

[1] 0.4072409

Thank you for your help.

Kind regards,

Isis

[vincentarelbundock]

@daSilva5

Yep, that makes sense! I agree that adding something to the brms vignette makes sense here.

Would you be able to provide a paragraph or two to describe the problem and show a simple example with mtcars?

I am not a brms expert and am afraid I'll write something wrong.