`R/weighted_posteriors.R`

`weighted_posteriors.Rd`

Extract posterior samples of parameters, weighted across models.
Weighting is done by comparing posterior model probabilities, via `bayesfactor_models`

.

weighted_posteriors(..., prior_odds = NULL, missing = 0, verbose = TRUE) # S3 method for data.frame weighted_posteriors(..., prior_odds = NULL, missing = 0, verbose = TRUE) # S3 method for stanreg weighted_posteriors( ..., prior_odds = NULL, missing = 0, verbose = TRUE, effects = c("fixed", "random", "all"), component = c("conditional", "zi", "zero_inflated", "all"), parameters = NULL ) # S3 method for brmsfit weighted_posteriors( ..., prior_odds = NULL, missing = 0, verbose = TRUE, effects = c("fixed", "random", "all"), component = c("conditional", "zi", "zero_inflated", "all"), parameters = NULL ) # S3 method for blavaan weighted_posteriors( ..., prior_odds = NULL, missing = 0, verbose = TRUE, effects = c("fixed", "random", "all"), component = c("conditional", "zi", "zero_inflated", "all"), parameters = NULL ) # S3 method for BFBayesFactor weighted_posteriors( ..., prior_odds = NULL, missing = 0, verbose = TRUE, iterations = 4000 )

... | Fitted models (see details), all fit on the same data, or a single |
---|---|

prior_odds | Optional vector of prior odds for the models compared to the first model (or the denominator, for |

missing | An optional numeric value to use if a model does not contain a parameter that appears in other models. Defaults to 0. |

verbose | Toggle off warnings. |

effects | Should results for fixed effects, random effects or both be returned? Only applies to mixed models. May be abbreviated. |

component | Should results for all parameters, parameters for the conditional model or the zero-inflated part of the model be returned? May be abbreviated. Only applies to brms-models. |

parameters | Regular expression pattern that describes the parameters that
should be returned. Meta-parameters (like |

iterations | For |

A data frame with posterior distributions (weighted across models) .

Note that across models some parameters might play different roles. For example,
the parameter `A`

plays a different role in the model `Y ~ A + B`

(where it is a main effect)
than it does in the model `Y ~ A + B + A:B`

(where it is a simple effect). In many cases centering
of predictors (mean subtracting for continuous variables, and effects coding via `contr.sum`

or
orthonormal coding via `contr.orthonorm`

for factors) can reduce this issue. In any case
you should be mindful of this issue.

See `bayesfactor_models`

details for more info on passed models.

Note that for `BayesFactor`

models, posterior samples cannot be generated from intercept only models.

This function is similar in function to `brms::posterior_average`

.

For `BayesFactor < 0.9.12-4.3`

, in some instances there might be
some problems of duplicate columns of random effects in the resulting data
frame.

Clyde, M., Desimone, H., & Parmigiani, G. (1996). Prediction via orthogonalized model mixing. Journal of the American Statistical Association, 91(435), 1197-1208.

Hinne, M., Gronau, Q. F., van den Bergh, D., and Wagenmakers, E. (2019, March 25). A conceptual introduction to Bayesian Model Averaging. doi: 10.31234/osf.io/wgb64

Rouder, J. N., Haaf, J. M., & Vandekerckhove, J. (2018). Bayesian inference for psychology, part IV: Parameter estimation and Bayes factors. Psychonomic bulletin & review, 25(1), 102-113.

van den Bergh, D., Haaf, J. M., Ly, A., Rouder, J. N., & Wagenmakers, E. J. (2019). A cautionary note on estimating effect size.

`bayesfactor_inclusion`

for Bayesian model averaging.

# \donttest{ if (require("rstanarm") && require("see")) { stan_m0 <- stan_glm(extra ~ 1, data = sleep, family = gaussian(), refresh = 0, diagnostic_file = file.path(tempdir(), "df0.csv") ) stan_m1 <- stan_glm(extra ~ group, data = sleep, family = gaussian(), refresh = 0, diagnostic_file = file.path(tempdir(), "df1.csv") ) res <- weighted_posteriors(stan_m0, stan_m1) plot(eti(res)) }#>#> Warning: Bayes factors might not be precise. #> For precise Bayes factors, it is recommended sampling at least 40,000 posterior samples.#>#> Error: Failed at retrieving data :( Please provide original model or data through the `data` argument## With BayesFactor if (require("BayesFactor")) { extra_sleep <- ttestBF(formula = extra ~ group, data = sleep) wp <- weighted_posteriors(extra_sleep) describe_posterior(extra_sleep, test = NULL) describe_posterior(wp$delta, test = NULL) # also considers the null }#>#>#>#>#> #> #> #>#> Summary of Posterior Distribution #> #> Parameter | Median | 95% CI #> ---------------------------------- #> Posterior | -0.12 | [-1.32, 0.07]## weighted prediction distributions via data.frames if (require("rstanarm")) { m0 <- stan_glm( mpg ~ 1, data = mtcars, family = gaussian(), diagnostic_file = file.path(tempdir(), "df0.csv"), refresh = 0 ) m1 <- stan_glm( mpg ~ carb, data = mtcars, family = gaussian(), diagnostic_file = file.path(tempdir(), "df1.csv"), refresh = 0 ) # Predictions: pred_m0 <- data.frame(posterior_predict(m0)) pred_m1 <- data.frame(posterior_predict(m1)) BFmods <- bayesfactor_models(m0, m1) wp <- weighted_posteriors(pred_m0, pred_m1, prior_odds = BFmods$BF[2] ) # look at first 5 prediction intervals hdi(pred_m0[1:5]) hdi(pred_m1[1:5]) hdi(wp[1:5]) # between, but closer to pred_m1 }#> Warning: Bayes factors might not be precise. #> For precise Bayes factors, it is recommended sampling at least 40,000 posterior samples.#>#> Warning: 'prior_odds = NULL'; Using uniform priors odds. #> For weighted data frame, 'prior_odds' should be specified as a numeric vector.#> Highest Density Interval #> #> Parameter | 95% HDI #> ---------------------------------- #> Mazda.RX4 | [ 6.57, 30.14] #> Mazda.RX4.Wag | [ 6.91, 31.09] #> Datsun.710 | [10.40, 34.08] #> Hornet.4.Drive | [ 9.77, 35.22] #> Hornet.Sportabout | [ 9.11, 32.21]# }