{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 8: Problems with the Error\n",
    "\n",
    "Load the packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy as sp\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generalized Least Squares\n",
    "\n",
    "Load the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>nhtemp</th>\n",
       "      <th>wusa</th>\n",
       "      <th>jasper</th>\n",
       "      <th>westgreen</th>\n",
       "      <th>chesapeake</th>\n",
       "      <th>tornetrask</th>\n",
       "      <th>urals</th>\n",
       "      <th>mongolia</th>\n",
       "      <th>tasman</th>\n",
       "      <th>year</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.66</td>\n",
       "      <td>-0.03</td>\n",
       "      <td>0.03</td>\n",
       "      <td>-0.66</td>\n",
       "      <td>0.33</td>\n",
       "      <td>-1.49</td>\n",
       "      <td>0.83</td>\n",
       "      <td>-0.12</td>\n",
       "      <td>1000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.63</td>\n",
       "      <td>-0.07</td>\n",
       "      <td>0.09</td>\n",
       "      <td>-0.67</td>\n",
       "      <td>0.21</td>\n",
       "      <td>-1.44</td>\n",
       "      <td>0.96</td>\n",
       "      <td>-0.17</td>\n",
       "      <td>1001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.60</td>\n",
       "      <td>-0.11</td>\n",
       "      <td>0.18</td>\n",
       "      <td>-0.67</td>\n",
       "      <td>0.13</td>\n",
       "      <td>-1.39</td>\n",
       "      <td>0.99</td>\n",
       "      <td>-0.22</td>\n",
       "      <td>1002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.55</td>\n",
       "      <td>-0.14</td>\n",
       "      <td>0.30</td>\n",
       "      <td>-0.68</td>\n",
       "      <td>0.08</td>\n",
       "      <td>-1.34</td>\n",
       "      <td>0.95</td>\n",
       "      <td>-0.26</td>\n",
       "      <td>1003</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.51</td>\n",
       "      <td>-0.15</td>\n",
       "      <td>0.41</td>\n",
       "      <td>-0.68</td>\n",
       "      <td>0.06</td>\n",
       "      <td>-1.30</td>\n",
       "      <td>0.87</td>\n",
       "      <td>-0.31</td>\n",
       "      <td>1004</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   nhtemp  wusa  jasper  westgreen  chesapeake  tornetrask  urals  mongolia  \\\n",
       "0     NaN -0.66   -0.03       0.03       -0.66        0.33  -1.49      0.83   \n",
       "1     NaN -0.63   -0.07       0.09       -0.67        0.21  -1.44      0.96   \n",
       "2     NaN -0.60   -0.11       0.18       -0.67        0.13  -1.39      0.99   \n",
       "3     NaN -0.55   -0.14       0.30       -0.68        0.08  -1.34      0.95   \n",
       "4     NaN -0.51   -0.15       0.41       -0.68        0.06  -1.30      0.87   \n",
       "\n",
       "   tasman  year  \n",
       "0   -0.12  1000  \n",
       "1   -0.17  1001  \n",
       "2   -0.22  1002  \n",
       "3   -0.26  1003  \n",
       "4   -0.31  1004  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "globwarm = pd.read_csv(\"data/globwarm.csv\")\n",
    "globwarm.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>nhtemp</td>      <th>  R-squared:         </th> <td>   0.476</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.446</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   15.47</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 25 Sep 2018</td> <th>  Prob (F-statistic):</th> <td>5.03e-16</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:55:20</td>     <th>  Log-Likelihood:    </th> <td>  50.995</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   145</td>      <th>  AIC:               </th> <td>  -83.99</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   136</td>      <th>  BIC:               </th> <td>  -57.20</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     8</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.2426</td> <td>    0.027</td> <td>   -8.980</td> <td> 0.000</td> <td>   -0.296</td> <td>   -0.189</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>wusa</th>       <td>    0.0774</td> <td>    0.043</td> <td>    1.803</td> <td> 0.074</td> <td>   -0.008</td> <td>    0.162</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>jasper</th>     <td>   -0.2288</td> <td>    0.078</td> <td>   -2.929</td> <td> 0.004</td> <td>   -0.383</td> <td>   -0.074</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>westgreen</th>  <td>    0.0096</td> <td>    0.042</td> <td>    0.229</td> <td> 0.819</td> <td>   -0.073</td> <td>    0.092</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>chesapeake</th> <td>   -0.0321</td> <td>    0.034</td> <td>   -0.943</td> <td> 0.347</td> <td>   -0.099</td> <td>    0.035</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>tornetrask</th> <td>    0.0927</td> <td>    0.045</td> <td>    2.057</td> <td> 0.042</td> <td>    0.004</td> <td>    0.182</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>urals</th>      <td>    0.1854</td> <td>    0.091</td> <td>    2.027</td> <td> 0.045</td> <td>    0.005</td> <td>    0.366</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>mongolia</th>   <td>    0.0420</td> <td>    0.046</td> <td>    0.917</td> <td> 0.361</td> <td>   -0.049</td> <td>    0.133</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>tasman</th>     <td>    0.1155</td> <td>    0.030</td> <td>    3.834</td> <td> 0.000</td> <td>    0.056</td> <td>    0.175</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>12.501</td> <th>  Durbin-Watson:     </th> <td>   0.817</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.002</td> <th>  Jarque-Bera (JB):  </th> <td>  17.388</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.490</td> <th>  Prob(JB):          </th> <td>0.000168</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 4.384</td> <th>  Cond. No.          </th> <td>    12.9</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:                 nhtemp   R-squared:                       0.476\n",
       "Model:                            OLS   Adj. R-squared:                  0.446\n",
       "Method:                 Least Squares   F-statistic:                     15.47\n",
       "Date:                Tue, 25 Sep 2018   Prob (F-statistic):           5.03e-16\n",
       "Time:                        15:55:20   Log-Likelihood:                 50.995\n",
       "No. Observations:                 145   AIC:                            -83.99\n",
       "Df Residuals:                     136   BIC:                            -57.20\n",
       "Df Model:                           8                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -0.2426      0.027     -8.980      0.000      -0.296      -0.189\n",
       "wusa           0.0774      0.043      1.803      0.074      -0.008       0.162\n",
       "jasper        -0.2288      0.078     -2.929      0.004      -0.383      -0.074\n",
       "westgreen      0.0096      0.042      0.229      0.819      -0.073       0.092\n",
       "chesapeake    -0.0321      0.034     -0.943      0.347      -0.099       0.035\n",
       "tornetrask     0.0927      0.045      2.057      0.042       0.004       0.182\n",
       "urals          0.1854      0.091      2.027      0.045       0.005       0.366\n",
       "mongolia       0.0420      0.046      0.917      0.361      -0.049       0.133\n",
       "tasman         0.1155      0.030      3.834      0.000       0.056       0.175\n",
       "==============================================================================\n",
       "Omnibus:                       12.501   Durbin-Watson:                   0.817\n",
       "Prob(Omnibus):                  0.002   Jarque-Bera (JB):               17.388\n",
       "Skew:                           0.490   Prob(JB):                     0.000168\n",
       "Kurtosis:                       4.384   Cond. No.                         12.9\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lmod = smf.ols(formula='nhtemp ~ wusa + jasper + westgreen + chesapeake + tornetrask + urals + mongolia + tasman',\n",
    "               data=globwarm).fit()\n",
    "lmod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot successive residuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X2MXNWZ5/Hv0005KeeFNtgJuO3GHq1xlImzGFqQEZokEBKDsguWA4FsoiVSMtYsQqshM95tlCiTSXZkT6wRWSlod6yd3WHyMhgI0/EMIM8QyGaE4og2DWFN4uCQxLhNAgE3UqY7uG0/+0dXNber71tV3ap7q+7vIyG6q27VPbe7/dxT5zznOebuiIhIuQzk3QAREek+BX8RkRJS8BcRKSEFfxGRElLwFxEpIQV/EZESUvAXESkhBX8RkRJS8BcRKaGz8m5AlJUrV/q6devyboYU2ImZkwCsWL4s55aIFMfBgwd/7e6rko4rbPBft24dExMTeTdDCuzeiecBuGF0bc4tESkOM/tFmuM07CMiUkIK/iIiJaTgLyJSQgr+IiIlpOAvIlJCCv4iIiWk4C8iUkIK/iIiJVTYRV4iIkHjk1Ps3n+Y49OzrB6qsmPLRrZuHs67WT1LwV9ECm98corb73+a2bnTAExNz3L7/U8D6AbQIg37iEjh7d5/eCHw183OnWb3/sM5taj3KfiLSOEdn55t6nFJpuAvIoW3eqja1OOSTMFfRApvx5aNVCuDix6rVgbZsWVjTi3qfZrwFZHCq0/qKtsnOwr+ItITtm4eVrDPUCbDPmZ2tZkdNrMjZjYWccxHzewZMztkZt/M4rwiItKatnv+ZjYI3Al8EDgGPG5m+9z9mcAxG4Dbgcvd/YSZva3d84qISOuy6PlfChxx9+fc/SRwN3BdwzF/ANzp7icA3P3FDM4rIiItyiL4DwPPB74/Vnss6ELgQjN7zMwOmNnVGZxXRERalMWEr4U85iHn2QC8H1gD/IuZvcvdpxe9kdl2YDvAyMhIBk0TEZEwWfT8jwFrA9+vAY6HHPNtd59z958Bh5m/GSzi7nvcfdTdR1etWpVB00REJEwWwf9xYIOZrTezZcBNwL6GY8aBKwDMbCXzw0DPZXBuERFpQdvB391PAbcC+4EfAfe4+yEz+6KZXVs7bD/wspk9AzwK7HD3l9s9t4iItCaTRV7u/iDwYMNjnw987cBnav+JiEjOVNtHRKSEFPxFREpIwV9EpIQU/EVESkjBX0SkhBT8RURKSMFfRKSEFPxFREpIwV9EpIQU/EVESkjBX0SkhBT8RURKSMFfRKSEFPxFREpIwV9EpIQU/EVESiiTzVxERKR145NT7N5/mOPTs6weqrJjy0a2bh7u6DkV/EVEcjQ+OcXt9z/N7NxpAKamZ7n9/qcBOnoDUPAXyUkevT0pnt37Dy8E/rrZudPs3n9YwV+k3+TV25PiOT4929TjWdGEr0gO4np7Ui6rh6pNPZ4VBX+RHOTV25Pi2bFlI9XK4KLHqpVBdmzZ2NHzKviL5CCv3p4Uz9bNw+zctonhoSoGDA9V2bltk7J9RPrRji0bF435Q7reniaJ+9PWzcNd/z0q+IvkoP4PvZlArkliyZKCv0hOmu3t5ZUSWET6BNQ+BX+RHqFJ4nn6BJSNTCZ8zexqMztsZkfMbCzmuOvNzM1sNIvzipSJJonnKU02G20HfzMbBO4ErgHeCXzMzN4ZctxbgP8M/KDdc4qUUV4pgUWjT0DZyKLnfylwxN2fc/eTwN3AdSHHfQn4MvDbDM4pUjp5pQQWjT4BZSOLMf9h4PnA98eAy4IHmNlmYK27/6OZ/UkG5xTpqqJMMOaRElg0rabJymJZBH8LecwXnjQbAO4APpn4Rmbbge0AIyMjGTRNpH2aYCyWVtJkZaksgv8xYG3g+zXA8cD3bwHeBXzXzADOA/aZ2bXuPhF8I3ffA+wBGB0ddUQKQCmWxaNPQO3LYsz/cWCDma03s2XATcC++pPu/qq7r3T3de6+DjgALAn8IkWlCUbpR20Hf3c/BdwK7Ad+BNzj7ofM7Itmdm277y+SN00wSj/KZJGXuz8IPNjw2Ocjjn1/FucU6RZNMEo/0gpfkQSaYJR+pOAvkoImGKXfqJ6/iEgJKfiLiJSQgr+ISAkp+IuIlJCCv4hICSn4i4iUkIK/iEgJKfiLiJSQgr+ISAkp+IuIlJDKO4h0WFF2ARMJUvAX6SDtAiZFpWEfkQ6K2wVMJE8K/iIdpF3ApKg07CMSo93x+tVDVaZCAr12AZO8qecvEqE+Xj81PYvz+nj9+ORU6vfYsWUj1crgose0C5gUgXr+0tOeOHqCrzz8bEcyaeLG69OeQ7uASVEp+EvPeuLoCe47eIy50w5kn0mT1Xi9dgGTItKwj/Ssh57+5ULgr8sykyZqXF7j9dIPFPylZ03PzoU+nlUmjcbrpZ9p2Ed61lC1EnoDyKpnrvF66Wfq+UvPumbTeVQGbdFj7fbMxyenuHzXI6wfe4DLdz0CwGNjV3LHjRcBcNveJ7l81yNNZfyIFJF6/tKzLh5ZAcD3fvLrTHrmUaUYJn7xCt86OKUSDdJXFPylp108soKd296dyXtFpXb+3Q+e57SHTywr+EuvymTYx8yuNrPDZnbEzMZCnv+MmT1jZj80s++Y2QVZnFckS1ETxY2BP+l4kV7QdvA3s0HgTuAa4J3Ax8zsnQ2HTQKj7v5u4D7gy+2eVyRrURPFg2ahjyvlU3pZFj3/S4Ej7v6cu58E7gauCx7g7o+6+0zt2wPAmgzOK5KpqNTOj122tq2Uz8ZJ5LjJ4maOFWlHFmP+w8Dzge+PAZfFHP8p4KEMziuSaXmHuNTO0QvOaSnls5l6/qr9L92URfAP+0wcOkhqZp8ARoH3RTy/HdgOMDIykkHTpJ91orxDVCmGVks0NFMfKItaQiJpZTHscwxYG/h+DXC88SAzuwr4LHCtu78W9kbuvsfdR919dNWqVRk0TfpZp8s7ZKGZ+kBlr/2vIa/uyiL4Pw5sMLP1ZrYMuAnYFzzAzDYDf8V84H8xg3OKdLy8QxaaqQ9U5lpCWZTPlua0Hfzd/RRwK7Af+BFwj7sfMrMvmtm1tcN2A28G7jWzJ81sX8TbiaQ2VK2EPt4YLPPsUTZTH6jMtYS03WX3ZbLIy90fBB5seOzzga+vyuI8IkHXbDpv0Zg/LA2WeU2iBncAG1pe4Q1nDfDq7FzsZHG7tYTa3XUsT2Uf8sqDVvhKV2UZoNKUd/izfzjU9UnUxhvOiZk5qpVB7rjxosRztjqx3OuZQtrusvsU/KVrOhGg4so7jE9OcWImfF4gLNA0vrbVm1QeWTu9nim0Y8vGRX8bUJ4hr7wo+EvXpA1Q7X46qL8+LsBb7biw9233JpXHEEavD5uofHb3KfhL16QJUO0G3sbXR3GI7BW3e5PKYwijH4ZNtN1ld6mev3RNmlTGdrM+wl4fpdnecthNKiw1MSxrx4Ar3tG5tStlzhSS1ij4S9ekCVDNBuQnjp5YlMaZNJYf1Gxefdqb1NbNw3zkkuFFS98d+NbBqY6lmW7dPMzObZsYHqpiwPBQlZ3bNqknLZE07FMyeaYDphnXjRq+GDBj/dgDi14TVt7BiKgt0iCuV5xm8jHpJvXoj19a0o5OT8Bq2ESaoeBfIkVIB0wKUGGBF16vqR9sc1h5B4clN4BqZZCPXDLMoz9+KfSmE3ZD3LltU0s3qfqng2Y/wfRyjr70JgX/Bv38j7AX0gEbPx0MmEXuohVV3sGZH/ZI8zuMuiHu3LaJx8aujGxn2E0qOK7fzARsJ27K/fx3LNlQ8A8oQs+4k7JIB+xGUAl+Olg/9kDoMcenZ6lWBkMnd4eqldjAHdTsDTF4/dXK4imz+rj+6AXnNJW3nvVNud//jiUbmvAN6Pf6Iu0WDsuj+FZcmyM22OLkqXTZPtDcDbHx+mfmziw5Jhi0007AZp2j3+9/x5INBf+AXl8ok6TddMA8gkpcm2dOhgf5mbkzqW9IzdwQ06aR1v9etm4e5rGxK/nZrg/z2NiVkb3urKt59vvfsWRDwT+g30vqBnujML83bT14pwmWeQSVuB50VFVPIPUNKSr3ft251SWVQNOmkTb795J1jn6//x1LNsw9TWJc942OjvrExERXzxm2OrRaGey7fOlWrzMqj354qJp6jL0djfMNI+dU+f5zr0Qeb5CY2TNz8lRo/Z/GjKHKoC3JLArT6t9LlnMpZfk7lnBmdtDdRxOPU/BfrAxZEq0G8TyDSti5K4OGAScTgnK9jUCq0g/NGqpWEss1d1sZ/o4lnIK/LGgMBFHDFwb8bNeHU7/X2dUKZjA9M8fQ8grudCwIRt2wli8bxJ3EgF4f6mpmBXBaP0/4mYl0U9rgrzH/PheWoRORJJNqTLg+iXnHjRfx2qkznJiZw5mvWT89O9exLKCoeYWZk6cXzQnEvT7t3ETc+zQa1ji69CgF/x7T7JaEYRkq9VWwQc1OMCZlvmSdBRR1YxqqVhZl1UQF49VD1dj3CE4of/w9I0smYCuDRmVg8U/NmL/RabNx6UVa5NVDWlm8E9XbbWYVbDPv2+wxaYUtmqoMGtdsOm/h+/HJKf71tVNLXhu8sYXNWXzh2t9dcu2jF5yzZMwcWNgnIDghXP89TPzilcgSEiJFo+DfQ1pZCRo1xt9uhk7c3EHwmKyEFYV774UrgdfnA8KKuq1YXuFP//3i4J5mIjSqBtHWzcOh8w+zc6f5xoGjS24IwbaLFImCfw9pJc9+x5aN7Lj3KebOvB4WKwPWdp33qAJsdfXedpZZJ40B+fb7f7ioqmdY6sKJmbmF4af669sNxnGfpoKKVjdJJEhj/j2k5cU7jQP8zcxoRmhcfLVieYWhamXRQiygo+Ugwqp6hmn1vFHzK818otGcgBSVUj17SCt59nkuzGr23M1+SlgXUfQtSjPXHPezhqVzB0n7CBRlkVW38/+13qD7lOrZh1rZramTJRmSMo+izjE1PbvkNa0UjYsr7xCmmWtOml9p/D2EZQiFvTZP3S7Ml0chQElPY/49ptkx605t7J0m8yhuUtgbXtPKZPY1m85bNOafpJlrTrpp1n8P9Z7tNw4c5exqhTdWBkLLRcS9Z7d0ez+HXtg/oszU8+9zndrYO02Fz7BzN6q/Ju5TQtSY+cUjK7j+kjULuf1JUxnrzk0f/NPMrzT2bKdn5/jt3BlWLA//RJJ082l2DUezul2YT9VFi03Bv891amPvuH/Y9SB2294neWNlYGEiOMpUbceuuOdv2/sknxt/eslzF4+s4LGxKxkeqibu3XvguRMJR7wegMNWQgdvmuOTU/zxPU+F3gDdafqG240hkm5X+1R10WLLJPib2dVmdtjMjpjZWMjzbzCzvbXnf2Bm67I4r6STtq58MyJXyy6vLApiJ2bmeO3UGe648aLI1bcGS7ZqbOTANw4cjQyGaXqTSecIBuD6Oes3gOBNs35c1Pu9OjvX9A23G3sldOpTYFHOJ81pe8zfzAaBO4EPAseAx81sn7s/EzjsU8AJd/83ZnYT8BfAje2eW/ITteJ2ulbrJ6gexKL2vU2bb+YQOV6cZtFZ0rBQVCmMxiyhpNIWq4eqTc/NdGOIJGyhXCezb7p9PmlOFhO+lwJH3P05ADO7G7gOCAb/64Av1L6+D/iqmZnH5JmemDnJvRPPZ9A86ZStm1fz0NO/ZHp2juXLBnlt7nRkIJ+anmXu9JlFrxmqViI3YY8yNT278Hfx+M9fr+X/3gtXJk/+2vzCsItHVkS+d9I5nzh6IvYmUxk03nvhykV/u08cPbHomq/ZdN6SNpwd8bM4u1rJ/N/BH121YeHrudNnOv7vrNvnk3SyCP7DQPC3eQy4LOoYdz9lZq8C5wK/Dh5kZtuB7QCr1qzPoGmSJE1ginLxyIqFY//8gR9FbqsIr6dlBl9Tf10zN4Co9M76e9avJYz7/PNR1xd1M6qf84mjJ7jv4LHIthlw/SVrFr1//TX1m9L07NzCewSPC8tcaqxdJJKlLIJ/2Kfpxu5XmmNw9z3AHphf5HXD6Nr2WyeRxienGJ88vjCEMT07x/jkcS5bf27TH83/y30/jHwuqngaQGVwIPUGK1HvU/87uWF0LTu3vRuA9WMPhH4KeXV2jqi/q6eOTfP1A0eXPP7v/u353DC6lq88/GzkJ4uoRVxhr5k77XzvJ79eaGu97ZetP1dDJNI1WQT/Y0DwX9Ma4HjEMcfM7CzgbCB6/z3pirR52GlWaUaNuQ+axU52bt08zMQvXllUFK3uTcsGqQwOtLRBTCvrGx798Uuxj8eNv0ddYzNj+UnzBFotK1nKIvg/Dmwws/XAFHAT8B8ajtkH3Ax8H7geeCRuvF+6Iy63vu5z40+nqlYZNpmbtqTBoz9+KbSXPrR8WapyDGFBMao9cZkmSYE6rkIqzJezaAzMWS2ya6Wct0ictoN/bQz/VmA/MAj8b3c/ZGZfBCbcfR/w18DXzOwI8z3+m9o9r7QvKjAZLKRUhvXIwz4dtJPZ0WymSz3YT03PsnzZIHOnzixULa0HxZ3bNrFz26bI9oTdMJICddQN5Yp3rIoMzK3chMJotaxkTYXdSmx8corb9j4Z2utO2vM2zX6/aTVTAC6s4FqY4GsbA/0V71jFtw5OLQnIH7lkOPTx4KeXsJtG/UYU1YYshmui5jCy/D1If0hb2E21fUps6+Zh/mjvk6HPJeWXhw1btBrk0vSOg739NOrtDxsuifo08+iPX4r9tADh4/K3JfwMs9hDoFM1mqS8FPxLbjghqEQNCzUOW7QzJp00ZJS2tx/W/qiFW2GOT8+2FKi7EZizGj4SqVNtn5KLW4If9pwBH3/PyJIA2Wp5gmAdIIA7brxoSQmKpBW1jYJBsZkVsq0G626UMehUjSYpL/X8+0A7Y8ppJmqT3nt8cipyOKYx+I5PTvFn/3AotOxx1EboSUM9lUHjTcvOCk0JjZvUDn4CaCdYd6uMQRbDRyJ1mvDtcWl39+pUjnjSkEzjxOuO+55KrL/fGJjj6v8MVSuRC8ii2lef3A3eYJQzL/1CE74lkSYFsJM54nFDMo296d37D6faeKXxiHp1zcae+tbNq7l4ZEXsNUT1yiF6UZdIGSj497g0OfKdzBFvZtVrOxUq69U1gwF87vSZVK9tHC7RgikRBf+elybTJM1K3qzPP1wra5zm2KCoIZ6wnP9Wq0NqwZSIsn16XppMk6gsluBK3mYEtxv819dOURlcXLcvavJ0x5aNDMQU1Y/aCD3rzBltLyiinn9PiZu0jZvM3bFlY+hK3rjNUeLaEBwymZ6dozJgrFheYXomuQDb4IBxpmHc/xPvGeG/bd208P3oBed0NHNGC6ZEFPwLK6kkQeM4ddKkZ6sreRuFDZnMnXGWLzuLyc9/aKHtYUXOoiZ8GydeO53SqAVTIgr+hRRWSTNtgbUoSSt500oaMombTO3EcMsTR0/wlYefbepTQl7bC6oksxSJgn/BjE9OhQb6uJIEce9VDzZDyytUBmyh+iW01tsdWl4JXaA1tHx+t6u4ydS0wy1pg2TjLlnNlpToZuBVhpEUjSZ8C2b3/sOpNzSH6J57PdhMTc/iMB+wbX5RVDvlAaLWBE7PzDE+ORXbu08zOd3Y7nqQDJuYfujpXy4ZRkpTUiIPrZa/EOkU9fwLJi4VspmSBKFj86edN73hLJ780w+13L5Xo/bHBW6//+nIjchXB1I/k0pJpE3DjNqrt4hZO8owkqJR8C+YQTNOR3SvP/6ekdQlCToVbOJy9WfnTvPGygDVymDkZGrScEsz7Y7acL1oWTvjk1MMRPxei9ZWKQ8F/4KJCvzAonTIJJ1KZwzLlAmanpnjjhsvanliM6ndwfmA6rJBBg2CIz9Fy9qpD2OF/V6L1lYpFwX/gonKyhluMmh3Kp2xHsT/+J6nInuy7UymxrW7cdJ05uRpBmvzGK1s8t4NUbWPkja2F+k0Bf+CySpodzKdsf4enby5hLX78l2PLAmkp5225zE6KWoY64y7Ar/kSsG/YLIM2p1MZ+z0zSXsfXpx0lSriaWoFPwLKMug3cmFRd3Olc8qkHZzsZVWE0tRKfj3sbQLi3pl5WlYIK0MWlOBtNuLrfJaTSySRMG/jyXlzIdtqdiNlaet3mwaA+nZ1QrXbDqvqXbmUc5Z2y9KESn497G4MfK47Rc7GQzb7XkHA2ljPf80N5VenDcQ6QQF/z4WN0Yet/0iRAfDxnpB7jSVZtmpnnfam4omYEXmqbZPH4urpZPU0w0LhmH1gqZn5xJr8AR1quedtnZO2vpC9c1qLt/1yKJrintOpJe01fM3s3OAvcA64OfAR939RMMxFwH/A3grcBr4c3ff2855JZ24ycbd+w9HlmmIykZJ+rSQ1IPvZJmDtDeVsJ/JFe9Yxe79h7lt75MMLa/wm9+eWqh+GvwEAfTVBLqUW7vDPmPAd9x9l5mN1b7/rw3HzAD/0d2fNbPVwEEz2+/u022eW1KImmyMKtMwVK3whWt/t6k8+zTHdLrMQVRBubOrlSWPBX8mjcNFYeWqg58gkoasVLpZekW7wf864P21r+8CvktD8Hf3nwS+Pm5mLwKrAAX/HLWSgphmA/aoHnynyxxYxN7AUY8ntatR3I0v+Jw2h5de0W7wf7u7vwDg7i+Y2dviDjazS4FlwE/bPK9koNkUxKSibnE9+E6XOZgO6bHHPZ7Urkb1m1rSZLGyiaRXJAZ/M3sYOC/kqc82cyIzOx/4GnCzu5+JOGY7sB1gZGSkmbfvOb04Ltz4aaGZbJ9OZ9m0+v5pPs0Eb2pJq3WVTSS9IjH4u/tVUc+Z2a/M7Pxar/984MWI494KPAB8zt0PxJxrD7AHYHR0tJkNrXpKL48Lt7pgqdNlDlp9/6hVw29adlbkTS3upq1yDtIr2h322QfcDOyq/f/bjQeY2TLg74G/dfd72zxfXyjjuHCnyxy0+v7Nvi7p5qdyDtIrzGM2D0l8sdm5wD3ACHAUuMHdXzGzUeAP3f3TZvYJ4P8AhwIv/aS7Pxn33qOjoz4xMdFy24ps/dgDofv0GvCzXR/udnN6Vn2F7w2ja3NuiUhxmNlBdx9NOq6tnr+7vwx8IOTxCeDTta+/Dny9nfP0m6hx4QEzxien+rKX2ItzHCL9TCt8cxC2yhTmt3BMs0q21zSuDE67GlhEOkfBPwdbNw+zc9smBkOS0MNKEvS6tKUXRKR7FPxzsnXzMGci5lv6LSdcue8ixaOqnl3UOO4dVZKg33LC281913yBSPYU/LskLLe/MmhUBmyhiBhknxNehMDZTu573JoIEWmdhn26JGzce+608+Y3nsXwUBUDhoeqmdS5qSvKRGt9jqOV69R8gUhnqOffJVHj29Mzc0x+/kMdOWeRFpO1ujI46uc2NT3Ljvt+yFC1QmVwQMNAIk1Sz79Losa3Ozm+Hxc4e2UjkqSfz/TsnNJGRVqg4N+iZnd0SrODVNbiAmev5NpHrYkI0jCQSPMU/FvQylh6O+PerUoKnL0QNBt/blGUNirSHI35t6DVsfRWx71bFSwyFlW2uBeCZvDndvmuR1pKGy1C1pNIkajn34JeWrS0dfMwj41dyXAOcw6d0MrwWVGynkSKRMG/BXlM3rYrjzmHTggOA8H8nsNJw2dKFxVZSsM+LejFDTv6qc58fRioXtI56Rp66ZOaSLco+LegVwNpt+ccikJbK4ospeDforIG0l7Ui5/URDpNwV/6Xq9+UhPpJAV/KQV9UhNZrG+Dv/K6RUSi9WXwjysDrBuAiEif5vkrr1tEJF5fBn/ldYuIxOvL4N+LK3BFRLqpL4N/v5QyEBHplL6c8FVet4hIvL4M/qC87l6nVF2Rzmpr2MfMzjGzfzazZ2v/XxFz7FvNbMrMvtrOOaX/qQSzSOe1O+Y/BnzH3TcA36l9H+VLwP9t83xSAkrVFem8doP/dcBdta/vAraGHWRmlwBvB/6pzfNJCShVV6Tz2g3+b3f3FwBq/39b4wFmNgD8JbCjzXNJSShVV6TzEoO/mT1sZv8v5L/rUp7jFuBBd38+xbm2m9mEmU289NJLKd9e+o1SdUU6LzHbx92vinrOzH5lZue7+wtmdj7wYshhvwf8vpndArwZWGZmv3H3JfMD7r4H2AMwOjrqaS9C+otSdUU6r91Uz33AzcCu2v+/3XiAu3+8/rWZfRIYDQv8IkFK1RXprHaD/y7gHjP7FHAUuAHAzEaBP3T3T7f5/n1P+ewikoe2gr+7vwx8IOTxCWBJ4Hf3vwH+pp1z9hOVnhaRvPRlbZ9eoXx2EcmLgn+OlM8uInlR8M+R8tlFJC8K/jlSPruI5KVvq3r2AuWzi0heFPxzpnx2EcmDhn1EREpIwV9EpIQU/EVESkjBX0SkhBT8RURKSMFfRKSEFPxFREpIwV9EpITMvZgbZpnZS8AvunjKlcCvu3i+btK19SZdW2/K+9oucPdVSQcVNvh3m5lNuPto3u3oBF1bb9K19aZeuTYN+4iIlJCCv4hICSn4v25P3g3oIF1bb9K19aaeuDaN+YuIlJB6/iIiJVTa4G9m55jZP5vZs7X/r4g59q1mNmVmX+1mG1uV5trM7CIz+76ZHTKzH5rZjXm0NS0zu9rMDpvZETMbC3n+DWa2t/b8D8xsXfdb2ZoU1/YZM3um9nv6jpldkEc7W5F0bYHjrjczN7PCZ8nUpbk2M/to7Xd3yMy+2e02xnL3Uv4HfBkYq309BvxFzLH/Hfgm8NW8253VtQEXAhtqX68GXgCG8m57xPUMAj8FfgdYBjwFvLPhmFuA/1n7+iZgb97tzvDargCW177+T/10bbXj3gJ8DzgAjObd7gx/bxuASWBF7fu35d3u4H+l7fkD1wF31b6+C9gadpCZXQK8HfinLrUrC4nX5u4/cfdna18fB14EEheG5ORS4Ii7P+fuJ4G7mb/GoOA13wd8wMysi21sVeK1ufuj7j5T+/YAsKbLbWxVmt8bwJeY77D8tpuNa1Oaa/sD4E53PwHg7i92uY2xyhz83+7uLwDU/v+2xgPMbAD4S2BHl9vWrsRrCzKzS5nvvfy0C21rxTDwfOD7Y7XHQo9x91PAq8DS6eOPAAACR0lEQVS5XWlde9JcW9CngIc62qLsJF6bmW0G1rr7P3azYRlI83u7ELjQzB4zswNmdnXXWpdCX+/ha2YPA+eFPPXZlG9xC/Cguz9ftE5kBtdWf5/zga8BN7v7mSza1gFhP/zGNLU0xxRR6nab2SeAUeB9HW1RdmKvrda5ugP4ZLcalKE0v7ezmB/6eT/zn9b+xcze5e7THW5bKn0d/N39qqjnzOxXZna+u79QC4BhH8l+D/h9M7sFeDOwzMx+4+6RE1fdksG1YWZvBR4APufuBzrU1CwcA9YGvl8DHI845piZnQWcDbzSnea1Jc21YWZXMX9jf5+7v9altrUr6dreArwL+G6tc3UesM/MrnX3ia61sjVp/yYPuPsc8DMzO8z8zeDx7jQxXpmHffYBN9e+vhn4duMB7v5xdx9x93XAnwB/W4TAn0LitZnZMuDvmb+me7vYtlY8Dmwws/W1dt/E/DUGBa/5euARr82yFVzitdWGRv4KuLZo48YJYq/N3V9195Xuvq72b+wA89dY9MAP6f4mx5mfrMfMVjI/DPRcV1sZo8zBfxfwQTN7Fvhg7XvMbNTM/leuLWtfmmv7KPBe4JNm9mTtv4vyaW682hj+rcB+4EfAPe5+yMy+aGbX1g77a+BcMzsCfIb5LKfCS3ltu5n/5Hlv7ffUGGQKKeW19aSU17YfeNnMngEeBXa4+8v5tHgprfAVESmhMvf8RURKS8FfRKSEFPxFREpIwV9EpIQU/EVESkjBX0SkhBT8RURKSMFfRKSE/j+Z+ZlmZtHF9QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(lmod.resid.iloc[:-1],lmod.resid.iloc[1:])\n",
    "plt.axhline(0,alpha=0.5)\n",
    "plt.axvline(0,alpha=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the correlation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.58333899],\n",
       "       [0.58333899, 1.        ]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.corrcoef(lmod.resid.iloc[:-1],lmod.resid.iloc[1:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the uncorrelated model and then iterate the fit. Produces result very similar to the correlation of the residuals. Different from the R correlation of 0.71"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.58221337])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "globwarm = globwarm.dropna()\n",
    "X = sm.add_constant(globwarm.iloc[:,1:9])\n",
    "gmod = sm.GLSAR(globwarm.nhtemp, X, rho=1)\n",
    "res=gmod.iterative_fit(maxiter=6)\n",
    "gmod.rho"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternative method that shows the iterative nature of the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AR coefficients: [0.]\n",
      "AR coefficients: [0.57494252]\n",
      "AR coefficients: [0.58194786]\n",
      "AR coefficients: [0.58220381]\n",
      "AR coefficients: [0.58221337]\n",
      "AR coefficients: [0.58221373]\n"
     ]
    }
   ],
   "source": [
    "gmod = sm.GLSAR(globwarm.nhtemp, X, rho=1)\n",
    "for i in range(6):\n",
    "     results = gmod.fit()\n",
    "     print(\"AR coefficients: {0}\".format(gmod.rho))\n",
    "     rho, sigma = sm.regression.yule_walker(results.resid,order=gmod.order)\n",
    "     gmod = sm.GLSAR(globwarm.nhtemp, X, rho)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in the data for the blocked example. Need to make variety categorical. For reasons unclear to me, the mixed effect modeling function throws an error if the response is called *yield*. We create the same variable but with a different name: *grams*."
   ]
  },
  {
   "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>yield</th>\n",
       "      <th>block</th>\n",
       "      <th>variety</th>\n",
       "      <th>grams</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>296</td>\n",
       "      <td>I</td>\n",
       "      <td>1</td>\n",
       "      <td>296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>357</td>\n",
       "      <td>II</td>\n",
       "      <td>1</td>\n",
       "      <td>357</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>340</td>\n",
       "      <td>III</td>\n",
       "      <td>1</td>\n",
       "      <td>340</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>331</td>\n",
       "      <td>IV</td>\n",
       "      <td>1</td>\n",
       "      <td>331</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>348</td>\n",
       "      <td>V</td>\n",
       "      <td>1</td>\n",
       "      <td>348</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   yield block variety  grams\n",
       "1    296     I       1    296\n",
       "2    357    II       1    357\n",
       "3    340   III       1    340\n",
       "4    331    IV       1    331\n",
       "5    348     V       1    348"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "oatvar = pd.read_csv(\"data/oatvar.csv\", index_col=0)\n",
    "oatvar['variety'] = oatvar['variety'].astype('category')\n",
    "oatvar['grams'] = oatvar['yield']\n",
    "oatvar.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can do mixed effect model but would need to calculate correlation from this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "       <td>Model:</td>       <td>MixedLM</td> <td>Dependent Variable:</td>   <td>grams</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <td>No. Observations:</td>   <td>40</td>          <td>Method:</td>         <td>REML</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "     <td>No. Groups:</td>       <td>5</td>          <td>Scale:</td>        <td>1336.9044</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <td>Min. group size:</td>     <td>8</td>        <td>Likelihood:</td>     <td>-170.6771</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <td>Max. group size:</td>     <td>8</td>        <td>Converged:</td>         <td>Yes</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <td>Mean group size:</td>    <td>8.0</td>            <td></td>               <td></td>     \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "        <td></td>        <th>Coef.</th>  <th>Std.Err.</th>    <th>z</th>   <th>P>|z|</th> <th>[0.025</th>  <th>0.975]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>    <td>334.400</td>  <td>21.040</td>  <td>15.894</td> <td>0.000</td> <td>293.162</td> <td>375.638</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.2]</th> <td>42.200</td>   <td>23.125</td>   <td>1.825</td> <td>0.068</td> <td>-3.124</td>  <td>87.524</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.3]</th> <td>28.200</td>   <td>23.125</td>   <td>1.219</td> <td>0.223</td> <td>-17.124</td> <td>73.524</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.4]</th> <td>-47.600</td>  <td>23.125</td>  <td>-2.058</td> <td>0.040</td> <td>-92.924</td> <td>-2.276</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.5]</th> <td>105.000</td>  <td>23.125</td>   <td>4.541</td> <td>0.000</td> <td>59.676</td>  <td>150.324</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.6]</th> <td>-3.800</td>   <td>23.125</td>  <td>-0.164</td> <td>0.869</td> <td>-49.124</td> <td>41.524</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.7]</th> <td>-16.000</td>  <td>23.125</td>  <td>-0.692</td> <td>0.489</td> <td>-61.324</td> <td>29.324</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.8]</th> <td>49.800</td>   <td>23.125</td>   <td>2.154</td> <td>0.031</td>  <td>4.476</td>  <td>95.124</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Group Var</th>    <td>876.492</td>  <td>21.576</td>     <td></td>      <td></td>       <td></td>        <td></td>    \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary2.Summary'>\n",
       "\"\"\"\n",
       "          Mixed Linear Model Regression Results\n",
       "==========================================================\n",
       "Model:             MixedLM  Dependent Variable:  grams    \n",
       "No. Observations:  40       Method:              REML     \n",
       "No. Groups:        5        Scale:               1336.9044\n",
       "Min. group size:   8        Likelihood:          -170.6771\n",
       "Max. group size:   8        Converged:           Yes      \n",
       "Mean group size:   8.0                                    \n",
       "----------------------------------------------------------\n",
       "              Coef.  Std.Err.   z    P>|z|  [0.025  0.975]\n",
       "----------------------------------------------------------\n",
       "Intercept    334.400   21.040 15.894 0.000 293.162 375.638\n",
       "variety[T.2]  42.200   23.125  1.825 0.068  -3.124  87.524\n",
       "variety[T.3]  28.200   23.125  1.219 0.223 -17.124  73.524\n",
       "variety[T.4] -47.600   23.125 -2.058 0.040 -92.924  -2.276\n",
       "variety[T.5] 105.000   23.125  4.541 0.000  59.676 150.324\n",
       "variety[T.6]  -3.800   23.125 -0.164 0.869 -49.124  41.524\n",
       "variety[T.7] -16.000   23.125 -0.692 0.489 -61.324  29.324\n",
       "variety[T.8]  49.800   23.125  2.154 0.031   4.476  95.124\n",
       "Group Var    876.492   21.576                             \n",
       "==========================================================\n",
       "\n",
       "\"\"\""
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mmod = smf.mixedlm(\"grams ~ variety\",oatvar,groups=oatvar['block']).fit()\n",
    "mmod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GEE model is also possible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>GEE Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>                <td>grams</td>        <th>  No. Observations:  </th>    <td>40</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                         <td>GEE</td>         <th>  No. clusters:      </th>     <td>5</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                    <td>Generalized</td>     <th>  Min. cluster size: </th>     <td>8</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th></th>                      <td>Estimating Equations</td> <th>  Max. cluster size: </th>     <td>8</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Family:</th>                     <td>Gaussian</td>       <th>  Mean cluster size: </th>    <td>8.0</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Dependence structure:</th>     <td>Exchangeable</td>     <th>  Num. iterations:   </th>     <td>2</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>                   <td>Tue, 25 Sep 2018</td>   <th>  Scale:             </th> <td>2213.400</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance type: </th>            <td>robust</td>        <th>  Time:              </th> <td>15:55:20</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "        <td></td>          <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>    <td>  334.4000</td> <td>    9.409</td> <td>   35.541</td> <td> 0.000</td> <td>  315.959</td> <td>  352.841</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.2]</th> <td>   42.2000</td> <td>   22.420</td> <td>    1.882</td> <td> 0.060</td> <td>   -1.743</td> <td>   86.143</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.3]</th> <td>   28.2000</td> <td>   32.653</td> <td>    0.864</td> <td> 0.388</td> <td>  -35.798</td> <td>   92.198</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.4]</th> <td>  -47.6000</td> <td>   17.136</td> <td>   -2.778</td> <td> 0.005</td> <td>  -81.186</td> <td>  -14.014</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.5]</th> <td>  105.0000</td> <td>   23.634</td> <td>    4.443</td> <td> 0.000</td> <td>   58.678</td> <td>  151.322</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.6]</th> <td>   -3.8000</td> <td>   14.750</td> <td>   -0.258</td> <td> 0.797</td> <td>  -32.709</td> <td>   25.109</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.7]</th> <td>  -16.0000</td> <td>   26.032</td> <td>   -0.615</td> <td> 0.539</td> <td>  -67.022</td> <td>   35.022</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>variety[T.8]</th> <td>   49.8000</td> <td>   34.238</td> <td>    1.455</td> <td> 0.146</td> <td>  -17.306</td> <td>  116.906</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>      0.0045</td> <th>  Kurtosis:          </th> <td>     -0.2204</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Centered skew:</th> <td>      0.4576</td> <th>  Centered kurtosis: </th> <td>      0.0147</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                               GEE Regression Results                              \n",
       "===================================================================================\n",
       "Dep. Variable:                       grams   No. Observations:                   40\n",
       "Model:                                 GEE   No. clusters:                        5\n",
       "Method:                        Generalized   Min. cluster size:                   8\n",
       "                      Estimating Equations   Max. cluster size:                   8\n",
       "Family:                           Gaussian   Mean cluster size:                 8.0\n",
       "Dependence structure:         Exchangeable   Num. iterations:                     2\n",
       "Date:                     Tue, 25 Sep 2018   Scale:                        2213.400\n",
       "Covariance type:                    robust   Time:                         15:55:20\n",
       "================================================================================\n",
       "                   coef    std err          z      P>|z|      [0.025      0.975]\n",
       "--------------------------------------------------------------------------------\n",
       "Intercept      334.4000      9.409     35.541      0.000     315.959     352.841\n",
       "variety[T.2]    42.2000     22.420      1.882      0.060      -1.743      86.143\n",
       "variety[T.3]    28.2000     32.653      0.864      0.388     -35.798      92.198\n",
       "variety[T.4]   -47.6000     17.136     -2.778      0.005     -81.186     -14.014\n",
       "variety[T.5]   105.0000     23.634      4.443      0.000      58.678     151.322\n",
       "variety[T.6]    -3.8000     14.750     -0.258      0.797     -32.709      25.109\n",
       "variety[T.7]   -16.0000     26.032     -0.615      0.539     -67.022      35.022\n",
       "variety[T.8]    49.8000     34.238      1.455      0.146     -17.306     116.906\n",
       "==============================================================================\n",
       "Skew:                          0.0045   Kurtosis:                      -0.2204\n",
       "Centered skew:                 0.4576   Centered kurtosis:              0.0147\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ind = sm.cov_struct.Exchangeable()\n",
    "gmod = smf.gee(\"grams ~ variety\", \"block\", oatvar, cov_struct = ind ).fit()\n",
    "gmod.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'The correlation between two observations in the same cluster is 0.336'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ind.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Weighted Least Squares\n",
    "\n",
    "Read in the French Presidential Election data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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>EI</th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "      <th>C</th>\n",
       "      <th>D</th>\n",
       "      <th>E</th>\n",
       "      <th>F</th>\n",
       "      <th>G</th>\n",
       "      <th>H</th>\n",
       "      <th>J</th>\n",
       "      <th>K</th>\n",
       "      <th>A2</th>\n",
       "      <th>B2</th>\n",
       "      <th>N</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Ain</th>\n",
       "      <td>260</td>\n",
       "      <td>51</td>\n",
       "      <td>64</td>\n",
       "      <td>36</td>\n",
       "      <td>23</td>\n",
       "      <td>9</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>4</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>105</td>\n",
       "      <td>114</td>\n",
       "      <td>17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Alpes</th>\n",
       "      <td>75</td>\n",
       "      <td>14</td>\n",
       "      <td>17</td>\n",
       "      <td>9</td>\n",
       "      <td>9</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>32</td>\n",
       "      <td>31</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ariege</th>\n",
       "      <td>107</td>\n",
       "      <td>27</td>\n",
       "      <td>18</td>\n",
       "      <td>13</td>\n",
       "      <td>17</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>57</td>\n",
       "      <td>33</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bouches.du.Rhone</th>\n",
       "      <td>1036</td>\n",
       "      <td>191</td>\n",
       "      <td>204</td>\n",
       "      <td>119</td>\n",
       "      <td>205</td>\n",
       "      <td>29</td>\n",
       "      <td>13</td>\n",
       "      <td>13</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>6</td>\n",
       "      <td>466</td>\n",
       "      <td>364</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Charente.Maritime</th>\n",
       "      <td>367</td>\n",
       "      <td>71</td>\n",
       "      <td>76</td>\n",
       "      <td>47</td>\n",
       "      <td>37</td>\n",
       "      <td>8</td>\n",
       "      <td>34</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>163</td>\n",
       "      <td>142</td>\n",
       "      <td>17</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     EI    A    B    C    D   E   F   G   H   J  K   A2   B2  \\\n",
       "Ain                 260   51   64   36   23   9   5   4   4   3  3  105  114   \n",
       "Alpes                75   14   17    9    9   3   1   2   1   1  1   32   31   \n",
       "Ariege              107   27   18   13   17   2   2   2   1   1  1   57   33   \n",
       "Bouches.du.Rhone   1036  191  204  119  205  29  13  13  10  10  6  466  364   \n",
       "Charente.Maritime   367   71   76   47   37   8  34   5   4   4  2  163  142   \n",
       "\n",
       "                    N  \n",
       "Ain                17  \n",
       "Alpes               5  \n",
       "Ariege              6  \n",
       "Bouches.du.Rhone   30  \n",
       "Charente.Maritime  17  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fpe = pd.read_csv(\"data/fpe.csv\",index_col=0)\n",
    "fpe.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the model with weights:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A    1.067130\n",
       "B   -0.105051\n",
       "C    0.245958\n",
       "D    0.926188\n",
       "E    0.249397\n",
       "F    0.755110\n",
       "G    1.972212\n",
       "H   -0.566217\n",
       "J    0.611642\n",
       "K    1.210658\n",
       "N    0.529353\n",
       "dtype: float64"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wmod = smf.wls(\"A2 ~ A + B + C + D + E + F + G + H + J + K +N -1\", data=fpe, weights = 1/fpe.EI ).fit()\n",
    "wmod.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the model without weights:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A    1.075148\n",
       "B   -0.124559\n",
       "C    0.257447\n",
       "D    0.904537\n",
       "E    0.670677\n",
       "F    0.782533\n",
       "G    2.165655\n",
       "H   -0.854288\n",
       "J    0.144419\n",
       "K    0.518133\n",
       "N    0.558272\n",
       "dtype: float64"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lmod = smf.wls(\"A2 ~ A + B + C + D + E + F + G + H + J + K +N -1\", data=fpe).fit()\n",
    "lmod.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiplying the weights by a constant makes no difference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A    1.067130\n",
       "B   -0.105051\n",
       "C    0.245958\n",
       "D    0.926188\n",
       "E    0.249397\n",
       "F    0.755110\n",
       "G    1.972212\n",
       "H   -0.566217\n",
       "J    0.611642\n",
       "K    1.210658\n",
       "N    0.529353\n",
       "dtype: float64"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wmod = smf.wls(\"A2 ~ A + B + C + D + E + F + G + H + J + K +N -1\", data=fpe, weights = 53/fpe.EI ).fit()\n",
    "wmod.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit a reduced model. Doesn't seem to be an offset option so need to modify the response:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "C    0.225773\n",
       "D    0.969977\n",
       "E    0.390204\n",
       "F    0.744240\n",
       "N    0.608539\n",
       "dtype: float64"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = fpe.A2 - fpe.A - fpe.G - fpe.K\n",
    "X = fpe.loc[:,[\"C\",\"D\",\"E\",\"F\",\"N\"]]\n",
    "wmod = sm.WLS(y, X, weights = 53/fpe.EI ).fit()\n",
    "wmod.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Contrained least squares - need some trickery to get the weights right."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = fpe.A2\n",
    "X = fpe.loc[:,[\"A\",\"B\",\"C\",\"D\",\"E\",\"F\",\"G\",\"H\",\"J\",\"K\",\"N\"]]\n",
    "weights = 1/fpe.EI\n",
    "Xw = (X.T * np.sqrt(weights)).T\n",
    "yw = y * np.sqrt(weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A    1.000\n",
       "B    0.000\n",
       "C    0.208\n",
       "D    0.969\n",
       "E    0.359\n",
       "F    0.743\n",
       "G    1.000\n",
       "H    0.367\n",
       "J    0.000\n",
       "K    1.000\n",
       "N    0.575\n",
       "dtype: float64"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = sp.optimize.lsq_linear(Xw, yw, bounds=(0, 1))\n",
    "pd.Series(np.round(res.x,3),index=lmod.params.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IRWLS example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "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>speed</th>\n",
       "      <th>dist</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>7</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7</td>\n",
       "      <td>22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>8</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   speed  dist\n",
       "0      4     2\n",
       "1      4    10\n",
       "2      7     4\n",
       "3      7    22\n",
       "4      8    16"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cars = pd.read_csv(\"data/cars.csv\")\n",
    "cars.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show that heteroscedascity is present"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFsBJREFUeJzt3X2MXFd9xvHvk40Dy5s2IUuI13Yd2uAScInRKA11hVAS5ABR7LpQQqF121QuUlChLyY2Qa1QiezUFS9tEWAlVEYNeVEwjpUAJsGJ2iJwWOMUJzFu3ACJ1ylZBMvrKtjOr3/MXWe9zL7NnZk7Z87zkayde+fO3OOru8+ee8659ygiMDOz3nda1QUwM7POcOCbmWXCgW9mlgkHvplZJhz4ZmaZcOCbmWXCgW9mlgkHvplZJhz4ZmaZOL3qAkx29tlnx9KlS6suhplZUvbt2/eDiBicbbuuCvylS5cyPDxcdTHMzJIi6Xtz2c5NOmZmmXDgm5llwoFvZpYJB76ZWSYc+GZmmeiqUTpmZinauX+ErbsPcXRsnIUD/WxYtYw1K4aqLtavcOCbmZWwc/8Im3YcYPzYCQBGxsbZtOMAQNeFvgPfspdK7cy609bdh06G/YTxYyfYuvtQ151HDnzLWkq1M+tOR8fG57W+Su60tazNVDszm4uFA/3zWl8lB75lLaXamXWnDauW0b+g75R1/Qv62LBqWUUlmp4D37KWUu3MutOaFUNsXrucoYF+BAwN9LN57fKubBJ0G75lbcOqZae04UP31s6se61ZMdSVAT9Vy2r4kvok7Zd0V7F8nqS9kh6VdJukM1q1L7NWSal2ZlZWK2v47wEOAi8qlm8APhIRt0r6JHA18IkW7s+sJVKpnZmV1ZIavqRFwJuBG4tlAZcAdxSbbAfWtGJfZmbWnFY16XwUeB/wTLH8YmAsIo4Xy0cAV6HMzCpUOvAlXQE8FRH7Jq9usGlM8/n1koYlDY+OjpYtjpmZTaMVNfyVwJWSvgvcSr0p56PAgKSJPoJFwNFGH46IbRFRi4ja4OCsUzKamVmTSgd+RGyKiEURsRS4CtgTEe8A7gPeUmy2Driz7L7MzKx57bzx6lrgryUdpt6mf1Mb92VmZrNo6Y1XEXE/cH/x+jHgolZ+v5mZNc+PVjAzy4QD38wsEw58M7NMOPDNzDLhwDczy4QD38wsEw58M7NMOPDNzDLhwDczy4QD38wsEw58M7NMOPDNzDLhwDczy4QD38wsEw58M7NMOPDNzDLhwDczy0RLZ7wys+61c/8IW3cf4ujYOAsH+tmwahlrVgxVXSzrIAe+WQZ27h9h044DjB87AcDI2DibdhwAcOhnpHSTjqTnSnpA0n9LeljSB4v150naK+lRSbdJOqN8cc2sGVt3HzoZ9hPGj51g6+5DFZXIqtCKNvyngUsi4tXAhcDlki4GbgA+EhHnAz8Crm7BvsysCUfHxue13npT6cCPup8ViwuKfwFcAtxRrN8OrCm7LzNrzsKB/nmtt97UklE6kvokPQg8BdwD/C8wFhHHi02OAA0bCiWtlzQsaXh0dLQVxTGzKTasWkb/gr5T1vUv6GPDqmUVlciq0JLAj4gTEXEhsAi4CHhFo82m+ey2iKhFRG1wcLAVxTGzKdasGGLz2uUMDfQjYGign81rl7vDNjMtHaUTEWOS7gcuBgYknV7U8hcBR1u5LzObnzUrhhzwmWvFKJ1BSQPF637gMuAgcB/wlmKzdcCdZfdlZmbNa0UN/1xgu6Q+6n9Abo+IuyQ9Atwq6UPAfuCmFuzLzMyaVDrwI+JbwIoG6x+j3p5vZmZdwM/SMTPLhAPfzCwTDnwzs0w48M3MMuHANzPLhAPfzCwTDnwzs0w48M3MMuHANzPLhAPfzCwTDnwzs0x4EnOzTOzcP8LW3Yc4OjbOwoF+Nqxa5sclZ8aBb9nLIQh37h9h044DJycyHxkbZ9OOAwA993+16blJx7I2EYQjY+MEzwbhzv0jVRetpbbuPnQy7CeMHzvB1t2HKiqRVcGBb1nLJQiPjo3Pa731Jge+ZS2XIFw40D+v9dabHPiWtVyCcMOqZfQv6DtlXf+CPjasWlZRiawKDnzLWi5BuGbFEJvXLmdooB8BQwP9bF673B22mSk9SkfSYuAzwEuBZ4BtEfExSWcBtwFLge8CfxARPyq7P7NWmgi8Xh+lA/X/ay/+v2zuFBHlvkA6Fzg3Ir4p6YXAPmAN8CfADyNii6SNwJkRce1M31Wr1WJ4eLhUeczMciNpX0TUZtuudJNORDwZEd8sXv8UOAgMAauB7cVm26n/ETAzs4q0tA1f0lJgBbAXOCcinoT6HwXgJa3cl5mZzU/LAl/SC4DPAe+NiJ/M43PrJQ1LGh4dHW1VcczMbIqWBL6kBdTD/uaI2FGs/n7Rvj/Rzv9Uo89GxLaIqEVEbXBwsBXFMTOzBkoHviQBNwEHI+LDk97aBawrXq8D7iy7LzMza14rHp62Evgj4ICkB4t17we2ALdLuhp4HHhrC/ZlZmZNKh34EfFfgKZ5+9Ky329mZq3hO23NzDLhwDczy4QnQDEzKymVSXQc+GZmJaQ0m5ibdMzMSkhpEh0HvplZCSlNouPANzMrYeB5C+a1vkoOfDOzEqZ7wnzJJ8+3hQPfzKyEH48fm9f6KjnwzcxKSGleZAe+mVkJKc2L7HH4ZmYlpDQvsgPfekYVdzumcoeltVcqE8Q78DsspYBIraydvtsxpTsszaAH2vB37h9h5ZY9nLfxblZu2cPO/SNVF2laEwExMjZO8GxAdGOZUyorVHO3Y0p3WJpB4oHvUGqflMoK1dztmNIdlmaQeOA7lNonpbJCNUPjUhqOZwaJB75DqX1SKitUMzQupeF4ZaXUdGrTa0ngS/q0pKckPTRp3VmS7pH0aPHzzFbsazKHUvukVFaod5JuXrucoYF+BAwN9LN57fK2dp5Wsc8qpNZ0atNTtOCBD5JeB/wM+ExEvKpY94/ADyNii6SNwJkRce1M31Or1WJ4eHjO+506SgLqodTNv3SpjXxJpazWPiu37GGkwVXz0EA/X914SQUlsqkk7YuI2mzbtWRYZkT8h6SlU1avBl5fvN4O3A/MGPjzldINDxNSGa8LaZXV2ie1plObXjvH4Z8TEU8CRMSTkl7Sjp04lMzaa+FAf8Mafrc2ndr0Ku+0lbRe0rCk4dHR0aqLY2ZTpNafY9NrZ+B/X9K5AMXPpxptFBHbIqIWEbXBwcE2FsfMmpFL53QO2tmkswtYB2wpft7Zxn2ZWRu56bQ3tGpY5i3A14Blko5Iupp60L9B0qPAG4plMzOrSKtG6bx9mrcubcX3m5lZeZV32pqZWWc48M3MMuHANzPLhAPfzCwTDnwzs0w48M3MMuHANzPLhAPfzCwT7Xy0giXOz8M3a79O/p4lH/gOpfaYOrnMxCxHgI+vWYt0+vcs6SYdT73WPqlNEG+Wok7/niVdw5/pYLkWWo5nObJW8BX4zDr9e5Z04DuU2sezHPWeToevmwVn1+nfs6SbdKY7KA6l8srOcrRz/wgrt+zhvI13s3LLHjezVayK5k83C86u07OJJR34OU291ukALTPLkftWuk8V4esr8Nl1ejaxpJt0Jg5Kr7cRVnVp3OwsR+5b6T5VhK+bBeemk7OJJR34kMfUa6kFqGt23aeK8N2watkpFRXo3SvwVCTdpJOL1ALUfSvt1UzzXhXNn578vPskX8PPQWqXxlXV7D6w8wC37H2CExH0Sbz9txfzoTXL27rPTmu2ea+q5s8crsBT0vbAl3Q58DGgD7gxIjyZ+TyldmlcRbh8YOcB/v3rj59cPhFxcrmXQr9M857D19oa+JL6gI8DbwCOAN+QtCsiHmnnfntNip3TnQ6XW/Y+Me36Xgr81Jr3rLu0u4Z/EXA4Ih4DkHQrsBpoGPiPjf6ct33qa20uUroWndnPojPrzTi3PPA4tzzw+CyfyMeJiGnX99I5taDvNH554pmG63vp/2nt0e5O2yFgctXrSLHuJEnrJQ1LGj527Fibi2OWtsVn9XOaTl13murrzWbT7hq+Gqw7pSoWEduAbQC1Wi1u+4vXtrlI1oumtuFPeOfFS2Zt0knteS+pldfa7/Z3zW27dgf+EWDxpOVFwNE279MS1myYTYT6fEfppPi8F3e+WrMU07R9tuTLpdOB/wEuBUaAbwB/GBEPN9q+VqvF8PBw28pj3W1q+EJ9NFI7x26v3LKn4ZDXoYF+vrrxkrbs06zVJO2LiNps27W1DT8ijgPvBnYDB4Hbpwt7Mz/vxay92j4OPyK+AHyh3fux9Pl5L9YK7uOYnh+tYF2jikcy5PTE1Rz4Sa0zc+Bb1/DzXqwsP4N/Zn6Wjk2r05fGft6LleU+mZk58K2h1J7Bbwbuk5mNm3SsIV8aW4rcJzMz1/CtIV8aWyvk0iyYCge+NeRL4/bKYeigmwW7j5t0rKGqLo07PVl7FXIZOuhmwe7jGn6HpVKzq+LSOMXn2jQjtTmKm+Vmwe7jwO+g1AKt05fGDsLeCkI3C3YfN+l0kC9xZ5ZTEM5nfapyGjGTSlOkA7+Dcgm0ZjkIeysIc7mLOaU+GTfpdJAvcWdWdrJ29490nxxGzKTUFOnA76CygdasHILQ/SNWlZSu3B34HeSRL7NrNghTqmVZb0npyt2B32G5jHzp9FVFSrUs6y1VXbk3w4Hf46oIwiquKlKqZeUklebEMlLqk3Hg97gqgrCKq4qUalllpRKiqTUnlpFKn0ypYZmS3irpYUnPSKpNeW+TpMOSDklaVa6Y1qwqhgCWuapodjyzhwB23xBA33fSfcrW8B8C1gKfmrxS0gXAVcArgYXAvZJeHhEnfvUrrJ2quNxs9qqibI0wlVpWGSl1TrtfpfuUCvyIOAggaepbq4FbI+Jp4DuSDgMXAV8rsz9rTqeDsNnmlZTCrCophaj7VbpPu+60HQKemLR8pFhnGWi2eSWlMKtKSncj53JHcUpmreFLuhd4aYO3rouIO6f7WIN1Mc33rwfWAyxZsmS24lgimrmqcI1wdmU6p8t09jbz2ZRGr+Ri1sCPiMua+N4jwOJJy4uAo9N8/zZgG0CtVmv4R8Hy4DuRZ9dsiJbpHynz2Rz6VVLSrmGZu4DPSvow9U7b84EH2rQv6xG+E3lumgnRMv0j7lvpHaUCX9LvAf8CDAJ3S3owIlZFxMOSbgceAY4D13iEjs1FLncid1qZ/hH3rfSOUp22EfH5iFgUEc+JiHMiYtWk966PiF+PiGUR8cXyRTVrvVzCrExnb0odxTYzPw/fspZLmJUZMVPms6lMDJILP1rBspZLR3GZ/pEqOoqtPRTRPQNjarVaDA8PV10My0ynw3dqEEL9j0yvPQpi5ZY9DYfZDg3089WNl1RQot4laV9E1GbbzjV8y547itsjl/6RlLgN36zDcgnCXPpHUuLAN+uwXILQj1boPg58sw7LJQhzeWR1StyGb10lpcccgJ8xMxs/WqG7OPCta6Q2jM/PmLHUuEnHukZqMySlVl4zB751jdRGr6RWXjMHvnWN1EavpFZeMwe+dY3URq9UVV4/n8aa5U5b6xqpjV7x8/stNX6WjllC/Hwaa2Suz9Jxk45ZQtxRbGW4SccsIZ7ovb1Su/FvvlzDN0tIah3bKZnoHxkZGyd4tn+klzrFSwW+pK2Svi3pW5I+L2lg0nubJB2WdEjSqpm+x8zmxs+naZ8cbqQr26RzD7ApIo5LugHYBFwr6QLgKuCVwELgXkkv90TmZuX5sQztkUP/SNlJzL8cEceLxa8Di4rXq4FbI+LpiPgOcBi4qMy+zMzaKYcb6VrZhv9nwBeL10PAE5PeO1KsMzPrSjn0j8zapCPpXuClDd66LiLuLLa5DjgO3DzxsQbbNxzwL2k9sB5gyZIlcyiymXVar49egfRu/GvGrIEfEZfN9L6kdcAVwKXx7F1cR4DFkzZbBByd5vu3AdugfuPVHMps1jVyCMKc7u7t9f6RsqN0LgeuBa6MiF9MemsXcJWk50g6DzgfeKDMvsy6TQ7D+CCP0Su5KNuG/6/AC4F7JD0o6ZMAEfEwcDvwCPAl4BqP0LFek0sQ5jB6JRelhmVGxG/M8N71wPVlvt+sm+UShL67t3f4TluzJuUwjA/yGL2SCwe+WZNyCULf3ds7/PA0syblMIxvQq+PXsmFA9+sBAehpcRNOmZmmXDgm5llwoFvZpYJB76ZWSYc+GZmmXDgm5llwoFvZpYJB76ZWSYc+GZmmXDgm5llwoFvZpYJB76ZWSYc+GZmmfDTMs2sbXKY5D0lDnwza4uJSd4n5v2dmOQdcOhXpFSTjqR/kPStYgLzL0taWKyXpH+WdLh4/zWtKa6ZpSKXSd5TUrYNf2tE/FZEXAjcBfxdsf6NwPnFv/XAJ0rux8wSk8sk7ykpFfgR8ZNJi88Honi9GvhM1H0dGJB0bpl9mVlacpnkPSWlR+lIul7SE8A7eLaGPwQ8MWmzI8W6Rp9fL2lY0vDo6GjZ4phZl8hlkveUzBr4ku6V9FCDf6sBIuK6iFgM3Ay8e+JjDb4qGqwjIrZFRC0iaoODg83+P8ysy6xZMcTmtcsZGuhHwNBAP5vXLneHbYVmHaUTEZfN8bs+C9wN/D31Gv3iSe8tAo7Ou3RmljRP8t5dyo7SOX/S4pXAt4vXu4A/LkbrXAz8OCKeLLMvMzMrp+w4/C2SlgHPAN8D3lWs/wLwJuAw8AvgT0vux8zMSioV+BHx+9OsD+CaMt9tZmat5WfpmJllwoFvZpYJ1VtfuoOkUep9Aa12NvCDNnxvr/Fxmhsfp9n5GM1Nq47Tr0XErOPauyrw20XScETUqi5Ht/Nxmhsfp9n5GM1Np4+Tm3TMzDLhwDczy0Qugb+t6gIkwsdpbnycZudjNDcdPU5ZtOGbmVk+NXwzs+z1fOBLulzSoWL2rY1Vl6cbSFos6T5JByU9LOk9xfqzJN0j6dHi55lVl7UbSOqTtF/SXcXyeZL2FsfpNklnVF3GqkkakHSHpG8X59VrfT6dStJfFb9vD0m6RdJzO30u9XTgS+oDPk59Bq4LgLdLuqDaUnWF48DfRMQrgIuBa4rjshH4SkScD3ylWDZ4D3Bw0vINwEeK4/Qj4OpKStVdPgZ8KSJ+E3g19ePl86kgaQj4S6AWEa8C+oCr6PC51NOBD1wEHI6IxyLil8Ct1GfjylpEPBkR3yxe/5T6L+cQ9WOzvdhsO7CmmhJ2D0mLgDcDNxbLAi4B7ig2yf44SXoR8DrgJoCI+GVEjOHzaarTgX5JpwPPA56kw+dSrwf+nGfeypWkpcAKYC9wzsRjrIufL6muZF3jo8D7qD8RFuDFwFhEHC+WfU7By4BR4N+Kpq8bJT0fn08nRcQI8E/A49SD/sfAPjp8LvV64M955q0cSXoB8DngvVPmJzZA0hXAUxGxb/LqBpvmfk6dDrwG+ERErAB+TsbNN40U/RergfOAhdTnAH9jg03bei71euB75q1pSFpAPexvjogdxervT0w2X/x8qqrydYmVwJWSvku9OfAS6jX+geKyHHxOQf337EhE7C2W76D+B8Dn07MuA74TEaMRcQzYAfwOHT6Xej3wvwGcX/SEn0G9k2RXxWWqXNEOfRNwMCI+POmtXcC64vU64M5Ol62bRMSmiFgUEUupnzt7IuIdwH3AW4rNfJwi/g94opgMCeBS4BF8Pk32OHCxpOcVv38Tx6ij51LP33gl6U3Ua2V9wKcj4vqKi1Q5Sb8L/CdwgGfbpt9PvR3/dmAJ9RP0rRHxw0oK2WUkvR7424i4QtLLqNf4zwL2A++MiKerLF/VJF1IvWP7DOAx6rPcnYbPp5MkfRB4G/VRcvuBP6feZt+xc6nnA9/MzOp6vUnHzMwKDnwzs0w48M3MMuHANzPLhAPfzCwTDnwzs0w48M3MMuHANzPLxP8D9e+G42o4u3UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "lmod = smf.ols(formula='dist ~ speed', data=cars).fit()\n",
    "plt.scatter(lmod.fittedvalues,lmod.resid)\n",
    "plt.axhline(0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept   -17.579095\n",
       "speed         3.932409\n",
       "dtype: float64"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lmod.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iterate once using absolute residuals as response for linear relationship in the SD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-3.8198808759124097, Intercept   -11.180437\n",
       " speed         3.541541\n",
       " dtype: float64)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gamma = np.mean(np.abs(lmod.resid) - cars.speed)\n",
    "lmod = smf.wls(formula='dist ~ speed', data=cars, weights=np.sqrt(1/(gamma+cars.speed))).fit()\n",
    "gamma, lmod.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADmhJREFUeJzt3X2IHPd9x/HPp5IMQhEoRmdXlqVca4T/aF1qd/HFqBRRUB2bEquPRA1NHGhUik0SXEydxNShWDRUrf8IKTYyMbGJqzY07tWlLuqVJjg1scidrPjsCGGlONZJQrpE+Ck5Gvvy7R+3p5xWeze7sw9z+933C4Ru5+ZufwzDe2d/O3PjiBAAIJefq3oAAIDuI+4AkBBxB4CEiDsAJETcASAh4g4ACRF3AEiIuANAQsQdABJaW9UTb968OUZHR6t6egAYSFNTUz+IiJGi9SqL++joqCYnJ6t6egAYSLa/38p6TMsAQELEHQASIu4AkBBxB4CEiDsAJETcASAh4g4ACVV2nntZ949P69CRU5qP0Bpbe8e26cE9N1Q9LABYVQYq7vePT+srz7928fF8xMXHBB4AfmagpmUOHTnV1nIAGFYDFff5iLaWA8CwGqi4r7HbWg4Aw2qg4r53bFtbywFgWA3UB6qLH5pytgwArMxR0Xx1rVYL/uQvALTH9lRE1IrWG6hpGQBAa4g7ACRE3AEgIeIOAAkRdwBIiLgDQELEHQASIu4AkBBxB4CEiDsAJETcASAh4g4ACRF3AEiIuANAQsQdABIi7gCQEHEHgISIOwAkRNwBICHiDgAJEXcASIi4A0BCxB0AEiqMu+1ttr9u+7jtl21/ssk6tv0F2ydtv2j7pt4MFwDQirUtrPOupD+PiKO2N0qasj0REd9dss5tknbU/41Jerj+PwCgAoVH7hFxNiKO1r9+S9JxSVsbVrtD0hOx4HlJm2xv6fpoAQAtaWvO3faopBslHWn41lZJp5Y8ntHlLwAAgD5pOe623yPpa5I+FRFvNn67yY9Ek9+xz/ak7cnZ2dn2RgoAaFlLcbe9TgthfzIinmqyyoykbUseXyvpTONKEXEwImoRURsZGSkzXgBAC1o5W8aSviTpeEQ8tMxqT0v6SP2smfdLeiMiznZxnACANrRytsxOSX8sadr2sfqyz0jaLkkR8YikZyTdLumkpB9L+lj3hwoAaFVh3CPif9R8Tn3pOiHprm4NCgDQGa5QBYCEiDsAJETcASAh4g4ACRF3AEiIuANAQsQdABIi7gCQEHEHgISIOwAkRNwBICHiDgAJEXcASIi4A0BCxB0AEiLuAJAQcQeAhIg7ACRE3AEgIeIOAAkRdwBIiLgDQELEHQASIu4AkBBxB4CEiDsAJETcASAh4g4ACRF3AEiIuANAQsQdABIi7gCQEHEHgISIOwAkRNwBICHiDgAJEXcASKgw7rYfs33e9kvLfH+X7TdsH6v/+8vuDxMA0I61LazzZUlflPTECut8MyJ+uysjAgB0rPDIPSKelXShD2MBAHRJt+bcb7H9Hdv/YfuXuvQ7AQAltTItU+SopPdFxNu2b5c0LmlHsxVt75O0T5K2b9/ehacGADTT8ZF7RLwZEW/Xv35G0jrbm5dZ92BE1CKiNjIy0ulTAwCW0XHcbf+8bde/vrn+O3/Y6e8FAJRXOC1j+5CkXZI2256R9ICkdZIUEY9I+n1Jf2b7XUlzkj4UEdGzEQMAChXGPSL2Fnz/i1o4VRIAsEpwhSoAJETcASAh4g4ACRF3AEiIuANAQsQdABIi7gCQEHEHgISIOwAkRNwBICHiDgAJEXcASIi4A0BCxB0AEiLuAJAQcQeAhIg7ACRE3AEgIeIOAAkRdwBIiLgDQELEHQASIu4AkBBxB4CEiDsAJETcASAh4g4ACRF3AEiIuANAQsQdABIi7gCQEHEHgISIOwAkRNwBICHiDgAJEXcASIi4A0BChXG3/Zjt87ZfWub7tv0F2ydtv2j7pu4PEwDQjlaO3L8s6QMrfP82STvq//ZJerjzYQEAOlEY94h4VtKFFVa5Q9ITseB5SZtsb+nWAAEA7evGnPtWSaeWPJ6pL7uM7X22J21Pzs7OduGpAQDNdCPubrIsmq0YEQcjohYRtZGRkS48NQCgmW7EfUbStiWPr5V0pgu/FwBQUjfi/rSkj9TPmnm/pDci4mwXfi8AoKS1RSvYPiRpl6TNtmckPSBpnSRFxCOSnpF0u6STkn4s6WO9GiwAoDWFcY+IvQXfD0l3dW1EAICOcYUqACRE3AEgIeIOAAkRdwBIiLgDQELEHQASIu4AkBBxB4CEiDsAJETcASAh4g4ACRF3AEiIuANAQsQdABIi7gCQEHEHgISIOwAkRNwBICHiDgAJEXcASIi4A0BCxB0AEiLuAJAQcQeAhIg7ACRE3AEgIeIOAAkRdwBIiLgDQELEHQASIu4AkNDaqgeQ2f3j0zp05JTmI7TG1t6xbXpwzw1VDwvAECDuPXL/+LS+8vxrFx/PR1x8TOAB9BrTMj1y6MiptpYDQDcR9x6Zj2hrOQB0E3HvkTV2W8sBoJtairvtD9g+Yfuk7fuafP9O27O2j9X//Un3hzpY9o5ta2s5AHRT4QeqttdI+ntJuyXNSPq27acj4rsNq/5TRNzdgzEOpMUPTTlbBkAVWjlb5mZJJyPifyXJ9j9KukNSY9zR4ME9N5SK+fgLp3Xg8AmdeX1O12xar3tvvV57btzagxECyKqVaZmtkpae4jFTX9bo92y/aPufbTP3UNL4C6f16aemdfr1OYWk06/P6dNPTWv8hdNVDw3AAGkl7s0+AWw85ePfJI1GxK9I+i9Jjzf9RfY+25O2J2dnZ9sb6ZA4cPiE5t6Zv2TZ3DvzOnD4REUjAjCIWon7jKSlR+LXSjqzdIWI+GFE/F/94aOSfq3ZL4qIgxFRi4jayMhImfGmd+b1ubaWA0Azrcy5f1vSDtu/IOm0pA9J+qOlK9jeEhFn6w8/KOl4V0c5RK7ZtF6nm4T8mk3rW/p55usBSC0cuUfEu5LulnRYC9H+akS8bPuvbH+wvtonbL9s+zuSPiHpzl4NOLt7b71e69etuWTZ+nVrdO+t1xf+LPP1ABY5KrpislarxeTkZCXPvdqVPfre+fn/bnrUv3XTej1332/2YqgA+sz2VETUitbjD4etQntu3FpqKoX5egCLiHsinczXf/jRb+m57124+HjndVfqyY/f0tXxAegf/rZMImXn6xvDLknPfe+CPvzot7o+RgD9wZF7IotTOe3O1zeGvWh5I25KAqw+xD2ZsvP1ZXFTEmB1YloGHeGmJMDqxJE7tPO6K5tOwey87srCn+3kpiS7H/qGXjn/o4uPd1y1QRP37Cr8OQDFiDv05MdvKX22zBq7aciLbkrSGHZJeuX8j7T7oW+0FHheGICVEXdIUunTHveObbtkzn3p8pU0hr1o+VKdvjAAw4C4oyNV3JSkmy8MHPEjK+KOjpW9KUm/dXrEz4VeGCTEHZXYcdWGpkfaO67a0LPn7OSIf6ULvYoCz3UAqAJxRyUm7tlVeoqkiheGshd6dXodAO8WUBZxR2XKznV38sLQbytdB1AU907eLQzK9kHvEHcMpDKhquKIv5PrAMq+W+CzBUhcoYohMnHPrstC3uoR7XIXdBVd6LXc+f5F1wF0olefLRQZ2z+h0fv+/eK/sf0TrQ0YPcGRO4ZK2amJshd6lb0OoCpl3y2M7Z/Qubd+csmyc2/9RGP7J3Tks7u7Nr5G3FZyecQdaFGZqYlOrgPo5M9C9Ftj2IuWd8PibSXn3pmX9LPbSkpqKfBlXxgG5fMM4g70WNnrAMq+W6jis4UqHDh84mLYF829M68Dh08URrrsC0Onn2f0850GcQdWsTLvFjo5m2iQ3i10clvJsi8MnXye0ek7jXYRdyChfn+2cPXGK5pOwVy98YpS42hFJ7eVrOJ+w5280yiDuAO4RJl3C0c+u/uyD1Wv3nhFTz9MvffW6y85EpZau62k1NkLQ1n9fkEh7gC6opchb6bsbSWl8i8MnXye0e8XFEcLF1P0Qq1Wi8nJyUqeGwD6fbZM45y7tPCC8te/e0Nb0zK2pyKiVrgecQeA/ujG2TKtxp1pGQDok37ewJ4/PwAACRF3AEiIuANAQsQdABIi7gCQEHEHgIQqO8/d9qyk71fw1Jsl/aCC5x0UbJ9ibKOVsX2KdbKN3hcRI0UrVRb3qtiebOUCgGHF9inGNloZ26dYP7YR0zIAkBBxB4CEhjHuB6sewCrH9inGNloZ26dYz7fR0M25A8AwGMYjdwBIb6jibvtV29O2j9ke+r83bPsx2+dtv7Rk2ZW2J2y/Uv//vVWOsWrLbKPP2T5d34+O2b69yjFWyfY221+3fdz2y7Y/WV/OfqQVt0/P96Ghmpax/aqkWkRwDq4k278h6W1JT0TEL9eX/Y2kCxHxedv3SXpvRPxFleOs0jLb6HOS3o6Iv61ybKuB7S2StkTEUdsbJU1J2iPpTrEfrbR9/lA93oeG6sgdl4qIZyU13ur+DkmP179+XAs74tBaZhuhLiLORsTR+tdvSTouaavYjyStuH16btjiHpL+0/aU7X1VD2aVujoizkoLO6akqyoez2p1t+0X69M2Qznl0Mj2qKQbJR0R+9FlGraP1ON9aNjivjMibpJ0m6S76m+5gXY9LOk6Sb8q6aykv6t2ONWz/R5JX5P0qYh4s+rxrDZNtk/P96GhintEnKn/f17Sv0i6udoRrUrn6vOEi/OF5ysez6oTEeciYj4ifirpUQ35fmR7nRbC9WREPFVfzH5U12z79GMfGpq4295Q/0BDtjdI+i1JL638U0PpaUkfrX/9UUn/WuFYVqXFaNX9joZ4P7JtSV+SdDwiHlryLfYjLb99+rEPDc3ZMrZ/UQtH69LCjcH/ISL2Vzikytk+JGmXFv5C3TlJD0gal/RVSdslvSbpDyJiaD9QXGYb7dLC2+mQ9KqkP12cXx42tn9d0jclTUv6aX3xZ7Qwrzz0+9EK22everwPDU3cAWCYDM20DAAME+IOAAkRdwBIiLgDQELEHQASIu4AkBBxB4CEiDsAJPT/8gm+Fgi8lfQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "weights=np.sqrt(1/(gamma+cars.speed))\n",
    "plt.scatter(cars.speed, weights)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing for Lack of fit \n",
    "\n",
    "Read in the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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>Fe</th>\n",
       "      <th>loss</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.01</td>\n",
       "      <td>127.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.48</td>\n",
       "      <td>124.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.71</td>\n",
       "      <td>110.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.95</td>\n",
       "      <td>103.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.19</td>\n",
       "      <td>101.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     Fe   loss\n",
       "0  0.01  127.6\n",
       "1  0.48  124.0\n",
       "2  0.71  110.8\n",
       "3  0.95  103.9\n",
       "4  1.19  101.5"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corrosion = pd.read_csv(\"data/corrosion.csv\")\n",
    "corrosion.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda/lib/python3.7/site-packages/scipy/stats/stats.py:1394: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=13\n",
      "  \"anyway, n=%i\" % int(n))\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>loss</td>       <th>  R-squared:         </th> <td>   0.970</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.967</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   352.3</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 25 Sep 2018</td> <th>  Prob (F-statistic):</th> <td>1.06e-09</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:55:20</td>     <th>  Log-Likelihood:    </th> <td> -31.890</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    13</td>      <th>  AIC:               </th> <td>   67.78</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    11</td>      <th>  BIC:               </th> <td>   68.91</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</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>  129.7866</td> <td>    1.403</td> <td>   92.524</td> <td> 0.000</td> <td>  126.699</td> <td>  132.874</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Fe</th>        <td>  -24.0199</td> <td>    1.280</td> <td>  -18.769</td> <td> 0.000</td> <td>  -26.837</td> <td>  -21.203</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 1.170</td> <th>  Durbin-Watson:     </th> <td>   2.535</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.557</td> <th>  Jarque-Bera (JB):  </th> <td>   0.958</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.551</td> <th>  Prob(JB):          </th> <td>   0.619</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.255</td> <th>  Cond. No.          </th> <td>    2.99</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:                   loss   R-squared:                       0.970\n",
       "Model:                            OLS   Adj. R-squared:                  0.967\n",
       "Method:                 Least Squares   F-statistic:                     352.3\n",
       "Date:                Tue, 25 Sep 2018   Prob (F-statistic):           1.06e-09\n",
       "Time:                        15:55:20   Log-Likelihood:                -31.890\n",
       "No. Observations:                  13   AIC:                             67.78\n",
       "Df Residuals:                      11   BIC:                             68.91\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept    129.7866      1.403     92.524      0.000     126.699     132.874\n",
       "Fe           -24.0199      1.280    -18.769      0.000     -26.837     -21.203\n",
       "==============================================================================\n",
       "Omnibus:                        1.170   Durbin-Watson:                   2.535\n",
       "Prob(Omnibus):                  0.557   Jarque-Bera (JB):                0.958\n",
       "Skew:                           0.551   Prob(JB):                        0.619\n",
       "Kurtosis:                       2.255   Cond. No.                         2.99\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lmod = smf.ols(formula='loss ~ Fe', data=corrosion).fit()\n",
    "lmod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c1bec50f0>]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4lFXax/HvnQKEGqpCAAHpSNOIQCBYgGADRFTUtaJYQNouKrtrW30V5X0pFlTsvWNAKQHEJYCABiK9g5TQQglFQj/vHxl3ASeQQGaeSfL7XFeuZM48M7n32ZFfnnLuY845REREThXmdQEiIhKaFBAiIuKXAkJERPxSQIiIiF8KCBER8UsBISIifikgRETELwWEiIj4pYAQERG/Irwu4FxUqFDB1ahRw+syRETylXnz5u1wzlU803b5OiBq1KhBSkqK12WIiOQrZrY+J9vpFJOIiPilgBAREb8UECIi4pcCQkRE/FJAiIiIXwoIERHxSwEhIiJ+FcqA2Ln/EP/6bil7Dx7xuhQRkZBVKANi1pqdvP/TOjoMm87Updu8LkdEJCQVyoDo3LQK3z4cR9niRbjvwxT6fpbKzv2HvC5LRCSkFMqAAGhaLZpxfdowoH1dJi7eQvth0xn7axrOOa9LExEJCYU2IACKRITRr30dxvdtywXlS9Dv81/p+UEKmzMyvS5NRMRzhTog/lD3vFJ881Br/nltA35as4OOw5P5ZO56jh/X0YSIFF4BCwgze9fMtpvZ4hPGnjWzhWb2q5lNNrMqvnEzs5fNbLXv+YsDVVd2wsOM+9rWYnL/djSpWoZ/fLuYW9+aw7odv5+0XWJqGnFDplHz8fHEDZlGYmpasEsVEQmKQB5BvA90OmVsqHOuiXOuGfA98KRv/Gqgju+rF/B6AOs6rerli/PJfZfx4o2NWbplL51GJDM6eQ1Hjx0nMTWNwWMWkZaRiQPSMjIZPGaRQkJECqSArQfhnEs2sxqnjO094WEJ4I9zOF2AD13WFeI5ZhZtZpWdc1sCVd/pmBm3XFqdy+tV4p+Ji3l+wnK+X7iFrXsOknnk2EnbZh45xtCkFXRtHuNFqSIiARP0axBm9j9mthG4nf8eQcQAG0/YbJNvzN/re5lZipmlpKenB7TW80oXY/Qdl/Dqbc1J253J9n3+b4XVRW0RKYiCHhDOuX8456oBnwB9fMPmb9NsXj/aORfrnIutWPGMK+adMzPjuiZVmDqwHVGR4X63qRIdFfA6RESCzcu7mD4FbvT9vAmodsJzVYHNQa/oNMqWKMIL3RpTJPzkXVYsIoxBCfU8qkpEJHCCGhBmVueEh52B5b6fxwF3+u5magns8er6w+l0bR7DS92bULlMsf+MlSgaQcVSRT2sSkQkMAJ2kdrMPgMuByqY2SbgKeAaM6sHHAfWAw/6Np8AXAOsBg4A9wSqrnPVtXnMfy5Iz127k8fHLOL2t+dyS2w1/n5tA8pERXpcoYhI3gjkXUy3+hl+J5ttHdA7ULXkpcTUNIYmrWBzRiZVoqPod1Ud1u74nbdmrOXHFdt5tutFJDQ63+syRUTOmWZS50JiahqDvl5w0jyIfyQuov75pUh8OI7yJYvywEfz6P3JfNKzueNJRCS/UEDkwjPfLeHIsZNvrjpyzPHMd0toXLUM4/rE8beOdZmydBsdhk9nzPxNav4nIvmWAiIXdh/wv8DQH+OR4WH0ubIOE/q1oVaFEgz8cgH3vP8LaZonISL5kAIiAGpXKsVXD7bmqesbMnftLjoOm85Hs39T8z8RyVcUELkQnc0dSv7Gw8OMe+JqMnlAPBdfUJYnxi6hx+g5rE3fH+gyRUTyhAIiF57u3IjIsJMnfUeGGU93bpTta6qVK86H97ZgaPcmLN+6l04jZ/D6v7Oa/4mIhDIFRC50bR7D0JuaEhMdhQEx0VEMvanpGRv1mRk3xVZj6sB2XFGvIi9OWk7XUbNYsnlPcAoXETkLlp/vsomNjXUpKSlel5FrExdt4YmxS9h94DAPtqvFI1fWoVg2fZ5ERPKamc1zzsWeaTsdQXjg6saVmTownq7NYnjtxzVc+/IM5q3f5XVZIiInUUB4JLp4Ef7v5qZ8cG8LDh45Tvc3ZvP0uCX8fuio16WJiAAKCM+1q1uRpAHx3NnyAj6Y/RsdhyeTvDKw61yIiOSEAiIElCwawTNdLuLLB1pRNDKMO9/9mb99tYCMA4e9Lk1ECjEFRAi5tEY5JvRty8OXX8i3qWm0H5bMxEUh1/VcRAoJBUSIKRYZzqOd6jO2dxyVShXloU/m89DH89i+76DXpYlIIROwdt9ybi6KKcPYPnGMTl7LyB9W8dOanfzz2gZ0v6QqZv5WaA2uU9ueD0qod8b5ICKSv+gIIoRFhofR+4raTOjbljqVSjLo64Xc+e7PbNx1wNO6ElPTGDxm0UltzwePWURiapqndYlI3lJA5AO1K5Xkywda8a8ujZi/fjcJI5J5f9Y6z5r/DU1aQeaRYyeNZR45xtCkFZ7UIyKBoYDIJ8LCjDtb1SBpQDyxNcrx9HdLufnN2azeHvzmf5uzaV+e3biI5E8KiHymatnifHDPpfzfTU1ZtX0/14ycwWs/ruZIEJv/VYmOytW4iORPCoh8yMy48ZKqTB3YjvYNKzE0aQVdXp3F4rTgNP8blFCPqFN6R0VFhjMooV5Qfr+IBIcCIh+rWKooo26/hDf+cgnp+w/R5bVZvDhpOQdPuT6Q17o2j+GFbo1P6mr7QrfGuotJpIBRN9cCYs+BI/zPhKV8mbKJWhVK8GL3Jlxao5zXZYlICFI310KmTPFIXurelI97XsbhY8e56Y3ZPDl2MfsD1PwvMTWNuCHTqPn4eOKGTNMtriIFkAKigGlTpwJJ/eO5J64GH81ZT8dh0/lxxfY8/R2aByFSOCggCqASRSN46vpGfP1ga4oXjeCe935h4Be/svv3vGn+p3kQIoWDAqIAu+SCsozv24ZHrqzNuAWb6TB8OuMXbuFcrztpHoRI4aCAKOCKRoTz1471GNenDZXLRNH70/k88NE8tu89++Z/mgchUjgoIAqJhlVK8+3DrRl8dX2mr0znqmHT+fKXjWd1NKF5ECKFgwKiEIkID+OBdhcysV9bGlQuzaPfLOSOd35mw87cNf/TPAiRwiFg8yDM7F3gOmC7c+4i39hQ4HrgMLAGuMc5l+F7bjDQEzgG9HXOJZ3pd2gexNk7ftzx6c8bGDJxOceOO/6WUI+7W9cgPMz7VuIiElihMA/ifaDTKWNTgIucc02AlcBgADNrCPQAGvleM8rMwpGACQsz/tLyAiYPiOeyWuV49vuldH/jJ1Zt2+d1aSISIgIWEM65ZGDXKWOTnXN/zNyaA1T1/dwF+Nw5d8g5tw5YDbQIVG3yX1Wio3jv7ksZcUszftvxO9e+PJOXf1jF4aPBa/4nIqHJy2sQ9wITfT/HABtPeG6Tb+xPzKyXmaWYWUp6enqASywczIyuzWOYMrAdCRedz7ApK+n86kwWbsrwujQR8ZAnAWFm/wCOAp/8MeRnM78XR5xzo51zsc652IoVKwaqxEKpQsmivHJrc966M5bdBw7T9bVZvDBhWcCb/4lIaAp6QJjZXWRdvL7d/fcK+Sag2gmbVQU2B7s2ydKh4XlMHtCOWy6txpvJa+k0Ipk5a3d6XZaIBFlQA8LMOgGPAZ2dcyfeWzkO6GFmRc2sJlAH+DmYtcnJykRF8kK3Jnx632Ucd9Bj9Bz+8e0i9h084nVpIhIkAQsIM/sMmA3UM7NNZtYTeBUoBUwxs1/N7A0A59wS4EtgKTAJ6O2c03mNENC6dgUm9W/LfW1q8tnPG+g4PJlpy7d5XZaIBIHWg5AcS92wm8e+WcjKbfvp2qwKT17fiHIlinhdlojkUijMg5ACpnn1snz/SFv6XVWH8Yu20H7YdMYt2HzOzf9EJDQpICRXikSEMaBDXb57pA3VykbR97NU7v9wHlv3nH3zPxEJTQoIOSv1zy/NmIfj+Mc1DZi5Op0Ow6bz2c8bdDQhUoAoIOSshYcZ98fXYlK/eBrFlGbwmEXc9tZc1u/83evSRCQPKCDknNWoUIJP72vJ8zc0ZnHaHhJGJPP2jLUcO66jCZH8TAEheSIszLjtsupMHhhP3IUVeG78Mrq9/hMrtqr5n0h+pYCQPFW5TBRv3xXLy7c2Z+OuA1z3ygxGTF2p5n8i+ZACQvKcmdG5aRWmDmzHNY0rM2LqKq5/ZSa/blTzP5H8RAEhAVOuRBFG9mjOO3fFsifzCN1GzeK575eSeViT5EXyAwWEBNxVDc5j8sB4erSoztsz15EwIpmf1uzwuiwROQMFhARF6WKRPH9DYz67vyVhBre9NZfBYxayV83/REKWAkKCqtWF5ZnYL54H4mvxxS8b6TBsOlOXqvmfSChSQEjQRRUJZ/A1DUjsHUfZ4kW478MUHvkslZ37D3ldmoicQAEhnmlSNZpxfdowsENdJi3Oav6XmJqmdh0iIUIBIZ4qEhFG36vqML5vWy4oX4L+X/xKzw9S2JyR6XVpIoWeAkJCQt3zSvHNQ6154rqGzF6zk47Dk/l4znqOq12HiGcUEBIywsOMnm1qktQ/nqbVyvDPxMXc+tYc1u1Q8z8RLyggJORUL1+cj3texks3NmHplr10GpHMm9PXcPSY2nWIBFOE1wWI+GNm3HxpNdrVq8g/ExfzwsTljF+0hRdvbEKDyqVJTE1jaNIKNmdkUiU6ikEJ9ejaPMbrskUKFK1JLSHPOceERVt5atxiMg4c4cr6lUhemc7BExoARkWG80K3xgoJkRzQmtRSYJgZ1zapzJQB7ejctAqTl247KRwAMo8cY2jSCo8qFCmYFBCSb5QtUYRhtzTL9nndGiuStxQQku/EREf5Ha+SzbiInB0FhOQ7gxLqERUZ/qfxqmWj2JOp5n8ieUUBIflO1+YxvNCtMTHRURhQpUwxrmpQiZT1u+kwbDpJS7Z6XaJIgaC7mKTAWLRpD49+s5BlW/ZybePKPN25ERVLFfW6LJGQk9O7mDQPQvKl7OZBjOsTx+jktYycuoqZq3fw5HUN6XZxDGbmdcki+Y5OMUm+k5iaxuAxi0jLyMQBaRmZDB6ziMTUNCLDw+h9RW0m9GtD7Uol+etXC7j7vV9I0x1OIrmmgJB8Z2jSCjKPnLyu9anzIGpXKsVXD7Ti6esb8stvu+g4bDofzv5Nzf9EciFgAWFm75rZdjNbfMLYTWa2xMyOm1nsKdsPNrPVZrbCzBICVZfkf9nNdzh1PCzMuDsuq/nfxReU5cmxS7hl9GzWpO8PRpki+V4gjyDeBzqdMrYY6AYknzhoZg2BHkAj32tGmdmf72MUIfv5DtmNVytXnA/vbcHQ7k1YsXUfV4+cwah/r+aImv+JnFbAAsI5lwzsOmVsmXPOXz+ELsDnzrlDzrl1wGqgRaBqk/zN3zyIqMhwBiXUy/Y1ZsZNsdWY+td2XFmvEi9NWkHX12axOG1PoMsVybdC5RpEDLDxhMebfGN/Yma9zCzFzFLS09ODUpyEllPnQcRER+W4UV+lUsV4445LeP32i9m29xBdXpvF0KTlHDzlmoaIhM5trv7uQfR7NdE5NxoYDVnzIAJZlISurs1jzqlz69WNK9PqwvI8N34Zr/24homLt/LSjU2IrVEuD6sUyd9C5QhiE1DthMdVgc0e1SKFRHTxIvzvTU358N4WHDpynJvenM3T45bw+6GjXpcmEhJCJSDGAT3MrKiZ1QTqAD97XJMUEvF1KzJ5QDx3tarBB7N/o+PwZJJX6vSlSCBvc/0MmA3UM7NNZtbTzG4ws01AK2C8mSUBOOeWAF8CS4FJQG/nnE4KS9CUKBrB050b8dUDrSgaGcad7/7M375aQMaBw16XJuIZ9WISOcXBI8d4Zdoq3pi+lrLFi/Bsl0Zc3biy12WJ5BmtKCdylopFhjMooT7j+sRxXumiPPTJfB78aB7b9x70ujSRoFJAiGSjUZUyjO0dx2Od6jNtxXbaD5vOVykbyc9H3SK5kaOAMLN+ZlbasrxjZvPNrGOgixPxWkR4GA9dfiET+7Wl3vmlGPT1Qu5892c27jrgdWkiAZfTI4h7nXN7gY5AReAeYEjAqhIJMRdWLMkXvVrxbJdGzF+/m4QRybw/a52a/0mBltOA+GMi2zXAe865Bfif3CZSYIWFGXe0qkHSgHgurVGOp79byk1vzmb19n1elyYSEDkNiHlmNpmsgEgys1KAOp1JoVS1bHHev+dSht3clDXp+7lm5ExenbZKzf+kwMnRba5mFgY0A9Y65zLMrBxQ1Tm3MNAFno5ucxWvpe87xNPfLWH8wi00qFyaod2bcFFMGa/LEjmtvL7NtRWwwhcOfwH+CagNphR6FUsV5bXbLubNOy5hx/6s5n9DJqr5nxQMOQ2I14EDZtYUeBRYD3wYsKpE8pmERuczdUA7ul9clTemr+GakTP4ed2uM79QJITlNCCOuqxzUV2Akc65kUCpwJUlkv+UKR7Ji92b8HHPyzh87Dg3vzmbJxIXs+/gEa9LEzkrOQ2IfWY2GLiDrB5K4UBk4MoSyb/a1KnA5AHx3BtXk4/nridheDI/rtjudVkiuZbTgLgFOETWfIitZC3mMzRgVYnkc8WLRPDk9Q35+sHWlCgawT3v/cLAL35l9+9q/if5R44CwhcKnwBlzOw64KBzTtcgRM7gkgvK8n3fNvS9sjbjFmym/bDpfL9ws9p1SL6Q01YbN5O1PsNNwM3AXDPrHsjCRAqKohHhDOxYj+8eaUOV6Cj6fJrKAx/NY5ua/0mIy+k8iAVAB+fcdt/jisBU51zTANd3WpoHIfnN0WPHeWfmOoZNWUmRiDD+eW0Dbo6thpkaE0jw5PU8iLA/wsFnZy5eK1JgJaamETdkGjUfH0/ckGkkpqaddvuI8DAeaHchk/rH06ByaR77ZhF/eWcuG3aq+Z+Enpz+Iz/JzJLM7G4zuxsYD0wIXFkioS8xNY3BYxaRlpGJA9IyMhk8ZtEZQwKgZoUSfH5/S57rehELNu4hYUQy78xcxzE1/5MQktOL1IOA0UAToCkw2jn3WCALEwl1Q5NWkHnKjOnMI8cYmrQiR68PCzP+0vICJg+Ip9WF5Xn2+6Xc+PpPrNym5n8SGnJ8msg5941zbqBzboBz7ttAFiWSH2zOyMzVeHaqREfxzl2xjOzRjPU7f+fal2fw8g+rOHxUzf/EW6cNCDPbZ2Z7/XztM7O9wSpSJBRViY7K1fjpmBldmsUwdWA7Ol1UmWFTVtL51Zks2JhxrmWKnLXTBoRzrpRzrrSfr1LOudLBKlIkFA1KqEdUZPhJY1GR4QxKqHfW71m+ZFFeubU5b90Zy+4Dh7lh1CxemLCMzMNq/ifBpzuRRM5S1+YxvNCtMTHRURgQEx3FC90a07V5zDm/d4eG5zFlYDtuubQabyav5eqRycxes/PcixbJhRzNgwhVmgchhcFPq3fw+JhFbNh1gNsuq87jV9endDG1QpOzl9fzIETEI61rVyCpfzz3t63J5z9voOOwZKYt3+Z1WVII6AhCJB/5dWMGj329kBXb9tGlWRWevK4hM1btYGjSCjZnZFIlOopBCfXy5DSXFFw5PYJQQIjkM4ePHmfUv1fz2o+rKRIRxuGjxzly7L//HUdFhufZtRApmHSKSaSAKhIRRv/2dfn+kbZ/CgfI3WQ9kdNRQIjkU/XOL/WncPhDbifrifijgBDJx2KymZRXqVTRIFciBVHAAsLM3jWz7Wa2+ISxcmY2xcxW+b6X9Y2bmb1sZqvNbKGZXRyoukQKEn+T9QB2HTjMW8lr1fxPzkkgjyDeBzqdMvY48INzrg7wg+8xwNVAHd9XL+D1ANYlUmD4m6z39PUNaVe3Iv8zYRndRs1ixVY1/5OzE9C7mMysBvC9c+4i3+MVwOXOuS1mVhn4t3Ounpm96fv5s1O3O9376y4mEf+cc3y/cAtPj1vC3oNHePjy2vS+ojZFInRWWUL3Lqbz/vhH3/e9km88Bth4wnabfGN/Yma9zCzFzFLS09MDWqxIfmVmXN+0ClMGtuPaxpUZ+cMqrntlBqkbdntdmuQjofLnhL/1Fv0e2jjnRjvnYp1zsRUrVgxwWSL5W7kSRRjRoznv3h3LvoNH6fb6Tzz7/VIOHD7qdWmSDwQ7ILb5Ti3h+/7HMqabgGonbFcV2Bzk2kQKrCvrn8fkAfHcfll13pm5jk4jZvDT6h1elyUhLtgBMQ64y/fzXcDYE8bv9N3N1BLYc6brDyKSO6WKRfJc18Z83qslYQa3vT2Xx79ZyJ7MI16XJiEqkLe5fgbMBuqZ2SYz6wkMATqY2Sqgg+8xZK1vvRZYDbwFPByoukQKu5a1yjOpfzwPtKvFlykb6Th8OlOWqvmf/Jl6MYkUYgs3ZfDo1wtZvnUf1zWpzNOdG1GhpCbZFXSheheTiISQJlWjGdenDX/tUJfJS7bRYdh0ElPTyM9/OEreUUCIFHJFIsJ45Ko6jO/bhhoVStD/i1+59/1f1M9JFBAikqXOeaX4+sHWPHldQ+as3UXH4cl8NGc9x9Wuo9BSQIjIf4SHGfe2qcnkAfE0qxbNE4mL6fHWHNbt+N3r0sQDCggR+ZNq5YrzUc8WvHRjE5Zt2UunEcm8MX0NR48d97o0CSIFhIj4ZWbcfGk1pg5sR7u6FRkycTk3jPqJpZv3el2aBIkCQkRO67zSxXjzjkt47baL2bInk86vzuT/Jq/g0NFjXpcmAaaAEJEzMjOubVKZKQPa0blZFV6ZtpprX57JvPVq/leQKSBEJMfKlijCsJub8f49l5J5+Bjd3/iJZ75bwu+H1PyvIFJAiEiuXV6vEkkD4rmj5QW8N+s3EkYkM2OV2u8XNAoIETkrJYtG8K8uF/HlA60oEh7GHe/8zKNfL2DPATX/KygUECJyTlrULMeEfm156PIL+WZ+Gu2HT2fS4q1elyV5QAEhIuesWGQ4j3Wqz9jecVQsWZQHP55H70/mk77vkNelyTlQQIhInrkopgxj+8QxKKEeU5Zto/2w6Xwzb5Oa/+VTCggRyVOR4WH0vqI2E/q2pXalkvz1qwXc9d4vbNp9wOvSJJcUECISELUrleSrB1rxTOdGpPy2i4ThyXw4+zc1/8tHFBAiEjBhYcZdrWuQ1D+eiy8oy5Njl3DL6NmsSd/vdWmSAwoIEQm4auWK8+G9Lfjfm5qyctt+rh45g1H/Xs0RNf8LaQoIEQkKM6P7JVWZMjCe9g0q8dKkFXR9bRaL0/Z4XZpkQwEhIkFVqVQxRt1+CW/85WK27T1El9dm8dKk5Rw8ouZ/oUYBISKe6HRRZX4Y2I5uzWMY9e81XPPyDFJ+2+V1WXICBYSI5KnE1DTihkyj5uPjiRsyjcTUtGy3LVM8kqE3NeXDe1tw6MhxbnpzNk+NXcx+Nf8LCQoIEckzialpDB6ziLSMTByQlpHJ4DGLThsSAPF1KzJ5QDx3tarBh3PWkzA8mekr1fzPawoIEckzQ5NWkHnKtYTMI8cYmrTijK8tUTSCpzs34usHW1EsMoy73v2Zv365gIwDhwNVrpyBAkJE8szmjMxcjftzyQXlGN+3LX2uqM3YX9NoP2w6ExZtyasSJRcUECKSZ8pEReZqPDvFIsP5W0I9xvaJ4/wyxXj4k/k8+NE8tu89mBdlSg4pIEQkz5jlbvxMGlUpQ+LDcTzWqT7TVmyn/bDpfJmyUc3/gkQBISJ5JiObxYKyG8+JiPAwHrr8Qib1a0v980vz6NcLufPdn9m4S83/Ak0BISJ5pkp0VK7Gc6NWxZJ83qslz3ZpxPz1u0kYkcx7s9ZxTM3/AsaTgDCzfma22MyWmFl/31g5M5tiZqt838t6UZuInL1BCfWIigw/aSwqMpxBCfXy5P3Dwow7WtVg8sB2tKhZjme+W8pNb/zE6u378uT95WRBDwgzuwi4H2gBNAWuM7M6wOPAD865OsAPvsciko90bR7DC90aExMdhQEx0VG80K0xXZvH5OnviYmO4r27L2X4LU1Zu+N3rhk5k1enrVLzvzxmwb7YY2Y3AQnOuft8j58ADgE9gcudc1vMrDLwb+fcaf/siI2NdSkpKQGvWURC1479h3hq3BLGL9xC/fNLMbR7UxpXLeN1WSHNzOY552LPtJ0Xp5gWA/FmVt7MigPXANWA85xzWwB83yv5e7GZ9TKzFDNLSU/XTEuRwq5CyaK8dtvFvHnHJez6/TBdR81iyEQ1/8sLQT+CADCznkBvYD+wFMgE7nHORZ+wzW7n3GmvQ+gIQkROtCfzCM+PX8YXKRupWaEEQ7o15rJa5b0uK+SE8hEEzrl3nHMXO+figV3AKmCb79QSvu/bvahNRPKvMlGRvNi9CZ/cdxlHjx/nltFzeCJxMfsOnv1ttoWZV3cxVfJ9rw50Az4DxgF3+Ta5CxjrRW0ikv/F1a5AUv94erapycdzs5r//bhcf3PmllfzIL4xs6XAd0Bv59xuYAjQwcxWAR18j0VEzkrxIhE8cV1DvnmoNSWKRnDP+78w4Itf2fW7mv/llCfXIPKKrkGISE4cOnqM135cw6gfV1MmKpJnujTi2saVsbPtAZLPhfQ1CBGRYCoaEc7ADnX57pE2xJSNos+nqfT6aB7b1PzvtBQQIlJoNKhcmjEPtebv19QneWU67YdN54tfNqj5XzYUECJSqESEh9Er/kKS+sfTsHJpHvtmEbe/PZcNO9X871QKCBEplGpUKMFn97fk+Rsas3DTHjqOmM7bM9aq+d8JFBAiUmiFhRm3XVadKQPjaX1hBZ4bv4wbX/+JldvU/A8UECIiVC4TxTt3xTKyRzM27DrAtS/PYOTUVRw+Wrib/0V4XYCISCgwM7o0i6FN7Qo8891Shk9dycTFW3jxxiY0rRZ95jfIhcTUNIYmrWBzRiZVoqMYlFAvzzve5gUdQYiInKB8yaK8fGtz3r4zlowDR7hh1Cyen7CMzMN50/wvMTWNwWMWkZaRiQPSMjIZPGYRialpefL+eUkBISLiR/uG5zF5YDw9WlRndPJaOo1MZvaanef8vkOTVpB5SqfZzCPHGJq04pzfO68pIEREslG6WCTP39CYT++/DIBb35rD4DFSboqvAAAIgklEQVSL2HsOzf82Z2TmatxLCggRkTNofWEFJvWLp1d8Lb74ZQMdhyXzw7JtZ/VegVy3O68pIEREciCqSDh/v6YBYx6Oo0xUJD0/SKHvZ6ns3H8oV+8T6HW785ICQkQkF5pVi+a7R9owoH1dJi7eQofhyYz9NS3H7TqCtW53XlA3VxGRs7Ry2z4e/Xohv27M4Kr6lXjuhouoXCb0ThWdSt1cRUQCrO55pfjmodb889oGzFqzg47Dkvl07gaOF5B2HQoIEZFzEB5m3Ne2FpP7t6Nx1TL8/dtF3Pb2HH7b8Xu2r0lMTSNuyDRqPj6euCHTQnIOBCggRETyRPXyxfnkvssY0q0xS9L2kjAimdHJazh67OR2HZooJyJSCJkZPVpUZ8rAdrStU5HnJyznxtd/YvnWvf/ZRhPlREQKsfPLFOOtOy/hlVubs2l3Jte9PJNhU1Zy6OixfDVRTs36REQCwMy4vmkV4mpX4Nnvl/LyD6uYtHgLFUoWJd3P3AlNlBMRKWTKlSjC8Fua8d7dl7Lv4FF27D9EeJidtI0myomIFGJX1K/E5AHx3N6yOseOu/+ERChPlFNAiIgESalikTzXtTFf9GpJ9XLFAWhbpwJX1K/kcWX+KSBERILsslrlmdivLQ+0q8WXKRvpMGw6k5ds9bqsP1FAiIh4oFhkOIOvbkBi7zjKlShCr4/m0efT+ezIZfO/QFJAiIh4qEnVrOZ/f+tYl8lLttF+2HS+Td2U4+Z/gaSAEBHxWGR4GH2urMOEfm2oVaEEA75YwD3v/0Kax3MjFBAiIiGidqVSfPVga566viFz1+6i47DpfDRnvWfN/xQQIiIhJDzMuCeuJpMHxNO8elmeSFxMj9FzWJu+P+i1eBIQZjbAzJaY2WIz+8zMiplZTTOba2arzOwLMyviRW0iIqGgWrnifNSzBS91b8LyrXu5euQM3pj+5+Z/gRT0gDCzGKAvEOucuwgIB3oALwLDnXN1gN1Az2DXJiISSsyMm2OrMXVgOy6vV5EhE5fTddQslm7ee+YX5wGvTjFFAFFmFgEUB7YAVwJf+57/AOjqUW0iIiGlUulivHlHLK/ffjFb9xyi86szeWfmuoD/3qAHhHMuDfhfYANZwbAHmAdkOOeO+jbbBPidd25mvcwsxcxS0tPTg1GyiEhIuLpxZaYOjKdLsxgu8M3EDiQvTjGVBboANYEqQAngaj+b+r1s75wb7ZyLdc7FVqxYMXCFioiEoOjiRfi/m5vSvuF5Af9dXpxiag+sc86lO+eOAGOA1kC075QTQFVgswe1iYiIjxcBsQFoaWbFzcyAq4ClwI9Ad982dwFjPahNRER8vLgGMZesi9HzgUW+GkYDjwEDzWw1UB54J9i1iYjIf3myopxz7ingqVOG1wItPChHRET80ExqERHxSwEhIiJ+KSBERMQvBYSIiPhlobAoxdkys3RgfZB/bQVgR5B/Z36g/ZI97Rv/tF/8C8Z+ucA5d8aZxvk6ILxgZinOuViv6wg12i/Z077xT/vFv1DaLzrFJCIifikgRETELwVE7o32uoAQpf2SPe0b/7Rf/AuZ/aJrECIi4peOIERExC8FRDbMrJOZrTCz1Wb2uJ/ni/rWzl7tW0u7RvCrDL4c7Je7zSzdzH71fd3nRZ3BZmbvmtl2M1uczfNmZi/79ttCM7s42DV6IQf75XIz23PC5+XJYNfoBTOrZmY/mtkyM1tiZv38bOP5Z0YB4YeZhQOvkbWQUUPgVjNreMpmPYHdzrnawHCy1tQu0HK4XwC+cM418329HdQivfM+0Ok0z18N1PF99QJeD0JNoeB9Tr9fAGac8Hn5VxBqCgVHgb865xoALYHefv5b8vwzo4DwrwWw2jm31jl3GPicrFXwTtSFrLWzIat9+VW+9S0Kspzsl0LJOZcM7DrNJl2AD12WOWQtkFU5ONV5Jwf7pVByzm1xzs33/bwPWMafl1n2/DOjgPAvBth4wmN/a2T/ZxvfWtp7yFrHoiDLyX4BuNF3SPy1mVULTmkhL6f7rjBqZWYLzGyimTXyuphg852ebg7MPeUpzz8zCgj//B0JnHq7V062KWhy8r/5O6CGc64JMJX/HmUVdoXx85IT88lq+9AUeAVI9LieoDKzksA3QH/n3N5Tn/bzkqB+ZhQQ/m0CTvzL198a2f/ZxreWdhkK/qH0GfeLc26nc+6Q7+FbwCVBqi3U5eQzVeg45/Y65/b7fp4ARJpZBY/LCgoziyQrHD5xzo3xs4nnnxkFhH+/AHXMrKaZFQF6AONO2WYcWWtnQ9Za2tNcwZ9Ucsb9cso50s5knVuVrP10p+/OlJbAHufcFq+L8pqZnf/HtTsza0HWv0k7va0q8Hz/m98BljnnhmWzmeefGU+WHA11zrmjZtYHSALCgXedc0vM7F9AinNuHFn/537kW0N7F1n/WBZoOdwvfc2sM1l3aewC7vas4CAys8+Ay4EKZraJrCV1IwGcc28AE4BrgNXAAeAebyoNrhzsl+7AQ2Z2FMgEehSCP7QA4oA7gEVm9qtv7O9AdQidz4xmUouIiF86xSQiIn4pIERExC8FhIiI+KWAEBERvxQQIiLil25zFckDZnYMWHTCUFfn3G8elSOSJ3Sbq0geMLP9zrmSXtchkpd0ikkkQMws3MyGmtkvvuaFD3hdk0hu6BSTSN6IOmFG7Drn3A1krRmyxzl3qZkVBWaZ2WTn3DrvyhTJOQWESN7IdM41O2WsI9DEzLr7Hpcha/EXBYTkCwoIkcAx4BHnXJLXhYicDV2DEAmcJLIa0UUCmFldMyvhcU0iOaYjCJHAeRuoAcz3tXdOB7p6WpFILug2VxER8UunmERExC8FhIiI+KWAEBERvxQQIiLilwJCRET8UkCIiIhfCggREfFLASEiIn79PzbGMRTuG1FVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(corrosion.Fe, corrosion.loss)\n",
    "plt.xlabel(\"Fe\")\n",
    "plt.ylabel(\"loss\")\n",
    "xr = np.array(ax.get_xlim())\n",
    "ax.plot(xr, lmod.params[0] + lmod.params[1] * xr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the ANOVA-style model to get the replicate fits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "corrosion['Fefac'] = corrosion['Fe'].astype('category')\n",
    "amod = smf.ols(formula='loss ~ Fefac', data=corrosion).fit()\n",
    "ax.scatter(corrosion.Fe, amod.fittedvalues, marker='x')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last line contains the answer. Can ignore annoying warnings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda/lib/python3.7/site-packages/scipy/stats/_distn_infrastructure.py:879: RuntimeWarning: invalid value encountered in greater\n",
      "  return (self.a < x) & (x < self.b)\n",
      "/anaconda/lib/python3.7/site-packages/scipy/stats/_distn_infrastructure.py:879: RuntimeWarning: invalid value encountered in less\n",
      "  return (self.a < x) & (x < self.b)\n",
      "/anaconda/lib/python3.7/site-packages/scipy/stats/_distn_infrastructure.py:1821: RuntimeWarning: invalid value encountered in less_equal\n",
      "  cond2 = cond0 & (x <= self.a)\n"
     ]
    },
    {
     "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>df_resid</th>\n",
       "      <th>ssr</th>\n",
       "      <th>df_diff</th>\n",
       "      <th>ss_diff</th>\n",
       "      <th>F</th>\n",
       "      <th>Pr(&gt;F)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>11.0</td>\n",
       "      <td>102.850233</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>6.0</td>\n",
       "      <td>11.781667</td>\n",
       "      <td>5.0</td>\n",
       "      <td>91.068566</td>\n",
       "      <td>9.275621</td>\n",
       "      <td>0.008623</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   df_resid         ssr  df_diff    ss_diff         F    Pr(>F)\n",
       "0      11.0  102.850233      0.0        NaN       NaN       NaN\n",
       "1       6.0   11.781667      5.0  91.068566  9.275621  0.008623"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sm.stats.anova_lm(lmod, amod)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Polynomial fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  173.81821489,  -919.82377689,  1784.60949242, -1540.83916465,\n",
       "         552.23222816,   -76.09724981,   129.27393903])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pc = np.polyfit(corrosion.Fe, corrosion.loss, 6)\n",
    "pc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VNXdx/HPbyb7QvYEEiBsCasoEHGrG4uioqDi2iJWK7ZFq7aPVaSLto+FSmsfbd1QqWIRrRtSQZEKgiKCAVS2hCUQSEJIWLLvmfP8MTcS4wBBMnMnye/9es0rM3funflyM+Q3555z7xFjDEoppVRLDrsDKKWU8k9aIJRSSnmkBUIppZRHWiCUUkp5pAVCKaWUR1oglFJKeaQFQimllEdaIJRSSnmkBUIppZRHAXYHOBXx8fGmV69edsdQSql2Zf369QeNMQknWq9dF4hevXqRmZlpdwyllGpXRCS3NevpISallFIeaYFQSinlkRYIpZRSHmmBUEop5ZEWCKWUUh55rUCIyFwRKRKRzc2W/VFEvhaRL0XkQxFJtpaLiDwpIjut54d7K1dzCzfmc96s5fR+cDHnzVrOwo35vnhbpZRqF7zZgngJGNdi2WxjzFBjzBnAe8DvrOWXAWnWbSrwjBdzAe7iMP3tTeSXVGOA/JJqpr+9SYuEUkpZvFYgjDGrgMMtlpU1exgONM13OgGYZ9w+B6JFpJu3sgHMXppNdX3jt5ZV1zcye2m2N99WKaXaDZ+fKCcijwK3AKXAxdbiFGBfs9XyrGX7vZWjoKT6pJYrpVRn4/NOamPMDGNMD2A+cJe1WDyt6ml7EZkqIpkikllcXPy9cyRHhwKNLZY2WsuVUkrZOYrpVeBa634e0KPZc92BAk8bGWPmGGMyjDEZCQknvJTIMS3kPuKo4GgdamR98M9ZJj//3q+plFIdiU8LhIikNXt4FZBl3V8E3GKNZjobKDXGeO3wEo2NJDjKeSDgNdyNF8PaoLuIk3LCpA4aW7YslFKq8/HmMNcFwBqgv4jkicjtwCwR2SwiXwOXAPdYqy8BcoCdwPOAd7/GO50sHPUxo50b+H3Ay4Awq+Fmil1dWDjqY3A6vfr2SinVHnitk9oYc5OHxS8eY10DTPNWFk8eWZLFkdqn2RMymUpC+EvDDbxTdx7Ri7cxMaOnL6MopZRf6rRnUh+pqmF9sLuhMs35Lj91LgIclFTX4a5XSinVuXXOAtHo7pCOk3IOmUh6177CT5yLucX5IeDg7//dbndCpZSyXbueMOh7czqplWAOGRhR+zTgJKPuab4ImsZBE8XjH0Gty/DTC/sSGRLo83i1DY2U1zRQVl1P2Tc/62l0GWLCgty38EBiwoIIC3Ii4mmUsFJKnZrOWSCAdRNWcf8bGzjaiHJybsPTzLp2GCE7D/LUil3MW5PLj85O5cfn9SIxMqRN37/RZcg7UsXOogp2FVdYPyvZVVxBSVV9q18nKMBBUpdghveMIaNXLCN7xZKWGIHDoUVDKXVqpD0fb8/IyDCnMuXowo35zF6aTUFJNcnRodx/aX8mDksB4Ou8Ep5bmcOSzfsJdDqYNKI7U8/vQ6/48JN+n/KaerIKy9laUMa2/e5bVmE5tQ2ub9aJjwiib0IEfRMj6NYlhKiwQLqEBNIlNICoUPd9EThSVc+RyjqOVNW571fVse9wFV/sOUJxeS0A0WGBZKTGcFbvOK48PZmuUW1b3JRS7ZuIrDfGZJxwvc5cIFpj98FK5qzK4a0NeTQ0uji3bzzJ0SHERwSTEBlMfIT7FhEcQFF5DQWlNewvqWZ/aQ0FJdXkHakmv9nlO2LCAhnYrQsDu3UhPSmCfokR9E2IIDos6JRyGmPYe7iKdbsP88Wew2TuOULOwUqcDmH0gERuPqsnF6QlaMtCKaUFoq0Vldfwz9V7WLW9mIMVtRysqKPR5XnfBTiEpC4hJEeH0C0qlP5dIxnYLZKB3brQtUuIz/oMcg9VsmDdPt7I3Mehyjp6xIZy08ieXDeiBwmRwT7JoJTyP1ogvMzlMpRU17uLRXkt5bUNJEYGkxwdSnxEME4/+qZe29DIh1sOMH9tLp/nHCbQKdxwZg/uG5NOXIQWCqU6Gy0QyqNdxRXM/XQ3r32xj7BAJ3eN6set5/UiOEDPHleqs9ACoY5rZ1E5f1qSxfKsInrEhvLguIFcflrXbx3+Ol4nvlKq/WptgeicJ8op+iVGMvfWM3nl9pGEBwUw7dUNXPfsGjbnlwI6455SSlsQCvc5GW9k7uMvH26npKqOX16Szr/W5FJQWsW3v0O4SIkOZ/WDo+yKqpRqA9qCUK3mdAg3juzJf395AZcMTuKxD7LpUpbN3MDZQNO5Gi4WB83gr1XT7YyqlPIhLRDqG9FhQTx183D+Omkoe00S99Tfza+cb9BUHAY7col11oLLdcLXUkq1f1og1LeICNdm9OBXl59BTznAXxtvYLxjLd3lIFtdqWwd/x9w6MdGqc6g016LSR1fXGQI2dKL+52v8beGSWys7ccBRzx/0eKgVKeh/9uVR7Pf38a7AQ8xLWARbwY9QiMOnK56/rBos93RlFI+ogVCfZfLxZyaXzLYkcsWVyoT6x5hduBzDJY9HK5u4PGlWbiOcZkRpVTHoQVCfZfDQa0zgi2uVK6oexRwMLn+QX4XMI/Rzi95csUu7vzXeipqG+xOqpTyIq8VCBGZKyJFIrK52bLZIpIlIl+LyDsiEt3suekislNEskXkUm/lUq2z96o3mGRmcfQj4uAmHmX8tVP4/ZWDWJ5VxDVPryb3UKWdMZVSXuTNFsRLwLgWy5YBQ4wxQ4HtwHQAERkE3AgMtrZ5WkT04kA2mjgshZnXnE5KdCgCpESHMvOa07l6eHd+fF5v5t02kgNltUx4ajXrc4/YHVcp5QVePZNaRHoB7xljhnh47mpgkjHmhyIyHcAYM9N6binwsDFmzfFeX8+ktlfuoUqmzF1HYVkNT908nNEDk+yOpJRqhfZwJvVtwPvW/RRgX7Pn8qxlyo+lxoXz5s/OpX9SJFNfWc/rX+y1O5JSqg3ZUiBEZAbQAMxvWuRhNY9NGxGZKiKZIpJZXFzsrYiqleIjgnn1jrP5Qb94HnhrE09+tIP2fH0vpdRRPi8QIjIFGA/80Bz9S5IH9Gi2WnegwNP2xpg5xpgMY0xGQkKCd8OqVgkPDuCFKRlcMzyFx5dt5zcLNx9ztj2lVPvh0zOpRWQc8ABwoTGmqtlTi4BXReRxIBlIA9b5Mps6NYFOB3+97nS6dgnh6Y93UVxey5M3DSMkUMcaKNVeeXOY6wJgDdBfRPJE5HbgH0AksExEvhSRZwGMMVuAfwNbgQ+AacaYRm9lU94hIvx63AAevnIQy7Yd4LaXvqBSz5VQqt3S+SCUV7yzMY//eeNrhnaP4qVbRxIVFmh3JKWUpT2MYlId2NXDuvP0D4ezJb+MG+asobi81u5ISqmTpAVCec2lg7vy4q0Z5B6q4obn1pBfUm13JKXUSdACobzq/LQEXrl9JMXltVz/7Bp2H9RLcyjVXmiBUF6X0SuWBVPPprq+keueXUN2YbndkZRSraAFQvnEkJQo/n3nOTgdcMOcNWzKK7U7klLqBLRAKJ/plxjBG3eeS0RwADc//znrcw/bHUkpdRxaIJRP9YwL4993nkNCZDCTX1zHZzsP2h1JKXUMWiCUzyVHh/LanWfTIyaMH7/0BSuyi+yOpJTyQAuEskViZAivTT2btKQIps7L5IPN++2OpJRqQQuEsk1MeBCv3nE2Q7tHM+3VjbyzMc/uSEqpZrRAKFt1CQlk3m0jOat3LPe9/hWvrNljdySllEULhLJdeHAAc289kzEDk/jtu1t4asVOuyMppfDx5b6VOpaQQCfP/Gg497/xFbOXZlNe08AD4/ojIizcmM/spdkUlFSTHB3K/Zf2Z+IwnXBQKW/TAqH8RqDTwePXn0F4cADPrtxFeU09I1JjmPHOZqrr3Vd/zy+pZvrbmwC0SCjlZVoglF9xOIT/nTiELqGBPPPxLt7ekE91fT3Nj4ZW19cze2m2FgilvEwLhPI7IsID4wYQGRLAYx9kc6Zs5wuTjrtIuFgcNIPyqlBgrc1JlerYtJNa+a2fX9CHXwQuZL1JJ13yaSoOgx25xDprweWyO6JSHZoWCOW/HA76TPwdMwL+Ra5JorccIJJqtppUto7/Dzj046uUN+n/MOXXJo7owe7hD/Fq0KMcMRFcU/cIf+/3IhNH9LA7mlIdnhYI5dcWrt/HzV9NZoRjB28GPUIwdazYso+Zi7faHU2pDs9rBUJE5opIkYhsbrbsOhHZIiIuEclosf50EdkpItkicqm3cql2xOVi0HtXMkhy2eJKZUzdY8wMfJ7esp8XPtnJG1/stTuhUh2aN1sQLwHjWizbDFwDrGq+UEQGATcCg61tnhYRpxezqfbA4eBwYzBbXKlcUfco4OCW+gf5Q8BLDJRc7n9rEzOXbKPRZexOqlSH5LVhrsaYVSLSq8WybeAextjCBOA1Y0wtsFtEdgIjgTXeyqfah1+FzSS/pJKj32UcXFf/O5Kjwpg8MInnVuWQfaCcJ24cRlRooJ1Rlepw/KUPIgXY1+xxnrVMdXL3X9qf0MBv/+EPDQzk1+MG8MeJQ/jT1afx6Y6DXP30anKKK2xKqVTH5C8F4jtNCsDjcQMRmSoimSKSWVxc7OVYym4Th6Uw85rTSIkORYCU6FBmXnPaN2dR33xWT+b/5CxKquqZ8NRqPtbJh5RqM/5yJnUe0HzcYnegwNOKxpg5wByAjIwMPfjcCUwclnLcy2qc1SeORXedxx3z1nPbS1/w4GUDuOP8Pp4OZSqlToK/tCAWATeKSLCI9AbSgHU2Z1LtSPeYMN762TmMG9KVPy3J4o55mRyurLM7llLtmjeHuS7A3cncX0TyROR2EblaRPKAc4DFIrIUwBizBfg3sBX4AJhmjGn0VjbVMYUFBfDUzcN5+MpBrNp+kMueWMVnuw7aHUupdkuMab9HaTIyMkxmZqbdMZQf2lJQyt0LNrL7YCU/v6gv945JJ9DpLw1mpewlIuuNMRknWk//x6gOaXByFO/d/QOuG9Gdp1bs4vrn1rDvcJXdsZRqV7RAqA4rLCiAxyadzpM3DWPngQouf+ITXl27F5eeWKdUq2iBUB3eVacns+Se8xmc0oWH3tnEdc+tIauwzO5YSvk9LRCqU+gRG8aCO87mL9edTk5xBeOf/JSZ72+jqq7B7mhK+S0tEKrTEBEmjejO8l9dxDXDU3huZQ6X/G0VK7L05DqlPNECoTqdmPAgHpt0Oq9PPZuQQCc/fukLbpm7jq/2ldgdTSm/ogVCdVpn9YljyS/OZ/plA9iUV8KEp1bzk5cz2Vqg/RNKgZ4HoRQA5TX1vLR6D3M+yaG8poErTuvGfWPT6JcYaXc0pdpca8+D0AKhVDOlVfW88GkOcz/dTXV9I5cM6soPz+7JeX3jcTiaXdvJGJDjPFbKj2mBUOoUHK6s4/lPcnht3V6OVNWTGhfGzSN7MmlEd+K+eJxd+/K5Jf9qCkprSI4KYV7KO/TtkQIXT7c7ulInpGdSK3UKYsODeGDcAD5/aDRP3HgGSV1CmPl+FufMXM5ta+Ip2LGRWyuex2C4vXIOfXNeYde+fHdLQqkOwl8u962UXwoOcDLhjBQmnJHCjgPlzF+7l5c/c7Gch4ikkisda4iXUp6uv5L5+VezWg8zqQ5EC4RSrZSWFMnDVw3mpc/2AIZxzi9Y3jiM/7jOBQyU1vDip7s5s1cMA7t10YsDqnZPC4RSJyklKoTbK+dwW8AHNAYIG00aj9dPYp2cxh/f2wpASKCDod2jGd4zhhGpMQzrGU18RLDNyVVH8eb6PAZ0jWRISpRX30cLhFInwxh3h3TOB8xtGMcfGibzu4BXeDX4T+zqM5nQKx9jw74SNuSWsH7vEV74JIdnV7r7JeIjgkhLjCQtKYK0pEjSEyNIT4okJjzI5n+Uak8aXYYH3/qaOy/sowVCKb8iQt8eKexiMi/mX42U1vBi+FQuTElwj2KKCSM5JozxQ5MBqKlvZFN+KV/tK2H7gXK2H6jg7Q35VNQevQZUVGggqXFh9IwNIzUujNS4cFJjw+gVH05iZHDrpk7VYbedRlF5DQ0uQ3J0qNffSwuEUifr4un0NebbHdJmlMc/yCGBTs7sFcuZvWKPrmoM+0tr2H6gnJ1FFew5VEnuoSq+zivl/c2FNDa7HHlYkJPUuHB6xbkLRq+4MPolRjKgayThwdZ/3xUzddhtJ1JQUg2gBUIpv9WyGJzEt3URITk6lOToUC7qn/it5+obXRSUVLPnUBW5hyrZc7CKPYcqyS4sZ9nWAzRYxUMEUmPDGNStC0kHQjivZAvXuyr4G5OsYbcfsIvJ9NWWRIeTX1IDQHctEEp1LoFOh/sQU1w4kPCt5xoaXRSU1JB9oJxt+8vYtr+MLfvLWHIonX/yPwCkyz72mQQerLudlXmXs0aLQ4fT1ILo1p4LhIjMBcYDRcaYIdayWOB1oBewB7jeGHNE3AdZnwAuB6qAW40xG7yVTan2KMDpoGdcGD3jwhg7KOmb5b0eXAwYpgcs4FPXEF5tHE0tQVBWy/XPruGC9HgmDkuhe0yYfeFVmykoqSYqNJCIYO9/v/fmQO2XgHEtlj0IfGSMSQM+sh4DXAakWbepwDNezKVUh5ISFcLvAl7hzoD3eCVoFl8F38Fkx4dEBDupaWjkr8u2c/5jK5j84loWf72f2oZGuyOrU5B/pNon/Q/gxRaEMWaViPRqsXgCcJF1/2XgY+ABa/k8474w1OciEi0i3Ywx+72VT6kO4RjDbv8Y9BK39mmk7+S/k1dSzRuZebyRuY9pr24gJiyQa4Z354Yze5CepFerbW/yS6rpHtPOC8QxJDX90TfG7BeRph66FGBfs/XyrGVaIJQ6nhMNuxWhe0wY941N5xej0/hkRzH/ztzHvDV7ePHT3Vx+Wld+fekAesWH2/0vUa1UUFLNyN6xJ16xDfhLJ7WnnjSPVz0Tkam4D0PRs2dPb2ZSqn1o5bBbp0O4qH8iF/VP5GBFLfPW5PL8qhyWbT3AD89K5Rej04jVk/b8WnlNPWU1DaT46BCTry8Wc0BEugFYP5smA84DejRbrztQ4OkFjDFzjDEZxpiMhIQET6so1fmc5LDb+Ihgfjk2nZX3X8SkET2Yt2YPFz62gqdW7KSmXvso/NX+UvcQV1/1Qfi6QCwCplj3pwDvNlt+i7idDZRq/4NS3pfYJYSZ15zG0nsv4Kw+scxems3Ff/mY5VkH7I6mPMj34Uly4MUCISILgDVAfxHJE5HbgVnAWBHZAYy1HgMsAXKAncDzwM+9lUsp9V1pSZG8MOVMXpt6NlGhgdz2UiaPLt5KXYPL7miqmfwj7gLhq0NM3hzFdNMxnhrtYV0DTPNWFqVU65zdJ46F087jT0u28fwnu1m35wj/uGkYPWL1HAp/UFBSTYBDSIj0zZWB9YL1SqlvCQl08ocJQ3jmh8PJKa7g8ic/4f1NesTXHxSUVNM1KgSnwzdnyLeqQIjIPSLSxeojeFFENojIJd4Op5Syz2WndWPJL86nT0IEP5u/gd8u3Kwd2DYrKKnx2eElaH0L4jZjTBlwCe4LxPyYo/0HSqkOqkdsGG/ceQ5TL+jDK5/n8qMX1lJWU293rE4rv6TaLwtEU3vmcuCfxpiv8HzuglKqgwkKcPDQ5QP5x83D+HJfCT96YS0lVXV2x+p0GhpdFJbV+GwEE7S+QKwXkQ9xF4ilIhIJ6PAGpTqR8UOTeW7yCLL2l3PjnM85WFFrd6ROpai8lkYfTRTUpLUF4nbcF9Y70xhTBQTiPsyklOpERg9M4sVbM9hzqJIbnltDoXXilvK+oxMFhfjsPVtbIM4Bso0xJSLyI+A3QKn3Yiml/NX5aQnMu+0sCktruP65NeQdqbI7UqfQdJKcP/ZBPANUicjpwK+BXGCe11IppfzayN6x/OsnZ1FSVcf1z65hz8FKuyN1eAUlvr3MBrS+QDRYJ7NNAJ4wxjwB6HWClerEhvWMYcHUs6lpcHHz859TVK6Hm7ypoKSa6LDAo3OR+0BrC0S5iEwHJgOLRcSJux9CKdWJDU6OYt5tIzlcVcedr6zX8yS8KL+kmuQo37UeoPUF4gagFvf5EIW452qY7bVUSql2Y0hKFH+7/gw27i1h+tubcB9sUG2toMR3M8k1aVWBsIrCfCBKRMYDNcYY7YNQSgHus67/55J03tmYz9Mf77I7TofkPknOdyOYoPWX2rgeWAdcB1wPrBWRSd4MppRqX6Zd3I8JZyQze2k2H2zWaze1pbKaesprGnzegmhtb8cM3OdAFAGISALwX+BNbwVTSrUvIsKfrx1K7qEq7nv9K7rHhDEkJcruWB3CfmsEU4qP5qJu0to+CEdTcbAcOoltlVKdREigkzm3jCAmLJA75mUyb3UO581aTu8HF3PerOUsXL/vxC+ivqPAxxMFNWntH/kPRGSpiNwqIrcCi3FP8qOUUt+SGBnC81MyOFJWzvz3lpJfUokB8ksqSV90BQf//p0pYdQJ5Nlwkhy0vpP6fmAOMBQ4HZhjjHnAm8GUUu3X4K6R/E/gm2Sbntzo+BhwsThoBoMkl8OHD4FLL+V2MgpKqgl0CgkRvpkoqEmrz7gwxrwFvOXFLEqpjsLh4NG6GxjjyOTfrot4K+gRBjty2eJKZXzdH9nt0CPUJ6NpoiCHjyYKanLc35KIlItImYdbuYiU+SqkUqr9iQoN5r+u4XTjEL+s/xmVJpgr6h4lKtS334I7ggIfzwPR5LgFwhgTaYzp4uEWaYzp4quQSqn2R4z7sNLjQc+w1yTyvw0/YnHQDMTo4aWTVVDi23kgmtjSzrOmMN0sIltE5F5rWayILBORHdbPGDuyKaXagMvFv1y/ZrAjlwiqMQgLGkdTaGL5l+vX2gdxEpomCvK7FoQ3iMgQ4A5gJO4O7/EikoZ7vomPjDFpwEfWY6VUe+RwUOuMYIsrlSvqHsU9AaXhV/U/5YAjEbQPotUO2DBRUBM7fksDgc+NMVXGmAZgJXA17ivFvmyt8zIw0YZsSqk2sveqN5hkZnH0z4xQQiSPx/1er9d0Euw6BwLsKRCbgQtEJE5EwnBPY9oDSDLG7AewfibakE0p1UYmDkth5jWnkxIdiuAewz/hjGQ255fxRmae3fHajYJvzoHw7XWY4CSGubYVY8w2EfkzsAyoAL4CGlq7vYhMBaYC9OzZ0ysZlVJtY+KwFCYOS/nmsctlKCqr5ZH/bOG8tHhbjqu3N/mdrAWBMeZFY8xwY8wFwGFgB3BARLoBWD+LjrHtHGNMhjEmIyEhwXehlVKnzOEQZl83lEZj+N/3ttodp10oKKkmJiyQsCCff5+3bRRTovWzJ3ANsABYBEyxVpkCvGtHNqWUd3WPCePuUWm8v7mQlduL7Y7j9/KP+H4eiCZ2DSV4S0S2Av8BphljjgCzgLEisgMYaz1WSnVAPzm/N73jw3l40RZqG3QWuuOx6xwIsO8Q0/nGmEHGmNONMR9Zyw4ZY0YbY9Ksn4ftyKaU8r7gACcPXzWY3QcreeGT3XbH8Wt2nUUNeslupZRNLkxP4LIhXfn78h3kHamyO45fKqupp7y2gWQbRjCBFgillI1+M34QgvBH7bD26OgQ1zBb3l8LhFLKNinRodw9uh9LtxxgRbbHgYud1sKN+dz8/FoAHl60hYUb832eQQuEUspWP/lBH/okuDusa+q1wxrcxWH625s4XFkHQHFFLdPf3uTzIqEFQillq6AAB49cNZjcQ1U8vyrH7jh+YfbSbKpbFMvq+kZmL832aQ4tEEop252flsDlp3XlHyt2su+wdli7+x5aXvHW9U2fhK9ogVBK+YXfXDEIgL9+6Ntvyf7ordBHWRw0A2i6qKF7bo23Qh/1aQ4tEEopv5AcHcptP+jNwi8L2FJQancc+7hc9IlsZLAjlwRKAcPioBkMduTSJ7LRp3NpaIFQSvmNn17Yl6jQQB77oBO3IhwOou9dQ3HkIA7RhV8432GwI5eSLgOIvneNT+fS0AKhlPIbUaGB3HVxP1ZuL+aznQftjmMfh4OaW5fhwkE3OQTg8+IAWiCUUn5m8jmpJEeFMOuDrM47sZDLxVdz7gSgj2M/ACX/d47Pp2rVAqGU8ishgU7uG5vO13mlLNlUaHcc33O5KPm/c8irdD+8oe43bHGlEl2W5fMioQVCKeV3rhnenfSkCGYvzaK+0bffmm3ncJBT7mStayDuUUwOrqh7lC2uVHLKndoHoZTq3JwO4YFxA9hzqIrXvthndxyfu7Z6BitcZwBiLXEXiWurZ/g0hxYIpZRfGjUgkZG9YnnivzuorG31rMQdQreoEI4WhyYOn88LoQVCKeWXRIQHLhvAwYpa5n7aueaMmHxO6neWhQY6uf/S/j7NoQVCKeW3RqTGcOngJJ5blcOhilq74/hMSoz78t6JkcEI7qvezrzmNCYOS/FpDt/Pgq2UUifh/ksHsGzrSp7+eBe/HT/I7jg+kbW/jACH8OkDowgKsO97vLYglFJ+rV9iBFcP686/Ps+lqLzG7jg+kV1YTp+EcFuLA9hUIETkPhHZIiKbRWSBiISISG8RWSsiO0TkdREJsiObUsr/3D2qHw0uw7Mfd47LgWcVltO/axe7Y/i+QIhICvALIMMYMwRwAjcCfwb+ZoxJA44At/s6m1LKP/WKD+fqYSnMX5tLUVnHbkWU1dSTX1LNgK6Rdkex7RBTABAqIgFAGLAfGAW8aT3/MjDRpmxKKT/0TStiZcduRWwvLAfonAXCGJMP/AXYi7swlALrgRJjTNNg5zzAt931Sim/lhoXzjWdoBWRZRWI/p2xQIhIDDAB6A0kA+HAZR5W9XiVLhGZKiKZIpJZXFzsvaBKKb9zl9WKeGblLrujeE12YTmRwQGk+PikOE/sOMQ0BthtjCk2xtQDbwPnAtHWISeA7kCBp42NMXOMMRnGmIyEhATfJFZK+YWjrYi9HOigrYjswnL6d41EpOWZ1L5nR4HYC5wtImHi3gOjga3ACmCStc4U4F0bsiml/Nzdo9JodBme+bijplj3AAAP5klEQVTjtSKMMWwrLPOLw0tgTx/EWtyd0RuATVaGOcADwC9FZCcQB7zo62xKKf/XMy6Ma4en8Oq6jteK2F9aQ3lNg190UINNo5iMMb83xgwwxgwxxkw2xtQaY3KMMSONMf2MMdcZYzrPefVKqZNy18VpuDpgKyL7mw5q+8+BAD2TWinVDrlbEd15dd1eCks7TivCn0YwgRYIpVQ7ddeoflYrYqfdUdpMVmEZyVEhRIUG2h0F0AKhlGqnesS6WxELvtjXYc6LaBrB5C+0QCil2q2fX9yXRpfhuVXt/+zq+kYXu4or/Kb/AbRAKKXasdS4cCacnsz8tbkcbOfzReQUV1LfaBjYTVsQSinVJqaN6kdtg4sXPmnfs85lFZYB/tNBDVoglFLtXN+ECMYPTeaVNXs4Ullnd5zvLauwnACH0Cc+wu4o39ACoZRq9+66uB+VdY3MXd1+WxHZheX0TYiwfZKg5vwniVJKfU/9u0YybnBXXlq9h9LqervjfC/ZheUM8KP+B9ACoZTqIO4e3Y/y2gZe/myP3VFOWmm1e5Igf+p/AC0QSqkOYnByFGMGJjJ39W4qahtOvIEf2X7AfyYJak4LhFKqw7h7VBolVfW8sibX7ignJcvPrsHURAuEUqrDOL1HNBekJ/D8JzlU1bWfVkR2YRmRIQEkR4XYHeVbtEAopTqUe0b343BlHa+u3Wt3lFbLLixngJ9MEtScFgilVIcyIjWWc/vG8ezKHKrrGu2Oc0LGGLL87BpMTbRAKKU6nHtGp3Gwopb5a/2/L6LAmiTI3/ofQAuEUqoDOqtPHOf1i+PZlbv8vi9iW4H7Ehv+NoIJtEAopTqo+8akc7Cizu9HNH2yo5iQQAenpUTZHeU7tEAopTqkjF6xnJ8Wz3Orcqj00/MijDEszy7iB/3iCQl02h3nO7RAKKU6rPvGpnO4so6X1+yxO4pHu4or2He4mosHJNodxSOfFwgR6S8iXza7lYnIvSISKyLLRGSH9TPG19mUUh3L8J4xXNQ/gTmrciiv8b9rNH20rQiAi/trgQDAGJNtjDnDGHMGMAKoAt4BHgQ+MsakAR9Zj5VS6pTcNyadkqp6v7xG0/KsIgZ0jSQ5OtTuKB7ZfYhpNLDLGJMLTABetpa/DEy0LZVSqsM4vUc0owckMmdVDmV+1Ioora4nM/cIowf6Z+sB7C8QNwILrPtJxpj9ANZP/91rSql25b6x6ZTVNPDPT/fYHeUbq7YX0+gyjPLT/gewsUCISBBwFfDGSW43VUQyRSSzuLjYO+GUUh3KkJQoxg5K4oVPc/xmvogVWUXEhAVyRg//7W61swVxGbDBGHPAenxARLoBWD+LPG1kjJljjMkwxmQkJCT4KKpSqr27d0wa5TUNvPip/bPONboMH28v5qL+iTgd/nX9pebsLBA3cfTwEsAiYIp1fwrwrs8TKaU6rMHJUYwb3JW5n+7mUEWtrVm+yivhcGWd3w5vbWJLgRCRMGAs8HazxbOAsSKyw3pulh3ZlFId168uSae6vpG//Xe7rTmWbyvC6RAuTPPvoyC2FAhjTJUxJs4YU9ps2SFjzGhjTJr187Ad2ZRSHVdaUiQ/Oqsnr67dS7Y1SY8dlmcVMaJnDFFhgbZlaA27RzEppZRP3TsmnYjgAP743laMMT5//8LSGrbuL2OUHw9vbaIFQinVqcSEB3HvmHQ+3XmQ5Vkex8J4VdN7+vPw1iZaIJRSnc7kc1LpkxDOo4u3Udfg8ul7L88qIiU6lLTECJ++7/ehBUIp1ekEOh385oqB5Bys5JXPfXc58Jr6RlbvPMjogYl+N72oJ1oglFKd0sX9Ezk/LZ4n/rudw5V1PnnPtbsPU13f6PfDW5togVBKdUoiwm/HD6KyrpH/89Gw1+XbDhAS6OCcPnE+eb9TpQVCKdVppSdF8sOzejJ/7V62H/DusFd/nxzIEy0QSqlO7d4x6YQHOb0+7NXfJwfyRAuEUqpTiw0P4p4x6Xyy4yCLvirw2vss/roQ8N/JgTzRAqGU6vRuOSeVjNQYHnp7E7uKK9r89feXVvPcql2MGZjot5MDeaIFQinV6QU6Hfz95mEEBTiYNn8DNfWNbfr6f3xvK40uw++vHNymr+ttWiCUUgroFhXK4zecQVZhOY/8Z0ubve7K7cUs2VTI3aP60SM2rM1e1xe0QCillOXi/on87KK+LFi3j4Ub80/59WrqG/n9u5vpEx/OHRf0aYOEvqUFQimlmvnV2HRG9orloXc2sbPo1PojnluZw55DVfxhwhCCA9rH0NbmtEAopVQzAU4HT940jJBAJ9Pmb6C67vv1R+QequSpj3cyfmg3fpAW38YpfUMLhFJKtdA1KoS/3XAG24vKeXjRyfdHGGN4eNEWgpwOfjt+kBcS+oYWCKWU8uDC9ASmXdSP1zP3Mev9LGobWt+SWLrlACuyi7lvbDpJXUK8mNK7tEAopdQx3DsmjRsyevDsyl2Mf/JTNu49csJtKmsb+MN/tjCgayRTzkn1QUrv0QKhlFLHEOB08OdJQ3npx2dSUdvAtc98xp+WbPN4nkRVXQMfbN7Pz+ZvoKC0hkevHkKAs33/iQ2w401FJBp4ARgCGOA2IBt4HegF7AGuN8acuFwrpZSXXdQ/kQ/vu4A/Lclizqoclm09wGOThtInPpyPthXx4dZCPtlxkNoGF1Ghgfx6XH9GpMbaHfuUiR1zsorIy8AnxpgXRCQICAMeAg4bY2aJyINAjDHmgeO9TkZGhsnMzPRBYqWUclu98yAPvPU1+SXVCOAykBIdythBSVwyOImRvWL9vuUgIuuNMRknXM/XBUJEugBfAX1MszcXkWzgImPMfhHpBnxsjOl/vNfSAqGUskNlbQPPf5KDy8Alg5IYnNylXcwQ16S1BcKOQ0x9gGLgnyJyOrAeuAdIMsbsB7CKRPu55KFSqlMJDw7g3jHpdsfwOjvaQQHAcOAZY8wwoBJ4sLUbi8hUEckUkczi4mJvZVRKqU7PjgKRB+QZY9Zaj9/EXTAOWIeWsH4WedrYGDPHGJNhjMlISEjwSWCllOqMfF4gjDGFwD4RaepfGA1sBRYBU6xlU4B3fZ1NKaXUUbYMcwXuBuZbI5hygB/jLlb/FpHbgb3AdTZlU0ophU0FwhjzJeCpB320r7MopZTyzL8H6yqllLKNFgillFIeaYFQSinlkS2X2mgrIlIM5J7iy8QDB9sgTlvTXK3nj5lAc50szXVyTiVXqjHmhOcJtOsC0RZEJLM1p5z7muZqPX/MBJrrZGmuk+OLXHqISSmllEdaIJRSSnmkBQLm2B3gGDRX6/ljJtBcJ0tznRyv5+r0fRBKKaU80xaEUkopjzpsgRCRcSKSLSI7rRnqWj4fLCKvW8+vFZFezZ6bbi3PFpFLfZzrlyKyVUS+FpGPRCS12XONIvKldVvk41y3ikhxs/f/SbPnpojIDus2peW2Xs71t2aZtotISbPnvLK/RGSuiBSJyOZjPC8i8qSV+WsRGd7sOW/uqxPl+qGV52sR+cyaj6XpuT0issnaV206C1crcl0kIqXNfle/a/bccX//Xs51f7NMm63PU6z1nFf2l4j0EJEVIrJNRLaIyD0e1vHd58sY0+FugBPYhXtyoiDcM9gNarHOz4Fnrfs3Aq9b9wdZ6wcDva3Xcfow18VAmHX/Z025rMcVNu6vW4F/eNg2FvcFF2OBGOt+jK9ytVj/bmCuD/bXBbgvUb/5GM9fDrwPCHA2sNbb+6qVuc5tej/gsqZc1uM9QLxN++si4L1T/f23da4W614JLPf2/gK6AcOt+5HAdg//F332+eqoLYiRwE5jTI4xpg54DZjQYp0JwMvW/TeB0SIi1vLXjDG1xpjdwE7r9XySyxizwhhTZT38HOjeRu99SrmO41JgmTHmsDHmCLAMGGdTrpuABW303sdkjFkFHD7OKhOAecbtcyBa3HOceHNfnTCXMeYz633Bd5+t1uyvYzmVz2Vb5/LVZ2u/MWaDdb8c2AaktFjNZ5+vjlogUoB9zR7n8d2d/M06xpgGoBSIa+W23szV3O24vyk0CRH3bHqfi8jENsp0MrmutZq0b4pIj5Pc1pu5sA7F9QaWN1vsrf11IsfK7c19dbJafrYM8KGIrBeRqTbkOUdEvhKR90VksLXML/aXiITh/kP7VrPFXt9f4j7sPQxY2+Ipn32+7JoPwts8zR7ecrjWsdZpzbbfV6tfW0R+hPuS6Bc2W9zTGFMgIn2A5SKyyRizy0e5/gMsMMbUishPcbe+RrVyW2/manIj8KYxprHZMm/trxOx47PVaiJyMe4C8YNmi8+z9lUisExEsqxv2L6wAfelHypE5HJgIZCGn+wv3IeXVhtjmrc2vLq/RCQCd0G61xhT1vJpD5t45fPVUVsQeUCPZo+7AwXHWkdEAoAo3M3N1mzrzVyIyBhgBnCVMaa2abkxpsD6mQN8jPvbhU9yGWMONcvyPDCitdt6M1czN9LiEIAX99eJHCu3N/dVq4jIUOAFYIIx5lDT8mb7qgh4h7Y7rHpCxpgyY0yFdX8JECgi8fjB/rIc77PV5vtLRAJxF4f5xpi3Paziu89XW3ey+MMNd8soB/chh6bOrcEt1pnGtzup/23dH8y3O6lzaLtO6tbkGoa7Yy6txfIYINi6Hw/soI067FqZq1uz+1cDn5ujHWO7rXwx1v1YX+Wy1uuPu9NQfLG/rNfsxbE7Xa/g252I67y9r1qZqyfuPrVzWywPByKb3f8MGOfDXF2bfne4/9DutfZdq37/3splPd/0xTHcF/vL+nfPA/7vOOv47PPVZjva3264e/q34/5jO8Na9gfc38oBQoA3rP8w64A+zbadYW2XDVzm41z/BQ4AX1q3Rdbyc4FN1n+STcDtPs41E9hivf8KYECzbW+z9uNO4Me+zGU9fhiY1WI7r+0v3N8m9wP1uL+13Q78FPip9bwAT1mZNwEZPtpXJ8r1AnCk2Wcr01rex9pPX1m/4xk+znVXs8/W5zQrYJ5+/77KZa1zK+5BK82389r+wn3YzwBfN/s9XW7X50vPpFZKKeVRR+2DUEopdYq0QCillPJIC4RSSimPtEAopZTySAuEUkopjzrqmdRK+ZSINOIecthkojFmj01xlGoTOsxVqTYgIhXGmAi7cyjVlvQQk1JeIiJOEZktIl9YFzm80+5MSp0MPcSkVNsIFZEvrfu7jTFX4z4zt9QYc6aIBAOrReRD476MvFJ+TwuEUm2j2hhzRotllwBDRWSS9TgK91VKtUCodkELhFLeI8DdxpildgdR6vvQPgilvGcp8DPr8s2ISLqIhNucSalW0xaEUt7zAu7LSW+wprMtBnw5s51Sp0SHuSqllPJIDzEppZTySAuEUkopj7RAKKWU8kgLhFJKKY+0QCillPJIC4RSSimPtEAopZTySAuEUkopj/4ffTr/ETT72JwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(corrosion.Fe, corrosion.loss)\n",
    "plt.xlabel(\"Fe\")\n",
    "plt.ylabel(\"loss\")\n",
    "ax.scatter(corrosion.Fe, amod.fittedvalues, marker='x')\n",
    "grid = np.linspace(0,2,50)\n",
    "ax.plot(grid,np.poly1d(pc)(grid))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Manually calculate R-squared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.99653135])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pc, rss, _, _, _ = np.polyfit(corrosion.Fe, corrosion.loss, 6, full=True)\n",
    "tss = np.sum((corrosion.loss-np.mean(corrosion.loss))**2)\n",
    "1-rss/tss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Robust Regression\n",
    "\n",
    "Load Galapagos data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "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>Species</th>\n",
       "      <th>Area</th>\n",
       "      <th>Elevation</th>\n",
       "      <th>Nearest</th>\n",
       "      <th>Scruz</th>\n",
       "      <th>Adjacent</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Baltra</th>\n",
       "      <td>58</td>\n",
       "      <td>25.09</td>\n",
       "      <td>346</td>\n",
       "      <td>0.6</td>\n",
       "      <td>0.6</td>\n",
       "      <td>1.84</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bartolome</th>\n",
       "      <td>31</td>\n",
       "      <td>1.24</td>\n",
       "      <td>109</td>\n",
       "      <td>0.6</td>\n",
       "      <td>26.3</td>\n",
       "      <td>572.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Caldwell</th>\n",
       "      <td>3</td>\n",
       "      <td>0.21</td>\n",
       "      <td>114</td>\n",
       "      <td>2.8</td>\n",
       "      <td>58.7</td>\n",
       "      <td>0.78</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Champion</th>\n",
       "      <td>25</td>\n",
       "      <td>0.10</td>\n",
       "      <td>46</td>\n",
       "      <td>1.9</td>\n",
       "      <td>47.4</td>\n",
       "      <td>0.18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Coamano</th>\n",
       "      <td>2</td>\n",
       "      <td>0.05</td>\n",
       "      <td>77</td>\n",
       "      <td>1.9</td>\n",
       "      <td>1.9</td>\n",
       "      <td>903.82</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           Species   Area  Elevation  Nearest  Scruz  Adjacent\n",
       "Baltra          58  25.09        346      0.6    0.6      1.84\n",
       "Bartolome       31   1.24        109      0.6   26.3    572.33\n",
       "Caldwell         3   0.21        114      2.8   58.7      0.78\n",
       "Champion        25   0.10         46      1.9   47.4      0.18\n",
       "Coamano          2   0.05         77      1.9    1.9    903.82"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gala = pd.read_csv(\"data/gala.csv\",index_col=0)\n",
    "gala.drop('Endemics', axis=1, inplace=True)\n",
    "gala.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Least squares fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>Species</td>     <th>  R-squared:         </th> <td>   0.766</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.717</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   15.70</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 25 Sep 2018</td> <th>  Prob (F-statistic):</th> <td>6.84e-07</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:55:21</td>     <th>  Log-Likelihood:    </th> <td> -162.54</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    30</td>      <th>  AIC:               </th> <td>   337.1</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    24</td>      <th>  BIC:               </th> <td>   345.5</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     5</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>    7.0682</td> <td>   19.154</td> <td>    0.369</td> <td> 0.715</td> <td>  -32.464</td> <td>   46.601</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Area</th>      <td>   -0.0239</td> <td>    0.022</td> <td>   -1.068</td> <td> 0.296</td> <td>   -0.070</td> <td>    0.022</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Elevation</th> <td>    0.3195</td> <td>    0.054</td> <td>    5.953</td> <td> 0.000</td> <td>    0.209</td> <td>    0.430</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Nearest</th>   <td>    0.0091</td> <td>    1.054</td> <td>    0.009</td> <td> 0.993</td> <td>   -2.166</td> <td>    2.185</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Scruz</th>     <td>   -0.2405</td> <td>    0.215</td> <td>   -1.117</td> <td> 0.275</td> <td>   -0.685</td> <td>    0.204</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Adjacent</th>  <td>   -0.0748</td> <td>    0.018</td> <td>   -4.226</td> <td> 0.000</td> <td>   -0.111</td> <td>   -0.038</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>12.683</td> <th>  Durbin-Watson:     </th> <td>   2.476</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.002</td> <th>  Jarque-Bera (JB):  </th> <td>  13.498</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 1.136</td> <th>  Prob(JB):          </th> <td> 0.00117</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 5.374</td> <th>  Cond. No.          </th> <td>1.90e+03</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 1.9e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                Species   R-squared:                       0.766\n",
       "Model:                            OLS   Adj. R-squared:                  0.717\n",
       "Method:                 Least Squares   F-statistic:                     15.70\n",
       "Date:                Tue, 25 Sep 2018   Prob (F-statistic):           6.84e-07\n",
       "Time:                        15:55:21   Log-Likelihood:                -162.54\n",
       "No. Observations:                  30   AIC:                             337.1\n",
       "Df Residuals:                      24   BIC:                             345.5\n",
       "Df Model:                           5                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept      7.0682     19.154      0.369      0.715     -32.464      46.601\n",
       "Area          -0.0239      0.022     -1.068      0.296      -0.070       0.022\n",
       "Elevation      0.3195      0.054      5.953      0.000       0.209       0.430\n",
       "Nearest        0.0091      1.054      0.009      0.993      -2.166       2.185\n",
       "Scruz         -0.2405      0.215     -1.117      0.275      -0.685       0.204\n",
       "Adjacent      -0.0748      0.018     -4.226      0.000      -0.111      -0.038\n",
       "==============================================================================\n",
       "Omnibus:                       12.683   Durbin-Watson:                   2.476\n",
       "Prob(Omnibus):                  0.002   Jarque-Bera (JB):               13.498\n",
       "Skew:                           1.136   Prob(JB):                      0.00117\n",
       "Kurtosis:                       5.374   Cond. No.                     1.90e+03\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 1.9e+03. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lsmod = smf.ols(formula='Species ~ Area + Elevation + Nearest + Scruz  + Adjacent', data=gala).fit()\n",
    "lsmod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the robust regression using default which is Huber T:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Robust linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>Species</td>     <th>  No. Observations:  </th> <td>    30</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                 <td>RLM</td>       <th>  Df Residuals:      </th> <td>    24</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>               <td>IRLS</td>       <th>  Df Model:          </th> <td>     5</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Norm:</th>                <td>HuberT</td>      <th>                     </th>    <td> </td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Scale Est.:</th>            <td>mad</td>       <th>                     </th>    <td> </td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Cov Type:</th>              <td>H1</td>        <th>                     </th>    <td> </td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>           <td>Tue, 25 Sep 2018</td> <th>                     </th>    <td> </td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>               <td>15:55:21</td>     <th>                     </th>    <td> </td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>        <td>20</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>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>     <td>    6.3626</td> <td>   12.366</td> <td>    0.515</td> <td> 0.607</td> <td>  -17.875</td> <td>   30.600</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Area</th>      <td>   -0.0061</td> <td>    0.014</td> <td>   -0.422</td> <td> 0.673</td> <td>   -0.034</td> <td>    0.022</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Elevation</th> <td>    0.2476</td> <td>    0.035</td> <td>    7.146</td> <td> 0.000</td> <td>    0.180</td> <td>    0.315</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Nearest</th>   <td>    0.3590</td> <td>    0.681</td> <td>    0.528</td> <td> 0.598</td> <td>   -0.975</td> <td>    1.693</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Scruz</th>     <td>   -0.1952</td> <td>    0.139</td> <td>   -1.404</td> <td> 0.160</td> <td>   -0.468</td> <td>    0.077</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Adjacent</th>  <td>   -0.0546</td> <td>    0.011</td> <td>   -4.774</td> <td> 0.000</td> <td>   -0.077</td> <td>   -0.032</td>\n",
       "</tr>\n",
       "</table><br/><br/>If the model instance has been used for another fit with different fit<br/>parameters, then the fit options might not be the correct ones anymore ."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                    Robust linear Model Regression Results                    \n",
       "==============================================================================\n",
       "Dep. Variable:                Species   No. Observations:                   30\n",
       "Model:                            RLM   Df Residuals:                       24\n",
       "Method:                          IRLS   Df Model:                            5\n",
       "Norm:                          HuberT                                         \n",
       "Scale Est.:                       mad                                         \n",
       "Cov Type:                          H1                                         \n",
       "Date:                Tue, 25 Sep 2018                                         \n",
       "Time:                        15:55:21                                         \n",
       "No. Iterations:                    20                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          6.3626     12.366      0.515      0.607     -17.875      30.600\n",
       "Area          -0.0061      0.014     -0.422      0.673      -0.034       0.022\n",
       "Elevation      0.2476      0.035      7.146      0.000       0.180       0.315\n",
       "Nearest        0.3590      0.681      0.528      0.598      -0.975       1.693\n",
       "Scruz         -0.1952      0.139     -1.404      0.160      -0.468       0.077\n",
       "Adjacent      -0.0546      0.011     -4.774      0.000      -0.077      -0.032\n",
       "==============================================================================\n",
       "\n",
       "If the model instance has been used for another fit with different fit\n",
       "parameters, then the fit options might not be the correct ones anymore .\n",
       "\"\"\""
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exog = gala.drop('Species',axis=1)\n",
    "exog = sm.add_constant(exog)\n",
    "rlmod = sm.RLM(gala.Species,exog).fit()\n",
    "rlmod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the small weights:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Espanola        0.679642\n",
       "Gardner1        0.661450\n",
       "Gardner2        0.850097\n",
       "Pinta           0.537700\n",
       "SanCristobal    0.414224\n",
       "SantaCruz       0.174601\n",
       "SantaMaria      0.307863\n",
       "dtype: float64"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wts = rlmod.weights\n",
    "wts[wts < 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L1 or LAD fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda/lib/python3.7/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>QuantReg Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>      <td>Species</td>     <th>  Pseudo R-squared:  </th> <td>  0.5064</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>             <td>QuantReg</td>     <th>  Bandwidth:         </th> <td>   56.17</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>          <td>Least Squares</td>  <th>  Sparsity:          </th> <td>   106.7</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>          <td>Tue, 25 Sep 2018</td> <th>  No. Observations:  </th>  <td>    30</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>              <td>15:55:21</td>     <th>  Df Residuals:      </th>  <td>    24</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th> </th>                      <td> </td>        <th>  Df Model:          </th>  <td>     5</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>const</th>     <td>    1.3145</td> <td>   16.754</td> <td>    0.078</td> <td> 0.938</td> <td>  -33.263</td> <td>   35.892</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Area</th>      <td>   -0.0031</td> <td>    0.020</td> <td>   -0.156</td> <td> 0.877</td> <td>   -0.044</td> <td>    0.037</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Elevation</th> <td>    0.2321</td> <td>    0.047</td> <td>    4.945</td> <td> 0.000</td> <td>    0.135</td> <td>    0.329</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Nearest</th>   <td>    0.1637</td> <td>    0.922</td> <td>    0.177</td> <td> 0.861</td> <td>   -1.739</td> <td>    2.067</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Scruz</th>     <td>   -0.1231</td> <td>    0.188</td> <td>   -0.654</td> <td> 0.520</td> <td>   -0.512</td> <td>    0.266</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Adjacent</th>  <td>   -0.0519</td> <td>    0.015</td> <td>   -3.349</td> <td> 0.003</td> <td>   -0.084</td> <td>   -0.020</td>\n",
       "</tr>\n",
       "</table><br/><br/>The condition number is large, 1.9e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                         QuantReg Regression Results                          \n",
       "==============================================================================\n",
       "Dep. Variable:                Species   Pseudo R-squared:               0.5064\n",
       "Model:                       QuantReg   Bandwidth:                       56.17\n",
       "Method:                 Least Squares   Sparsity:                        106.7\n",
       "Date:                Tue, 25 Sep 2018   No. Observations:                   30\n",
       "Time:                        15:55:21   Df Residuals:                       24\n",
       "                                        Df Model:                            5\n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          1.3145     16.754      0.078      0.938     -33.263      35.892\n",
       "Area          -0.0031      0.020     -0.156      0.877      -0.044       0.037\n",
       "Elevation      0.2321      0.047      4.945      0.000       0.135       0.329\n",
       "Nearest        0.1637      0.922      0.177      0.861      -1.739       2.067\n",
       "Scruz         -0.1231      0.188     -0.654      0.520      -0.512       0.266\n",
       "Adjacent      -0.0519      0.015     -3.349      0.003      -0.084      -0.020\n",
       "==============================================================================\n",
       "\n",
       "The condition number is large, 1.9e+03. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l1mod = sm.QuantReg(gala.Species, exog).fit()\n",
    "l1mod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is not the same as least trimmed squares but as close as we can come with `sm.RLM`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Robust linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>Species</td>     <th>  No. Observations:  </th> <td>    30</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                 <td>RLM</td>       <th>  Df Residuals:      </th> <td>    24</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>               <td>IRLS</td>       <th>  Df Model:          </th> <td>     5</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Norm:</th>              <td>TrimmedMean</td>   <th>                     </th>    <td> </td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Scale Est.:</th>            <td>mad</td>       <th>                     </th>    <td> </td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Cov Type:</th>              <td>H1</td>        <th>                     </th>    <td> </td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>           <td>Tue, 25 Sep 2018</td> <th>                     </th>    <td> </td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>               <td>15:55:21</td>     <th>                     </th>    <td> </td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>         <td>5</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>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>     <td>    2.5164</td> <td>    7.657</td> <td>    0.329</td> <td> 0.742</td> <td>  -12.490</td> <td>   17.523</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Area</th>      <td>    0.0023</td> <td>    0.009</td> <td>    0.253</td> <td> 0.800</td> <td>   -0.015</td> <td>    0.020</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Elevation</th> <td>    0.2172</td> <td>    0.021</td> <td>   10.124</td> <td> 0.000</td> <td>    0.175</td> <td>    0.259</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Nearest</th>   <td>    0.6670</td> <td>    0.421</td> <td>    1.583</td> <td> 0.113</td> <td>   -0.159</td> <td>    1.493</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Scruz</th>     <td>   -0.2229</td> <td>    0.086</td> <td>   -2.589</td> <td> 0.010</td> <td>   -0.392</td> <td>   -0.054</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Adjacent</th>  <td>   -0.0443</td> <td>    0.007</td> <td>   -6.261</td> <td> 0.000</td> <td>   -0.058</td> <td>   -0.030</td>\n",
       "</tr>\n",
       "</table><br/><br/>If the model instance has been used for another fit with different fit<br/>parameters, then the fit options might not be the correct ones anymore ."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                    Robust linear Model Regression Results                    \n",
       "==============================================================================\n",
       "Dep. Variable:                Species   No. Observations:                   30\n",
       "Model:                            RLM   Df Residuals:                       24\n",
       "Method:                          IRLS   Df Model:                            5\n",
       "Norm:                     TrimmedMean                                         \n",
       "Scale Est.:                       mad                                         \n",
       "Cov Type:                          H1                                         \n",
       "Date:                Tue, 25 Sep 2018                                         \n",
       "Time:                        15:55:21                                         \n",
       "No. Iterations:                     5                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          2.5164      7.657      0.329      0.742     -12.490      17.523\n",
       "Area           0.0023      0.009      0.253      0.800      -0.015       0.020\n",
       "Elevation      0.2172      0.021     10.124      0.000       0.175       0.259\n",
       "Nearest        0.6670      0.421      1.583      0.113      -0.159       1.493\n",
       "Scruz         -0.2229      0.086     -2.589      0.010      -0.392      -0.054\n",
       "Adjacent      -0.0443      0.007     -6.261      0.000      -0.058      -0.030\n",
       "==============================================================================\n",
       "\n",
       "If the model instance has been used for another fit with different fit\n",
       "parameters, then the fit options might not be the correct ones anymore .\n",
       "\"\"\""
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ltsmod = sm.RLM(gala.Species, exog, M=sm.robust.norms.TrimmedMean()).fit()\n",
    "ltsmod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results come with standard errors so no need to do bootstrapping.\n",
    "\n",
    "Fit the Galapagos model without Isabela:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda/lib/python3.7/site-packages/statsmodels/regression/linear_model.py:719: RuntimeWarning: divide by zero encountered in log\n",
      "  llf += 0.5 * np.sum(np.log(self.weights))\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>WLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>Species</td>     <th>  R-squared:         </th> <td>   0.871</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>WLS</td>       <th>  Adj. R-squared:    </th> <td>   0.845</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   32.53</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 25 Sep 2018</td> <th>  Prob (F-statistic):</th> <td>6.14e-10</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:55:21</td>     <th>  Log-Likelihood:    </th> <td>    -inf</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    30</td>      <th>  AIC:               </th> <td>     inf</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    24</td>      <th>  BIC:               </th> <td>     inf</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     5</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>   22.5861</td> <td>   13.120</td> <td>    1.722</td> <td> 0.098</td> <td>   -4.492</td> <td>   49.664</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Area</th>      <td>    0.2957</td> <td>    0.061</td> <td>    4.884</td> <td> 0.000</td> <td>    0.171</td> <td>    0.421</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Elevation</th> <td>    0.1404</td> <td>    0.049</td> <td>    2.885</td> <td> 0.008</td> <td>    0.040</td> <td>    0.241</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Nearest</th>   <td>   -0.2552</td> <td>    0.706</td> <td>   -0.361</td> <td> 0.721</td> <td>   -1.713</td> <td>    1.203</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Scruz</th>     <td>   -0.0901</td> <td>    0.147</td> <td>   -0.614</td> <td> 0.545</td> <td>   -0.393</td> <td>    0.213</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Adjacent</th>  <td>   -0.0650</td> <td>    0.012</td> <td>   -5.433</td> <td> 0.000</td> <td>   -0.090</td> <td>   -0.040</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>13.716</td> <th>  Durbin-Watson:     </th> <td>   2.459</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.001</td> <th>  Jarque-Bera (JB):  </th> <td>  17.482</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 1.081</td> <th>  Prob(JB):          </th> <td>0.000160</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 6.052</td> <th>  Cond. No.          </th> <td>1.66e+03</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 1.66e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            WLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                Species   R-squared:                       0.871\n",
       "Model:                            WLS   Adj. R-squared:                  0.845\n",
       "Method:                 Least Squares   F-statistic:                     32.53\n",
       "Date:                Tue, 25 Sep 2018   Prob (F-statistic):           6.14e-10\n",
       "Time:                        15:55:21   Log-Likelihood:                   -inf\n",
       "No. Observations:                  30   AIC:                               inf\n",
       "Df Residuals:                      24   BIC:                               inf\n",
       "Df Model:                           5                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     22.5861     13.120      1.722      0.098      -4.492      49.664\n",
       "Area           0.2957      0.061      4.884      0.000       0.171       0.421\n",
       "Elevation      0.1404      0.049      2.885      0.008       0.040       0.241\n",
       "Nearest       -0.2552      0.706     -0.361      0.721      -1.713       1.203\n",
       "Scruz         -0.0901      0.147     -0.614      0.545      -0.393       0.213\n",
       "Adjacent      -0.0650      0.012     -5.433      0.000      -0.090      -0.040\n",
       "==============================================================================\n",
       "Omnibus:                       13.716   Durbin-Watson:                   2.459\n",
       "Prob(Omnibus):                  0.001   Jarque-Bera (JB):               17.482\n",
       "Skew:                           1.081   Prob(JB):                     0.000160\n",
       "Kurtosis:                       6.052   Cond. No.                     1.66e+03\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 1.66e+03. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wts = gala.index != 'Isabela'\n",
    "limod = smf.wls(formula='Species ~ Area + Elevation + Nearest + Scruz  + Adjacent', data=gala, weights=wts).fit()\n",
    "limod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show example with multiple outliers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "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>temp</th>\n",
       "      <th>light</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4.37</td>\n",
       "      <td>5.23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.56</td>\n",
       "      <td>5.74</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.26</td>\n",
       "      <td>4.93</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.56</td>\n",
       "      <td>5.74</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.30</td>\n",
       "      <td>5.19</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   temp  light\n",
       "0  4.37   5.23\n",
       "1  4.56   5.74\n",
       "2  4.26   4.93\n",
       "3  4.56   5.74\n",
       "4  4.30   5.19"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "star = pd.read_csv(\"data/star.csv\")\n",
    "star.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XlwHOd55/HvM7gvAiQOHuABggdEgjooQaJleRNF3ood2VGU26pSZZ31RpVUNnFVEnutimsr600qdpSqHLUpH+tU4o1i2YqkyKcsH5KcSJGUkKZsGuApAiRuzOAgiBuYefePGQxmQBwDYoCZ7vl9qlDCdDdmHkKYBy+efvppc84hIiL+Esh0ACIikn5K7iIiPqTkLiLiQ0ruIiI+pOQuIuJDSu4iIj6k5C4i4kNK7iIiPqTkLiLiQ/mZeuGamhrX0NCQqZcXEfGkU6dOhZxztasdl7Hk3tDQwMmTJzP18iIinmRmV1I5TmUZEREfUnIXEfEhJXcRER9SchcR8SEldxERH1JyFxHxoYy1QmaL509388SL5+kZmWRXVQkfeU8TDx+vz3RYIiLrktPJ/fnT3Tz+3BkmZ8MAdI9M8vhzZwCU4EXE03K6LPPEi+fjiX3e5GyYJ148n6GIRETSI6eTe8/I5Jq2i4h4RU4n911VJWvaLiLiFTmd3D/yniZKCvKStpUU5PGR9zRlKCIRkfTI6ROq8ydN1S0jIn6T08kdogleyVxE/CanyzIiIn6l5C4i4kNK7iIiPqTkLiLiQ0ruIiI+pOQuIuJDSu4iIj6k5C4i4kNK7iIiPqTkLiLiQ0ruIiI+pOQuIuJDOTs4TPdOFRE/y8nkrnuniojf5WRZRvdOFRG/y8nkrnuniojfpZTczazKzJ4xs3NmdtbM7l2038zsr83skpn9yMzu3Jhw00P3ThURv0t15f5XwLecc7cAtwNnF+3/GeBQ7OMx4NNpi3AD6N6pIuJ3q55QNbMtwE8AHwRwzs0AM4sO+zng/znnHPBGbKW/0znXm+Z400L3ThURv0ulW6YRCAJ/Z2a3A6eADzvnxhOOqQc6Ex53xbZlZXIH3TtVRPwtlbJMPnAn8Gnn3HFgHPjYomNsia9zizeY2WNmdtLMTgaDwTUHKyIiqUkluXcBXc65N2OPnyGa7Bcfsyfh8W6gZ/ETOec+55xrcc611NbW3ky8IiKSglWTu3OuD+g0s/mzje8G2hYd9lXg12JdM+8ArmVrvV1EJBekeoXq7wD/aGaFwGXg183sNwGcc58Bvgk8CFwCJoBf34BYRUQkRSkld+fcW0DLos2fSdjvgN9OY1wiIrIOOXmFqoiI3+Xk4LBEmg4pIn6U08ld0yFFxK9yuiyj6ZAi4lc5ndw1HVJE/Cqnk7umQ4qIX+V0ctd0SBHxq5w+oarpkCLiVzmd3EHTIUXEn3K6LCMi4ldK7iIiPqTkLiLiQ0ruIiI+lPMnVEUkN3z8+TM89WYnYefIM+ORE3v444dvzXRYGybnk7sGh4n438efP8OTb1yNPw47F3/s1wSf02WZ+cFh3SOTOBYGhz1/ujvToYlIGj31ZueatvtBTid3DQ4TyQ1h59a03Q9yOrlrcJhIbsgzW9N2P8jp5K7BYSK54ZETe9a03Q9yOrlrcJhIbvjjh2/l0Xfsja/U88x49B17fXsyFcBchmpOLS0t7uTJkxl57UTqlhERLzGzU865ltWOy/lWSA0OExE/yumyjIiIX+X8yl1EZCVeLd0quYuILGP+Qsf562HmL3QEsj7BqywjIrIML1/oqOQuIrIML1/oqLKMiM9ke4042+NLtKuqhO4lErkXLnT03Mp9YGKAFztepOt6F5nq0RfJVtk+DC/b41vMyxc6em7l/kbvG/zhq38IQGVRJc3VzQsfNc1sL92O+XhehMhKVqoRZ8PqONvjW2w+Jq/8pZHIc8n9vQ3v5UDVAVpDrbQNttE62Mrf/fjvmHNzAGwr3hZP9PNJv7a0NsNRi2yObK8RZyK+9ZaBvHqho+eSe2FeYTxpz5uam+LC8IV4sm8dbOW1H71GxEUAqCup42jNUZqrmzlaHf1vdUl1pv4JIhsm22vEmx2fl1sZ1yul5G5mHcB1IAzMLZ5rYGb3A18B2mObnnPOfSJ9Ya6sOL+Y22pv47ba2+LbJmYnuDB8IZrsQ9GE//3O7+OI1ul3lO1IKukcrT5KVXHVZoUssiE+8p6mpGQG2VUj3uz4vFYGSqe1rNx/yjkXWmH/vzrn3r/egNKltKCUO+ru4I66O+LbxmfHOTt4Nr66bxts43tXvxffX19en1TSOVJ9hC2FWzIRvshNyfYa8WbHl+1lqo3kubLMepQVlNGyo4WWHQt/eIzOjC4k/NgK/9tXvh3fv2/LvngpZz7hlxWUZSJ8kZRke414M+PL9jLVRko1uTvg22bmgM865z63xDH3mtkPgR7gD5xzrekKciNtKdzCiZ0nOLHzRHzbyNQIbYNttA210Rpq5a2Bt3ih/QUADGN/5f6FhF/TTNPWJkoLSjP1TxCRZWR7mWojpTTP3cx2Oed6zKwO+A7wO865f0nYvwWIOOfGzOxB4K+cc4eWeJ7HgMcA9u7de9eVK1fS9e/YcIOTg0knbNtCbQxMDgAQsACNlY1JJZ2mbU0U5RVlOGoRb9mIC5y8dNFUKlKd577mm3WY2R8BY865P1/hmA6gZaUafbbcrGM9BiYGFhJ+rKQzNDUEQL7lc3DrwYUOnZpmDlcdpiCvIMNRi2TWcsl2cWcLRFfZf/oLt3o6Gadb2m7WYWZlQMA5dz32+U8Dn1h0zA6g3znnzOweole+Dt5c6N5RV1pHXWkd9++5HwDnHP0T/fFE3zrYynevfpdnLz4LQEGggMNbDyet8BurGikIKOFLblipNTGXO1s2Qio19+3AP8eu+swHvuic+5aZ/SaAc+4zwC8Bv2Vmc8Ak8AGXg7MBzIwdZTvYUbaDd+97NxBN+N1j3fHunNbBVl5of4GnLzwNQFFeEU1bm+Kr++bqZhorG8kL5K30UiKetFICz+XOlo2Q8/dQzYSIi9B1vSupnNM22MbE3AQAJfkl3LLtlqSSTsOWBgLmuVFAIkn2f+wbLJVxjOU7W/LMiDjni3p5OmxYzT1dcjm5LyXiInSMdiSNVTg3dI7JuegPe1lBGUe2HUkq6eyp2KM5OuIp933ypSUTeH0scS+uuS+mGryfk/v0GEyPQsVO8Hlim4vM0X6tPb7Cbxts49zQOWYiMwBUFFYk9eA31zSzq2yXEr5krdVOmiaebA2YEV4iP9VXlfDaxx7YzLCzin+Te9tX4Olfg8JyqD4A1Yeg5hBUH4z+d9sBKCpPf8BZYjYyy9sjbyet8M8Pn2cuEh2cVlVUdUPC16RMySaptiauVMJp/+T7NjzObOXb5N7bcY6e//gK9ZEetk5eoXDkbWykExJ/DCp2Qc3BhMR/KPq4cg/48ETlTHiGiyMXkxL+xeGLhF10daRJmeJFK5VwtHJPQytktnl9qJzfO3UMOAZAeVE+TdX53LVlhFuLBmi0HnbOdbFlvIP8M8/A9LWFL84riq32DyYk/diqv8S7Q8NWmpSZeNL2tZ6lJ2XOn7jVpEzJJrl8dWk6eG7lPheO0DMyRfvgOO3BMdpD47QPTtAeGqN7eJJIwj9nW2kBd2yb5a6yEEeL+tnneqid7qRsrJ3AcAe4hBM3ZbULK/zEFf/WfeCTC48mZic4P3w+qQ+/41pHfFLmzrKd8RX+0W1HNSlTNsVKZZql9kH2DkbbDL4ty6xkei5M59AEl4PjdAyORxN/7KN/dDrp2N0VeZyoGuX20hBN+X3sjnSzbfIKxaOXsYmE668C+bB1/8IKP3HVX1bj+ZO6YzNjnB06Gy/ntA22cWV0YSyEJmXKRlrrVam6ijVHk/tKxqfn6BgcpyMUXeVfDo3TEUv8wxOz8eMCBkeqwpzYMsStxUEOBHrZOddF1XgH+dfasfDMwpMWV954Qrf6EGxrhILiTfu3pdtSkzK7xxbucalJmZIua62rqw7v45r7zSoryqd5VyXNuypv2DcyMUN7KLbaD45zOTTOvw9u5cvtOxifORo/rjgP7t46Rkt5iKOFAzTQQ910J+WXXibvh08lPKNB1d7kk7nzvwQ80MK50qTM+XLO6YHTN0zKXDw4rSTf/2NVZX3WelWqrmJNXc4k95VUlRZyfG8hx/duTdrunCM4Nk17MFbeiSX/bwzu4G+6JpiZi8SPrSmc4Z2VIxwvC9KU38+eSDfVQ1coufJv2OzEwpMu18JZfRAKs3f1W1VcxTvr38k7698Z3xaaDC2Uc0JtvN77Ol+7/DUgOinzQNWBpBO2mpQpi6113nouz2dfq5wpy6RbOOLovTaZVNdvj5V6OocnCcfO7BoRDpdc58SWIW4rDnIwr49dc11UTVyhYKwLS2zh3FK/KPF7r4VzYGIgaaSCJmXKSlRzXzvV3DNoZi5C1/DEkom/59pU/LgiZrizfIiW8kGOFUXLPNtnO6kY6yBvZnThCT3cwrnUpMzWwVauxVpUNSlT1jpvfb3z2b0+313JPUtNzoS5MrRQ258/qdsxOE5obP5kraPWRrm7YpC7ykLcUhAt89RMX6VkrBPzeAtn4qTM+ZJO22Ab12evA7FJmduakm5gvr9yvyZlyrr5YeWv5O5Bo1Oz8WSf1M4ZHOf6dHS8QAFz7A8McM+WIW4viZZ56sNdbJ24QsH00MKTLW7hnE/61QezsoUz4iJ0Xu9MWuGfHTybNCnzyLYjSaOR923Zp0mZsiZ+6LZRcvcR5xyD4zN0hKKr/fZFK/6p2eiJ3UrGuKWgn7sromWe/dbDjtkuKsavEoh4r4UzHAlzZfRKUv3+7OBZpsLR0pYmZcpa+WFejVohfcTMqCkvoqa8iJaGbUn7IhFH3+hUUm3/bGicb4bGuTo0wVzEESBCvQU5VjQQL/PsG+uhZvB7lExlbwtnXiCPxqpGGqsa+dkDPwssPSnzqXNPaVJmBqSrdv3x58/w1JudhJ0jz4xHTuzhjx++dQMizq1uG63cfWwuHKFrOLmjp2MwWvLpuTaJc1DGJPutl9tLgtxREuRQXh/1kW62Tl4lP5zwJsjiFs7ESZnzJZ0LwxeSJmUmdug0V2tS5nqlq3b98efP8OQbV2/Y/ug79m5IglfNfRMouWfW1GyYq4mjGoKxPv7QOMHr04BjB0M0Bnq5szTIseJgdCjbbBflU71LtHAuGs2Q4RbOmfAMF4cvLnTohFq5NHIpPimzurg6OkMnYZWvSZmpS1ft+sDj31xyZnueGW//6YPrinE5udIto7JMjiouyOPw9goOb6+4Yd/Y9Fy8pj//8Wrsv9cmZylihgbr42Cgj+OlQY6E+9nX30Nd5w8onLu+8EQZbOEszCuMrtJrkidlzg9Om6/hv9r96rKTMptrmtlWvG25l/C09Sa4dF0pulRiX2l7Ojx8vN5TyfxmKbnLDcqL8jlWX8mx+htHNQyPzySd1D0dGue52OPJ2TlqGKXRejic38fthUGaxvrYPXqaqnPfIJDhFs7i/GJur72d22tvj29balLm9zu/v/SkzNgqv7Loxu+LlywuTXSPTPL4c2cAUk566apd5y1zt6W8myyZeX1Vnk5K7rImW8sKuauskLv23TiqYeD6NJeDC7X9F4PjfHZwnCuD4xCeZa/1s9/6OFLQz7G5AQ4Ge9nZ8zVK54YXnmiTWzhLC0o5Xnec43XH49uSJmXGkv53r343vn93+e6kG58cqT5CReGNfwFlqydePH/DfUonZ8M88eL5lBNhumatP3Jiz5I190dO7FnT80B6fmn5iWrusuHCEUfPyGR0xR8co2NwIn4BV9fwBOVujAPWS6P10lzUz5HCfhropW6mizy3MLEzky2c16avcXbobNJoBa9OykxXO2C2dcv4oYc9FTqhKp4wP4O/PTaKef6/HaEJBkYnqLdgPPEfKxrgcH4feyLdVM6F4s/hLIBV7tn0Fs7hqeGFSZmhVtqG2ugb7wOye1KmX5OgH3rYU6ETquIJRfl5HKyr4GBdBbA9ad/iGfyvhSZ4MhRd+U9PXWN/LOkfDPRy9HqQg+Md7Lr8KoWRhfk9rrAc26AWzq3FW7mv/j7uq78vvm0tkzKbq5s5vO3wpk/K9Ovt63Kphz0VWrmLJ12bmI21bo7F2jjnV/zjlE8P0BiIJv5DgV6OFvbTYL3UzPUntXC6LfXYJrRwJk7KnF/lD09HzzPkWz6Hth5K6sE/VHVowydlpqOkkm0nL/3Qw54KlWUkJyXO4O8YTB7O1js4wq5wD42xFf/h/F5uye9nj+umNDK+8Bx5RbHV/sa0cDrn6BvvS0r2rYOtjMYmgXphUma2JtKV6vfZ9svoZim5iywSiTh6YjP45+f0dMRO8k6M9NHgeuIr/qb8Pg7m9bEj3EceCwksUlpLoCb9LZzOObrGuuLlnPmTtmOzY0D2TcrMxrr9Sr9wgKz8ZXQzlNxF1mA2HImd2E0e1XB14BoF16/GVvvRVX9TQT+N1kNl5Fr8610gH7d1fzTxp6mFc6lJmW2DbUzORZNqJidlZuPJy5V+4QBZ98voZumEqsgaFOQFaKwtp7G2/IZ9iTP42wfH+WKs5BMK9lM1cSWa+AM9NA700jR0hj3uOxSw0MIZLqokUHMIW2MLZ8AC7Nuyj31b9vFgY/RS/MRJmfMlnWcuPMOTZ58EopMykwanVTezu2J32ufoZOPJy5u5atbP915VchdZRUlhHrfs2MItO7bcsC9xBv/50DjfCo1zNTjK1OAVts90RhP/XA8HJns51PMiNW5hCqezAOGK3eTVHo4m/hRaOJeblHn52uWkHvwvnv3ihk7KTGfHTbpq4av9wsm2X0YbTWUZkQ3gnGNofCZ645WEk7p9wSCBobepD/dEV/vWy4FALwesl2Km418fzi/DVR9YSPxrbOGcDc9yaeRS0gr/4sjFtE7KTFfHTbpq4aq5LzpOyV1kc0Uijv7rU8m3WgyOcT10laKRy+yjhwOx+v6BvF52ESKQUOGeKdtJoOYQ+bVru5F6qpMyE1f4NSU1G/q9SPeJ2ZV+4ahbZpMouYvcKD6DPzaGuWNwnK6BIeZCb1M+1s5+ovX9A7EVfwUT8a8NBwqZrdpPQd1h8moPp9TCmTgpc76kc/na5YVJmaV18WQ/v8pfPClzPUkzG0/MZru0nlA1sw7gOhAG5hY/sUX/lvsr4EFgAvigc+4Haw1aJNfl5wVoqCmjoaaMn0oqX/9EfAZ/e2icN0LjfCk4xuBAN4HBi2yduhq9YjfYS2PoJPsC3yCfSPyrZ4qqCW87SNGOplgrZ3TFX7x135KTMs8NnUta4b/c+XJ8f+KkzJHh7fzt92aZnI5eZbvWYV3ZeGLWL9ZyQvWnnHOhZfb9DHAo9nEC+HTsvyKSJkvP4I8m5cQZ/N8IjXM1eI2J/rfJH77EjtlOGud6aZzs5UDP81TbaPyrI5bPVMVerOYQxdubsNpDlFYf4s7qg9xZdzx+Und+UmZiW+b8pMz8Riib2UZ4ajfhyXpmpnbzZ98Op5Tc/ToKIRukVJaJrdxblkvuZvZZ4BXnoq0AZnYeuN8517vcc6osI7I55mfwx0/qDvQRHrhA8bXL7I50xa/Y3W99FNpc/Otm8iuYqTpAXt1hinc0xU7sLrRwXpu+xl2f+nsCJV0EirvJK+4iULgwvrlhS8NCl05NM0e2HaG0oPSG+DbzHqp+kO4+dwd828wc8Fnn3OcW7a8HOhMed8W2LZvcRWRz3DiDvwn4yaQZ/D8YHOe54CijfZchdInysXb2zXWzv7+PxuB32dn2dPz5IgSYLN2FbT3An1DOmZE63nZ3cDnyMwzkFRIo7mHr1n4a90xyqv8U32z/JhCdlNlY2ZjUoXOhs4JnT3XHb9gRdo5nT3XTsm+bJ092ZpNUk/t9zrkeM6sDvmNm55xz/5Kwf6n+qRv+JDCzx4DHAPbu3bvmYEUkfcyM7VuK2b6lmHsPVMe2HgOSZ/B/OzROd3+Qqf4L5A1domqig/3Xe2kcu8L7rZdfLVho4RxzxXTM7aTQDtPg7qDwlv9EqKKWNjdN67W3aR1s5d96/i0+KRMXwHbXUTxVHyvr7GZyeseabhwiS1tzt4yZ/REw5pz784RtKsuI5IjoDP7YjJ7gGG+1tTHWc5Y9kYUWzkbrpd5CBGwhv4wVbWe6spFA7SEmt+/hXFEBH321lfHiYQLF3QTyo8PbnAsQmd7Br9x676ZOyvSKtLVCmlkZEHDOXY99/h3gE865byUc8z7gvxPtljkB/LVz7p6VnlfJXWRjZLKfe2JmLjZ/f5zOgUGu91yAwYuUjrazY64zfuOVLbbQwjntCrjsdvCDvDpOFZZzsSjAcNkMbssIozPRG64XBApo2toUT/ZHq49yoOoA+YHcu8g+ncm9Efjn2MN84IvOuT8xs98EcM59JtYK+X+A9xJthfx159yKmVvJXST9snUULyTM4A9ep7+3i6m+C0z3naNq8kp8tb/P+sm3aAunAy4UbeOHlTs5X17OpULHhblhxsLRMlC2TcrcLLqISSQHZeMo3tX88w+6+NS3ztM3OkVNsfHT9VPUh7soHLlM1UQH+4iWempslAhwNT+fM0XFnC7dxtmSYi7lzTIVG8tcklfMkeojHE2Yhb9ZkzI3i6ZCiuSgm5mMuF7rLQP9/J27+fk7dy+5b34G/7nQBN29PYx1n8OFLlAy2s4916/yq/Swx/roLYTWoiJaCws5Mz7GP/Wf5slYm0dZoIijVQdo3n4nzbW3b9ikzGyj5C7iI5t9xefiMtBar1BdTSBg7N5ayu6tpXCoBrgtvm9+Bv8bwVGCXZeY6jtH/eDbNA10UDNzlfyiAfqLJ2ktKqRtYpR/HGxl1mKjkcmnqWg7t9Uc4bb6d9C8+13sLF/fpMxso7KMiI9sds09W8tA8zP4O/sGGO48y1z/BWz4IpNzl7ieP0CoeJzzRflcLCxgLpbQKyPQ6Mo4VFzPsZpj3NP4k+zacw9WdOOM/0xSWUYkB80n8M3qlslEGSgVSTP47ziYtG90apaO4Bi9XZcZ7jpD39ApgtOXCAUG6C28zrMz53m69wL0Pkf1XJjDs7DPVdBQtJdbam6jYe9dbN17lEBVem+knm5auYvITdvMlftGt3jOz+A/393Nmcsv0R46Sc/UZXosxED+LC5Wsambm6N5eoZbZsLUhyvZWbiX8qominc2sW1fM5W7j2AlW1d+sXVQt4yIbLjNKgNlusVzbHqc17t+yL93vM6F4A/onOogmDCAbefsHMdmZjg6PUPz9DS7Z0uYLtzL5Jb9WM0hynYdoaahmYodB9d1I3VQcheRTbIZF01lY23/+sx1zg2d40zwx5zsfYtzoTMEZ4Px/XWzxtHpae6cGad5eoYj0zOUuAD9eTvpv+XXuPOX/8dNva5q7iKyKR4+Xr/hq+dsrO1XFFZw9467uXvH3fzX2BDLa9PXaBtsi9/45MehH/PK+MIUlrpICfunjbvCA9y5wfEpuYtI1vPKTT0qiyq5d9e93Lvr3vi2oamhaMJPmIV/T9PGz6tXcheRrOflm3psK97Gu+rfxbvq3xXfFo6EV/iK9FByF5Gst9ktnhttM+bfKLmLiCdsRm3fT/wzTUdEROKU3EVEfEjJXUTEh5TcRUR8SMldRMSHlNxFRHxIyV1ExIeU3EVEfEgXMYnkkM2Y4JhNr5vLlNwlLfTmzX4bfb/TbHvdXKeyjKzb/Ju3e2QSx8Kb9/nT3ZkOTRI88eL5pMFbAJOzYZ548bwvXzfXaeUu67bSm1crs+yxUTPRV/urLRtnsecCrdxl3fTm9YblZp+vZyZ6Kn+1bcTryuqU3GXd9Ob1ho+8p4mSguRRs+udiZ5KyWUjXldWp+Qu66Y3rzc8fLyeP/2FW6mvKsGI3n90vTeYTuWvto14XVmdau6ybn67kYKfpXsmulduf5eLlNwlLXQjhdyUyu3v1AqZGSrLiMhNS6XkolbIzNDKXUTWZbW/2tRNlRlauYvIhlI3VWYouYvIhlI3VWaknNzNLM/MTpvZ15fY90EzC5rZW7GP/5beMEXEqx4+Xs8v3lVPnhkAeWb84l06Ab/R1rJy/zBwdoX9X3bO3RH7+Pw64xIRn3j+dDfPnuom7BwAYed49lS3Zg9tsJSSu5ntBt4HKGmLyJqoWyYzUl25/yXwUSCywjG/aGY/MrNnzGzPUgeY2WNmdtLMTgaDwbXGKiIepG6ZzFg1uZvZ+4EB59ypFQ77GtDgnLsN+C7whaUOcs59zjnX4pxrqa2tvamARcRb1C2TGams3O8DHjKzDuBLwANm9mTiAc65QefcdOzh/wXuSmuUIuJZ6pbJjFWTu3PucefcbudcA/AB4CXn3KOJx5jZzoSHD7HyiVcRySEaHJYZN32Fqpl9AjjpnPsq8Ltm9hAwBwwBH0xPeCKS7XSLxexkLtaetNlaWlrcyZMnM/LaIpIei4eCQbTkkrgyT+UYSZ2ZnXLOtax2nK5QFfGI5093c98nX2L/x77BfZ98KSv6xFNpc1QrZGZocJiIB6RrbG66SyiptDmqFTIztHIX8YB0rH5Tud/pWqXS5qhWyMxQchfxgHSsfjeiPJJKm6NaITNDZRkRD0jH7ew2ojySyi0WdRvGzFC3jIgHLNVxYoAj2jeeSrK875MvLfkLor6qhNc+9kCaI5aNom4ZER9JvBAIFhI7pF47V3kktyi5i3jEw8free1jD1BfVcLiv7dTqZ3rStHcopq7iMesp3a+2v1OxT+0chfxGLUWSiqU3EU8RrVzSYXKMiIes57WQg35yh1K7iIedDO183SNMBBvUFlGJEdogFduUXIXyREa4JVblNxFckRlScGatou3KbmL5AiztW0Xb1NyF8kRIxOza9ou3qZuGREPupmWxnRMlrxZasHcfFq5i3jMzd50I1MXP23ETUJkdUruIh5zsy2NmRocphbMzFBZRsRjvDY4TC2YmaGVu4jHeG1wmNfi9QsldxGP8drgMK/F6xcqy4h4jNfuSeq1eP1C91AVEfEQ3UNVRCSHKbmLiPiQkruIiA/phKrkBF3+LrlGyV18T3cgklwUBJXAAAAFzUlEQVSksoz4ni5/l1yUcnI3szwzO21mX19iX5GZfdnMLpnZm2bWkM4gRdZDl79LLlrLyv3DwNll9n0IGHbOHQT+AvjUegMTSRdd/i65KKXkbma7gfcBn1/mkJ8DvhD7/Bng3Wa6v4tkB13+Lrko1ROqfwl8FKhYZn890AngnJszs2tANRBad4Qi66TL3yUXrZrczez9wIBz7pSZ3b/cYUtsu2GugZk9BjwGsHfv3jWEKbI+mRh1K5JJqZRl7gMeMrMO4EvAA2b25KJjuoA9AGaWD1QCQ4ufyDn3Oedci3Oupba2dl2Bi4jI8lZduTvnHgceB4it3P/AOffoosO+CvwX4HXgl4CXXKYmkonkAF2UJau56YuYzOwTwEnn3FeBvwX+wcwuEV2xfyBN8YnIIrooS1KxpuTunHsFeCX2+f9M2D4F/HI6AxORpa10UZaSu8zTFaoiHqOLsiQVSu4iHqOLsiQVSu4iHqOLsiQVmgop4jG6KEtSoeQu4kG6KEtWo7KMiIgPKbmLiPiQkruIiA8puYuI+JCSu4iIDym5i4j4kGVqeKOZBYErN/nlNXj3RiCKPTMUe2Z4NfZsjnufc27VmekZS+7rYWYnnXMtmY7jZij2zFDsmeHV2L0adyKVZUREfEjJXUTEh7ya3D+X6QDWQbFnhmLPDK/G7tW44zxZcxcRkZV5deUuIiIryNrkbmbFZvbvZvZDM2s1s/+1wrG/ZGbOzLLi7HaqsZvZr5hZW+yYL252nEtJJXYz22tmL5vZaTP7kZk9mIlYl2JmebG4vr7EviIz+7KZXTKzN82sYfMjXN4qsf9e7GflR2b2PTPbl4kYl7NS7AnHZNX7dN5qsWfj+zQV2Tzydxp4wDk3ZmYFwKtm9oJz7o3Eg8ysAvhd4M1MBLmMVWM3s0PA48B9zrlhM6vLVLCLpPJ9/zjwtHPu02Z2FPgm0JCBWJfyYeAssGWJfR8Chp1zB83sA8CngF/dzOBWsVLsp4EW59yEmf0W8Gd4J/ZsfZ/OWzb2LH6friprV+4uaiz2sCD2sdQJgv9N9Ad9arNiW02Ksf8G8DfOueHY1wxsYojLSjF2x8IboRLo2aTwVmRmu4H3AZ9f5pCfA74Q+/wZ4N1mZpsR22pWi90597JzbiL28A1g92bFtpoUvu+Qhe9TSCn2rHyfpiJrkzvE/1x6CxgAvuOce3PR/uPAHufcsn8KZspqsQOHgcNm9pqZvWFm7938KJeWQux/BDxqZl1EV+2/s8khLucvgY8CkWX21wOdAM65OeAaUL05oa1qtdgTfQh4YWPDWZMVY8/m9ymrf9+z9n26mqxO7s65sHPuDqKrlHvM7Nj8PjMLAH8B/H6m4lvJSrHH5AOHgPuBR4DPm1nV5ka5tBRifwT4e+fcbuBB4B9i/z8yxszeDww4506tdNgS2zLeLpZi7PPHPgq0AE9seGApWC32bH6fpvh9z9r36WqyOrnPc86NAK8Aib81K4BjwCtm1gG8A/hqtp2sWSZ2gC7gK865WedcO3Ce6A9R1lgh9g8BT8eOeR0oJjqLI5PuAx6K/Sx8CXjAzJ5cdEwXsAfAzPKJlpSGNjPIZaQSO2b2n4E/BB5yzk1vbojLWi32bH6fpvozk9Xv02U557LyA6gFqmKflwD/Crx/heNfIXrCyROxE02YX4h9XkO0XFDtkdhfAD4Y+/wI0Zq7ZTr2hPjuB76+xPbfBj4T+/wDRE8KZzzeFGM/DrwNHMp0jGuNfdExWfM+TfH7npXv01Q+snnlvhN42cx+BPwH0drv183sE2b2UIZjW00qsb8IDJpZG/Ay8BHn3GCG4k2USuy/D/yGmf0QeIpoos94eWMpi+L+W6DazC4Bvwd8LHORrW5R7E8A5cA/mdlbZvbVDIa2Ko+8T5fkkffpqnSFqoiID2Xzyl1ERG6SkruIiA8puYuI+JCSu4iIDym5i4j4kJK7iIgPKbmLiPiQkruIiA/9fxX9I7QFc+O3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gs1 = smf.ols(formula='light ~ temp', data=star).fit()\n",
    "gs2 = smf.RLM(star.light, sm.add_constant(star.temp), data=star).fit()\n",
    "gs3 = smf.RLM(star.light, sm.add_constant(star.temp), data=star,M=sm.robust.norms.TrimmedMean(c=0.1) ).fit()\n",
    "plt.scatter(star.temp, star.light)\n",
    "plt.plot([3.4, 4.7], [gs1.params[0] + gs1.params[1]*3.4, gs1.params[0] + gs1.params[1]*4.7])\n",
    "plt.plot([3.4, 4.7], [gs2.params[0] + gs2.params[1]*3.4, gs2.params[0] + gs2.params[1]*4.7])\n",
    "plt.plot([3.4, 4.7], [gs3.params[0] + gs3.params[1]*3.4, gs3.params[0] + gs3.params[1]*4.7])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Does not work. This class of estimators cannot do the job for this problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/json": {
       "Software versions": [
        {
         "module": "Python",
         "version": "3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)]"
        },
        {
         "module": "IPython",
         "version": "6.5.0"
        },
        {
         "module": "OS",
         "version": "Darwin 17.7.0 x86_64 i386 64bit"
        },
        {
         "module": "pandas",
         "version": "0.23.4"
        },
        {
         "module": "numpy",
         "version": "1.15.1"
        },
        {
         "module": "matplotlib",
         "version": "2.2.3"
        },
        {
         "module": "seaborn",
         "version": "0.9.0"
        },
        {
         "module": "scipy",
         "version": "1.1.0"
        },
        {
         "module": "patsy",
         "version": "0.5.0"
        },
        {
         "module": "statsmodels",
         "version": "0.9.0"
        }
       ]
      },
      "text/html": [
       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)]</td></tr><tr><td>IPython</td><td>6.5.0</td></tr><tr><td>OS</td><td>Darwin 17.7.0 x86_64 i386 64bit</td></tr><tr><td>pandas</td><td>0.23.4</td></tr><tr><td>numpy</td><td>1.15.1</td></tr><tr><td>matplotlib</td><td>2.2.3</td></tr><tr><td>seaborn</td><td>0.9.0</td></tr><tr><td>scipy</td><td>1.1.0</td></tr><tr><td>patsy</td><td>0.5.0</td></tr><tr><td>statsmodels</td><td>0.9.0</td></tr><tr><td colspan='2'>Tue Sep 25 15:55:21 2018 BST</td></tr></table>"
      ],
      "text/latex": [
       "\\begin{tabular}{|l|l|}\\hline\n",
       "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n",
       "Python & 3.7.0 64bit [Clang 4.0.1 (tags/RELEASE\\_401/final)] \\\\ \\hline\n",
       "IPython & 6.5.0 \\\\ \\hline\n",
       "OS & Darwin 17.7.0 x86\\_64 i386 64bit \\\\ \\hline\n",
       "pandas & 0.23.4 \\\\ \\hline\n",
       "numpy & 1.15.1 \\\\ \\hline\n",
       "matplotlib & 2.2.3 \\\\ \\hline\n",
       "seaborn & 0.9.0 \\\\ \\hline\n",
       "scipy & 1.1.0 \\\\ \\hline\n",
       "patsy & 0.5.0 \\\\ \\hline\n",
       "statsmodels & 0.9.0 \\\\ \\hline\n",
       "\\hline \\multicolumn{2}{|l|}{Tue Sep 25 15:55:21 2018 BST} \\\\ \\hline\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "Software versions\n",
       "Python 3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)]\n",
       "IPython 6.5.0\n",
       "OS Darwin 17.7.0 x86_64 i386 64bit\n",
       "pandas 0.23.4\n",
       "numpy 1.15.1\n",
       "matplotlib 2.2.3\n",
       "seaborn 0.9.0\n",
       "scipy 1.1.0\n",
       "patsy 0.5.0\n",
       "statsmodels 0.9.0\n",
       "Tue Sep 25 15:55:21 2018 BST"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%load_ext version_information\n",
    "%version_information pandas, numpy, matplotlib, seaborn, scipy, patsy, statsmodels"
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}