# Simple example that shows the advantages of Bayesian Model Averaging (BMA)

I'm incorporating a Bayesian Model Averaging (BMA) approach in my research and will soon give a presentation about my work to my colleagues. However, BMA isn't really that well-known in my field, so after presenting them with all the theory and before actually applying it to my problem, I want to present a simple, yet instructive example on why BMA works.

I was thinking about a simple example with two models one can choose from, but the true data generating model (DGM) is somewhere in between and the evidence doesn't really favor any one of them. So if you choose one and continue from them on, you would ignore model uncertainty and make an error, but BMA, although the true model is not part of the model set, at least gives correct posterior density of the parameter of interest. For instance, there are two weather forecasts each day (A and B) and one wants to predict the weather best, so in classical statistics you would first try to find the best forecaster between the two, but what if the truth is somewhere in between (that is, sometimes A is right, sometimes B). But I couldn't formalize it. Something like that but I'm very open to ideas. I hope this question is specific enough!

In the literature, I haven't found any nice examples from what I have read so far:

• Kruschke (2011), while a great introduction to Bayesian statistics, doesn't really focus on BMA and the coin toss example he has in chapter 4 is great for introducing Bayesian statistics, but doesn't really convince a fellow researcher to use BMA. ("Why again do I have three models, one saying the coin is fair and two saying it's biased in either direction?")
• All the other stuff I read (Koop 2003, Koop/Poirier/Tobias (2007), Hoeting et al. (1999) and tons of others) are great references, but I haven't found a simple toy example in them.

But maybe I just missed a good source here.

So does anyone have a good example he or she uses to introduce BMA? Maybe by even showing the likelihoods and posteriors because I think that would be quite instructive.

• A short update: I just came across this presentation which refers to Freedman's paradox in section 2. There is a short example in which 39 random covariates are simulated and if one just looks for the best model, one eventually finds significant covariates. Model averaging is apparantly a cure for that problem. I'm not posting a solution with code here because frankly, I don't know how the figures there are derived. Oct 10, 2013 at 15:25
• (Continued) What exactly are they averaging over? The best parameter? All parameters (I think that would only make sense in this specific example). Still, I think the charts in combination with the hint to Freedman's paradox is quite helpful. Maybe it helps some. Oct 10, 2013 at 15:26

I did something similar recently. Not so much trying to convince others, but doing a small project which allowed me to get a little taste of BMA. What I did was to generate a dataset with a binary response, three independent variables which had an effect on the response and seven variables which did not have any effect on the response. I then compared the BMA results to the frequentist estimates in logistic regression. I think that at least in this case the BMA approach appears to be quite good. If you want to make it more accessible you can always name the variables or something instead of calling them the generic $X$ and $y$.

The R code I used for this is presented below. Hope it can inspire you!

# The sample size
n <- 100

# The 'true' coefficient vector
Beta <- cbind(c(-1.5, 0.45, -3))

# Generate the explanatory variables which have an effect on the outcome
set.seed(1)
X <- cbind(rnorm(n, 0, 1), rnorm(n, 4, 2), rnorm(n, 0.5, 1))

# Convert this into probabilities
prob <- 1/(1+exp(-X %*% Beta))

# Generate some uniform numbers. If the elements are smaller than the corresponding elements in the prob vector, then return 1.
set.seed(2)
runis <- runif(n, 0, 1)
y <- ifelse(runis < prob, 1, 0)

X <- cbind(X, rpois(n, 3))        # Redundant variable 1 (x4)
X <- cbind(X, rexp(n, 10))        # Redundant variable 2 (x5)
X <- cbind(X, rbeta(n, 3, 10))    # Redundant variable 3 (x6)
X <- cbind(X, rbinom(n, 10, 0.5)) # Redundant variable 4 (x7)
X <- cbind(X, rpois(n, 40))       # Redundant variable 5 (x8)
X <- cbind(X, rgamma(n, 10, 20))  # Redundant variable 6 (x9)
X <- cbind(X, runif(n, 0, 1))     # Redundant variable 7 (x10)

# The BMA
library(BMA)
model <- bic.glm(X, y,  glm.family="binomial", factor.type=FALSE, thresProbne0 = 5, strict = FALSE)

# The frequentist model
model2 <- glm(y~X, family = "binomial")

old.par <- par()
par(mar=c(3,2,3,1.5))
plot(model, mfrow=c(2,5))
par(old.par)

summary(model)
summary(model2)

• That's a nice example, so +1 for it. However, as you already pointed out, it does not really help in convincing others to use BMA. Actually, I run it and it even took some time to convince me that BMA is better here than the classical approach: the best model is not the true model (the best model only includes x2 and x3) and the parameters for model2 aren't that much off, at least for the relevant parameters. However, it shows some significant parameters x5 and x6 that shouldn't be there and BMA does a great job in telling you that this is not significant, so this is a plus for BMA. Oct 7, 2013 at 10:28

A great resource for this is:
Bayesian Model Averaging with BMS by Stefan Zeugner (2012)

It is using the R-package BMS, more info can be found here:
http://bms.zeugner.eu/

Two hands-on tutorials for reproducing real-world examples with the package can be found here:

A more general motivational and current introduction to Bayesian methods is the following paper:

The Time Has Come: Bayesian Methods for Data Analysis in the Organizational Sciences by John K. Kruschke, Herman Aguinis and Harry Joo

• Thanks for the links, but they are not really what I'm after. I actually knew and used the package before (it's great) and I agree that their documentation is really instructive. But again, the authors intention is not to convince someone (hopefully in less than 5 minutes) why they should use BMA, but given that they want to use it, how to do it with their package. So they start with the attitude example and if you scroll through your first link, there really isn't any table or figure where you would scream: "Geez, I'm happy to have used BMA!" Sep 23, 2013 at 10:34
• Continued: Just to be clear, this is of course not in any way a critique of their documentation: it's not their intention in the first place. Maybe to give an example I'm after. Let's assume you want to explain the problem of outliers in a linear regression. You would probably start with somethink like the charts here. Of course, the problem in real data will never be that easy. It's gonna be hard how you define an outlier etc. But with such a chart, everyone knows what's going on. Sep 23, 2013 at 10:36
• @Christoph_J: Do you know this paper: indiana.edu/~kruschke/articles/KruschkeAJ2012.pdf - it is not so much about BMA but about convincing someone to use Bayesian methods in the first place - perhaps this is something that is helpful for you :-) Sep 23, 2013 at 12:19
• I didn't know that one and it is a really good introductory paper to Bayesian methods, so +1 for the link. Thanks. Sep 23, 2013 at 18:53
• @Christoph_J: I edited the post accordingly: Your +1 for the answer did not work, it is still 0 (?!?) so could you please do it again - Thank you :-) Sep 24, 2013 at 11:28