How to Calculate Bayes Factor and Credible Interval in R

What is a Bayes Factor?

The Bayes Factor is a measure of evidence in favor of one hypothesis over another. It quantifies how much more likely the data is under one hypothesis compared to another. A Bayes Factor greater than 1 indicates evidence in favor of the alternative hypothesis, while a value less than 1 indicates evidence in favor of the null hypothesis.

Where:

• P(Data|H1) is the probability of the data given the alternative hypothesis
• P(Data|H0) is the probability of the data given the null hypothesis

Bayes Factors are interpreted on a logarithmic scale:

• 0 to 1: Negative evidence against H1
• 1 to 3: Anecdotal evidence
• 3 to 10: Substantial evidence
• 10 to 100: Strong evidence
• 100+: Very strong evidence

What is a Credible Interval?

A credible interval is a Bayesian equivalent of a confidence interval. It represents the range of values that contain a specified probability of the parameter of interest. Unlike confidence intervals, credible intervals are based on the posterior distribution of the parameter.

Where:

• θ is the parameter of interest
• θₗ and θᵤ are the lower and upper bounds of the interval
• α is the credible level (typically 0.05 for 95% credible intervals)

Credible intervals provide a direct probability statement about the parameter, making them more intuitive than confidence intervals.

Calculating Bayes Factor in R

In R, you can calculate Bayes Factors using the BayesFactor package. Here's a basic example:

library(BayesFactor)

# Example data
data <- c(2.1, 2.0, 2.2, 2.1, 2.3, 2.0, 2.2, 2.1, 2.3, 2.2)

# Calculate Bayes Factor for mean difference
bf <- ttestBF(data, mu = 2.1, nullInterval = c(2.0, 2.2))
print(bf)

The ttestBF function performs a t-test and returns the Bayes Factor. The nullInterval parameter specifies the range of values for the null hypothesis.

You can also calculate Bayes Factors for other tests using functions like regressionBF for linear regression or proportionBF for proportion tests.

Calculating Credible Interval in R

To calculate credible intervals in R, you can use the rstanarm package for Bayesian regression or the brms package for more advanced Bayesian models. Here's an example using rstanarm:

library(rstanarm)

# Example data
data <- data.frame(y = c(2.1, 2.0, 2.2, 2.1, 2.3, 2.0, 2.2, 2.1, 2.3, 2.2))

# Fit Bayesian model
fit <- stan_glm(y ~ 1, data = data, prior = normal(2.1, 0.1))

# Get credible interval
summary(fit)$summary[, "5%":"95%"]

This code fits a simple Bayesian model and extracts the 95% credible interval for the mean. The prior parameter specifies the prior distribution for the parameter.

For more complex models, you can use the brms package which provides a more flexible framework for Bayesian modeling.

Example Calculation

Let's calculate both the Bayes Factor and Credible Interval for a simple example where we want to test if the mean of a dataset is different from 2.1.

Bayes Factor Calculation

Using the BayesFactor package:

library(BayesFactor)

# Sample data
data <- c(2.1, 2.0, 2.2, 2.1, 2.3, 2.0, 2.2, 2.1, 2.3, 2.2)

# Calculate Bayes Factor
bf <- ttestBF(data, mu = 2.1, nullInterval = c(2.0, 2.2))
print(bf)

The output might show a Bayes Factor of approximately 0.33, indicating weak evidence against the null hypothesis that the mean is between 2.0 and 2.2.

Credible Interval Calculation

Using the rstanarm package:

library(rstanarm)

# Fit Bayesian model
fit <- stan_glm(y ~ 1, data = data.frame(y = data), prior = normal(2.1, 0.1))

# Get credible interval
ci <- summary(fit)$summary[, "5%":"95%"]
print(ci)

The output might show a 95% credible interval of approximately [2.05, 2.25], indicating that we're 95% confident the true mean falls within this range.