Python_ANOVA/One_Way_Python_ANOVA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "\n",
    "datafile = 'https://vincentarelbundock.github.io/Rdatasets/csv/datasets/PlantGrowth.csv'\n",
    "df = pd.read_csv(datafile)\n",
    "\n",
    "#Create a boxplot\n",
    "df.boxplot('weight', by='group', figsize=(12, 8))\n",
    "\n",
    "ctrl = df['weight'][df.group == 'ctrl']\n",
    "\n",
    "grps = pd.unique(df.group.values)\n",
    "d_data = {grp:df['weight'][df.group == grp] for grp in grps}\n",
    "\n",
    "k = len(pd.unique(df.group))  # number of conditions\n",
    "N = len(df.values)  # conditions times participants\n",
    "n = df.groupby('group').size()[0] #Participants in each condition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Anova usig Scipy\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.846087862380136 0.0159099583256229\n"
     ]
    }
   ],
   "source": [
    "from scipy import stats\n",
    "\n",
    "F, p = stats.f_oneway(d_data['ctrl'], d_data['trt1'], d_data['trt2'])\n",
    "\n",
    "print(F, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Degrees of Fredoom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Calculate the degrees of freedom\n",
    "DFbetween = k - 1\n",
    "DFwithin = N - k\n",
    "DFtotal = N - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sum of Squares Between, WIthin, and Total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Between\n",
    "SSbetween = (sum(df.groupby('group').sum()['weight']**2)/n) \\\n",
    "    - (df['weight'].sum()**2)/N\n",
    "\n",
    "# Within\n",
    "sum_y_squared = sum([value**2 for value in df['weight'].values])\n",
    "SSwithin = sum_y_squared - sum(df.groupby('group').sum()['weight']**2)/n\n",
    "\n",
    "# Total\n",
    "SStotal = sum_y_squared - (df['weight'].sum()**2)/N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Means Square Errors and F-value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "MSbetween = SSbetween/DFbetween\n",
    "MSwithin = SSwithin/DFwithin\n",
    "\n",
    "F = MSbetween/MSwithin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Obtaining the p-value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = stats.f.sf(F, DFbetween, DFwithin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate Effect sizes: Eta Squared and Omega Squared:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sum_sq</th>\n",
       "      <th>df</th>\n",
       "      <th>F</th>\n",
       "      <th>PR(&gt;F)</th>\n",
       "      <th>EtaSq</th>\n",
       "      <th>OmegaSq</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Group</th>\n",
       "      <td>3.76634</td>\n",
       "      <td>2</td>\n",
       "      <td>4.84609</td>\n",
       "      <td>0.01591</td>\n",
       "      <td>0.264148</td>\n",
       "      <td>0.204079</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Residual</th>\n",
       "      <td>10.49209</td>\n",
       "      <td>27</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            sum_sq  df        F   PR(>F)     EtaSq   OmegaSq\n",
       "Group      3.76634   2  4.84609  0.01591  0.264148  0.204079\n",
       "Residual  10.49209  27                                      "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eta_sqrd = SSbetween/SStotal\n",
    "om_sqrd = (SSbetween - (DFbetween * MSwithin))/(SStotal + MSwithin)\n",
    "\n",
    "results = {'sum_sq':[SSbetween, SSwithin],\n",
    "           'df':[DFbetween, DFwithin],\n",
    "           'F':[F, ''],\n",
    "            'PR(>F)':[p, ''],\n",
    "          'EtaSq':[eta_sqrd, ''],\n",
    "          'OmegaSq':[om_sqrd, '']}\n",
    "\n",
    "aov_table1 = pd.DataFrame(results, \n",
    "                          index=['Group', 'Residual'])\n",
    "    \n",
    "aov_table1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statsmodels One-Way ANOVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            sum_sq    df         F   PR(>F)\n",
      "group      3.76634   2.0  4.846088  0.01591\n",
      "Residual  10.49209  27.0       NaN      NaN\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "from statsmodels.formula.api import ols\n",
    "\n",
    "mod = ols('weight ~ group',\n",
    "                data=df).fit()\n",
    "                \n",
    "aov_table = sm.stats.anova_lm(mod, typ=2)\n",
    "print(aov_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculating Eta-Squared"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            sum_sq    df         F   PR(>F)     EtaSq\n",
      "group      3.76634   2.0  4.846088  0.01591  0.264148\n",
      "Residual  10.49209  27.0       NaN      NaN       NaN\n"
     ]
    }
   ],
   "source": [
    "esq_sm = aov_table['sum_sq'][0]/(aov_table['sum_sq'][0] + aov_table['sum_sq'][1])\n",
    "\n",
    "aov_table['EtaSq'] = [esq_sm, 'NaN']\n",
    "print(aov_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pairwise Comparisons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>t</th>\n",
       "      <th>P&gt;|t|</th>\n",
       "      <th>Conf. Int. Low</th>\n",
       "      <th>Conf. Int. Upp.</th>\n",
       "      <th>pvalue-hs</th>\n",
       "      <th>reject-hs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>trt1-ctrl</th>\n",
       "      <td>-0.371</td>\n",
       "      <td>0.278782</td>\n",
       "      <td>-1.330791</td>\n",
       "      <td>0.194388</td>\n",
       "      <td>-0.943013</td>\n",
       "      <td>0.201013</td>\n",
       "      <td>0.194388</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>trt2-ctrl</th>\n",
       "      <td>0.494</td>\n",
       "      <td>0.278782</td>\n",
       "      <td>1.771996</td>\n",
       "      <td>0.087682</td>\n",
       "      <td>-0.078013</td>\n",
       "      <td>1.066013</td>\n",
       "      <td>0.167675</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>trt2-trt1</th>\n",
       "      <td>0.865</td>\n",
       "      <td>0.278782</td>\n",
       "      <td>3.102787</td>\n",
       "      <td>0.004459</td>\n",
       "      <td>0.292987</td>\n",
       "      <td>1.437013</td>\n",
       "      <td>0.013318</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            coef   std err         t     P>|t|  Conf. Int. Low  \\\n",
       "trt1-ctrl -0.371  0.278782 -1.330791  0.194388       -0.943013   \n",
       "trt2-ctrl  0.494  0.278782  1.771996  0.087682       -0.078013   \n",
       "trt2-trt1  0.865  0.278782  3.102787  0.004459        0.292987   \n",
       "\n",
       "           Conf. Int. Upp.  pvalue-hs  reject-hs  \n",
       "trt1-ctrl         0.201013   0.194388      False  \n",
       "trt2-ctrl         1.066013   0.167675      False  \n",
       "trt2-trt1         1.437013   0.013318       True  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pair_t = mod.t_test_pairwise('group')\n",
    "pair_t.result_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}