modelsummary: bayesian multinomial logit

Vincent Arel-Bundock
2020-12-14

This notebook shows how to define a custom tidy method in order to use the modelsummary package to summarize bayesian multinomial logit estimated using the brms package.

To begin, we load packages, simulate data, and estimate a model:

library(brms)
library(modelsummary)
library(parameters)
library(dplyr)
 
N <- 1000
dat = list(
  Y   = factor(sample(c('pos', 'neg', 'neutral', 'dk'), replace = TRUE, size=N)),
  X1  = sample(c(0,1), size=N, replace=TRUE),
  X2  = rnorm(N, mean=10, sd=2),
  X3  = rnorm(N, mean=10, sd=2),
  X4  = rnorm(N, mean=10, sd=2),
  grp = sample(1:10, N, replace=TRUE)
)

mod = brm(
  Y ~ X1 + (1 |r| grp), 
  family = categorical(link="logit", refcat="neutral"),
  data   = dat,
  iter   = 10,
  cores  = 4)

We want to summarize this model in a table where each column corresponds to a different level of the dependent variable (“pos”, “neg”, “neutral”, “dk”). The challenge is that brms stores the name of each level as part of the term name, whereas, the modelsummary_wide needs that information in a separate column (the coef_group argument).

To see this, we can use the get_estimates function. This function is a simple wrapper which tries to extract parameter estimates sequentially using the parameters, broom, and broom.mixed packages, and then returns them in format that modelsummary can understand. Here, we have:

get_estimates(mod)
               term    estimate   pd rope.percentage       rhat
1  b_mudk_Intercept -0.05299578 0.50            0.25 153.119207
3 b_muneg_Intercept -0.32360280 0.75            0.00  10.969310
5 b_mupos_Intercept -0.12091592 1.00            0.25   1.541910
2         b_mudk_X1  0.19955939 0.75            0.00   3.068074
4        b_muneg_X1  0.19452704 0.70            0.00   2.287111
6        b_mupos_X1  0.25312436 0.75            0.25 138.679745

Notice that the term names include the “dk”, “neg”, and “pos” strings.

One workaround is to assign a new class called “custom” to the model object, and to define a tidy method to pre-process the information extracted from this class. Inside that method works in X steps:

  1. Use the parameters::parameters function to extract estimates. We use this function because it supports bayesian multi-level models, whereas the broom does not.
  2. Use the parameters::standardize_names to, well…, standardize column names.
  3. Detect group names using regular expressions and create a new column.
  4. Clean up variable labels by removing group identifiers from the term labels.
class(mod) = c("custom", class(mod))

tidy.custom <- function(x, conf.level=.95, ...) {
  out = parameters::parameters(x, ci=conf.level, verbose=FALSE)
  out = parameters::standardize_names(out, style="broom")
  out = out %>%
        mutate(group = case_when(grepl("dk", term) ~ "Don't Know",
                                 grepl("neg", term) ~ "Negative",
                                 grepl("pos", term) ~ "Positive"),
               term = gsub("dk|neg|pos", "", term))
  return(out)
}

And now we can use the modelsummary_wide function to summarize results:

modelsummary_wide(mod, statistic = "conf.int")
Model 1 Don’t Know Model 1 Negative Model 1 Positive
b_mu_Intercept -0.053 -0.324 -0.121
[-0.558, 0.253] [-0.515, 0.181] [-1.201, -0.091]
b_mu_X1 0.200 0.195 0.253
[-0.585, 0.867] [-0.217, 0.418] [-0.085, 0.726]
Num.Obs. 1000
Log.Lik. -1.433
elpd -1441.872
elpd.se 3.519
looic 2883.744
looic.se 7.037
waic 2890.144

The beautiful thing about your new class is that it can now apply to any model. Imagine that you have 2 models to summarize and that you want to display them side by side. Estimate your new model, but don’t forget to assign the custom class:

mod2 = brm(
  Y ~ X1 + X2 + X3 + (1 |r| grp), 
  family = categorical(link="logit", refcat="neutral"),
  data   = dat,
  iter   = 10,
  cores  = 4)
Warning: There were 20 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: The largest R-hat is Inf, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat
class(mod2) = c("custom", class(mod2))

models = list(mod, mod2)

