ANOVA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run for interactive plots with jupyter lab\n",
    "%matplotlib widget          \n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import statsmodels.api as sm\n",
    "\n",
    "import matplotlib\n",
    "# matplotlib.use('nbagg')    # Doesn't work with jupyter lab\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DescribeResult(nobs=30, minmax=(0.04995345894608716, 0.9648400471483856), mean=0.5558004906644833, variance=0.0770657116798012, skewness=-0.23634961375889468, kurtosis=-1.0681272762337057)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Creating artificial data for regression\n",
    "\n",
    "np.random.seed(1)       # For reproducibility\n",
    "X = np.random.random((30,2))\n",
    "b = [0.2,0.3]\n",
    "e = np.random.random(30)\n",
    "\n",
    "y = np.dot(X,b) + e\n",
    "\n",
    "stats.describe(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.182</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.121</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   3.002</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 22 May 2020</td> <th>  Prob (F-statistic):</th>  <td>0.0665</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>19:14:27</td>     <th>  Log-Likelihood:    </th> <td> -3.4473</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    30</td>      <th>  AIC:               </th> <td>   12.89</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    27</td>      <th>  BIC:               </th> <td>   17.10</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    0.5371</td> <td>    0.120</td> <td>    4.490</td> <td> 0.000</td> <td>    0.292</td> <td>    0.783</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>        <td>    0.2963</td> <td>    0.177</td> <td>    1.674</td> <td> 0.106</td> <td>   -0.067</td> <td>    0.659</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>        <td>    0.2649</td> <td>    0.192</td> <td>    1.379</td> <td> 0.179</td> <td>   -0.129</td> <td>    0.659</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 4.062</td> <th>  Durbin-Watson:     </th> <td>   2.074</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.131</td> <th>  Jarque-Bera (JB):  </th> <td>   1.944</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.310</td> <th>  Prob(JB):          </th> <td>   0.378</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 1.918</td> <th>  Cond. No.          </th> <td>    4.94</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.182\n",
       "Model:                            OLS   Adj. R-squared:                  0.121\n",
       "Method:                 Least Squares   F-statistic:                     3.002\n",
       "Date:                Fri, 22 May 2020   Prob (F-statistic):             0.0665\n",
       "Time:                        19:14:27   Log-Likelihood:                -3.4473\n",
       "No. Observations:                  30   AIC:                             12.89\n",
       "Df Residuals:                      27   BIC:                             17.10\n",
       "Df Model:                           2                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept      0.5371      0.120      4.490      0.000       0.292       0.783\n",
       "x1             0.2963      0.177      1.674      0.106      -0.067       0.659\n",
       "x2             0.2649      0.192      1.379      0.179      -0.129       0.659\n",
       "==============================================================================\n",
       "Omnibus:                        4.062   Durbin-Watson:                   2.074\n",
       "Prob(Omnibus):                  0.131   Jarque-Bera (JB):                1.944\n",
       "Skew:                          -0.310   Prob(JB):                        0.378\n",
       "Kurtosis:                       1.918   Cond. No.                         4.94\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "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",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>x1</th>\n",
       "      <td>0.229370</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.801840</td>\n",
       "      <td>0.105704</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x2</th>\n",
       "      <td>0.155665</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.901506</td>\n",
       "      <td>0.179227</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Residual</th>\n",
       "      <td>2.210333</td>\n",
       "      <td>27.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            sum_sq    df         F    PR(>F)\n",
       "x1        0.229370   1.0  2.801840  0.105704\n",
       "x2        0.155665   1.0  1.901506  0.179227\n",
       "Residual  2.210333  27.0       NaN       NaN"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#### Not very useful for general OLS. Needs work to explain with a better example ####\n",
    "\n",
    "# ANOVA\n",
    "# Need to use a different api to make the anova call\n",
    "\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "y2 = np.reshape(y,(np.size(y),1))                   # Convert np vector to np array\n",
    "y2X = np.concatenate((y2,X),axis=1)                # Stack the columns horizontally\n",
    "df = pd.DataFrame(y2X,columns=['y','x1','x2'])      # Create a data frame\n",
    "\n",
    "Reg1 = smf.ols(formula='y ~ x1 + x2', data=df)\n",
    "fit1 = Reg1.fit()\n",
    "fit1.summary()\n",
    "sm.stats.anova_lm(fit1,typ=2)         # two-way ANOVA table"
   ]
  }
 ],
 "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": 4
}