I’ll follow @Bob_Carpenter’s suggestion to throw all the draws together, and I’ll use the example code you provided to show how to do this with the beta parameters.
First, it helps to see how to do it with one model. The betas from one model look like this:
> str(anes_2012_models[[1]]$beta)
num [1:10, 1:1000] 0.0398 -0.0173 0.087 -0.1759 0.31 ...
The dimension for iterations (1:1000) is last, which is the opposite of how {tidybayes} and {posterior} expect matrices of draws to be structured. So, we have to transpose this matrix with t() to move that dimension to the front. Given the transposed matrix, you can use either tidybayes::tidy_draws() or posterior::as_draws_df() to convert into a data frame of draws; e.g.:
> posterior::as_draws_df(t(anes_2012_models[[1]]$beta))
# A draws_df: 1000 iterations, 1 chains, and 10 variables
...1 ...2 ...3 ...4 ...5 ...6 ...7 ...8
1 0.040 -0.0173 0.08698 -0.176 0.310 0.0821 -0.0243 -0.055
2 -0.073 -0.0982 -0.04095 -0.229 -0.669 0.1062 -0.0385 -0.115
3 0.091 -0.0411 0.04597 -0.182 0.301 0.2498 -0.1360 -0.059
4 -0.077 -0.0164 0.04660 0.072 0.089 -0.0632 0.0055 -0.342
5 -0.067 0.0141 -0.00810 0.213 0.803 -0.0118 -0.0289 -0.048
6 -0.071 -0.0266 0.09155 0.046 -0.148 -0.0174 0.0114 -0.066
7 -0.045 -0.0702 0.00095 -0.197 -0.200 0.0243 0.0332 -0.044
8 -0.098 0.0047 0.04084 -0.520 0.262 -0.0028 -0.0028 0.012
9 -0.189 -0.0459 -0.06672 -0.195 -0.238 0.2645 -0.1209 -0.346
10 0.100 -0.0017 0.08825 0.156 -0.169 0.0696 0.0050 -0.159
# ... with 990 more draws, and 2 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
As for throwing all the draws together from each model, there are plenty of ways to do it, depending on your comfort with different data structures. If you are comfortable with manipulating matrices in R, I’d probably just bind all the matrices together upfront and then convert to a data frame. Something like this:
> betas_matrix = do.call(cbind, lapply(anes_2012_models, `[[`, "beta"))
> str(betas_matrix)
num [1:10, 1:10000] 0.0398 -0.0173 0.087 -0.1759 0.31 ...
(If you haven’t seen something like lapply(anes_2012_models, "[[", "beta") before, this says to create a new list that is the result of doing x[["beta"]] for every x in anes_2012_models.)
Given betas_matrix, you then do tidybayes::tidy_draws(t(betas_matrix)) or posterior::as_draws_df(t(betas_matrix)) to get the data frame version.
Alternatively, you could convert to a list of data frames first, and then bind them together. If you do that, I would use posterior::bind_draws() instead of rbind() to do it, because you need to re-number the .draw column in order for the result to be usable by some packages (e.g., {bayesplot}). So while you could do do.call(rbind, betas_dfs) as the second step and have it mostly work, the slightly more complicated second step below is safer:
> betas_dfs = lapply(anes_2012_models, function(m) posterior::as_draws_df(t(m$beta)))
> betas_df = do.call(posterior::bind_draws, c(betas_dfs, list(along = "draw")))
> betas_df
# A draws_df: 10000 iterations, 1 chains, and 10 variables
...1 ...2 ...3 ...4 ...5 ...6 ...7 ...8
1 0.040 -0.0173 0.08698 -0.176 0.310 0.0821 -0.0243 -0.055
2 -0.073 -0.0982 -0.04095 -0.229 -0.669 0.1062 -0.0385 -0.115
3 0.091 -0.0411 0.04597 -0.182 0.301 0.2498 -0.1360 -0.059
4 -0.077 -0.0164 0.04660 0.072 0.089 -0.0632 0.0055 -0.342
5 -0.067 0.0141 -0.00810 0.213 0.803 -0.0118 -0.0289 -0.048
6 -0.071 -0.0266 0.09155 0.046 -0.148 -0.0174 0.0114 -0.066
7 -0.045 -0.0702 0.00095 -0.197 -0.200 0.0243 0.0332 -0.044
8 -0.098 0.0047 0.04084 -0.520 0.262 -0.0028 -0.0028 0.012
9 -0.189 -0.0459 -0.06672 -0.195 -0.238 0.2645 -0.1209 -0.346
10 0.100 -0.0017 0.08825 0.156 -0.169 0.0696 0.0050 -0.159
# ... with 9990 more draws, and 2 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
Then you should be able to use posterior::summarise_draws() on the result:
> posterior::summarise_draws(betas_df)
# A tibble: 10 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ...1 -0.00789 -0.00482 0.0911 0.0709 -0.164 0.138 1.00 1946. 9397.
2 ...2 -0.00818 -0.00400 0.0899 0.0687 -0.158 0.134 1.00 3010. 8798.
3 ...3 0.0328 0.0227 0.0865 0.0703 -0.0962 0.186 1.00 1007. 7968.
4 ...4 -0.121 -0.0857 0.265 0.213 -0.609 0.267 1.00 8863. 9647.
5 ...5 -0.0142 -0.00835 0.285 0.224 -0.491 0.444 1.01 562. 7651.
6 ...6 0.0745 0.0544 0.116 0.0962 -0.0832 0.296 1.00 1742. 7226.
7 ...7 -0.0325 -0.0234 0.0549 0.0451 -0.137 0.0434 1.01 3517. 8856.
8 ...8 -0.212 -0.166 0.255 0.225 -0.695 0.121 1.00 1011. 6938.
9 ...9 -0.000628 -0.000369 0.00484 0.00384 -0.00911 0.00709 1.00 9439. 7883.
10 ...10 -0.0246 -0.0166 0.0822 0.0661 -0.170 0.100 1.00 640. 9614.
Or pass it into {bayesplot} functions:
bayesplot::mcmc_areas(betas_df)
Or make some visualizations with tidybayes (here we go through gather_variables() to get a long-format data frame):
betas_df |>
gather_variables() |>
ggplot(aes(x = .value, y = .variable)) +
# normalize density within groups b/c variance is quite different across variables
# otherwise you can't see the densities for anything but `...9`
stat_halfeye(normalize = "groups")

