5_MLR_Determination_Confidence_fkingfisher.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "aa13cd0f-de9a-4a9a-9177-7c53d9aea186",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Exercise 3.12, taken from Montgomery, Douglas C., et al. <u>Introduction to Linear Regression Analysis</u>, John Wiley & Sons, Incorporated, 2012.\n",
    "\n",
    "A chemical engineer studied the effect of the amount of surfactant and time on clathrate formation. Clathrates are used as cool storage media. Table B.8 summarizes the experimental results.\n",
    "\n",
    "**a. Fit a multiple linear regression model relating clathrate formation to these regressors.**\n",
    "\n",
    "Given:\n",
    "\n",
    " * y = Clathrate Formation\n",
    " * x1 = level of surfactant\n",
    " * x2 = time of reaction\n",
    " \n",
    "The model equation of clathrate formation is:\n",
    "\n",
    "    y_hat = beta_0 + beta_1(x1) + beta_2(x2) + beta_12(x1)(x2)\n",
    "\n",
    "The multiple linear regression model is:\n",
    "\n",
    "    y_hat = 12.5013 + 256.7374(x1) + 0.0988(x2) + 0.7613(x1)(x2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "cd90d453-fe1f-44f2-a453-51c7779907ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "const     12.501284\n",
       "x1       256.737401\n",
       "x2         0.098792\n",
       "x1_x2      0.761270\n",
       "dtype: float64"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "\n",
    "# Import clathrate formation data\n",
    "dataset = pd.read_csv(\"https://raw.githubusercontent.com/forestwhite/Regression/main/data/data-table-B8.csv\")\n",
    "# display(dataset.head(8))\n",
    "\n",
    "# Option 1: Add column(s) for interaction terms, uses explicit index names because python doesn't gracefully handle indices\n",
    "# dataset['x1x2'] = dataset['x1']*dataset['x2']\n",
    "# Create Xs and y; for option 1, include term 'x1x2'\n",
    "# Xb = sm.add_constant(dataset[['x1','x2','x1x2']])\n",
    "# y1 = dataset['y']\n",
    "# Generate a linear model \n",
    "# model = sm.OLS(y1, Xb)\n",
    "\n",
    "# Option 2: use scikit-learn's Polynomial features to generate interaction terms\n",
    "# Note: sklearn a.k.a. scikit-learn is not always installed and configured \n",
    "# to operate in jupyter environment. Command line \"conda install scikit-learn\" in \n",
    "# the jupyter environment is often required, along with updates to other libraries.\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "# Create Xs and y;\n",
    "X = dataset[['x1','x2']]\n",
    "y1 = dataset['y']\n",
    "\n",
    "# degree = 2 indicates the second order interactions\n",
    "interaction = PolynomialFeatures(degree=2, interaction_only=True, include_bias=False)\n",
    "\n",
    "X_inter = interaction.fit_transform(X)\n",
    "\n",
    "# The linear model generates sequential x# labels, which hides which interactions are in which variable\n",
    "#X_inter_const = sm.add_constant(X_inter)\n",
    "#model = sm.OLS(y1, X_inter_const).fit()\n",
    "# interaction.get_feature_names_out() is a key for the interactions, used below to \n",
    "#print(interaction.get_feature_names_out())\n",
    "\n",
    "# create a dictionary that maps X using the itertools.count() function, simlar to an enumeration\n",
    "from itertools import count\n",
    "x_to_feature = dict(zip(('x{}'.format(i) for i in count()), X.columns))\n",
    "\n",
    "# substitute x# with column names or combinations of column names from the dataset,\n",
    "# replaces spaces with underscores to ease column references\n",
    "\n",
    "import re\n",
    "# necessary so `x11` enters the loop before `x1`, preventing duplicate x1 match\n",
    "feature_keys = list(x_to_feature.keys())[::-1]\n",
    "features = []\n",
    "\n",
    "for feature in interaction.get_feature_names_out():\n",
    "    feature = re.sub(' ', '_', feature)\n",
    "    features.append(feature)\n",
    "\n",
    "X_all_inter = pd.DataFrame(X_inter, columns=features)\n",
    "X_all_inter_const = sm.add_constant(X_all_inter)\n",
    "\n",
    "lm312 = sm.OLS(y1, X_all_inter_const).fit()\n",
    "display(lm312.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b0ddae9-7fe7-4246-8a4e-0e18da05c52b",
   "metadata": {
    "tags": []
   },
   "source": [
    "**b. Test for significance of regression. What conclusions can you draw?**\n",
    "\n",
    "Besides a few outliers, the variance is constant/random for this data set, as demonstrated in the residual vs. fitted plot.\n",
    "\n",
    "Additionally, the data set demonstrates normality as demonstrated by the normal probability plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "a82ebcb5-6eb6-48d1-abe2-aebb34206c1e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAAF1CAYAAAAJAjeKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAj10lEQVR4nO3df3RkeVnn8ffT04EUECroNDGkGWlRIJJV8GQ4KNgiv8wgMIrigOKCrht/rAq7ogLqgV1R0R0F97irRmFR+TkCKsoQZlCwdRUkjANkCCBMD8xkQqZHTJGBilTsZ/+o6ph0p3vSqVupH/f9Oienb+6t3Pv0d+6kP/Wt594bmYkkSZJUJoe6XYAkSZJ00AzBkiRJKh1DsCRJkkrHECxJkqTSMQRLkiSpdAzBkiRJKh1DsCSdR0R8X0Rcd4Ht742IHyrgOI+LiNva3U8nRcTvRMQvdHD/D4qIjIjDnTqGJG1nCJY0ECLiloioR8RdEfHZiHhtRNynnX1m5usz88lF1dgPIuJlEdFojeOZr5/JzB/JzF9sveac0N76udd1p2pJuniGYEmD5GmZeR/gEcAjgRd3t5y+9ebMvM+2r1/rdkGSVDRDsKSBk5mfBd5FMwwDEBGPjoi/i4i1iPhQRDxu27bnRcTNEbEeEScj4vu2rf/bba97UkR8LCJqEfFbQGzbtmMm9OyP9yPiByJiqXWMmyPih89Xf0T8bEQst1778Yh4wi6veXRrxvuSbeu+MyI+3Fp+VEQsRMTnI2I1In7j4kbxnOO9NiJeHhH3Bt4JPGDbTPH3Ai8Brmp9/6HWz1Qj4tURsdL6+7z8TL0RcUlEXB0Rd0bEzcC3t1OfJF0sQ7CkgRMRR4ErgE+2vp8A3gG8HPgy4IXAWyPiSCvU/S/giswcAb4JuHGXfV4KvBX4eeBS4FPAYy6irDuApwL3BX4AeGVEfMMux3ko8OPA5a16vg245ezXZeb7gC8Aj9+2+nuBN7SWfxP4zcy8L/Bg4JqLqPW8MvMLNMf29m0zxW8Afpl/n0H++tbL/wDYBL6a5sz8k4EzPdT/meZ4PBKYBr67iPokaa8MwZIGyZ9GxDpwK83Q+dLW+ucA12bmtZl5OjOvBxaAp7S2nwamIqKSmSuZedMu+34K8NHMfEtmNoBXAZ/da2GZ+Y7M/FQ2/TVwHfDNu7z034B7Al8bEUOZeUtmfuo8u30j8GyAiBhp1fjG1rYG8NURcWlm3tUKzXv1Pa0Z8zNfD7iIn6VVzxjNsPyCzPxCZt4BvBJ41pljAK/KzFsz83PAr1zsMSSpHYZgSYPkO1qzp48DHkZzxhbgK4Fnbg92wGOB8dbM5lXAjwArEfGOiHjYLvt+AM1wDUBm5vbv705EXBER74uIz7WO/5Rt9W3JzE8CLwBeBtwREW+6QAh9A/CMiLgn8Azghsz8dGvbfwIeAnwsIj4QEU/da63ANZk5uu3r9ov42TO+EhiiOaZnxvx3gfu3tu8YT+DTSNIBMgRLGjitmdbXAle3Vt0K/NFZwe7emfmK1uvflZlPAsaBjwG/t8tuV4AHnvkmImL79zRbE+617fuv2Pbae9JspbgaGMvMUeBatvUUn1X/GzLzsTSDZAK/ep7XfZRmeLyCna0QZOY/ZeazaYbOXwXe0mr9KELuYd2twL8Cl24b8/tm5sNb23eMJ3BZQbVJ0p4YgiUNqlcBT4qIRwCvA54WEd/WuiBruHWbr6MRMRYRT28FxH8F7qLZknC2dwAPj4hntC52+0m2BV2afcTHI+KyiKiy884U96DZ4nAK2IyIK2j2x54jIh4aEY9vBecNoH6ees54Q6uW48Afb9vPcyLiSGaeBtZaqy+0n4uxCnx56++5fd2DIuIQQGau0Gz5+PWIuG9EHIqIB0fEt7Refw3wk63/BvcDXlRQbZK0J4ZgSQMpM08Bfwj8QmbeClxJ8w4Gp2jOUv40zd+Bh4CfAm4HPgd8C/Bju+zvTuCZwCuAfwa+Bvh/27ZfD7wZ+DDwQeAvtm1bpxlUrwH+heas7dvPU/o9W8e4k2bP8f1bdZ/PG2m2f/xVq8YzZoCbIuIumhfJPSszNwBad3DYrR95TzLzY63j3rytZ/hMAP/niLihtfwfab4B+CjNv/dbaM62Q3O2/V3Ah4AbgLfttx5J2o9otrVJkiRJ5eFMsCRJkkrHECxJkqTSMQRLkiSpdAzBkiRJKh1DsCRJkkrncDcOOjMzk/Pz8904tCRJkspl1wcTdWUm+M4777z7F0mSJEkdYjuEJEmSSscQLEmSpNIxBEuSJKl0DMGSJEkqHUOwJEmSSscQLEmSpNIxBEuSJKl0DMGSJEkqHUOwJEmSSqcrj02WJElSZyyt1JhfXGV5rc7EaIWZqTEmx6vdLqvnOBMsSZI0IJZWasydOEmt3mC8Okyt3mDuxEmWVmrdLq3nGIIlSZIGxPziKtXKENXKEIcitpbnF1e7XVrPMQRLkiQNiOW1OiPDO7tdR4YPs7xW71JFvcsQLEmSNCAmRiusb2zuWLe+scnEaKVLFfUuQ7AkSdKAmJkao1ZvUKs3OJ25tTwzNdbt0nqOIViSJGlATI5XmT1+jGpliJXaBtXKELPHj3l3iF14izRJkqQBMjleNfTugTPBkiRJKh1DsCRJkkrHECxJkqTSMQRLkiSpdLwwTpIkSedYWqkxv7jK8lqdidEKM1NjA3XBnTPBkiRJ2mFppcbciZPU6g3Gq8PU6g3mTpxkaaXW7dIKYwiWJEnSDvOLq1QrQ1QrQxyK2FqeX1ztdmmFMQRLkiRph+W1OiPDO7tmR4YPs7xW71JFxSskBEfEf42ImyJiMSLeGBHDRexXkiRJB29itML6xuaOdesbm0yMVrpUUfHaDsERMQH8JDCdmVPAJcCz2t2vJEmSumNmaoxavUGt3uB05tbyzNRYt0srTFHtEIeBSkQcBu4F3F7QfiVJknTAJserzB4/RrUyxEptg2pliNnjxwbq7hBt3yItM5cj4mrgM0AduC4zrzv7dRExC8wCXHbZZe0eVpIkSR00OV4dqNB7tiLaIe4HXAkcAx4A3DsinnP26zJzLjOnM3P6yJEj7R5WkiRJ2rci2iGeCJzMzFOZ2QDeBnxTAfuVJEmSOqKIEPwZ4NERca+ICOAJwFIB+5UkSZI6ou0QnJnvB94C3AB8pLXPuXb3K0mSJHVK2xfGAWTmS4GXFrEvSZIkqdN8YpwkSZJKxxAsSZKk0jEES5IkqXQMwZIkSSodQ7AkSZJKxxAsSZKk0jEES5IkqXQMwZIkSSodQ7AkSZJKxxAsSZKk0jEES5IkqXQMwZIkSSodQ7AkSZJKxxAsSZKk0jEES5IkqXQMwZIkSSodQ7AkSZJKxxAsSZKk0jEES5IkqXQMwZIkSSodQ7AkSZJKxxAsSZKk0jEES5IkqXQMwZIkSSodQ7AkSZJKxxAsSZKk0jEES5IkqXQMwZIkSSqdw0XsJCJGgd8HpoAEfjAz/76IfUuSJKk/La3UmF9cZXmtzsRohZmpMSbHq90uCyhuJvg3gfnMfBjw9cBSQfuVJElSH1paqTF34iS1eoPx6jC1eoO5EydZWql1uzSggBAcEfcFjgOvBsjML2XmWrv7lSRJUv+aX1ylWhmiWhniUMTW8vziardLA4qZCf4q4BTwfyPiHyPi9yPi3me/KCJmI2IhIhZOnTpVwGElSZLUq5bX6owM7+y8HRk+zPJavUsV7VRECD4MfAPw25n5SOALwIvOflFmzmXmdGZOHzlypIDDSpIkqVdNjFZY39jcsW59Y5OJ0UqXKtqpiBB8G3BbZr6/9f1baIZiSZIkldTM1Bi1eoNavcHpzK3lmamxbpcGFBCCM/OzwK0R8dDWqicAH213v5IkSepfk+NVZo8fo1oZYqW2QbUyxOzxYz1zd4hCbpEG/ATw+oi4B3Az8AMF7VeSJEl9anK82jOh92yFhODMvBGYLmJfkiRJUqf5xDhJkiSVjiFYkiRJpWMIliRJUukUdWGcJEnSgVtaqTG/uMryWp2J0QozU2M9eyGWeoszwZIkqS8trdSYO3GSWr3BeHWYWr3B3ImTLK3Uul2a+oAzwSXmu2dJUj+bX1ylWhmiWhkC2PpzfnHVf890twzBfardAHvm3XO1MrTj3XMv3cRakqQLWV6rM14d3rFuZPgwy2v1LlWkfmII7kNFBFjfPatf+QmGpDMmRivU6o2tf8MA1jc2mRitdLEq9Qt7gvvQ9gB7KGJreX5xdc/7WF6rMzK88z2Q757V6+z/k7TdzNQYtXqDWr3B6cyt5ZmpsW6Xpj5gCO5DRQTYidEK6xubO9b57lm9rog3gJIGx+R4ldnjx6hWhlipbVCtDNnWpz2zHaIPFfHxz8zUGHMnTgLNAL2+sUmt3uCqy48WXq9UFPv/JJ1tcrxq6NW+OBPch4r4+Md3z+pHfoIhSSqKM8F96EyA3X5x0FWXH73oAOu7Z/UbP8GQJBUlMvPADzo9PZ0LCwsHflxJ/c+7Q0iSLlLsttKZYEl9xU8wJElFsCdYkiRJpeNMsNriR9OSJKkfOROsffPBBZIkqV8ZgrVvPrhAkiT1K0Ow9s1HL0uSpH5lCNa++eACSZLUrwzB2rcinlwnSZLUDYZg7ZuPXpYkSf3KW6SpLT64QHvl7fQkSb3EmWBJHeft9CRJvcYQLKnjvJ2eJKnXGIIldZy305Mk9RpDsKSO83Z6kqReU1gIjohLIuIfI+IvitqnpMHg7fQkSb2myJng5wNLBe5P0oDwdnqSpF5TyC3SIuIo8O3ALwH/rYh9Shos3k5PktRLipoJfhXwM8DpgvYnSZIkdUzbITgingrckZkfvJvXzUbEQkQsnDp1qt3DSpIkSftWxEzwY4CnR8QtwJuAx0fE685+UWbOZeZ0Zk4fOXKkgMNKkiRJ+9N2CM7MF2fm0cx8EPAs4K8y8zltVyZJkiR1iPcJliRJUukUcneIMzLzvcB7i9ynJEmSVDRngiVJklQ6hmBJkiSVjiFYkiRJpWMIliRJUukYgiVJklQ6hmBJkiSVjiFYkiRJpWMIliRJUukYgiVJklQ6hmBJkiSVjiFYkiRJpWMIliRJUukYgiVJklQ6hmBJkiSVjiFYkiRJpXO42wVIAEsrNeYXV1leqzMxWmFmaozJ8Wq3y5IkSQPKmWB13dJKjbkTJ6nVG4xXh6nVG8ydOMnSSq3bpUmSpAFlCFbXzS+uUq0MUa0McShia3l+cbXbpUmSpAFlCFbXLa/VGRne2ZkzMnyY5bV6lyqSJEmDzp5gdd3EaIVavUG1MrS1bn1jk4nRSherGhz2W0uSdC5ngtV1M1Nj1OoNavUGpzO3lmemxrpdWt+z31qSpN0ZgtV1k+NVZo8fo1oZYqW2QbUyxOzxY85WFsB+a0mSdmc7hHrC5HjV0NsBy2t1xqvDO9bZby1JkjPB0kCbGK2wvrG5Y5391pIkGYKlgWa/tSRJuzMESwPMfmtJknZnT7A04Irot/Y2a5KkQeNMsKQL8jZrkqRBZAiWdEHeZk2SNIjaDsER8cCIeE9ELEXETRHx/CIKk9QbfKy1JGkQFdETvAn8VGbeEBEjwAcj4vrM/GgB+5bUZT7Wujzs/e4cx1bqPW3PBGfmSmbe0FpeB5aAiXb3K6k3eJu1crD3u3McW6k3FdoTHBEPAh4JvH+XbbMRsRARC6dOnSrysJI6yNuslYO9353j2Eq9qbBbpEXEfYC3Ai/IzM+fvT0z54A5gOnp6SzquJI6z8daDz4fsd05jq3UmwqZCY6IIZoB+PWZ+bYi9ilJOjg+YrtzHFupNxVxd4gAXg0sZeZvtF+SJOmg2fvdOY6t1JuKmAl+DPD9wOMj4sbW11MK2K8k6YDY+905jq3UmyLz4Ntzp6enc2Fh4cCPK0mSpNKJ3Vb6xDhJkiSVjiFYkiRJpWMIliRJUukYgiVJklQ6hmBJkiSVjiFYkiRJpWMIliRJUukYgiVJklQ6h7tdQNksrdSYX1xlea3OxGiFmakxnxokSZJ0wJwJPkBLKzXmTpykVm8wXh2mVm8wd+IkSyu1bpcmSZJUKobgAzS/uEq1MkS1MsShiK3l+cXVbpcmSZJUKobgA7S8VmdkeGcHysjwYZbX6l2qSJIkqZzsCT5AE6MVavUG1crQ1rr1jU0mRitdrEqSJPUKrx06OM4EH6CZqTFq9Qa1eoPTmVvLM1Nj3S5NkiR1mdcOHSxngg/Q5HiV2ePHdrzDu+ryo77Dkw6YMy2SetH2a4eArT/nF1f9HdUBhuADNjle9USWuujMTEu1MrRjpmX2+DH/35TUVctrdcarwzvWee1Q59gOIalUvEuLpF41MVphfWNzxzqvHeocZ4IllYozLZ1lq4m0fzNTY8ydOAk0fy+tb2xSqze46vKjXa5sMDkTLKlUnGnpHC/qkdpz5tqhamWIldoG1cqQrVod5EywpFJxpqVzvKhHap/XDh0cZ4IllYozLZ3jA4Ek9RNngiWVjjMtneEDgST1E2eCJUmF8IFAkvqJM8EaGF6VLnWXDwSS1E8iMw/8oNPT07mwsHDgx9Xg2v4AhO0XO9nrKUlS6cVuK22H0EDwAQiSJOliGII1ELwqXZIkXQxDsAaCD0CQJEkXo5AQHBEzEfHxiPhkRLyoiH1KF8Or0iVJ0sVoOwRHxCXA/wauAL4WeHZEfG27+5Uuhg9AkCRJF6OIW6Q9CvhkZt4MEBFvAq4EPlrAvov15sedu+6h3wOP+DFofBHe9pRztz/8eTD1PPjinfDn333u9q//UXjYVfD5W+Gd33/u9umfggc/DT73cbj+h8/d/uifh698ItxxI7znBeduf+wvw8Q3wfLfwd++5Nzt3/oquP8j4NPvhve9/NztT/pd+LKHwqf+HBZ+/dztV/wR3PeB8LE3w4d++9ztT3sL3OtSWHwt3PTac7c/41oYuhfc+H/g49ecu/2q9zb//MDVcPNf7Nx2uALf9c7m8t//InzmL3dur3w5PP2tzeW/eTHc/vc7t48chae8rrn8nhcweceNTJ7Z9jngIw+B8bnm99fNwr98YufP3/8RzfEDuPY5sH7bzu0P+Eb45l9pLr/9u6D+zzu3X/YE+MZfaC6/9QrYPKv/+KueCpe/sLnsuXfu9gE697jjxp3b7/cQeLLnnuee5945PPf689w7s98BU0Q7xARw67bvb2ut2yEiZiNiISIWTp06VcBhJUmSpP1p+z7BEfFM4Nsy84da338/8KjM/Inz/Yz3CZYkSdIB6dh9gm8DHrjt+6PA7QXsV5IkSeqIIkLwB4CviYhjEXEP4FnA2wvYryRJktQRbV8Yl5mbEfHjwLuAS4DXZOZNbVcmSZIkdUgRd4cgM68Fri1iX5IkSVKn+cQ4SZIklY4hWJIkSaVjCJYkSVLpGIIlSZJUOoZgSZIklY4hWJIkSaVjCJYkSVLpGIIlSZJUOoZgSZIklU4hT4yTJEnqV0srNeYXV1leqzMxWmFmaozJ8Wq3y1KHORMsSZJKa2mlxtyJk9TqDcarw9TqDeZOnGRppdbt0tRhhmBJklRa84urVCtDVCtDHIrYWp5fXO12aeowQ7AkSSqt5bU6I8M7u0NHhg+zvFbvUkU6KPYESz3KHjVJ6ryJ0Qq1eoNqZWhr3frGJhOjlS5WpYPgTLDUg+xRk6SDMTM1Rq3eoFZvcDpza3lmaqzbpanDDMFSD7JHTZIOxuR4ldnjx6hWhlipbVCtDDF7/JifvJWA7RBSD1peqzNeHd6xzh41SeqMyfGqobeEnAmWetDEaIX1jc0d6+xRkySpOIZgqQfZoyZJUmcZgqUeZI+aJEmdZU+w1KPsUZMkqXOcCZYkSVLpGIIlSZJUOrZDSGfxSW2SJA0+Z4KlbXxSmyRJ5WAIlrbxSW2SJJWDIVjaZnmtzsjwzi4hn9QmSdLgMQRL2/ikNkmSyqGtEBwR/zMiPhYRH46IP4mI0YLqkrrCJ7VJklQO7c4EXw9MZebXAZ8AXtx+SVL3+KQ2SZLKoa1bpGXmddu+fR/w3e2VI3WfT2pTv/G2fpJ08YrsCf5B4J3n2xgRsxGxEBELp06dKvCwklRe3tZPkvbnbkNwRLw7IhZ3+bpy22t+DtgEXn++/WTmXGZOZ+b0kSNHiqlekkrO2/pJ0v7cbTtEZj7xQtsj4rnAU4EnZGYWVZgk6e4tr9UZrw7vWOdt/STp7rXVExwRM8DPAt+SmV8spiRJ0l5NjFao1RtUK0Nb6wbhtn72OUvqtHZ7gn8LGAGuj4gbI+J3CqhJkrRHg3hbP/ucJR2Edu8O8dVFFSJJunhnbuu3fdb0qsuP9vWs6fY+Z2Drz/nF1b7+e0nqLW2F4LLx4zlJvWjQbutnn7Okg+Bjk/fIj+ck6WD4+HJJB8EQvEfehkiSDsYg9jlL6j2G4D1aXqszMryze8SP5ySpeD6+XNJBsCd4jwb1NkTqDPvHpfYMWp+zpN7jTPAe+fGc9sr+cUmSep8heI/8eE57Zf+4JEm9z3aIi+DHc9oLb+8kSVLvcyZYKpi3d5IkqfcZgqWC2T8uSVLvMwRLBbN/XJKk3mdPsNQB9o9LktTbnAmWJElS6RiCJUmSVDqGYEmSJJWOIViSJEmlYwiWJElS6RiCJUmSVDqGYEmSJJWOIViSJEml48MyJEnSgVtaqTG/uMryWp2J0QozU2M+ZEgHyplgSZJ0oJZWasydOEmt3mC8Okyt3mDuxEmWVmrdLk0lYgiWJEkHan5xlWpliGpliEMRW8vzi6vdLk0lYgiWJEkHanmtzsjwzo7MkeHDLK/Vu1SRysgQLEmSDtTEaIX1jc0d69Y3NpkYrXSpIpWRIViSJB2omakxavUGtXqD05lbyzNTY90uTSViCJYkSQdqcrzK7PFjVCtDrNQ2qFaGmD1+zLtD6EB5izRJknTgJserhl51VSEzwRHxwojIiLi0iP1JkiRJndR2CI6IBwJPAj7TfjmSJElS5xUxE/xK4GeALGBfkiRJUse1FYIj4unAcmZ+qKB6JEmSpI672wvjIuLdwFfssunngJcAT97LgSJiFpgFuOyyyy6iREnqTUsrNeYXV1leqzMxWmFmaswLfSSpT0Tm/roYIuI/AH8JfLG16ihwO/CozPzshX52eno6FxYW9nVcSeoFSys15k6cpFoZYmT4MOsbm9TqDW/zJEm9J3Zbue9bpGXmR4D7b+094hZgOjPv3O8+JalfzC+uUq0MUa0MAWz9Ob+4agiWpD7gfYIlaR+W1+qMV4d3rBsZPszyWr1LFWnQ2X4jFauwEJyZDypqX5LU6yZGK9Tqja0ZYID1jU0mRisXtR+DjfZie/vNeHWYWr3B3ImTtt9IbfCxyZK0DzNTY9TqDWr1Bqczt5Znpsb2vI8zwaZWb+wINksrtQ5Wrn60vf3mUMTW8vziardLk/qWIViS9mFyvMrs8WNUK0Os1DaoVoYuelbOYKO9Wl6rMzK888Nb22+k9tgTLEn7NDlebeujaPuKtVdFtd9I+nfOBEtSl0yMVljf2NyxzmCj3RTRfiNpJ0OwJHWJwUZ7VUT7jaSd9v2wjHb4sAxJavLuEJLUccU+LEOS1L52+4olSftjO4QkSZJKxxAsSZKk0jEES5IkqXQMwZIkSSodQ7AkSZJKxxAsSZKk0jEES5IkqXQMwZIkSSodQ7AkSZJKxxAsSZKk0jEES5IkqXQMwZIkSSodQ7AkSZJKxxAsSZKk0jEES5IkqXQMwZIkSSodQ7AkSZJK53C3C5AkSQdnaaXG/OIqy2t1JkYrzEyNMTle7XZZ0oFzJliSpJJYWqkxd+IktXqD8eowtXqDuRMnWVqpdbs06cAZgiVJKon5xVWqlSGqlSEORWwtzy+udrs06cAZgiVJKonltTojwzs7IUeGD7O8Vu9SRVL32BMsSdJ5DFr/7MRohVq9QbUytLVufWOTidFKF6uSuqPtmeCI+ImI+HhE3BQRv1ZEUUVbWqnxyus/wQv/+EO88vpP2PskSbpbg9g/OzM1Rq3eoFZvcDpza3lmaqzbpUkHrq0QHBHfClwJfF1mPhy4upCqCjSIv8QkSZ03iP2zk+NVZo8fo1oZYqW2QbUyxOzxY309uy3tV7vtED8KvCIz/xUgM+9ov6Ribf8lBmz9Ob+46v/0kqTzWl6rM14d3rFuEPpnJ8er/vsn0X47xEOAb46I90fEX0fE5ed7YUTMRsRCRCycOnWqzcPunRcBSJL2Y2K0wvrG5o519s9Kg+NuQ3BEvDsiFnf5upLmTPL9gEcDPw1cExGx234ycy4zpzNz+siRI4X+JS7EX2KSpP2wf1YabHcbgjPziZk5tcvXnwG3AW/Lpn8ATgOXdrroi+EvMUnSftg/Kw22dnuC/xR4PPDeiHgIcA/gznaLKtKZX2Lbb3Fz1eVH/SUmSbpb9s9Kg6vdEPwa4DURsQh8CXhuZmb7ZRXLX2KSJEnarq0QnJlfAp5TUC2SJEnSgfCxyZIkSSodQ7AkSZJKp92eYEmSetLSSm3HRdEzU2NeHyJpizPBkqSBs7RSY+7ESWr1BuPVYWr1BnMnTrK0Uut2aZJ6hCFYkjRw5hdXqVaGqFaGOBSxtTy/uNrt0iT1CEOwJGngLK/VGRne2fE3MnyY5bV6lyqS1GsMwZKkgTMxWmF9Y3PHuvWNTSZGK12qSFKvMQRLkgbOzNQYtXqDWr3B6cyt5ZmpsW6XJqlHGIIlSQNncrzK7PFjVCtDrNQ2qFaGmD1+zLtDSNriLdIkSQNpcrxq6JV0Xs4ES5IkqXScCZYkSRfFB5HsznHpL84ES5KkPfNBJLtzXPqPIViSJO2ZDyLZnePSfwzBkiRpz3wQye4cl/5jCJYkSXvmg0h257j0H0OwJEnaMx9EsjvHpf8YgiVJ0p75IJLdOS79JzLzwA86PT2dCwsLB35cSZIklU7sttKZYEmSJJWOIViSJEmlYwiWJElS6RiCJUmSVDqGYEmSJJWOIViSJEmlYwiWJElS6RiCJUmSVDqGYEmSJJWOIViSJEml05XHJkfEKeDTB35guBS4swvHLQPHtnMc285xbDvHse0cx7ZzHNvO6ebY3pmZM2ev7EoI7paIWMjM6W7XMYgc285xbDvHse0cx7ZzHNvOcWw7pxfH1nYISZIklY4hWJIkSaVTthA81+0CBphj2zmObec4tp3j2HaOY9s5jm3n9NzYlqonWJIkSYLyzQRLkiRJ5QnBETETER+PiE9GxIu6Xc8giYhbIuIjEXFjRCx0u55+FhGviYg7ImJx27ovi4jrI+KfWn/er5s19qvzjO3LImK5de7eGBFP6WaN/SoiHhgR74mIpYi4KSKe31rvudumC4yt526bImI4Iv4hIj7UGtv/3lrveduGC4xrz52zpWiHiIhLgE8ATwJuAz4APDszP9rVwgZERNwCTGem91ZsU0QcB+4C/jAzp1rrfg34XGa+ovUG7n6Z+bPdrLMfnWdsXwbclZlXd7O2fhcR48B4Zt4QESPAB4HvAJ6H525bLjC234PnblsiIoB7Z+ZdETEE/C3wfOAZeN7u2wXGdYYeO2fLMhP8KOCTmXlzZn4JeBNwZZdrks6RmSeAz521+krgD1rLf0DzH0BdpPOMrQqQmSuZeUNreR1YAibw3G3bBcZWbcqmu1rfDrW+Es/btlxgXHtOWULwBHDrtu9vw18iRUrguoj4YETMdruYATSWmSvQ/AcRuH+X6xk0Px4RH261S/ixZ5si4kHAI4H347lbqLPGFjx32xYRl0TEjcAdwPWZ6XlbgPOMK/TYOVuWEBy7rOvJdyV96jGZ+Q3AFcB/aX3sLPWD3wYeDDwCWAF+vavV9LmIuA/wVuAFmfn5btczSHYZW8/dAmTmv2XmI4CjwKMiYqrLJQ2E84xrz52zZQnBtwEP3Pb9UeD2LtUycDLz9tafdwB/QrP9RMVZbfUFnukPvKPL9QyMzFxt/bI+Dfwenrv71ur9eyvw+sx8W2u1524Bdhtbz91iZeYa8F6afauetwXZPq69eM6WJQR/APiaiDgWEfcAngW8vcs1DYSIuHfrYg0i4t7Ak4HFC/+ULtLbgee2lp8L/FkXaxkoZ/6ha/lOPHf3pXUhzKuBpcz8jW2bPHfbdL6x9dxtX0QciYjR1nIFeCLwMTxv23K+ce3Fc7YUd4cAaN2K41XAJcBrMvOXulvRYIiIr6I5+wtwGHiDY7t/EfFG4HHApcAq8FLgT4FrgMuAzwDPzEwv8LpI5xnbx9H8aC6BW4AfPtMLqL2LiMcCfwN8BDjdWv0Smr2rnrttuMDYPhvP3bZExNfRvPDtEpqTgtdk5v+IiC/H83bfLjCuf0SPnbOlCcGSJEnSGWVph5AkSZK2GIIlSZJUOoZgSZIklY4hWJIkSaVjCJYkSVLpGIIlSZJUOoZgSZIklY4hWJIkSaXz/wGdjNt7xKYIvAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtkAAAFzCAYAAADrIhWLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABHuElEQVR4nO3deZyVc//H8denEQ1iUKIhJclODCJue2Ut5bbd9iUhobSMcnP/Ki2j0KJNKFtFSYiE276W0iJDu2YilanUtE3f3x/fM+6pZjkzc865zpl5Px+PHnPOda5zXZ+u+3L37tvn+n7NOYeIiIiIiEROlaALEBERERGpaBSyRUREREQiTCFbRERERCTCFLJFRERERCJMIVtEREREJMIUskVEREREImy3oAuIhho1ari6desGXYaIiIiIVBTOwapV8NtvsGUL7LknMzZuXOWcq1nY7hUyZNetW5fp06cHXYaIiIiIJLqNG2HkSOjXD7KzoXFjeOQRuPhirEqVpUV9rUKGbBERERGRclm/HoYOhf79YeVKOOccGDMGzj8fzEr8ukK2iIiIiEi+nBwYPBiefBLWrIGmTaF7dzj77FIdRiFbRERERGT1anjqKRg4ENatg8svh27d4PTTy3S4wGYXMbNDzey/ZjbfzOaZ2f2F7GNmNtDMFpjZbDM7OYhaRURERKSC+v136NwZDjsMevaEiy6CmTNh8uQyB2wIdiR7G9DROfe9mVUHZpjZNOfcjwX2uRhoEPp1OjA09FNEREREpOyysiAjA0aMgM2b4dpr4eGH4dhjI3L4wEK2c24FsCL0er2ZzQdSgYIhuwUwxjnngK/NLMXMDg59V0RERESkdJYuhT594LnnIC8PbrwR0tPhyCMjepq46Mk2s7pAI+CbnT5KBX4t8H55aNsuIdvM2gBtAOrUqROVOkVEREQkQS1YAI8/Di++6GcHue026NIF6tWLyukCX/HRzPYGJgAPOOfW7fxxIV9xhR3HOTfCOZfmnEurWbPQOcFFREREpLKZPx9uuAEaNoRXX4V77oFFi2DYsKgFbAh4JNvMquID9svOuYmF7LIcOLTA+0OA7FjUJiIiIiIJ7Icf/IOMEybAnntChw7QsSMcdFBMTh9YyDYzA0YB851zA4rYbTLQzszG4h94XKt+bBERERHJN2lmFhlTM8nOyaV2SjK9Dsnl3NdG+NlB9tnHP8z4wANQo0ZM6wpyJLsJcCMwx8xmhbY9DNQBcM4NA6YAlwALgI3ArbEvU0RERETi0aSZWaRPnEPu1jxOWf4j7ceP5ZzF37Nln33Z/T//gfbtISUlkNqCnF3kcwrvuS64jwPujU1FIiIiIpJIMt77iZMWfE/7L8dyxrI5rNpzX/qccwsfnteaaf++LNDa4mJ2ERERERGRsDkHU6fy9JAOpGXN5/e996fH+XfwyonNyd29GpYbdIEK2SIiIiKV0s69zJ2aNaRlo9Sgyyre9u3w1lv+gcbp0zlk3wPpftHdvHbCRWzebfe/d6udkhxgkV7gU/iJiIiISGzl9zJn5eTigKycXNInzmHSzKygSytcXh6MHw+NGkHLlrBmDTz7LN+8/zUTTr9ih4CdXDWJTs0aBldriEK2iIiISCWTMTWT3K15O2zL3ZpHxtTMgCoqwrZt8NJLcNxxcM01fvnzMWMgMxNuv50Wp9Wjd6vjSU1JxoDUlGR6tzo+Lkbk1S4iIiIiUslk5xTetFzU9pjbssWvzNi7NyxcCMcfD+PGQevWkJS0w64tG6XGRajemUayRURERCqZonqWA+9l3rQJhg6FBg3gjjv89HuTJsGsWXD11bsE7HimkC0iIiJSyXRq1pDkqjsG1kB7mTduhKeegvr1/bLnqakwZQp89x20aAFVEi+yql1EREREpJLJb68IfHaR9evhmWegf3/44w8491zfJnLeeWDFLqcS9xSyRURERCqhQHuZc3Jg0CA/er1mDTRrBt27w1lnBVNPFChki4iIiEhsrFrlg/WgQbBuHVx+uQ/Xp50WdGURp5AtIiIiItH122++JWToUN9/3bo1dOsGJ50UdGVRo5AtIiIiEmMJudpiWWRlQb9+MGKEn5bv2mt9uD7mmKArizqFbBEREZEYyl9tMX8xmPzVFoGKE7SXLIE+feD55/1S6DfeCOnpfmq+SiLx5kMRERERSWAJs9piWfzyC9x2mw/Tzz/vX//yCzz3XKUK2KCRbBEREZEdRLuVI+5XWyyLH3+EXr1g7FjYfXc/13WnTnDIIUFXFhiFbBEREZGQWLRy1E5JJquQQB34aotlMWsW9OwJEyfCnntCx47+V61aQVcWOLWLiIiIiITEopUj7lZbLItvv4UrroBGjWDaNHj4Yd+H3a+fAnaIRrJFREREQmLRyhE3qy2WxeefQ48e8P77sP/+8H//B/fdBykpQVcWdxSyRUREREJi1coR6GqLpeUcfPSRD9effAI1a0LfvnD33VC9etDVxS21i4iIiIiEVIhWjkhxDt59F5o0gQsvhJ9/hief9G0hnTsrYJdAI9kiIiIiIQndyhEp27fD5Mn+gcYZM6BOHXjmGbj1VqhWLejqEoZCtoiIiEgBCdXKEUl5eTBhgg/Xc+ZA/fowahTccIOflk9KRe0iIiIiIpXZtm3w4otw3HFwzTWwdat//9NPfjEZBewyUcgWERERqYy2bIFnn4WGDeGmm3yYHj8e5s71o9e7qeGhPHT1REREJGFFe3XGCmnTJr/Med++sGwZpKXBgAFw+eVQReOvkaKQLSIiIgkpFqszVigbN8Lw4ZCRAStWwJln+vfNmoFZ0NVVOPrrioiIiCSkWKzOWCGsXw99+kDdutChAxx1FHz4oV9YpnlzBewo0Ui2iIiIJKRYrM6Y0HJyYOBAeOop+PNPP2L9yCN+3muJOoVsERERSUixWp0x4axa5ReNGTwY1q2DK66A7t3h1FODrqxSCbRdxMyeM7OVZja3iM/PNbO1ZjYr9Ovfsa5RRERE4pNWZ9zJb79Bp06+LaR3bz9yPWsWvPmmAnYAgh7JfgEYDIwpZp/PnHOXxaYcERERSRRanTFk+XLo1w9GjvTT8l13HTz8MBxzTNCVVWqBhmzn3KdmVjfIGkRERCRxVdrVGQEWL/YPND7/PDjn57pOT4cjjgi6MiExZhc5w8x+MLN3zezYonYyszZmNt3Mpv/xxx+xrE9EREQkdn75BW69FRo0gBdegNtv99tGjVLAjiNBt4uU5HvgMOfcX2Z2CTAJaFDYjs65EcAIgLS0NBezCkVERERiYd486NULxo3zqzO2a+d7sFMr6Uh+nIvrkWzn3Drn3F+h11OAqmZWI+CyRERERGJn5kxo3RqOOw4mT4aHHoIlS/zUfArYcSuuR7LN7CDgd+ecM7PT8H8pWB1wWSIiIiLR9+230KMHvP027LOPn4bvgQfggAOCrkzCEGjINrNXgXOBGma2HHgUqArgnBsGXAXcbWbbgFzgWuecWkFERESk4vrsMx+up02D/ff3r9u1g5SUoCuTUgh6dpHrSvh8MH6KPxEREalAJs3M0tR7BTkHH33kA/Unn8CBB0LfvnD33VC9etDVSRnEdbuIiIiIVDyTZmaRPnEOuVvzAMjKySV94hyAyhe0nYN33/Xh+uuvoXZt32t9552w555BVyflENcPPoqIiEjFkzE18++AnS93ax4ZUzMDqigA27fDG29AWhpceimsWAFDh8KiRXD//QrYFYBCtoiIiMRUdk5uqbZXKHl5fgq+E0+EVq1g7Vo/v/Uvv0DbtrDHHkFXKBGikC0iIiIxVTsluVTbK4Rt22DMGDj2WLj2Wh+2X3oJfvoJbrsNqlYNukKJMIVsERERialOzRqSXDVph23JVZPo1KxhQBVF0ZYtMHIkNGwIN98M1arBa6/B3Lnwr3/Bbno8rqLS/7IiIiISU/kPN1bo2UU2bfJtIH37wq+/+t7rJ5+Eyy8Hs6CrkxhQyBYREZGYa9kotWKF6nwbNsDw4ZCRAb/9Bk2a+JHspk0VrisZhWwRERGR8lq3Dp55Bvr3h1Wr4Pzz4ZVX4NxzFa4rKYVsERERkbL6808YOBCeftq/bt4cHnkEzjwz6MokYArZIiIiIqW1ahUMGACDB8P69dCiBXTv7nuvRVDIFhEREQnfb7/BE0/4hWNyc+Gqq3y4PuGEoCuTOKOQLSIiIlKSX3+Ffv38Q4xbt8L118PDD8PRRwddmcQphWwRERGRoixeDL17wwsvgHN+ruuuXeGII4KuTOKcQraIiIjIzn7+GR5/3K/KmJQEd9wBXbrAYYcFXZkkCIVsERERkXxz50KvXjB+POyxB9x3H3TqBLVrB12ZJBiFbBEREZGZM6FnT5g4EfbeGx56CDp0gFq1gq5MEpRCtoiIiFRe33wDPXrAO+/Avvv6Oa7vvx8OOCDoyiTBKWSLiIhI5fPpp37keto02H9//7pdOx+0RSJAIVtEREQqB+fgww/9yPWnn8KBB/pp+e6+27eIiESQQraIiIhUbM7BlCk+XH/zDaSm+mXQ77wTkpODrk4qqCpBFyAiIiISFdu3+wcZTzkFLrvMr9Y4bBgsXAjt2ytgS1QpZIuIiEjFkpcHY8fCiSdC69awfj089xz88gvcdZefmk8kyhSyRUREpGLYuhVGj4ZjjoHrrvMj2S+/DPPnw623QtWqQVcolYh6skVERCSxbdniw3Xv3n4Z9BNPhNdeg1atoIrGEyUYuvNEREQkMeXmwuDBUL8+tGkDNWrA5Ml+YZmrrlLAlkBpJFtEREQSy4YN/gHGJ57wDzOedRaMGgUXXQRmQVcnAihki4iISKJYtw6GDIEBA2DVKjj/fHj1VTjnHIVriTsK2SIiIhLf/vzTz2v99NOQkwMXXwzdu8OZZwZdmUiRFLJFREQkPv3xhx+1HjLET8PXsiV06wZpaUFXJlKiQJ8IMLPnzGylmc0t4nMzs4FmtsDMZpvZybGuUURERGJsxQro2BHq1oW+ff3I9Q8/wBtvKGBLwgj6sdsXgObFfH4x0CD0qw0wNAY1iYiISBB+/RXatYN69XxrSOvW8OOPMG4cnHBC0NWJlEqg7SLOuU/NrG4xu7QAxjjnHPC1maWY2cHOuRWxqVBERESibtEi6NMHXngBnIObb4b0dD81n0iCivee7FTg1wLvl4e27RKyzawNfrSbOnXqxKQ4ERERKYfMTHj8cb8q4267wZ13QpcuoD/HpQIIul2kJIXNx+MK29E5N8I5l+acS6tZs2aUyxIREZEymzvXL3t+9NF+Zcb27f1o9pAhCthSYcT7SPZy4NAC7w8BsgOqRURERMrj+++hZ0//AOPee0PnztChAxx4YNCViURcvI9kTwZuCs0y0hhYq35sERGRBPP113DppXDKKfDRR/Dvf8PSpb4PWwFbKqhAR7LN7FXgXKCGmS0HHgWqAjjnhgFTgEuABcBG4NZgKhUREZFS++QTP3L9wQdwwAHQqxfcey/su2/QlYlEXdCzi1xXwucOuDdG5YiIiEh5OedDdY8e8NlnUKsWZGRA27a+RUSkkoj3nmwRERFJBM7BO+/4ketvvoHUVBg4EO64A5KTg65OJObivSdbRERE4tn27TBxou+3vvxy+P13GDYMFi6E++5TwJZKSyFbRERESi8vD1591a/E2Lo1/PUXPP88/Pwz3HUX7LFH0BWKBEohW0RERMK3datfmfHoo+H6632byCuvwPz5cMstULVq0BWKxAX1ZIuIiEjJNm+G0aOhd29YsgROOglefx2uvBKqaMxOZGf6r0JERESKlpsLgwbBEUf4NpADD4S33vILy7RurYAtUgSNZIuIiMiu/vrLP8D4xBP+YcazzoJRo+Cii8As6OpE4p5CtoiIiPzPunUweDAMGACrV8MFF8C4cXDOOUFXJpJQFLJFREQE1qyBp5/2c1vn5MAll0D37nDGGUFXJpKQFLJFREQqs5Ur4cknYcgQWL8eWrb04fqUU4KuTCShKWSLiIhURitW+OXOhw2DTZvg6quhWzc4/vigKxOpEBSyRUREKpNly6BfP3j2Wdi2Df71L0hPh6OOCroykQpFIVtERKQyWLTIz3E9erR/f/PN0LUr1K8fbF0iFZRCtoiISEX200/w+ON+VcbddoM2baBzZ6hTJ+jKRCo0hWwREZGKaM4c6NULxo+H5GS4/37o2BFq147I4SfNzCJjaibZObnUTkmmU7OGtGyUGpFji1QECtkiIiIVyYwZ0LMnTJoEe+8NXbpAhw5Qs2bETjFpZhbpE+eQuzUPgKycXNInzgFQ0BYJ0VqoIiIiFcFXX8Gll0JaGnz8MTz6KCxd6vuwIxiwATKmZv4dsPPlbs0jY2pmRM8jksg0ki0iIpLIPvkEevSADz+EAw7wLSL33gv77hu1U2bn5JZqu0hlpJFsERGRROMcvP8+/OMfcO65MHcuPPEELFkCDz8c1YANUDsluVTbRSojhWwREZFE4Ry89RY0bgzNmsHixTBokP/ZsaPvwY6BTs0aklw1aYdtyVWT6NSsYUzOL5II1C4iIiIS77Zvhzfe8A80zpoFdevC8OF+rus99oh5OfkPN2p2EZGiKWSLiIjEq7w8GDfO91n/+CMceSS88AJcfz1UrRpoaS0bpSpUixRD7SIiIiLxZutWH6aPPtove27mF5P58Uc/eh1wwBaRkmkkW0REJF5s3uzDdZ8+/iHGRo1gwgRo2RKqaFxMJJHov1gREZGg5ebCwIFQvz60bQu1asHbb/uFZVq1UsAWSUAayRYREQnKX3/BsGF++r3ff4ezz4bnn4cLL/QtIiKSsBSyRUREYm3tWhg8GJ58Elav9qF6/Hg/77WIVAgK2SIiIrGyZg08/bRvDcnJ8cugd+sGZ5wRSDmTZmZpGj6RKFHIFhERibaVK2HAABgyxLeIXHkldO8OJ58cWEmTZmaRPnEOuVvzAMjKySV94hwABW2RCAj0SQoza25mmWa2wMy6FvL5uWa21sxmhX79O4g6RUREyiQ7Gx580C8e068fXHYZzJkDEycGGrDBLySTH7Dz5W7NI2NqZkAViVQsJY5km1l9YLlzbrOZnQucAIxxzuWU58RmlgQMAS4ClgPfmdlk59yPO+36mXPusvKcS0REJKaWLvWhetQo2LYNbrgB0tOhYfwsO56dk1uq7SJSOuGMZE8A8szsCGAUUA94JQLnPg1Y4Jxb5JzbAowFWkTguCIiIsFYuBDuuAOOOAJGjoSbboKff/ZzX8dRwAaonZJcqu0iUjrhhOztzrltwJXAU865B4GDI3DuVODXAu+Xh7bt7Awz+8HM3jWzY4s6mJm1MbPpZjb9jz/+iEB5IiIiYfrpJx+oGzaEl17yc10vWAAjRsDhhwddXaE6NWtIctWkHbYlV02iU7P4+suASKIK58HHrWZ2HXAzcHloWyTWcy1sAlC30/vvgcOcc3+Z2SXAJKBBYQdzzo0ARgCkpaXtfBwREZHImzMHevaE116D5GS4/3546CE4OBJjUeEryywh+Z9rdhGR6AgnZN8KtAV6OecWm1k94KUInHs5cGiB94cA2QV3cM6tK/B6ipk9Y2Y1nHOrInB+ERGRspkxA3r0gDffhOrVoWtX/4BjzZoxL6U8s4S0bJSqUC0SJSW2i4QeROyCH1XGObfYOdcnAuf+DmhgZvXMbHfgWmBywR3M7CAzv+SVmZ0Wqnd1BM4tIiJSel9+CZdcAmlp8Mkn8OijsGQJPP54IAEbNEuISLwqMWSb2eXALOC90PuTzGxysV8KQ6jPux0wFZgPjHfOzTOztmbWNrTbVcBcM/sBGAhc65xTK4iIiMSOc/Dxx3DBBdCkCXz3nQ/VS5fCY4/B/vsHWp5mCRGJT+G0izyGnwnkYwDn3KxQy0i5OeemAFN22jaswOvBwOBInEtERKRUnIP33/c9159/DgcdBP3781bjy+nz6a9kP/5ZWH3Mpe2XLu3+tVOSySokUGuWEJFghTO7yDbn3Nqdtmk0WUREKibn4K234PTToXlz3w4yaBAsWsSk866h87sLycrJxfG//udJM7MKPVR+v3S09gfNEiISr8IJ2XPN7HogycwamNkg4Mso1yUiIhJb27fD669Do0ZwxRWwapWfgm/hQmjXDpKTS93/HO39wT+82LvV8aSmJGNAakoyvVsdrwcaRQIWTrvIfUA3YDPwKr6Hukc0ixIREYmZbdtg3Djo1Qvmz4cjj/SLx1x/PVTdccba0vY/R3t7Ps0SIhJ/wpldZKNzrptz7lTnXFro9aZYFCciIhI1W7fC88/z1+EN4IYbyPxjA/++thuTXp4GN9+8S8CG0q+SGO3tIhK/igzZZvaWmU0u6lcsixQREYmYzZth2DBo0ABuu41lW5O468qHaX7bYMYcdgbpb/5YZA90afufo72/iMSv4tpFnohZFSIiIiUozawbhe7bcD8YORL69YPsbDj9dB46tw2v1zoB7H+LEOf3QBd27NKukhjt/UUkfllFnHY6LS3NTZ8+PegyREQkQnZe1RD8CG9hD/jtvO+eW3K5dfZ73DfzTaqtWQX/+Ac88ghccAH10qcUOl2WAYv7XBrF35GIVARmNsM5l1bYZ0WOZJvZeOfc1WY2h0Km7HPOnRDBGkVERIpU3KwbO4fs/H2rb97AzTPe4vbv3mS/Tev59ohTOO2NCT5kh2iOaRGJluLaRe4P/bwsFoWIiIgUpTSzbmxc8TsPTp/MrTPeYp/NG/ig/qkMOeMaZqUexeICARt8D3RhI+TqgRaR8ioyZDvnVoRe3uOc61LwMzPrC3TZ9VsiIiKRF9aI88qV0L8/XwwbxJ5bcnn3yDMZfOY1zKtVH/DzR+9MPdAiEi3hzJN9EbsG6osL2SYiIhIVxY44Z2XBE0/A8OGweTOrm17OtXWaM3u/Q3fdtxCaY1pEoqG4KfzuDvVjNzSz2QV+LQZmx65EERGp7Apb1fCpxvvRcmQvOPxwv+z51VfD/Pkc+u4kbmt7uVZAFJFAFTm7iJntC+wH9Aa6FvhovXNuTQxqKzPNLiIiUoEtWAC9e8OYMX7qvVtvhS5dfNgWEYmhMs0u4pxbC6wFrjOzJKBWaP+9zWxv59yyqFQrIiJSmPnz4fHH4ZVX/GqMbdtC585w6KElf1dEJMZK7Mk2s3bAY8DvwPbQZgdoCj8REYm+2bOhZ094/XVIToYHH4SOHeHgg4OuTESkSOE8+PgA0NA5tzrKtYiISCVT7CqO06f7cP3mm1C9OqSnwwMPQM2agdYsIhKOcEL2r/i2ERERkYjZeWXGrJxc0ifOYf8fpvOPccPgvfcgJQUeewzat4f99gu0XhGR0ggnZC8CPjazd4DN+RudcwOiVpWIiFR4O6zi6BxnLJvDfV+O5cxls6FGDf9w4z33wD77BFuoiEgZhBOyl4V+7R76JSIiUm7ZObngHOcs/p52X47j1KwfWbnXfvQ873a6v/U07LVX0CWKiJRZiSHbOfefWBQiIiKViHNcnT2T6z4Yw0krfiGrek0euagt409oSo0a+9JdAVtEElw4s4vUBDoDxwLV8rc7586PYl0iIlIRbd8OEyZAz570nT2bZSkH0aX5fUw87ny2JlUtdmVGEZFEEk67yMvAOOAyoC1wM/BHNIsSEZHEtvOsIZ0vqE+LzM/8PNfz50PDhjB6NDOPPofPP1zItpxcUneeXUREJIGFE7IPcM6NMrP7nXOfAJ+Y2SfRLkxERBJTwVlDquZt5cxP3+ekfq/Dn9lw3HEwdixcdRUkJdECaHHqYUGXLCISceGE7K2hnyvM7FIgGzgkeiWJiEgiy5iaSV7uJm6YM422X7/GIev+YE6t+qTf8B96j+4OVaoEXaKISNSFE7J7mtm+QEdgELAP8GBUqxIRkcS0cSNNPxjLXd9M4KC/1vB97YZ0b3oPHx+ehpnRWwFbRCqJcGYXeTv0ci1wXnTLERGRhLR+PQwdCv378+jKlXx96HF0uLQDXx52IpgBUDslOeAiRURiJ5zZRZ4H3M7bnXO3RaUiERFJHDk5MHgwPPkkrFkDTZvy2dV30WbJnv9baAY0a4iIVDrhtIu8XeB1NeBKfF+2iIhUVqtXw1NPwcCBsG4dXHYZdO8Op5/O2UDvnWYX0awhIlLZhNMuMqHgezN7FfggahWJiEj8+v136N8fnnkGNmyA1q2hWzdo1GiH3Vo2SlWoFpFKrSxPoDQA6kTi5GbW3MwyzWyBmXUt5HMzs4Ghz2eb2cmROK+IiJRSVhY88ADUq+dD9hVXwNy58PrruwRsEREJryd7Pb4n20I/fwO6lPfEZpYEDAEuApYD35nZZOfcjwV2uxgf6hsApwNDQz9FRCQWli6FPn3guecgLw9uvBHS0+HII4OuTEQkroXTLlI9Suc+DVjgnFsEYGZjgRZAwZDdAhjjnHPA12aWYmYHO+dWRKkmEREBWLDAr8744ot+dpDbboMuXfxItoiIlKjYkG1mycC/gGNCm6YDrzvntkTg3KnArwXeL2fXUerC9kkFdgnZZtYGaANQp05EullERCqf+fOhVy949VXYfXe4+27o3BkO0RpkIiKlUWRPtpkdD8wHzgaWAEuBZsAXoRHlnuU8txWybeepAsPZx290boRzLs05l1azZs1yliYiUsn88AP8859w7LEwaRJ06ACLF/vZQxSwRURKrbiR7IHAnc65aQU3mtmFwFxgXjnPvRw4tMD7Q9h1asBw9hERkbL67jvo2RMmT4bq1X2/9YMPQo0aQVcmIpLQiptd5OCdAzaAc+4DYCt+vuzy+A5oYGb1zGx34Fpg8k77TAZuCs0y0hhYq35sEZEI+OILaN4cTjsNPvsM/vMf/5Bjr14K2CIiEVDcSHYVM9vDObe54EYzqwZsdc5tLM+JnXPbzKwdMBVIAp5zzs0zs7ahz4cBU4BLgAXARuDW8pxTRKRScw7++1/o0QM+/hhq1vQzh9x9N+yzT9DViYhUKMWF7DHABDNr55xbAmBmdfFtJC9G4uTOuSn4IF1w27ACrx1wbyTOJSJSaTkH773n20K+/BIOPhgGDIA2bWCvvYKuTkSkQioyZDvneoZGmj81sz1DmzcATzjnBsWkOhERKbvt2+Gtt3y4nj4dDj0Uhgzx0/FVq/b3bpPCXAI93P1ERKSEKfycc4OBwWZWPfR+fUyqEhGRssvLgwkTfH/17Nlw+OEwciTcdJOflq+ASTOzSJ84h9yteQBk5eSSPnEOwA4BOtz9RETEC2tZdefcegVsEZE4t20bvPQSHHccXHMNbN4MY8ZAZibccccuARsgY2rm38E5X+7WPDKmZpZpPxER8cIK2SIiEse2bIFRo+Coo/yy51Wr8m2foZx90yDqzdufJk98yqSZWYV+NTsnN6zt4e4nIiJeicuqi4hIfNi5J7rLuYdxxcz3/Qwhy5bBKafAG28w6dBTSJ80j9ytfnHe4lo7aqckk1VIUK6dklym/URExCtxJNvM9jSzR8xsZOh9AzO7LPqliYhIvvye6KycXPbYuolm08ZyetPT4Z57IDUVpkzxC8u0bEnGtF/Cbu3o1KwhyVWTdtiWXDWJTs0almk/ERHxwhnJfh6YAZwRer8ceA14O1pFiYjIjjKmZlLlr/W0nTmFO757gxob1/JVnePpeU06g0Z0ALO/9y1Na0f+yHZJs4aEu5+IiHjhhOz6zrlrzOw6AOdcrlmB/zcXEZHoysmh9TvPcev0yey3aT2f1m3EwCbXMv2QYzFg0E7/l1za1o6WjVLDCsvh7iciIuGF7C1mlgw4ADOrD2wu/isiIlJuq1bBU0/BoEF0WLeOaUecxuAzruGH2v9r0SgsOHdq1nCH6fZArR0iIrEWTsh+FHgPONTMXgaaALdEsygRkUrtt9+gf38YOhQ2bIDWrfnvVW1oP297WMFZrR0iIsEzv3J5CTuZHQA0Bgz42jm3KtqFlUdaWpqbPn160GWIiJROVhb06wcjRvhp+a69Fh5+GI49FtCKiyIi8cbMZjjn0gr7rMiRbDM7eadNK0I/65hZHefc95EqUESkUluyxE/D9/zzfin0G2+E9HRo0GCH3dQTLSKSOIprF+lfzGcOOD/CtYiIVC6//AK9e7P9xRfZ5ozXjr+AiU1v5MZrz6FlA4VpEZFEVmTIds6dF8tCREQqkmJbO378EXr1grFjyatalZcbXcozaVfy2z41/MdFLBwjIiKJo8QHH82sGnAPcBZ+BPszYJhzblOUaxMRSUj5C8fkP6SYv+LivpnzOO/1ETBxIuy5J3TsSMuqpzEnb8cZQvIXjlHIFhFJXOHMLjIGWA8MCr2/DngR+Ge0ihIRSWQZUzN3mAXkxOxM2n01jvMWfAv77OMfZnzgAahRg7ld3yn0GEUtKCMiIokhnJDd0Dl3YoH3/zWzH6JVkIhIossPyGnL59H+i7H8Y8lMcqrtzYCz/kWHtwZDSsrf+5Z24RgREUkM4YTsmWbW2Dn3NYCZnQ58Ed2yRESir7RT4oW1v3Ncvmo+178/msa/zmXVnvvS55xbeLHRJaTUOoAOBQI2aOEYEZGKKpyQfTpwk5ktC72vA8w3szmAc86dELXqRESipKi+aSj8gcMS93cO3nsPevRg4FdfsXLv/fm/8+/klZOasalqNS0cIyJSyZS4GI2ZHVbc5865pRGtKAK0GI2IlKRJn48KbdNITUnmi667zlBa1P6H7LMHnx/zF/TsCTNmQJ060LUrkxs1pe/HSxWcRUQqsDItRpPPObfUzPYDDi24vxajEZFEVtSDheFur7I9j4szv+S+r8bBH0ugfn149lm/kMzuu3MFcEXj+hGuWkREEkU4U/j1AG4BFuKn8AMtRiMiCa60Dxzm75+0PY8rfvyEe78azxFrlrOkZh148UW/BPpu4XTgiYhIZRDOnwhXA/Wdc1uiXYyISHmU5kHG0j5w2Pn8ekz/v6e544txHJbzG/Nr1uWBVg9zbvpd1E2rE5Xfj4iIJK5wQvZcIAVYGd1SRETKrrQPMob9wOGmTfDcc7To25cWy5bxU+qRtDn/Dn485Rweuvho9VmLiEihwnnwMQ14Ex+2N+dvd85dEd3Syk4PPopUPqV9kLFEGzfC8OGQkQErVsCZZ8Ijj0CzZmAWgYpFRCTRlevBR2A00BeYA2yPZGEiIpFS2gcZi7R+PQwZAgMGwB9/wLnnwksvwXnnKVyLiEjYwgnZq5xzA6NeiYhIOZR75cScHBg4EJ56Cv78049Yd+8OZ50V0TpFRKRyqBLGPjPMrLeZnWFmJ+f/inplIiKl0KlZQ5KrJu2wLayVE1etgm7d4LDD4NFH4eyz4dtv/cIyCtgiIlJG4YxkNwr9bFxgW7mm8DOz/YFxQF1gCXC1c+7PQvZbAqwH8oBtRfW8iIiUeuXE336D/v1h6FDff926tR+5PvHEGFYtIiIVVYkPPkblpGb9gDXOuT5m1hXYzznXpZD9lgBpzrlVpTm+HnwUkSItXw79+sHIkbBlC1x3HTz8MBxzTNCViYhIginvg4+Y2aXAsUC1/G3Ouf8rR00tgHNDr0cDHwO7hGwRkYhZvBj69IHnnwfn4KaboGtXaNAg6MpERKQCCmfFx2HAnsB5wLPAVcC35TxvLefcCgDn3AozO7CI/Rzwvpk5YLhzbkQ5zysilc3PP0Pv3n5VxqQkuP126NIF6tYNujIREanAwhnJPtM5d4KZzXbO/cfM+gMTS/qSmX0AHFTIR91KUV8T51x2KIRPM7OfnHOfFnG+NkAbgDp1tPqaSKU3bx706gXjxsHuu0O7dtCpE6Rq8RgREYm+cEJ2/pxYG82sNrAaqFfSl5xzFxb1mZn9bmYHh0axD6aI1SSdc9mhnyvN7A3gNKDQkB0a5R4Bvie7pPpEpIKaORN69oSJE2GvvaBjR/+rVq2gKxMRkUoknCn83jazFCAD+B4/G8ir5TzvZODm0Oub8StK7sDM9jKz6vmvgab4VSdFRHb17bdw+eVw8snwwQd+ppClS/1DjgrYIiISYyWOZDvneoReTjCzt4Fqzrm15TxvH2C8md0OLAP+CRAaKX/WOXcJUAt4w/wKa7sBrzjn3ivneUWkovnsM+jRA6ZNg/3396/btYOUlKArExGRSqzIkG1mpwK/Oud+C72/CWgNLDWzx5xza8p6UufcauCCQrZnA5eEXi8CNGGtiOzKOfjoIx+oP/kEDjwQ+vaFu++G6tWDrk5ERKTYdpHhwBYAM/sHfvR5DLCWUO+ziEhMOQdTpsCZZ8KFF8Ivv/hl0Bcvhs6dFbBFRCRuFNcuklRgtPoaYIRzbgK+bWRW1CsTEcm3fTu8+aZ/oPH77/0S6EOHwi23QLVqJX5dREQk1oobyU4ys/wQfgHwUYHPwlrERkSkXPLy/BR8J54IrVrB2rUwapQfwW7bVgFbRETiVnFh+VXgEzNbhZ/G7zMAMzsC3zIiIhId27bBK6/A449DZiYcfTS89BJccw3spr/ji4hI/CvyTyvnXC8z+xA4GHjfOZc/93QV4L5YFCcilcyWLTB6tF/+fNEiOOEEGD8eWreGKuHMOCoiIhIfih0Scs59Xci2n6NXjohUSps2+TaQvn3h118hLQ2efBIuu0zhWkREEpL+3VVEgrNhAwwfDhkZ8Ntv0KQJjBwJTZuCnyNfREQkISlki0jsrVsHzzwD/fvDqlVw3nm+B/vccxWuRUSkQlDIFpHY+fNPGDgQnn7av27e3C9/3qRJ0JWJiIhElEK2iETfH3/4HuvBg2H9emjRArp1g1NPLfZrk2ZmkTE1k+ycXGqnJNOpWUNaNkqNUdEiIiJlp5AtItGzYoVvCRk6FHJz4aqrfLg+8cQSvzppZhbpE+eQuzUPgKycXNInzgFQ0BYRkbinx/ZFJPJ+/RXuuw/q1fMj2K1awbx5fjq+MAI2QMbUzL8Ddr7crXlkTM2MRsUiIiIRpZFsEYmcxYuhd2944QVwDm66CdLT4YgjSn2o7JzcUm0XERGJJwrZIlJ+P//sV2d86SVISoI77oAuXeCww8p8yNopyWQVEqhrpySXp1IREZGYULuIiJTd3Llw3XV+2fPx432LyKJFfnq+cgRsgE7NGpJcNWmHbclVk+jUrGG5jisiIhILGskWkdKbORN69oSJE2HvveGhh6BDB6hVK2KnyH+4UbOLiIhIIlLIFpHwffMN9OgB77wD++4LjzwC998PBxxQ4lfLMh1fy0apCtUiIpKQ1C4iIiX79FO46CJo3Jgtn3/B8Ka3ccItw2my57lMWrapxK/nT8eXlZOL43/T8U2amRX92kVERAKgkWwRKZxz8OGHfuT600/hwAOZe383bt7jFFbb7gCsC3Pu6uKm49NItYiIVEQayRaRHTnn20HOOMOPXi9c6JdBX7yYuw46/++AnS+cuas1HZ+IiFQ2Ctki4m3f7h9kPOUUuOwy+O03v1LjwoXQvj3suWeZw3JR0+5pOj4REamoFLJFKru8PBg71q/E2Lo1rF8Pzz0Hv/wCbdvCHnv8vWtZw7Km4xMRkcpGIVukstq6FUaPhmOO8XNd5+X5xWTmz4dbb4WqVXf5SlnDcstGqfRudTypKckYkJqSTO9Wx6sfW0REKiw9+ChS2Wze7MN1nz5+GfQTT4TXXoNWraBK8X/vLs/c1ZqOT0REKhOFbJHKIjcXRo2Cvn1h+XI49VT/QONll4FZ2IdRWBYRESmZQrZIRbdhAwwbBk884R9mbNIEnn0WmjYtVbgWERGR8Clki1RU69bBkCEwYACsWgXnnw+vvgrnnKNwLSIiEmUK2SIVzZ9/+jaQp5+GnBy4+GLo3h3OPDPoykRERCoNhWyRiuKPP/yo9ZAhfhq+li2hWzdISwu6MhERkUonkCn8zOyfZjbPzLabWZEJwMyam1mmmS0ws66xrFEkYaxYAR07Qt26/qHGiy+GH36AN94oMmBPmplFkz4fUa/rOzTp8xGTZmbFtmYREZEKLqiR7LlAK2B4UTuYWRIwBLgIWA58Z2aTnXM/xqZEkTj3668+VD/7LGzbBtdfD+npcPTRxX5t0sws0ifOIXdrHgBZObmkT5wDoFlDREREIiSQkWzn3HznXGYJu50GLHDOLXLObQHGAi2iX51InFu0CNq0gfr1YfhwuOEGyMyEMWNKDNjg57jOD9j5crfmkTG1pP8kRUREJFzx3JOdCvxa4P1y4PSAahEJXmYmPP44vPwyJCXBnXdCly5Qp06pDpOdk1uq7SIiIlJ6UQvZZvYBcFAhH3Vzzr0ZziEK2eaKOV8boA1AnVKGDpG4Nncu9OwJ48dDtWrQvj089BDUrl2mw9VOSSarkEBdOyW5vJWKiIhISNRCtnPuwnIeYjlwaIH3hwDZxZxvBDACIC0trcgwLpIwvv/eh+s33oC994bOnaFDBzjwwHIdtlOzhjv0ZAMkV02iU7OG5a1YREREQuK5XeQ7oIGZ1QOygGuB64MtSSQGvv4aevSAKVNg333h3//2o9cHHBCRw+c/3JgxNZPsnFxqpyTTqVlDPfQoIiISQYGEbDO7EhgE1ATeMbNZzrlmZlYbeNY5d4lzbpuZtQOmAknAc865eUHUKxITn3ziR64/+MAH6p49oV07H7R3MmlmVrlCcstGqQrVIiIiUWTOVbzOirS0NDd9+vSgyxApmXM+VPfoAZ99BrVq+X7rtm19i0ghdp6CD3y7R+9Wxys4i4iIxJCZzXDOFbooRSBT+IlUes7B22/DGWdA06Z+Wr6nn4bFi33ILiJgg6bgExERSQQK2SKxtH07TJgAJ58Ml18Ov/8Ow4bBwoW+7zq55Bk+NAWfiIhI/FPIFomFvDx49VU44QS46irYsAGefx5+/hnuugv22CPsQxU11Z6m4BMREYkfCtki0bR1K7zwgl+J8frrfZvIyy/D/Plwyy1QtWqpD9mpWUOSqybtsE1T8ImIiMSXeJ7CTyRxbd4Mo0dD796wZAmcdBK8/jpceSVUKd/fbTUFn4iISPxTyBaJpNxcePZZ6NcPli+H006DQYPg0kvBClvEtGw0BZ+IiEh8U8gWiYS//vIPMD7xhH+Y8ayzYNQouOiiiIZrERERSQwK2SLlsW4dDB4MAwbA6tVwwQUwbhycc07QlYmIiEiAFLJFymLNGj+v9cCBkJMDl1wC3bv7ea9FRESk0lPIFimNlSvhySdhyBBYvx5atvTh+pRTgq5MRERE4ohCtkg4VqyAjAzfd71pE1x9NTz8sJ/3WkRERGQnCtkixVm2DPr29Q8xbtsG//oXpKfDUUcFXZmIiIjEMYVskcIsWuTnuB492r+/+Wbo2hXq1w+2LhEREUkICtkiBf30Ezz+OLzyCuy2G7RpA507Q506QVcmIiIiCUQhWwRgzhzo1QvGj4dq1aB9e3joIahdO+jKREREJAEpZEvlNmMG9OwJkybB3ntDly7w4INw4IFBVyYiIiIJTCFbKqevvvLhesoU2Hdf+Pe/4f77Yf/9g65MREREKgCFbKk8nINPPvHh+sMP4YADfIvIvff6oB1lk2ZmkTE1k+ycXGqnJNOpWUNaNkqN+nlFREQk9hSypeJzDqZNgx494PPPoVYteOIJuOsu3yJShEiG4kkzs0ifOIfcrXkAZOXkkj5xDoCCtoiISAVUJegCRKLGOXjrLWjcGJo1g8WL/TLoixdDx44lBuz0iXPIysnF8b9QPGlmVplKyZia+XfAzpe7NY+MqZllOp6IiIjEN4VsqXi2b4cJE+Dkk+GKK/xS6MOHw8KFcN99kJxc4iEiHYqzc3JLtV1EREQSm0K2VBx5eX5+6+OPh6uugg0b4IUX4Oef/XzXe+wR9qEiHYprpxQe7IvaLiIiIolNPdmScHbule58/uG0+PFjv4jML7/Ascf6sH311ZCUVKZz1E5JJquQQF3WUNypWcMderIBkqsm0alZwzIdT0REROKbRrIloRTsla66bSvnfDyRU5o2hltvhb328m0is2fDddeVOWCDD8XJVXf8fnlCcctGqfRudTypKckYkJqSTO9Wx+uhRxERkQpKI9mSUDKmZrJ940Zu+WEqd30zgYP/Ws2sg4/k6Rb3kfHCw2AWkfPkh99ITrnXslGqQrWIiEgloZAtieOvv7h06svc+d1Eam7I4ZtDjqXTJQ/wed2TMDMyIhSw8ykUi4iISFkpZEv8W7sWBg+GJ5/k4dWr+eywk2h3xTV8U+f4v3fRA4QiIiISTxSyJX6tWQNPPeXntl67Fi69lE+uakPbhbvrAUIRERGJa3rwUeLPypXQtSscdphfpfH882HGDHj7bc655Qo9QCgiIiJxL5CRbDP7J/AYcDRwmnNuehH7LQHWA3nANudcWqxqlABkZ0NGhl84ZtMmuOYaePhhP+91AeqVFhERkXgXVLvIXKAVMDyMfc9zzq2Kcj0SpKVLoV8/GDUKtm2Df/3Lh+uGagERERGRxBRIyHbOzQewCM8GIQlm4ULo3RtGj/ZT791yi28TOfzwoCsTERERKZd4f/DRAe+bmQOGO+dGBF2QRMBPP/nVGV95BXbbDe66Czp3hjp1gq5MREREJCKiFrLN7APgoEI+6uacezPMwzRxzmWb2YHANDP7yTn3aRHnawO0AaijsBafZs+GXr3gtdcgORnuvx8eeggOPjjoykREREQiKmoh2zl3YQSOkR36udLM3gBOAwoN2aFR7hEAaWlprrznlgiaMcPPEvLmm1C9um8JefBBqFkz6MpEREREoiJu20XMbC+ginNufeh1U+D/Ai5LgEkzs8JbbvzLL6FnT3j3XUhJgUcfhfbtYf/9Y16ziIiISCwFMk+2mV1pZsuBM4B3zGxqaHttM5sS2q0W8LmZ/QB8C7zjnHsviHrlfybNzCJ94hyycnJxQFZOLukT5zBpZpbfwTn4+GO44AJo0gS++873Xy9dCo89poAtIiIilYI5V/E6K9LS0tz06YVOvS3l1KTPR2Tl5O6yPXXfanxx8jY/cv3551CrFnTqBG3bwl57BVCpiIiISHSZ2Yyi1nGJ23YRiU/ZOwds57hg4be0/3IsrPgFDjkEBg2C22/3DzeKiIiIVEIK2VIqtVOSycrJxdx2mmd+yX1fjeOYlYvJ2u9gGDECbroJ9tgj6DJFREREAqWQLaXS+YL6fN5rCG0+H0uD1b+ycP9Uul7egTO630fqaXWDLk9EREQkLihkS3i2boWXXqLF44/TYsECFtWqy32Xd2JW44voePExtChsdhERERGRSkohW4q3eTM8/zz06eNnCGnUCCZM4PCWLRlUJZDJaURERETinkK2FG7jRhg5Evr1g+xsOP10GDIELrkEzIKuTkRERCSuKWTLjv76C4YOhSeegJUr4R//gNGj/bzXCtciIiIiYVHIFm/tWj/13pNPwpo1cNFF0L27D9kiIiIiUioK2ZXd6tXw9NMwcKAP2pde6sN148ZBVyYiIiKSsBSyK6uVK6F/f3jmGd8i0qqVD9eNGgVdmYiIiEjCU8iubLKyfL/18OGwaRNccw106wbHHRd0ZSIiIiIVhkJ2ZbF0KfTtC6NGQV4e3HADpKdDw4ZBVyYiIiJS4ShkV3QLFkDv3jBmjJ8d5JZboGtXOPzwoCsTERERqbAUsiuq+fPh8cfhlVegalVo2xY6d4ZDDw26MhEREZEKTyG7opk9G3r2hNdfh+RkePBB6NgRDj446MpEREREKg2F7Ipi+nQfrt98E6pX9y0hDz4INWsGXZmIiIhIpaOQnei+/BJ69ID33oOUFHjsMWjfHvbbL+jKRERERCothexE5Bx8/LEP1//9L9So4R9uvOce2GefoKsTERERqfQUshOJczB1qm8L+eILOOggv6DMXXfBXnsFXZ2IiIiIhChkJwLn4K23/Mj19Ol+hpDBg+H226FataCrExEREZGdVAm6ACnG9u3w2mtw0knQogWsXg0jRvi5r++9VwFbREREJE4pZMejbdvgpZf8UudXX+2XPx89Gn7+Ge68E3bfPegKRURERKQYCtnxZMsWeO45OPpouPFGSEqCsWPhxx/hpptgN3X3iIiIiCQCpbZ4sHmzD9d9+sCyZXDyyTBxom8RqaK/B4mIiIgkGoXsIG3cCCNHQr9+kJ0NjRvD0KFw8cVgFnR1IiIiIlJGCtlBWL/eh+n+/WHlSjjnHBgzBs4/X+FaREREpAJQyI6lnBwYNAieegrWrIGLLoJHHoGzzw66MhERERGJIIXsWFi92gfrgQNh3Tq47DLo3h1OPz3oykREREQkChSyo+n3331LyDPPwIYN0KqVD9eNGgVdmYiIiIhEUSBTV5hZhpn9ZGazzewNM0spYr/mZpZpZgvMrGuMyyy7rCx44AGoV8+H7CuugLlzmdR9IE2m/km9ru/QpM9HTJqZFXSlIiIiIhIFQc0PNw04zjl3AvAzkL7zDmaWBAwBLgaOAa4zs2NiWmVpLV0Kd98Nhx/ulz2/5hqYPx9eeYVJW1JInziHrJxcHJCVk0v6xDkK2iIiIiIVUCAh2zn3vnNuW+jt18Ahhex2GrDAObfIObcFGAu0iFWNpbJgAdx2GxxxBIwaBbfcAr/8As8/D0ceCUDG1Exyt+bt8LXcrXlkTM0MoGARERERiaZ46Mm+DRhXyPZU4NcC75cD8fWk4Pz50KsXvPqqX+r87ruhc2c4ZNe/M2Tn5BZ6iKK2i4iIiEjiilrINrMPgIMK+aibc+7N0D7dgG3Ay4UdopBtrpjztQHaANSpU6fU9ZbKDz9Az54wYQIkJ0OHDtCxIxxU2G/Xq52STFYhgbp2SnI0KxURERGRAEQtZDvnLizuczO7GbgMuMA5V1h4Xg4cWuD9IUB2MecbAYwASEtLKzKMl8t33/lwPXkyVK8O6enw4INQo0aJX+3UrCHpE+fs0DKSXDWJTs0aRqVUEREREQlOIO0iZtYc6AKc45zbWMRu3wENzKwekAVcC1wfoxJ39MUX0KMHTJ0K++0H//kP3Heffx2mlo1SAd+bnZ2TS+2UZDo1a/j3dhERERGpOILqyR4M7AFMM7+M+NfOubZmVht41jl3iXNum5m1A6YCScBzzrl5MavQOfjvf324/vhjP1rduzfccw/ss0+ZDtmyUapCtYiIiEglEEjIds4dUcT2bOCSAu+nAFNiVVfopPDee74t5Msv4eCDYcAAaNMG9torpqWIiIiISGKKh9lF4sP27fDWWz5cT58Ohx4KQ4b4qfmqVQu6OhERERFJIEEtRhM/8vJg/Hi/1HnLlrBmDYwc6ee+vuceBWwRERERKbXKG7K3bYOXXoLjjvMrM27eDGPGQGYm3HGHn/daRERERKQMKl/I3rLFr8p41FFw441QtSqMGwfz5vn3u6mDRkRERETKp/Ikyk2b4LnnoG9fWLYMTj4Z3ngDrrgCqlS+v2uIiIiISPRU/JC9cSOMGAEZGZCdDWecAcOGQfPmYIUtKikiIiIiUj4VN2SvXw/PPAP9+8Mff8A55/ie6/PPV7gWERERkaiqmCF7xQqoW9fPFNK0KXTvDmefHXRVIiIiIlJJVMxm5OxsaNIEvvnGL4WugC0iIiIiMWTOuaBriDgz+wNYGuPT1gBWxficFYGuW+npmpWNrlvZ6LqVja5b2ei6lZ6uWdlE6rod5pyrWdgHFTJkB8HMpjvn0oKuI9HoupWerlnZ6LqVja5b2ei6lY2uW+npmpVNLK5bxWwXEREREREJkEK2iIiIiEiEKWRHzoigC0hQum6lp2tWNrpuZaPrVja6bmWj61Z6umZlE/Xrpp5sEREREZEI00i2iIiIiEiEKWSXkZllmNlPZjbbzN4ws5Qi9mtuZplmtsDMusa4zLhiZv80s3lmtt3Minyi18yWmNkcM5tlZtNjWWM8KsV1071WgJntb2bTzOyX0M/9ithP9xsl3z/mDQx9PtvMTg6izngSxjU718zWhu6tWWb27yDqjDdm9pyZrTSzuUV8rnttJ2FcM91rhTCzQ83sv2Y2P/Tn6P2F7BO1+00hu+ymAcc5504AfgbSd97BzJKAIcDFwDHAdWZ2TEyrjC9zgVbAp2Hse55z7iRNSwSEcd10rxWqK/Chc64B8GHofVEq9f0W5v1zMdAg9KsNMDSmRcaZUvw391no3jrJOfd/MS0yfr0ANC/mc91ru3qB4q8Z6F4rzDago3PuaKAxcG8s/79NIbuMnHPvO+e2hd5+DRxSyG6nAQucc4ucc1uAsUCLWNUYb5xz851zmUHXkWjCvG6613bVAhgdej0aaBlcKXEvnPunBTDGeV8DKWZ2cKwLjSP6b66MnHOfAmuK2UX32k7CuGZSCOfcCufc96HX64H5QOpOu0XtflPIjozbgHcL2Z4K/Frg/XJ2/R9XduWA981shpm1CbqYBKF7bVe1nHMrwP8fLXBgEfvpfgvv/tE9tqNwr8cZZvaDmb1rZsfGprSEp3utbHSvFcPM6gKNgG92+ihq99tukThIRWVmHwAHFfJRN+fcm6F9uuH/OeLlwg5RyLYKPZ1LONcsDE2cc9lmdiAwzcx+Cv0tvsKKwHWrdPcaFH/dSnGYSne/FSKc+6dS3mPFCOd6fI9fcvkvM7sEmIT/J2kpnu610tO9Vgwz2xuYADzgnFu388eFfCUi95tCdjGccxcW97mZ3QxcBlzgCp8LcTlwaIH3hwDZkasw/pR0zcI8Rnbo50ozewP/z7IVOvRE4LpVunsNir9uZva7mR3snFsR+qe/lUUco9Ldb4UI5/6plPdYMUq8HgX/MHfOTTGzZ8yshnNuVYxqTFS610pJ91rRzKwqPmC/7JybWMguUbvf1C5SRmbWHOgCXOGc21jEbt8BDcysnpntDlwLTI5VjYnIzPYys+r5r4Gm+Af/pHi613Y1Gbg59PpmYJd/EdD99rdw7p/JwE2hJ/EbA2vz23EqqRKvmZkdZGYWen0a/s/c1TGvNPHoXisl3WuFC12TUcB859yAInaL2v2mkeyyGwzsgf/nZYCvnXNtzaw28Kxz7hLn3DYzawdMBZKA55xz84IrOVhmdiUwCKgJvGNms5xzzQpeM6AW8Ebomu4GvOKcey+wouNAONdN91qh+gDjzex2YBnwTwDdb7sq6v4xs7ahz4cBU4BLgAXARuDWoOqNB2Fes6uAu81sG5ALXFvEv3pWKmb2KnAuUMPMlgOPAlVB91pRwrhmutcK1wS4EZhjZrNC2x4G6kD07zet+CgiIiIiEmFqFxERERERiTCFbBERERGRCFPIFhERERGJMIVsEREREZEIU8gWEREREYkwhWwRkVIyswPMbFbo129mlhV6nWNmP8a4lpZmdkyB9/9nZqVe3MjM6ppZoXOEm9mxZvaRmf1sZgvN7D9mFvE/P4r7vZjZx2aWFulziohEi0K2iEgpOedWO+dOcs6dBAwDngy9PgnYHunzmVlxaxq0BP4Ops65fzvnPojguZPxizX0cc4dCRyPXxXz/kido4CWRPH3IiISSwrZIiKRlWRmI81snpm9HwqpmFl9M3vPzGaY2WdmdlRo+2Fm9qGZzQ79rBPa/oKZDTCz/wJ9C/u+mZ0JXAFkhEbS64e+d1XoGKea2Zdm9oOZfWtm1UMj1p+Z2fehX2eW8Pu5HvjCOfc+QGiF23ZAp9A5HjOzh/J3NrO5ZlY39HpSqN55ZtamwD5/mVmvUF1fm1mtkn4vBZlZUzP7KlT/a2a2d2h7HzP7MXQtnyj9/3QiIpGjkC0iElkNgCHOuWOBHKB1aPsI4D7n3CnAQ8Azoe2DgTHOuROAl4GBBY51JHChc65jYd93zn2JH2XuFBpZX5j/RfNLfY8D7nfOnQhciF8JbiVwkXPuZOCanc5XmGOBGQU3hM6TbGYpJXz3tlC9aUB7MzsgtH0v/Cq5JwKfAncW93spyMxqAN1D1+VkYDrQwcz2B64Ejg1dy54l1CYiElVaVl1EJLIWO+dmhV7PAOqGRlrPBF4LLeEOsEfo5xlAq9DrF4F+BY71mnMur4TvF6UhsMI59x2Ac24dgJntBQw2s5OAPHyQL44BhS0NbIVs21l7M7sy9PpQ/F9AVgNbgLdD22cAF4VxrHyN8S0lX4Suxe7AV8A6YBPwrJm9U+D4IiKBUMgWEYmszQVe5wHJ+H81zAn1bZekYKDdEPpZmu/nKyocPwj8DpwYOu6mEo4zD/jHDgc2OxxY5ZzLMbNt7PivotVC+5yLHz0/wzm30cw+zv8M2Oqcy68tj9L9WWTANOfcdbt8YHYacAFwLb6l5fxSHFdEJKLULiIiEmWhUeTFZvZPAPNODH38JT4UAvwL+LyU318PVC/ktD8Btc3s1NB3qoceoNwXP8K9HbgRSCqh/JeBswrM8pGMbzF5NPT5EuDk0GcnA/VC2/cF/gwF7KPwI9AlKer3UtDXQBMzOyJ0zj3N7MjQaP++zrkpwAP4h1BFRAKjkC0iEhv/Am43sx/wo8MtQtvbA7ea2Wx86C1q1o6ivj8W6GRmM82sfv7Ozrkt+J7rQaHvTMOPJD8D3GxmX+NbRTZQDOdcLv6BxG5m9jOwCv8g5MuhXSYA+5vZLOBu4OfQ9veA3UK/rx74cFySQn8vO9XzB3AL8Gro2F8DR+HD+duhbZ/gR+xFRAJj//sXOxERkeKZWUtgAHCec25pwOWIiMQthWwRERERkQhTu4iIiIiISIQpZIuIiIiIRJhCtoiIiIhIhClki4iIiIhEmEK2iIiIiEiEKWSLiIiIiESYQraIiIiISIT9PwPOvCxV98eJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Import clathrate formation data\n",
    "dataset = pd.read_csv(\"https://raw.githubusercontent.com/forestwhite/Regression/main/data/data-table-B8.csv\")\n",
    "\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "# Create Xs and y;\n",
    "X = dataset[['x1','x2']]\n",
    "y1 = dataset['y']\n",
    "\n",
    "# degree = 2 indicates the second order interactions\n",
    "interaction = PolynomialFeatures(degree=2, interaction_only=True, include_bias=False)\n",
    "X_inter = interaction.fit_transform(X)\n",
    "\n",
    "# create a dictionary that maps X using the itertools.count() function, simlar to an enumeration\n",
    "from itertools import count\n",
    "x_to_feature = dict(zip(('x{}'.format(i) for i in count()), X.columns))\n",
    "\n",
    "# substitute x# with column names or combinations of column names from the dataset,\n",
    "# replaces spaces with underscores to ease column references\n",
    "\n",
    "import seaborn as sns\n",
    "# necessary so `x11` enters the loop before `x1`, preventing duplicate x1 match\n",
    "feature_keys = list(x_to_feature.keys())[::-1]\n",
    "features = []\n",
    "\n",
    "for feature in interaction.get_feature_names_out():\n",
    "    feature = re.sub(' ', '_', feature)\n",
    "    features.append(feature)\n",
    "X_all_inter = pd.DataFrame(X_inter, columns=features)\n",
    "X_all_inter_const = sm.add_constant(X_all_inter)\n",
    "lm312 = sm.OLS(y1, X_all_inter_const).fit()\n",
    "\n",
    "# credit to Jeff Macaluso for demonstrating model adequecy tests\n",
    "# see https://jeffmacaluso.github.io/post/LinearRegressionAssumptions/\n",
    "\n",
    "# Calculating residuals for the plot\n",
    "predictions = lm312.predict()\n",
    "df_results = pd.DataFrame({'Actual': y1, 'Predicted': predictions})\n",
    "df_results['Residuals'] = abs(df_results['Actual']) - abs(df_results['Predicted'])\n",
    "\n",
    "# Plotting the residuals\n",
    "plt.subplots(figsize=(12, 6))\n",
    "ax = plt.subplot(111)  # To remove spines\n",
    "plt.scatter(x=df_results.index, y=df_results.Residuals, alpha=0.5)\n",
    "plt.plot(np.repeat(0, df_results.index.max()), color='darkorange', linestyle='--')\n",
    "ax.spines['right'].set_visible(False)  # Removing the right spine\n",
    "ax.spines['top'].set_visible(False)  # Removing the top spine\n",
    "plt.title('Residuals vs. Fitted')\n",
    "plt.show()\n",
    "\n",
    "# create Q-Q plot to test normality using statsmodels\n",
    "fig, ax = plt.subplots(figsize=(12, 6))\n",
    "fig = sm.qqplot(dataset['y'], fit=True, line='45', ax = ax)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "028faa54-4f4a-45ac-abcc-5ecdfba2208c",
   "metadata": {
    "tags": []
   },
   "source": [
    "**c. Use t tests to assess the contribution of each regressor to the model. Discuss your findings.**\n",
    "\n",
    "A very low p-value (2.318e-13) and relatively high R-squared coefficient (0.8518) support that the model is overall significant and a reasonably good fit for actual population from which this sample is taken. For each individual predictor variable, we have:\n",
    "\n",
    " * (Intercept): t = 6.602 -> Pr(>|t|) = 0.000 < 0.05, therefore significant\n",
    " * x1: t = 3.482 -> Pr(>|t|) = 0.001 < 0.05, therefore significant\n",
    " * x2: t = 8.281 -> Pr(>|t|) = 0.000 < 0.05, therefore significant\n",
    " * x1:x2: t = 1.492 -> Pr(>|t|) = 0.146 > 0.05, therefore not significant "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "b2caab6c-d127-404a-85bd-e96c34287949",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "' 0.000'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "' 0.001'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "' 0.000'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "' 0.146'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "\n",
    "# Import clathrate formation data\n",
    "dataset = pd.read_csv(\"https://raw.githubusercontent.com/forestwhite/Regression/main/data/data-table-B8.csv\")\n",
    "\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "# Create Xs and y;\n",
    "X = dataset[['x1','x2']]\n",
    "y1 = dataset['y']\n",
    "\n",
    "# degree = 2 indicates the second order interactions\n",
    "interaction = PolynomialFeatures(degree=2, interaction_only=True, include_bias=False)\n",
    "X_inter = interaction.fit_transform(X)\n",
    "\n",
    "# create a dictionary that maps X using the itertools.count() function, simlar to an enumeration\n",
    "from itertools import count\n",
    "x_to_feature = dict(zip(('x{}'.format(i) for i in count()), X.columns))\n",
    "\n",
    "# substitute x# with column names or combinations of column names from the dataset,\n",
    "# replaces spaces with underscores to ease column references\n",
    "\n",
    "import seaborn as sns\n",
    "# necessary so `x11` enters the loop before `x1`, preventing duplicate x1 match\n",
    "feature_keys = list(x_to_feature.keys())[::-1]\n",
    "features = []\n",
    "\n",
    "import re\n",
    "for feature in interaction.get_feature_names_out():\n",
    "    feature = re.sub(' ', '_', feature)\n",
    "    features.append(feature)\n",
    "X_all_inter = pd.DataFrame(X_inter, columns=features)\n",
    "X_all_inter_const = sm.add_constant(X_all_inter)\n",
    "\n",
    "lm312 = sm.OLS(y1, X_all_inter_const).fit()\n",
    "display(lm312.summary().tables[1].data[1][4])\n",
    "display(lm312.summary().tables[1].data[2][4])\n",
    "display(lm312.summary().tables[1].data[3][4])\n",
    "display(lm312.summary().tables[1].data[4][4])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5890d6fd-e11b-43e9-9f38-d543119e61fa",
   "metadata": {
    "tags": []
   },
   "source": [
    "**d. Calculate R2 and RAdj 2 for this model. Compare these values to the R2 and RAdj 2 for the simple linear regression model relating clathrate formation to time. Discuss your results.**\n",
    "\n",
    "The computed r-squared and adjusted r-squared values for both models is summarized in the table below:\n",
    "\n",
    "| model   | r-squared           | adj. r-squared     |\n",
    "|---------|---------------------|--------------------|\n",
    "| y~x1+x2 | 0.8415008375400974  | 0.8318948276940428 |\n",
    "| y~x2    | 0.46761604551975844 | 0.4519576939173985 |\n",
    "\n",
    "R-squared values are also known as coefficients of determination because this value indicates how much the response fit estimate, y_hat, is determined  by the model predictor variables. In this case, time (x2) alone determines between 45.2% and 46.8% of the y_hat fit, whereas time (x2) and surfuctant (x1) together determine between 83.2% and 84.2% of the y_hat fit. In this way, the 2 variable model is much better at estimating y values than the single variable model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "8e5db4e8-fa47-459e-9d6b-242380f70907",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8415008375400974"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "0.8318948276940428"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "0.46761604551975844"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "0.4519576939173985"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "\n",
    "# Import clathrate formation data\n",
    "dataset = pd.read_csv(\"https://raw.githubusercontent.com/forestwhite/Regression/main/data/data-table-B8.csv\")\n",
    "\n",
    "# Create Xs and y;\n",
    "X = dataset[['x1','x2']]\n",
    "y1 = dataset['y']\n",
    "\n",
    "X1 = sm.add_constant(X)\n",
    "result = sm.OLS(y1, X1).fit()\n",
    "#print dir(result)\n",
    "display(result.rsquared, result.rsquared_adj)\n",
    "# 0.8415008375400974 0.8318948276940428\n",
    "\n",
    "X = dataset[['x2']]\n",
    "X2 = sm.add_constant(X)\n",
    "result = sm.OLS(y1, X2).fit()\n",
    "#print dir(result)\n",
    "display(result.rsquared, result.rsquared_adj)\n",
    "# 0.46761604551975844 0.4519576939173985"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cad2d31-8bcb-47e3-9797-53d00057cf86",
   "metadata": {
    "tags": []
   },
   "source": [
    "**e. Find a 95% CI for the regression coefficient for time for both models in part d. Discuss any differences.**\n",
    "\n",
    "From the linear models, the confidence intervals for beta_2 (time coefficient) are:\n",
    "\n",
    " * surfuctant+time: 0.08862340997573717 < beta_2 < 0.12924542726171231\n",
    " * time alone: 0.06136762210397215 < beta_2 < 0.13403371454369178\n",
    " \n",
    "When time is the only predictor, the 95% confidence interval for beta_2 coefficient is wider than the 95% confidence interval for beta_2 when there are 2 predictor variables, no matter which method used to calculate the interval. This implies that the model with 2 predictors is more precise than the model with only the time predictor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "99012edc-174a-4a70-8eba-0af37bccab6e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.08862340997573717"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "0.12924542726171231"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "0.06136762210397215"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "0.13403371454369178"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "\n",
    "# Import clathrate formation data\n",
    "dataset = pd.read_csv(\"https://raw.githubusercontent.com/forestwhite/Regression/main/data/data-table-B8.csv\")\n",
    "\n",
    "# Create Xs and y for 2-parameter model\n",
    "X = dataset[['x1','x2']]\n",
    "y1 = dataset['y']\n",
    "\n",
    "X1 = sm.add_constant(X)\n",
    "result = sm.OLS(y1, X1).fit()\n",
    "\n",
    "display(result.t_test('x2').conf_int(alpha=0.05)[0,0])\n",
    "display(result.t_test('x2').conf_int(alpha=0.05)[0,1])\n",
    "\n",
    "# Create x and y for 2-parameter model\n",
    "X = dataset[['x2']]\n",
    "X2 = sm.add_constant(X)\n",
    "result = sm.OLS(y1, X2).fit()\n",
    "\n",
    "display(result.t_test('x2').conf_int(alpha=0.05)[0,0])\n",
    "display(result.t_test('x2').conf_int(alpha=0.05)[0,1])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}