2. Linear Regression/4. Linear Regression with Statsmodels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Linear Regression**: Human Brain Weight 🧠"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This dataset was compiled using a medical study conducted on a group of people. This dataset shows a few variations of head sizes and masses of brains, it also consists additional gender and age group columns. The goal of the study is to predict the masses of the brains using the data, and later compare it with the actual masses of brains mentioned in the dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Linear Regression with Statsmodels**\n",
    "\n",
    "$\\qquad$ <span style=\"color:gray\"><b>0.</b> Settings </span><br>\n",
    "$\\qquad$ <span style=\"color:gray\"><b>1.</b> Dataset </span><br>\n",
    "$\\qquad$ <span style=\"color:gray\"><b>2.</b> Data Preprocessing </span><br>\n",
    "$\\qquad$ <span style=\"color:gray\"><b>3.</b> Data Preparation </span><br>\n",
    "$\\qquad\\qquad$ <span style=\"color:gray\"><b>3.1.</b> One-Hot-Encoding </span><br>\n",
    "$\\qquad\\qquad$ <span style=\"color:gray\"><b>3.2.</b> Ordinal-Encoding </span><br>\n",
    "$\\qquad\\qquad$ <span style=\"color:gray\"><b>3.3.</b> Update </span><br>\n",
    "$\\qquad\\qquad$ <span style=\"color:gray\"><b>3.4.</b> Split in train and validation </span><br>\n",
    "$\\qquad$ <span style=\"color:gray\"><b>4.</b> Linear Regression with Statsmodels </span><br>\n",
    "$\\qquad\\qquad$ <span style=\"color:gray\"><b>4.1.</b> Basic </span><br>\n",
    "$\\qquad\\qquad$ <span style=\"color:gray\"><b>4.2.</b> Advanced </span><br>\n",
    "$\\qquad\\qquad$ <span style=\"color:gray\"><b>4.3.</b> Fancy </span><br>\n",
    "$\\qquad\\qquad$ <span style=\"color:gray\"><b>4.4.</b> GLM </span><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **0.** Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Libraries\n",
    "from statsmodels.graphics.regressionplots import plot_leverage_resid2\n",
    "from sklearn.preprocessing import OrdinalEncoder, OneHotEncoder\n",
    "from statsmodels.stats.outliers_influence import OLSInfluence\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "from statsmodels.formula.api import glm\n",
    "from statsmodels.formula.api import ols\n",
    "\n",
    "import statsmodels.formula.api as smf\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "import pandas as pd  \n",
    "import numpy as np\n",
    "\n",
    "np.random.seed(42)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **1.** Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "    DATASET INFORMATIONS\n",
    "\n",
    "    |--------------|------------|-------|--------------------------------------------------------------|\n",
    "    | Name         | Data Type  | Meas. | Description                                                  |\n",
    "    |--------------|------------|-------|--------------------------------------------------------------|\n",
    "    | Gender       | nominal    |       | 1 represents Male , 2 represents Female                      |\n",
    "    | Age Range    | nominal    |       | 1 represents >18 years of age, 2 represents <18 years of age |\n",
    "    | Head Size    | continuous | cm^3  | Head volume                                                  |\n",
    "    | Brain Weight | continuous | grams | Mass of brains                                               |\n",
    "    |--------------|------------|-------|--------------------------------------------------------------|  \n",
    "\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Gender</th>\n",
       "      <th>Age Range</th>\n",
       "      <th>Head Size(cm^3)</th>\n",
       "      <th>Brain Weight(grams)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4512</td>\n",
       "      <td>1530</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3738</td>\n",
       "      <td>1297</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4261</td>\n",
       "      <td>1335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3777</td>\n",
       "      <td>1282</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4177</td>\n",
       "      <td>1590</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Gender  Age Range  Head Size(cm^3)  Brain Weight(grams)\n",
       "0       1          1             4512                 1530\n",
       "1       1          1             3738                 1297\n",
       "2       1          1             4261                 1335\n",
       "3       1          1             3777                 1282\n",
       "4       1          1             4177                 1590"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Import the dataset\n",
    "data = pd.read_csv('./dataset.csv')\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['Gender', 'Age Range', 'Head Size(cm^3)', 'Brain Weight(grams)'], dtype='object')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['Gender', 'Age Range', 'Head Size', 'Brain Weight'], dtype='object')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.columns = ['Gender', 'Age Range', 'Head Size', 'Brain Weight']\n",
    "data.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 237 entries, 0 to 236\n",
      "Data columns (total 4 columns):\n",
      " #   Column        Non-Null Count  Dtype\n",
      "---  ------        --------------  -----\n",
      " 0   Gender        237 non-null    int64\n",
      " 1   Age Range     237 non-null    int64\n",
      " 2   Head Size     237 non-null    int64\n",
      " 3   Brain Weight  237 non-null    int64\n",
      "dtypes: int64(4)\n",
      "memory usage: 7.5 KB\n"
     ]
    }
   ],
   "source": [
    "data.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **2.** Data Preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Gender          0\n",
       "Age Range       0\n",
       "Head Size       0\n",
       "Brain Weight    0\n",
       "dtype: int64"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Null elements\n",
    "data.isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are no missing values in this dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **3.** Data Preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = data[['Gender', 'Age Range', 'Head Size']]\n",
    "Y = data['Brain Weight']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Gender       int64\n",
       "Age Range    int64\n",
       "Head Size    int64\n",
       "dtype: object"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is not good, the 'Age Range' and the 'Gender' are not numbers that should be considered as numbers, they are categories and should therefore be dummies-converted. Actually to 'Age Range' we can apply an Ordinal-Encoding (since it makes sense to give growth meaning), while to 'Gender' we apply a One-Hot-Encoder."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **3.1.** One-Hot-Encoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "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>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   0  1\n",
       "0  1  0\n",
       "1  1  0\n",
       "2  1  0\n",
       "3  1  0\n",
       "4  1  0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## ONE-HOT-ENCODING: 'Gender'\n",
    "\n",
    "# Encoding\n",
    "gender_encoded = OneHotEncoder().fit_transform(X[['Gender']]);\n",
    "\n",
    "# Transform in a DataFrame\n",
    "gender_encoded = pd.DataFrame(gender_encoded.toarray()).astype(int)\n",
    "\n",
    "# Show results\n",
    "gender_encoded.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    1\n",
       "1    1\n",
       "2    1\n",
       "3    1\n",
       "4    1\n",
       "Name: Gender, dtype: int64"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# What is the association?\n",
    "X['Gender'].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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>Male</th>\n",
       "      <th>Female</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>232</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>233</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>234</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>235</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>236</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>237 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Male  Female\n",
       "0       1       0\n",
       "1       1       0\n",
       "2       1       0\n",
       "3       1       0\n",
       "4       1       0\n",
       "..    ...     ...\n",
       "232     0       1\n",
       "233     0       1\n",
       "234     0       1\n",
       "235     0       1\n",
       "236     0       1\n",
       "\n",
       "[237 rows x 2 columns]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# It is possible to conclude:\n",
    "#   1 -> Male   -> [1, 0]\n",
    "#   2 -> Female -> [0, 1]\n",
    "\n",
    "gender_encoded.columns = ['Male', 'Female']\n",
    "gender_encoded"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **3.2.** Ordinal-Encoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "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>0</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   0\n",
       "0  0\n",
       "1  0\n",
       "2  0\n",
       "3  0\n",
       "4  0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## ORDINAL-ENCODING: 'Age Range' \n",
    "\n",
    "# Encoding\n",
    "age_encoded = OrdinalEncoder().fit_transform(X[['Age Range']]);\n",
    "\n",
    "# Transform in a DataFrame\n",
    "age_encoded = pd.DataFrame(age_encoded).astype(int)\n",
    "\n",
    "# Show results\n",
    "age_encoded.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "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>0</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   0\n",
       "0  1\n",
       "1  1\n",
       "2  1\n",
       "3  1\n",
       "4  1"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Perform a replacement to have: 0 == <18, 1 == >18\n",
    "age_encoded = age_encoded.replace({0:1, 1:0})\n",
    "age_encoded.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    1\n",
       "1    1\n",
       "2    1\n",
       "3    1\n",
       "4    1\n",
       "Name: Age Range, dtype: int64"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# What is the association?\n",
    "X['Age Range'].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It is possible to conclude:\n",
    "#   1 -> >18 -> 1\n",
    "#   2 -> <18 -> 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **3.3.** Update"
   ]
  },
  {
   "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>Age Range</th>\n",
       "      <th>Head Size</th>\n",
       "      <th>Male</th>\n",
       "      <th>Female</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>4512</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>3738</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>4261</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>3777</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>4177</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>232</th>\n",
       "      <td>0</td>\n",
       "      <td>3214</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>233</th>\n",
       "      <td>0</td>\n",
       "      <td>3394</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>234</th>\n",
       "      <td>0</td>\n",
       "      <td>3233</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>235</th>\n",
       "      <td>0</td>\n",
       "      <td>3352</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>236</th>\n",
       "      <td>0</td>\n",
       "      <td>3391</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>237 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Age Range  Head Size  Male  Female\n",
       "0            1       4512     1       0\n",
       "1            1       3738     1       0\n",
       "2            1       4261     1       0\n",
       "3            1       3777     1       0\n",
       "4            1       4177     1       0\n",
       "..         ...        ...   ...     ...\n",
       "232          0       3214     0       1\n",
       "233          0       3394     0       1\n",
       "234          0       3233     0       1\n",
       "235          0       3352     0       1\n",
       "236          0       3391     0       1\n",
       "\n",
       "[237 rows x 4 columns]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Add 'Gender' encoded\n",
    "X = pd.concat([X, gender_encoded], axis=1)\n",
    "\n",
    "# Drop the old 'Gender'  \n",
    "X = X.drop('Gender', axis=1) \n",
    "\n",
    "# Substitute 'Age Range' with 'age_encoded'\n",
    "X['Age Range'] = age_encoded\n",
    "\n",
    "# Check\n",
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **3.4.** Split in train and validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(189, 4)\n",
      "(48, 4)\n",
      "(189,)\n",
      "(48,)\n"
     ]
    }
   ],
   "source": [
    "X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.2, random_state=5)\n",
    "print(X_train.shape)\n",
    "print(X_test.shape)\n",
    "print(Y_train.shape)\n",
    "print(Y_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **4.** Linear Regression with **Statsmodels**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **4.1.** Basic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "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>const</th>\n",
       "      <th>Age Range</th>\n",
       "      <th>Head Size</th>\n",
       "      <th>Male</th>\n",
       "      <th>Female</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>4512</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>3738</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>4261</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>3777</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>4177</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   const  Age Range  Head Size  Male  Female\n",
       "0    1.0          1       4512     1       0\n",
       "1    1.0          1       3738     1       0\n",
       "2    1.0          1       4261     1       0\n",
       "3    1.0          1       3777     1       0\n",
       "4    1.0          1       4177     1       0"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# In statsmodel we have to manually add the intercept\n",
    "X_with_intercept = sm.add_constant(X)\n",
    "X_with_intercept.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Divide in train and test (the new dataset, with intercept)\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(X_with_intercept, Y, test_size = 0.2, random_state=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit regression model\n",
    "model = sm.OLS(Y, X_with_intercept)\n",
    "model = model.fit();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:           Brain Weight   R-squared:                       0.653\n",
      "Model:                            OLS   Adj. R-squared:                  0.648\n",
      "Method:                 Least Squares   F-statistic:                     146.0\n",
      "Date:                Thu, 16 Jun 2022   Prob (F-statistic):           2.94e-53\n",
      "Time:                        18:54:45   Log-Likelihood:                -1345.7\n",
      "No. Observations:                 237   AIC:                             2699.\n",
      "Df Residuals:                     233   BIC:                             2713.\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const        255.2074     35.975      7.094      0.000     184.329     326.085\n",
      "Age Range     23.9684      9.481      2.528      0.012       5.290      42.647\n",
      "Head Size      0.2442      0.015     16.212      0.000       0.215       0.274\n",
      "Male         138.8753     21.339      6.508      0.000      96.833     180.917\n",
      "Female       116.3321     15.902      7.315      0.000      85.001     147.663\n",
      "==============================================================================\n",
      "Omnibus:                        7.989   Durbin-Watson:                   1.922\n",
      "Prob(Omnibus):                  0.018   Jarque-Bera (JB):                8.255\n",
      "Skew:                           0.357   Prob(JB):                       0.0161\n",
      "Kurtosis:                       3.571   Cond. No.                     1.29e+19\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The smallest eigenvalue is 1.9e-29. This might indicate that there are\n",
      "strong multicollinearity problems or that the design matrix is singular.\n"
     ]
    }
   ],
   "source": [
    "# Inspect the results\n",
    "print(model.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:\n",
      "--------------------------\n",
      "const        255.207360\n",
      "Age Range     23.968445\n",
      "Head Size      0.244212\n",
      "Male         138.875307\n",
      "Female       116.332053\n",
      "dtype: float64\n",
      "\n",
      "R2:\n",
      "--------------------------\n",
      "0.6528346387526656\n"
     ]
    }
   ],
   "source": [
    "print('Parameters:')\n",
    "print('--------------------------')\n",
    "print(model.params)\n",
    "print('\\nR2:')\n",
    "print('--------------------------')\n",
    "print(model.rsquared)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R2 Score: 0.6119769433954168\n",
      "Mean Absolute Error: 57.41370465780287\n",
      "Mean Square Error: 5433.661065204953\n",
      "Root Mean Square Error: 73.71337100692759\n"
     ]
    }
   ],
   "source": [
    "# Evaluation\n",
    "Y_pred_train = model.predict(X_train)\n",
    "\n",
    "# Metrics\n",
    "r2   = r2_score(Y_train, Y_pred_train)\n",
    "mae  = mean_absolute_error(Y_train, Y_pred_train)\n",
    "mse  = mean_squared_error(Y_train, Y_pred_train)\n",
    "rmse = np.sqrt(mse)\n",
    "\n",
    "print(f'R2 Score: {r2}')\n",
    "print(f'Mean Absolute Error: {mae}')\n",
    "print(f'Mean Square Error: {mse}')\n",
    "print(f'Root Mean Square Error: {rmse}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **4.2.** Advanced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "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>const</th>\n",
       "      <th>age_range</th>\n",
       "      <th>head_size</th>\n",
       "      <th>male</th>\n",
       "      <th>female</th>\n",
       "      <th>brain_weight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>4158</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1432</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>174</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>3615</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>3458</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1225</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>198</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3493</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1240</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>217</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3170</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1178</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>118</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>4114</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1362</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>230</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3685</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>189</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3703</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>206</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3402</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1220</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3478</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1270</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>189 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     const  age_range  head_size  male  female  brain_weight\n",
       "34     1.0          1       4158     1       0          1432\n",
       "174    1.0          1       3615     0       1          1350\n",
       "55     1.0          1       3458     1       0          1225\n",
       "198    1.0          0       3493     0       1          1240\n",
       "217    1.0          0       3170     0       1          1178\n",
       "..     ...        ...        ...   ...     ...           ...\n",
       "118    1.0          0       4114     1       0          1362\n",
       "230    1.0          0       3685     0       1          1350\n",
       "189    1.0          0       3703     0       1          1350\n",
       "206    1.0          0       3402     0       1          1220\n",
       "99     1.0          0       3478     1       0          1270\n",
       "\n",
       "[189 rows x 6 columns]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_advanced = pd.concat([X_train, Y_train], axis=1, join=\"inner\")\n",
    "data_advanced.columns = ['const', 'age_range', 'head_size', 'male', 'female', 'brain_weight']\n",
    "data_advanced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:           brain_weight   R-squared:                       1.000\n",
      "Model:                            OLS   Adj. R-squared:                  1.000\n",
      "Method:                 Least Squares   F-statistic:                 6.184e+29\n",
      "Date:                Thu, 16 Jun 2022   Prob (F-statistic):               0.00\n",
      "Time:                        18:58:07   Log-Likelihood:                 4950.2\n",
      "No. Observations:                 189   AIC:                            -9890.\n",
      "Df Residuals:                     184   BIC:                            -9874.\n",
      "Df Model:                           4                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "================================================================================\n",
      "                   coef    std err          t      P>|t|      [0.025      0.975]\n",
      "--------------------------------------------------------------------------------\n",
      "Intercept    -2.224e-13   3.87e-13     -0.575      0.566   -9.86e-13    5.41e-13\n",
      "const         1.331e-13   3.87e-13      0.344      0.731    -6.3e-13    8.96e-13\n",
      "age_range    -1.298e-13   1.56e-13     -0.834      0.405   -4.37e-13    1.77e-13\n",
      "head_size     2.577e-16   3.49e-16      0.738      0.462   -4.32e-16    9.47e-16\n",
      "male         -7.184e-14   2.49e-13     -0.288      0.774   -5.64e-13     4.2e-13\n",
      "female       -9.374e-14   1.69e-13     -0.555      0.580   -4.27e-13     2.4e-13\n",
      "brain_weight     1.0000   1.02e-15   9.77e+14      0.000       1.000       1.000\n",
      "==============================================================================\n",
      "Omnibus:                       16.338   Durbin-Watson:                   0.033\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               11.131\n",
      "Skew:                          -0.467   Prob(JB):                      0.00383\n",
      "Kurtosis:                       2.264   Cond. No.                     6.55e+19\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The smallest eigenvalue is 6.57e-31. This might indicate that there are\n",
      "strong multicollinearity problems or that the design matrix is singular.\n"
     ]
    }
   ],
   "source": [
    "# Fit regression model\n",
    "results = smf.ols('brain_weight ~ const + age_range + head_size + male + female + brain_weight', data=data_advanced).fit()\n",
    "\n",
    "# Inspect the results\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.17081575,  0.28949481, -0.40392374, -1.76000874, -0.4298968 ,\n",
       "         1.29561558,  1.62414217],\n",
       "       [-1.80547454,  0.07133669,  1.79590552,  2.1731347 ,  0.1571729 ,\n",
       "        -2.57656716,  0.30772795],\n",
       "       [-0.74447704, -0.0994916 ,  0.74053131,  0.89607959, -0.53502283,\n",
       "        -0.66887857, -0.88822894],\n",
       "       [-1.60307423, -1.22349671,  0.1161628 ,  1.08563457, -1.09475496,\n",
       "        -2.17832207,  0.31930361],\n",
       "       [-1.89119075,  1.61537366, -3.64550613, -0.70950834,  1.67638495,\n",
       "        -3.76176773,  2.85450996]])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Influence test:\n",
    "#       once created, an object of class OLSInfluence holds attributes and methods \n",
    "#       that allow users to assess the influence of each observation. \n",
    "#       For example, we can compute and extract the first few rows of DFbetas by:\n",
    "\n",
    "test_class = OLSInfluence(results)\n",
    "test_class.dfbetas[:5, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAGDCAYAAAA2xlnwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABDbklEQVR4nO3deXxU9dn//9dFCBAgJixpKwmboCAoQg1oa11wQ+oCIlpRLCotP630bm9vqRtLRKxY/OrtXVxrraJWREHEgqIoKrWVioICIhYUlGBZFAQkERKu3x/nJE6GBGYgMxOG9/PxyIM5n/OZM9eZCbnmfM5nMXdHRERE0ke9VAcgIiIitUvJXUREJM0ouYuIiKQZJXcREZE0o+QuIiKSZpTcRURE0oySu4jslZm5mXUMHz9gZqNq+fiXm9nfa/OYNbzOHmOPPM/9fJ0iM3tif49zIEjWZyfxqZ/qAOTAZ2argF+4+5xUx3IgMrN2wKfAi+7+04jyJ4AV7l6UotCq5e5XpTqGfXUgxy4SD125ywHPzDJSHUMtOc7Mfry/BzGztP7Snu7nl0oWUF5IA/oQJWHMrJ6Z3WBmK83sSzObYmbNw30vmtnwqPrvm9mA8HFnM3vFzL4ys+VmdlFEvUfN7H4zm2Vm3wC9zexsM1toZlvM7HMzK4o69s/NbHUYxygzW2Vmp+8tzmrOaZmZnROxXd/MNpjZD82skZk9ER5js5m9Y2bfj+Mt+wNw2x7ez1+a2YrwPZlhZq0i9rmZXWNm/wb+bWanmNkaM/udma03sy/MrL+Z/dTMPg6PcVPE83uZ2T/DuL8ws4lm1qCGOB41s3Hh4xfMbFvEzy4zuzzct6fPsEV4DlvM7F9Ahz2cd7vw/Iaa2WfAa2H5leHnscnMZptZ27DczOzu8Ly3mNliMzsqOvZwe0R4vmvN7Mqo133dzH4RsV2l+dnM7gl/17aY2btmdmIN8cf8e2Fm15tZsZltDd+z08LyrDD2TWb2YRj3mojnVbmdEPUZNTOzv4W/p5vCxwVR53mbmb0FbAcOq63PTlJHyV0S6ddAf+BkoBWwCbg33PcUMKiiopl1AdoCM82sCfAK8Ffge8DFwH1hnQqXECTCbODvwDfAz4Fc4GzgajPrH3Hs+4BLgUOBHCA/xjijVYkb6ANsdPf3gCHhsVsDLYCrgJKa357d3AccYeGXjkhmdipwO3BReA6rgclR1foDxwEV79MPgEYE5zoa+BMwGDgWOBEYZWbtw7rlwH8DLYEfAacBv9pbwO5+rrs3dfemwIXAf4BXY/gM7wVKw3O5MvzZm5OBI4E+ZtYPuAkYAOQB8wg+G4AzgZOAIwg+j4uAL6MPZmZnAdcBZwCHA7u973vxDtAdaE5wns+YWaNq6sX0e2FmnYDhQE93zyb43VoV7h5DkEQ7hOVD4oizHvAXgv9fbcLXnhhV5zJgGMH/pw3U/mcnyebu+tHPfv0Q/AE6vZryZcBpEduHAjsJ+npkEyTktuG+24BHwsc/A+ZFHetBYEz4+FFg0l5i+l/g7vDxaOCpiH2NgR0VMe8pzmqO2xHYCjQOt58ERoePrwT+AXSL8/1rB3j4vvwKeDssfwIoCh//GfhDxHOahjG2C7cdODVi/ykEf8Qzwu3ssM5xEXXeBfrXENNvgecith3oGPH+j4uqfwSwHvjJ3j5DICOMvXPEvt8Df9/L+3NYRNmLwNCI7XoEV51tgVOBj4HjgXpRx6qMHXgEGB91DpHn+TpBX5KK/ZfXFGO4fxNwTPi4CHgint+L8HdrPcGXjMyofZ8AZ0VsDwPWVPf51PQZRezrDmyK2H4dGBuxXWufnX5S96Mrd0mktsBzYVPkZoIkWg583923AjMJrgoguBp+MuJ5x1U8L3zupQRXohU+j3whMzvOzOaGTY9fE1wdtQx3t4qs7+7bqXolV2Oc0Sfk7ivC/eeaWWPgPIIrHIDHgdnA5LCZ9w9mlhnD+xTpYeD7ZnZuVHkrgqv1iji2hecQ2QLxedRzvnT38vBxxZXiuoj9JQRfEjCzI8Lm2v+Y2RaCP9gtiYGZ5QDPAyPdvaLZek+fYR7BF5nIeFezd5H12wL3RBz7K8CAfHd/jeDK9F5gvZk9ZGaHVHO8Kr8XMcZQycyuC28LfB3GkEP171lMvxfh79ZvCb4YrDezyfbdrZd9jtXMGpvZgxbcltoCvAnkWtW+KtHvbW1/dpJkSu6SSJ8Dfd09N+KnkbsXh/ufAgaZ2Y8Imo/nRjzvjajnNXX3qyOOHb2c4V+BGUBrd88BHiD4Yw/wBRB5jzGLoHk01jijVTTN9wM+DP8o4+473f0Wd+8C/Bg4h+BWQczcfQdwC3BrRPwAawn+6FacQ5PwHCJj3J8lHu8HPgIOd/dDCJq8bc9PCforELz3c939oYhde/oMNwBlBM3UFdrEEGPk+X0O/H9Rx89y938AuPv/ufuxBLcojgBGVHO8L/YSwzcErTwVKr9chvfXf0fQ5N/M3XOBr6nmPYvn98Ld/+ruPyH4rB24I8ZYt9cUK/A/QCeCVptDCG5ZEBVr9Htb25+dJJmSu9SWzLDjUMVPfYIEe5t919EpL7xXWmEWwR+xscDT7r4rLP8bwb3ny8wsM/zpaWZH7uH1s4Gv3L3UzHoR3JOv8CzBlfaPLegkVkTVP2x7izPaZIL7ulfz3VU7ZtbbzI4Or4i2EDRf7qr+EHv0OMGXnbMiyp4CrjCz7mbWkODKer67r9qH41cnmyDmbWbWmeDcYnEb0AT4TVR5jZ9h2JowDSgKryq7EN89ZAg+sxvNrCsErQdmdmH4uGfYkpNJkKBLqf5zmAJcbmZdwlaYMVH7FwEDwhg7AkMj9mUTJLkNQH0zGw1U1zoQ8++FmXUys1PDz7eUoGWlot6U8HybWdAZ7tfVxHqJmWWEfQlOjoq1BNhsQUfR6POMlujPTpJAyV1qyyyCPyAVP0XAPQRX0y+b2VbgbYIOXwC4+7cEfyhOJyJJhk32ZxI02a8l6KR1B9BwD6//K2Bs+DqjCf4YVhxvKcEfw8kEV0DbCO5tfhtW2WOc0dz9C+CfBFdhT0fs+gHBF4ktBE33bxAk6orJUx7YQ/yRxy8Pz6F5RNkcYBQwNTyHDnx3S6M2XEfwhWgrQce7p/dcvdIggnvbm+y7HvOXxvAZDie4JfAfgvvDf4knWHd/Ljze5LCpeQnQN9x9SHgOmwiajL8EJlRzjBcJ+ma8BqwI/410N0HfjHXAY3x32wiCZvaXCO7tryZIxtG3RSrU+HsRpSEwHthI8L58D7gx3HdL+DqfAi9X8/zfAOcCmwma0KdH7PtfICs87tth3DVK9GcnyWHu+9OSJ3LgMbOmBH8ED3f3T1McjkjczOwUgg57BXupKgcpXbnLQcHMzg2bEZsAdwKL+W6YkYhIWlFyl4NFP4ImxrUEY5ovdjVbiUiaUrO8iIhImtGVu4iISJpJaHI3s7MsmJd4hZndUM3+hmb2dLh/vgWrYxEOvXjMgjmhl5nZjbsdXERERKqVsNWVwjGd9xLM27wGeMfMZrj7hxHVhhJMg9jRzC4mGG7xM4I5qhu6+9Hh+NMPzeypPY3pbdmypbdr1y5BZyMiIlL3vPvuuxvdPS+6PJFLJ/YiWIv6EwAzm0w4o1dEnX4E46EhGAc60cyMYLakJuFEKFkEY0237OnF2rVrx4IFC2r1BEREROoyM6t2+t9ENsvnU3VShzVUnQe7Sh13LyOYvrEFQaL/hmCyjs+AO939qwTGKiIikjbqaoe6XgQLd7QC2gP/Y2aHRVcys2FmtsDMFmzYsCHZMYqIiNRJiUzuxVRdXKCAqotcVKkTNsHnEEwVeQnwUrjgwnrgLaAw+gXc/SF3L3T3wry83W451JqJEydSWFhIw4YNufzyy6vsmzJlCkceeSTZ2dl06dKF6dOnJywOERGRWCQyub8DHG5m7cPFOi4mmL870gy+W3RgIPBaOLHIZwRrMlesfnU8wYpVKdGqVStGjhzJlVdeWaW8uLiYwYMHc9ddd7FlyxYmTJjAJZdcwvr161MUqYiISAKTe3gPfTjBAgvLgCnuvtTMxprZeWG1PwMtzGwFcC1QMVzuXqCpmS0l+JLwF3f/IFGx7s2AAQPo378/LVq0qFK+Zs0acnNz6du3L2bG2WefTZMmTVi5cmWKIhUREUlsb3ncfRbBamGRZaMjHpcSDHuLft626srrmsLCQo488khmzJjB2WefzQsvvEDDhg3p1q1bqkMTEZGDWEKTe7rLyMjg5z//OZdccgmlpaU0aNCAZ555hiZNmqQ6NBEROYjV1d7yB4Q5c+bwu9/9jtdff50dO3bwxhtv8Itf/IJFixalOjQRETmIKbnvh0WLFnHSSSdRWFhIvXr16NmzJ8cddxxz5sxJdWgiInIQU3KPQVlZGaWlpZSXl1NeXk5paSllZWX07NmTefPmVV6pL1y4kHnz5umeu4iIpJSSewzGjRtHVlYW48eP54knniArK4tx48Zx8sknU1RUxMCBA8nOzuaCCy7gpptu4uOPP65xXPyrr75K586dady4Mb1792b16u9mDrzuuus4/PDDyc7OpnPnzkyaNCnJZyoiIukgbdZzLyws9Loyt/y0adOoV68es2fPpqSkhEcffRSAjRs30qFDBx5++GHOPfdcRo0axbx583j77bcBGDNmDIMGDeKII47gnXfe4ayzzmLmzJn8+Mc/TuHZiIhIXWVm77r7bpO86co9AWoaFz9t2jS6du3KhRdeSKNGjSgqKuL999/no4+C+XluueUWOnfuTL169TjuuOM48cQT+ec//5mKUxARkQOYknsSLV26lGOOOaZyu0mTJnTo0IGlS5fuVrekpIR33nmHrl27JjNEERFJA0ruSbRt2zZycnKqlOXk5LB169bd6l511VUcc8wx9OnTJ1nhiYhImtAkNknUtGlTtmypuiz9li1byM7OrlI2YsQIlixZwty5cwmWtxcREYmdkvs+mL6wmAmzl7N2cwmtcrMY0acT/XtEL1W/u65du/LYY49Vbn/zzTesXLmyStP7mDFjePHFF3njjTc45JBDEhK/iIikNzXLx2n6wmJunLaY4s0lOFC8uYQbpy1m+sLvVrOtaVz8+eefz5IlS5g6dSqlpaWMHTuWbt260blzZwBuv/12/vrXvzJnzpzdOuOJiIjESsk9ThNmL6dkZ3mVspKd5UyYvbxyu6Zx8Xl5eUydOpWbb76ZZs2aMX/+fCZPnlz5vJtuuonPPvuMjh070rRpU5o2bcrvf//7pJ2biIikB41zj1P7G2ZS3TtmwKfjz07464uIiFTQOPda0io3K65yERGRZFNyj9OIPp3IysyoUpaVmcGIPp1SFJGIiEhV6i0fp4pe8fvSW15ERCQZlNz3Qf8e+fuVzPd1KJ2IiEgslNyTrGIoXUWP+4qhdIASvIiI1Ardc0+yWIbSiYiI7A8l9yRbu7kkrnIREZF4KbknmYbSiYhIoim5J5mG0omISKKpQ12SaSidiIgkmpJ7CuzvUDoREZE9UbO8iIhImlFyFxERSTNK7iIiImlGyV1ERCTNKLmLiIikGSV3ERGRNKPkLiIikmaU3EVERNKMkruIiEiaUXIXERFJM0ruIiIiaUbJXUREJM0kNLmb2VlmttzMVpjZDdXsb2hmT4f755tZu7D8UjNbFPGzy8y6JzJWERGRdJGw5G5mGcC9QF+gCzDIzLpEVRsKbHL3jsDdwB0A7v6ku3d39+7AZcCn7r4oUbGKiIikk0ReufcCVrj7J+6+A5gM9Iuq0w94LHz8LHCamVlUnUHhc0VERCQGiUzu+cDnEdtrwrJq67h7GfA10CKqzs+Ap6p7ATMbZmYLzGzBhg0baiVoERGRA12d7lBnZscB2919SXX73f0hdy9098K8vLyEx/Pvf/+bRo0aMXjwYADmzp3L0UcfTW5uLi1atOD888+nuLg44XGIiIjsSSKTezHQOmK7ICyrto6Z1QdygC8j9l9MDVftqXDNNdfQs2fPyu0uXbowe/ZsNm/ezNq1azn88MO5+uqrd3vet99+y9ChQ2nbti3Z2dl0796dF198cbd6Y8eOxcyYM2dOQs9DRETSWyKT+zvA4WbW3swaECTqGVF1ZgBDwscDgdfc3QHMrB5wEXXkfvvkyZPJzc3ltNNOqyz7/ve/T6tWrSq3MzIyWLFixW7PLSsro3Xr1rzxxht8/fXXjBs3josuuohVq1ZV1lm5ciXPPPMMhx56aELPQ0RE0l/Cknt4D304MBtYBkxx96VmNtbMzgur/RloYWYrgGuByOFyJwGfu/sniYoxVlu2bGH06NHcddddu+377LPPyM3NJSsrizvvvJPf/e53u9Vp0qQJRUVFtGvXjnr16nHOOefQvn173n333co611xzDXfccQcNGjRI6LmIiEj6q5/Ig7v7LGBWVNnoiMelwIU1PPd14PhExherUaNGMXToUAoKCnbb16ZNGzZv3sxXX33Fn/70Jzp37rzX461bt46PP/6Yrl27AvDMM8/QsGFDfvrTn9Z67CIicvBJaHJPB4sWLWLOnDksXLhwj/WaN2/OkCFDOOaYYyguLqZ+/erf2p07d3LppZcyZMgQOnfuzNatW7npppt45ZVXEhG+iIgchJTc9+L1119n1apVtGnTBoBt27ZRXl7Ohx9+yHvvvVelbllZGevXr2fLli00b958t2Pt2rWLyy67jAYNGjBx4kQAioqKuOyyy2jXrl3Cz0VERA4OFvZfO+AVFhb6ggULav2427dvZ8uWLZXbd955J6tWreL+++9n3rx5dO3alcMPP5wvv/ySa665hhUrVuyW9AHcnSuvvJJVq1Yxa9YssrKyAOjevTtr1qypvNLfsGEDOTk5XH/99Vx//fW1fj4iIpI+zOxddy+MLteV+140btyYxo0bV243bdqURo0akZeXR3FxMf/zP//D+vXryc7O5pRTTuG5556r9jhXX301y5YtY86cOZWJHeDVV19l586dlds9e/bkrrvuom/fvok7KRERSWu6ck+C1atX065dOxo2bFjlXvyDDz7IpZdeWqVuu3btePjhhzn99NOTHaaIiBxgdOWeQm3btiXWL1GRY99FRET2RZ2eflZERETip+QuIiKSZtQsvw+mLyxmwuzlrN1cQqvcLEb06UT/HtEL3omIiKSGknucpi8s5sZpiynZWQ5A8eYSbpy2GCDmBK8vByIikkhqlo/ThNnLKxN7hZKd5UyYvTym51d8OSjeXILz3ZeD6Qu1VKyIiNQOJfc4rd1cEld5tP39ciAiIrI3Su5xapWbFVd5tP39ciAiIqk3efJkjjzySJo0aUKHDh2YN28eb7/9NmeccQbNmzcnLy+PCy+8kC+++CIl8Sm5x2lEn05kZWZUKcvKzGBEn04xPX9/vxyIiEhqvfLKK1x//fX85S9/YevWrbz55pscdthhbNq0iWHDhrFq1SpWr15NdnY2V1xxRUpiVIe6OFV0fNvXDnEj+nSq0iEP4vtyICIiqTVmzBhGjx7N8ccHq5Ln5+dX+bfC8OHDOfnkk5MeHyi575P+PfL3uXf7/n45EBGR1CkvL2fBggWcd955dOzYkdLSUvr378+ECROqrBsC8Oabb9K1a9eUxKm55UVERGK0du1a8vPzOfbYY3nhhRfIzMykX79+nHLKKdx2222V9T744ANOOeUUnn/+eU488cSExVPT3PK65y4iIhKjiqvzX//61xx66KG0bNmSa6+9llmzZlXWWbFiBX379uWee+5JaGLfEyV3ERGRGDVr1oyCggLMrLIs8vHq1as5/fTTGTVqFJdddlkqQgSU3EVERPZo4sSJFBYW0rBhQy6//HKuuOIK/vjHP7J+/XqmT5/O4MGDWbJkCT/+8Y858cQTGT58OFdddRXFxcX069eP5s2bU1BQwAMPPJC0mNWhTkREZA9atWrFyJEjmT17NiUlJYwaNYqNGzfSsWNHtm3bRp8+fXjqqac4++yz+fzzzykqKqKoqIiSkhLq1avH9u3b+fDDD+nduzedOnWid+/eCY9ZV+4iIiJ7MGDAAPr370+LFi0AyMzM5L777uPOO+/k+OOP58UXXyQ3N5eXX36ZRo0asWDBAv7zn/+wa9cu1q5dS2ZmJscccwwDBw7kkUceSUrMSu4iIiL7YOnSpRxzzDGV2xWz1S1dupSKkWiRI9LcnSVLliQlNiV3ERGRfbBt2zZycnKqlOXk5LB161ays7M54YQTuPXWWyktLeW9995j6tSpbN++PSmxKbmLiIjsg6ZNm7Jly5YqZVu2bCE7OxuAJ598kk8//ZTWrVtz9dVXM3jwYAoKCpISm5J7NaJ7Rkbavn07v/rVr2jZsiU5OTmcdNJJlfvmzp1L7969ycnJoV27dskNWkREkqpr1668//77ldvffPMNK1eurJyVrm3btvztb39jw4YNzJ8/n40bN9KrV6+kxKbe8tWI7hkZadiwYZSVlbFs2TKaN2/OokWLKvc1adKEK6+8kkGDBvH73/8+yVGLiEgilJWVUVZWRnl5OeXl5ZSWllK/fn3OP/98RowYwfV3PsSbJW34cObDZLRsx0cl2XQGli1bRkFBAQ0bNmTKlCm8/PLLLFu2LCkxK7lXY8CAAQAsWLCANWvWVJZ/9NFHzJgxgzVr1nDIIYcAcOyxx1bu79WrF7169WLOnDnJDVhERBJm3Lhx3HLLLZXbTzzxBGPGjKGoqIhr73iQ348cwc6v19Pg0CNoefZ13DhtMQCr3pjNbbfdxvbt2+nRowcvvfQSeXl5SYlZyT0O//rXv2jbti1jxozh8ccf59BDD6WoqIgLLrgg1aGJiEiCVIxbr87Lm7/Hob+oOjlNyc5yJsxezls3/Jbf/va3iQ+wGrrnHoc1a9awZMkScnJyWLt2LRMnTmTIkCFJa2YREZG6Ze3mkrjKk0XJPQ5ZWVlkZmYycuRIGjRowMknn0zv3r15+eWXUx2aiIikQKvcrLjKk0XJPQ7dunXbrSxywQARETm4jOjTiazMjCplWZkZjOjTKUURBZTcq1FWVkZpaWmVnpFlZWWcdNJJtGnThttvv52ysjLeeust5s6dS58+fQDYtWsXpaWl7Ny5E3entLSUHTt2pPhsREQkUfr3yOf2AUeTn5uFAfm5Wdw+4Gj698hPaVxK7tUYN24cWVlZjB8/nieeeIKsrCzGjRtHZmYmzz//PLNmzSInJ4df/vKXTJo0ic6dOwPw5ptvkpWVxU9/+lM+++wzsrKyaNy4Md27d+fFF18E4O233+aMM86gefPm5OXlceGFF/LFF19Uvnbfvn1p2rRp5U+DBg04+uijU/I+iIjIgcki5709kBUWFvqCBQtSHQYQTGQwYcIELr/8ctq0acOsWbMYNGgQixcvZtmyZZWrCNWvX5/hw4ezdu1aXnrppWqPdcopp3DqqacyevToJJ+FiIjszfSFxdw4bTElO8sry7IyM5J29W5m77p74W7lSu7J0a1bN8aMGbPbsLn33nuPk08+ma1bt+72nFWrVtGhQwdWrlypGe9EROqgE8a/RnE1PePzc7N464ZTE/76NSX3hDbLm9lZZrbczFaY2Q3V7G9oZk+H++ebWbuIfd3M7J9mttTMFptZo0TGmkjr1q3j448/rpySMNKbb75ZbTnApEmTOPHEE5XYRUTqqINuKJyZZQD3An2BLsAgM+sSVW0osMndOwJ3A3eEz60PPAFc5e5dgVOAnYmKNZF27tzJpZdeypAhQyrvzVf44IMPGDt2LBMmTKj2uZMmTdptbnsREak7DsahcL2AFe7+ibvvACYD/aLq9AMeCx8/C5xmwdiyM4EP3P19AHf/0t3LOcDs2rWLyy67jAYNGjBx4sQq+1asWEHfvn255557OPHEE3d77t///nf+85//MHDgwGSFKyIicaqrQ+ESOf1sPvB5xPYa4Lia6rh7mZl9DbQAjgDczGYDecBkd/9D9AuY2TBgGECbNm1q/QRqMn1hMRNmL2ft5hJa5WYxok+n3TpOuDtDhw5l3bp1zJo1i8zMzMp9q1ev5vTTT2fUqFFcdtll1b7GY489xoABA2jatGlCz0VERPZdxd/+veWEZKurc8vXB34C9AS2A6+GnQZejazk7g8BD0HQoS4ZgUX3jCzeXFK5SEDkh3n11VezbNky5syZQ1bWd80zxcXFnHrqqQwfPpyrrrqq2tcoKSlhypQpPPfccwk8ExERqQ39e+SnPJlHS2SzfDHQOmK7ICyrtk54nz0H+JLgKv9Nd9/o7tuBWcAPExhrzCbMXl5lyAN8t0hAhdWrV/Pggw+yaNEifvCDH1SOWX/yySd5+OGH+eSTTygqKqoynj3S9OnTyc3NpXfv3kk5JxERSS8JGwoXJuuPgdMIkvg7wCXuvjSizjXA0e5+lZldDAxw94vMrBnwKsHV+w7gJeBud59Z0+slayhc+xtmUt07ZsCn489O+OuLiIhUSPpQOHcvA4YDs4FlwBR3X2pmY83svLDan4EWZrYCuBa4IXzuJuAugi8Ei4D39pTYk6mu9owUERGpoEls4pTq2YhEREQq1HTlXlc71NVZdbVnpIiISAUl932wvz0jYxlKJyIisq+U3JMs1qF0IiIi+0pLviZZLEPpRERE9oeSe5LV1UUGREQkfSi5J5mG0omISKIpuSdZXV1kQERE0oc61CWZhtKJiEiiKbmnQF1cZEBERNKHmuVFRETSjJJ7NSZOnEhhYSENGzbk8ssvr7Jv+/bt/OpXv6Jly5bk5ORw0kknVe6bMGECRx11FNnZ2bRv354JEyYkOXIRERE1y1erVatWjBw5ktmzZ1NSUnWI2rBhwygrK2PZsmU0b96cRYsWVe5zdyZNmkS3bt1YuXIlZ555Jq1bt+biiy9O8hmIiMjBTAvH7MHIkSNZs2YNjz76KAAfffQRvXr1Ys2aNRxyyCF7ff5//dd/4e788Y9/rNW4REREIAVLvqajf/3rX7Rt25YxY8bQsmVLjj76aKZOnVptXXdn3rx5dO3aNclRiojIwU7JPQ5r1qxhyZIl5OTksHbtWiZOnMiQIUNYtmzZbnWLiorYtWsXV1xxRQoiFRGRg5mSexyysrLIzMxk5MiRNGjQgJNPPpnevXvz8ssvV6k3ceJEJk2axMyZM2nYsGGKohURkYOVOtTFoVu3bruVmVmV7UceeYTx48fz5ptvUlBQkKzQREREKunKvRplZWWUlpZSXl5OeXk5paWllJWVcdJJJ9GmTRtuv/12ysrKeOutt5g7dy59+vQB4Mknn+Smm27ilVde4bDDDkvxWYiIyMFKyb0a48aNIysri/Hjx/PEE0+QlZXFuHHjyMzM5Pnnn2fWrFnk5OTwy1/+kkmTJtG5c2cg6F3/5Zdf0rNnT5o2bUrTpk256qqrUnw2IiJysNFQOBERkQOUhsKJiIgcJJTcRURE0oySu4iISJrRULh9MH1hsdZjFxGROkvJPU7TFxZz47TFlOwsB6B4cwk3TlsMoAQvIiJ1gprl4zRh9vLKxF6hZGc5E2YvT1FEIiIiVSm5x2nt5pK4ykVERJJNyT1OrXKz4ioXETmYffvttwwdOpS2bduSnZ1N9+7defHFFyv3v/rqq3Tu3JnGjRvTu3dvVq9eXbmvuLiYfv360bx5cwoKCnjggQdScQoHJCX3OI3o04mszIwqZVmZGYzo0ylFEYmI1F1lZWW0bt2aN954g6+//ppx48Zx0UUXsWrVKjZu3MiAAQO49dZb+eqrrygsLORnP/tZ5XMHDx5M+/btWbduHTNnzuSmm25i7ty5KTybA4dmqNsH6i0vIrLvunXrxpgxY/jyyy959NFH+cc//gHAN998Q8uWLVm4cCEFBQVkZ2ezfv168vLyABg2bBglJSU8/vjjqQy/Tqlphjr1lt8H/XvkK5mLiOyDdevW8fHHH9O1a1fuv/9+jjnmmMp9TZo0oUOHDixdupT8/OBvbOQFqLuzZMmSpMd8IFKzfApMX1jMCeNfo/0NMzlh/GtMX1ic6pBERBJu586dXHrppQwZMoTOnTuzbds2cnJyqtTJyclh69atZGdnc8IJJ3DrrbdSWlrKe++9x9SpU9m+fXuKoj+wKLknWcU4+eLNJTjfjZNXgheRdLZr1y4uu+wyGjRowMSJEwFo2rQpW7ZsqVJvy5YtZGdnA8Ey2p9++imtW7fm6quvZvDgwRQUFCQ99gORknuSaZy8iBxs3J2hQ4eybt06pk6dSmZmJgBdu3bl/fffr6z3zTffsHLlSrp27QpA27Zt+dvf/saGDRuYP38+GzdupFevXik5hwONknuSaZy8iBxsrr76apYtW8YLL7xAVtZ3w4bPP/98lixZwtSpUyktLWXs2LF069aNzp07A7Bs2TK2bt3Kjh07eOKJJ3j55Ze59tprU3UaBxQl9yTTOHkROZisXr2aBx98kEWLFvGDH/yApk2b0rRpU5588kny8vKYOnUqN998M82aNWP+/PlMnjy58rmzZ8/msMMOo1mzZjzwwAO89NJLlT3nZc8SOhTOzM4C7gEygIfdfXzU/obAJOBY4EvgZ+6+yszaAcuAirbqt939qj29VjKHwu2P6LnpIRgnf/uAo9UDX0RE4pL0oXBmlgHcC5wBrAHeMbMZ7v5hRLWhwCZ372hmFwN3ABUzGKx09+6Jii9VKhK4xsmLiEiiJHKcey9ghbt/AmBmk4F+QGRy7wcUhY+fBSaamSUwpjpB4+RFRCSREnnPPR/4PGJ7TVhWbR13LwO+BlqE+9qb2UIze8PMTkxgnCIiImmlrs5Q9wXQxt2/NLNjgelm1tXdqwyINLNhwDCANm3apCBMERGRuieRyb0YaB2xXRCWVVdnjZnVB3KALz3o5fctgLu/a2YrgSOAKj3m3P0h4CEIOtQl4iRERCSxtF5H7Utks/w7wOFm1t7MGgAXAzOi6swAhoSPBwKvububWV7YIQ8zOww4HPgkgbGKiEgKaNbOxIgpuZtZYzMbZWZ/CrcPN7Nz9vSc8B76cGA2wbC2Ke6+1MzGmtl5YbU/Ay3MbAVwLXBDWH4S8IGZLSLoaHeVu38V57mJiEgdp1k7EyPWZvm/AO8CPwq3i4FngL/t6UnuPguYFVU2OuJxKXBhNc+bCkyNMTYRETlAadbOxIi1Wb6Du/8B2Ang7tuBtB+yJiIiiaVZOxMj1uS+w8yyAAcwsw6EHd5ERET21Yg+ncjKzKhSlpWZwYg+nVIUUXqItVl+DPAS0NrMngROAC5PVFAiInJw0KydiRHz3PJm1gI4nqA5/m1335jIwOJ1oMwtLyIiUlv2a255M/th+PCL8N82ZpYDrA57xYuIiEgdEWuz/H3AD4EPCK7cjwKWAjlmdrW7v5yg+ERERCROsXaoWwv0cPdCdz8W6EEwqcwZwB8SFZyIiIjEL9bkfoS7L63YCJdt7Vyx4puIiIjUHbE2yy81s/uByeH2z4APzawh4dh3ERERqRtiTe6XA78CfhtuvwVcR5DYe9d6VGlOiySIiEgixZTc3b0E+H/hT7RttRpRmqtYJKFiLuWKRRIAJXgREakVsS4cc7iZPWtmH5rZJxU/iQ4uHWmRBBERSbRYO9T9BbgfKCNohp8EPJGooNKZFkkQEZFEizW5Z7n7qwQz2q129yLg7MSFlb60SIKIiCRarMn9WzOrB/zbzIab2flA0wTGlbZ6d86Lq1xERCResSb33wCNgf8CjgUGA0MSFVQ6m/vRhrjKRURE4rXX3vJmlgH8zN2vI+gZf0XCo0pjuucuIiKJttcrd3cvB36ShFgOCrrnLiIiiRZrs/xCM5thZpeZ2YCKn4RGlqZ0z11E5MA3efJkjjzySJo0aUKHDh2YN28eAFOmTOHII48kOzubLl26MH369JTEF+sMdY2AL4FTI8ocmFbrEaU53XMXETmwvfLKK1x//fU8/fTT9OrViy++CFZDLy4uZvDgwTz//POcddZZzJo1iwsvvJBVq1bxve99L6kxxjpDne6zR9if6WN1z11E5MA2ZswYRo8ezfHHHw9Afn7w93/+/Pnk5ubSt29fAM4++2yaNGnCypUrk57cY52h7ggze9XMloTb3cxsZGJDq5sqpo8t3lyC8930sdMXFsf0fN1zFxE5cJWXl7NgwQI2bNhAx44dKSgoYPjw4ZSUlFBYWMiRRx7JjBkzKC8vZ/r06TRs2JBu3bolPc5Y77n/CbiRcAU4d/8AuDhRQdVl+zt97Ig+ncjKzKhSlpWZwYg+nWotRhERSYx169axc+dOnn32WebNm8eiRYtYuHAh48aNIyMjg5///OdccsklNGzYkEsuuYQHH3yQJk2aJD3OWJN7Y3f/V1RZWW0HcyDY32b1/j3yuX3A0eTnZmFAfm4Wtw84WovGiIgcALKyglbWX//61xx66KG0bNmSa6+9llmzZjFnzhx+97vf8frrr7Njxw7eeOMNfvGLX7Bo0aKkxxlrct9oZh0IOtFhZgOBLxIWVR2mZnURkYNXs2bNKCgowMwqyyoeL1q0iJNOOonCwkLq1atHz549Oe6445gzZ07S44w1uV8DPAh0NrNignXdr0pUUHXZ/jar7+89exERSa0rrriCP/7xj6xfv55NmzZx9913c84559CzZ8/KpnqAhQsXMm/evJTcc491KNxqdz/dzJoA9dx9ayKDqssqms/3tbf8nu7Zq2leRKTuGzVqFBs3buSII46gUaNGXHTRRdx88800atSIoqIiBg4cyLp168jLy+Omm27izDPPTHqM5u57r2T2GfAS8DTwmsfypCQrLCz0BQsWpDqMvWp/w0yqe/MM+HS8FtoTEZHYmdm77l4YXR5rs3xnYA5B8/ynZjbRzDQl7T7QPXsREUm0mJK7u2939ynuPgDoARwCvJHQyNLUiD6ddnvT64XlIiIitSHWK3fM7GQzuw94l2A62osSFlUaW7D6K3ZFle0Ky0VE5MBQ09zyDz/8MB07dqRp06acddZZrF27NiXxxTpD3SqCHvLzgKPd/SJ3n5rAuNLWU/M/j6tcRETqloq55f/yl7+wdetWbrrvaf571lp+cMntXP3bEfzX+If46quvaN++PYMGDUpJjLH2lu/m7lsSGslBoryGvog1lYuISN0SObf89IXF3Pn3jZR4U0pWvEOjI07g/vd30KbDBkaNGkV+fj4rV66kQ4cOSY0x1mb5H2hu+dqRETHxQSzlIiJSd0TPLf+zU7pTPGsiu3Z+G9bwyuHNFQPLlixZkvQ4Nbd8kg06rnVc5SIiUndEzy3//SH3sGPdJ3z9z6dpdNixbP/o7+xY/ylrNmxm7NixmBnbt29PepyaWz7JCts2J/oa3cJyERGp26Lnlm/d6gdk9+xPycoFZLXrTu5PLmHD9N/zxUNDadeuHdnZ2RQUFCQ9Ts0tn2S3vLB0t0lsPCwXEZG6LXpu+RF9OtGw/ndTkmf/8Bw6XvMIU95cwgUXXEBZWRlHHXVU0uNM6NzyZnaWmS03sxVmdkM1+xua2dPh/vlm1i5qfxsz22Zm18UYZ523afvOuMpFRKRuiZxb/uR2jTlkxct8v+uP8bId5Jb+h9+ffxQ/bFHOsGHD+M1vfkOzZs2SHmOsk9h84u6nA3lAZ3f/CXD+np5jZhnAvUBfoAswyMy6RFUbCmxy947A3cAdUfvvAl6MJUYREZFkGDVqFD179uSII47gyCOP5PQTj+Pfs/7E+yN7s+u1exh8Yid69erFj370I2699daUxBjrUDgA3P2biM1rgf/dQ/VewAp3/wTAzCYD/YAPI+r0A4rCx88CE83M3N3NrD/wKRD5miIiIimVmZnJfffdx3333VelvFGjRnzwwQcpiqqqmGeoq8bexm7lA5Ezs6wJy6qt4+5lwNdACzNrClwP3LLHAMyGmdkCM1uwYcOGeGJPmdyszLjKRURE4rU/yT2Rs64UAXe7+7Y9BuD+kLsXunthXl5eAsOpPUXnda12bvmi87qmIhwREUlDe2yWN7OtVJ/EDdjbMmbFQOTg7YKwrLo6a8ysPpADfAkcBww0sz8AucAuMyt194l7ec0DQkaGsavcq2yLiIjUlj0md3fP3o9jvwMcbmbtCZL4xcAlUXVmAEOAfwID+W6t+BMrKphZEbAtXRL7hNnL2Vle9fvSznJnwuzl9O8RfddCRETquukLi5kwezlrN5fQKjeLEX06pfzveVwd6uLh7mVmNhyYDWQAj7j7UjMbCyxw9xnAn4HHzWwF8BUHwax3azeXxFUuIiJ11/SFxdw4bTElO8sBKN5cwo3TFgOkNMEnLLkDuPssYFZU2eiIx6XAhXs5RlFCgkuRVrlZFFeTyFvl7u0uh4iI1DUTZi+vTOwVKuaWT2Vy358OdbIPRvTpRFZmRpWyrMwMRvTplKKIRERkX9XV1lgl9yTr3yOfC47Nr1wFLsOMC47NT/n9GRERiV9Nra6pbo1Vct8H0xcWc8L412h/w0xOGP8a0xdGDwLY83OnvltcuX57uTtT3y2O6xgiIlI31NXWWCX3OFV0nijeXILzXeeJWJPznu7PiIjIgaV/j3xuH3A0+blZGJCfm8XtA45OeWtsQjvUpaP97TxRV+/PiIjIvunfo+7dWtWVe5yq6+m+p/JodfX+jIiIpA8l9zhVdISLtTxaXb0/IyIi6UPN8nGq6AgXa3m0iqabujabkYiIpA8l9zhlmFWbyGO9coe6eX9GRETSh5rl47S/V+4iIiKJpiv3ONXGlXtdXGRARETSh5J7nPb3yr2uLjIgIiLpQ83ycdrf3vKaxEZERBJNyT1O+3vlrklsREQk0ZTc45Rfw2QzNZVH0yQ2IiKSaErucerdOS+u8miaxEZERBJNyT1OMz/4Iq7yaHV1kQEREUkf6i0fp03bd8ZVXh1NYiMiIomk5J4CGucuIiKJpOQep6zMepTs3FVteSymLyxmxDPvs3NX0Lu+eHMJI555H9A4dxERqR265x6nRlGd4fZWHq1oxtLKxF5h5y6naMbS/Y5NREQElNzjtrmGe+s1le9Wr6SG59dQLiIiEi8l9zhpnLqIiNR1Su5x2t9x6s0aZ8ZVLiIiEi8l9zj175HPBcfmV84ln2HGBcfGPrRtzLldycyoOg99ZoYx5tyutR6riIgcnJTc4zR9YTFT3y2unEu+3J2p7xYzfWFxTM/v3yOfCQOPqTKJzYSBx6invIiI1BoNhYvTnlZ1izVBaxIbERFJJF25x0mruomISF2n5B4n9ZYXEZG6Tsk9Tvu7KpyIiEiiKbnH6W/vV7/6W03lIiIiyabkHifNMCciInWdesungFaFExGRRFJyT7LpC4u5cdriyuF0xZtLuHHaYkCrwomISO1Qs3yS7WmcvIiISG1Qco+TxVkeTePkRUQk0RKa3M3sLDNbbmYrzOyGavY3NLOnw/3zzaxdWN7LzBaFP++b2fmJjDMejRtUv257TeXRNE5eREQSLWHJ3cwygHuBvkAXYJCZdYmqNhTY5O4dgbuBO8LyJUChu3cHzgIeNLM60T9g+47yuMqjaZy8iIgkWiKv3HsBK9z9E3ffAUwG+kXV6Qc8Fj5+FjjNzMzdt7t7WVjeCPAExhmX3BqWZq2pPNrcjzbEVS4iIhKvRCb3fODziO01YVm1dcJk/jXQAsDMjjOzpcBi4KqIZF/JzIaZ2QIzW7BhQ3KSo9fwNaOm8mjFNdxbr6lcREQkXnW2Q527z3f3rkBP4EYza1RNnYfcvdDdC/PyktOsvb+T2FSsAx9ruYiISLwSmdyLgdYR2wVhWbV1wnvqOcCXkRXcfRmwDTgqYZHGoV4NObim8mjlNVzi11QuIiISr0Qm93eAw82svZk1AC4GZkTVmQEMCR8PBF5zdw+fUx/AzNoCnYFVCYw1ZrtqyME1lUfLr6FXfE3lIiIi8UpYcg/vkQ8HZgPLgCnuvtTMxprZeWG1PwMtzGwFcC1QMVzuJ8D7ZrYIeA74lbtvTFSsyaTe8iIikmgJHV7m7rOAWVFloyMelwIXVvO8x4HHExnbvsrNyqz2/npuVmy95fe0qty4/kfvV2wiIiJQhzvU1VVF53WNqzyaVpUTEZFEU3KP04LVX8VVLiIikmxK7nF68u3P4ioXERFJNiX3ONXUKV4D2UREpK5QchcREUkzSu5xalLD6m81lUc7oUPzuMpFRETipeQep9vOP5qMqOnoMuoZt50f2zC2J3/5o90S+QkdmvPkL39UazGKiMjBrU4so3og6d8jn2cWfMZbK7/rHX98+2b07xG9Jk7NlMhFRCSRdOUep5HTF1dJ7ABvrfyKkdMXx3yM6QuLOWH8a7S/YSYnjH+N6Qujp9wXERHZd0rucXpq/udxlUebvrCYG6ctpnhzCU6w1OuN0xYrwYuISK1Rco/T/q7qNmH2ckp2llcpK9lZzoTZy/c7NhEREVByT7rizSVxlYuIiMRLyT3JMqz6hd9rKhcREYmXknuS7W+zvoiIyN4ouSeZrtxFRCTRlNyTTFfuIiKSaEruSZafmxVXuYiISLyU3JOsd+e8uMpFRETipeSeZDM/+CKuchERkXgpuSfZpu074yoXERGJl5K7iIhImlFyFxERSTNK7nGq6Q2L9Y2sV8Nw9prKRURE4qXkHqddcZbvVq+G4ew1lYuIiMRLyV1ERCTNKLmLiIikGSV3ERGRNKPkLiIikmaU3EVERNKMkruIiEiaUXIXERFJM0ruIiIiaUbJXUREJM0ouYuIiKQZJXcREZE0o+SeZFbDAjE1lYuIiMQrocndzM4ys+VmtsLMbqhmf0MzezrcP9/M2oXlZ5jZu2a2OPz31ETGmUyZNSz/VlO5iIhIvBKW3M0sA7gX6At0AQaZWZeoakOBTe7eEbgbuCMs3wic6+5HA0OAxxMVZ7LtKK9++beaykVEROKVyCv3XsAKd//E3XcAk4F+UXX6AY+Fj58FTjMzc/eF7r42LF8KZJlZwwTGKiIikjYSmdzzgc8jtteEZdXWcfcy4GugRVSdC4D33P3bBMUpIiKSVuqnOoA9MbOuBE31Z9awfxgwDKBNmzZJjExERKTuSuSVezHQOmK7ICyrto6Z1QdygC/D7QLgOeDn7r6yuhdw94fcvdDdC/Py8mo5/MTIzcqMq1xERCReiUzu7wCHm1l7M2sAXAzMiKozg6DDHMBA4DV3dzPLBWYCN7j7WwmMMenOOebQuMpFRETilbDkHt5DHw7MBpYBU9x9qZmNNbPzwmp/BlqY2QrgWqBiuNxwoCMw2swWhT/fS1SsyTT3ow1xlYuIiMQroffc3X0WMCuqbHTE41LgwmqeNw4Yl8jYUqV4c0lc5SIiIvHSDHUiIiJpRsldREQkzSi5i4iIpBkldxERkTSj5C4iIpJmlNxFRETSjJK7iIhImlFyFxERSTNK7iIiImlGyV1ERCTNKLnHKSuz+respnIREZFkU0aK0+0Duu32ptULy2PRIMPiKhcREYlXQheOSUf9e+QDMGH2ctZuLqFVbhYj+nSqLN+bHeUeV7mIiEi8lNz3Qf8e+TEncxERkWRTs7yIiEiaUXIXERFJM0ruIiIiaUbJXUREJM0ouYuIiKQZJXcREZE0o+QuIiKSZpTcRURE0oySu4iISJpRchcREUkzSu4iIiJpRsldREQkzSi5J1luVmZc5SIiIvFSck+yovO6klmv6trtmfWMovO6pigiERFJN1ryNcn2dz14ERGRvVFyTwGtBy8iIomkZnkREZE0o+QuIiKSZpTcRURE0oySu4iISJpRchcREUkzSu4iIiJpRsldREQkzSQ0uZvZWWa23MxWmNkN1exvaGZPh/vnm1m7sLyFmc01s21mNjGRMYqIiKSbhCV3M8sA7gX6Al2AQWbWJaraUGCTu3cE7gbuCMtLgVHAdYmKT0REJF0l8sq9F7DC3T9x9x3AZKBfVJ1+wGPh42eB08zM3P0bd/87QZIXERGROCQyuecDn0dsrwnLqq3j7mXA10CLWF/AzIaZ2QIzW7Bhw4b9DFdERCQ9HNAd6tz9IXcvdPfCvLy8VIcjIiJSJyQyuRcDrSO2C8KyauuYWX0gB/gygTGJiIikvUQm93eAw82svZk1AC4GZkTVmQEMCR8PBF5zd09gTCIiImkvYUu+unuZmQ0HZgMZwCPuvtTMxgIL3H0G8GfgcTNbAXxF8AUAADNbBRwCNDCz/sCZ7v5houIVERFJFwldz93dZwGzospGRzwuBS6s4bntEhmbiIhIurJ0aQU3sw3A6lTHsQ9aAhtTHUQCpON5peM5QXqeVzqeE+i8DiTJOqe27r5bj/K0Se4HKjNb4O6FqY6jtqXjeaXjOUF6nlc6nhPovA4kqT6nA3oonIiIiOxOyV1ERCTNKLmn3kOpDiBB0vG80vGcID3PKx3PCXReB5KUnpPuuYuIiKQZXbmLiIikGSX3FDGzR8xsvZktSXUstcXMWpvZXDP70MyWmtlvUh1TbTCzRmb2LzN7PzyvW1IdU20xswwzW2hmf0t1LLXFzFaZ2WIzW2RmC1IdT20xs1wze9bMPjKzZWb2o1THtD/MrFP4GVX8bDGz36Y6rtpgZv8d/q1YYmZPmVmjpMegZvnUMLOTgG3AJHc/KtXx1AYzOxQ41N3fM7Ns4F2g/4E+s6CZGdDE3beZWSbwd+A37v52ikPbb2Z2LVAIHOLu56Q6ntoQzm5Z6O5pNW7azB4D5rn7w+GU3o3dfXOKw6oVZpZBsNbIce5+IM5XUsnM8gn+RnRx9xIzmwLMcvdHkxmHrtxTxN3fJJhyN224+xfu/l74eCuwjN2X+T3geGBbuJkZ/hzw34rNrAA4G3g41bHInplZDnASwZTduPuOdEnsodOAlQd6Yo9QH8gKF0RrDKxNdgBK7pIQZtYO6AHMT3EotSJsvl4ErAdecfd0OK//BX4H7EpxHLXNgZfN7F0zG5bqYGpJe2AD8JfwNsrDZtYk1UHVoouBp1IdRG1w92LgTuAz4Avga3d/OdlxKLlLrTOzpsBU4LfuviXV8dQGdy939+4ESxf3MrMD+laKmZ0DrHf3d1MdSwL8xN1/CPQFrglvgR3o6gM/BO539x7AN8ANqQ2pdoS3GM4Dnkl1LLXBzJoB/Qi+kLUCmpjZ4GTHoeQutSq8Jz0VeNLdp6U6ntoWNoXOBc5KcSj76wTgvPD+9GTgVDN7IrUh1Y7wygl3Xw88B/RKbUS1Yg2wJqLF6FmCZJ8O+gLvufu6VAdSS04HPnX3De6+E5gG/DjZQSi5S60JO579GVjm7nelOp7aYmZ5ZpYbPs4CzgA+SmlQ+8ndb3T3gnD1xYuB19w96VcXtc3MmoSdOQmbrc8EDvgRKe7+H+BzM+sUFp0GHNAdVSMMIk2a5EOfAcebWePwb+JpBP2PkkrJPUXM7Cngn0AnM1tjZkNTHVMtOAG4jOAqsGJ4y09THVQtOBSYa2YfAO8Q3HNPm6Fjaeb7wN/N7H3gX8BMd38pxTHVll8DT4a/h92B36c2nP0XfgE7g+DqNi2ErSvPAu8BiwnybNJnq9NQOBERkTSjK3cREZE0o+QuIiKSZpTcRURE0oySu4iISJpRchcREUkzSu4iCWBmbmb/L2L7OjMrSnIMr5tZYfh4VsVY/f043imJWD3OzP5RQ/mjZjZwH4+5x1jD8cdUfCYWqqasu5n9M1zh6wMz+9m+xCOSbPVTHYBImvoWGGBmt+/L6mRmVt/dy2orGHdPynwD+xK3uyd99i7gt2a2hWBq0NuAN4Cu1ZStAn7u7v82s1bAu2Y2O80WbZE0pCt3kcQoI5i44r+jd5hZOzN7LbwSfNXM2oTlj5rZA2Y2H/hDuH2/mb1tZp+EV6OPhGt5PxpxvPvNbIHtYa15C9Y4b2lmV0VMMPSpmc0N958ZXqG+Z2bPhOsDYGZnWbB++HvAgBqOfbmZzTCz14BXw1niHjGzf4WLnPQL63UNyxaF5354WL4t/NfMbKKZLTezOcD3ouMPHxea2evh415h3AvN7B8RM7hFxndyxDkvNLNsd78byAP+C3jJ3V+uoexjd/83gLuvJVg4KK+690GkLlFyF0mce4FLLViuM9IfgcfcvRvwJPB/EfsKgB+7+7XhdjPgRwRfEmYAdwNdgaPNrHtY52Z3LwS6ASebWbeaAnL3B8IFcHoSzFd+V5g0RwKnhwuuLACuNbNGwJ+Ac4FjgR/s4Vx/CAx095OBmwmms+0F9AYmhDORXQXcE75+Yfj6kc4HOgFdgJ8T23zcHwEnhoupjKb6WduuA64JX/dEoMTMfkOwytr/AWeZ2RnVlUUexMx6AQ2AlTHEJZJSapYXSRB332JmkwiuBEsidv2I766CHwf+ELHvGXcvj9h+wd3dzBYD69x9MYCZLQXaAYuAiyxY2rQ+wVS5XYAP9hLePQQJ+AULVojrArwV3nZuQDA1cmeCBTD+Hb7mE0BNS6i+4u5fhY/PJFiU5rpwuxHQJjzmzRasIz+t4rgRTgKeCs9/bdgSsDc5wGNhK4ADmdXUeYvgS8yT4euuMbP/C9/XIncvCu+3z6mmjPDcDyX4rIa4e7otkStpSFfuIon1v8BQINa1t7+J2v42/HdXxOOK7fpm1p7gyvS0sCVgJkEyrZGZXQ60BSqa8I0gOXcPf7q4e7xrHUTGbcAFEcdr4+7L3P2vBEt7lgCzzOzUOI5fxnd/ryLP71ZgrrsfRdDCsNu5u/t44BdAFsEXmM4ezrvt7kXhv15dGYCZHULwvt7s7m/HEbNIyii5iyRQeDU7hSDBV/gHwUpsAJcC8/bjJQ4hSKxfm9n3CZbPrJGZHUvwZWBwxBXo28AJZtYxrNPEzI4gaPJuZ2YdwnqDYoxpNvDriitfM+sR/nsY8Im7/x/wPMFthEhvAj8zs4zwSrl3xL5VBLcGAC6IKM8BisPHl9dwzh3cfbG730Gw8E/nGM+jYq3x54BJ7v5srM8TSTUld5HE+39Ay4jtXwNXWLC612XAb/b1wO7+PrCQIBH/laAJek+GA80JVrlbZGYPu/sGgsT4VBjTP4HO7l5K0Aw/M+xQtz7GsG4laB7/ILx9cGtYfhGwxMwWAUcBk6Ke9xzwb4KlTCeFcVS4BbjHzBYAkbct/gDcbmYLqfk242/NbEl4bjuBF2M8j4qYTwIuj+iU1z2O54ukhFaFExERSTO6chcREUkzSu4iIiJpRsldREQkzSi5i4iIpBkldxERkTSj5C4iIpJmlNxFRETSjJK7iIhImvn/AXB9QKLjEuOLAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Useful information on leverage can also be plotted:\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "fig = plot_leverage_resid2(results, ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **4.3.** Fancy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_fancy = {\n",
    "    'x1': X_train['Head Size'],\n",
    "    'x2': X_train['Male'],\n",
    "    'x3': X_train['Female'],\n",
    "    'x4': X_train['Age Range'],\n",
    "    'y': Y_train}\n",
    "\n",
    "# Fit any form of regression model\n",
    "results = ols('y ~ x1 + np.sin(x1) + I(x1-5) + x2 + x3 + np.exp(x4)', data=data_fancy).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.614\n",
      "Model:                            OLS   Adj. R-squared:                  0.606\n",
      "Method:                 Least Squares   F-statistic:                     73.23\n",
      "Date:                Thu, 16 Jun 2022   Prob (F-statistic):           5.07e-37\n",
      "Time:                        19:03:52   Log-Likelihood:                -1080.4\n",
      "No. Observations:                 189   AIC:                             2171.\n",
      "Df Residuals:                     184   BIC:                             2187.\n",
      "Df Model:                           4                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept     27.0967      4.521      5.994      0.000      18.178      36.016\n",
      "x1            67.8632     11.293      6.009      0.000      45.583      90.143\n",
      "np.sin(x1)    -1.5710      7.669     -0.205      0.838     -16.702      13.560\n",
      "I(x1 - 5)    -67.6204     11.310     -5.979      0.000     -89.935     -45.306\n",
      "x2            20.5428      7.787      2.638      0.009       5.179      35.906\n",
      "x3             6.5540      5.675      1.155      0.250      -4.642      17.750\n",
      "np.exp(x4)    13.1935      6.457      2.043      0.042       0.453      25.934\n",
      "==============================================================================\n",
      "Omnibus:                        7.167   Durbin-Watson:                   2.232\n",
      "Prob(Omnibus):                  0.028   Jarque-Bera (JB):                7.105\n",
      "Skew:                           0.386   Prob(JB):                       0.0287\n",
      "Kurtosis:                       3.554   Cond. No.                     4.93e+19\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The smallest eigenvalue is 2.06e-30. This might indicate that there are\n",
      "strong multicollinearity problems or that the design matrix is singular.\n"
     ]
    }
   ],
   "source": [
    "# Inspect the results\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **4.4.** GLM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "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>const</th>\n",
       "      <th>age_range</th>\n",
       "      <th>head_size</th>\n",
       "      <th>male</th>\n",
       "      <th>female</th>\n",
       "      <th>brain_weight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>4158</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1432</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>174</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>3615</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>3458</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1225</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>198</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3493</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1240</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>217</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3170</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1178</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>118</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>4114</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1362</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>230</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3685</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>189</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3703</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>206</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3402</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1220</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>3478</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1270</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>189 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     const  age_range  head_size  male  female  brain_weight\n",
       "34     1.0          1       4158     1       0          1432\n",
       "174    1.0          1       3615     0       1          1350\n",
       "55     1.0          1       3458     1       0          1225\n",
       "198    1.0          0       3493     0       1          1240\n",
       "217    1.0          0       3170     0       1          1178\n",
       "..     ...        ...        ...   ...     ...           ...\n",
       "118    1.0          0       4114     1       0          1362\n",
       "230    1.0          0       3685     0       1          1350\n",
       "189    1.0          0       3703     0       1          1350\n",
       "206    1.0          0       3402     0       1          1220\n",
       "99     1.0          0       3478     1       0          1270\n",
       "\n",
       "[189 rows x 6 columns]"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_glm = pd.concat([X_train, Y_train], axis=1, join=\"inner\")\n",
    "data_glm.columns = ['const', 'age_range', 'head_size', 'male', 'female', 'brain_weight']\n",
    "data_glm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:           brain_weight   No. Observations:                  189\n",
      "Model:                            GLM   Df Residuals:                      185\n",
      "Model Family:                Gaussian   Df Model:                            3\n",
      "Link Function:               identity   Scale:                          5520.6\n",
      "Method:                          IRLS   Log-Likelihood:                -1080.4\n",
      "Date:                Thu, 16 Jun 2022   Deviance:                   1.0213e+06\n",
      "Time:                        19:04:48   Pearson chi2:                 1.02e+06\n",
      "No. Iterations:                     3   Pseudo R-squ. (CS):             0.7894\n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept    262.0164     42.121      6.221      0.000     179.460     344.573\n",
      "head_size      0.2426      0.018     13.737      0.000       0.208       0.277\n",
      "male         138.1129     24.881      5.551      0.000      89.347     186.879\n",
      "female       123.9035     18.699      6.626      0.000      87.253     160.553\n",
      "age_range     22.5626     11.054      2.041      0.041       0.896      44.229\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "# Formula-based\n",
    "formula = 'brain_weight ~ head_size + male + female + age_range'\n",
    "model = glm(formula=formula, data=data_glm, family=sm.families.Gaussian()).fit()\n",
    "print(model.summary())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.10 64-bit (system)",
   "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.9.10"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "9805bccb659ffcb6952b8b29440036fac091f8170d71a5a402ad9a8f2988209f"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}