Bayesian Model Averaging (BMA)
Martin Schaffhauser
2024-08-09
1 Preface
This document compiles multiple sections from different RMarkdown files.
2 Introduction
Currently the introduction is not written.
3 Bayesian Model Averaging (BMA)
Bayesian Model Averaging …
BMA can be done two ways in R.
3.1 BMA with bayestestR package
Explanation of Bayesian Model Averaging using the bayestestR package
A total of 4 models was ran to demonstrate the process of averaging the posterior distributions of parameters over multiple models both manually and via the bayestestR package.
{r run_models_bayestestR, eval=knitr::opts_chunk$get(‘run_bma_example_models’)}
model_1_bayestestR <-
stan_glm(
formula = death_in_hospital ~ per_type + age + sex +
ICU_admission_MPI_geq26,
family = binomial("logit"),
data = data,
refresh = 0,
chains = 8,
iter = 20000,
warmup = 3000,
diagnostic_file = "SAVED_OBJECTS/MODEL_DIAGNOSTIC_FILES/bma_example_model_1_bayestestR.csv",
seed = 140187,
cores = 8,
prior = list(
dist = "normal",
df = NA,
location = c(0, 0.0593464944351461, 0, 2.01356879752913),
scale = c(
0.821156881008601,
0.0236417370439603,
0.6931472,
0.549877737532996
),
autoscale = FALSE
),
prior_intercept = list(
dist = "normal",
df = NA,
location = -1.12816447279787,
scale = 0.915343324434467,
autoscale = FALSE
)
)
model_2_bayestestR <-
stan_glm(
formula = death_in_hospital ~ per_type + age + sex,
family = binomial("logit"),
data = data,
refresh = 0,
chains = 8,
iter = 20000,
warmup = 3000,
diagnostic_file = "SAVED_OBJECTS/MODEL_DIAGNOSTIC_FILES/bma_example_model_2_bayestestR.csv",
seed = 291221,
cores = 8,
prior = list(
dist = "normal",
df = NA,
location = c(0, 0.0593464944351461, 0),
scale = c(0.821156881008601,
0.0236417370439603, 0.6931472),
autoscale = FALSE
),
prior_intercept = list(
dist = "normal",
df = NA,
location = -1.12816447279787,
scale = 0.915343324434467,
autoscale = FALSE
)
)Averaging the posterior estimations over multiple models can be done with the weighted_posteriors-function of the bayestestR package.
averaged_posterior_bayestestR <-
bayestestR::weighted_posteriors(
model_1_bayestestR,
model_2_bayestestR,
verbose = T)# save and load afert estimating marginal likelihood via bridge sampling because somehow the bridgesampler cannot handle loaded model data - it must run before hand and be in RAM this way
save(model_1_bayestestR, model_2_bayestestR, averaged_posterior_bayestestR, file = "REPORTS/COMPILED_REPORT/bma_example_models_bayestestR.RData")load("REPORTS/COMPILED_REPORT/bma_example_models_bayestestR.RData")summary(model_1_bayestestR)##
## Model Info:
## function: stan_glm
## family: binomial [logit]
## formula: death_in_hospital ~ per_type + age + sex + ICU_admission_MPI_geq26
## algorithm: sampling
## sample: 136000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 159
## predictors: 5
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) -5.4 1.0 -6.7 -5.3 -4.1
## per_typePOP 0.2 0.3 -0.3 0.2 0.6
## age 0.1 0.0 0.0 0.1 0.1
## sexm 0.1 0.4 -0.3 0.1 0.6
## ICU_admission_MPI_geq26over 1.4 0.4 0.9 1.3 1.8
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.3 0.0 0.3 0.3 0.4
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 127326
## per_typePOP 0.0 1.0 169602
## age 0.0 1.0 147170
## sexm 0.0 1.0 155510
## ICU_admission_MPI_geq26over 0.0 1.0 144514
## mean_PPD 0.0 1.0 152455
## log-posterior 0.0 1.0 61575
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).summary(model_2_bayestestR)##
## Model Info:
## function: stan_glm
## family: binomial [logit]
## formula: death_in_hospital ~ per_type + age + sex
## algorithm: sampling
## sample: 136000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 159
## predictors: 4
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) -4.7 0.9 -5.9 -4.7 -3.5
## per_typePOP 0.2 0.3 -0.2 0.2 0.6
## age 0.1 0.0 0.0 0.1 0.1
## sexm -0.2 0.3 -0.7 -0.2 0.2
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.3 0.0 0.3 0.3 0.4
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 126476
## per_typePOP 0.0 1.0 155661
## age 0.0 1.0 136246
## sexm 0.0 1.0 150537
## mean_PPD 0.0 1.0 147197
## log-posterior 0.0 1.0 61239
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).Switching parameters to zero in models where they do not occur is done automatically. The output are averaged posterior distributions of all parameters weighted by the posterior model probabilities.
Here are the means of all parameter posterior distributions as a way of summarizing the output:
# summarise averaged posterior retrieved via bayestestR
averaged_posterior_bayestestR %>%
as.data.frame() %>%
summarise(
across(2:5, \(x) mean(x, na.rm=T))
) %>% exp()## per_typePOP age sexm ICU_admission_MPI_geq26over
## 1 1.204773 1.056488 0.9954626 2.31036The posterior model probabilities can also be accessed:
posterior_weights_automatic <-
attr(averaged_posterior_bayestestR, "weights") %>%
mutate(
prop_weight = (weights / sum(weights))
)
posterior_weights_automatic## Model weights prop_weight
## 1 per_type + age + sex + ICU_admission_MPI_geq26 84311 0.6199338
## 2 per_type + age + sex 51689 0.38006623.1.1 Multimodal posteriors
When multiple models are being considered in forming an averaged posterior over those models and weighing the parameters according to posterior model probabilities it is quite likely to retrieve multimodel posterior distributions for those parameters that are present in some models while not present in others. The different peaks in the posterior correspond to different models. The parameter Bayes factor estimation via the bayestestR package cannot handle comparing multimodal posterior distributions with the averaged priors over the prior model specifications as it will not be able to converge using the its internal algorithms. -> can this be achieved by calculating the BFs manually? So in cases of including or not including several parameters the for the parameter BF estimation only variables of interest that are present in all models are viable for this estimation. The other parameters only work as controlling variables to control confounding effects taking uncertainty with regards to the real data generating process (model specification) into account.
Error converging: Error in oldlogspline(x) : * no convergence -> Troubleshooting it seems like considering multiple models in an averaged prior and posterior leads to wonky distributions with multiple spikes in parameter space -> this corresponds to different values the parameter can take on in different contraints from differing model formulas of corse - so there might be a spike at some OR values for one model and another spike somewhere else. this seems to make parameter bayes factor calculation impossible via the bayestestR package -> can manual calculation help?
See the multimodal posterior distributions for parameters that occur in some models but not in others:
bayestestR::hdi(averaged_posterior_bayestestR, ci=.95) %>% plot()
3.2 BMA manually
Explanation of Bayesian Model Averaging - manual calculation
A total of 4 models was ran to demonstrate the process of averaging the posterior distributions of parameters over multiple models both manually and via the bayestestR package.
model_1_manual <-
stan_glm(
formula = death_in_hospital ~ per_type + age + sex +
ICU_admission_MPI_geq26,
family = binomial("logit"),
data = data,
refresh = 0,
chains = 8,
iter = 20000,
warmup = 3000,
diagnostic_file = "SAVED_OBJECTS/MODEL_DIAGNOSTIC_FILES/bma_example_model_1_manual.csv",
seed = 291221,
cores = 8,
prior = list(
dist = "normal",
df = NA,
location = c(0, 0.0593464944351461, 0, 2.01356879752913),
scale = c(
0.821156881008601,
0.0236417370439603,
0.6931472,
0.549877737532996
),
autoscale = FALSE
),
prior_intercept = list(
dist = "normal",
df = NA,
location = -1.12816447279787,
scale = 0.915343324434467,
autoscale = FALSE
)
)
model_2_manual <-
stan_glm(
formula = death_in_hospital ~ per_type + age + sex,
family = binomial("logit"),
data = data,
refresh = 0,
chains = 8,
iter = 20000,
warmup = 3000,
diagnostic_file = "SAVED_OBJECTS/MODEL_DIAGNOSTIC_FILES/bma_example_model_2_manual.csv",
seed = 140187,
cores = 8,
prior = list(
dist = "normal",
df = NA,
location = c(0, 0.0593464944351461, 0),
scale = c(0.821156881008601,
0.0236417370439603, 0.6931472),
autoscale = FALSE
),
prior_intercept = list(
dist = "normal",
df = NA,
location = -1.12816447279787,
scale = 0.915343324434467,
autoscale = FALSE
)
)save(model_1_manual, model_2_manual, file = "REPORTS/COMPILED_REPORT/bma_example_models_manual.RData")load("REPORTS/COMPILED_REPORT/bma_example_models_manual.RData")summary(model_1_manual)##
## Model Info:
## function: stan_glm
## family: binomial [logit]
## formula: death_in_hospital ~ per_type + age + sex + ICU_admission_MPI_geq26
## algorithm: sampling
## sample: 136000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 159
## predictors: 5
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) -5.4 1.0 -6.7 -5.3 -4.1
## per_typePOP 0.2 0.4 -0.3 0.2 0.6
## age 0.1 0.0 0.0 0.1 0.1
## sexm 0.1 0.3 -0.3 0.1 0.6
## ICU_admission_MPI_geq26over 1.4 0.4 0.9 1.3 1.8
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.3 0.0 0.3 0.3 0.4
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 121943
## per_typePOP 0.0 1.0 165390
## age 0.0 1.0 140899
## sexm 0.0 1.0 154339
## ICU_admission_MPI_geq26over 0.0 1.0 144545
## mean_PPD 0.0 1.0 154191
## log-posterior 0.0 1.0 58701
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).summary(model_2_manual)##
## Model Info:
## function: stan_glm
## family: binomial [logit]
## formula: death_in_hospital ~ per_type + age + sex
## algorithm: sampling
## sample: 136000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 159
## predictors: 4
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) -4.7 0.9 -5.9 -4.7 -3.5
## per_typePOP 0.2 0.3 -0.2 0.2 0.6
## age 0.1 0.0 0.0 0.1 0.1
## sexm -0.2 0.3 -0.7 -0.2 0.2
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.3 0.0 0.3 0.3 0.4
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 129576
## per_typePOP 0.0 1.0 161214
## age 0.0 1.0 137277
## sexm 0.0 1.0 156652
## mean_PPD 0.0 1.0 146860
## log-posterior 0.0 1.0 61331
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).