Model selection in a bayesian framework: Which model is more likely?

Suppose we have two competing models named model 0 and model 1. A natural question to ask is "which model is more likely to be have generated the data?". An appealing feature of the bayesian framework is that we can answer this question.

To answer this question we need to calculate the probability that model 0 is the true model and the probability that model 1 is the true model.

I.e. we want to calculate

for k = 0,1.

where m = k means model k is the true model.

How do we calculate these probabilities?

Using Bayes' theorem,

where the probability of the data is given by

And since,

to find each probability we simply evaluate an integral (in practice this normally requires us to carry out monte carlo integration).

Once we have these probabilities we can also calculate the Bayes Factors/Posterior Odds given by

What I really like about this approach to model selection is that it doesn't matter whether the models are nested or not, the approach is always the same. In a frequentist/classical setting testing non-nested models is a lot tougher to handle. Another great property is that as the sample size tends to infinity, with probability 1, this approach will choose the best model. In a frequentist setting even as the sample size tends to infinity we'll still be choosing the wrong models with a probability at least as large as the chosen significance level.

However, a problem with this approach is that it requires us to evaluate integrals which can be computationally expensive. But this problem is becoming less of an issue as computers become cheaper and more powerful.

Example

The following R code considers two bayesian linear regression models with a vague prior. To carry out the monte carlo integration of the posteriors we make use of the gibbs sampler (discussed in my previous post here) to obtain samples from the non-standard distributions.

See R code below

# Model selection example in a bayesian framework

# M0: y_t = theta1 + theta2 * x1^(1/3) + e_t

# M1: y_t = phi1 + phi2 * x1^(2/3) + e_t

# where e_t ~ iidN(0,1)

# want the posterior probabilities of the two models

# postprob_mj = prob(y|model j true) * priorprob_mj / p(y)

# need to calculate the integral of p(y|theta,h)p(theta,h)

# and the integral of p(y|phi,h)p(phi,h)

# recall we chose the vague improper prior p(theta,h) = p(phi,h) = 1/h

# need to sample from the posterior to get an approximation of the integral

# # # # load required package # # # #

library(mvtnorm) 

# # # # generate data # # # #

x1 <- runif(50, min = 0, max = 10)

y <- 1 + 2 * x1^(2/5) + rnorm(50, sd = 0.5)

# # # # # # # # Model 0 # # # # # # # 

Z <- cbind(1,x1^(1/3))

ZZ <- crossprod(Z,Z)

ZZI <- solve(ZZ)

theta.ols <- ZZI %*% t(Z) %*% y

resid_m0 <- y - Z %*% theta.ols

n <- length(y)

k_m0 <- ncol(Z)

v_m0 <- n - k_m0

s2_m0 <- sum(resid_m0^2)/v_m0

# # # # # # # # Model 1 # # # # # # # 

W <- cbind(1,x1^(2/3))

WW <- crossprod(W,W)

WWI <- solve(WW)

phi.ols <- WWI %*% t(W) %*% y

resid_m1 <- y - W %*% phi.ols

n <- length(y)

k_m1 <- ncol(W)

v_m1 <- n - k_m1

s2_m1 <- sum(resid_m1^2)/v_m1

# # # # define gibbs sampler function # # # #

gibbs <- function(nrDraws,m,v,s2,X) {

h_sample <- rgamma(nrDraws, v/2, v*s2/2) 

beta_sample <- matrix(0, nrow = nrDraws, ncol = ncol(X))

XX <- crossprod(X,X)

XXI <- solve(XX)

for(i in 1:nrDraws) {

beta_sample[i,] <- rmvnorm(1,m,(1/h_sample[i]) * XXI)

}

sample <- cbind(beta_sample,h_sample)

return(sample)

}

# # # get sample for model 0 using gibbs sampler # # #

nrDraws <- 10000

m0_sample <- gibbs(nrDraws,theta.ols,v_m0,s2_m0,Z)

# # # get sample for model 1 using gibbs sampler # # #

nrDraws <- 10000

m1_sample <- gibbs(nrDraws,phi.ols,v_m1,s2_m1,W)

# # # # define posterior density function # # # #

f <- function(beta,h,m,v,s2,X) {

XXI <- solve(crossprod(X,X))

x <- dmvnorm(beta,m,1/sqrt(h) * XXI) * dgamma(h, v/2, v*s2/2)

out <- prod(x)

return(out)

}

# # # # monte carlo integration of p(y|m_0) # # # #

post_m0 <- rep(NA,nrow(m0_sample))

for(i in 1:nrow(m0_sample)) {

post_m0[i] <-f(m0_sample[i,1:2],m0_sample[i,3],theta.ols,v_m0,s2_m0,Z)

}

int_m0 <- sum(post_m0) / nrow(m0_sample)

# # # # monte carlo integration of p(y|m_1) # # # #

post_m1 <- rep(NA,nrow(m1_sample))

for(i in 1:nrow(m1_sample)) {

post_m1[i] <-f(m1_sample[i,1:2],m1_sample[i,3],phi.ols,v_m1,s2_m0,W)

}

int_m1 <- sum(post_m1) / nrow(m1_sample)

# # # calculate posterior model probabilities # # #

priorprob_m0 <- 1/2

priorprob_m1 <- 1/2

py <- priorprob_m0 * int_m0 + priorprob_m1 * int_m1

postprob_m0 <- int_m0 * priorprob_m0 / py

postprob_m1 <- int_m1 * priorprob_m1 / py

# # #  calculate Bayes Factor/Posterior Odds  # # #

BF <- postprob_m1/postprob_m0