{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparing three or more means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "import scipy.stats\n",
    "import statsmodels.stats.multicomp\n",
    "\n",
    "age_shf = pd.read_pickle(\"../data/processed/age_shf\")\n",
    "age_ebn = pd.read_pickle(\"../data/processed/age_ebn\")\n",
    "age_oxf = pd.read_pickle(\"../data/processed/age_oxf\")\n",
    "\n",
    "# we want just a sample; a population is more likely to be statistically significant!\n",
    "age_shf = age_shf.sample(n = 1000, replace = True, random_state = 42)\n",
    "age_ebn = age_ebn.sample(n = 1000, replace = True, random_state = 42)\n",
    "age_oxf = age_oxf.sample(n = 1000, replace = True, random_state = 42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous section we compared the means of two different groups to see if they were statistically significantly different.\n",
    "Specifically we tested to see if the mean age of residents of Eastbourne is greater than the mean age of residents of Sheffield, and it was.\n",
    "This suggests that, overall, the population of Eastbourne is older than the population of Sheffield, and so there might be some truth to the caricature that Eastbourne is a popular retirement destination.\n",
    "\n",
    "But what if we want to compare the means of more than two groups?\n",
    "Could we perform multiple $t$--tests?\n",
    "Well, no.\n",
    "It has been argued that the process of performing multiple $t$--tests is laborious because you need to perform one for each pair of means.\n",
    "These days, with the automation provided by computers, we no longer need to do the tests manually so this is no longer relevant; performing multiple $t$--tests in a loop is now a trivial task with `R` or `Python`.\n",
    "\n",
    "But, there's a more important reason we shouldn't perform multiple $t$--tests, even though we can.\n",
    "By performing multiple $t$--tests we increase the risk of making a [Type I error](08-hypothesis-testing.html#Errors-interpreting-the-results).\n",
    "Remember we can never *prove* things with statistics, only calculate how likely something is to happen by chance alone.\n",
    "\n",
    "By convention we use a threshold of 0.05, or less than a 5% chance that the results we observe are by chance.\n",
    "But at 5% there's a 1 in 20 chance we commit a Type I error.\n",
    "So if we performed 20 pairs of $t$--tests we would be almost guaranteed to make at least one Type I error!\n",
    "XKCD, as always, neatly sums this up.\n",
    "\n",
    "![Significant](https://imgs.xkcd.com/comics/significant.png)\n",
    "\n",
    "So instead of making multiple $t$--tests we perform an analysis of variance (ANOVA).\n",
    "\n",
    "## Example\n",
    "\n",
    "For this example I use the ages of a sample of 1,000 people from Sheffield and Eastbourne as before, but to compare an additional mean I include a sample of 1,000 people from Oxford.\n",
    "Oxford is one of the youngest populations overall for any city in the UK, perhaps because of the student population, so this test will see if there is a difference between the ages of Oxford, Sheffield, and Eastbourne.\n",
    "Our null hypothesis ($H_0$) is therefore that there is no statistically significant difference between the mean ages of residents in Oxford, Sheffield, and Eastbourne.\n",
    "\n",
    "## Assumptions\n",
    "\n",
    "The assumptions of ANOVA are, not surprisingly, very similar to those of a $t$--test:\n",
    "\n",
    "- Measurements/observations are independent of each other (i.e. from different people)\n",
    "- Measurements are on a continuous scale\n",
    "- If the sample size is small the variables should have a normal distribution.\n",
    "\n",
    "For our example the observations are independent (the ages of people in Sheffield do not affect the ages of people in Eastbourne or Oxford).\n",
    "Age is a continuous measure (even though it's measured here on a discrete scale we can treat it as continuous).\n",
    "We have sample sizes of 1,000 per group so do not need to check for normality, but if we did we would plot a series of [QQ plots](11-comparing-two-independent-groups.html#Assumptions) as we did for $t$--tests.\n",
    "\n",
    "## Equality of variance\n",
    "\n",
    "As with the $t$--test it is ideal if the variances of the three (or more) groups are equal, so we are looking for *equality of variance*.\n",
    "However, as a rule of thumb [ANOVA is robust to differences in variance of up to 4x](https://stats.stackexchange.com/a/56994/35674).\n",
    "That is, if the largest variance is no larger than 4x the smallest variance the ANOVA is still reliable.\n",
    "\n",
    "First, let's try a Levene test for equality of variance as we might not even need to worry about this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scipy.stats.levene(\n",
    "    age_oxf.C_AGE_NAME, age_shf.C_AGE_NAME, age_ebn.C_AGE_NAME\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the Levene statistic is statistically significant ($p << 0.01$) so we must reject the null that there is no difference between the variances, and assume that the variances are not equal.\n",
    "So is the largest variance less than 4x the smallest variance?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_var = max(age_shf.C_AGE_NAME.var(), age_ebn.C_AGE_NAME.var(), age_oxf.C_AGE_NAME.var())\n",
    "min_var = min(age_shf.C_AGE_NAME.var(), age_ebn.C_AGE_NAME.var(), age_oxf.C_AGE_NAME.var())\n",
    "print(max_var / min_var)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The maximum variance is therefore less than twice as large as the minimum variance, so our rule of thumb holds for this data set and we can interpret the ANOVA as reliable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ANOVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scipy.stats.f_oneway(\n",
    "    age_shf, age_ebn, age_oxf\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ANOVA result is statistically significant ($p << 0.01$) so there is a statistically significant difference between the mean ages of people in Sheffield, Eastbourne, and Oxford overall.\n",
    "Now let's run some post--hoc tests to see which group(s) actually differ.\n",
    "First we should transform the data into 'long' form, so one row per obvservation.\n",
    "So that we know which observation is from which city I add a '`city`' variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "age_shf = age_shf.assign(city = \"Sheffield\")\n",
    "age_ebn = age_ebn.assign(city = \"Eastbourne\")\n",
    "age_oxf = age_oxf.assign(city = \"Oxford\")\n",
    "ages = age_shf.append(age_ebn)\n",
    "ages = ages.append(age_oxf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the post--hoc tests:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "posthoc = statsmodels.stats.multicomp.MultiComparison(\n",
    "    ages.C_AGE_NAME, ages.city\n",
    ")\n",
    "\n",
    "posthoc = posthoc.tukeyhsd()\n",
    "print(posthoc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The post--hoc tests suggest that there is a difference in ages between people in Oxford and Eastbourne, and Sheffield and Eastbourne, but *not* between Oxford and Sheffield.\n",
    "This would indicate that people in Eastbourne are older, on average, than people in Oxford or in Sheffield, but that people in Oxford are not younger than people in Sheffield overall."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}