05_15_inference_b.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/COGS108/Lectures-Fa23/blob/master/Week_05/05_15_inference_b.ipynb)\n",
    "\n",
    "## Inference II: Building on Regression\n",
    "- Association between variables: **correlation** (& causation)\n",
    "- Visualizing multple linear regression\n",
    "- Comparison of means: **t-test**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "# Import seaborn and apply its plotting styles\n",
    "import seaborn as sns\n",
    "sns.set(style=\"white\", font_scale=2)\n",
    "\n",
    "# import matplotlib\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "# set plotting size parameter\n",
    "plt.rcParams['figure.figsize'] = (17, 7)\n",
    "\n",
    "# Statmodels & patsy\n",
    "import patsy\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "from scipy.stats import pearsonr, norm, ttest_ind\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "#improve resolution\n",
    "%config InlineBackend.figure_format ='retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlation: Pearson correlation ($r$)\n",
    "\n",
    "- measures the strength of association between two variables \n",
    "- takes values [-1,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Positive Correlation \n",
    "- $r$ > 0\n",
    "- as $r$ approaches 1, x and y are highly correlated "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "\n",
    "# 1000 random integers between 0 and 50\n",
    "x = np.random.randint(0, 50, 1000)\n",
    "  \n",
    "# Positive Correlation with some noise\n",
    "y = x + np.random.normal(0, 10, 1000)\n",
    "\n",
    "np.corrcoef(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x, y, s=50);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Negative Correlation \n",
    "- $r$ < 0\n",
    "- as $r$ approaches -1, x and y are highly negatively correlated "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.randint(0, 50, 1000)\n",
    "\n",
    "# Negative Correlation with some noise\n",
    "y = 100 - x + np.random.normal(0, 5, 1000)\n",
    "\n",
    "np.corrcoef(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x, y, s=50);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# No Correlation \n",
    "- $r \\approx 0$\n",
    "- as $r$ approaches 0, x and y are not correlated "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.randint(0, 50, 1000)\n",
    "y = np.random.randint(0, 50, 1000)\n",
    "\n",
    "np.corrcoef(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.scatterplot(x, y, s=50);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clicker Question #1\n",
    "\n",
    "What would you expect Pearson correlation for the following relationship to be?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(888)\n",
    "\n",
    "# 1000 random integers between 0 and 50\n",
    "x = np.random.randint(0, 50, 200)\n",
    "  \n",
    "# Positive Correlation with some noise\n",
    "y = x + np.random.normal(0, 10, 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.scatterplot(x, y, s=50);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- A) -0.8\n",
    "- B) -0.4\n",
    "- C) 0\n",
    "- D) 0.4\n",
    "- E) 0.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.corrcoef(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression Reminder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$outcome = \\beta_0 + \\beta_1*predictor$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\beta_1 \\approx r $ if ...\n",
    "\n",
    "- $predictor$ and $outcome$ each have a standard deviation of 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# inputs are mean, sd, and number of observations\n",
    "x = np.random.normal(0, 1, 100)\n",
    "y = np.random.normal(0, 1, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# look at the data\n",
    "sns.scatterplot(x, y, s=50);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# regression as we've done previously\n",
    "outcome, predictors = patsy.dmatrices('y ~ x')\n",
    "mod = sm.OLS(outcome, predictors)\n",
    "res = mod.fit() \n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## plot model for data\n",
    "sns.scatterplot(x, y, s=50);\n",
    "\n",
    "# Generate and plot the model fit line\n",
    "xs = np.arange(x.min(), x.max())\n",
    "ys = res.params[1] * xs + res.params[0]\n",
    "plt.plot(xs, ys, '--k', linewidth=4, label='Model');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('model beta: ' , res.params[1])\n",
    "print('corr: ', pearsonr(x,y)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Data\n",
    "\n",
    "Returning to the regression question we looked at last class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in file; specify that it is tab-separated file\n",
    "df = pd.read_csv('https://raw.githubusercontent.com/shanellis/datasets/master/index.txt', sep=\"\\t\")\n",
    "df = df[df['ViolCrime'] < 50]\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clicker Question #2\n",
    "\n",
    "What is most correlated with Viol Crime?\n",
    "\n",
    "- A) PovPct\n",
    "- B) Brth15to17\n",
    "- C) Brth18to19\n",
    "- D) TeenBrth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.corr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that each variable is perfectly correlated with itself (b/c you're comparing the exact same values!)\n",
    "\n",
    "`TeenBrth` is highly correlated with `Brth15to17` and `Brth18to19`....which makes sense as the latter are subsets of the former. \n",
    "\n",
    "Last lecture we saw that `PovPct` plays a role in `TeenBrth`, so it makes sense that these values are somewhat correlated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run model\n",
    "outcome, predictors = patsy.dmatrices('ViolCrime ~ TeenBrth + PovPct + Brth18to19 + Brth15to17', df)\n",
    "mod = sm.OLS(outcome, predictors)\n",
    "res = mod.fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run model\n",
    "outcome, predictors = patsy.dmatrices('TeenBrth ~ PovPct', df)\n",
    "mod = sm.OLS(outcome, predictors)\n",
    "res = mod.fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run model\n",
    "outcome, predictors = patsy.dmatrices('TeenBrth ~ ViolCrime + PovPct', df)\n",
    "mod = sm.OLS(outcome, predictors)\n",
    "res = mod.fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run model\n",
    "outcome, predictors = patsy.dmatrices('TeenBrth ~ ViolCrime + PovPct ', df)\n",
    "mod = sm.OLS(outcome, predictors)\n",
    "res = mod.fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### linear regression\n",
    "$\\beta_1 = 1.93$  \n",
    "\n",
    "#### correlation\n",
    "$r = 0.67$  \n",
    "$R^2 =  0.45$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$R-squared$ ($R^2$) is equal to Pearson's correlation ($r$) squared.\n",
    "\n",
    "**$R-squared$** is the percentage of the the $outcome$ variable variation that is explained by a linear model. It takes a value between 0 and 100%.\n",
    "\n",
    "Here, 45.3% of the variation in Teen Birth ($outcome$) is explained by Poverty Percentage ($predictor$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot the model fit line\n",
    "\n",
    "# Plot the orginal data (as before)\n",
    "sns.scatterplot(x='PovPct', y='TeenBrth', data=df, alpha=0.3, s=100);\n",
    "\n",
    "# Generate and plot the model fit line\n",
    "xs = np.arange(df['PovPct'].min(), df['PovPct'].max())\n",
    "ys = 1.9292 * xs + 16.7423\n",
    "plt.plot(xs, ys, '--k', linewidth=4, label='Model')\n",
    "\n",
    "plt.xlabel('Poverty Percentage')\n",
    "plt.ylabel('Teen Birth Rate')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Returning to multiple linear regression\n",
    "\n",
    "- Visualizing the relationship\n",
    "- Correlation coefficients\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.histplot(df['ViolCrime'], bins=10)\n",
    "plt.xlabel('Violent Crime');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['ViolCrime'].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "~Chat with your neighbor.~\n",
    "\n",
    "What is your interpretation of these data?\n",
    "\n",
    "Should this variable be used as is in our regression model? \n",
    "- If yes, why is that your conclusion?\n",
    "- If no, why not? What could we do to use this variable?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clicker Question #3\n",
    "\n",
    "Should this variable be used in our regression model as is? \n",
    "- A) Yes \n",
    "- B) No\n",
    "- C) I don't know"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add a new column 'Violent'\n",
    "df['Violent'] = df['ViolCrime'].apply(lambda x: 0 if x < 6.7 else 1)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.boxplot(x='Violent', y='TeenBrth', data=df);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outcome, predictors = patsy.dmatrices('TeenBrth ~ PovPct + Violent', df)\n",
    "mod = sm.OLS(outcome, predictors)\n",
    "res = mod.fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may have noted above that there is an $R-squared$ and an $Adj R-squared$ in the model output.\n",
    "\n",
    "Every time a predictor is added to a model, $R-Squared$ will increase, making you *think* you should add all the predictors.\n",
    "\n",
    "$Adj R-Squared$ is adjusted to account for the number of predictors in a model (applies in multiple regression), allowing multiple models to be compared more accurately. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### One predictor (Poverty Percentage)\n",
    "$R^2 = 0.45$  \n",
    "$Adj R^2 = 0.44$\n",
    "\n",
    "\n",
    "\n",
    "#### Two predictors (Violent and Poverty Percentage)\n",
    "$R^2 = 0.61$  \n",
    "$Adj R^2 = 0.59$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparison of the adjusted correlation coefficients suggests that even after adjusting for the fact that we're adding a second predictor, the second model explains more of the variance in the outcome variable, Teen Birth Rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clicker Question #4\n",
    "\n",
    "How are you doing? \n",
    "\n",
    "- A) totally got it\n",
    "- B) kinda got it\n",
    "- C) neutral\n",
    "- D) kinda confused\n",
    "- E) super confused"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing multiple regression models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# relationship between originalpredictor and outcome;\n",
    "sns.scatterplot(x='PovPct', y='TeenBrth', \n",
    "                data=df, s=200);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# color by Violent\n",
    "sns.scatterplot(x='PovPct', y='TeenBrth', hue='Violent', \n",
    "                data=df, s=200 );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot with regression line\n",
    "    sns.lmplot(x=\"PovPct\", y=\"TeenBrth\", hue=\"Violent\", \n",
    "           data=df, ci=None, aspect=2.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ Teen Birth = 19.8 + (1.30  * Poverty Percentage) + (10.5 * Violent) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clicker Question #5\n",
    "\n",
    "How are you doing? \n",
    "\n",
    "- A) totally got it\n",
    "- B) kinda got it\n",
    "- C) neutral\n",
    "- D) kinda confused\n",
    "- E) super confused"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cautions about correlation\n",
    "\n",
    "- Anscombe's Quartet: A Cautionary Tale\n",
    "- Correlation != Causation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Code in this example taken from [here](https://matplotlib.org/gallery/specialty_plots/anscombe.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5])\n",
    "y1 = np.array([8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68])\n",
    "y2 = np.array([9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74])\n",
    "y3 = np.array([7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73])\n",
    "x4 = np.array([8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8])\n",
    "y4 = np.array([6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit(x):\n",
    "    return 3 + 0.5 * x\n",
    "\n",
    "xfit = np.array([np.min(x), np.max(x)])\n",
    "\n",
    "plt.subplot(221)\n",
    "plt.plot(x, y1, 'ks', xfit, fit(xfit), 'r-', lw=2)\n",
    "plt.axis([2, 20, 2, 14])\n",
    "plt.setp(plt.gca(), xticklabels=[], yticks=(4, 8, 12), xticks=(0, 10, 20))\n",
    "plt.text(3, 12, 'I', fontsize=20)\n",
    "\n",
    "plt.subplot(222)\n",
    "plt.plot(x, y2, 'ks', xfit, fit(xfit), 'r-', lw=2)\n",
    "plt.axis([2, 20, 2, 14])\n",
    "plt.setp(plt.gca(), xticks=(0, 10, 20), xticklabels=[],\n",
    "         yticks=(4, 8, 12), yticklabels=[], )\n",
    "plt.text(3, 12, 'II', fontsize=20)\n",
    "\n",
    "plt.subplot(223)\n",
    "plt.plot(x, y3, 'ks', xfit, fit(xfit), 'r-', lw=2)\n",
    "plt.axis([2, 20, 2, 14])\n",
    "plt.text(3, 12, 'III', fontsize=20)\n",
    "plt.setp(plt.gca(), yticks=(4, 8, 12), xticks=(0, 10, 20))\n",
    "\n",
    "plt.subplot(224)\n",
    "xfit = np.array([np.min(x4), np.max(x4)])\n",
    "plt.plot(x4, y4, 'ks', xfit, fit(xfit), 'r-', lw=2)\n",
    "plt.axis([2, 20, 2, 14])\n",
    "plt.setp(plt.gca(), yticklabels=[], yticks=(4, 8, 12), xticks=(0, 10, 20))\n",
    "plt.text(3, 12, 'IV', fontsize=20)\n",
    "\n",
    "# verify the stats\n",
    "pairs = (x, y1), (x, y2), (x, y3), (x4, y4)\n",
    "for x, y in pairs:\n",
    "    print('mean=%1.2f, std=%1.2f, r=%1.2f' % (np.mean(y), np.std(y),\n",
    "          np.corrcoef(x, y)[0][1]))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlation != Causation\n",
    "\n",
    "Measures a relationship, but neither correlation nor linear regression can say one variable _causes_ the other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](img/maine_margarine.png)\n",
    "\n",
    "Source: http://www.tylervigen.com/spurious-correlations "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](img/pool_cage.png)\n",
    "\n",
    "Source: http://www.tylervigen.com/spurious-correlations \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Difference in means: t-test\n",
    "\n",
    "The t-test tests for a difference in means between groups."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Student's t-test assumptions:\n",
    "\n",
    "- Data are continuous\n",
    "- Normally distributed\n",
    "- Large enough sample size\n",
    "- Equal variance b/w groups\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$H_0: \\bar x = \\bar y$\n",
    "\n",
    "$H_a: \\bar x \\ne \\bar y$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.groupby(['Violent'])['TeenBrth'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$H_0: \\bar x_{ViolentTeenBirthRate} = \\bar y_{NotViolentTeenBirthRate}$\n",
    "\n",
    "$H_a: \\bar x_{ViolentTeenBirthRate} \\ne \\bar y_{NotViolentTeenBirthRate}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "group1 = df.where(df.Violent==0).dropna()['TeenBrth']\n",
    "group2 = df.where(df.Violent==1).dropna()['TeenBrth']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ttest_ind(group1, group2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The p-value here is $1.08\\times 10^{-7}$ ...which means that observing this difference in means due to random chance alone is very unlikely. \n",
    "\n",
    "We reject the null in favor of the alternative hypothesis, concluding that the means between the groups are not the same and that violence plays a role on Teen Birth Rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clicker Question #6\n",
    "\n",
    "To understand the relationship between height and gender in this classroom, what would you use?\n",
    "\n",
    "- A) linear regression\n",
    "- B) correlation\n",
    "- C) t-test\n",
    "- D) something else"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differences in means by ANOVA\n",
    "\n",
    "Parametric test assumptions\n",
    "\n",
    "- Population distributions are normal\n",
    "- Samples have equal variances\n",
    "- Independence\n",
    "\n",
    "Hypothesis\n",
    "\n",
    "$H_0: \\bar x_1 = \\bar x_2 = \\bar x_3 \\ldots $\n",
    "\n",
    "$H_a: $ at least one group mean is different \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regions = pd.read_csv('https://raw.githubusercontent.com/cphalpert/census-regions/master/us%20census%20bureau%20regions%20and%20divisions.csv')\n",
    "regions = regions.rename({'State':'Location'},axis=1).set_index('Location')\n",
    "regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.set_index('Location').join(regions['Region'])\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So how would you run an ANOVA checking for a difference in TeenBrth across regions?\n",
    "\n",
    "this is left as an exercise to the reader. But there are two major options for you.  If you want to manually group the categorical variable ($\\bar x_1 = \\bar x_2 = \\bar x_3$) you can use scipy:\n",
    "https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.f_oneway.html\n",
    "\n",
    "If you want to build a model that is more precise or custom, with automatic grouping by a variable inside the dataframe then you can use statsmodels:\n",
    "https://www.statsmodels.org/dev/anova.html"
   ]
  },
  {
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}