modelsummary_wide(models, statistic = "conf.int")
Model 1 Don’t Know Model 1 Negative Model 1 Positive Model 2 Don’t Know Model 2 Negative Model 2 Positive
b_mu_Intercept -0.053 -0.324 -0.121 4.692 0.004 -1.274
[-0.558, 0.253] [-0.515, 0.181] [-1.201, -0.091] [0.961, 13.269] [-8.559, 11.849] [-16.773, 8.865]
b_mu_X1 0.200 0.195 0.253 -0.997 -0.979 -0.411
[-0.585, 0.867] [-0.217, 0.418] [-0.085, 0.726] [-1.569, -0.500] [-1.141, -0.030] [-1.579, 0.633]
b_mu_X2 -0.707 -0.310 0.055
[-1.065, 0.130] [-0.886, 0.099] [-0.636, 1.309]
b_mu_X3 0.119 0.389 0.058
[-0.535, 0.849] [-0.475, 1.226] [-0.466, 0.762]
Num.Obs. 1000 1000
Log.Lik. -1.433 -3.117
elpd -1441.872 -2647.271
elpd.se 3.519 56.412
looic 2883.744 5294.542
looic.se 7.037 112.825
waic 2890.144 6391.041

By default, the columns have useful and explicit titles, but you may not want to see the “Model X” label repeated over and over. Instead, we can assign empty model names " " and " " (varying the number of spaces) and use kableExtra to add spanning column labels:

library(kableExtra)

models = list(
  " " = mod,
  "  " = mod2)

modelsummary_wide(models, statistic = "conf.int") %>%
  add_header_above(c(" " = 1, "First model" = 3, "Second model" = 3))
First model
Second model
Don’t Know Negative Positive Don’t Know Negative Positive
b_mu_Intercept -0.053 -0.324 -0.121 4.692 0.004 -1.274
[-0.558, 0.253] [-0.515, 0.181] [-1.201, -0.091] [0.961, 13.269] [-8.559, 11.849] [-16.773, 8.865]
b_mu_X1 0.200 0.195 0.253 -0.997 -0.979 -0.411
[-0.585, 0.867] [-0.217, 0.418] [-0.085, 0.726] [-1.569, -0.500] [-1.141, -0.030] [-1.579, 0.633]
b_mu_X2 -0.707 -0.310 0.055
[-1.065, 0.130] [-0.886, 0.099] [-0.636, 1.309]
b_mu_X3 0.119 0.389 0.058
[-0.535, 0.849] [-0.475, 1.226] [-0.466, 0.762]
Num.Obs. 1000 1000
Log.Lik. -1.433 -3.117
elpd -1441.872 -2647.271
elpd.se 3.519 56.412
looic 2883.744 5294.542
looic.se 7.037 112.825
waic 2890.144 6391.041

You can also stack models vertically like this:

modelsummary_wide(
  list(mod, mod2), 
  statistic="conf.int",
  stacking="vertical")
Don’t Know Negative Positive
Model 1 b_mu_Intercept -0.053 -0.324 -0.121
[-0.558, 0.253] [-0.515, 0.181] [-1.201, -0.091]
Model 1 b_mu_X1 0.200 0.195 0.253
[-0.585, 0.867] [-0.217, 0.418] [-0.085, 0.726]
Model 2 b_mu_Intercept 4.692 0.004 -1.274
[0.961, 13.269] [-8.559, 11.849] [-16.773, 8.865]
Model 2 b_mu_X1 -0.997 -0.979 -0.411
[-1.569, -0.500] [-1.141, -0.030] [-1.579, 0.633]
Model 2 b_mu_X2 -0.707 -0.310 0.055
[-1.065, 0.130] [-0.886, 0.099] [-0.636, 1.309]
Model 2 b_mu_X3 0.119 0.389 0.058
[-0.535, 0.849] [-0.475, 1.226] [-0.466, 0.762]
Num.Obs. 1000
Log.Lik. -3.117
elpd -2647.271
elpd.se 56.412
looic 5294.542
looic.se 112.825
Model 1 elpd -1441.872
Model 1 elpd.se 3.519
Model 1 logLik -1.433
Model 1 looic 2883.744
Model 1 looic.se 7.037
Model 1 nobs 1000
Model 1 sigma 1.000
Model 1 waic 2890.144
waic 6391.041