Calculating probabilities after bayesian regression

You could try to automate the process of selecting variable combinations that you want to predict at. Here's an example with some built-in data. To create the prediction data frame pred, we get three quantiles for each of the numeric columns and all of the unique categories for the categorical columns. We then combine these into a data frame containing every combination of these values and run predict to get the probabilities.

library(arm)
library(vcd) # For Baseball data
library(tidyverse)

# Get a data subset for modeling
d = Baseball %>% 
  select(league86, league87, atbat86, homer86, posit86) %>% 
  filter(posit86 %in% c("1B","2B","C")) 

m1 = bayesglm(league87 ~ atbat86 + homer86 + posit86 + league86, data=d, 
              family=binomial,
              drop.unused.levels=FALSE)

# Create data frame with desired combinations of independent variables then add predictions
pred = d %>% 
  select(-league87) %>% 
  map(
    ~if(is.numeric(.x)) quantile(.x, prob=c(0.25, 0.5, 0.75)) else unique(.x)
  ) %>% 
  expand.grid() %>% 
  mutate(predicted = predict(m1, ., type="response")) 

head(pred)
#>   league86 atbat86 homer86 posit86  predicted
#> 1        N  255.00       4      2B 0.89825874
#> 2        A  255.00       4      2B 0.07354152
#> 3        N  371.00       4      2B 0.92800690
#> 4        A  371.00       4      2B 0.10385790
#> 5        N  515.75       4      2B 0.95385433
#> 6        A  515.75       4      2B 0.15672025

Created on 2021-05-05 by the reprex package (v1.0.0)