5. Bonferroni Correction

tags: #statistics/inferential/anova

If you're testing n different hypotheses, the simplest way to avoid inflating the type I error rate is to apply a Bonferroni Correction by dividing the critical p-value by the number of pairwise comparisons to be made: $$\textrm{Bonferroni Correction} = \frac{\alpha}{\textrm{Number of Tests}}$$
For example, if $\alpha =0.05$ and you have 5 groups in your ANOVA, there will be 10 pairwise comparisons. Dividing 0.05 by 10 gives a corrected significance level of 0.005:

If the p-value for a pairwise comparison is less than 0.005, you can reject the null hypothesis.

Bonferroni Correction in Python

To apply the Bonferroni correction in Python for an ANOVA, you can use the multipletests() function from the statsmodels.stats.multitest module.

Two-way ANOVA:

To apply the Bonferroni Correction to a Two-way ANOVA in Python:

import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.multitest import multipletests

# load data
data = pd.read_csv('data.csv')

# create model
model = sm.formula.ols('Y ~ A + B + A:B', data=data).fit()

# calculate ANOVA table
anova_table = sm.stats.anova_lm(model, typ=2)

# apply Bonferroni correction to p-values
pvalues = anova_table['PR(>F)']
corrected_pvalues = multipletests(pvalues, alpha=0.05, method='bonferroni')[1]
# multipletests returns a list of values, we are [1] to index for the second element which are the pvals_corrected that is returned 

# add corrected p-values to the ANOVA table
anova_table['Bonferroni-corrected p-value'] = corrected_pvalues

print(anova_table)

One-way ANOVA:

To perform a Bonferroni Correction on a one-way ANOVA:

import statsmodels.api as sm
from statsmodels.formula.api import ols

# assume your data is stored in a pandas dataframe called 'df' 
# with the dependent variable 'y' and the independent variable 'x'

model = ols('y ~ x', data=df).fit()
anova_table = sm.stats.anova_lm(model, type=2)

# apply Bonferroni correction to p-values
pvalues = anova_table['PR(>F)']
corrected_pvalues = multipletests(pvalues, alpha=0.05, method='bonferroni')[1]

# add corrected p-values to the ANOVA table
anova_table['Bonferroni-corrected p-value'] = corrected_pvalues

print(anova_table)