ANOVA_implementation.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Notebook for implementation of ANOVA one-way tests in pure Python\n",
    "##### Brian Henn, Insight Data Science, August 2018\n",
    "\n",
    "References: \n",
    "https://en.wikipedia.org/wiki/F-test#One-way_ANOVA_example\n",
    "https://www.marsja.se/four-ways-to-conduct-one-way-anovas-using-python/\n",
    "http://www.stat.purdue.edu/~ghobbs/STAT_512/Lecture_Notes/ANOVA/Topic_20.pdf\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import stuff\n",
    "\n",
    "#import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(30, 3)\n",
      "group\n",
      "ctrl    10\n",
      "trt1    10\n",
      "trt2    10\n",
      "Name: weight, dtype: int64\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5,0.98,'')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAEaCAYAAAA8ObDuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHUpJREFUeJzt3XuYXXV97/H3hwRKNCheMIdWUXuzQSxapj5VcsoElKOVPtpqq2m1pxKb9hxL7UULFam1NRVqH3s5bY8nNrb0YlApaAs2ajEDpooablYI+FgviKBIFSSKmMD3/LHWwGaYyexZM3v2nsz79Tz72bPX9bv3ZGU+s+a7fitVhSRJkqS5OWjYBUiSJElLkUFakiRJ6sAgLUmSJHVgkJYkSZI6MEhLkiRJHRikJUmSpA4M0pIOeElWJqkk/7YA29qZZN9C1DVK+5IkzZ1BWtJAJTmuDbGXzzD/59r5leSJ08xfleTbSb6V5LsGX/HCWsgQPwxJHprktCTvT3JLkruT7Eny6STvSPLiJAcPu05JGoaVwy5A0gHvKuDrwFiSh1XVN6bMPxEoIO3XW6fMPx74LuCDVXV3lwKqal+StcA3u6y/XCVZB7wT+G7gy8AO4EZgBfB44CRgA3A58IwhlSlJQ2OQljRQVXVvkgngp4ATgH+ZssiJwATww0wfpE9sny+ZZx3Xz2f95SbJDwPbgVXAGcBbqmrvlGVWAj8NbFr8CiVp+GztkLQYJkPwib0TkzwBeGI7/1Jg/TTrzhikk/x8kokkt7ftH9cleW2SQ6YsN2N7RZLvTnJukq8muSvJVUlemuRZ7Tqvm+4NJTk4yeuSfKZtd7gxyZt62xySvAKYDJ8n9bSwzLjdGfZ1aJI/TPL5dl+fSXJW7/tM8qi2/huSZIbtbG/3fWwfu/0L4KHAG6vqnKkhGpoz/VX1LuA5U/Zz32eX5MeSvC/J19ppj+1Z7keTXNh+9ne37+8vkvy3aWqfsV88ySvabb90yvSb2s/q8CR/leTm9t/JtUleOdPnJEn98oy0pMXwofb5pCnTT+qZfwfw00mOrqrrAJI8DBgDbgeu7F0xybnAL9C0Gpzfrv9MYDNwYpL/UVX37K+oNrB9FDiK5qz45cCRwBbg/bO8p/No2hm2A3cCz6M5c/to4JfaZa4E/gA4C/gc8Hc96182y/Z7/RPw1PZ5H/AC4PeB49qvqar/SvJu4GXAOE0bRu97fTzwbOBjVXXN/naW5EnAf6dphfnj2YqrqpkuiFwH/C7Ne90KPIb2F4skLwDeRdPWcz7N9/FHgVcCz09yfFXdONu++/BdNP++VgPvaF//DM0vCj8IvGoB9iFpuaoqHz58+Bj4A/gScC9wRM+0f6QJoSuBJ9OEql/tmf+T7bQLpmzrFe30dwGHTpn3B+28V/ZMW9lO+7cpy57bTt88ZfqPAN9p571uyryd7fSPA4/omb4a+CxN0D1itn33+ZlN7ut64PCe6ava/RewoWf6j7XTzptmW29s5728j/2e2i67o+P3+lnt+gVsnGb+w2j65vcBz5wy78x2vfdN81nsm2F/k/8eXjpl+k3t9EuBQ3qmP5rmF5uaun8fPnz4mMvD1g5Ji2UHzQWFve0b64EPV9MicC1wKw9s/5ipreNVNEH3FVX17Snz3kBzBvvn91dMkkOBF9MEuj/snVdVV9KE/P357ar6es86e4BtNBfiHTfLunP1hqq6vWdfdwGvbV+e2jP9cpqLO38qyRGT09te5lNpztq/s4/9TbZWfGm6mUl+M8nvTXkcNc2iu6pqas87NP3yhwPvqKqPTJn3RzRnp5+b5Hv6qLUfZ1TVdyZfVNVtNH+5AHj5Au1D0jJka4ekxXIJTbg9EXhXO4rGkcCf9CwzATw7yUFVdS/3B+n7epuTHAYcA3wF+M0Z2ly/DaydpZ61NH/m/0hVTTeax07gF/ez/q5ppn2xfX7ELPueq0unmXYZzRn+p02Z/n9pWlNeThNKoTmzfyTwf6rqW33sb/JDrRnm/yYwNeT+G00A7vXxGdb/kfb5Q1NnVNXeJDuBn6NpZ5k2zM/Bd2hadqaaaJ+nfn6S1DeDtKTFMnlW+aQpz71hagL4WeBpSb4APAX4UlXd0LPMI9vnNcDr97O/2W5k8vD2+SszzJ9pOsA97Rnomfa5YpZ9z9WtUydU1XeSfJ3738ekfwTeDGxK8uaqKuCX23lb+tzfLe3ztGeEq6r3gsHzaM7sT+fLM0yfrPmWGeZPTj98PzX269b2M5hqsrapn58k9c3WDkmLopoLx/4T+P4kj6M523w7TSvCpMkL5E5sH+HBbR13tM+fqKrs5zHbTUImx7NeM8P8maYPw2OmTmhH7HgE938eALRnnP8O+D6akUImLzL896r6VJ/7+/f2+elJVneueuYz2pM1P2h0jtaRU5aD5ux7kkz3c2t/gfsxM4zOMbnvO6aZJ0l9MUhLWkyTofhZNGNKX9q2cAD3jfV8C/cH6d51Jpe5HbgBeEqS+ZyxvA64G3hqkodOM3/dPLbda/L9zecs9QnTTPtxmv/Dr5pm3l+1z79MM4LIQcD/63dn7V8APkwz/N1vzanS/kzWPD51Rjt84PFTloOml/0gpj9LPraffR1CcxHmVJP7nu7zk6S+GKQlLabJNo7foGnR2DHNMhM0Q6+d3L6e7kYsbwEOBbYmedCf5pM8Msl+e1/bixTfTXNW97W989p193uxYr/aXxS+TjPEXle/2/tLQ5JV3H+B5N9Ms8/raT7H59PcLOVrNO91Lk6jGf7urCSvzjS3AW/PDh82x+0CXEDz14iXJvnRKfN+i+auidurqrc/erLf+pd6F05yMs1wdvtz9pQxtx9NMzoITPP5SVK/7JGWtJg+RPPn/qf0vJ5qB81tp58I3DAlTAFQVVuSHEcTEk9I8gGaC90eCXwvTRB/G/Crs9Tz2zRnJl+b5Jk0Y0ofSdPzezHNGM33zrh2/y4BXpTkvTRnQPcBE1W1s49176EZVu9TSXrHkf5e4L00I4VM569o3tsRwJ9MM7rJflXVNUmeSzPKx5tpLuzcQXNB5Qqa1ogTaW4f/lmaoeb63fY3kmxst/3hdvzrL9KcWX42cDPwv6astpUmZJ/V/qKzG/ghmpvBXAi8cIbd3UQT9j+V5J9pLjB9UVv/n08zaogk9c0gLWnRVNVXk/wHze3AbwOm69ntPUs9423Bq+qXk1xM077wbJo+2f8CvkAzWsU/9FHPLUmeQXN29ydoWgCub7e5lyawfmPmLfTtNJoAfBJwCs1fA8+iGRmkHy+kubByA03Q/1L7+uwZLqSDJlx+neaMe78XGT5AVX04yQ/SDJ33kzQtOY+g+Wxupmn/uAB4T+/wcn1u+4Ik64DfAZ5LM7b0LTSjjvxBVd0yZfkvJzmB5ns7TjN04idoPtMnMXOQvpsm8L+JZiSQR9H06m8G/nIuNUvSVJn5/2BJWr6SnENzxvpZVTVjoB9VSX6Appf80qqa7tbrB7wkNwHfrqrvH3Ytkg5M9khLWtaSfPc0046luVX1bfR/1njUvIZm1JO/GHYhknSg6qu1I8lvcP8tWP+D5hazc+q3k6QRdXWS3TRtJt8CfpCmzeMgmttb3z3M4uaiHepuA02rw/8ErqRp8ZAkDcCsZ6TbW7T+GjBWVcfQXGTykkEXJkmL5K00N+X4OeDXgWcC24H1VdXP7bRHyQ/Q9AL/DPB+4IW9wwtKkhbWrD3SbZC+HDiW5qKb99Bc6fyBwZcnSZIkjaZZz0i3Q0/9Mc3QUrcAdxiiJUmStNzN2iOd5BE0g/o/kWYA/XcneWlV/cOU5TbRjOnKqlWrjnvc4x43gHK12O69914OOshrUqVR4nEpjSaPzQPHpz/96duq6ojZluvnYsNnAZ+rqq8CJLmApofwAUG6qrbQjlU6NjZWu3btmnPRGj0TExOMj48PuwxJPTwupdHksXngSPKFfpbr59emG4EfS/KQJKEZ/H73fIqTJEmSlrp+eqQ/BpxPM4zSf7TrdLpLliRJknSg6Gsc6ap6Pc3taCVJkiThnQ0lSZKkTgzSkiRJUgcGaUmSJKkDg7QkSZLUgUFakiRJ6sAgLUmSJHVgkJYkSZI6MEhLkiRJHRikJUmSpA4M0pIkSVIHBmlJkiSpA4O0JEmS1IFBWpIkSerAIC1JkiR1YJCWJEmSOjBIS5IkSR0YpCVJkqQODNKSJElSB7MG6SRPSnJ1z+MbSX59MYqTJEmSRtXK2RaoqhuApwIkWQF8CbhwwHVJkiRJI22urR0nAf9ZVV8YRDGSJEnSUjHXIP0SYNsgCpEkSZKWklRVfwsmhwA3A0+uqq9MM38TsAlgzZo1x5133nkLWaeGZM+ePaxevXrYZUjq4XEpjSaPzQPH+vXrr6iqsdmWm0uQfj7wyqo6ebZlx8bGateuXX1tV6NtYmKC8fHxYZchqYfHpTSaPDYPHEn6CtJzae3YgG0dkiRJEtBnkE7yEODZwAWDLUeSJElaGmYd/g6gqr4FPGrAtUiSJElLhnc2lCRJkjowSEuSJEkdGKQlSZKkDgzSkiRJUgcGaUmSJKkDg7QkSZLUgUFakiRJ6sAgLUmSJHVgkJYkSZI6MEhLkiRJHRikJUmSpA4M0pIkSVIHBmlJkiSpg5XDLkCSJGmYkgy7hAeoqmGXoD55RlqSJC1rVbUgj8efftGCbEdLh0FakiRJ6sAgLUmSJHVgkJYkSZI66CtIJzk8yflJrk+yO8kzBl2YJOmBtm3bxjHHHMNJJ53EMcccw7Zt24ZdkiQta/2O2vFnwPaqelGSQ4CHDLAmLQCvQJYOLNu2bePMM89k69at3HPPPaxYsYKNGzcCsGHDhiFXJ0nL06xnpJM8DPhxYCtAVX2nqm4fdGGaH69Alg4smzdvZuvWraxfv56VK1eyfv16tm7dyubNm4ddmiQtW/20dnwv8FXgb5JcleSvkzx0wHVJknrs3r2bdevWPWDaunXr2L1795AqkiT109qxEvgR4LSq+liSPwPOAM7qXSjJJmATwJo1a5iYmFjgUjUsfi+l4TvqqKM49dRT2blzJzfeeCNHHXUU69at46ijjvIYlUaIx+Py0k+Qvgm4qao+1r4+nyZIP0BVbQG2AIyNjdX4+PhC1ahh2n4xfi+l4TvllFN461vfyjnnnMPRRx/Nddddx+mnn86v/MqveIxKo8KfmcvOrEG6qr6c5ItJnlRVNwAnAdcNvjRJ0qQdO3Zw+umn8/a3v53du3ezdu1aTj/9dN7znvcMuzRJWrb6HbXjNOAf2xE7Pgu8fHAlSZKm2r17N1dddRVvfOMbmZiYYHx8nL179/KmN71p2KVJ0rLV1zjSVXV1VY1V1Q9X1Quq6uuDLkySdL+1a9eyc+fOB0zbuXMna9euHVJFkiTvbChJS8CZZ57Jxo0b2bFjB/v27WPHjh1s3LiRM888c9ilSdKy1W9rhyRpiCZvunLaaafd1yO9efNmb8YiSUNkkJakJWLDhg1s2LDhvh5pSXDsGz7AHXftHXYZ93nCGRcPuwQevupgrnn9ycMuY1kwSEuSpCXrjrv28vmznzfsMgBG5pfcUQjzy4U90pIkSVIHBmlJkiSpA4O0JEmS1IFBWpIkSerAIC1JkiR1YJCWJEmSOnD4O0mStGQdtvYMnnLuGcMu437nDrsAOGwtwGgMCXigM0hLkqQl687dZzuO9BSOI714bO2QJEmSOjBIS5IkSR3Y2jGCjn3DB7jjrr3DLuM+o/AnooevOphrXn/ysMuQJEm6j0F6BN1x1177vaYYhTAvSRpNI/UzYvvwa3n4qoOHXcKyYZCWJElL1qiceIIm0I9SPRo8e6QlSZKkDvo6I53k88CdwD3AvqoaG2RRy51jYj6YY2JKkqRRM5fWjvVVddvAKtF9HBPzwUaq/02SJAlbOyRJkqRO+g3SBXwgyRVJNg2yIEmSJGkp6Le14/iqujnJY4APJrm+qi7rXaAN2JsA1qxZw8TExMJWusyMyue3Z8+ekallVOqQhm2UjktJD+Sxubz0FaSr6ub2+dYkFwJPBy6bsswWYAvA2NhYjUJf7ZK1/eKR6EuG0emRHqXPRBq2kTkuJT2QP6uWnVmDdJKHAgdV1Z3t1ycDvz/wyiRJkhZBkoXb1jnz30ZVzX8jWhT9nJFeA1zY/iNbCbyjqrYPtCpJGkHHvuED3HHX3nlt4wvnnLJA1czf40+/aN7bePiqg7nm9ScvQDXS8CxUcPWvRcvPrEG6qj4LHLsItUjSSLvjrr3zH5ry7Pn/wB6lH9YOTSlpOXP4O0mSJKkDg7QkSZLUgUFakiRJ6sAgLUmSJHVgkJYkSZI6MEhLkiRJHRikJUmSpA76ukW4Ft9Ijc26ffi1PHzVwcMuQeKwtWfwlHPPGHYZjXOHXUDjsLUA8xxbW5KWKIP0CJr3DR8W0BPOuHik6pGG6c7dZ4/E8eANWSRpNNjaIUmSJHVgkJYkSZI6MEhLkiRJHRikJUmSpA4M0pIkSVIHBmlJkiSpA4O0JEmS1IHjSEvSHIzMuMkjcKMk8GZJkpY3g/QBKsnCbeuc+W+jqua/EWnIRuFmLOCNkiRpVPTd2pFkRZKrklw0yIK0MKpqQR47duxYkO1IkiQdaObSI/0qYPegCpEkSZKWkr6CdJLHAs8D/nqw5UiSJElLQ79npP8U+G3g3gHWIkmSJC0Zs15smOQU4NaquiLJ+H6W2wRsAlizZg0TExMLVaOGaM+ePX4vpRHkcSmNHn9mLj+Z7UKwJG8CXgbsAw4FHgZcUFUvnWmdsbGx2rVr10LWqSGZmJhgfHx82GVI6uGoHdJo8mfmgSPJFVU1Nttys7Z2VNXvVNVjq+oJwEuAD+0vREuSJEnLgXc2lCRJkjqY0w1ZqmoCmBhIJZIkSdIS4hlpSZIkqQODtCRJktSBQVqSJEnqwCAtSZIkdWCQliRJkjowSEuSJEkdzGn4O0nS/CRZmO2cM/9tzHZnW0nS/nlGWpIWUVXN+7Fjx44F2Y4kaX4M0pIkSVIHBmlJkiSpA4O0JEmS1IFBWpIkSerAIC1JkiR1YJCWJEmSOjBIS5IkSR0YpCVJkqQODNKSJElSBwZpSZIkqYNZg3SSQ5N8PMk1Sa5N8obFKEySJEkaZSv7WOZu4MSq2pPkYGBnkn+tqssHXJskSZI0smYN0lVVwJ725cHtowZZlCRJkjTq+uqRTrIiydXArcAHq+pjgy1LkiRJGm1pTjj3uXByOHAhcFpVfWrKvE3AJoA1a9Ycd9555y1knRqSPXv2sHr16mGXIamHx6U0mjw2Dxzr16+/oqrGZltuTkEaIMnrgW9W1R/PtMzY2Fjt2rVrTtvVaJqYmGB8fHzYZUjq4XEpjSaPzQNHkr6CdD+jdhzRnokmySrgWcD18y9RkiRJWrr6GbXjSODcJCtogve7quqiwZYlSZIkjbZ+Ru34JPC0RahFkiRJWjK8s6EkSZLUgUFakiRJ6sAgLUmSJHVgkJYkSZI6MEhLkiRJHRikJUmSpA4M0pIkSVIHBmlJkiSpA4O0JEmS1IFBWpIkSerAIC1JkiR1YJCWJEmSOjBIS5IkSR0YpCVJkqQODNKSJElSBwZpSZIkqQODtCRJktSBQVqSJEnqYNYgneRxSXYk2Z3k2iSvWozCJEmSpFG2so9l9gG/VVVXJjkMuCLJB6vqugHXJkmSJI2sWc9IV9UtVXVl+/WdwG7gewZdmCRJkjTK5tQjneQJwNOAjw2iGEmSJGmp6Ke1A4Akq4F/An69qr4xzfxNwCaANWvWMDExsVA1aoj27Nnj91IaMR6X0mjy2Fx+UlWzL5QcDFwEvL+q3jLb8mNjY7Vr164FKE/DNjExwfj4+LDLkNTD41IaTR6bB44kV1TV2GzL9TNqR4CtwO5+QrQkSZK0HPTTI3088DLgxCRXt4+fGHBdkiRJ0kibtUe6qnYCWYRaJEmSpCXDOxtKkiRJHRikJUmSpA4M0pIkSVIHBmlJkiSpA4O0JEmS1IFBWpIkSerAIC1JkiR1YJCWJEmSOjBIS5IkSR0YpCVJkqQODNKSJElSBwZpSZIkqQODtCRJktSBQVqSJEnqwCAtSZIkdWCQliRJkjowSEuSJEkdGKQlSZKkDmYN0knenuTWJJ9ajIIkSZKkpaCfM9J/CzxnwHVIkiRJS8qsQbqqLgO+tgi1SJIkSUuGPdKSJElSBysXakNJNgGbANasWcPExMRCbVpDtGfPHr+X0ojxuJRGk8fm8rNgQbqqtgBbAMbGxmp8fHyhNq0hmpiYwO+lNFo8LqXR5LG5/NjaIUmSJHXQz/B324CPAk9KclOSjYMvS5IkSRpts7Z2VNWGxShEkiRJWkps7ZAkSZI6MEhLkiRJHRikJUmSpA4M0pIkSVIHBmlJkiSpA4O0JEmS1IFBWpIkSerAIC1JkiR1YJCWJEmSOjBIS5IkSR0YpCVJkqQODNKSJElSBwZpSZIkqQODtCRJktSBQVqSJEnqwCAtSZIkdWCQliRJkjowSEuSJEkd9BWkkzwnyQ1JPpPkjEEXJUmSJI26WYN0khXAXwLPBY4GNiQ5etCFSZIkSaOsnzPSTwc+U1WfrarvAOcBzx9sWZIkSdJo6ydIfw/wxZ7XN7XTJEmSpGVrZR/LZJpp9aCFkk3ApvblniQ3zKcwjYxHA7cNuwhJD+BxKY0mj80Dx+P7WaifIH0T8Lie148Fbp66UFVtAbb0VZqWjCS7qmps2HVIup/HpTSaPDaXn35aOz4B/ECSJyY5BHgJ8M+DLUuSJEkabbOeka6qfUl+FXg/sAJ4e1VdO/DKJEmSpBHWT2sHVfU+4H0DrkWjyXYdafR4XEqjyWNzmUnVg64blCRJkjQLbxEuSZIkdWCQ1oMkecH+7l6Z5PeSvHoxa5KWiySHJ/nf/c5Lsj3J7UkuWpwKpeVpLsdmkqcm+WiSa5N8MsmLF69SLSaDtKbzAprbwT9Ikr766iV1djjwoB/WSVbMMO/NwMsWoS5puZvLsfkt4Beq6snAc4A/TXL4olSpRWUoEkl+AXg1zY12bgWOA05I8jrghcBW4CPA8Tj0oTRoZwPfl+RqYC+wB7gFeCrwyZ55H6yq11TVJUnGh1attHzM6dicXKmqbk5yK3AEcPvil61BMkgvc0meDJwJHF9VtyV5JPAW4KKqOr9dBuDwqjqhff17QypXWg7OAI6pqqe2Afni9vXnkjxhct4Q65OWq07HZpKnA4cA/7mItWqR2NqhE4Hzq+o2gKr62gzLvXPxSpLU4+NV9blhFyHpQWY9NpMcCfw98PKqundxytJiMkgrNC0ds/nmoAuRNC2PPWk07ffYTPIwmrPWr6uqyxenJC02g7QuAX42yaMA2taOO4HDhlqVtHzt7/jz2JSGp+9jM8khwIXA31XVuxehNg2JPdLLXFVdm2QzcGmSe4CrgLcBb0vya8CLhlqgtMxU1X8l+fcknwLuAr4yw7x/rarXJPkw8EPA6iQ3ARur6v3DqV46cM3l2ASuAX4ceFSSX2wX+8Wqunqx69ZgeWdDSZIkqQNbOyRJkqQODNKSJElSBwZpSZIkqQODtCRJktSBQVqSJEnqwCAtSZIkdWCQlqQRlsTx/iVpRBmkJWmIkpyV5PokH0yyLcmrk0wk+cMklwKvSvL4JJck+WT7fFS77t8meVHPtva0z+NJLktyYZLrkrw1if/fS9IC80yHJA1JkjHghcDTaP4/vhK4op19eFWd0C73LzS3Gj43yanAnwMvmGXzTweOBr4AbAd+Gjh/wd+EJC1jnqGQpOFZB7y3qu6qqjuBf+mZ986er58BvKP9+u/b9Wbz8ar6bFXdA2zrcx1J0hwYpCVpeLKfed/cz7xqn/fR/j+eJMAh0ywz02tJ0jwZpCVpeHYCP5nk0CSrgefNsNxHgJe0X/98ux7A54Hj2q+fDxzcs87Tkzyx7Y1+cc86kqQFYo+0JA1JVX0iyT8D19D0Mu8C7phm0V8D3p7kNcBXgZe3098GvDfJx4FLeOBZ7I8CZwNPAS4DLhzIm5CkZSxV/rVPkoYlyeqq2pPkITSBd1NVXTnPbY4Dr66qUxaiRknS9DwjLUnDtSXJ0cChwLnzDdGSpMXjGWlJkiSpAy82lCRJkjowSEuSJEkdGKQlSZKkDgzSkiRJUgcGaUmSJKkDg7QkSZLUwf8HwRq+WPE62fUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# load test dataset from web\n",
    "data_url = 'https://vincentarelbundock.github.io/Rdatasets/csv/datasets/PlantGrowth.csv'\n",
    "plants = pd.read_csv(data_url)\n",
    "\n",
    "# print out some basic info about this dataframe\n",
    "print(plants.shape)\n",
    "print(plants.groupby('group')['weight'].count())\n",
    "\n",
    "# make a boxplot by the group\n",
    "fig = plt.figure()\n",
    "fig.set_size_inches([12,4])\n",
    "ax = fig.subplots(1)\n",
    "plants.boxplot('weight', by='group', ax=ax)\n",
    "ax.set_ylim([0,8])\n",
    "ax.set_title('Weight by Group', fontsize=20)\n",
    "plt.suptitle('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10 10 10]\n"
     ]
    }
   ],
   "source": [
    "grps = pd.unique(plants.group.values)\n",
    "d_data = {grp:plants['weight'][plants.group == grp] for grp in grps}\n",
    " \n",
    "k = len(pd.unique(plants.group))  # number of conditions\n",
    "N = len(plants.values)  # conditions times participants\n",
    "n = plants.groupby('group').size().values #Participants in each condition\n",
    "\n",
    "print(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.492090000000001\n",
      "0.38859592592592596\n",
      "3.766340000000002\n",
      "1.883170000000001\n"
     ]
    }
   ],
   "source": [
    "# now let's manually compute the ANOVA parameters to conduct a one-way test between the group and weight\n",
    "\n",
    "# first, sum of squared errors of the model (SSE)\n",
    "# i.e., the sum of squared differences between each point and its group mean\n",
    "\n",
    "for i, group in enumerate(plants['group'].values):\n",
    "    #print(group, i)\n",
    "    plants.loc[i, 'group_mean'] = plants[plants['group'] == group].weight.values.mean()\n",
    "    \n",
    "SSE = 0\n",
    "for grp in grps:\n",
    "    SSE = SSE + sum(\n",
    "        (plants[plants['group'] == grp].group_mean.values - plants[plants['group'] == grp].weight.values)**2)\n",
    "\n",
    "print(SSE)\n",
    "\n",
    "# now compute MSE\n",
    "df_model = N - k\n",
    "MSE = SSE/df_model\n",
    "\n",
    "print(MSE)\n",
    "\n",
    "# next, let's compute the treatment squared errors (differences between groups)\n",
    "SSTr = 0\n",
    "y_barbar = plants.weight.values.mean()\n",
    "for i, grp in enumerate(grps):\n",
    "    SSTr = SSTr + n[i]*(y_barbar - plants[plants['group'] == grp].group_mean.values[0])**2\n",
    "    \n",
    "#print(y_barbar)\n",
    "print(SSTr)\n",
    "\n",
    "\n",
    "#SSbetween = (sum(plants.groupby('group').sum()['weight']**2)/n[0]) - (plants['weight'].sum()**2)/N\n",
    "#print(SSbetween)\n",
    "\n",
    "# now compute MSTr\n",
    "df_groups = k - 1\n",
    "MSTr = SSTr/df_groups\n",
    "print(MSTr)\n",
    "\n",
    "#SSbetween = (sum(plants.groupby('group').sum()['weight']**2)/n[0]) - (plants['weight'].sum()**2)/N\n",
    "#print(SSbetween)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.846087862380138\n",
      "0.01590995832562281\n"
     ]
    }
   ],
   "source": [
    "# now compute the F-value and the p-value using the scipy F-distribution\n",
    "\n",
    "F = MSTr/MSE\n",
    "print(F)\n",
    "\n",
    "p = 1 - stats.f.cdf(F, df_groups, df_model)\n",
    "print(p)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:insight]",
   "language": "python",
   "name": "conda-env-insight-py"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}