{ "cells": [ { "cell_type": "markdown", "id": "565fd595", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Linear regression in Python" ] }, { "cell_type": "markdown", "id": "0e698b4b", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Goals of this lecture\n", "\n", "- Implementing linear regression in `statsmodels`. \n", "- Interpreting model summaries:\n", " - Interpreting *coefficients* for **continuous** data. \n", " - Interpreting *coefficients* for **categorical** data. " ] }, { "cell_type": "markdown", "id": "52fa3450", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "##### Libraries" ] }, { "cell_type": "code", "execution_count": 1, "id": "cc9e5aa0", "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/anaconda3/lib/python3.9/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.1\n", " warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n" ] } ], "source": [ "## Imports\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import statsmodels.api as sm\n", "import statsmodels.formula.api as smf\n", "import seaborn as sns\n", "import scipy.stats as ss" ] }, { "cell_type": "code", "execution_count": 2, "id": "23f12fbc", "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "%matplotlib inline\n", "%config InlineBackend.figure_format = 'retina' # makes figs nicer!" ] }, { "cell_type": "markdown", "id": "59ec0fe4", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Regression in Python\n", "\n", "There are many packages for implementing **linear regression** in Python, but we'll be focusing on [`statsmodels`](https://www.statsmodels.org/stable/index.html)." ] }, { "cell_type": "markdown", "id": "89424277", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Loading our first dataset" ] }, { "cell_type": "code", "execution_count": 3, "id": "a1633468", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
EducationSeniorityIncome
021.586207113.10344899.917173
118.275862119.31034592.579135
212.068966100.68965534.678727
\n", "
" ], "text/plain": [ " Education Seniority Income\n", "0 21.586207 113.103448 99.917173\n", "1 18.275862 119.310345 92.579135\n", "2 12.068966 100.689655 34.678727" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_income = pd.read_csv(\"data/models/income.csv\")\n", "df_income.head(3)" ] }, { "cell_type": "markdown", "id": "f72ef247", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Exploratory visualization\n", "\n", "Does `Education` correlate with `Income`?" ] }, { "cell_type": "code", "execution_count": 4, "id": "f8051a3a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAw0AAAILCAYAAACqz3BJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABYlAAAWJQFJUiTwAABEBklEQVR4nO3de3xcd33n/9fHM6Ob6yQy0NLa/bl2i00WWrCNuyU4cQytS7lXhN1sL+AEE5dly9K4rbtAu2mXXsKvW+5LHWpstgFCmxq626a02ybxBQg1likthFvlGJwW2GCRiEiWZqTv/jFHRpGl8Uie0dxez8djHsc653vO+UrHI533nO8lUkpIkiRJ0nyWNboCkiRJkpqboUGSJElSRYYGSZIkSRUZGiRJkiRVZGiQJEmSVJGhQZIkSVJFhgZJkiRJFRkaJEmSJFVkaJAkSZJUkaFBkiRJUkWGBkmSJEkVGRokSZIkVZRvdAU6UUScAi4DHmhwVSRJktS+fgB4JKW09lIPZGhojMt6e3tXXnnllSsbXRFJkiS1p/vvv5+xsbGaHMvQ0BgPXHnllStPnDjR6HpIkiSpTW3evJnBwcEHanEs+zRIkiRJqsjQIEmSJKkiQ4MkSZKkigwNkiRJkioyNEiSJEmqyNAgSZIkqSJDgyRJkqSKDA2SJEmSKjI0SJIkSarI0CBJkiSpIkODJEmSpIpaLjRExHUR8Y6IOBoRj0REiojbL7LPVRFxV0ScjYjRiPhMRLwuInIV9nlFRPx9RHw7Ih6OiHsj4gW1/44kSZLUrkYnSoxNTFKcnGJsYpLRiVKjq7Qo+UZXYBHeCDwN+DZwBnhypcIR8WLgz4BzwIeAs8ALgbcAzwJeNsc+vw/syY7/HqALuB743xHxiymld9bqm5EkSVL7GStOMjJWZN+RIQ4NnmF4tEh/X4GBTavZvW0dK3oK9Bbm/fy66bTckwbgl4D1wGXAqysVjIjLKN/0TwLXppRemVL6FeDpwCeA6yLi+ln7XEU5MPwz8CMppV9KKb0G2Ew5cPx+RPxATb8jSZIktY2x4iSDp4e5+s33sP/YKYZHiwAMjxbZf+wUV996DydPDzNWnGxwTavXcqEhpXRPSulLKaVURfHrgCcAd6SUPjXjGOcoP7GAC4PHL2TL304pDc/Y5wHgXUA3cMMiqy9JkqQ2N3KuyI0HjzNemppz+3hpihsOHmfkXHGJa7Z4LRcaFujZ2fKjc2w7AowCV0VEd5X7/NWsMpIkSdJ5oxMl9h0emjcwTBsvTXHb4SHGWqSPQ7uHhg3Z8ouzN6SUSsApyv061gFExHJgFfDtlNK/znG8L2XL9dWcPCJOzPXiIv0wJEmS1JqC4NDgmarKHjr5IBD1rVCNtGJH6IW4PFs+PM/26fVXLLK8JEmS5jE6USII8rmgNJlIJPq62vv2M5+L830YLubsoxPkc4aGVjB9larpHzFTVeVTSpvnPGn5acOmBZ5TkiSpJbTbyEELUZpM9PcVqgoOK5d3UZpMtMKPot2bJ00/Gbh8nu2XzSp3sfIXexIhSZLU0dpx5KCFSCQGNq2uquzAxlUs/LPrxmj30PCFbHlBH4SIyANrgRIwBJBSehR4EPiuiPjeOY73pGx5QR8JSZIktefIQQvR15Vn97Z1dOcr32Z355dx07Z19LZIc612Dw13Z8vnzrHtGqAP+HhKabzKfX5qVhlJkiRl2nXkoIVa0VPgwM4t8waH7vwyDuzcwoqewhLXbPHaPTTcCTwEXB8Rz5heGRE9wJuyL989a58/zJZviIj+Gfv8APAaYBw4UK8KS5Iktap2HTlooXoLOTau6efo3u3s2rqWlcu7gHIfhl1b13J073Y2rulvqX4drfE8ZIaIeAnwkuzLJ2bLZ0bEwezfD6WUfhkgpfRIRLyKcni4NyLuoDyr84soD8d6J/ChmcdPKX08Iv4AuBn4TETcCXQB/x5YCfxiNtGbJEmSZmjXkYMWo7eQo7eQ4+Yd69mzY8P5EaQgtUyTpJlar8bwdOAVs9aty14Ap4Ffnt6QUvpIRGwD3gC8FOgBvkw5FLx9rpmlU0p7IuIzwH8CbgKmgEHg/08p/UVNvxtJkqQ20a4jB12KmUPMtvL32nKhIaV0C3DLAvf5GPC8Be7zPuB9C9lHkiSpk02PHLT/2KmLlm2lkYPU/n0aJEmStETadeQgGRokSZJUQ+04cpAMDZIkSaqhdhw5SC3Yp0GSJKmWRidKBHF+dJtEekznVS1cu40cJEODJEnqUGPFSUbGiuw7MsShwTMMjxbp7yswsGk1u7etY0VPwU/DL1G7jBwkQ4MkSepAY8VJBk8Pc+PB44+ZvXh4tMj+Y6e4/b7THNi5xWY0UsY+DZIkqeOMnCteEBhmGi9NccPB44ycq26iMqndGRokSVJHGZ0ose/w0LyBYdp4aYrbDg8xNlFaoppJzcvQIEmSOkoQHBo8U1XZQycfBKK+FZJagKFBkiR1lHwuGB6trtnR2UcnyOcMDZKhQZIkdZTSZKK/r7qJxVYu78qGCpU6m6FBkiR1lERiYNPqqsoObFwFGBokQ4MkSeoofV15dm9bR3e+8m1Qd34ZN21b52RkEoYGSZLUgVb0FDiwc8u8waE7v4wDO7ewoqe6ZkxSuzM0SJKkjtNbyLFxTT9H925n19a1rFzeBZT7MOzaupaje7c7sZs0g8/bJElSR+ot5Ogt5Lh5x3r27NhAPhdZp+dkkyRpFt8RkiSpo/XNCAg+WJDmZvMkSZIkSRUZGiRJkiRVZGiQJEmSVJGhQZIkSVJFdoSWJEktYXSiRBDnRzlKpMd0YpZUP77TJElSUxsrTjIyVmTfkSEODZ5heLRIf1+BgU2r2b1tHSt6Cs6nINWZoUGSJDWtseIkg6eHufHgccZLU+fXD48W2X/sFLffd5oDO7c4EZtUZ/ZpkCRJTWvkXPGCwDDTeGmKGw4eZ+RccYlrJnUWQ4MkSWpKoxMl9h0emjcwTBsvTXHb4SHGJkpLVDOp8xgaJElSUwqCQ4Nnqip76OSDQNS3QlIHMzRIkqSmlM8Fw6PVNTs6++gE+ZyhQaoXQ4MkSWpKpclEf1+hqrIrl3dRmkx1rpHUuQwNkiSpKSUSA5tWV1V2YOMqwNCwWKMTJcYmJilOTjE2Mcmo/UM0i0OuSpKkptTXlWf3tnXcft/pip2hu/PLuGnbOnqd6G3BnAND1fJJgyRJaloregoc2LmF7vzctyzd+WUc2LmFFT3VNWPSd0zPgXH1m+9h/7FT5/uPTM+BcfWt93Dy9DBjxckG11TNwNAgSZKaVm8hx8Y1/Rzdu51dW9eycnkXUO7DsGvrWo7u3e7EbovkHBhaCJ/jSZKkptZbyNFbyHHzjvXs2bGBfC6yTs/JJkmLtNA5MPbsWL+on/XoRIkgzl+zRKLPa9aSvGqSJKklzLzZ9MHCpVnoHBh7dmxY0PHtK9F+DA2SJEkdpp5zYEz3lZjd9Gm6r8Tt953mwM4tNitrMfZpkCRJ6jD1nAPDvhLtydAgSZLUYeo1B8ZC+0qMOR9Ey7B5kiRJUoep1xwYM/tKPOX7LuPnfmwNP37ld/Nd3QW+PV7k/3zuG7z/k6f57L88sqi+EmocQ4MkSVIHmp4D44Z5mhItZg6MfC545FyJ337JU7nqhx7P//zEA7z5o58/3xH6JRtX8c6f2cTHv/wQv/G/PrugvhJqrEjJKdeXWkSc2LRp06YTJ040uiqSJKmDjRUnGTlX5LbDQxw6+SBnH51g5fIuBjau4qZFjHI0NjHJX3zmX7iir4v/9IHBecPIu352E2cfneCFP/J99HbZGbpeNm/ezODg4GBKafOlHssnDZIkSR2q5nNgROKZP/g4nvPfD1fsCP2a9w/yd3u2EeGH163CjtCSJEkdrq8rT29XjkJuGb1duUVPmpcSHPjYA1V1hD7wsQewwUvrMDRIkiR1uNGJEmMTkxQnpxibmGR0kaMaLWTSuA+ffBCwT0OrsHmSJElSh6r1zM31nDROjeWTBkmSpA40PXPz1W++h/3HTp2/2Z+eufnqW+/h5OlhxoqTVR+znpPGqbEMDZIkSR2oHjM312vSODWeoUGSJKnD1Gvm5ulJ47rzlW8xFzppnBrP0CBJktRhFtJh+dACOyxPTxo3X3BYzKRxajxDgyRJUoepZ4fl3kKOjWv6Obp3O7u2rmXl8i6g3Idh19a1HN27nY1r+hfUwVqN5zMhSZKkDjPdYbma4DDdYXkh9/g1nzRODeeTBkmSpA6zVB2WazVpnBrP0CBJktRh7LCshTI0SJIkdSA7LGshDA2SJEkdyA7LWgifNUmSJHUoOyyrWv5vkCRJ6nB9MwKCDxY0F5snSZIkSarI0CBJkiSpIkODJEmSpIoMDZIkSZIqMjRIkiRJqsjQIEmSJKkiQ4MkSZKkigwNkiRJkioyNEiSJEmqyNAgSZIkqSJDgyRJkqSKDA2SJEmSKjI0SJIkSarI0CBJkiSpIkODJEmSpIoMDZIkSZIqMjRIkiRJqsjQIEmSJKkiQ4MkSZKkigwNkiRJkioyNEiSJEmqyNAgSZIkqSJDgyRJkqSKDA2SJEmSKuqY0BARz4+Iv4mIMxExFhFDEfGnEfHMecpfFRF3RcTZiBiNiM9ExOsiIrfUdZckSZIaqSNCQ0TcCvwFsAn4KPA2YBB4MfCxiPi5WeVfDBwBrgE+DLwL6ALeAtyxdDWXJEnSUhidKDE2MUlxcoqxiUlGJ0qNrlJTyTe6AvUWEU8Efhn4OvAjKaVvzNi2Hbgb+C3g9mzdZcB7gEng2pTSp7L1v56VvS4irk8pGR4kSZJa3FhxkpGxIvuODHFo8AzDo0X6+woMbFrN7m3rWNFToLdgQ5NOeNKwhvL3+cmZgQEgpXQPMAI8Ycbq67Kv75gODFnZc8Absy9fXdcaS5Ikqe7GipMMnh7m6jffw/5jpxgeLQIwPFpk/7FTXH3rPZw8PcxYcbLBNW28TggNXwImgB+NiMfP3BAR1wArgL+dsfrZ2fKjcxzrCDAKXBUR3XWoqyRJkpbIyLkiNx48znhpas7t46Upbjh4nJFzxSWuWfNp++ZJKaWzEbEX+APgcxHxEeCbwA8CLwL+D7B7xi4bsuUX5zhWKSJOAU8B1gH3Vzp3RJyYZ9OTF/I9SJIkqbZGJ0rsOzw0b2CYNl6a4rbDQ+zZsZ7erra/dZ5XJzxpIKX0VmCAckh6FfBrwMuArwIHZzVbujxbPjzP4abXX1HzikqSNA87aUq1FQSHBs9UVfbQyQeBWNR52uW92xFxKSJ+Ffgd4O3AO4GvUf60/3eB90fE01NKv1rt4bJluljBlNLmeepzgvJITpIkVWQnTak+8rk434fhYs4+OkE+t7DQ0G7v3bZ/0hAR1wK3Av8rpXRzSmkopTSaUhoEfhp4ENgTEeuyXaafJFx+wcHKLptVTpKkurCTplQ/pclEf1+hqrIrl3dRmrzo58XnteN7t+1DA/CCbHnP7A0ppVHg7yn/HDZmq7+QLdfPLh8ReWAtUAKGal5TSZJmsJOmVD+JxMCm1VWVHdi4iioamZzXju/dTggN06McPWGe7dPrJ7Ll3dnyuXOUvQboAz6eUhqvTfUkSbrQQjtpjrVoO2mpUfq68uzeto7ufOXb4e78Mm7atq7qTtDt+t7thNBwNFveFBGrZm6IiJ8CngWcAz6erb4TeAi4PiKeMaNsD/Cm7Mt317XGkqSOt1SdNKVOtqKnwIGdW+YNDt35ZRzYuYUVPdU1Y4L2fe92QkfoOynPw/DjwP0R8WHKHaGvpNx0KYBfSyl9EyCl9EhEvCrb796IuAM4S3l41g3Z+g8t+XchSeoo9e6kKQl6Czk2runn6N7t3HZ4iEMnH+TsoxOsXN7FwMZV3LSIDsvt+t5t+9CQUpqKiOcBrwGup9z5uY9yELgLeHtK6W9m7fORiNgGvAF4KdADfBm4OStffaM2SZIWYbqTZjU3H9OdNFtoIBapafQWcvQWcty8Yz17dmwgn4us03Na1LwM7frebfvQAJBSKgJvzV7V7vMx4Hl1qpIkSRVNd9Lcf+zURcsutJOmpAv1zQgIl3IT367v3U7o0yBJUsupVydNSfXVru9dQ4MkSU2qHp00JdVfO753DQ2SJDWpmZ00d21dy8rlXUC5HfSurWs5unc7G9f0t9SsslInaMf3btind+lFxIlNmzZtOnHiRKOrIklqEaMTJYK45E6akpZWI9+7mzdvZnBwcDCltPlSj+VvG0mSWkCtOmlKWlrt8t61eZIkSZKkigwNkiRJkioyNEiSJEmqyNAgSZIkqSJDgyRJkqSKDA2SJEmSKnLIVUmSWsDssd4T6TFDOXYCfwZS4/hOkySpiY0VJxkZK7LvyBCHBs8wPFqkv6/AwKbV7N62jhU9hZaaVXYx/BlIjWdokCSpSY0VJxk8PcyNB48zXpo6v354tMj+Y6e4/b7THNi5hY1r+tv2ptmfgdQc7NMgSVKTGjlXvOBmeabx0hQ3HDzOyLniEtds6fgzkJqDoUGSpCY0OlFi3+GheW+Wp42Xprjt8BBjE6UlqtnS8WcgNQ9DgyRJTSgIDg2eqarsoZMPAlHfCjWAPwOpeRgaJElqQvlcMDxaXZObs49OkM+13w2zPwOpeRgaJElqQqXJRH9foaqyK5d3UZpMda7R0vNnIDUPQ4MkSU0okRjYtLqqsgMbVwHtd8Psz0BqHoYGSZKaUF9Xnt3b1tGdr/ynuju/jJu2raO3DSc582cgNQ9DgyRJTWpFT4EDO7fMe9PcnV/GgZ1bWNFTXROeVuTPQGoOhgZJkppUbyHHxjX9HN27nV1b17JyeRdQbr+/a+taju7d3vaTmvkzkJpDpGT7v6UWESc2bdq06cSJE42uiiSpRYxOlAiCfC6yDr/pkpvjzD5mItHXxE186vEzkNrZ5s2bGRwcHEwpbb7UY/lOkySpBcy8mb/UD9XHipOMjBXZd2SIQ4NnGB4t0t9XYGDTanZvW8eKnkJTfnJfy5+BpIUxNEiS1EHGipMMnh7mxoPHHzPT8vBokf3HTnH7fac5sHOLTX4kPYZ9GiRJ6iAj54oXBIaZxktT3HDwOCPnqptUbSmNTpQYm5ikODnF2MQkoxOlRldJ6hg+aZAkqUOMTpTYd3ho3sAwbbw0xW2Hh9izY31T9Blo1eZUUjvxSYMkSR0iCA4Nnqmq7KGTDwJR3wpVYbo51dVvvof9x04xPFp+AjLdnOrqW+/h5OlhxoqTDa6p1N4MDZIkdYh8Ls7fdF/M2UcnyOcaHxpauTmV1E4MDZIkdYjSZKK/r7pJ0FYu78qGNW2chTanGrOPg1Q3hgZJkjpEIjGwaXVVZQc2rgIaGxpasTmV1K4MDZIkdYi+rjy7t62jO1/5z393fhk3bVvX8E7QrdicSmpXhgZJkmqo2YcFXdFT4MDOLfMGh+78Mg7s3MKKnuqaMdVTqzWnktpZ48dRkySpDbTKsKC9hRwb1/RzdO92bjs8xKGTD3L20QlWLu9iYOMqbmqiuk43p9p/7NRFyzZDcyqpnUVKvsGWWkSc2LRp06YTJ040uiqSpBqYb5bladOf3jfbLMujEyWCIJ+L7FP61PAmSbN9Y+QcV996T8XO0N35ZRzdu53vXtGzhDWTmt/mzZsZHBwcTCltvtRj2TxJkqRL1KrDgvZ15entylHILaO3K9d0gQFaqzmV1M4MDZIkXQKHBa2vmc2pdm1dy8rlXUC5D8OurWs5und70z3BkdpR832kIElSC1nosKB7dmyoc43aT28hR28hx8071rNnx4ambk4ltSvfaZIkXQKHBV06fTMCgg8WpKVl8yRJki6Bw4JK6gSGBkmSLkGrzbIsSYthaJAk6RK02izLkrQYhgZJki6Rw4Kq1TX7TOZqPD/ukCTpErXSLMvSTK0yk7kazxmhG8AZoSWpfbXCLMsStO5M5qqeM0JLktSkWmGWZQladyZzNYahQZIkqcM4k7kWytAgSZLUYRY6kzk4KWGnMzRIkiR1GGcy10IZGiRJqiGHrlQrcCZzLZS9syRJqgGHrlQrmZ7JfP+xUxct60zmAp80SJJ0yaaHrrz6zfew/9ip880+hkeL7D92iqtvvYeTp4cZK042uKZSmTOZa6EMDZIkXSKHrlQrciZzLYShQZKkS+DQlWpVM2cy37V1LSuXdwHlPgy7tq7l6N7tTuym83zWJEnSJVjo0JV7dmyoc42k6vUWcvQWcty8Yz17dmxwJnPNy/8NkiRdAoeuVDvomxEQfLCgudg8SZKkS+DQlZI6gaFBkqRLMD10ZTUculJSqzI0SJJ0CRy6UlInMDRIknSJHLpSUrurS2iIiBdGxB0R8Q8R8eUZ66+MiF+NiFX1OK8kSY3g0JWS2l1Nn5FGRAAHgZ/LVo0BvTOKDAO/AwRway3PLUlSIzl0paR2VusnDf8R+HngALAS+P2ZG1NKXwM+Bjy/xueVJKkp9HXl6e3KUcgto7crZ2CQ1BZqHRpeCfwD8KqU0sPMPUTEl4C1NT6vJEmSpDqpdWjYANyTUqo0ntw3gCfU+LySJEmS6qTWoaEE9FykzCrg2zU+ryRJkqQ6qXVo+BxwbdYh+gIR0QM8GzhZ4/NKkiRJqpNah4Y/Bp4MvCUiHnPsiMgBfwB8H+URliRJkiS1gFoP6bAPeBHwWuBlwAhARNwJ/BjlwPDnKaX31/i8kiRJkuqkpk8aUkqTwAuA3wK6gPWU52QYAPqA/0Y5TEiSJElqETUfPDqlVAJuiYjfpBwaHgc8DHw+CxWSJEmSWkjdZpzJhl39Qr2OL0mSJGlp1LojtCRJkqQ2U/MnDRGxGvgl4OnAaqAwR7GUUvrBWp9bkiRJUu3VNDRExLXAXZQneCsBX8+WFxSt5XklSZIk1U+tnzS8GcgBLwc+kFKaqvHxJUmSJC2xWoeGHwY+mFK6vcbHlSRJktQgte4IPQycrfExJUmSJDVQrUPDXwDbanxMSZIkSQ1U69DweuDyiHhXRCyv8bElSZIkNUBN+zSklB6KiOcCnwReHhFfpDwb9BxF03Nqee5qRMTVwOuAq4CVlJtS/SPw1pTSXbPKXgW8EfgxyqNBfRl4L/AOZ7aWJElSJ6n1kKtPAe4B+rNVG+cpmmp53mpExBuB/wY8RLkZ1b8Cj6dcx2spDxU7XfbFwJ8B54APUQ4XLwTeAjwLeNkSVl2SJElqqFqPnvQHwOOA3wDeB/xLM3wqHxEvoxwY/hYYSCmNzNpemPHvy4D3AJPAtSmlT2Xrfx24G7guIq5PKd2xVPWXJEmSGqnWfRqeCRxKKb0ppfTVJgkMy4BbgVHgZ2YHBoCUUnHGl9cBTwDumA4MWZlzlJsrAby6fjWWJEmSmkutnzRMAA/U+JiX6ipgLXAnMBwRzweeSrnp0d+nlD4xq/yzs+VH5zjWEcrh46qI6E4pjdepzpIkSVLTqHVouBf40Rof81JtyZZfBwYpT0B3XkQcAa5LKf3fbNWGbPnF2QdKKZUi4hTwFGAdcH+lE0fEiXk2Pbm6qkuSJEmNV+vmSb8K/JuI+LWIiBofe7G+O1v+AtAL/DiwgvLThr8GrgH+dEb5y7PlXKM+zVx/RU1rKUmSJDWpWj9peCPwT8BvA6+KiE8z/5Crr6zxueeTy5ZB+YnCP2RffzYifpryE4VtEfHMOZoqzWU6DF10BKiU0uY5D1B+ArGpinNJkiRJDVfr0LBzxr/XZq+5JGCpQsNwthyaERjKlUhpLCL+OqvLjwKf4Dsh53Lmdlm2nO9JhCRJktRWah0a5gsJjfSFbPmtebZPh4reGeWfAawHHtMnISLylL/HEjBU01pKkiRJTarWM0KfruXxauQI5Zv8J0VEV0ppYtb2p2bLB7Ll3cDPAs8FPjir7DVAH3DEkZMkSZLUKWrdEbrppJQeojyr8+WUJ507LyJ+AvhJyk2NpodYvZPyrNHXR8QzZpTtAd6UffnuOldbktSiRidKjE1MUpycYmxiktGJUqOrJEmXrNbNkwCIiB8DdgEbKY8y9DDlpj4HUkofr8c5L+Jm4N8Cb4iIa4C/B9YAP0155udXpZS+BZBSeiQiXkU5PNwbEXcAZ4EXUR6O9U7KIUSSpPPGipOMjBXZd2SIQ4NnGB4t0t9XYGDTanZvW8eKngK9hdzFDyRJTajmoSEi3gT8F74zytC0pwM3RsStKaXX1/q8laSUvhER/5by6E4/DfwYMAL8JfC7KaX7ZpX/SERsA94AvBToAb5MOXy8PaV00ZGTJEmdY6w4yeDpYW48eJzx0tT59cOjRfYfO8Xt953mwM4tbFzTb3CQ1JJq2jwpIl4GvB74CuUnDesodzBel339FWBvRPy7Wp63Gimlsymlm1NKa1NKXSmlx6WUXjw7MMwo/7GU0vNSSv0ppd6U0g+nlN6SUppc6rpLkprbyLniBYFhpvHSFDccPM7IueIS10ySaqPWfRp+kfLMy1tSSu9NKT2QUhrPlu+lPDvz/wVeU+PzSpLUEKMTJfYdHpo3MEwbL01x2+EhxuzjIKkF1To0PA24M+t8fIFs/Z9SbqokSVLLC4JDg2eqKnvo5INc2HpXkppfrUNDHhi9SJlR6tQBW5KkpZbPBcOj1TU7OvvoBPmcoUFS66l1aPgy8IKImPO42frnAf9c4/NKktQQpclEf1+hqrIrl3dRmnQsDUmtp9ah4YPAlcCfR8STZm6IiB+kPFzpvwE+UOPzSpLUEInEwKbVVZUd2LgKMDRIaj21Dg1/QHkG5ucD90fEVyLikxFxGvgC8BLgY1k5SZJaXl9Xnt3b1tGdr/wntTu/jJu2raO3yxa6klpPTUNDSmkC+AnK8xucAlZTHjHp+7Ov3wA8JysnSVJbWNFT4MDOLfMGh+78Mg7s3MKKnuqaMUlSs6n5xx0ppSLwu8DvRsR3AZcDD6eUvl3rc0mS1Ax6Czk2runn6N7t3HZ4iEMnH+TsoxOsXN7FwMZV3OSM0JJaXF2fkWZBwbAgSWp7vYUcvYUcN+9Yz54dG8jnIuv0nGySJKnl1XpG6M0R8RsR8T3zbH9itv3ptTyvJEnNoq8rT29XjkJuGb1dOQODpLZQ647Qe4BdwDfm2f514JXAzTU+ryRJkqQ6qXVoeCZwT0ppzvHksvV3A8+q8XklSZIk1UmtQ8MTgTMXKfMvwPfW+LySJEmS6qTWoWEUeMJFyjwBGK/xeSVJkiTVSa1Dw6eBF2dDrV4gIi4DXpyVkyRJktQCah0abqP8JOH/RMSPzNwQEU8D/gZ4fFZOkiRJUguo6ThwKaUPRcRPAS8HTkbE14EHgVXA9wABvC+l9MFanleS2tXoRIkgzo/5n0j0OYSnJGmJ1WNG6J0R8XHgF4GnUO4cDfBPwNtTSn9U63NKUrsZK04yMlZk35EhDg2eYXi0SH9fgYFNq9nt7MKSpCVWl4+rUkq3AbdFRB9wBfCtlNJoPc4lSe1mrDjJ4Olhbjx4nPHS1Pn1w6NF9h87xe33nebAzi1sXNNvcJAkLYla92l4jJTSaErpXwwMklS9kXPFCwLDTOOlKW44eJyRc8UlrpkkqVPVNTRIkhZmdKLEvsND8waGaeOlKW47PMTYRGmJaiZJ6mQ1Dw0RsS0i/iIivhERxYiYnOPlXzlJmkMQHBq82ByZZYdOPkh5fAlJkuqrpn0aIuL5wEeAHPAV4AuAAUGSqpTPBcOj1TU7OvvoBPmcoUGSVH+17gh9C1AEnp9S+psaH1uS2l5pMtHfV6gqOKxc3kVpMmFfaElSvdW6edJTgQ8ZGCRpcRKJgU2rqyo7sHEVkOpbIUmSqH1o+DZwtsbHlKSO0deVZ/e2dXTnK/967s4v46Zt6+h1ojdJ0hKodWj4O+CZNT6mJHWUFT0FDuzcMm9w6M4v48DOLazoKSxxzSRJnarWoWEv8IMR8caIsHeeJC1CbyHHxjX9HN27nV1b17JyeRdQ7sOwa+taju7d7sRukqQlVevn2v8V+Czwm8CNEfFp4FtzlEsppVfW+NyS1DZ6Czl6Czlu3rGePTs2kM8FpckEJJskSZKWXK3/8uyc8e8fyF5zSYChQZIuom9GQPDBgiSpUWodGtbW+HiSJEmSGqymoSGldLqWx5MkSZLUeLXuCC1JkiSpzRgaJEmSJFV0yc2TImJyEbullJLDf0iSJEktoBY37ouZj8E5HCRJkqQWccmhIaVkEydJkiSpjXnDL0mSJKkiQ4MkSZKkigwNkiRJkioyNEiSJEmqyNAgSZIkqSJDgyRJkqSKDA2SJEmSKjI0SJIkSarI0CBJkiSpIkODJEmSpIoMDZIkSZIqMjRIkiRJqsjQIEmSJKkiQ4MkSZKkigwNkiRJkioyNEiSJEmqyNAgSZIkqaJ8oyugpTE6USII8rmgNJlIJPq6vPySJEm6OO8a29xYcZKRsSL7jgxxaPAMw6NF+vsKDGxaze5t61jRU6C3kGt0NSVJktTEDA1tbKw4yeDpYW48eJzx0tT59cOjRfYfO8Xt953mwM4tbFzTb3CQJEnSvOzT0MZGzhUvCAwzjZemuOHgcUbOFZe4ZpIkSWolhoY2NTpRYt/hoXkDw7Tx0hS3HR5ibKK0RDWTJElSqzE0tKkgODR4pqqyh04+CER9KyRJkqSWZWhoU/lcMDxaXbOjs49OkM8ZGiRJkjQ3Q0ObKk0m+vsKVZVdubyL0mSqc40kSZLUqgwNbSqRGNi0uqqyAxtXAYYGSZIkzc3Q0Kb6uvLs3raO7nzlS9ydX8ZN29bR60RvkiRJmoehoY2t6ClwYOeWeYNDd34ZB3ZuYUVPdc2YJEmS1JkMDW2st5Bj45p+ju7dzq6ta1m5vAso92HYtXUtR/dud2I3SZIkXZRtUtpcbyFHbyHHzTvWs2fHBvK5yDo9J5skSZIkqSreNXaIvhkBwQcLkiRJWgibJ0mSJEmqyNAgSZIkqSJDgyRJkqSKDA2SJEmSKjI0SJIkSarI0CBJkiSpIodc7RCjEyWCOD9PQyI9ZhhWSZIkaT7eNba5seIkI2NF9h0Z4tDgGYZHi/T3FRjYtJrd29axoqfgjNCSJEmqyNDQxsaKkwyeHubGg8cZL02dXz88WmT/sVPcft9pDuzcwsY1/QYHSZIkzcs+DW1s5FzxgsAw03hpihsOHmfkXHGJayZJkqRWYmhoU6MTJfYdHpo3MEwbL01x2+EhxiZKS1QzSZIktRpDQ5sKgkODZ6oqe+jkg0DUt0KSJElqWYaGNpXPBcOj1TU7OvvoBPmcoUGSJElz68jQEBE/HxEpe+2ap8xVEXFXRJyNiNGI+ExEvC4iWqLHcGky0d9XqKrsyuVdlCZTnWskSZKkVtVxoSEivh94B/DtCmVeDBwBrgE+DLwL6ALeAtyxBNW8ZInEwKbVVZUd2LgKMDRIkiRpbh0VGiIigAPAN4E/nKfMZcB7gEng2pTSK1NKvwI8HfgEcF1EXL80NV68vq48u7etoztf+RJ355dx07Z19DrRmyRJkubRUaEBeC3wbOAG4NF5ylwHPAG4I6X0qemVKaVzwBuzL19dz0rWyoqeAgd2bpk3OHTnl3Fg5xZW9FTXjEmSJEmdqWNCQ0RcCfwe8LaU0pEKRZ+dLT86x7YjwChwVUR017iKNddbyLFxTT9H925n19a1rFzeBZT7MOzaupaje7c7sZskSZIuqiPapEREHvhj4CvA6y9SfEO2/OLsDSmlUkScAp4CrAPuv8h5T8yz6ckXqUPN9BZy9BZy3LxjPXt2bCCfi6zTc7JJkiRJkqrSKXeNvwFsBLamlMYuUvbybPnwPNun119Rg3otmb4ZAcEHC5IkSVqItg8NEfGjlJ8u/PeU0idqcchsedHhhlJKm+ep0wlgUw3qIkmSJNVdW/dpmNEs6YvAr1e52/SThMvn2X7ZrHKSJElSW2vr0AB8F7AeuBI4N2NCtwT816zMe7J1b82+/kK2XD/7YFkIWQuUgKG61lySJElqEu3ePGkc2D/Ptk2U+zkcoxwUppsu3Q38LPBc4IOz9rkG6AOOpJTGa15bSZIkqQm1dWjIOj3vmmtbRNxCOTS8L6X0RzM23QncClwfEe+YnqshInqAN2Vl3l23SkuSJElNpq1Dw2KklB6JiFdRDg/3RsQdwFngRZSHY70T+FADqyhJkiQtqXbv07AoKaWPANsoT+b2UuAXgSJwM3B9SumiIydJkiRJ7aJjnzSklG4Bbqmw/WPA85aqPpIkSVKz8kmDJEmSpIoMDZIkSZIqMjRIkiRJqsjQIEmSJKkiQ4MkSZKkigwNkiRJkioyNEiSJEmqyNAgSZIkqSJDgyRJkqSKDA2SJEmSKjI0SJIkSaoo3+gKSJrb6ESJIMjngtJkIpHo6/ItK0mSlp53IFKTGStOMjJWZN+RIQ4NnmF4tEh/X4GBTavZvW0dK3oK9BZyja6mJEnqIIYGqYmMFScZPD3MjQePM16aOr9+eLTI/mOnuP2+0xzYuYWNa/oNDpIkacnYp0FqIiPnihcEhpnGS1PccPA4I+eKS1wzSZLUyQwNUpMYnSix7/DQvIFh2nhpitsODzE2UVqimkmSpE5naJCaRBAcGjxTVdlDJx8Eor4VkiRJyhgapCaRzwXDo9U1Ozr76AT5nKFBkiQtDUOD1CRKk4n+vkJVZVcu76I0mepcI0mSpDJDg9QkEomBTaurKjuwcRVgaJAkSUvD0CA1ib6uPLu3raM7X/lt2Z1fxk3b1tHrRG+SJGmJGBqkJrKip8CBnVvmDQ7d+WUc2LmFFT3VNWOSJEmqBUOD1ER6Czk2runn6N7t7Nq6lpXLu4ByH4ZdW9dydO92J3aTJElLzvYNUpPpLeToLeS4ecd69uzYQD4XWafnZJMkSZLUEN6BSE2qb0ZA8MFCbY1OlAjifCBLpMf8vCVJ0mP5V1JSxxgrTjIyVmTfkSEODZ5heLRIf1+BgU2r2b1tHSt6Cjb9kiRpDoYGSR1hrDjJ4Olhbjx4nPHS1Pn1w6NF9h87xe33nebAzi32GZEkaQ52hJbUEUbOFS8IDDONl6a44eBxRs5VNyu3JEmdxNAgqe2NTpTYd3ho3sAwbbw0xW2HhxibKC1RzSRJag2GBkltLwgODZ6pquyhkw8CUd8KSZLUYgwNktpePhcMj1bX7OjsoxPkc4YGSZJmMjRIanulyUR/X3WzaK9c3pXNiyFJkqYZGiS1vURiYNPqqsoObFwFGBokSZrJ0CCp7fV15dm9bR3d+cq/8rrzy7hp2zpn3pYkaRZDg6SOsKKnwIGdW+YNDt35ZRzYuYUVPdU1Y5IkqZMYGiR1hN5Cjo1r+jm6dzu7tq5l5fIuoNyHYdfWtRzdu92J3SRJmofP4CV1jN5Cjt5Cjpt3rGfPjg3kc5F1ek42SZIkqQL/SkrqOH0zAoIPFiRJujibJ0mSJEmqyCcNklQDoxMlgjjf5CmRHvNEQ5KkVuZfNEm6BGPFSUbGiuw7MsShwTMMjxbp7yswsGk1u7etY0VPwc7VkqSWZ2iQpEUaK04yeHqYGw8eZ7w0dX798GiR/cdOcft9pzmwc4ujMkmSWp59GiRpkUbOFS8IDDONl6a44eBxRs4Vl7hmkiTVlqFBukSjEyXGJiYpTk4xNjHJ6ESp0VXSEhidKLHv8NC8gWHaeGmK2w4PMeb/C0lSC7N5krRItmXvbEFwaPBMVWUPnXyQPTs21LlGkiTVj6FBWgTbsiufC4ZHq2t2dPbRCfK5qHONJEmqH5snSYtgW3aVJhP9fYWqyq5c3pXNPC1JUmsyNEgLZFt2ASQSA5tWV1V2YOMqwNAgSWpdhgZpgRbalh1sltKO+rry7N62ju585V+j3fll3LRtHb1O9CZJamGGBmmBbMuuaSt6ChzYuWXe4NCdX8aBnVtY0VNdMyZJkpqVoUFaINuya1pvIcfGNf0c3budXVvXsnJ5F1C+7ru2ruXo3u12hpcktQWfl0sLNN2Wff+xUxcta1v29tdbyNFbyHHzjvXs2bGBfC6yoJhskiRJahs+aZAWyLbsmktfV57erhyF3DJ6u3Jed0lSWzE0SItgW3ZJktRJDA3SItiWXZIkdRKfn0uLZFt2SZLUKbyzkS5R34yA4IMFSZLUjmyeJEmSJKkiQ4MkSZKkigwNkiRJkioyNEiSJEmqyNAgSZIkqSJDgyRJkqSKDA2SJEmSKjI0SJIkSarI0CBJkiSpIkODJEmSpIoMDZIkSZIqMjRIkiRJqijf6ApIrW50okQQ5HNBaTKRSPR1+daSJEntwzsbaZHGipOMjBXZd2SIQ4NnGB4t0t9XYGDTanZvW8eKngK9hVyjqylJknTJDA3SIowVJxk8PcyNB48zXpo6v354tMj+Y6e4/b7THNi5hY1r+g0OkiSp5dmnQVqEkXPFCwLDTOOlKW44eJyRc8UlrpkkSVLtGRqkBRqdKLHv8NC8gWHaeGmK2w4PMTZRWqKaSZIk1YehQVqgIDg0eKaqsodOPghEfSskSZJUZ/ZpUMeo1ShH+VwwPFpds6Ozj06QzxkaJElSazM0qO3VepSj0mSiv69QVXBYubyL0mTCvtCSJKmV2TxJbW16lKOr33wP+4+dOn+jPz3K0dW33sPJ08OMFSerPmYiMbBpdVVlBzauAtJiqi5JktQ0DA1qa/UY5aivK8/ubevozld++3Tnl3HTtnX0OtGbJElqcYYGta16jnK0oqfAgZ1b5g0O3fllHNi5hRU9hQXVWZIkqRkZGtS26jnKUW8hx8Y1/Rzdu51dW9eycnkXUO7DsGvrWo7u3e7EbpIkqW3YbkJtq96jHPUWcvQWcty8Yz17dmw4PyoTJJskSZKkttL2Txoi4nERsSsiPhwRX46IsYh4OCKORcQrI2LOn0FEXBURd0XE2YgYjYjPRMTrIsKPjlvE9ChH1Zge5Wgx+rry9HblKOSW0duVMzBIkqS20/ahAXgZ8B7g3wKfBN4K/BnwVOCPgD+JiMd8xBwRLwaOANcAHwbeBXQBbwHuWKqK69I4ypEkSVJtdEJo+CLwImB1SulnU0r/JaV0I/Bk4KvAS4GB6cIRcRnlkDEJXJtSemVK6VeApwOfAK6LiOuX+HvQIjjKkSRJUm20fWhIKd2dUvrfKaWpWeu/Bvxh9uW1MzZdBzwBuCOl9KkZ5c8Bb8y+fHX9aqxacpQjSZKkS9fpH61O95KdOdbms7PlR+cofwQYBa6KiO6U0nilg0fEiXk2PXlBtdSizRzl6LbDQxw6+SBnH51g5fIuBjau4qZFzAgtSZLUaTo2NEREHnh59uXMgLAhW35x9j4ppVJEnAKeAqwD7q9rJVUTjnIkSZJ0aTr5jun3KHeGviul9Ncz1l+eLR+eZ7/p9Vdc7AQppc1zrc+eQGyqrprNbXSiRBDnb8QTib4mvRGfWS8fLEiSJFWvOe/u6iwiXgvsAT4P/PxCd8+WHT3UzlhxkpGxIvuODHFo8AzDo0X6+woMbFrNbpv8SJIktZWOCw0R8RrgbcDngOeklM7OKjL9JOFy5nbZrHIdZ6w4yeDpYW48eJzx0nf6lw+PFtl/7BS333eaAzu3OCOyJElSm2j70ZNmiojXAe8E/gnYno2gNNsXsuX6OfbPA2spd5weqlM1m97IueIFgWGm8dIUNxw8zsi56mZjliRJUnPrmNAQEXspT872acqB4RvzFL07Wz53jm3XAH3Axy82clK7Gp0ose/w0LyBYdp4aYrbDg8xNlGqWE6SJEnNryNCQ0T8OuWOzycoN0l6qELxO4GHgOsj4hkzjtEDvCn78t31qmuzC4JDg2eqKnvo5IN8pwuIFmp0osTYxCTFySnGJiYZNYBJkqQGafs+DRHxCuC3KM/wfBR4bcQFN7IPpJQOAqSUHomIV1EOD/dGxB3AWcqzSm/I1n9oaWrffPK5YHi0umZHZx+dIJ8zNCyUncwlSVKzafvQQLkPAkAOeN08ZQ4DB6e/SCl9JCK2AW8AXgr0AF8GbgbenlLq2JGTSpOJ/r5CVcFh5fIuSpPJ4U0XwE7mkiSpGbV986SU0i0ppbjI69o59vtYSul5KaX+lFJvSumHU0pvSSlNNuDbaBqJxMCm1VWVHdi4ig4fmXbB7GQuSZKaUduHBtVWX1ee3dvW0Z2v/F+nO7+Mm7atc8blBbCTuSRJalaGBi3Yip4CB3ZumTc4dOeXcWDnFlb0FJa4Zq3NTuaSJKlZGRq0YL2FHBvX9HN073Z2bV3LyuVdQLkPw66tazm6d7tt7hfBTuaSJKlZ2XZEi9JbyNFbyHHzjvXs2bGBfC4oTSYg2SRpkexkLkmSmpVPGnRJ+rry9HblKOSW0duVMzBcAjuZS5KkZmVokJqEncwlSVKzMjRITcRO5pIkqRkZGqQmYidzSZLUjGzfIDUZO5lLkqRm4x2I1KT6ZgQEHyxIkqRGsnmSJEmSpIoMDZIkSZIqMjRIkiRJqsjQIEmSJKkiQ4MkSZKkigwNkiRJkioyNEiSJEmqyNAgSZIkqSJDgyRJkqSKDA2SJEmSKjI0SJIkSaooUkqNrkPHiYhv9vb2rrzyyisbXRVJkiS1qfvvv5+xsbGzKaXHXeqxDA0NEBGngMuAB5b41E/Olp9f4vOqel6j5uc1an5eo+bnNWp+XqPWcLHr9APAIymltZd6IkNDB4mIEwAppc2Nrovm5jVqfl6j5uc1an5eo+bnNWoNS3md7NMgSZIkqSJDgyRJkqSKDA2SJEmSKjI0SJIkSarI0CBJkiSpIkdPkiRJklSRTxokSZIkVWRokCRJklSRoUGSJElSRYYGSZIkSRUZGiRJkiRVZGiQJEmSVJGhQZIkSVJFhoY2ERHXRcQ7IuJoRDwSESkibr/IPldFxF0RcTYiRiPiMxHxuojILVW9O8lCrlFEPCki9kbE3RHx1YiYiIivR8SfR8T2pa57p1jM+2jW/vuzfVJE/FA969qpFvm7LiLiFRFxb/b7biwiTkXEn0TE+qWqe6dY6DWKiO6IeE1E/H1EPBQR346I+yPi7RGxZinr3gki4nERsSsiPhwRX87eDw9HxLGIeGVEzHlv6D3D0lnoNVqqe4Z8rQ6khnsj8DTg28AZ4MmVCkfEi4E/A84BHwLOAi8E3gI8C3hZPSvboRZyjf4b8O+BzwF3Ub4+G4AXAS+KiP+cUnp7favbkRb0PpopIl4I3Jjt+111qZ1g4b/reoA/BV4AfAH4ADACfB9wNbAe+GId69uJqr5GEZEH/o7y353PAx8ExoEtwC8CL4+Iq1JKn6t3pTvIy4B3A/8K3AN8BfgeYAD4I+CnIuJlacbsv94zLLmFXqOluWdIKflqgxewHXgSEMC1QAJun6fsZcA3KP9ifsaM9T3Ax7N9r2/099RurwVeo53AxjnWbwMmsmv3vY3+ntrttZBrNGu/JwBfA+4A7s32+6FGfz/t+FroNQLelZX5HWDZHNsLjf6e2u21wN91L8u2/+3s6wP8ZrbtvY3+ntrpBTyb8g3/7J/3EynfnCbgpTPWe8/Q/NdoSe4ZbJ7UJlJK96SUvpSy/yUXcR3lm5w7UkqfmnGMc5Q/IQJ4dR2q2dEWco1SSgdTSifnWH+Y8k1pF3BV7WvZ2Rb4Pprptmz5mlrXSY+1kGsUET8I/AJwHHhDSmlqjuMV61DNjrbA99G6bPmXc1yfP8+WT6hd7ZRSujul9L9n/7xTSl8D/jD78toZm7xnWGILvUZLdc9g86TO9Oxs+dE5th0BRoGrIqI7pTS+dNVSlaZvckoNrYUAiIidwEuAn04pfTMiGlshzfQfKPfdex9wWdaE7PuBbwJ3p5S+3MjKCYDPZsufioi3zbpJekG2/NslrlMnm+vvi/cMzWWh9wA1u2cwNHSmDdnygna8KaVSRJwCnkL5E6D7l7JiqizrFPgcyr+kjzS4Oh0vux5vo9z04iMNro4utCVbXg78M/C4GdtSRLwbeG1KaXLJa6ZpfwkcotxW+x8j4m8pN6fYDGwF3gG8s3HV6xxZ/5KXZ1/ODAjeMzSJCtdovvI1vWeweVJnujxbPjzP9un1V9S/KqpWRHQD7we6gVtSSsMNrlJHy0aveB/lzp6vbXB1NLfvzpa/BXwK+GFgBeU/ov8M/Efg1xtTNQFkTZiuA26hfHP6WuCXKfeLOAJ8wFC3ZH4PeCpwV0rpr2es956hecx3jS5Qj3sGQ4PmMt2+YqHtulUn2ZB2f0x5lIoPAb/f2BoJ+CXKncxeZYBrWtNDQf4r5eZj/5RS+nZK6W7KN6pTwM0R0dWwGna4bHSrD1EOCq8BvpfyTerzgDXAkWzkHtVRRLwW2EN5BKufX+ju2dJ7hjpayDWq1z2DoaEzTX8qcPk82y+bVU4NlL35b6c8ysifAD+3iI66qqGIeBLw28CBlNJdja6P5jUd5j6aUhqbuSGl9A/AKcpPHq5c6orpvF+j/LvtDSmlfSmlr6WUHkkp/RXlYFeg3ARQdRIRr6H8M/4csD2ldHZWEe8ZGqyKazSzbN3uGQwNnekL2fKCSY2y9nJrKXeYGVrKSulC2fX4IHA95fHlfyalZAfoxnsK5Ue+N8yYzC1FRKL89AHgS9m6lzSslpr+XfetebZPh4re+ldF85ju7HzP7A1ZsDsLrImIx83erksXEa+j3GfknyjfjH5tjmLeMzRQlddoumxd7xnsCN2Z7gZ+Fngu5f9cM10D9AFHHAWhsbImE38CvBj4n8ANcw0ZqYZ4ANg/z7bnUx5L+0+BR7Kyaoy/ozxB2FNnb8ja+z4p+/KBJayTHqs7W14wrGp2jaY/xZ5Yshp1iIjYS7mN/KeBn0gpPTRPUe8ZGmQB12hJ7hl80tCZ7gQeAq6PiGdMr8zalr4p+/LdjaiYyrI/lh+m/Obfj4GhqaSUPp1S2jXXi+98Kvf6bN2nG1jVTvdXlD/9/MmI+IlZ236dcnOLw5U+uVPdHc2Wr89+7810C+UPN4+nlEaWtFZtLiJ+nfLN6AngOZVuRvGeoSEWco2W6p4hbBrdHrImEC/Jvnwi8JOU/1hO/0J+KKX0y7PK30l5Svg7KD8CfhHl0SvuBP6d7eZrayHXKCIOUJ7h8SHgfzB3B7N7U0r31q3CHWih76N5jnEv5SZKT3IegNpbxO+6rcDfUJ7c6MPAacpDsV4D/F9ga0rpgqEktXgL/F23CrgPWE35ic9HgTHKHTh/NPv3c1JKn1ia2re/iHgFcBCYpDyk7Vx9ER5IKR2csc9L8J5hySz0Gi3ZPcNc00T7ar0X5U9kUoXXA3Ps8yzgLsrteseAf6Q8Ikyu0d9PO74Wco0oz+BYqWyiPIRaw7+vdnot5n00xzGmr90PNfr7acfXIn/X/RvKI4h8g3Izl68C+4DVjf5+2vG10GtEuWnS71Me4/9cdo1OAweAJzf6+2m3VxXXJ1G+wZy9n/cMTXqNluqewScNkiRJkiqyT4MkSZKkigwNkiRJkioyNEiSJEmqyNAgSZIkqSJDgyRJkqSKDA2SJEmSKjI0SJIkSarI0CBJkiSpIkODJEmSpIoMDZIkSZIqMjRIkiRJqsjQIElalIi4NiJSRNzS6Lpcioi4NyJSo+shSc3M0CBJHSS7yb/Y69pG17OWIuJg9n39QKPrIkmtKt/oCkiSGuI3K2x7YKkq0SReDvQ1uhKS1MwMDZLUgVJKtzS6Ds0ipfSVRtdBkpqdzZMkSRVFxPdExP6I+HpEjEXEpyPiFRXKPxARD8yz7Zb5mkBFxJMj4r3Z/uMR8Y2IOBoRr55V7iURcXtEfDEiHo2Ib0fEiYh4bUQsm1U2AdN1PTWjCdYDM8rM2achIpZFxC9ExPHsHI9m/3717PNMnys71uMj4raI+Nfs+/hsRNww389LklqBTxokSfOKiMcBHwfWAcey1/cCfwj8TQ3P83zgT4Fu4KPAB4ErgKcBvwq8e0bx3wOmgE8CDwKXA88G3gZsAX5+RtnfBF6SHedtwLey9d/i4v4Y+Bngq8AfAQn4aeB/AFuBn51jnyuAjwETwJ1AD3Ad8N6ImEopva+K80pS0zE0SFIHqjDi0bmU0u/N+Pp3KQeGt6aUfmnG/u8EPlGjujwe+ADlv0nPTikdnrV99axdnp9S+udZZZYBB4CXR8Q7U0qfhHIzrKwD9NOy7+GBKuv0HygHhpPANSmlb2fr3wgcBn4mIv4ypfSBWbs+DdgP7E4pTWb7vAX4DLAXMDRIakmGBknqTP91nvUPU/4kn4goUP40fQS4ZWahlNKnIuL9fKfpz6V4BXAZ8PbZgSE715lZX//zHGWmIuJtlDs1/yTlpxCX4sZs+WvTgSE7z6MRsRf4W2AX5bAz0yhw83RgyPb5XER8DLgmIlaklEYusW6StOTs0yBJHSilFPO8rphR7MmURxX6dErp4TkOc2+NqvNj2fKvqikcEY+LiN+LiM9kfQ1S1ifhRFZkVQ3qtIlyE6h759h2GJgENs6x7UsppUfmWP/VbHlFDeomSUvOJw2SpPlcni2/Ps/2r9XoPFdkywcvVjAirgCOA2uBvwf+J3AWKGXH+c+U+0VcqsuBsymlidkbUkqliHgI+O459vvWPMcrZctcDeomSUvO0CBJms/004XvmWf7E+dZPwV0zbPtijnWfStbrgL+8SJ12kU5MPzm7GFjI+KZlENDLTwMrIyIQkqpOOs8eeDxwFxPFCSpLdk8SZI0n89TbqP/9Ii4fI7t186z3zDwPVmfiNmeMce6+7LlT1VRpx/Kln82x7Zt8+wz3b9gIZ/yn6T8N/KaObZdkx1rcAHHk6SWZmiQJM0p+4T9/cAKZnWEjohnMPeQo1BuNpQHHjM3QUTsBJ41R/n3Uf7U/tURccFN+qzRkx7IltfOKrMR+C/z1Oeb2fL/m2f7XN6bLX83Is7PFp39e3p0qf0LOJ4ktTSbJ0lSB6ow5CrAR1JKn87+/XrgOcDrsqAwPU/DvwfuAl40x/7voBwY3h0Rz6HcCfhpwFXAXwAvmFk4pfRQRPwM5XkN7omIv6I8ROllwI8A30+5SRKU+zD8CvDWiNgOfAl4UnbMQ1m9Zvu7bJ/3RMSdwLeBb6WU3jnfDyCl9IGIeDHw74DPRsRHKM/T8JKsLn+SUnr/fPtLUrsxNEhSZ5pvyFUof5r/aTh/Q/8s4HeAF1JuXvQF4NVZuQtCQzbE6I/P2KcEHAWeCQwwKzRk+/xlFkr2Ug4pOyg3c/o85bkipsv9S0RcTfnT/q2Uh1f9PPAfKQ+DekFoSCn9dUTsAV4F/BLl/hangXlDQ+Y/UB4p6UZgd7bufuC/89jJ5iSp7UVKqdF1kCRJktTE7NMgSZIkqSJDgyRJkqSKDA2SJEmSKjI0SJIkSarI0CBJkiSpIkODJEmSpIoMDZIkSZIqMjRIkiRJqsjQIEmSJKkiQ4MkSZKkigwNkiRJkioyNEiSJEmqyNAgSZIkqSJDgyRJkqSKDA2SJEmSKjI0SJIkSarI0CBJkiSpov8H1suZrjO+/4IAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 261, "width": 390 }, "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sns.scatterplot(data = df_income, x = \"Education\", y = \"Income\")" ] }, { "cell_type": "markdown", "id": "0f36cba6", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Building a regression model\n", "\n", "We can build a regression model using [`statsmodels.formula.api.ols`](https://www.statsmodels.org/dev/generated/statsmodels.formula.api.ols.html), here imported as `smf`.\n", "\n", "```python\n", "smf.ols(data = df_name, formula = \"Y ~ X\").fit()\n", "```" ] }, { "cell_type": "code", "execution_count": 5, "id": "0e32d29c", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "statsmodels.regression.linear_model.RegressionResultsWrapper" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mod = smf.ols(data = df_income, formula = \"Income ~ Education\").fit()\n", "type(mod)" ] }, { "cell_type": "markdown", "id": "1a3835c7", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Inspecting output" ] }, { "cell_type": "code", "execution_count": 6, "id": "4d4c89ce", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: Income R-squared: 0.812
Model: OLS Adj. R-squared: 0.805
Method: Least Squares F-statistic: 120.8
Date: Fri, 10 Feb 2023 Prob (F-statistic): 1.15e-11
Time: 13:09:06 Log-Likelihood: -115.90
No. Observations: 30 AIC: 235.8
Df Residuals: 28 BIC: 238.6
Df Model: 1
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
Intercept -41.9166 9.769 -4.291 0.000 -61.927 -21.906
Education 6.3872 0.581 10.990 0.000 5.197 7.578
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 0.561 Durbin-Watson: 2.194
Prob(Omnibus): 0.756 Jarque-Bera (JB): 0.652
Skew: 0.140 Prob(JB): 0.722
Kurtosis: 2.335 Cond. No. 75.7


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: Income R-squared: 0.812\n", "Model: OLS Adj. R-squared: 0.805\n", "Method: Least Squares F-statistic: 120.8\n", "Date: Fri, 10 Feb 2023 Prob (F-statistic): 1.15e-11\n", "Time: 13:09:06 Log-Likelihood: -115.90\n", "No. Observations: 30 AIC: 235.8\n", "Df Residuals: 28 BIC: 238.6\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept -41.9166 9.769 -4.291 0.000 -61.927 -21.906\n", "Education 6.3872 0.581 10.990 0.000 5.197 7.578\n", "==============================================================================\n", "Omnibus: 0.561 Durbin-Watson: 2.194\n", "Prob(Omnibus): 0.756 Jarque-Bera (JB): 0.652\n", "Skew: 0.140 Prob(JB): 0.722\n", "Kurtosis: 2.335 Cond. No. 75.7\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mod.summary() ### Lots of stuff here!" ] }, { "cell_type": "markdown", "id": "185b467d", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Extracing specific information" ] }, { "cell_type": "code", "execution_count": 7, "id": "fbdddaa9", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Intercept -41.916612\n", "Education 6.387161\n", "dtype: float64" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## coefficients\n", "mod.params" ] }, { "cell_type": "code", "execution_count": 8, "id": "8840119c", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Intercept 9.768949\n", "Education 0.581172\n", "dtype: float64" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## standard errors\n", "mod.bse" ] }, { "cell_type": "code", "execution_count": 9, "id": "7848964a", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Intercept 1.918257e-04\n", "Education 1.150567e-11\n", "dtype: float64" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## p-values for coefficients\n", "mod.pvalues" ] }, { "cell_type": "markdown", "id": "17f2e935", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Interpreting our model" ] }, { "cell_type": "markdown", "id": "e7565924", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### The linear equation\n", "\n", "Recall that the linear equation is written:\n", "\n", "$Y = \\beta_0 + \\beta_1 * X_1 + \\epsilon$\n", "\n", "How would these terms **map onto** the coefficients below?" ] }, { "cell_type": "code", "execution_count": 10, "id": "7b2367ae", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Intercept -41.916612\n", "Education 6.387161\n", "dtype: float64" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Coefficients\n", "mod.params" ] }, { "cell_type": "markdown", "id": "26d09cb0", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Rewriting the linear equation\n", "\n", "We can *insert* our coefficients into the equation.\n", "\n", "$Y = -41.92 + 6.39 * X_1 + \\epsilon$\n", "\n", "Which can then be used to generate a prediction for a given value of $X$." ] }, { "cell_type": "code", "execution_count": 11, "id": "752bf6bd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "85.88\n" ] } ], "source": [ "x = 20\n", "y = -41.92 + 6.39 * x\n", "print(y) ### our predicted value of Y, given X!" ] }, { "cell_type": "markdown", "id": "9872fdaf", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Check-in\n", "\n", "Write a function called `predict_y(x)` which takes in a value `x` and outputs a prediction based on these learned coefficients." ] }, { "cell_type": "code", "execution_count": 12, "id": "ac1726c5", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "### Your code here" ] }, { "cell_type": "markdown", "id": "4f71e16e", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Solution" ] }, { "cell_type": "code", "execution_count": 13, "id": "9aa89674", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def predict_y(x):\n", " return -41.92 + 6.39 * x" ] }, { "cell_type": "code", "execution_count": 14, "id": "9ca3bf95", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "21.979999999999997" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "predict_y(10)" ] }, { "cell_type": "code", "execution_count": 15, "id": "ebb3176c", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "149.77999999999997" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "predict_y(30)" ] }, { "cell_type": "markdown", "id": "f4f3f750", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Understanding the *intercept*\n", "\n", "> The **intercept** term ($\\beta_0$) is the predicted value of $\\hat{Y}$ when $X = 0$.\n", "\n", "What does this `Intercept` value mean here?" ] }, { "cell_type": "code", "execution_count": 16, "id": "f3a90872", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Intercept -41.916612\n", "Education 6.387161\n", "dtype: float64" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Coefficients\n", "mod.params" ] }, { "cell_type": "markdown", "id": "768fbb77", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### The intercept and the linear equation\n", "\n", "If $X = 0$, the linear equation reduces to:\n", "\n", "$Y = -41.92 + \\epsilon$\n" ] }, { "cell_type": "code", "execution_count": 17, "id": "57321217", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "-41.92" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "predict_y(0) ### predicted value when x = 0" ] }, { "cell_type": "markdown", "id": "332ad301", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Understanding the *slope*\n", "\n", "> For a **continuous** variable, the *slope* is the predicted change in $Y$ for every 1-unit change in $X$.\n", "\n", "What does this *slope* term mean here?" ] }, { "cell_type": "code", "execution_count": 18, "id": "7aae42d1", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Intercept -41.916612\n", "Education 6.387161\n", "dtype: float64" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Coefficients\n", "mod.params" ] }, { "cell_type": "markdown", "id": "617ca5e5", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Check-in: more practice with `statsmodels`\n", "\n", "Build a regression model predicting `Income` from `Seniority`. What are the `params`? What do they mean?" ] }, { "cell_type": "code", "execution_count": 19, "id": "db6ec0f4", "metadata": {}, "outputs": [], "source": [ "### Your code here" ] }, { "cell_type": "markdown", "id": "9e43fb4b", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Solution" ] }, { "cell_type": "code", "execution_count": 20, "id": "86e85ca7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Intercept 39.158326\n", "Seniority 0.251288\n", "dtype: float64" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mod_seniority = smf.ols(data = df_income, formula = \"Income ~ Seniority\").fit()\n", "mod_seniority.params" ] }, { "cell_type": "markdown", "id": "51306431", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Check-in\n", "\n", "Why were *these* parameters chosen as opposed to some other $\\beta_0$ and $\\beta_1$?" ] }, { "cell_type": "code", "execution_count": 21, "id": "a3017795", "metadata": {}, "outputs": [], "source": [ "### Your response here" ] }, { "cell_type": "markdown", "id": "95e9c71f", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Solution\n", "\n", "Linear regression finds the set of parameters $\\beta$ that **minimizes the residual sum of squares (RSS)**.\n", "\n", "I.e., it finds the **line of best fit**!" ] }, { "cell_type": "markdown", "id": "e51d77bd", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Other relevant output\n", "\n", "`statsmodels` also gives us:\n", "\n", "- A **standard error** associated with each coefficient. \n", "- A **t-value** associated with each coefficient. \n", "- A **p-value** associated with each coefficient. \n", "- A **confidence interval** associated with each coefficient." ] }, { "cell_type": "markdown", "id": "097a6e9f", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Standard error\n", "\n", "- These **standard errors** can be used to construct a **confidence interval** around our parameter estimate. \n", "- Assumption: parameters are **estimates**, which have some amount of *normally-distributed* error." ] }, { "cell_type": "code", "execution_count": 22, "id": "464374df", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Intercept 9.768949\n", "Education 0.581172\n", "dtype: float64" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mod.bse ## Standard error associated with each coefficient" ] }, { "cell_type": "markdown", "id": "30b6a93b", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "Can *report* standard error as follows:\n", "\n", "> `Education` was positively related with `Income`, $[\\beta = 6.34, SE = 0.58]$." ] }, { "cell_type": "markdown", "id": "d1aff9d1", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### t-value and p-value\n", "\n", "The **standard error** can also be used to calculate $t$-statistics for each coefficient, which can in turn be used to estimate a $p$-value to test for significance." ] }, { "cell_type": "code", "execution_count": 23, "id": "25e65cc3", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Intercept -4.290801\n", "Education 10.990148\n", "dtype: float64" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## t-values = params / bse\n", "mod.params / mod.bse" ] }, { "cell_type": "code", "execution_count": 24, "id": "d602accb", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Intercept -4.290801\n", "Education 10.990148\n", "dtype: float64" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## double-checking our work\n", "mod.tvalues" ] }, { "cell_type": "code", "execution_count": 25, "id": "801a3b6b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Intercept 1.918257e-04\n", "Education 1.150567e-11\n", "dtype: float64" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Check for significance\n", "mod.pvalues" ] }, { "cell_type": "markdown", "id": "87a6b208", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Check-in\n", "\n", "Knowing what you know about **standard error of the mean**, how do you think our *sample size* affects the standard error for our coefficients?" ] }, { "cell_type": "code", "execution_count": 26, "id": "70d2fb38", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "### Your response here" ] }, { "cell_type": "markdown", "id": "fa19df16", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Solution\n", "\n", "A larger sample ($n$) results in a **lower standard error**.\n", "\n", "*Coming up soon*: We'll discuss how to actually *calculate* standard error for our coefficients." ] }, { "cell_type": "markdown", "id": "049d5e88", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Confidence interval\n", "\n", "Finally, the `conf_int` function `return`s a **confidence interval** for all of our parameters.\n", "\n", "- This is calculated using the **standard error**.\n", "- By default, this assumes $\\alpha = .05$, i.e., a $95\\%$ CI. \n", " - Crucial assumption is that distribution of **sample statistics** ($\\beta$) is **normal**. " ] }, { "cell_type": "code", "execution_count": 27, "id": "56998da8", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01
Intercept-61.927397-21.905827
Education5.1966857.577637
\n", "
" ], "text/plain": [ " 0 1\n", "Intercept -61.927397 -21.905827\n", "Education 5.196685 7.577637" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## alpha = .05\n", "mod.conf_int()" ] }, { "cell_type": "code", "execution_count": 28, "id": "5329fa43", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01
Intercept-68.910782-14.922442
Education4.7812327.993091
\n", "
" ], "text/plain": [ " 0 1\n", "Intercept -68.910782 -14.922442\n", "Education 4.781232 7.993091" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## alpha = .01\n", "mod.conf_int(alpha = .01)" ] }, { "cell_type": "markdown", "id": "cd09ad45", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Check-in\n", "\n", "What is the $95\\%$ confidence interval for our coefficients for the model using `Seniority` to predict `Income`?" ] }, { "cell_type": "code", "execution_count": 29, "id": "4098f4f6", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "### Your code here" ] }, { "cell_type": "markdown", "id": "1ca0152c", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Solution\n", "\n", "These could be reported as follows:\n", "\n", "> `Seniority` was positively related to `Income`, $[\\beta = 0.25]$, $95\\%$ $CI = [0.09, 0.41]$." ] }, { "cell_type": "code", "execution_count": 30, "id": "0fc03181", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01
Intercept21.71421256.60244
Seniority0.0907760.41180
\n", "
" ], "text/plain": [ " 0 1\n", "Intercept 21.714212 56.60244\n", "Seniority 0.090776 0.41180" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mod_seniority.conf_int()" ] }, { "cell_type": "markdown", "id": "ae39777b", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Regression with categorical predictors" ] }, { "cell_type": "markdown", "id": "e909d3da", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### What are *categorical* variables?\n", "\n", "> A **categorical** (or **qualitative**) variable takes on one of several discrete values.\n", "\n", "Examples:\n", "\n", "- `spam` or `not spam`. \n", "- `male` or `female`. \n", "- `right-handed` vs. `left-handed`.\n", "- `smoker` vs. `non-smoker`. \n", "- `treatment` vs. `placebo`. \n", "- `Condition A` vs. `Condition B`.\n", "\n", "**Note**: some variables (e.g., `handnedness` or `gender`) are treated either *categorically* or *continuously*." ] }, { "cell_type": "markdown", "id": "064fbb62", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Example dataset: *Stroop task*\n", "\n", "> In psychology, the [Stroop effect](https://en.wikipedia.org/wiki/Stroop_effect) is the delay in reaction time between congruent and incongruent stimuli." ] }, { "cell_type": "code", "execution_count": 31, "id": "58160ddc", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ConditionRT
0Congruent12.079
1Congruent16.791
2Congruent9.564
\n", "
" ], "text/plain": [ " Condition RT\n", "0 Congruent 12.079\n", "1 Congruent 16.791\n", "2 Congruent 9.564" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_stroop = pd.read_csv(\"data/models/stroop.csv\")\n", "df_stroop.head(3)" ] }, { "cell_type": "markdown", "id": "82b9b02c", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Building a linear model\n", "\n", "*How do we interpret these parameters?*" ] }, { "cell_type": "code", "execution_count": 32, "id": "18fc9b2e", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Intercept 14.051125\n", "Condition[T.Incongruent] 7.964792\n", "dtype: float64" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mod_stroop = smf.ols(data = df_stroop, formula = \"RT ~ Condition\").fit()\n", "mod_stroop.params" ] }, { "cell_type": "markdown", "id": "a2073c3c", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Interpreting the slope for categorical variables\n", "\n", "$\\Large Y = \\beta_0 + \\beta_1X_1 + \\epsilon$\n", "\n", "- $\\beta_0 = 14.05$\n", "- $\\beta_1 = 7.96$\n", "\n", "**What is $X_1$ here?**" ] }, { "cell_type": "code", "execution_count": 33, "id": "8ce50fcb", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Intercept 14.051125\n", "Condition[T.Incongruent] 7.964792\n", "dtype: float64" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mod_stroop.params" ] }, { "cell_type": "markdown", "id": "29d2d434", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Using indicator variables\n", "\n", "> A **indicator variable** (alternatively: \"dummy variable\") represents the possible values of a categorical variable with different numbers, e.g., `0` vs. `1` for a binary variable.\n", "\n", "In our case, this variable might take the form:\n", "\n", "- if `Condition == Congruent` --> 0\n", "- if `Condition == Incongruent` --> 1\n", "\n", "This is also called **dummy coding** or **treatment coding**.\n" ] }, { "cell_type": "markdown", "id": "39e316ab", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Interpreting our equation with indicator variables (pt. 1)\n", "\n", "Our equation might thus look like this:\n", "\n", "$Y = 14.05 + 7.96X_{incongruent=1} + \\epsilon$\n", "\n", "- What happens when $X$ is `Congruent`? \n", "- What happens when $X$ is `Incongruent`?" ] }, { "cell_type": "markdown", "id": "ee38dd0d", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### When `Condition == Congruent`\n", "\n", "When `Condition == Congruent`, the equation reduces to the `Intercept`.\n", "\n", "$Y = 14.05 + 7.96 *0 + \\epsilon = 14.05 + \\epsilon$\n", "\n", "**Note**: This is actually *identical* to the `mean` of the `Congruent` condition!" ] }, { "cell_type": "code", "execution_count": 34, "id": "150c9298", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
RT
Condition
Congruent14.051125
Incongruent22.015917
\n", "
" ], "text/plain": [ " RT\n", "Condition \n", "Congruent 14.051125\n", "Incongruent 22.015917" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_stroop.groupby(\"Condition\").mean()" ] }, { "cell_type": "markdown", "id": "dde232d5", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### When `Condition == Incongruent`\n", "\n", "When `Condition == Incongruent`, the equation is the `Intercept` plus the `slope`.\n", "\n", "$Y = 14.05 + 7.96 * 1 + \\epsilon = 14.05 + 7.96 + \\epsilon$\n", "\n", "**Note**: The resulting value is the `mean` of the `Incongruent` condition." ] }, { "cell_type": "code", "execution_count": 35, "id": "3cfd3b8e", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "22.01" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "14.05 + 7.96" ] }, { "cell_type": "code", "execution_count": 36, "id": "bcd5d68a", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
RT
Condition
Congruent14.051125
Incongruent22.015917
\n", "
" ], "text/plain": [ " RT\n", "Condition \n", "Congruent 14.051125\n", "Incongruent 22.015917" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_stroop.groupby(\"Condition\").mean()" ] }, { "cell_type": "markdown", "id": "765d601a", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Check-in\n", "\n", "What does our `slope` parameter reflect then?" ] }, { "cell_type": "markdown", "id": "10cbcf95", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Solution\n", "\n", "In univariate regression with a categorical, binary predictor, the slope reflects the **difference** in means between the two *levels* of that predictor." ] }, { "cell_type": "code", "execution_count": 40, "id": "2f380acb", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvwAAAILCAYAAACUz0guAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABYlAAAWJQFJUiTwAABJx0lEQVR4nO3dd3jd5PnG8fvxjrP3JnuRhAw7ZUOAUnbZlJUEShJmKavlB4WSsltoyyojCYUk7E0ps2VvsDMgZJG99x6O1/v7Q3IxRkrs2JZ95O/nunwp1vtI5zmByLd1XknmnBMAAACAeEqq6QYAAAAAVB8CPwAAABBjBH4AAAAgxgj8AAAAQIwR+AEAAIAYI/ADAAAAMUbgBwAAAGKMwA8AAADEGIEfAAAAiDECPwAAABBjBH4AAAAgxgj8AAAAQIyl1HQDic7MFkhqJGlhDbcCAACAeOssabNzrktFNiLwV16jevXqNevTp0+zmm4EAAAA8TVz5kzt2LGjwtsR+CtvYZ8+fZrl5ubWdB8AAACIsaysLE2ePHlhRbdjDj8AAAAQYwR+AAAAIMYI/AAAAECMEfgBAACAGCPwAwAAADFG4AcAAABijMAPAAAAxBiBHwAAAIgxAj8AAAAQYwR+AAAAIMYI/AAAAECMEfgBAACAGCPwAwAAADGWUIHfzJqb2Ugze9nM5prZDjPbZGafmNkFZrbb92Nmj5qZ87+6R9E3AAAAUFNSarqBCjpd0kOSVkh6X9JiSa0lnSJpvKRjzOx055wL2tjMTpD0a0lbJTWIpGMAAFDrbc8v1LzV236yvlur+spMS7S4BPxYov0fPEfSLyW97pwrLllpZtdL+krSqfLC/4tlNzSzlpLGSXpWUhtJh0bRMAAAqP3mrd6mEx745CfrX7vsIPXv0LgGOgKqTkJN6XHOveece6102PfXr5T0sP/t0JDNx/rLS6upPQAAAKDWSbQz/LtS4C8Lyw6Y2XmSTpJ0snNunZlF2BYAAABQc2IR+M0sRdJw/9u3yox1knSvpCecc69U4jVyQ4Z67+k+AQAAgOqWUFN6duFOSf0kveGce7tkpX/XngnyLtK9vIZ6AwAAAGpMwp/hN7PLJV0taZakYWWGr5R3ce5xzrkNlXkd51xWyOvnShpcmX0DAAAA1SWhz/Cb2aXypuvMkHSYc259qbEekm6T9Jhz7o0aahEAAACoUQkb+M3sCkkPSJouL+yvLFPSV1K6pPNLPWjLmZnTD7fk/N5fd1JUfQMAAABRSsgpPWZ2rbx5+1MlHemcWxtQtlDSoyG7OE7evfifl7TZrwUAAABiJ+ECv5ndKOlmSbmSflF6Gk9pzrmpkkaG7OMDeYH/eufc3OrpFAAAAKh5CRX4zWyEvLBfJOljSZcH3FN/oXPu8YhbAwAAAGqlhAr8krr4y2RJV4TUfCjp8SiaAQAAAGq7hAr8zrkxksZUwX6GVnYfAAAAQCJI2Lv0AAAAANg9Aj8AAAAQYwR+AAAAIMYI/AAAAECMEfgBAACAGCPwAwAAADFG4AcAAABijMAPAAAAxBiBHwAAAIgxAj8AAAAQYwR+AAAAIMYI/AAAAECMEfgBAACAGCPwAwAAADFG4AcAAABijMAPAAAAxBiBHwAAAIgxAj8AAAAQYwR+AAAAIMYI/AAAAECMEfgBAACAGCPwAwAAADFG4AcAAABijMAPAAAAxBiBHwAAAIgxAj8AAAAQYwR+AAAAIMYI/AAAAECMEfgBAACAGCPwAwAAADFG4AcAAABijMAPAAAAxBiBHwAAAIgxAj8AAAAQYwR+AAAAIMYI/AAAAECMEfgBAACAGCPwAwAAADFG4AcAAABijMAPAAAAxBiBHwAAAIgxAj8AAAAQYwR+AAAAIMYI/AAAAECMEfgBAACAGCPwAwAAADFG4AcAAABijMAPAAAAxBiBHwAAAIgxAj8AAAAQYwR+AAAAIMYI/AAAoE5zzmnWys2BY0XORdwNUPVSaroBAACAmuCc0ytTl2n8xwv03fLgwD9qQo4uOLiLzjugszJSkyPuEKgaBH4AAFDnFBQV69oXvtFLU5btsm7N1p26881Zeue7lRo/Yoia1U+LqEOg6jClBwAA1CnOOf2+HGG/tMmLN2rEP7/S9vzCauwMqB4JFfjNrLmZjTSzl81srpntMLNNZvaJmV1gZkll6nuY2bVm9p6ZLTGzfDNbZWavmtlhNfU+AABAzXlx8jK9XIGwX+LbZZv05zdnVUNHQPVKqMAv6XRJ4yTtK+lLSfdIelFSP0njJT1nZlaq/hZJd0pqLekNSX+V9Kmk4yS9Z2aXR9Y5AACocc45jf94/h5v/2zOEm3aUVCFHQHVL9EC/xxJv5TUwTl3jnPuOufcryX1lrRE0qmSTilV/5akwc65vs65C/36UyQdIalA0l1m1jbi9wAAAGpIzqINmrVyyx5vn1dQrBdzl1ZhR0D1S6jA75x7zzn3mnOuuMz6lZIe9r8dWmr94865KQH7+VDSB5LSJB1QbQ0DAIBa5a3pKyu/j+8qvw8gSgkV+Hej5PO18l5NU9F6AACQ4FZtzqv0PlZXwT6AKMXitpxmliJpuP/tW+Wo7yRvWs92SR+V8zVyQ4Z6l2d7AABQ84qKK/8grcIq2AcQpVgEfnkX5vaT9IZz7u1dFZpZuqQnJaVL+r1zbkME/QEAgFqgaRXcR5978SPRJHzg9++0c7WkWZKG7aY2WdIkSQdKelbS3eV9HedcVsg+cyUNLu9+AABAzdm/a3M99eXiSu8DSCQJPYffzC6VdK+kGZIOc86t30VtsqQn5N3a8zlJ5zrn+EwOAIA6pF5qsmz3Zbt09r57VUkvQFQS9gy/mV0h6e+Spks6wjm3ehe1KZKekhf2n5I03DlXFEWfAACg5jnn9OgnC3T7GzNVmbN9h/VqqU7N61dZX0AUEjLwm9m18ubtT5V0pHNu7S5q0+Sd0T9R0kRJ55e9rScAAIivnYVFuuHl6Xq+kvfPb5qZqptO6FtFXQHRSbjAb2Y3SrpZUq6kX+xmGk+6pJckHSvpUUmjCfsAANQda7fu1EWTcpWzqHL36GiamarHzv+ZOrfg7D4ST0IFfjMbIS/sF0n6WNLlZj+ZibfQOfe4/+eH5YX9tZKWSfpjQP0HzrkPqqllAABQQ2Ys36xRE3O0bOOOwPEkk/q3b6zvlm9S4S5OBx7co4VuPrGfuhD2kaASKvBL6uIvkyVdEVLzoaTHy9S3kPTHXez3g0r2BQAAapG3pq/UVc9N1fb84Ev2Gmak6IGzB+vQni21ekue7vvv93oi4O49D54zWMf2b1vd7QLVKqECv3NujKQxFagfWl29AACA2sc5pwfem6u//mdOaE2XFvU1fkS2urVsIElq1TBDvxqyV2Dg79g0s9p6BaKSUIEfAAAgzI78Iv3+xW/02rTloTUH92ihB84arMaZqRF2BtQsAj8AAEh4KzfladTEHH27bFNozfkHdtYfju2jlOSEfgwRUGEEfgAAkNCmLtmo0RNztHrLzsDxlCTTLSf101k/44FZqJsI/AAAIGG9OnWZfvfCN8oPuc1Os/ppeuicwdq3a/OIOwNqDwI/AABIOMXFTne9M1sPfTAvtKZX64YaPyJbHZtx4S3qNgI/AABIKFt3FuqKZ6bqvzNXhdb8vE9r3XPmQDVIJ+oA/CsAAAAJY8n67Ro5IUezV20JrblkaDdd84teSkr6ycM2gTqJwA8AABLCl/PX6aIncrVhe0HgeFpKku46bR+dOLB9xJ0BtRuBHwAA1HpPf7VYN74yXYXFLnC8VcN0jR2erYEdm0TbGJAACPwAAKDWKiwq1q2vz9Tjny0MrdmnQ2ONHZatNo0zomsMSCAEfgAAUCtt2l6gy56erI+/Xxtac8KAdrrrtH2UkZocYWdAYiHwAwCAWmfemq0aOSFHC9ZuC6255hc9delh3WXGxbnArhD4AQBArfLhnDW67KnJ2pJXGDiemZasv50xUEf3axNxZ0BiIvADAIBawTmnxz5dqFtfn6GQa3PVvkk9jRuerb3bNYq2OSCBEfgBAECNyy8s1o2vTNezOUtCa4Z0bqqHzs1SiwbpEXYGJD4CPwAAqFFrt+7UxU/k6uuFG0JrzsjuoFtO6qf0FC7OBSqKwA8AAGrMzBWbNXJCjpZt3BE4nmTSH47bW78+sDMX5wJ7iMAPAABqxNvfrdSVz07V9vyiwPGGGSl64OzBOrRny4g7A+KFwA8AACLlnNODH8zTXW/PDq3p3DxT40cMUfdWDSLsDIgnAj8AAIhMXkGRfv/CN/rXtOWhNQd1b6EHzh6kJplpEXYGxBeBHwAARGLlpjyNnpSjb5ZuCq0574DOuuG4PkpJToqwMyDeCPwAAKDaTV2yUaMn5mj1lp2B4ylJpptP7Kez990r4s6A+CPwAwCAavXq1GX6/QvfaGdhceB408xUPXRulvbr2jzizoC6gcAPAACqRXGx01//M1v/eH9eaE3P1g306Igh6tgsM8LOgLqFwA8AAKrc1p2FuvLZqfrPjFWhNT/v00p//9VANcxIjbAzoO4h8AMAgCq1ZP12jZqYo1krt4TWXDy0m675RS8lJ/EwLaC6EfgBAECV+XL+Ol385GSt35YfOJ6WkqQ/n9pfJw/qEHFnQN1F4AcAAFXi2a8X64ZXpqugyAWOt2yYrrHDsjRor6YRdwbUbQR+AABQKYVFxbrtjZl67NOFoTX92zfW2OFZatu4XnSNAZBE4AcAAJWwaUeBLntqsj7+fm1ozXH7tNXdpw1QvbTkCDsDUILADwAA9sj8NVs1ckKO5q/dFlpz9ZE9ddnh3WXGxblATSHwAwCACvtozhpd9tRkbc4rDByvl5qsv/9qgI7u1zbizgCUReAHAADl5pzT458t1C3/nqHi4Gtz1b5JPY0dnqW+7RpH2xyAQAR+AABQLvmFxbrpX9P19FdLQmuyOzXVw8Oy1KJBeoSdAdgVAj8AANitdVt36uInJ+urBetDa07P6qBbT+6n9BQuzgVqEwI/AADYpVkrN2vkhBwt3bAjcDzJpOuP7aMLDurCxblALUTgBwAAod75bqWufHaqtuUXBY43TE/RfWcP0mG9WkXcGYDyIvADAICfcM7pwQ/m6e53ZsuFXJzbuXmmxo/IVvdWDaNtDkCFEPgBAMCP5BUU6doXv9GrU5eH1hzYvbn+cfZgNclMi7AzAHuCwA8AAP5n1eY8jZ6Yo2lLN4XWDN+/k248fm+lJidF2BmAPUXgBwAAkqRvlm7UqIk5WrV5Z+B4SpLpTyf21Tn7doq4MwCVQeAHAAD617Tl+t3z07SzsDhwvElmqh46J0v7d2secWcAKovADwBAHVZc7PT3/87R/e/NDa3p2bqBxg8for2aZ0bYGYCqQuAHAKCO2razUFc9N1Vvf7cqtOaI3q10z5kD1TAjNcLOAFQlAj8AAHXQ0g3bNXJCjmat3BJac+GhXfX7o3orOYmHaQGJjMAPAEAd8/XC9bpoUq7WbcsPHE9LTtKdp/bXKYM7RNwZgOpA4AcAoA559uvFuuGV6SooCn6aVosG6Ro7PEuD92oacWcAqguBHwCAOqCwqFh3vDlLj36yILSmX/tGGjssW+2a1IuwMwDVjcAPAEDMbdpRoN88PUUfzVkTWnNc/7a6+/QBqpeWHGFnAKJA4AcAIMYWrN2mCyZ8rflrtoXWXPnznrr8iO4y4+JcII4I/AAAxNTH36/RpU9O1ua8wsDxeqnJ+tsZA3RM/7YRdwYgSgR+AABixjmnCZ8t1C2vz1RRcfDFue0aZ2jciGz1bdc44u4ARI3ADwBAjOQXFuumf32np79aHFozeK8memRYtlo2TI+wMwA1hcAPAEBMrN+Wr4ufyNWXC9aH1pw6uINuP6Wf0lO4OLe0bq3q67XLDgpcDyS6hAr8ZtZc0smSjpPUX1J7SfmSvpX0mKTHnHPFAdsdIOkGSftJypA0V9I/Jd3vnCuKpnsAAKrP7JVbNHLi11qyfkfguJl0/TF9NPLgLlycGyAzLUX9OzC9CfGUUIFf0umSHpK0QtL7khZLai3pFEnjJR1jZqc75/43YdHMTpT0oqQ8Sc9KWi/pBEl/l3Sgv08AABLWf2as0hXPTNG2/OBzWA3TU3TfWYN0WO9WEXcGoDZItMA/R9IvJb1e+ky+mV0v6StJp8oL/y/66xtJGiepSNJQ51yOv/5GSe9JOs3MznTOPRPpuwAAoAo45/TQh/N019uz5YKvzVWn5pkaPzxbPVo3jLY5ALVGUk03UBHOufecc6+VnbbjnFsp6WH/26Glhk6T1FLSMyVh36/PkzfFR5Iurr6OAQCoHnkFRbry2an6y1vhYf+Abs31yiUHEvaBOi7RzvDvSoG/LH2z4cP95VsB9R9J2i7pADNLd87trM7mAACoKqs352n0pFxNXbIxtGbYfp30xxP2VmpyQp3bA1ANYhH4zSxF0nD/29Lhvpe/nFN2G+dcoZktkNRXUldJM3fzGrkhQ70r1i0AAHvu26WbNGpijlZuzgscT04yjfllXw3br1PEnQGorWIR+CXdKamfpDecc2+XWl9yuf2mkO1K1jeppr4AAKgyr01brt+9ME15BT+5IZ0kqUlmqh48Z7AO6NYi4s4A1GYJH/jN7HJJV0uaJWlYRTf3lyGzH3/gnMsKef1cSYMr+LoAAJRbcbHTPf+do/vemxta071VAz06IludmnPfeAA/ltCB38wulXSvpBmSjnDOlX3SSMkZ/LAb6zYqUwcAQK2yPb9QVz07TW99tzK05vDerXTvmQPVMCM1ws4AJIqEvZLHzK6Q9ICk6ZIO8+/UU9Zsf9kzYPsUSV3kXeQ7v5raBABgjy3dsF2nPvT5LsP+hYd01bjh2YR9AKESMvCb2bXyHpw1VV7YXx1S+p6/PDpg7BBJmZI+4w49AIDaJmfhep30j081c8XmwPG05CTdffoAXXdsHyUn8eRcAOESLvD7D826U1KuvGk8a3dR/oKktZLONLPsUvvIkHSr/+1D1dUrAAB74rmcJTpr3BdauzU/cLxFg3Q9PXo/nZbVIeLOACSihJrDb2YjJN0s78m5H0u63OwnZzUWOucelyTn3GYzGyUv+H9gZs9IWi/vab29/PXPRtM9AAC7VlTsdMcbMzX+kwWhNX3bNdK44dlq16RehJ0BSGQJFfjlzbmXpGRJV4TUfCjp8ZJvnHOvmNmhkv4g6VRJGZLmSrpK0n3OhT2fEACA6GzOK9BvnpqiD+esCa05tn8b3X36AGWmJdqPbwA1KaGOGM65MZLG7MF2n0o6tqr7AQCgKixYu00jJ3yteWu2hdZc8fMeuvzwHkpivj6ACkqowA8AQNx8OnetLnlysjbtKAgcz0hN0l9PH6jj9mkbcWcA4oLADwBADXDOadIXi/Sn12aoqDh4dmnbxhkaNzxb/dqHPU4GAHaPwA8AQMQKiop107++01NfLg6tGbRXEz0yLEutGmZE2BmAOCLwAwAQofXb8nXxE7n6ckHZh8P/4JTB7XX7yf2VkZocYWcA4orADwBAROas2qILJnytJet3BI6bSdcd01ujDu6qgNtOA8AeIfADABCBd2eu0m+fmaqtOwsDxxukp+jeMwfqiD6tI+4MQNwR+AEAqEbOOT3y0Xz9+a1ZCnvyy17NMjV+RLZ6tm4YbXMA6gQCPwAA1SSvoEjXv/StXpqyLLRmv67N9NA5WWpaPy3CzgDUJQR+AACqwerNeRo9KVdTl2wMrTl3v7100wl9lZqcFF1jAOocAj8AAFVs+rJNGjUxRys25QWOJyeZxpywt4bt3znaxgDUSQR+AACq0OvfrNDVz09VXkFx4Hjjeql68JzBOrB7i4g7A1BXEfgBAKgCxcVO97z7ve579/vQmm4t6+vREUPUuUX9CDsDUNcR+AEAqKTt+YW6+rlpenP6ytCaob1a6r6zBqlRRmqEnQEAgR8AgEpZtnGHRk3I0YwVm0NrRh/SVdce3VvJSTxMC0D0CPwAAOyh3EXrdeGkXK3dmh84npacpNtO7qfTsztG3BkA/IDADwDAHnghd6muf+lb5RcFX5zbokGaHhmWpaxOzSLuDAB+jMAPAEAFFBU73fnmTI37eEFoTZ+2jTR+RLbaN6kXYWcAEIzADwBAOW3OK9Bvn56i92evCa05pl8b/fWMAcpM40csgNqBoxEAAOWwcO02jZyYo7mrt4bWXH5ED11xRA8lcXEugFqEwA8AwG58NnetLn5ysjbtKAgcz0hN0t2nD9Dx+7SLuDMA2D0CPwAAuzDp84Ua89oMFRW7wPE2jTI0bni2+ndoHHFnAFA+BH4AAAIUFBXrT699pye+WBxaM7BjE40dlqVWjTIi7AwAKobADwBAGRu25euSJyfr8/nrQmtOGdRet5/SXxmpyRF2BgAVR+AHAKCU71dt0QUTcrR4/fbAcTPp2qN768JDusqMi3MB1H4EfgAAfO/NWqXLn56qrTsLA8frpyXrvrMG6Yg+rSPuDAD2HIEfAFDnOec07uP5uuPNWXLB1+aqY7N6Gj98iHq1aRhtcwBQSQR+AECdlldQpOtf/lYvTV4WWrNvl2Z66NwsNaufFmFnAFA1CPwAgDpr9ZY8XTgpV1MWbwytOXvfvTTmhL5KS0mKrjEAqEIEfgBAnTR92SaNmpijFZvyAseTk0w3nbC3hu3XiYtzASQ0Aj8AoM5549sVuuq5qcorKA4cb1wvVQ+eM1gHdm8RcWcAUPUI/ACAOqO42Om+977XPf/9PrSma8v6enTEEHVpUT/CzgCg+hD4AQB1wvb8Ql3z/DS98e3K0JpDe7bUfWcNUuN6qRF2BgDVi8APAIi95Rt3aNTEHH23fHNozciDuui6Y/soOYn5+gDihcAPAIi13EUbdOGkXK3dujNwPDXZdNvJ/XVGdseIOwOAaBD4AQCx9WLuUl330rfKLwq+OLdFgzQ9fG6Wsjs3i7gzAIgOgR8AEDtFxU5/eWuWHvlofmhNn7aNNG54ljo0zYywMwCIHoEfABArW/IK9Ntnpuq9WatDa47q21p/O2Og6qfzYxBA/HGkAwDExqJ12zRyQo6+X701tObyw7vrip/3VBIX5wKoIwj8AIBY+GzeWl3y5GRt3F4QOJ6ekqS7Tx+gEwa0i7gzAKhZSZXdgZm9Z2bDq6IZAAD2xBNfLNLwR78KDfttGmXohYsOIOwDqJOq4gz/UEkfVMF+AACokIKiYt382gxN+mJRaM2Ajk00bliWWjXKiLAzAKg9mNIDAEhIG7fn65InJ+uzeetCa04a2E53nrqPMlKTI+wMAGoXAj8AIOHMXb1FF0zI0aJ12wPHzaTfHdVLFx/aTWZcnAugbiPwAwASyvuzVuvyp6doy87CwPH6acm698xB+vnerSPuDABqp6oK/EMreAbFOeduqaLXBgDUAc45jf94gW5/c6acC67p0LSeHh0xRL3aNIy2OQCoxaoq8B8q7+Ld3XGSzF8S+AEA5bKzsEh/eHm6XshdGlrzsy7N9PC5WWpWPy3CzgCg9quqwP+h/wUAQJVas2WnLnoiV7mLNoTWnPWzjvrTL/spLaXSd5sGgNipqsD/gXPu5iraFwAAkqTpyzZp9MQcLd+UFzienGS68bg+GnFAZy7OBYAQXLQLAKiV3vx2ha56bpp2FBQFjjfKSNE/zhmsg3u0jLgzAEgsBH4AQK3inNN9787V3/87J7Sma8v6Gj88W11bNoiwMwBITDUS+M2spXNuTU28NgCg9tqRX6RrXpim179ZEVpzSM+Wuv+sQWpcLzXCzgAgcVVF4J8gaWp5Cs2ssaRrJV0mqVEVvDYAICZWbNqhURNzNH3Z5tCaXx/YRdcf21spyVycCwDlVenA75w7X5LMrJOkLEkFkr5yzq0qqTGzDElXSrpGUlNJwY9GBADUSZMXb9CFk3K1ZsvOwPHUZNNtJ/XXGUM6RtwZACS+KpnSY2b3SbpE3j32JSnfzK52zj1oZkPlfQrQQdJOSfdKuqMqXhcAkPhemrxU//fSt8ovLA4cb14/TQ8Py9KQzs0i7gwA4qHSgd/MRsibolMsaaa80N9L0n1mtk3SI5KS/eWtzrnllX1NAEDiKyp2+svbs/TIh/NDa3q3aajxI7LVoWlmhJ0BQLxUxSTI8yTlSzrYOdfPOddX0uGSiiQ9KmmlpMHOuUuqIuyb2Wlmdr+ZfWxmm83MmdkTu6hPN7NLzewrM1trZlvNbKaZ3edPQwIARGxLXoFGT8zZZdj/xd6t9eLFBxD2AaCSqmJKzz6SXnbOfV6ywjn3kZm9Iuk0Sb92zn1bBa9T4gZJAyRtlbRUUu+wQjNLkfSupAMlzZL0tLxpRUMk/UbScDM7wDk3owr7AwDswuJ12zVy4teas2praM1lh3XXVUf2VFISD9MCgMqqisDfWNLcgPXf+8vPA8Yq40p5QX+upEMlvb+L2pPlhf13Jf3COfe/CaJm9idJf5R3IfGvq7hHAECAz+et0yVP5mrD9oLA8fSUJP3ltH104sD2EXcGAPFVFVN6kuTdmaesAklyzu2ogtf4H+fc+865751zrhzlXf3l66XDvu9Vf8kjGgEgAk9+uUjDHv0yNOy3bpSu5y7cn7APAFWsqh68VZ7wXRO+85fHmNm9ZUL/8f7yvxH3BAB1SkFRsW759wxN/HxRaM2ADo01dni2WjfKiLAzAKgbqirwjzGzMUEDZlYUsNo556J4yu/rkl6SdIqkb83sv/IuMM6SdJCk+yU9UJ4dmVluyFDoNQQAUNdt3J6vS5+arE/nrgutOXFgO/351H2UkZocYWcAUHdUVeiu6FVVkVyF5ZxzZnaavLn6N0rau9Twu5Kecs4F/UICAKikuau3auSEr7VwXfizFn93VC9dMrSbzLg4FwCqS1U8abfWPt/cf8LvREnHSLpU3rz97fIu5L1P0kdmdrpz7tXwvXicc1khr5EraXCVNQ0AMfD+7NW6/Kkp2rKzMHA8My1Z9/xqoH7Rt03EnQFA3RPFtJqa9H+STpf0W+fcI6XWv+mf+Z8q78m/uw38AIDdc87p0U8W6PY3Zqo45OquDk3rafyIbPVu0yja5gCgjop74C+5MPcnt+50zk0zs/WSOplZc+dc+ARTAMBu7Sws0h9enq4XcpeG1vysczM9dO5gNW+QHmFnAFC3xT3wl/xE+cmtN80sXVLJ6aX8yDoCgBhau3WnLpyUq9xFG0JrzhzSUTef2E9pKbV2JigAxFLcA//HkvpJut7MPnXO7Sw1Nkbe+//aObelJpoDgDiYsXyzRk3M0bKNwY9dSTLpxuP31nkHdObiXACoAQkX+M3sJEkn+d+WXO21v5k97v95rXPuGv/Pt0k6QdIRkmaZ2VuSdsi7aPdn/p9/W/1dA0A8vTV9pa58dqp2FATf8KxhRor+cfZgHdKTZxwCQE1JuMAvaaCkEWXWddUPT9VdJOkaSXLOLTOzwZKulXScpPPlPRl4haTHJf3ZOTer+lsGgHhxzun+9+bqb/+ZE1rTtUV9jRuRrW4tG0TYGQCgrIQL/M65MfKm45S3fo28XwCu2V0tAGD3duQX6XcvTNO/v1kRWnNwjxZ64KzBapyZGmFnAIAgCRf4AQA1Z+WmPI2amKNvl20KrTn/wM76w7F9lJLMxbkAUBsQ+AEA5TJl8QaNnpSrNVt2Bo6nJptuObGfzvzZXhF3BgDYFQI/AGC3XpmyTL9/8RvlFxYHjjern6aHzhmsfbs2j7gzAMDuEPgBAKGKi53ueme2HvpgXmhN7zYNNW54tjo2y4ywMwBAeRH4AQCBtu4s1BXPTNF/Z64OrTly79b6+68GqkE6P04AoLbiCA0A+Ikl67dr5IQczV4V/lzCSw/rpquP7KWkJB6mBQC1GYEfAPAjX8xfp4ufyNWG7QWB42kpSbrrtH104sD2EXcGANgTBH4AwP889eVi/fHV6SosdoHjrRqma+zwbA3s2CTaxgAAe4zADwBQYVGxbn19ph7/bGFozT4dGmvssGy1aZwRXWMAgEoj8ANAHbdpe4EufWqyPpm7NrTmhAHtdNdp+ygjNTnCzgAAVYHADwB12NzVWzVqYo4WrN0WWvO7o3rpkqHdZMbFuQCQiAj8AFBHfThnjS57arK25BUGjmemJetvZwzU0f3aRNwZAKAqEfgBoI5xzumfny7Uba/PUMi1uWrfpJ7Gj8hWn7aNom0OAFDlCPwAUIfsLCzSja9M13M5S0NrhnRuqofOzVKLBukRdgYAqC4EfgCoI9Zu3amLn8jV1ws3hNb8Krujbjmpn9JSkiLsDABQnQj8AFAHzFyxWSMn5GjZxh2B40km3XDc3jr/wM5cnAsAMUPgB4CYe/u7lbry2ananl8UON4wI0UPnD1Yh/ZsGXFnAIAoEPgBIKacc/rH+3N19ztzQmu6tKiv8SOy1a1lgwg7AwBEicAPADGUV1Ck373wjV6btjy05qDuLfSPswercWZqhJ0BAKJG4AeAmFm5KU+jJ+Xom6WbQmvOO6Czbjiuj1KSuTgXAOKOwA8AMTJ1yUaNnpij1Vt2Bo6nJJluPrGfzt53r4g7AwDUFAI/AMTEq1OX6XcvfKP8wuLA8aaZqXro3Czt17V5xJ0BAGoSgR8AElxxsdPd78zWgx/MC63p1bqhxo/IVsdmmRF2BgCoDQj8AJDAtu4s1JXPTtV/ZqwKrfl5n1a658xBapDOIR8A6iKO/gCQoJas366RE3I0e9WW0JpLhnbTNb/opaQkHqYFAHUVgR8AEtCX89fp4icna/22/MDxtJQk/eXUfXTSoPYRdwYAqG0I/ACQYJ75arFufHW6Copc4HjLhukaOyxLg/ZqGnFnAIDaiMAPAAmisKhYt70xU499ujC0pn/7xho7PEttG9eLrjEAQK1G4AeABLBpe4Eue3qyPv5+bWjN8fu01V2nDVC9tOQIOwMA1HYEfgCo5eat2apRE3I0f+220Jqrj+ypyw7vLjMuzgUA/BiBHwBqsY/mrNGlT03WlrzCwPF6qcn6+68G6uh+bSLuDACQKAj8AFALOef02KcLdevrM1QcfG2u2jepp7HDs9S3XeNomwMAJBQCPwDUMvmFxfrjq9P1zNdLQmuyOzXVw8Oy1KJBeoSdAQASEYEfAGqRdVt36uInJuurhetDa07P6qBbT+6n9BQuzgUA7B6BHwBqiZkrNmvkhBwt27gjcDzJpOuP7aMLDurCxbkAgHIj8ANALfDOdyt1xbNTtT2/KHC8YXqK7j97kIb2ahVxZwCAREfgB4Aa5JzTgx/M093vzJYLuTi3c/NMjR8xRN1bNYi2OQBALBD4AaCG5BUU6doXv9GrU5eH1hzYvbn+cfZgNclMi7AzAECcEPgBoAas2pyn0RNzNG3pptCaEft30g3H763U5KQIOwMAxA2BHwAiNm3JRo2elKNVm3cGjqckmf50Yl+ds2+niDsDAMQRgR8AIvTq1GX6/QvfaGdhceB408xUPXhOlvbv1jzizgAAcUXgB4AIFBc7/e0/c/TA+3NDa3q2bqDxw4dor+aZEXYGAIg7Aj8AVLNtOwt15bNT9c6MVaE1R/RupXvOHKiGGakRdgYAqAsI/ABQjZas365RE3M0a+WW0JqLDu2m3x3VS8lJPEwLAFD1CPwAUE2+WrBeFz2Rq/Xb8gPH01KS9OdT++vkQR0i7gwAUJcQ+AGgGjz79WLd8Mp0FRQFP02rZcN0jR2WpUF7NY24MwBAXUPgB4AqVFhUrNvfmKV/frogtKZf+0YaNzxbbRvXi7AzAEBdReAHgCqyaUeBfvP0FH00Z01ozXH7tNXdpw1QvbTkCDsDANRlBH4AqALz12zVyIk5mr9mW2jNVUf21G8O7y4zLs4FAESHwA8AlfTx92t06ZOTtTmvMHC8Xmqy/nbGAB3Tv23EnQEAQOAHgD3mnNPjny3Ura/PVFFx8MW57RpnaNyIbPVt1zji7gAA8BD4AWAP5BcW66Z/TdfTXy0Jrcnq1FQPn5ullg3TI+wMAIAfI/ADQAWt35avi57I1VcL1ofWnJbVQbed3E/pKVycCwCoWQR+AKiAWSs3a+SEHC3dsCNwPMmk647po5EHd+HiXABArZBU0w1UlJmdZmb3m9nHZrbZzJyZPbGbbczMRpjZB2a23sx2mNkCM3vOzHpG1TuAxPafGat06oOfhYb9hukpevS8IRp1SFfCPgCg1kjEM/w3SBogaaukpZJ676rYzDIkPS/peEmzJT0laYukdpIOltRT0pxq7BdAgnPO6cEP5unud2bLBV+bq07NM/XoiGx1b9Uw2uYAANiNRAz8V8oL+nMlHSrp/d3U/1Ve2L9D0g3OueLSg2aWWh1NAoiHvIIi/d+L3+iVqctDaw7o1lwPnjNYTTLTIuwMAIDySbjA75z7X8Df3UfmZtZN0kWSvpb0B+d+em7OOVdQ1T0CiIfVm/M0alKupi3ZGFozfP9OuvH4vZWanHAzJAEAdUTCBf4KOkvedQoTJDUysxMkdZS0TtJ7zrm5NdkcgNrrm6UbNXpirlZuzgscT0kyjfllX527X6eIOwMAoGLiHviH+MvGkuZJal5qzJnZQ5Iud84V7W5HZpYbMrTLawgAJJ7Xpi3XNc9P087C4sDxJpmpevCcwTqgW4uIOwMAoOLi/hl0K395s6QcSf0lNZR0hLxfAC6RdGPNtAagtikudvrrO7P1m6enhIb9Hq0a6NVLDyTsAwASRtzP8Jc88WaFpJOdcyX30nvPzE6TNFnSVWZ2u3Muf1c7cs5lBa33z/wPrqqGAdSMbTsLddVzU/X2d6tCaw7v3Ur3njlQDTO41h8AkDjiHvg3+Mu3SoV9SZJzbpqZLZDUTVIfSdOibg5A7bB0w3aNnJCjWSu3hNZceEhX/f7o3kpO4v76AIDEEvfAP1vSLyRtDBkv+YWgXiTdAKh1chau14WTcrVuW/CHfGnJSbrjlP46NatDxJ0BAFA14j6H/11/2a/sgJmlS+rhf7swqoYA1B7Pfb1EZ437IjTst2iQrqdH70fYBwAktLif4X9T0nxJR5nZkc65/5Qau1He3Xs+dM6trJHuANSIwqJi3fHmLD36yYLQmr7tGmnc8Gy1a8IHgACAxJZwgd/MTpJ0kv9tG3+5v5k97v95rXPuGklyzuWb2QhJ70h608xelrRI3u06D5G0RtLoaDoHUBtszivQb56aog/nrAmtObZ/G919+gBlpiXcIRIAgJ9IxJ9mAyWNKLOuq/8leYH+mpIB59wnZpYt6SZJh0lqImmVpLGSbnHOLa3mfgHUEgvWbtPICV9r3pptoTVX/LyHLj+8h5K4OBcAEBMJF/idc2MkjangNjMk/ao6+gGQGD75fq0ufWqyNu0oCBzPSE3S384YqGP7t424MwAAqlfCBX4AqAjnnCZ+vkg3/3uGiopdYE27xhkaOzxb/do3jrg7AACqH4EfQGzlFxbrpn99p6e/WhxaM3ivJnp4WJZaNcyIsDMAAKJD4AcQS+u35eviJ3L15YL1oTWnDG6vO07pr/SU5NAaAAASHYEfQOzMWbVFF0z4WkvW7wgcN5OuO6a3Rh3cVWZcnAsAiDcCP4BY+e+MVfrtM1O0Lb8ocLxBeoruO2ugDu/dOuLOAACoGQR+IGa25xdq3uqf3nayW6v6sb6vvHNOD384X395e5Zc8LW56tQ8U+OHZ6tH64bRNgcAQA2K709/oI6at3qbTnjgk5+sf+2yg9S/QzzvQpNXUKTrXvpWL09ZFlqzf9fmevCcwWpaPy3CzgAAqHkEfgAJbfXmPI2elKupSzaG1py731666YS+Sk1Oiq4xAABqCQI/gIT17dJNGjUxRys35wWOJyeZxpywt4bt3znaxgAAqEUI/AAS0r+/Wa5rnp+mvILiwPHG9VL10DmDdUD3FhF3BgBA7ULgB5BQioud7vnvHN333tzQmu6tGmj88Gx1blE/ws4AAKidCPwAEsb2/EJd9ew0vfXdytCaw3q11L1nDVKjjNQIOwMAoPYi8ANICMs27tDICTmauWJzaM3oQ7rq2qN7KzmJh2kBAFCCwA+g1stdtF4XTsrV2q35geNpyUm6/ZT+Oi2rQ8SdAQBQ+xH4AdRqz+cs0R9enq78ouCLc1s0SNMjw7KU1alZxJ0BAJAYCPwAaqWiYqc735ypcR8vCK3Zu20jjRuRrfZN6kXYGQAAiYXAD6DW2ZxXoMufnqIPZq8JrTmmXxv99YwBykzjMAYAwK7wkxJArbJw7TZdMOFrzVuzLbTmt0f00G+P6KEkLs4FAGC3CPwAao3P5q7VxU9O1qYdBYHjGalJ+uvpA3XcPm0j7gwAgMRF4AdQK0z6fKHGvDZDRcUucLxt4wyNG56tfu0bR9wZAACJjcAPoEYVFBVrzL++05NfLg6tGbRXEz0yLEutGmZE2BkAAPFA4AdQYzZsy9fFT+bqi/nrQ2tOGdRet5/SXxmpyRF2BgBAfBD4AdSIOau2aOSEHC1evz1w3Ez6v6N7a/QhXWXGxbkAAOwpAj+AyL03a5Uuf3qqtu4sDByvn5as+84apCP6tI64MwAA4ofADyAyzjmN/Wi+7nxrllzwtbnaq1mmxo/IVs/WDaNtDgCAmCLwA4hEXkGRrn/pW700ZVlozX5dm+nBc7LUrH5ahJ0BABBvBH4A1W71ljxdOClXUxZvDK05Z9+9NOaXfZWanBRdYwAA1AEEfgDVavqyTRo1MUcrNuUFjicnmcacsLeG7d852sYAAKgjCPwAqs3r36zQ1c9PVV5BceB443qpevCcwTqwe4uIOwMAoO4g8AOocsXFTve++73ufff70JpuLetr/Igh6tKifoSdAQBQ9xD4AVSp7fmFuvq5aXpz+srQmkN7ttT9Zw9So4zUCDsDAKBuIvADqDLLNu7QqAk5mrFic2jNqIO76P+O6aPkJB6mBQBAFAj8AKpE7qL1unBSrtZuzQ8cT0023XZyf52R3THizgAAqNsI/AAq7YXcpbr+pW+VXxR8cW6LBml6+NwsZXduFnFnAACAwA9gjxUVO/35rVka+9H80Jo+bRtp3PAsdWiaGWFnAACgBIEfwB7Zklegy5+eovdnrwmtObpvG/31jAGqn86hBgCAmsJPYQAVtnDtNo2cmKO5q7eG1lx+RA9dcUQPJXFxLgAANYrAD6BCPpu7Vpc8NVkbtxcEjmekJunu0wfo+H3aRdwZAAAIQuAHUG6TvlikMf/6TkXFLnC8TaMMjRuerf4dGkfcGQAACEPgB7BbBUXFuvm1GZr0xaLQmgEdm2jcsCy1apQRYWcAAGB3CPwAdmnDtnxd8uRkfT5/XWjNyYPa645T+isjNTnCzgAAQHkQ+AGE+n7VFo2cmKNF67YHjptJvz+qty46tKvMuDgXAIDaiMAPIND7s1brN09P0dadhYHj9dOSde+Zg/TzvVtH3BkAAKgIAj+AH3HOadzH83XHm7Pkgq/NVcdm9TR++BD1atMw2uYAAECFEfgB/M/OwiJd/9J0vTh5aWjNvl2a6aFzs9SsflqEnQEAgD1F4AcgSVq9JU8XTcrV5MUbQ2vO+tle+tMv+yotJSm6xgAAQKUQ+AFo+rJNGj0xR8s35QWOJyeZ/nj83hq+fycuzgUAIMEQ+IE67o1vV+jq56ZpR0FR4HijjBQ9eE6WDurRIuLOAABAVSDwA3WUc073vTtXf//vnNCari3r69ERQ9SlRf0IOwMAAFWJwA/UQTvyi3TN89P0+rcrQmsO6dlS9581SI3rpUbYGQAAqGoEfiBG1mzZqWdzFgeOfT5/rfq0bajVW3Zq1MQcfbd8c+h+Ljioi647prdSkrk4FwCAREfgB2Igr6BIt/x7hp7LWaKCouCb59/+xiyN/XC+8ouKtTkv+GFaqcmm207qrzOGdKzOdgEAQIQI/ECC25JXoPMe+1q5izbstnbttvzQseb10/TwsCwN6dysKtsDAAA1jMAPJLDiYqdLn5pSrrC/K73bNNT4Ednq0DSzijoDAAC1BYEfSGD//naFPpqzplL7+MXerfX3Xw1U/XQOBwAAxFHCXZFnZqeZ2f1m9rGZbTYzZ2ZPVGD7R/1tnJl1r85egeo28bOFldr+yD6t9PC5WYR9AABiLBF/yt8gaYCkrZKWSupd3g3N7ARJv/a3bVAt3QERmb1yi3IqOZVn/fYCJSXx5FwAAOIs4c7wS7pSUk9JjSRdXN6NzKylpHGSnpWUWz2tAdH5euH6Su9jyuINKiwqroJuAABAbZVwgd85975z7nvnXPC9B8ON9ZeXVnVPQE3YtKOg0vsodgq9RScAAIiHRJzSU2Fmdp6kkySd7JxbZ8YUBiS+1OSq+f84LSXhfu8HAAAVEPvAb2adJN0r6Qnn3CuV2E/YNKByX0MAVKW2jetVeh8NM1JUPy25CroBAAC1VaxP7ZlZkqQJ8i7SvbyG2wGq1NBeLZVZybB+/D5txSdeAADEW9zP8F8p6VBJxznnKnU7E+dcVtB6/8z/4MrsG9gTDTNSdfKg9nryy8V7vI9z9+tUhR0BAIDaKLZn+M2sh6TbJD3mnHujpvsBqsOvD+qyx3P5D+reQn3bNa7ijgAAQG0T28Avqa+kdEnnl3rQljMzJ++svyR97687qca6BCqhW8sGuuOUfSq8Xfsm9fS3Xw2oho4AAEBtE+cpPQslPRoydpykNpKel7TZrwUS0mlZHVTsnK5/6VsVFu/+brXdWzXQY+cNUauGGRF0BwAAalpsA79zbqqkkUFjZvaBvMB/vXNuboRtAdXijOyO6teuscZ+NE///mZFYPBv2SBd5x3YWecd0Fn102P7Tx8AAJSRcD/1/ek3J/nftvGX+5vZ4/6f1zrnrom4LaDG7d2uke45c5BOzeqgYY9+9ZPxcSOyNbBjk+gbAwAANSrhAr+kgZJGlFnX1f+SpEWSCPyos5rUSwtcn8ztNwEAqJMS7qJd59wY55zt4qtzOfYx1K9lOg8AAABiLeECPwAAAIDyI/ADAAAAMUbgBwAAAGKMwA8AAADEGIEfAAAAiDECPwAAABBjBH4AAAAgxgj8AAAAQIwR+AEAAIAYI/ADAAAAMUbgBwAAAGKMwA8AAADEGIEfAAAAiDECPwAAABBjBH4AAAAgxgj8AAAAQIwR+AEAAIAYI/ADAAAAMUbgBwAAAGKMwA8AAADEGIEfAAAAiDECPwAAABBjBH4AAAAgxgj8AAAAQIwR+AEAAIAYI/ADAAAAMUbgBwAAAGKMwA8AAADEGIEfAAAAiDECPwAAABBjBH4AAAAgxgj8AAAAQIwR+AEAAIAYI/ADAAAAMUbgBwAAAGKMwA8AAADEGIEfAAAAiDECPwAAABBjBH4AAAAgxgj8AAAAQIyl1HQDAKpWt1b19dplBwWuBwAAdQ+BH4iZzLQU9e/QuKbbAAAAtQRTegAAAIAYI/ADAAAAMUbgBwAAAGKMwA8AAADEGIEfAAAAiDECPwAAABBjBH4AAAAgxgj8AAAAQIwR+AEAAIAYI/ADAAAAMUbgBwAAAGKMwA8AAADEGIEfAAAAiDECPwAAABBj5pyr6R4Smpmtq1evXrM+ffrUdCsAAACIsZkzZ2rHjh3rnXPNK7Idgb+SzGyBpEaSFtZwK0BZvf3lrBrtAgASC8dO1GadJW12znWpyEYEfiCmzCxXkpxzWTXdCwAkCo6diCPm8AMAAAAxRuAHAAAAYozADwAAAMQYgR8AAACIMQI/AAAAEGPcpQcAAACIMc7wAwAAADFG4AcAAABijMAPAAAAxBiBHwAAAIgxAj8AAAAQYwR+AAAAIMYI/AAAAECMEfiBEGbW28zuN7PpZrbJzPLNbLmZvW5mF5hZRk33GEdm1tnMnJk9XtO9APiB/++Sh/ckMDM7z//veF5N94JopdR0A0BtZGZ/lHSTvF+Kv5A0QdJWSa0lDZU0XtLFkrJrqEUAAIByIfADZZjZ9ZL+JGmJpNOdc18G1Bwv6eqoewMAAKgopvQApZhZZ0ljJBVIOjYo7EuSc+7fko4us+0ZZvaRP/1nh5l9a2bXmVl6wOss9L8yzewuM1tsZjvNbK6ZXWtmFrCNmdlvzWyGmeWZ2TIze8DMGpfsr0z9/z66NbOjzewDvzdX8l53NXXGrw/8+N7MjjKzN8xsrd/3PP99NKnMezWzMZIW+N+OKJlCwEfQQO1U+jji//kZ/7iQZ2Y5/smRsG1/ZWbvmtl6v36hmT1tZtll6tLN7P/M7Bsz225mm83sYzM7o6r68Y+j95jZUr92lpldZWZdg46T/v6dP/4bv7cdZvaBP77LqTP+2AcB61PM7BIz+8J/n9vNbIqZXWZmSWVqK/Re/dd7zP/2sTLH185BfSI+OMMP/Nj5klIlPeOcm76rQufczpI/m9ntkq6TtFbSU/Km/xwj6XZJR5nZkc65gjK7SJX0jqR2kt6UVCjpJEl3SsqQ9ylDaf+QN41ouaSxkvIl/VLSz/x9ld1/idPk/XLypqSHJXXe1fvaHfOmO/1J0npJ/5a0WtI+kq6RdKyZ7e+c21xms/K+1w8kNZH0W0nTJL1Sah9TK9M3gGrVSdJXkuZLmiSpmaRfSXrVzH7unHu/pND/Jf8xSSPkHTNfkrRGUgdJh0maLSnHr02T9LakQyXNkncczJR3XHvWzAY6566vZD8Zkt6TNFjSFElPSmos6Q+SDt7N+77Xr3ld0huSinZTH8rMUiW9JukoeX8HT0nKk/d3cr+kfSUNC9i0vO/1cUkbJZ0o6VX9+Ji6cU/7RoJwzvHFF1/+l6R3JTlJIyuwzf7+NosltSm1PkXewdtJur7MNgv99W9IqldqfSt5B96NklJLrT/Yr58tqUmp9WmSPvLHFpZ5jfP89cWSjg7ou7M//njI+/rAO0T8aN1h/jafle6jzOv9vZLvdZd98cUXXzXz5f+7dGXWlfx7dZJuKjN2VMm//TLrR/vrv5LUuMxYsqS2pb6/rtTxI6XU+lalji0HVLKfG/31T0uyUus7yvtF5CfHI3nh2UlaJqlLwN9VyfHwvF38XX5QZt0Yf/39kpLL/J086o+dWMn3usu++IrvF1N6gB9r6y+XVmCbX/vLW51zK0tWOucK5c3zL5Y0MmTby51zO0pts1remZfGknqVqhvhL29zzm0sVZ8v7wfirrzqnHurHO+jPC73l6NK9+H38ri8M0bnhG1bzvcKIPEsknRr6RXOubflnQj5WZna3/jLC51zm8psU+ScW1Fq1a/lBdSr/GNqSd1qSbf43wYdXyvSzwh5x+nrnHOuVP0SSfcE7Lu0vzjnFuymZrf86TqXSVop6Urn3P8+KfD/fLW8v4eg42tF3ivqKKb0AD9WMp+8IreeG+wv3ys74JybY2ZLJXUxsyZlQvIm59zcgP0t8ZdNS60b5C8/Caj/Qt4UmTBf7WKsovaXN3XodDM7PWA8TVJLM2vunFtXan1F3iuAxDO1dEgtZYm844YkyczqS+onaZVzbsqudmhmDSV1l7TMOTcroKTkmDsoYKy8/TSS1E3SEufcwoD6oGNuaVV1fO0pqbmk7yXdYD+9jEuSdkjqE7C+XO8VdRuBH/ix5ZJ6y5tLWl6N/eWKkPEVkvby6zaWWr8xqFg/hPfkgNdYVbbYOVdkZuvKri9l5S7GKqq5vOPGTbupayCpdE8bQ+qC3iuAxLMxZH2hfnyDkCb+clk59lmeY2vpfe5JP4385U+OrbtZX6Kqjq/N/WUP7fr42iBg3caQ2rLvFXUY/yMAP1ZyNueICmxT8pF0m5DxtmXq9kTJRbCtyw6YWbJ++GERJOzTimJ/GfaLf5OAdZskbXDO2W6+Fu2iHwB110Z/2b4ctTV6bN3N+hIVPr4G3c1MP7yHl3dzbO2ym36AQAR+4Mcekzdl5VQz23tXhfbD7TZLPpYeGlDTXd6nBQvKznmvoJLXOChgbD/t2ad1G/xlx7ID/sfcPQO2+UJSUzPruwevV14lH01z1h+IGefcNknTJbU2s6CpOKVrt0iaJ6m9mfUIKDnMX06uRD+b5d3dpn3IrSmDjrnlEXp8VfADG2fJ+2VoP/9uPdWF42sdReAHSvHncI6RNxf99bL3gy5hZiW3uZSkf/rLG8ysZamaZEl3y/t39mglW5voL/9gZiUfc5fcsu72Pdmh/8N0lqQDS/9y4/f9N0n1Ajb7u78cZ2btyg6aWX0z229P+illg7yzZntVcj8Aaqf7/OUjpY9nknfxqpm1LbXqn/KurbrLPzaV1LWQd3edkprKmCjvOH2H2Y+eC9JR0hV7uM8ceWf5zzazzFL7bCbpL2WL/QuS75f3qcV9ZvaT46+Ztd3diahyKJlqyfG1jmEOP1CGc+52MyuZp/61mX0m7+C9Vd7Hu4fIm2eZ49d/ZmZ/kfR7SdPN7AVJ2+Tdh7+fvGlCd1Wypw/NbKy829l9Z2Yvyvsk4gR5HwUv1w8fIVfEXfJ+GfnUzJ7XD/d8TpV3H/wBZfp418z+T9Idkr43szfkPSirgbx7QR8q7/3+6KFkFeGc22pmX0o62MyelDRH3lmpfznnvtnT/QKoNcbLO3M+XN5x5FV5t79sJ+lweQF+jF97t7xj6YmSpvnHnExJp8u7NedfnHO7u7B2d/4i77kgZ0rqZWbvyLt+4Ax5tz0+SRU8vjrnVvjHr2GSpprZ6/KuFzjW32fQpxu3yDvmXiTpBDN7T961Dq3k/cw5UN6zAWZU7O39yOeStku6wv/lo+QahfvL3jEJ8ULgBwI45272A/Al8gLw+fIeELVO3q0n/yzpiVL115rZFHm3VRsuLzDPk3SDpL/6t8+srIvlnZG/UN4PhHWSXpZ0vbzbiM6r6A6dc//0z2hdJe/WdBvk3SrzekkvhmzzZzP7VN4tOg+S94N4k7wfTGPlPSymsobJ+zThaElnyTvDt1QSgR9IcP6tL0f4wXq0vGCdLu8i3I8l/atUbb6ZHSnvGHW2vFt6Fso7IXGFc+7pKuhnh5kdJulmeQ/0ulLeiYzb/X5O0g9z/StilLxAfZakS+XdJvM+eSdafvKUYOdcgZmdJOlceffLP17eyZQ1fj83ynso2B5zzm0ws1PlndA6X1J9f+gJVe5aCNRyVuqWswASkD+3dY68pwOfVdP9AEBcmNkoeScyLnLOPVLT/QB7ijn8QIIwszb+w1lKr8vUDw+GeTnypgAgBkKuSeoo76x6oaR/R94UUIWY0gMkjisknWVmH8j76LuNvNuHdpB3AfHzNdYZACS2F/274+TKu1tOZ3lTajLlPYG3PM8NAGotpvQACcLMjpB0jaSBkprJO+s0R96c+XuccwU11x0AJC4zu0TetUM95F2wu1Xe7ZAfcM69VJO9AVWBwA8AAADEGHP4AQAAgBgj8AMAAAAxRuAHAAAAYozADwAAAMQYgR8AAACIMQI/AAAAEGMEfgAAACDGCPwAgMiZ2UIzW1hm3Xlm5szsvAruy/lPoAYABCDwA0BMmFlvM7vfzKab2SYzyzez5Wb2upldYGYZNd3jngj65QAAUH48aRcAYsDM/ijpJnkncr6Q9LWkrZJaSxoqqaukXOdcdk31WFpJgHfOdS61rrGktpJWOOc27aq2zL56S9runFtcbQ0DQAJLqekGAACVY2bXS/qTpCWSTnfOfRlQc7ykq6PurSL8kL9pt4U/3W5WNbQDALHBlB4ASGBm1lnSGEkFko4NCvuS5Jz7t6Sjy2x7hpl95E//2WFm35rZdWaWHvA6C/2vTDO7y8wWm9lOM5trZteamQVsY2Z2mZl9Z2Z5ZrbMzB7wz+QHvZcfzeE3s6Fm5iR1ktTJHyv5erzUdoFz+M2ssZndYWaz/dffYGZvm9nPA2qH+vsZY2YD/WlQG81su5l9aGYHBPUMAImAM/wAkNjOl5Qq6Rnn3PRdFTrndpb82cxul3SdpLWSnpI3/ecYSbdLOsrMjnTOFZTZRaqkdyS1k/SmpEJJJ0m6U1KGvE8ZSrtH0uWSVkgaK++XkhMl7SspTVL+bt7bQn+fV5TaX4mpu9rQzJpI+lTS3vKmN90jqYWkMyS9Y2YXO+ceCdg0W9LvJX0uabykvSSdKuldMxvonJu9m54BoNZhDj8AJDAze1fS4ZJGOefGl3Ob/SV9Jm8K0M+ccyv99SmSXpZ0vKQ/OOduL7XNQnln2t+UdKpzboe/vpWkOX5Zy5JfEvwz4p9Kmue/xnp/fYak9yXtJ2lRmTn850l6TNL5zrnHy7z2rubwO0kfOueGllr3iKTR8n7RuMj5P+zMrIekHHm/oPRyzi301w/1+1LA618o6WFJDznnLgnqAQBqM6b0AEBia+svl1Zgm1/7y1tLwr4kOecK5c3zL5Y0MmTby0vCvr/NakmvSmosqVepuvP95W0lYd+vz5P3yUK1MbNUSefK+9TiOlfqzJZz7ntJ98n7hGF4wOaflg77vn/K+zTjZ9XSMABUMwI/ACS2krnzFfm4drC/fK/sgHNujrxfHrr402JK2+ScmxuwvyX+smnAa3wYUP+xvABdXXpLypQ0rfQvG6WUvO9BAWM5ZVf4n1qs0o/fHwAkDAI/ACS25f6yQwW2KblodkXI+IoydSU2htSXhPfkgNdYVbbYOVckad2uW6yU8r6/JgFjG0O2KdSP3x8AJAwCPwAktk/85REV2Kbk1pdtQsbblqnbEyXbti47YGbJkppXYt/lfe3qfH8AkDAI/ACQ2B6Td/ebU81s710Vlrrd5hR/OTSgpru8TwsWOOc2VqKvyf7y0ICxg1Wxu8QVqWJn12dL2i5poJkFTcM5zF9ODhgDgNgh8ANAAvPvMjNG3kWor5tZ4JN0zexoeXfYkbyLUCXpBjNrWaomWdLd8n42PFrJ1h73l38ws2alXiND0h0V3Nc6SS3NrF55ip1z+ZKelNRA0s2lx8ysm7xbhRZImlTBPgAgIXEffgBIcM652/1bat4k6Wsz+0zexadb5U2pOURSye0o5Zz7zMz+Iu9+89PN7AVJ2+Tdh7+fvGlCd1Wyp0/N7H5Jvyn1GiX34d+g8Pn1Qd6VNETSW2b2kaSd8i7IfW0X2/yfvE8SLjOzIfJuuVlyH/6Gki5zzi2o4NsCgIRE4AeAGHDO3Wxmz0u6RN6UlfPl3Wt+nbyHVP1Z0hOl6q81symSLpN3e8pUeffMv0HSX/2z5JX1W3n36L9U0oV+Ly9Lul7StArs51Z5F9ieIOlAedN7JkgKDfzOufX+8wauk3SKpKsk7ZD0laS7nHPvVPC9AEDC4sFbAAAAQIwxhx8AAACIMQI/AAAAEGMEfgAAACDGCPwAAABAjBH4AQAAgBgj8AMAAAAxRuAHAAAAYozADwAAAMQYgR8AAACIMQI/AAAAEGMEfgAAACDGCPwAAABAjBH4AQAAgBgj8AMAAAAxRuAHAAAAYozADwAAAMQYgR8AAACIsf8HicJudsnr+nAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 261, "width": 382 }, "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sns.pointplot(data = df_stroop, x = \"Condition\", y = \"RT\")" ] }, { "cell_type": "markdown", "id": "e065b841", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Check-in\n", "\n", "Load the `tips` dataset from `seaborn` and fit a linear model predicting `tip` from `time`. \n", "\n", "- Inspecting the model summary.\n", "- How should you interpret the resulting coefficients?\n", "\n" ] }, { "cell_type": "code", "execution_count": 41, "id": "0d814367", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "df_tips = sns.load_dataset(\"tips\")\n", "### Your code here" ] }, { "cell_type": "markdown", "id": "4f39a6fe", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Solution\n", "\n", "The intercept reflects the mean `tip` for `Lunch`, and the slope reflects the *difference* in `tip` amount between `Lunch` and `Dinner`." ] }, { "cell_type": "code", "execution_count": 42, "id": "13a6207a", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Intercept 2.728088\n", "time[T.Dinner] 0.374582\n", "dtype: float64" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mod_tips = smf.ols(data = df_tips, formula = \"tip ~ time\").fit()\n", "mod_tips.params" ] }, { "cell_type": "code", "execution_count": 43, "id": "83885766", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
tip
time
Lunch2.728088
Dinner3.102670
\n", "
" ], "text/plain": [ " tip\n", "time \n", "Lunch 2.728088\n", "Dinner 3.102670" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_tips[['time', 'tip']].groupby(\"time\").mean()" ] }, { "cell_type": "markdown", "id": "98c633c6", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Variables with $>2$ levels?\n", "\n", "Many categorical variables have $>2$ **levels**:\n", "\n", "- `Assistant Professor` vs. `Associate Professor` vs. `Full Professor.\n", "- Conditions `A`, `B`, and `C`. \n", "- `Europe` vs. `Asia` vs. `Africa` vs. `Americas` vs. `Oceania`.\n", "\n", "By default, `statsmodels` will use **dummy coding**, choosing the *alphabetically first* level as the **reference level**." ] }, { "cell_type": "markdown", "id": "cb849c37", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Check-in\n", "\n", "- Read in the dataset at `data/housing.csv`. \n", "- Then, build a linear model predicting `median_house_value` from `ocean_proximity`. \n", "- Inspect the model `summary()`. How should you interpret the results?" ] }, { "cell_type": "code", "execution_count": null, "id": "5998452b", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "### Your code here" ] }, { "cell_type": "markdown", "id": "8b48282f", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Solution" ] }, { "cell_type": "code", "execution_count": 51, "id": "8fc82cb1", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "Intercept 240084.285464\n", "ocean_proximity[T.INLAND] -115278.893463\n", "ocean_proximity[T.ISLAND] 140355.714536\n", "ocean_proximity[T.NEAR BAY] 19128.026326\n", "ocean_proximity[T.NEAR OCEAN] 9349.691963\n", "dtype: float64" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_housing = pd.read_csv(\"data/housing.csv\")\n", "mod_housing = smf.ols(data = df_housing, formula = \"median_house_value ~ ocean_proximity\").fit()\n", "mod_housing.params" ] }, { "cell_type": "code", "execution_count": 52, "id": "7f4e918b", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
median_house_value
ocean_proximity
<1H OCEAN240084.285464
INLAND124805.392001
ISLAND380440.000000
NEAR BAY259212.311790
NEAR OCEAN249433.977427
\n", "
" ], "text/plain": [ " median_house_value\n", "ocean_proximity \n", "<1H OCEAN 240084.285464\n", "INLAND 124805.392001\n", "ISLAND 380440.000000\n", "NEAR BAY 259212.311790\n", "NEAR OCEAN 249433.977427" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_housing[['ocean_proximity', 'median_house_value']].groupby(\"ocean_proximity\").mean()" ] }, { "cell_type": "markdown", "id": "d3538fd1", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### How to *code* your variables?\n", "\n", "> **Coding** refers to the approach taken to representing the different levels of a categorical variable in a regression model.\n", "\n", "There are a [variety of different approaches to contrast coding](https://stats.oarc.ucla.edu/r/library/r-library-contrast-coding-systems-for-categorical-variables/), some of which are summarized below.\n", "\n", "|Approach|Description|What is the `Intercept`?|What is the `Slope`?|\n", "|--------|-----------|-------|----------|\n", "|Dummy coding|Choose one level as **reference**.|`mean` of reference group.|Difference between that level and `mean` of reference group.|\n", "|Deviation coding|Compare each level to **grand mean**.|`mean` of all groups.|Difference between that level and grand mean.|\n", "\n", "And many more!\n", "\n", "By default, `statsmodels` uses **dummy coding**." ] }, { "cell_type": "markdown", "id": "16be39fd", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Conclusion\n", "\n", "- The `statsmodels` package can be used to build **linear regression models**. \n", "- These models produce many useful *summary* features, including:\n", " - The fit $\\beta$ parameters (**intercept** and **slope**). \n", " - The **standard error** and **p-value** of those parameters.\n", "- Each $\\beta$ parameter has a particular interpretation with respect to $\\hat{Y}$." ] } ], "metadata": { "celltoolbar": "Slideshow", "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.9.12" } }, "nbformat": 4, "nbformat_minor": 5 }