{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lab 3-3: Analysis of Variance (ANOVA)\n",
"\n",
"For this example, we are interested in how different fertilizer treatments have affected plant growth (as measured by plant height).\n",
"\n",
"Also note, you can find a nice online tutorial on [ANOVA using python here](https://reneshbedre.github.io/blog/anova.html). ANOVA is one of the most commonly used stats tools, so many video and online resources exist. If you find an especially good one, please let me know!"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import scipy.stats as stats\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# load data file (note that we are specifying a separator other than the default comma in read_csv, the symbol for \"tab\" is \"\\t\")\n",
"df = pd.read_csv(\"../data/ANOVA_fertilizer_treatment.txt\", sep=\"\\t\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we look at the dataframe we just loaded, we can see that we have one control group (no fertilizer) and three different fertilizer treatments."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Using boxplot, we can start to visually see differences between the treatment groups\n",
"df.boxplot(column=['Control', 'F1', 'F2', 'F3'], grid=False)\n",
"\n",
"# Add labels\n",
"plt.xlabel('Fertilizer Treatment')\n",
"plt.ylabel('plant height (cm)')\n",
"plt.title('Plant Height Measurements\\nFor Each Treatment Group');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**State our null and alternative hypothesis:**\n",
"\n",
"$H_0$: All groups have the same central mean\n",
"\n",
"$H_1$: The means are different from each other\n",
"\n",
"We want 95% confidence, so choose $\\alpha = 0.05$\n",
"\n",
"In this case, we perform a one-way (also called one-factor) ANOVA to determine whether our null hypothesis ($H_0$) is true or not. We can use the [scipy.stats.f_oneway()](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.f_oneway.html) function for this. This function will return the F-statistic, and a p-value. **We can reject the null hypothesis if $p < \\alpha$**"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"F-statistic = 27.46\n",
"p = 2.711994408537598e-07\n"
]
}
],
"source": [
"# stats f_oneway functions takes the groups as input and returns an F and P-value\n",
"fvalue, pvalue = stats.f_oneway(df['Control'], df['F1'], df['F2'], df['F3'])\n",
"\n",
"# print the results\n",
"print(\"F-statistic = {}\".format( np.round(fvalue,2)))\n",
"print(\"p = {}\".format( pvalue ))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our p-value is smaller than our chosen $\\alpha$, so in this case we can reject the null hypothesis."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"Take a look at the lecture notes to recall how the F-statistic is calculated.\n",
"\n",
"The p-value for this test is determined by looking up that F-statistic in the [F-distribution](https://en.wikipedia.org/wiki/F-distribution). The p-value is this example much less than 0.05, so we can reject the null. However, if we want to know more, such as which of the groups are actually different from which other groups, we need more information."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"We are going to use functions from the [statsmodels](https://www.statsmodels.org/stable/index.html) python package for the examples below."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"import statsmodels.api as sm\n",
"from statsmodels.formula.api import ols"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code will produce an ANOVA table with our results using the [statsmodels.stats.anova.anova_lm()](https://www.statsmodels.org/dev/generated/statsmodels.stats.anova.anova_lm.html) function."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
sum_sq
\n",
"
df
\n",
"
F
\n",
"
PR(>F)
\n",
"
\n",
" \n",
" \n",
"
\n",
"
C(treatments)
\n",
"
251.440000
\n",
"
3.0
\n",
"
27.464773
\n",
"
2.711994e-07
\n",
"
\n",
"
\n",
"
Residual
\n",
"
61.033333
\n",
"
20.0
\n",
"
NaN
\n",
"
NaN
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" sum_sq df F PR(>F)\n",
"C(treatments) 251.440000 3.0 27.464773 2.711994e-07\n",
"Residual 61.033333 20.0 NaN NaN"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# reshape the d dataframe suitable for statsmodels package \n",
"df_reshaped = pd.melt(df.reset_index(), id_vars=['index'], value_vars=['Control', 'F1', 'F2', 'F3'])\n",
"\n",
"# replace column names\n",
"df_reshaped.columns = ['index', 'treatments', 'value']\n",
"\n",
"# Ordinary Least Squares (OLS) model\n",
"model = ols('value ~ C(treatments)', data=df_reshaped).fit()\n",
"anova_table = sm.stats.anova_lm(model, typ=2)\n",
"\n",
"# display the results table\n",
"anova_table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compare the above to the table we made in class. Both of these python packages are valid ways to solve the basic problem."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"However, we need to apply Tukey's test to tell which groups are different from which other groups. [Read the documentation for the pairwise_tukeyhsd function here](https://www.statsmodels.org/stable/generated/statsmodels.stats.multicomp.pairwise_tukeyhsd.html)."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"from statsmodels.stats.multicomp import pairwise_tukeyhsd"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Multiple Comparison of Means - Tukey HSD, FWER=0.05 \n",
"=====================================================\n",
" group1 group2 meandiff p-adj lower upper reject\n",
"-----------------------------------------------------\n",
"Control F1 7.6 0.0 4.7771 10.4229 True\n",
"Control F2 4.8667 0.0006 2.0437 7.6896 True\n",
"Control F3 8.2 0.0 5.3771 11.0229 True\n",
" F1 F2 -2.7333 0.0599 -5.5563 0.0896 False\n",
" F1 F3 0.6 0.9324 -2.2229 3.4229 False\n",
" F2 F3 3.3333 0.0171 0.5104 6.1563 True\n",
"-----------------------------------------------------\n"
]
}
],
"source": [
"# perform multiple pairwise comparison (Tukey HSD)\n",
"m_comp = pairwise_tukeyhsd(endog=df_reshaped['value'], groups=df_reshaped['treatments'], alpha=0.05)\n",
"\n",
"# display the results table\n",
"print(m_comp)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The \"reject\" column in this table shows us that all of the fertilizer treatments are different from the control (reject=True, we can reject the null hypothesis), and that treatments F2 and F3 are different from each (reject=True) other but not from F1 (reject=False).\n",
"\n",
"Look at the boxplot graph above again. The boxplots shows the interquartile ranges and plus and minus 1.5 times the interquartile range -- these are related to the 95% confidence range but do not show it exactly. Note that the lower and upper ranges shown in the Tukey table above suggest a real difference if the 95% confidence range does not include 0. (If it includes 0, then there's a greater than 5% chance that there's really no difference at all.)\n",
"\n",
"This is the extent of what you need to complete the homework problem. I encourage you to read the details about each of these functions, review the lecture notes, or other resources you can find online. Chapter 7 of the [USGS book](https://pubs.usgs.gov/tm/04/a03/tm4a3.pdf) covers in much more detail one-factor ANOVA as well as several more tests for comparing multiple groups of samples."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.13"
}
},
"nbformat": 4,
"nbformat_minor": 4
}