Bayes factors can be used as part of a Bayesian approach to null hypothesis testing. However, Bayes factors can be highly sensitive to choice of prior. Also, rejecting a point null hypothesis (e.g., \(H_0:\theta = 0.5\)) is hardly ever an interesting or meaningful conclusion. Rather, the posterior distribution provides all relevant information to make decisions about practically meaningful issues. Credible intervals are usually much more useful, and provide much more information about how large is the parameter/difference/effect than just “evidence to reject the null hypothesis”.
This RMarkdown file provides a template for you to fill in. Read the file from start to finish, completing the parts as indicated. Some code is provided for you. Be sure to run this code, and make sure you understand what it’s doing. Some blank “code chunks” are provided; you are welcome to add more (Code > Insert chunk or CTRL-ALT-I) as needed. There are also a few places where you should type text responses. Feel free to add additional text responses as necessary.
You can run individual code chunks using the play button. You can use objects defined in one chunk in others. Just keep in mind that chunks are evaluated in order. So if you call something x in one chunk, but redefine x as something else in another chunk, it’s essential that you evaluate the chunks in the proper order.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. When you are finished
You’ll need a few R packages, which you can install using install.packages
library(ggplot2)
library(dplyr)
library(knitr)
library(janitor)
knitr::opts_chunk$set(
echo = TRUE,
warning = FALSE,
message = FALSE
)(Continuing a problem from a previous assignment.)
Does uniform color give athletes an advantage over their competitors? To investigate this question, Hill and Barton (Nature, 2005) examined the records in the 2004 Olympic Games for four combat sports: boxing, tae kwon do, Greco-Roman wrestling, and freestyle wrestling. Competitors in these sports were randomly assigned to wear either a red or a blue uniform. The researchers investigated whether competitors wearing one color won significantly more often than those wearing the other color. They collected data on a total of 457 matches.
Let the parameter \(\theta\) represent the probability that a competitor wearing red wins a match. Let \(Y\) be the number of matches, out of 457, in which the competitor wearing red defeats the competitor wearing blue.
Let’s first consider a null hypothesis of \(H_0:\theta = 0.5\) and an alternative hypothesis of \(H_a:\theta \neq 0.5\). Generally in a Bayesian analysis a parameter \(\theta\) is a random variable. However, if the null hypothesis here is true then \(\theta\) is just a constant.
The interpretation of the Bayes factor as the ratio of the posterior and prior odds breaks down for a point null hypothesis like \(H_0:\theta = 0.5\) because if \(\theta\) has a continuous distribution then the probability that \(\theta\) equals 0.5 is 0, so the odds are 1/0. (Or 0/1, but then both posterior and prior odds would be 0, so the ratio would be 0/0.)
Instead, we use the interpretation of the Bayes factor as the ratio of the likelihoods.
Compute the likelihood of observing the sample data if the null hypothesis is true. (You could use simulation, but you don’t have to if the null hypothesis is true.)
# Enter your code here.
# Red beats blue 248 out of 457 times
#Use binomial likelihood function since binomial assumptions hold given sample data
lh = pbinom(248,size = 457, prob = .5)
lh## [1] 0.969395If the alternative hypothesis is true, then \(\theta\) can potentially take any value other than 0.5, so it can truly be considered a random variable. Assume a Uniform(0, 1) = Beta(1, 1) prior for \(\theta\). Use simulation to approximate the likelihood of observing the sample data if the alternative hypothesis is true.
# Enter your code here.
Nrep = 1000000
theta = rbeta(Nrep, 1, 1)
y = rbinom(Nrep, 457, theta)
lha = sum(y != 248) / Nrep
lha## [1] 0.997823Compute the Bayes Factor in favor of the alternative hypothesis as the ratio of the likelihoods from the previous parts. Which hypothesis does the Bayes Factor favor?
# Enter your code here.
bf = lha/lh
bf## [1] 1.029325Now suppose instead we had assumed a Beta(10, 10) prior for \(\theta\). Use simulation to approximate the likelihood of observing the sample data if the alternative hypothesis is true.
# Enter your code here.
Nrep = 1000000
theta = rbeta(Nrep, 10, 10)
y = rbinom(Nrep, 457, theta)
lha2 = sum(y != 248) / Nrep
lha2## [1] 0.992865Compute the Bayes Factor in favor of the alternative hypothesis as the ratio of the likelihood from the previous part and the likelihood under the null hypothesis from above Which hypothesis does the Bayes Factor favor?
# Enter your code here.
bf2 = lha2/lh
bf2## [1] 1.024211Now suppose instead we had assumed a Beta(200, 200) prior for \(\theta\). Use simulation to approximate the likelihood of observing the sample data if the alternative hypothesis is true.
# Enter your code here.
Nrep = 1000000
theta = rbeta(Nrep, 200, 200)
y = rbinom(Nrep, 457, theta)
lha3 = sum(y != 248) / Nrep
lha3## [1] 0.988225Compute the Bayes Factor in favor of the alternative hypothesis as the ratio of the likelihood from the previous part and the likelihood under the null hypothesis from above Which hypothesis does the Bayes Factor favor?
# Enter your code here.
bf3 = lha3/lh
bf3## [1] 1.019424Is the Bayes Factor sensitive to the choice of prior? Consider each of the prior distributions. Which prior distribution places the highest prior probability on \(\theta\) being close to 0.5? Far from 0.5? Which prior distribution leads to the largest Bayes Factor in favor of the alternative hypothesis? Does this make sense?
TYPE YOUR RESPONSE HERE. The Bayes Factor is sensitive to the prior when there is more confidence placed on the prior. In the cases where “low” confidence is placed on the prior, such as the Beta(1,1) prior, there is much more placement of confidence placed on what we observed in the sample data, this leads to the highest observed Bayes factor of 1.029331 in favor of alternative hypothesis. Similarly, the Beta(10,10) prior, there is much more placement of confidence placed on what we observed in the sample data, this leads to the highest observed Bayes factor of 1.024437 in favor of alternative hypothesis; similar to the Beta(1,1) case. The Bayes factor in favor of the alternative when the prior is Beta(200,200) is 0.988301/ 0.969395 = 1.019503. Observing 248 wins by the team wearing red in 457 flips is 1.019503 times more likely under the alternative hypothesis being true than under the alternative hypothesis being true. This makes sense as we are placing much more confidence on our prior assumption that the true probability of winning when wearing red is not equal to .5, but centered near .5.
A point null hypothesis is hardly ever a realistic hypothesis. Do we really think that \(\theta\) is exactly equal to 0.500000000? A more practical and realistic hypothesis can be formed by considering a range of values that are “practically equivalent” to 0.5. What this range of values is depends on the context and should be determined before collecting data.
In this context, we might say that if the probability of winning by the competitor wearing red is between 0.49 and 0.51, then any color advantage is small enough to be considered unimportant. That is, we might consider a null hypothesis of \(H_0: 0.49 \le \theta \le 0.51\) and an alternative that \(\theta\) is outside of [0.49, 0.51].
Assume a Uniform(0, 1) = Beta(1, 1) prior distribution for \(\theta\). Compute the prior odds in favor of the alternative hypothesis, the posterior odds in favor of the alternative hypothesis, and the Bayes Factor in favor of the alternative hypothesis. (You can compute probabilities using pbeta; you don’t need simulation.)
# Enter your code here.
alpha_prior = 1
beta_prior = 1
prior = pbeta(0.51, alpha_prior, beta_prior) - pbeta(0.49, alpha_prior, beta_prior)
prior_odds = (1 - prior)/prior
# data
n = 457
y = 248
# posterior
alpha_post = alpha_prior + y
beta_post = beta_prior + n - y
post = pbeta(0.51, alpha_post, beta_post) - pbeta(0.49, alpha_post, beta_post)
post_odds = (1 - post)/post
BF = post_odds/prior_odds
BF## [1] 0.2740815Assume a Beta(10, 10) prior distribution for \(\theta\). Compute the prior odds in favor of the alternative hypothesis, the posterior odds in favor of the alternative hypothesis, and the Bayes Factor in favor of the alternative hypothesis. (You can compute probabilities using pbeta; you don’t need simulation.)
# Enter your code here.
alpha_prior = 10
beta_prior = 10
prior = pbeta(0.51, alpha_prior, beta_prior) - pbeta(0.49, alpha_prior, beta_prior)
prior_odds = (1 - prior)/prior
# data
n = 457
y = 248
# posterior
alpha_post = alpha_prior + y
beta_post = beta_prior + n - y
post = pbeta(0.51, alpha_post, beta_post) - pbeta(0.49, alpha_post, beta_post)
post_odds = (1 - post)/post
BF = post_odds/prior_odds
BF## [1] 0.931962Assume a Beta(200, 200) prior distribution for \(\theta\). Compute the prior odds in favor of the alternative hypothesis, the posterior odds in favor of the alternative hypothesis, and the Bayes Factor in favor of the alternative hypothesis. (You can compute probabilities using pbeta; you don’t need simulation.)
# Enter your code here.
alpha_prior = 200
beta_prior = 200
prior = pbeta(0.51, alpha_prior, beta_prior) - pbeta(0.49, alpha_prior, beta_prior)
prior_odds = (1 - prior)/prior
# data
n = 457
y = 248
# posterior
alpha_post = alpha_prior + y
beta_post = beta_prior + n - y
post = pbeta(0.51, alpha_post, beta_post) - pbeta(0.49, alpha_post, beta_post)
post_odds = (1 - post)/post
BF = post_odds/prior_odds
BF ## [1] 1.805793Is the Bayes Factor sensitive to the choice of prior? Consider each of the prior distributions. Which prior distribution places the highest prior probability on the alternative hypothesis being true? Which prior distribution leads to the largest Bayes Factor in favor of the alternative hypothesis? Does this make sense?
TYPE YOUR RESPONSE HERE. Yes, the Bayes Factor is sensitive in the choice of prior. The prior that places the highest prior probability on the alternative being true is the Beta(200,200) prior because it stretches the probability the least of all three priors. Thus, this Beta(200,200) prior also results in the largest Bayes Factor in favor of the alternative hypothesis. This makes sense as in the case of Beta(1, 1) this prior gives various prior probabilities of theta for which the likelihood of 248/457 is basically 0, which brings the average likelihood down.
An alternative to hypothesis testing is to just use credible intervals to determine ranges of plausible values for the parameters, and to assess any hypotheses or regions of practical equivalence in light of the credible values.
Assume a Uniform(0, 1) = Beta(1, 1) prior for \(\theta\). Compute a posterior 98% credible interval for \(\theta\) given the sample data.
# Enter your code here.
theta = seq(0, 1, 0.0001) # the grid is just for plotting
# prior
alpha_prior = 1
beta_prior = 1
prior = dbeta(theta, alpha_prior, beta_prior)
# data
n = 457
y = 248
# likelihood
likelihood = dbinom(y, n, theta)
# posterior
alpha_post = alpha_prior + y
beta_post = beta_prior + n - y
posterior = dbeta(theta, alpha_post, beta_post)
qbeta(c(0.01, 0.99), alpha_post, beta_post)## [1] 0.4882458 0.5961777Assume a Beta(10, 10) prior for \(\theta\). Compute a posterior 98% credible interval for \(\theta\) given the sample data.
# Enter your code here.
theta = seq(0, 1, 0.0001) # the grid is just for plotting
# prior
alpha_prior = 10
beta_prior = 10
prior = dbeta(theta, alpha_prior, beta_prior)
# data
n = 457
y = 248
# likelihood
likelihood = dbinom(y, n, theta)
# posterior
alpha_post = alpha_prior + y
beta_post = beta_prior + n - y
posterior = dbeta(theta, alpha_post, beta_post)
qbeta(c(0.01, 0.99), alpha_post, beta_post)## [1] 0.4876717 0.5935857Assume a Beta(200, 200) prior for \(\theta\). Compute a posterior 98% credible interval for \(\theta\) given the sample data.
# Enter your code here.
theta = seq(0, 1, 0.0001) # the grid is just for plotting
# prior
alpha_prior = 200
beta_prior = 200
prior = dbeta(theta, alpha_prior, beta_prior)
# data
n = 457
y = 248
# likelihood
likelihood = dbinom(y, n, theta)
# posterior
alpha_post = alpha_prior + y
beta_post = beta_prior + n - y
posterior = dbeta(theta, alpha_post, beta_post)
qbeta(c(0.01, 0.99), alpha_post, beta_post)## [1] 0.4830347 0.5623168Are the posterior 98% credible intervals highly sensitive to the choice of prior? Based on the credible intervals, do you have evidence that \(\theta\) lies outside of [0.49, 0.51]?
TYPE YOUR RESPONSE HERE. After observing these three credible intervals we can see that the lower bounds of the 98% credible intervals only change slightly when altering the prior distribution. Additionally, the upper bound only changes slightly when going from Beta(1,1) to Beta(10,10) and again only slightly changes when utilizing a prior of Beta(200,200). Thus, because of the sample size of the data, the posterior 98% credible intervals do not appear to be highly sensitive to the choice of prior. Based on these credible intervals, we do not have evidence to say that theta lies outside of [0.49, 0.51] because all three intervals contain the values [0.49, 0.51] inside their respective posterior 98% credible intervals. Thus we do not have evidence that \(\theta\) lies outside of [0.49, 0.51].
Suppose instead that there were \(y=2480\) successes in a sample of size \(n=4570\).
Assume a Beta(200, 200) prior for \(\theta\). Compute a posterior 98% credible interval for \(\theta\) given this sample data.
# Enter your code here.
theta = seq(0, 1, 0.0001) # the grid is just for plotting
# prior
alpha_prior = 200
beta_prior = 200
prior = dbeta(theta, alpha_prior, beta_prior)
# data
n = 4570
y = 2480
# likelihood
likelihood = dbinom(y, n, theta)
# posterior
alpha_post = alpha_prior + y
beta_post = beta_prior + n - y
posterior = dbeta(theta, alpha_post, beta_post)
qbeta(c(0.01, 0.99), alpha_post, beta_post)## [1] 0.5227674 0.5556570Based on this sample, would you have evidence that \(\theta\) lies outside of [0.49, 0.51]?
TYPE YOUR RESPONSE HERE. Based on this credible interval, we do have evidence to say that theta lies outside of [0.49, 0.51] because the interval does not contain the values in the range [0.49, 0.51] inside the respective posterior 98% credible interval after updating the data to include a much larger sample size. Thus we do have evidence to saythat \(\theta\) lies outside of [0.49, 0.51].
TYPE YOUR RESPONSE HERE. 1) One important concept from this lab is that when using Bayes Factor model comparison, where the parameter has to be exactly equal to the hypothesized value, there tends to be an extreme sensitivity of the model to the choice of prior corresponding to the alternative hypothesis. 2) I am still wondering if we can do this model comparison procedure but mixing in different prior distributions (normal, poisson, etc.) along with other posterior distributions, or do they have to have their own “dual” distributions like Beta-Binomial and Normal-Normal?