{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Multivariate Thinking\n", "> A Summary of lecture \"Exploratory Data Analysis in Python\", via datacamp\n", "\n", "- toc: true \n", "- badges: true\n", "- comments: true\n", "- author: Chanseok Kang\n", "- categories: [Python, Datacamp]\n", "- image: images/brfss-logreg.png" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Limits of simple regression\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "from empiricaldist import Pmf, Cdf\n", "from scipy.stats import linregress\n", "import statsmodels.formula.api as smf" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "brfss_original = pd.read_hdf('./dataset/brfss.hdf5', 'brfss')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using StatsModels" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LinregressResult(slope=0.06988048092105248, intercept=1.5287786243362973, rvalue=0.11967005884864361, pvalue=1.3785039162157718e-238, stderr=0.002110976356332355)\n", "Intercept 1.528779\n", "INCOME2 0.069880\n", "dtype: float64\n" ] } ], "source": [ "# Run regression with linregress\n", "subset = brfss_original.dropna(subset=['INCOME2', '_VEGESU1'])\n", "xs = subset['INCOME2']\n", "ys = subset['_VEGESU1']\n", "res = linregress(xs, ys)\n", "print(res)\n", "\n", "# Run regression with StatsModels\n", "results = smf.ols('_VEGESU1 ~ INCOME2', data=brfss_original).fit()\n", "print(results.params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Multiple regression" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "gss = pd.read_hdf('./dataset/gss.hdf5', 'gss')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot income and education" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Income (1986 $)')" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEGCAYAAACkQqisAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df5xddX3n8dcbMjqjhhmQOB0TMFFiBVNBGDDVNY8ULA5sK9qqYLsmtWxntVBru27FtVsUXVe2rax0FTdKasAf/LC1IMXEGMVoV35MlB/BKImESsJ0gIZcIkx0Qj7943wnPTO5d+6dO3PunTvzfj4e9zHnfs/33Pu5Z3585vvjfI8iAjMzs3oc0ewAzMysdTmJmJlZ3ZxEzMysbk4iZmZWNycRMzOr27xmB9Boxx57bCxevLjZYZiZtYwtW7Y8HhELyu2bc0lk8eLFDAwMNDsMM7OWIemfK+1zd5aZmdXNScTMzOrmJGJmZnVzEjEzs7o5iZiZWd3m3OwsM7O5ZNtgifVbh9i9d5iFXR30LevmxJ7OaXt9t0TMzGapbYMl1mzeSWl4hJ7OdkrDI6zZvJNtg6Vpew8nETOzWWr91iE6O9ro7GjjCOnQ9vqtQ9P2Hk4iZmaz1O69w8xvHztqMb99Hrv3Dk/beziJmJnNUgu7Oti3/8CYsn37D7Cwq2Pa3sNJxMxslupb1k1peITS8AgHIw5t9y3rnrb3cBIxM5ulTuzppH/FEjo72hgs7aezo43+FUumdXaWp/iamc1iJ/Z0TmvSGM8tETMzq5uTiJmZ1c1JxMzM6uYkYmZmdXMSMTOzujmJmJlZ3ZxEzMysbk4iZmZWNycRMzOrm5OImZnVrdAkIqlL0pcl/UjSNkm/KukYSRslbU9fj051JelKSTsk3Svp1NzrrE71t0tanSs/TdJ96ZgrJanIz2NmZmMV3RL5BLA+Il4GnAxsAy4BNkXEUmBTeg5wDrA0PfqBqwAkHQNcCrwKOAO4dDTxpDr9ueP6Cv48ZmaWU1gSkXQUsAK4GiAifhERe4HzgHWp2jrgjWn7POCayNwOdEnqAV4PbIyIPRHxBLAR6Ev7joqI70VEANfkXsvMzBqgyJbIi4HHgL+V9ANJn5X0XKA7IgYB0tcXpPoLgYdzx+9KZROV7ypTfhhJ/ZIGJA089thjU/9kZmYGFJtE5gGnAldFxCuBp/j3rqtyyo1nRB3lhxdGrImI3ojoXbBgwcRRm5lZzYq8n8guYFdE3JGef5ksiQxJ6omIwdQl9Wiu/nG54xcBj6TylePKb0vli8rUNzObVbYNlli/dYjde4dZ2NVB37LuQu8RMhmFtUQi4l+AhyX9cio6C/ghcDMwOsNqNXBT2r4ZWJVmaS0HSqm7awNwtqSj04D62cCGtG+fpOVpVtaq3GuZmc0K2wZLrNm8k9LwCD2d7ZSGR1izeSfbBkvNDg0o/s6GfwR8QdKzgAeBd5AlrhskXQj8FHhLqnsrcC6wA3g61SUi9kj6MHBXqndZROxJ2+8CPgd0AF9LDzOzWWP91iE6O9ro7GgDOPR1/dahGdEaKTSJRMTdQG+ZXWeVqRvARRVeZy2wtkz5ALBsimGamc1Yu/cO09PZPqZsfvs8du8dblJEY/mKdTOzGWxhVwf79h8YU7Zv/wEWdnU0KaKxiu7OMjMz6h8c71vWzZrNO4GsBbJv/wFKwyOcf/qiKkc2hlsiZmYFm8rg+Ik9nfSvWEJnRxuDpf10drTRv2LJjBgPAbdEzMwKN9XB8RN7OmdM0hjPLREzs4Lt3jvM/Pax/7PPpMHxqXASMTMr2EwfHJ8KJxEzs4L1LeumNDxCaXiEgxGHtvuWdTc7tClzEjEzK9hMHxyfCg+sm5k1wEweHJ8Kt0TMzKxuTiJmZlY3JxEzM6ubk4iZmdXNScTMzOrmJGJmZnVzEjEzs7o5iZiZWd2cRMzMrG5OImZmVjcnETMzq5uTiJmZ1c1JxMzM6uYkYmZmdXMSMTOzujmJmJlZ3XxTKjOzGm0bLLF+6xC79w6zsKuDvmXds/JGU5PhloiZWQ22DZZYs3knpeERejrbKQ2PsGbzTrYNlpodWlM5iZiZ1WD91iE6O9ro7GjjCOnQ9vqtQ80OramcRMzMarB77zDz28eOAMxvn8fuvcNNimhmKDSJSHpI0n2S7pY0kMqOkbRR0vb09ehULklXStoh6V5Jp+ZeZ3Wqv13S6lz5aen1d6RjVeTnMbO5a2FXB/v2HxhTtm//ARZ2dTQpopmhES2RX4uIUyKiNz2/BNgUEUuBTek5wDnA0vToB66CLOkAlwKvAs4ALh1NPKlOf+64vuI/jpnNRX3LuikNj1AaHuFgxKHtvmXdzQ6tqZoxO+s8YGXaXgfcBrwvlV8TEQHcLqlLUk+quzEi9gBI2gj0SboNOCoivpfKrwHeCHytYZ/EzFrKVGZXndjTSf+KJWOOP//0RXN+dlbRSSSAr0sK4P9FxBqgOyIGASJiUNILUt2FwMO5Y3elsonKd5UpP4ykfrIWC8cff/xUP5OZtaDR2VWdHW1jZlf1r1gyqUQy15PGeEV3Z70mIk4l66q6SNKKCeqWG8+IOsoPL4xYExG9EdG7YMGCajGb2Szk2VXFKDSJRMQj6eujwFfIxjSGUjcV6eujqfou4Ljc4YuAR6qULypTbmZ2GM+uKkZhSUTScyXNH90Gzga2AjcDozOsVgM3pe2bgVVpltZyoJS6vTYAZ0s6Og2onw1sSPv2SVqeZmWtyr2WmdkYnl1VjJrGRCT1Aq8FXggMkyWDb4wOdlfQDXwlzbqdB3wxItZLugu4QdKFwE+Bt6T6twLnAjuAp4F3AETEHkkfBu5K9S7Lve+7gM8BHWQD6h5UN7Oy+pZ1s2bzTiBrgezbf4DS8Ajnn76oypE2EWWToSrslH4PeDewE9hC1vXUDrwUeA1ZMvkfEfHTwiOdJr29vTEwMNDsMMysCbz2VX0kbcldpjFGtZbIc8kGx8t2Gko6hez6jJZJImY2d3l21fSbMIlExCer7L97esMxM7NWUnVgXdIL0sA4kjokfUDSx0ZnWJmZ2dxVy+ys64Dnp+0PAScATwBfLCooMzNrDRMmkbTY4UuAlWn7fGAA+BfgRZJWSXpF8WGamdlMVG1g/TayKb3bgE5gCPgq2dXiF6f9c/uOLGZmc1i1gfV/lvQJ4BagDVgVET+VdDzweCtN7TUzs+lX9WLDiLhK0rXAwYh4OhX/K/C2QiMzM7MZr6Yr1iPiZ+OeP1VMOGZm1kp8e1wzM6ubk4iZmdVtUklE0lHpvuZHV69tZmazXbXrRD4v6di0/XrgfuBy4G5Jb5noWDMzm/2qDayfHBGPp+1LgddGxEMpsWwCbiw0OjMzm9GqdWcdIemotH2QtFpvSixF35/dzMxmuGqJ4EPAtyR9Evgn4EZJNwFnAuuLDs7MzGa2ales3yDpB8B/JrsR1TzgV4EvRcSGBsRnZmYzWC1XrG8H3teAWMzMrMVUm531YklrJX1E0vMkfUbSVkk3SlrcmBDNzGymqjaw/jngLuBnwO3Aj4FzyMZD1hYamZmZzXjVksj8iLgqIj4GHBURfxURD0fE1YAvODQzm+OqJZGDkl4q6XTgOZJ6ASSdABxZeHRmZjajVRtY/zOym1AdBN4IvF/SycBRwB8UHJuZmc1w1ab4bgJ+OVf03XS1+hMR8UyhkZmZ2YxXdYqvpOcBfcBxwAFgO/D1guMyM7MWUG2K71uBb5ElkYuBM4C3ky3A+IriwzMzs5msWkvkz4HlEfF06sb6QkS8PiWQTwOvLjxCMzObsarNzhIwnLafAl4AEBH3kg2um5nZHFYtidwKrJf038nGQW4EkHQMWYKpStKRkn4g6Zb0fImkOyRtl3S9pGel8men5zvS/sW513h/Kv9xuq/JaHlfKtsh6ZJJfG4zM5sGEyaRiHgf8AngF8BlEfHRtGsvcGqN7/HHwLbc88uBKyJiKfAEcGEqv5Bs1tcJwBWpHpJOAi4AXk42NvOplJiOBD5JdgX9ScDbUl0zM2uQqrfHjYhb05XqG3NlByPi59WOlbQI+I/AZ9NzkS0j/+VUZR3Z9ScA56XnpP1npfrnAddFxM8jYiewg2yA/wxgR0Q8GBG/AK5Ldc3MrEEmdY/1PEn31VDt/5BdsHgwPX8+sDciDqTnu4CFaXsh8DBA2l9K9Q+VjzumUrmZmTXIhLOzJP1WpV3AL1U59jeARyNii6SVuePGiyr7KpWXS4BRpgxJ/UA/wPHHHz9B1GZmNhnVpvheD3yB8n+c26sc+xrgDZLOTXWPImuZdEmal1obi4BHUv1dZBc07pI0D+gE9uTKR+WPqVQ+RkSsAdYA9Pb2lk00ZmY2edWSyL3AX0XE1vE7JL1uogMj4v3A+1PdlcB7I+J3Jd0IvJlsDGM1cFM65Ob0/Htp/zcjIiTdDHxR0seBFwJLgTvJWihLJS0BdpMNvv9O1U9sZmbTploSeQ/wZIV9b6rzPd8HXCfpI8APgKtT+dXAtZJ2kLVALgCIiPsl3QD8kGzZlYtG1+2SdDGwgWxF4bURcX+dMZmZWR0UMbd6d3p7e2NgYKDZYZiZtQxJWyKit9y+agPrfw58KiL2VNh/JvCciLhl6mGamVW3bbDE+q1D7N47zMKuDvqWdXNiT2ezw5qzqnVn3Qd8VdJ+4PvAY2SD5EuBU4BvAB+tfLiZ2fTZNlhizeaddHa00dPZTml4hDWbd9K/YokTSZNUu5/ITcBNkpaSzbbqIRsj+TzQHxHDEx1vZjad1m8dorOjjc6ONoBDX9dvHXISaZKq9xMBiIjtZPcRMTObsnq7pHbvHaanc+zVBfPb57F7r/+fbZa6r1g3M6vHaJdUaXhkTJfUtsFS1WMXdnWwb/+BMWX79h9gYVdHUeFaFU4iZtZQ+S6pI6RD2+u3DlU9tm9ZN6XhEUrDIxyMOLTdt6y7AZFbOU4iZtZQu/cOM799bE96rV1SJ/Z00r9iCZ0dbQyW9tPZ0eZB9SaraUxE0kuBq4DuiFiW7mz4hoj4SKHRmdmss7Crg9LwyKFBcZhcl9SJPZ1OGjNIrS2Rz5AtYTICh+5seEFRQZnZ7OUuqdml1iTynIi4c1zZgbI1zcwm4C6p2aWm7izgcUkvIa3mK+nNwGBhUZnZrOYuqdmj1iRyEdlS6i+TtBvYCfynwqIyM7OWUOvFhg8Cr5P0XOCIiNhXbFhmZtYKap2d1QWsAhYD87Jbn0NEvLuwyMzMbMartTvrVuB2sgUZD1apa2Zmc0StSaQ9Iv600EjMzKzl1JpErpX0B8AtwM9HCyvdZ8TMZj/f18Og9utEfgH8Jdn9z7ekh28PaDZHTWURRZtdam2J/ClwQkQ8XmQwZtYafF8PG1VrS+R+4OkiAzGz1jGVRRRtdqm1JfIMcLekbzF2TMRTfM3moKkuomizR61J5B/Sw8yMvmXdrNm8E8haIPv2H6A0PML5py9qcmTWaLVesb5O0rOAl6aiH0fESHFhmdlMNrqIYn521vmnL/J4yBxU6xXrK4F1wEOAgOMkrY6IzcWFZmYzmRdRNKi9O+uvgbMj4sdw6CZVXwJOKyowMzOb+WqdndU2mkAAIuIBoG2C+mZmNgfU2hIZkHQ1cG16/rtkFxyamdkcVmsSeRfZPUXeTTYmshn4VFFBmZlZa6g1icwDPhERHweQdCTw7MKiMrPCee0rmw61jolsAvJXEXUA35joAEntku6UdI+k+yV9KJUvkXSHpO2Srk9Th5H07PR8R9q/OPda70/lP5b0+lx5XyrbIemSGj+L2Zznta9sutSaRNoj4mejT9L2c6oc83PgzIg4GTgF6JO0HLgcuCIilgJPABem+hcCT0TECcAVqR6STgIuAF4O9AGfknRkag19EjgHOAl4W6prZlXk1746Qjq0vX7rULNDsxZTaxJ5StKpo08knQZMuEhOZEYTT1t6BHAm8OVUvg54Y9o+Lz0n7T9L2S0UzwOui4ifR8ROYAdwRnrsiIgHI+IXwHWprplV4bWvbLrUOibyHuBGSY+k5z3A+dUOSq2FLcAJZK2GnwB7I+JAqrILWJi2FwIPA0TEAUkl4Pmp/Pbcy+aPeXhc+asqxNEP9AMcf/zx1cI2m/W89pVNl5paIhFxF/AysllafwicGBFVp/hGxDMRcQqwiKzlcGK5aumrKuybbHm5ONZERG9E9C5YsKBa2GazXt+ybkrDI5SGRzgYcWi7b1l3s0OzFlNrdxbA6cArgFeSjT+sqvXAiNgL3AYsB7okjbaAFgGjrZtdwHEAaX8nsCdfPu6YSuVmVsXo2ledHW0MlvbT2dFG/4olnp1lk1br2lnXAi8B7iZbFh6y//qvmeCYBcBIROyV1AG8jmyw/FvAm8nGMFYDN6VDbk7Pv5f2fzMiQtLNwBclfRx4IbAUuJOsJbJU0hJgN9ng++/U+LnN5jyvfWXTodYxkV7gpIgo211UQQ+wLo2LHAHcEBG3SPohcJ2kjwA/AK5O9a8mu5f7DrIWyAUAEXG/pBuAHwIHgIsi4hkASRcDG4AjgbURcf8k4jMzsylSLXlB0o3AuyNisPiQitXb2xsDA749vJlZrSRtiYjecvtqbYkcC/xQ0p2MvbPhG6YhPjMza1G1JpEPFhmEmZm1plrvbPjtogMxM7PWM2ESkbSP8tdeiOyi9KMKicrMzFrChEkkIuY3KhAzM2s9k7nY0MzMbAwnETMzq5uTiJmZ1c1JxMzM6lbrdSLWBL59qZnNdG6JzFC+famZtQInkRnKty81s1bg7qwZavfeYXo628eU+falNp67PK3Z3BKZoRZ2dbBv/4ExZb59qeW5y9NmAieRGcq3L7Vq3OVpM4GTyAzl25daNbv3DjO/fWyPtLs8rdE8JjKD+falNpGFXR2Uhkfo7Gg7VOYuT2s0t0TMWpS7PG0mcBIxa1Hu8rSZwN1ZZi3MXZ7WbE4iZk3maz2slbk7y6yJfK2HtTonEbMm8rUe1uqcRMyayNd6WKtzEjFrIi9vY63OScSsiXyth7U6JxGzJvK1HtbqPMXXrMl8rYe1ssJaIpKOk/QtSdsk3S/pj1P5MZI2Stqevh6dyiXpSkk7JN0r6dTca61O9bdLWp0rP03SfemYKyWpqM9jZmaHK7I76wDwXyPiRGA5cJGkk4BLgE0RsRTYlJ4DnAMsTY9+4CrIkg5wKfAq4Azg0tHEk+r0547rK/DzmJnZOIUlkYgYjIjvp+19wDZgIXAesC5VWwe8MW2fB1wTmduBLkk9wOuBjRGxJyKeADYCfWnfURHxvYgI4Jrca5mZWQM0ZExE0mLglcAdQHdEDEKWaCS9IFVbCDycO2xXKpuofFeZcrOG89IlNlcVPjtL0vOAvwPeExFPTlS1TFnUUV4uhn5JA5IGHnvssWohm02Kly6xuazQJCKpjSyBfCEi/j4VD6WuKNLXR1P5LuC43OGLgEeqlC8qU36YiFgTEb0R0btgwYKpfSizcbx0ic1lRc7OEnA1sC0iPp7bdTMwOsNqNXBTrnxVmqW1HCilbq8NwNmSjk4D6mcDG9K+fZKWp/dalXsts0nbNljiio0P8N4b7+GKjQ/U3JLw0iU2lxXZEnkN8HbgTEl3p8e5wMeAX5e0Hfj19BzgVuBBYAfwGeAPASJiD/Bh4K70uCyVAbwL+Gw65ifA1wr8PDaLTaVLykuX2FxW2MB6RHyX8uMWAGeVqR/ARRVeay2wtkz5ALBsCmGaAWO7pIBDX9dvHao6QN63rJs1m3cCWQtk3/4DlIZHOP/0RRMeZzYbeNkTM6bWJeWlS2wu87InZmRdUqXhkUMtEJhcl5SXLrG5yknEZo2pXKvhLimz+rg7y2aFqV6r4S4ps/q4JWLTrhlXb09lYHyUu6TMJs8tEZtWzbp629dqmDWHk4hNq2Zdve1rNcyaw0nEplWzWgS+zaxZcziJ2LRqVovAA+NmzeGB9Rp4me/aNXOqrAfGzRrPLZEqvMz35LhFYDa3uCVSxXRMHZ1r3CIwmzucRKrYvXeYns72MWWeOlocdx2atRZ3Z1XhqaON465Ds9bjlkgVc3VNpVa96tzMGsstkSrm4kCxrzo3s1q5JVKDuTZQ3KwWwVSXYzezxnMSKVgrDhQ3azLBXO06NGtl7s4qUKsOFPuqczOrlVsiBWrVgWJfdW5mtXJLpECtOlDsFoGZ1cotkQI1c6B4qmMxbhGYWS3cEilQs5Ynb9WxGDNrPU4iBWpWt1CzbgxlZnOPu7MK1oxuIa/3ZWaN4pbILOT1vsysUZxEZiHfKtbMGsVJZBbyFF0zaxSPicxSnqJrZo1QWEtE0lpJj0ramis7RtJGSdvT16NTuSRdKWmHpHslnZo7ZnWqv13S6lz5aZLuS8dcKUlFfRYzMyuvyO6szwF948ouATZFxFJgU3oOcA6wND36gasgSzrApcCrgDOAS0cTT6rTnztu/HuZmVnBCksiEbEZ2DOu+DxgXdpeB7wxV35NZG4HuiT1AK8HNkbEnoh4AtgI9KV9R0XE9yIigGtyr2VmZg3S6IH17ogYBEhfX5DKFwIP5+rtSmUTle8qU16WpH5JA5IGHnvssSl/CDMzy8yU2VnlxjOijvKyImJNRPRGRO+CBQvqDNHMzMZr9OysIUk9ETGYuqQeTeW7gONy9RYBj6TylePKb0vli8rUr2rLli2PS/rnuqKHY4HH6zy2SI5rchzX5DiuyZmNcb2o0o5GJ5GbgdXAx9LXm3LlF0u6jmwQvZQSzQbgo7nB9LOB90fEHkn7JC0H7gBWAX9TSwARUXdTRNJARPTWe3xRHNfkOK7JcVyTM9fiKiyJSPoSWSviWEm7yGZZfQy4QdKFwE+Bt6TqtwLnAjuAp4F3AKRk8WHgrlTvsogYHax/F9kMsA7ga+lhZmYNVFgSiYi3Vdh1Vpm6AVxU4XXWAmvLlA8Ay6YSo5mZTc1MGVhvFWuaHUAFjmtyHNfkOK7JmVNxKWsEmJmZTZ5bImZmVjcnETMzq5uTSBmS+iT9OC3ueEmZ/c+WdH3af4ekxQ2I6ThJ35K0TdL9kv64TJ2VkkqS7k6Pvyg6rvS+D6XFMO+WNFBmf8UFNguM6Zdz5+FuSU9Kes+4Og05X5NZjLTMsWUXIC0wrr+U9KP0ffqKpK4Kx074PS8grg9K2p37Xp1b4dgJf3cLiOv6XEwPSbq7wrFFnq+yfxsa9jMWEX7kHsCRwE+AFwPPAu4BThpX5w+BT6ftC4DrGxBXD3Bq2p4PPFAmrpXALU04Zw8Bx06w/1yyKdgClgN3NOF7+i/Ai5pxvoAVwKnA1lzZ/wYuSduXAJeXOe4Y4MH09ei0fXTBcZ0NzEvbl5eLq5bveQFxfRB4bw3f5wl/d6c7rnH7/xr4iyacr7J/Gxr1M+aWyOHOAHZExIMR8QvgOrIFIvPyC0l+GThLKnYp+ogYjIjvp+19wDYmWC9shqm0wGajnAX8JCLqXalgSmJyi5HmlV2AtMi4IuLrETF6b+XbGbsyRENUOF+1qOV3t5C40u//W4EvTdf71WqCvw0N+RlzEjlcpUUfy9ZJv3Al4PkNiQ5I3WevJLtaf7xflXSPpK9JenmDQgrg65K2SOovs7+Wc1qkC6j8y92M8wWVFyPNa/Z5+30qX8Rb7XtehItTN9vaCl0zzTxfrwWGImJ7hf0NOV/j/jY05GfMSeRwtSzuOKkFIKeTpOcBfwe8JyKeHLf7+2RdNieTLQPzD42ICXhNRJxKdl+YiyStGLe/mefrWcAbgBvL7G7W+apVM8/bB4ADwBcqVKn2PZ9uVwEvAU4BBsm6jsZr2vkC3sbErZDCz1eVvw0VDytTNqlz5iRyuEqLQZatI2ke0El9ze9JkdRG9kPyhYj4+/H7I+LJiPhZ2r4VaJN0bNFxRcQj6eujwFfIuhXyajmnRTkH+H5EDI3f0azzlQyNdulp7GKkeU05b2lw9TeA343UcT5eDd/zaRURQxHxTEQcBD5T4f2adb7mAb8FXF+pTtHnq8Lfhob8jDmJHO4uYKmkJem/2AvIFojMG11IEuDNwDcr/bJNl9TnejWwLSI+XqHOL42OzUg6g+z7+68Fx/VcSfNHt8kGZreOq3YzsEqZ5aQFNouMK6fif4jNOF85+Z+h/GKkeRuAsyUdnbpvzk5lhZHUB7wPeENEPF2hTi3f8+mOKz+G9qYK71fL724RXgf8KCJ2ldtZ9Pma4G9DY37Gipgt0OoPstlED5DN9PhAKruM7BcLoJ2se2QHcCfw4gbE9B/Impn3Anenx7nAO4F3pjoXA/eTzUq5HXh1A+J6cXq/e9J7j56vfFwCPpnO531Ab4O+j88hSwqdubKGny+yJDYIjJD953ch2RjaJmB7+npMqtsLfDZ37O+nn7MdwDsaENcOsj7y0Z+x0VmILwRuneh7XnBc16afnXvJ/jj2jI8rPT/sd7fIuFL550Z/pnJ1G3m+Kv1taMjPmJc9MTOzurk7y8zM6uYkYmZmdXMSMTOzujmJmJlZ3ZxEzMysbk4iNqtJekZjV/MttyrzSkm3TPP7rpT06tzzd0paNU2v3TPd8U7iva+TtLQZ720zU2H3WDebIYYj4pQmvO9K4GfA/weIiE9P42v/KdlV24WQdGREPFNh91XAnwF/UNT7W2txS8TmpHTfiR9J+i7ZkhWj5R+U9N7c861pUTskrUoLAN4j6dpU9pvK7inzA0nfkNSd6r8T+JPU+nlt/nUlnSLpdv37PTuOTuW3Sbpc0p2SHpD02grh/zawPh3zHUmHkqSkf5L0inSV9FpJd6XYzkv7F6djvp8er07lK5Xdk+KLwH3p+H9Mn3WrpPPTW3wHeF1a6sPMScRmvY5x3VnnS2on+0/+N8lWX/2lai+ibIXfDwBnRrZg4+hNwb4LLI+IV5ItPf5nEfEQ8Gngiog4JSK+M+7lrgHeFxGvILsK+9LcvnkRcQbwnnHlo3EsAZ6IiJ+nos8Cv5f2vRR4dkTcm2L9ZkScDvwa8JdpyY1HgV+PbDHA84Ercy9/BtnV1CeRLQf+SEScHBHLSEkrsrWrdgAnVztnNjc4idhsN5z+kI8+rgdeBuyMiO2RLdnw+Rpe50zgyxHxOEBEjC64ue5t7QIAAAIoSURBVAjYIOk+4L8BEy4nL6kT6IqIb6eidWQ3Oxo1unjeFmBxmZfoAR7LPb8R+I20AN/vky3BAdkaSJcou9PebWRL9RwPtAGfSfHeSHbzolF3RsTOtH0fWYvjckmvjYhSrt6jZMt6mDmJ2JxVab2fA4z9vWhPX1XhmL8B/m9E/ArwX3L16zXawniG8mOWw/n3iGyRxI1kNyB6K/DFXLy/nUuex0fENuBPgCGylkQv2R0ARz2Ve90HgNPIksn/0thbB7enOMycRGxO+hGwRNJL0vO35fY9RHYLVJTdC35JKt8EvFXS89O+Y1J5J7A7befvT72P7FalY6T/6J/IjXe8Hfj2+HoTeIDDWyifJeuWuivXQtoA/FFuleJX5uIdTN1Sbye7pexhJL0QeDoiPg/8FemcJC8lW0jQzEnEZr3xYyIfi4j9QD/wj2lgPX/b3L8DjkndQO8i+6NNRNwP/E/g25LuAUaX3P4gcKOk7wCP517nq8CbRgfWx8W0mmyM4l6ymyxdVuuHiYingJ9IOiFXtgV4EvjbXNUPk3Vd3Stpa3oO8ClgtaTbyZLBU5T3K8Cd6Tx8APgIgKRusi7CRi3lbzOcV/E1azGS3gScFhF/np6/kGzc42WphVHke/8J8GREXF3k+1jrcEvErMVExFfIut1IFzDeQTarqtAEkuwlmwxgBrglYmZmU+CWiJmZ1c1JxMzM6uYkYmZmdXMSMTOzujmJmJlZ3f4NhVxTd7NZtRUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Group by educ\n", "grouped = gss.groupby('educ')\n", "\n", "# Compute mean income in each group\n", "mean_income_by_educ = grouped['realinc'].mean()\n", "\n", "# Plot mean income as a scatter plot\n", "plt.plot(mean_income_by_educ, 'o', alpha=0.5)\n", "\n", "# Label the axes\n", "plt.xlabel('Education (years)')\n", "plt.ylabel('Income (1986 $)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Non-linear model of education" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "gss['age2'] = gss['age'] ** 2" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Intercept -23241.884034\n", "educ -528.309369\n", "educ2 159.966740\n", "age 1696.717149\n", "age2 -17.196984\n", "dtype: float64\n" ] } ], "source": [ "# Add a new column with educ squared\n", "gss['educ2'] = gss['educ'] ** 2\n", "\n", "# Run a regression model with educ, educ2, age and age2\n", "results = smf.ols('realinc ~ educ + educ2 + age + age2', data=gss).fit()\n", "\n", "# Print the estimated parameters\n", "print(results.params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualizing regression results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Making predictions" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 12182.344976\n", "1 11993.358518\n", "2 11857.672098\n", "3 11775.285717\n", "4 11746.199374\n", "dtype: float64\n" ] } ], "source": [ "# Run a regression model with educ, educ2, age, and age2\n", "results = smf.ols('realinc ~ educ + educ2 + age + age2', data=gss).fit()\n", "\n", "# Make the DataFrame\n", "df = pd.DataFrame()\n", "df['educ'] = np.linspace(0, 20)\n", "df['age'] = 30\n", "df['educ2'] = df['educ'] ** 2\n", "df['age2'] = df['age'] ** 2\n", "\n", "# Generate and plot the predictions\n", "pred = results.predict(df)\n", "print(pred.head())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualizing predictions" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEGCAYAAACkQqisAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhU1fnA8e8bEkiAEFZDAJEoiwFUhIj7hqIRrUvdUCtQrbSKtmr7U6y1tmpbtVZb96JSwQUQXMAtFEFALIqgKIGAhD0QA4KEAAlkeX9/3BMcwiQZQu5MJnk/zzPP3Dn33DtvJsubc86954iqYowxxtRGTKQDMMYYE70siRhjjKk1SyLGGGNqzZKIMcaYWrMkYowxptZiIx1AuLVv3167desW6TCMMSZqLFq06HtV7RBsX6NLIt26dWPhwoWRDsMYY6KGiKyrap91ZxljjKk1SyLGGGNqzZKIMcaYWmt0YyLBlJSUkJubS3FxcaRDqVfi4+Pp0qULcXFxkQ7FGFNPWRIBcnNzSUxMpFu3bohIpMOpF1SVrVu3kpubS2pqaqTDMcbUU9adBRQXF9OuXTtLIAFEhHbt2lnrzBhTLUsijiWQA9lnYoypiSURY4xp6NbOg8+eAx+W/rAkUo+8/fbbiAjLly+vs3M+//zzHHPMMfTr14/TTjuNZcuW7dv3t7/9je7du9OrVy+mT59eZ+9pjKk/VqxZR+HrP2fLrKd5avo3ZOcV1On5LYnUIxMmTOC0005j4sSJdXbOa6+9liVLlrB48WLuuusu7rzzTgCWLVvGxIkTWbp0KZmZmdxyyy2UlZXV2fsaYyIve9N2it66jeZ7tzK910Ns3RvLmLlr6jSRWBKpJ3bu3Mmnn37KSy+9tF8SKS8v55ZbbqFPnz5cdNFFDBkyhClTpgCwaNEizjzzTAYMGMD5559PXl7eAedt1arVvu1du3btG+eYOnUqQ4cOpVmzZqSmptK9e3cWLFjg81dpjAmnDbNepF/hHOYf8Su2tOpDUkIcSQlxZGbl19l72CW+lX04Gr5bUrfn7HgMXPBwtVXeeecdMjIy6NmzJ23btuXLL7+kf//+vPXWW6xdu5YlS5awefNm0tLSuOGGGygpKeG2225j6tSpdOjQgUmTJnHvvfcyduzYA879zDPP8Pjjj7N3715mzZoFwMaNGznppJP21enSpQsbN26s26/bGBM5W1dx5urHWN9qAAs7X7+vODE+lo3bi+rsbawlUk9MmDCBoUOHAjB06FAmTJgAwLx587jyyiuJiYmhY8eOnH322QCsWLGCrKwsBg8eTL9+/XjooYfIzc0Neu5Ro0axatUqHnnkER566CHAuw+kMrsay5gGoqwE3rwRjYljStf7UGmyb1dhcSmdWyfU2VtZS6SyGloMfti6dSuzZs0iKysLEaGsrAwR4dFHHw36xx68JNCnTx/mz58f8vsMHTqUm2++GfBaHhs2bNi3Lzc3l06dOh3aF2KMqR8+/its+orvBz/Phtw2JBWVkBgfS2FxKQVFJVx9Qpc6eytridQDU6ZMYdiwYaxbt461a9eyYcMGUlNTmTdvHqeddhpvvvkm5eXl5OfnM3v2bAB69erFli1b9iWRkpISli5desC5V65cuW/7/fffp0ePHgBcfPHFTJw4kT179rBmzRpWrlzJwIED/f9ijTH+WvMJzHsCjr+eLqdew8gzUklKiCOvoJikhDhGnpFKWkpSnb2dtUTqgQkTJjB69Oj9yi6//HJef/11nnnmGWbOnEnfvn3p2bMnJ554IklJSTRt2pQpU6bw61//moKCAkpLS7n99tvp06fPfud5+umn+eijj4iLi6NNmzaMGzcOgD59+nDVVVfRu3dvYmNjeeaZZ2jSpAnGmCi2exu8/UtoeyRkeL0qaSlJdZo0KpOquksaqvT0dK28KFV2djZpaWkRiqhmO3fupGXLlmzdupWBAwfy6aef0rFjx7C8d33/bIwxjipMHg7L34cb/wudB9TZqUVkkaqmB9tnLZEocNFFF7F9+3b27t3LfffdF7YEYoyJIotehmVT4dw/1WkCqYklkShQMQ5ijDFBbc6GzNFw5Nlwym/C+tY2sO40tm69UNhnYkwUKCmCyT+HZolw2b8hJrx/1i2J4C2+tHXrVvujGaBiPZH4+PhIh2KMqc7038OWbLjseUhMDvvbW3cW3j0Tubm5bNmyJdKh1CsVKxsaY+qpZVNh4Vg45Tbofm5EQrAkAsTFxdnqfcaY6LJ9PUy7DTr1h0F/jFgY1p1ljDHRpqwU3rwJysvhipcgtmnEQrGWiDHGRJs5D8OGz+Dyl7wbCyPIWiLGGBNNVs+GuY9Bv5/BMVdEOhp/k4iItBaRKSKyXESyReRkEWkrIjNEZKV7buPqiog8KSI5IvKNiPQPOM9wV3+liAwPKB8gIkvcMU+KTUNrjGnICvO9bqz2PWHIo5GOBvC/JfIvIFNVjwaOA7KB0cBMVe0BzHSvAS4AerjHSOA5ABFpC9wPnAgMBO6vSDyuzsiA4zJ8/nqMMSYyysvgrV/AnkK4ahw0bRHpiAAfk4iItALOAF4CUNW9qroduAQY56qNAy5125cA49XzGdBaRFKA84EZqrpNVX8AZgAZbl8rVZ2v3g0e4wPOZYwxDcvcx2DNXBjydzis/sxn52dL5EhgC/AfEflKRF4UkRZAsqrmAbjnw1z9zsCGgONzXVl15blByg8gIiNFZKGILLR7QYwxUWfNXJj9Nzh2KBz/s0hHsx8/k0gs0B94TlWPB3bxY9dVMMHGM7QW5QcWqo5R1XRVTe/QoUP1URtjTH2yczOlk29ga8IR/H7vCJ74aCXZeQWRjmofP5NILpCrqp+711Pwkkq+64rCPW8OqH94wPFdgE01lHcJUm6MMQ1DeRk7J/yc8qICJnV7kHZt21JQVMKYuWvqTSLxLYmo6nfABhHp5YrOAZYB04CKK6yGA1Pd9jRgmLtK6ySgwHV3TQfOE5E2bkD9PGC621coIie5q7KGBZzLGGOi3yeP03LjPDIPv5M9bY8mRoSkhDiSEuLIzMqPdHSA/zcb3ga8JiJNgdXAz/ES1xsiciOwHrjS1f0AGALkALtdXVR1m4g8CHzh6j2gqtvc9s3Ay0AC8KF7GGNM9FvzCcz+K18mncvKzpft9x9/YnwsG7cXRSy0QL4mEVVdDARbDeucIHUVGFXFecYCY4OULwT6HmKYxhhTvxR+B1NugLZH8b/uf6BwTxlJCT+mkcLiUjq3TohggD+yaU+MMSYMsvMKyMzKZ+P2Ijq3TiCjb3Lwtc/LSr0EsncnDJ/GOWWdGDN3DeC1QAqLSykoKuHqE+rHDNs27YkxxvgsO6+AMXPXUFBUQkpSfPWD47MehHWfwk/+BYelkZaSxMgzUklKiCOvoJikhDhGnpEaPAFFgLVEjDHGZ5lZ+fsGxIF9z5lZ+fsng+UfwKf/hPQb4Nir9hWnpSTVm6RRmbVEjDHGZxu3F5EYv///7AcMjm9bA+/8ClL6wfl/C3OEtWdJxBhjfNa5dQKFxaX7le03OF5SDJPdnQ9XjYO46FmW2pKIMcb4LKNvMgVFJRQUlVCuum87o69bEz1zNOR9DZf9G9p0i2isB8uSiDHG+KzawfHFE2DRf+C0O6DXBZEO9aDZwLoxxoRB0MHxvG/gvduh2+lw9h8iE9ghspaIMcZEwu5tMOln0LwdXPEfaBKd/9NHZ9TGGBPNysvhrZGwYxPckAkto3d2cUsixhgTbnMehpwZcOHj0CXYzFDRw7qzjDEmnFZkwpxHoN913k2FUc6SiDHGhMvWVV43VspxcOE/QIKtrRddLIkYY0w47N0Fk66HmBi46hWIqx+z8B4qGxMxxhi/qcK0X8PmZfCzN6HNEZGOqM5YS8QYY/z2v6cgawoM+gN0P2A5pahmScQYY/yUMxM+uh96XwKn/zbS0dQ5SyLGGOOXbau9BaY6pMElzzaIgfTKLIkYY4wf9uyEidd520Nfg2YtIxuPT2xg3Rhj6poqvHMzbFnuDaS3TY10RL6xJGKMMXXtk39A9jQ47yE4alCko/GVdWcZY0xd+nY6zHoIjrkSTr410tH4zpKIMcbUlS0r4M1fQMdj4CdPNsiB9MqsO8sYY0KUnVdAZlY+G7cX0bl1Ahl9k39cI2T3NpgwFGKbwdDXoWnzyAYbJtYSMcaYEGTnFTBm7hoKikpISYqnoKiEMXPXkJ1XAGUlMHkEFOTC1a9B68MjHW7YWEvEGGNCkJmVT1JCHEkJcQD7njOz8kkreQHWzIFLn4OuJ0YyzLCzJGKMMSHYuL2IlKT4/coS42M5fPVEyBvjDaL3uzZC0UWOr91ZIrJWRJaIyGIRWejK2orIDBFZ6Z7buHIRkSdFJEdEvhGR/gHnGe7qrxSR4QHlA9z5c9yxDX8UyxgTEZ1bJ1BYXLpfWbstC7jsu39C98Ew+IEIRRZZ4RgTOVtV+6lqxfJdo4GZqtoDmOleA1wA9HCPkcBz4CUd4H7gRGAgcH9F4nF1RgYcl+H/l2OMaYwy+iZTUFRCQVEJ5arEbF/LFavvpTQpFa54CWKaRDrEiIhEd9YlwFluexwwG7jblY9XVQU+E5HWIpLi6s5Q1W0AIjIDyBCR2UArVZ3vyscDlwIfhu0rMcZElWqvrqpBWkoSI89IJTMrn++3fs8d60eTEAtNrn8D4kM7R0Pkd0tEgf+KyCIRGenKklU1D8A9H+bKOwMbAo7NdWXVlecGKTfGmANUe3VViNJSkrhj0JH8pfRx2hevo8nV46HdUT5GXf/53RI5VVU3ichhwAwRWV5N3WDjGVqL8gNP7CWwkQBdu3atPmJjTINU7dVVIbZGAJj+e8iZAT/5Fxx5Vt0HGmV8bYmo6ib3vBl4G29MI991U+GeN7vquUDgxdVdgE01lHcJUh4sjjGqmq6q6R06dDjUL8sYE4U2bi8iMX7//5sT42PZuL0o9JMseAEW/Nu7EmvAiLoNMEqFlEREJF1E7hCRv4vIAyJylRvwru6YFiKSWLENnAdkAdOAiiushgNT3fY0YJi7SuskoMB1d00HzhORNm5A/TxguttXKCInuauyhgWcyxhj9hPs6qrC4lI6tw5xrfOVH8GHd0GvIY32Sqxgqk0iIjJCRL4E7gESgBV4LYfT8LqnxolIVf1DycA8EfkaWAC8r6qZwMPAYBFZCQx2rwE+AFYDOcALwC0AbkD9QeAL93igYpAduBl40R2zChtUN8ZUofLVVRXbGX2Taz44f5l3R3pyH/jpC432SqxgxLsYqoqdIqOAsaoatL0nIv2Adqo606f46lx6erouXLgw0mEYYyKgVldn7dwML5wDZXvhplmQ1Piu3xGRRQG3aeyn2oF1VX2mhv2LDyUwY4wJp7SUpIMbRC8pgonXwq4tcMOHjTKB1KTGq7PclVW7VHWXiCQAdwKJwL8qLtU1xpgGp7wc3hoJuV/AVa9Ap+MjHVG9FMrA+kSgndv+M9Ad+AF43a+gjDEm4mbc9+PqhL0vjnQ09VZNA+vDgaOAs9z21cBC4DvgCBEZJiLH+h+mMcaE0YIXYP7TcMJNjWJ1wkNRU3fWbKAIyAaSgHzgXbwb/W51+0O/3dMYY+q7FR96l/L2vAAueKRRrE54KGoaWF8nIv8C3gPigGGqut5d1vu9qq4PR5DGGBMWG7+EKTdAynGNelLFg1HjwLqqPicirwDlqrrbFW8FrvE1MmOMCacf1sHrV0Pz9nDNJGjaItIRRYWQ5s5S1Z2VXu/yJxxjjImAoh/gtSuhbA+MeA8SQ7gB0QC2sqExprErKYaJ18G21XD929ChV6QjiiqWRIwxjVd5Gbx1E6z7FC5/CVJPj3REUeegZvEVkVZuSdo2Ndc2xph6TBUyR3v3gpz/VzjmikhHFJVquk/kVRFp77bPB5YCjwCLReTKMMRnjDH+mPcELBjj3Qdy8qhIRxO1aurOOk5Vv3fb9wOnq+pal1hmApN9jc4YY/yw+HWY+Wc45koY/GCko4lqNXVnxYhIK7ddDqwHcInFxlOMMdFn5Ucw9VZvVcJLnoUYv1cJb9hqSgR/Bj4WkWeAT4HJIjIVGARk+h2cMcbUqY1fwhvDILm3N6libNNIRxT1arpj/Q0R+Qr4BdDT1T8ZmKCq08MQnzHG1I3vV8JrV0CLdnDdmxDfquZjTI1CuWN9JXB3GGIxxhh/FOTC+EtBYuD6d+xmwjpU09VZR4rIWBF5SERaisgLIpIlIpNFpFt4QjTGmEOwayu8chns2QE/exPaHRXpiBqUmkaUXsZb13wn8BneGusX4I2HjPU1MmOMOVR7Cr0urB/WwTUTvIkVTZ2qKYkkqupzqvow0EpVH1PVDar6EmA3HBpj6q/SPTDpZ5D3NVz5MnQ7LdIRNUg1JZFyEekpIicAzUUkHUBEugM2R7Ixpn6qmM5k9Wy45Gk4ekikI2qwahpYvwtvEapy4FLgHhE5DmgF3ORzbMYYc/BU4f07YdlUOO8v0O/aSEfUoNV0ie9MIHBKy3nubvUfVLXM18iMMeZgqcJ//wCLXobT7oBTbGlbv9V4ia+ItAQygMOBUmAl8F+f4zLGmIM359Ef10Y/5/5IR9Mo1HSJ71XAx3hJ5FZgIHA93gSMx/ofnjHGhGj+MzD7r3DctXDBo7Y2epjU1BL5A3CSqu523Vivqer5LoE8D5zie4TGGFOTRS/D9N9D2sVw8VM2H1YY1fRJC1DktncBhwGo6jd4g+vGGBNZ30yGd2+H7oO9haWa2Nyw4VTTp/0BkCkic/BuMpwMICJt8RKMMcZEzvL34e1fwhGnwtU2oWIkVNsSUdW7gX8Be4EHVPWvbtd2oH8obyAiTUTkKxF5z71OFZHPRWSliEwSkaauvJl7neP2dws4xz2ufIVbHKuiPMOV5YjI6IP4uo0x0S7nI5g8Ajr1g2snQlxCpCNqlGrsOFTVD9yd6jMCyspVdU+I7/EbIDvg9SPAE6raA/gBuNGV34h36XB34AlXDxHpDQwF+uAN8D/rElMT4Bm8FlJv4BpX1xjT0K2eDROvgw694Lop0Cwx0hE1WrUefRKRJSHU6QJcCLzoXgveWiRTXJVxeDcxAlziXuP2n+PqXwJMVNU9qroGyMG7SmwgkKOqq1V1LzDR1TXGNGRr58HrQ6HtUXD9VGjeNtIRNWrVjomIyE+r2gV0DOH8/8S7673i34R2wHZVLXWvc4HObrszsAFAVUtFpMDV74w3+SNBjtlQqfzEKr6OkcBIgK5du4YQtjGmXlr/Gbx2FbTuCsOmemuDmIiqaWB9EvAaoEH2xVd3oIhcBGxW1UUiclZFcZCqWsO+qsqDtaKCxYmqjgHGAKSnpwetY4yp5zZ8Aa9eAa1SYPg0aNkh0hEZak4i3wCPqWpW5R0icm4Nx54KXCwiQ/ASTiu8lklrEYl1rZEuwCZXPxfvrvhcEYkFkoBtAeUVAo+pqtwY05Bs/BJevRxatIfh70JiKB0hJhxqGhO5HdhRxb7LqjtQVe9R1S6q2g1vYHyWql6Hdwf8Fa7acGCq257mXuP2z1JVdeVD3dVbqUAPYAHeOic93NVeTd17TKvh6zHGRJu8r71FpRKSvATSqlOkIzIBapqA8ZNq9i2s5XveDUwUkYeAr4CXXPlLwCsikoPXAhnq3mepiLwBLMObu2tUxeSPInIrMB1vWvqxqrq0ljEZY+qjTV95y9o2S4Th70Hrw2s+xoSVeP/sV7FT5A/As6q6rYr9g4DmqvqeT/HVufT0dF24sLb5zxgTNhsXeS2QZkkw4j1oc0SkI2q0RGSRqqYH21fTmMgS4F0RKQa+BLbgjW/0APoBHwF/rfpwY4yphdxFP3ZhjXjfuxrLyc4rIDMrn43bi+jcOoGMvsmkpSRFMNjGraY71qeq6qnAr4CleN1GO4BXgYGqeoeqbvE/TGNMo5G7EF65FJq3gREfHJBAxsxdQ0FRCSlJ8RQUlTBm7hqy8woiGHDjFtJMZaq6Em8dEWOM8c+GBfDKT72rsEa8B0ld9tudmZVPUkIcSQlxAPueM7PyrTUSITbdpTEm7IJ2SZVke5fxtkz2rsJK6nzAcRu3F5GStP8taonxsWzcXnRAXRMeNum+MSasgnVJzfrwTcrHX+rd/zHivaAJBKBz6wQKi0v3KyssLqVza5t8MVIsiRhjwiqwSypGhOOKFvDLDXezrWknbwykmvtAMvomU1BUQkFRCeWq+7Yz+iaH8SswgUJKIiLSU0RmikiWe32su/zXGGMOysbtRSTGez3p3b+fxcXLf8fW5kfy5OH/hMTqk0FaShIjz0glKSGOvIJikhLiGHlGqo2HRFCoYyIvAP8H/Bu8lQ1F5HXgIb8CM8Y0TJ1bJ1BQVMKJhR9x/so/811ib8anPkabxNBm401LSbKkUY+E2p3VXFUXVCorDVrTGGOqkdE3mbRNb5Gx8n5yW/Xj5SOfIL8k3rqkolSoLZHvReQo3Cy5InIFkOdbVMaYBitt3QTSvnuMNW1O4d8d/0xyyyRG2g2DUSvUJDIKbyr1o0VkI7AG+JlvURljGh5VmP0wzHkYjr6I1CvG8nBss0hHZQ5RqDcbrgbOFZEWQIyqFvobljGmQSkvh8zRsODf0O86+MmT0MRuU2sIQvouikhrYBjQDYj1Vq0FVf21b5EZYxqGshKYOgq+mQQn3wqDH4QYu7ugoQj1X4EP8JaoXQKU+xeOMaZBKSmCySPg20wYdB+c/luQYIuVmmgVahKJV9U7fY3EGNOwFBfAhGtg3f/gwn/ACb+IdETGB6EmkVdE5CbgPWBPRWFV64wYYxq+aqdkL8yH166Azcvg8hfhmCuqP5mJWqF2TO4F/g7MBxa5h63sZEwjVe2U7FtXwUuDYWsOXDPREkgDF2pL5E6gu6p+72cwxpjoUNWU7F/Nn0Vazm+9SsPfgy4DIhWiCZNQWyJLgd1+BmKMiR6B819VOKZoAVcs+RU0bQk3zrAE0kiE2hIpAxaLyMfsPyZil/ga0whVzH9V0QLpnf8ug3P+wpYW3Um+8d0aJ1I0DUeoSeQd9zDGGDL6JjNm7hpQ5Zytr3H6+mdZ0SIdrhpPsiWQRiXUO9bHiUhToKcrWqGqJf6FZYypz9JSkhh52uGUTPstx+a/TXaHDLjkadK6dIh0aCbMQr1j/SxgHLAWEOBwERmuqnP9C80YU28V7yBt1k2QPxNOu5O0QffZXeiNVKjdWf8AzlPVFeAtUgVMAGzkzJjGpiAXXr8aNmd7c2ANGB7piEwEhZpE4ioSCICqfisicT7FZIypr/K+9hLInp1w3WTofk6kIzIRFmoSWSgiLwGvuNfX4d1waIxpLL79rzcPVkIbuHE6JPeJdESmHgg1idyMt6bIr/HGROYCz/oVlDGmHlGFBS9A5t2Q3BeufQNapUQ6KlNPhJpEYoF/qerjACLSBLDVZIyJYtXOfVWhrAQ+vAsWjoWeF3jzYDVrGZmATb0U6uUUM4GEgNcJwEfVHSAi8SKyQES+FpGlIvJnV54qIp+LyEoRmeQuHUZEmrnXOW5/t4Bz3ePKV4jI+QHlGa4sR0RGh/i1GNPoVTv3VYXd2+DVn3oJ5NTfwNDXLIGYA4SaROJVdWfFC7fdvIZj9gCDVPU4oB+QISInAY8AT6hqD+AH4EZX/0bgB1XtDjzh6iEivYGhQB8gA3hWRJq41tAzwAVAb+AaV9cYU4PAua9iRPZtZ2blexW2fAsvngPrP4NLn4PBD0BMk8gGbeqlUJPILhHpX/FCRAYARdUdoJ6KxBPnHgoMAqa48nHApW77Evcat/8c8ZZQvASYqKp7VHUNkAMMdI8cVV2tqnuBia6uMaYGwea+SoyPZeP2IsiZCS+eC8U7vEkU+10boShNNAh1TOR2YLKIbHKvU4CrazrItRYWAd3xWg2rgO2qWuqq5AKd3XZnYAOAqpaKSAHQzpV/FnDawGM2VCo/sYo4RgIjAbp27VpT2MY0eJXnvgIoLCrhwt1T4bV/Qoej4dqJ0Np+X0z1Qp325AsRORrohXd11vJQpj1R1TKgn1uj/W0gLVg19xxszUytpjxYK0qDlKGqY4AxAOnp6UHrGNOY7Jv7Cq8FUrR7F0PWPMLAgkzoNQR+OgaaJUY4ShMNQm2JAJwAdHPHHC8iqOr4UA5U1e0iMhs4CWgtIrGuNdIFqGjd5AKHA7kiEgskAdsCyisEHlNVuTGmGmkpSYw8I5XMrHx2blnLTZv+SMed2XDWPXDGXTaFiQlZqHNnvQIcBSzGmxYevP/6q0wiItIBKHEJJAE4F2+w/GPgCrwxjOHAVHfINPd6vts/S1VVRKYBr4vI40AnoAewAK+F0kNEUoGNeIPv1nlrTIjSUpJI25MFi2+BkmIYOgGOHhLpsEyUCbUlkg70VtWD6QpKAca5cZEY4A1VfU9ElgETReQh4CvgJVf/Jby13HPwWiBDAVR1qYi8ASwDSoFRrpsMEbkVmA40Acaq6tKDiM+YxqviBsLp90CbVBjxOnToWfNxxlQioeQFEZkM/FpV8/wPyV/p6em6cKEtD28asZIieP+3sPg16JnhjX/EJ9V8nGm0RGSRqqYH2xdqS6Q9sExEFrD/yoYX10F8xphw2boK3hgO+UvgzLvhzNE2/mEOSahJ5E9+BmGMCYPs9+Cdm72bBq+bAj0GRzoi0wCEeonvHL8DMcb4pKwUZj0An/4LOh0PV46DNkdEOirTQFSbRESkkOD3XgjeTemtfInKGFM3CvNhyg2wbh6k3wAZD0OszZ1q6k61SURV7W4jY6LV6tnw1khv+pLLxsBxNU4yYcxBO5ibDY0x0aCsFOY8DHMfg/Y94Pp3INnmJjX+sCRiTENSsBHe/AWs/x/0+xkMeRSatoh0VKYBsyRiTEOx4kPv6quyEvjpC3DsVZGOyDQClkSMiXale+CjP8Fnz0LHY+HKl6HdUZGOyjQSlkTqsZCWLzWN2+ZsePMm7+bBgb+E8x60q69MWNmtqvVUSMuXmsZLFT4fA2POgsI8uGaiN/5hCcSEmbVE6qnA5UuBfc+ZWfnWGmnsCvNh6ijImcHOroN4LfkuVn7TnM7rv7XWqgk7a4nUU9UuX2oarxUfwnMnw9pP+HAjfeoAABT9SURBVO7UB7mv+R/JK2tlrVUTMZZE6qnOrRMoLC7dr6ywuJTOrRMiFJGJqOIdMO02mDAUWnWCkXOYQAZJzZuSlBBHjMi+lmtmVn6kozWNiCWReiqjbzIFRSUUFJVQrrpvO6NvcqRDM+G2eg48dwp89Sqcejv8YiYcdrS1Vk29YEmknqpYvjQpIY68gmKSEuIYeUaq9Xc3Jnt3wfu/g/EXewPmN0yHwX/eN3hurVVTH9jAej2WlpJkSaOxWjffu3Hwh7Vw0i0w6D5o2ny/Khl9kxkzdw3gtUAKi0spKCrh6hO6RCBg01hZS8SY+mTvLsi8B/5zAWg5jHgfMv52QAIBa62a+sFaIsbUF6tmwbu/ge3rIf1GGPwANGtZ7SHWWjWRZknEmAhbsWYdpR/eQ5/N77Mtvis7Lp5Mt/7nRTosY0JiScSYSFEl99PX6fTxvbQo28FnnUfwUYfhbF0Zw8iUAmthmKhgScSYSNi+AT68iy4rPmBTQi+m9nyKLS170QIojSmxmQlM1LAkYkw4lZXA/GdgziMAvJd8M98eeT0SE7evit3rYaKJJRFjwmXtp/D+nbBlOfS6EC54mJVfFLOjqISkgFs77F4PE03sEl9j/Lbre3j7Znh5COzd7c24e83r0LqrzUxgop61RIzxS1kpLBwLH//Fu//jtDvhjP/b756Pins9AteNufqELjYeYqKGJRFj/LB6Nnw4GrZkQ+qZcMGjcNjRQavavR4mmvnWnSUih4vIxyKSLSJLReQ3rrytiMwQkZXuuY0rFxF5UkRyROQbEekfcK7hrv5KERkeUD5ARJa4Y54UEfHr6zEmJD+shYnXwfhLoGQ3XP0aDJtaZQIxJtr5OSZSCvxWVdOAk4BRItIbGA3MVNUewEz3GuACoId7jASeAy/pAPcDJwIDgfsrEo+rMzLguAwfvx5jqranEGY+AE8P9O48H3QfjFoAaReB/W9jGjDfurNUNQ/Ic9uFIpINdAYuAc5y1cYBs4G7Xfl4VVXgMxFpLSIpru4MVd0GICIzgAwRmQ20UtX5rnw8cCnwoV9fkzEHKCuFL8fB7Idh12Y45ipvpt1WnSIdmTFhEZYxERHpBhwPfA4kuwSDquaJyGGuWmdgQ8Bhua6suvLcIOXG+E8VVnwAM+6HrSvZ3fEEph71MIvKutP5851k9LU7zk3j4PslviLSEngTuF1Vd1RXNUiZ1qI8WAwjRWShiCzcsmVLTSEbU73cRfDyhTDxWgA2nP8i97b+Oytij7Zlak2j42sSEZE4vATymqq+5YrzXTcV7nmzK88FDg84vAuwqYbyLkHKD6CqY1Q1XVXTO3TocGhflGm8Ni+HN4bBi4Pg+2/hwn/ALfOZsvM4W6bWNFq+dWe5K6VeArJV9fGAXdOA4cDD7nlqQPmtIjIRbxC9wHV3TQf+GjCYfh5wj6puE5FCETkJr5tsGPCUX1+Pafiy8wr2u18jo2+y1yW1bY035rHkDYhrDmfcBaf+GpolArBxexEpSfH7ncumLjGNhZ9jIqcC1wNLRGSxK/s9XvJ4Q0RuBNYDV7p9HwBDgBxgN/BzAJcsHgS+cPUeqBhkB24GXgYS8AbUbVDd1Ep2XgFj5q4hKSFuX5fUpJmf85u4t2mz4g2IiYWTR8Gpd0CLdvsd27l1AgVFJSQl/Dj/lU1dYhoLP6/OmkfwcQuAc4LUV2BUFecaC4wNUr4Q6HsIYRoDQGZW/r5uqBZ7NjMobzzHfPc2Qjmk/xzO+B0kdgx6rC1Taxozu2PdGLwuqV7NtnFiznh6b34XoZxlHYYwLel67r3w/GqPtalLTGNmScSYrasYseXvpG3+AJUYlh32E77oMpwN2mG/Lqrq2NQlprGyJGIajCoHxquS9zV8+iQsfYveMXHMa3MpCztfj7bqZF1SxoTIkohpEIINjI+Zu4aRZ6Tun0hUYdVML3msmQNNW8LJo4g5+TY67IxHsvLZZF1SxoTMkoipcwfdIqgDgQPjwL7nfcvMlu6FrDfhf0/B5qWQmALn/hkGjICE1gCkJWJJw5iDZEnE1KmQWwR1rKp7NQq+3wRzp8IXY6FwExzWGy59DvpeAbFNfYvHmMbCkoipUzW2CHxS+V6N5MKl9MmdSJ8fZsKKEjjybLj4Keh+js2qa0wdsiRi6lSk7t7O6JvM2NnLOXbrJ5y89S1Sdi6lOKY5O/pcR9szR0GHnr6+vzGNlSURU6cicvd2/jLSFo/j4bUTaLKngPymXZl15P/R+cyf0+sIm9jZGD9ZEglBJAaKo1XY7t7esxOWvgWLxsHGhdCkKU3SfgL9h5GceibJ1mVlTFiIN9tI45Genq4LFy4MuX7gQHHgH0W/B4qjmW9Jt7wc1s2DbybB0ndg707ocDT0Hw7HXn3AnFbGmLohIotUNT3YPmuJ1CBSA8XRrM7v3s5f6iWOJVNgx0Zomgh9LvWSR5cTbKDcmAiyJFIDm+Y7vCpaMUWbV3H63k85YedM4rcu82bR7X4unPcg9BoCcTZDrjH1gSWRGtg03+GTs3wxq2a+ytU759KpaAUAa+N7E3/qA3Q85Tpo0T7CERpjKrMkUoPGOs13WC4mUIXvlsC3mbBsKt3zs+gObEo8hjndfkNOu0HeJIjEcYclEGPqJUsiNWiM03z7etd5SRGsmQsrPoRvp3t3kQMcfiJTO95Kfufz2BWfsq96oqp1HRpTj1kSCcGhDBRH4+XBdXoxgaq3Hvnq2bBqFqyeA6VF3sSHR50NPe+F7oMhMZnVM771ug4DDreuQ2PqN0siPorUPFKH6pAvJtiR582Qu3q29yjM88rbdIP+w6Dn+dDtNIhttt9hjbXr0JhoZknER9F6efBBXUygCttWw/r5sG6+97xtlbcvoS0ceSYceZb3aNOt2vdtjF2HxkQ7SyI+itbLg6ttERQXeIs5bfoKcr+A9Z/Bri3egQltoOvJMGC4lzSSj4GYmIN6b1sh0JjoYknER5G8PPhQxmIqWgSzv1qBbspmAGs4PnYNSVOyYGvOjxVbH+Hdu9H1JC95tOtx0EnDGBPdLIn4KFJ9/Ac1FqMKu7d5XVBblsPmbNi8jLTN2aTtzP+xXmIn6HQ8HDfUe0453qYZMcZYEgnZ52MgvhW06gytOnnPcfHVHhKpPv7KYzHtm5bQrngTS+dlkdZdvDGMfY81sKfgx4NjE+Cwo70WxmFp0CENOvaFxI6+xmyMiU6WREJRXg7Tfw/lJfuXN2/nJZPEjt52xaNFe2jeHhLakNasJWkDWkDTTtC0xaFN11FeDiW7Ye8ub/LBvbtgzw7YvdU9tsHubZyenUM7KSBx72Za7tlMfNnOH8+xFJAm0LortE2FLunQ9kjv0aGX10UV06T2MRpjGhVLIqGIiYHR62HHJtiR6z0XbPS2CzbCznzIXwa7v4fS4urPJTEQ1wKaxEJMnDcnVJNY7zkmFrQcyku9hFFe+uOjtNhLIDVp2pIe0oqdTZLYHt+VDUnp7Gx6GJulHeUtU7jy7IFeAmkSV/O5jDGmBpZEQtW0ObTv7j2qs3e3l0x2b4WiH1xrYadrObjWw95dUFbiEkQJlJf9+DqmiddSiIn1titexyV4LZl9j5YQ1xyaJbqWTzvvktq4eDZVNX39qanQzq58MsbUHUsida1pc2ja1ftvP0LsfgtjTLhYEmmg7H4LY0w4+HZRv4iMFZHNIpIVUNZWRGaIyEr33MaVi4g8KSI5IvKNiPQPOGa4q79SRIYHlA8QkSXumCdFbGUiY4wJNz/vDHsZyKhUNhqYqao9gJnuNcAFQA/3GAk8B17SAe4HTgQGAvdXJB5XZ2TAcZXfyxhjjM98SyKqOhfYVqn4EmCc2x4HXBpQPl49nwGtRSQFOB+YoarbVPUHYAaQ4fa1UtX56i0SPz7gXMYYY8Ik3HNUJKtqHoB7PsyVdwY2BNTLdWXVlecGKQ9KREaKyEIRWbhly5ZD/iKMMcZ46stER8HGM7QW5UGp6hhVTVfV9A4dOtQyRGOMMZWFO4nku64o3PNmV54LHB5QrwuwqYbyLkHKjTHGhFG4k8g0oOIKq+HA1IDyYe4qrZOAAtfdNR04T0TauAH184Dpbl+hiJzkrsoaFnAuY4wxYSLeuLQPJxaZAJwFtAfy8a6yegd4A+gKrAeuVNVtLhE8jXeF1W7g56q60J3nBuD37rR/UdX/uPJ0vCvAEoAPgds0hC9GRLYA62r5ZbUHvq/lsX6yuA6OxXVwLK6D0xDjOkJVg44F+JZEGiIRWaiq6ZGOozKL6+BYXAfH4jo4jS2u+jKwbowxJgpZEjHGGFNrlkQOzphIB1AFi+vgWFwHx+I6OI0qLhsTMcYYU2vWEjHGGFNrlkSMMcbUmiWRIEQkQ0RWuGnmRwfZ30xEJrn9n4tItzDEdLiIfCwi2SKyVER+E6TOWSJSICKL3eOPfsfl3netm5Z/sYgsDLK/yqn+fYypV8DnsFhEdojI7ZXqhOXzOphlEYIcG3QpBB/j+ruILHffp7dFpHUVx1b7Pfchrj+JyMaA79WQKo6t9nfXh7gmBcS0VkQWV3Gsn59X0L8NYfsZU1V7BDyAJsAq4EigKfA10LtSnVuA5932UGBSGOJKAfq77UTg2yBxnQW8F4HPbC3Qvpr9Q/BuCBXgJODzCHxPv8O7YSrsnxdwBtAfyAooexQY7bZHA48EOa4tsNo9t3HbbXyO6zwg1m0/EiyuUL7nPsT1J+B3IXyfq/3dreu4Ku3/B/DHCHxeQf82hOtnzFoiBxoI5KjqalXdC0zEm6o+UOCU9lOAc9xd975R1TxV/dJtFwLZVDNzcT1T1VT/4XIOsEpVaztTwSHRg1sWIVDQpRD8jEtV/6uqpe7lZ+w/R11YVPF5hSKU311f4nK//1cBE+rq/UJVzd+GsPyMWRI5UFXTzwet437hCoB2YYkOcN1nxwOfB9l9soh8LSIfikifMIWkwH9FZJGIjAyyP5TP1E9DqfqXOxKfF1S9LEKgSH9uN+C1IIOp6Xvuh1tdN9vYKrpmIvl5nQ7kq+rKKvaH5fOq9LchLD9jlkQOFMo08wc1FX1dEpGWwJvA7aq6o9LuL/G6bI4DnsKbqywcTlXV/ngrVI4SkTMq7Y/k59UUuBiYHGR3pD6vUEXyc7sXKAVeq6JKTd/zuvYccBTQD8jD6zqqLGKfF3AN1bdCfP+8avjbUOVhQcoO6jOzJHKgqqafD1pHRGKBJGrX/D4oIhKH90Pymqq+VXm/qu5Q1Z1u+wMgTkTa+x2Xqm5yz5uBt/G6FQKF8pn65QLgS1XNr7wjUp+XU9WyCIEi8rm5wdWLgOvUdZxXFsL3vE6par6qlqlqOfBCFe8Xqc8rFvgpMKmqOn5/XlX8bQjLz5glkQN9AfQQkVT3X+xQvKnqAwVOaX8FMKuqX7a64vpcXwKyVfXxKup0rBibEZGBeN/frT7H1UJEEiu28QZmsypVq2qq/3Co8j/ESHxeAapaFiFQ0KUQ/AxKRDKAu4GLVXV3FXVC+Z7XdVyBY2iXVfF+ofzu+uFcYLmq5gbb6ffnVc3fhvD8jPlxtUC0P/CuJvoW70qPe13ZA3i/WADxeN0jOcAC4MgwxHQaXjPzG2CxewwBfgX8ytW5FViKd1XKZ8ApYYjrSPd+X7v3rvi8AuMS4Bn3eS4B0sP0fWyOlxSSAsrC/nnhJbE8oATvP78b8cbQZgIr3XNbVzcdeDHg2Bvcz1kO3hIJfseVg9dHXvEzVnEVYifgg+q+5z7H9Yr72fkG749jSuW43OsDfnf9jMuVv1zxMxVQN5yfV1V/G8LyM2bTnhhjjKk1684yxhhTa5ZEjDHG1JolEWOMMbVmScQYY0ytWRIxxhhTa5ZETIMmImWy/2y+wWZlPktE3qvj9z1LRE4JeP0rERlWR+dOqet4D+K9J4pIj0i8t6mfYiMdgDE+K1LVfhF437OAncD/AFT1+To89514d237QkSaqGpZFbufA+4CbvLr/U10sZaIaZTcuhPLRWQe3pQVFeV/EpHfBbzOcpPaISLD3ASAX4vIK67sJ+KtKfOViHwkIsmu/q+AO1zr5/TA84pIPxH5TH5cs6ONK58tIo+IyAIR+VZETq8i/MuBTHfMJyKyL0mKyKcicqy7S3qsiHzhYrvE7e/mjvnSPU5x5WeJtybF68ASd/z77mvNEpGr3Vt8ApzrpvowxpKIafASKnVnXS0i8Xj/yf8Eb/bVjjWdRLwZfu8FBqk3YWPFomDzgJNU9Xi8qcfvUtW1wPPAE6raT1U/qXS68cDdqnos3l3Y9wfsi1XVgcDtlcor4kgFflDVPa7oRWCE29cTaKaq37hYZ6nqCcDZwN/dlBubgcHqTQZ4NfBkwOkH4t1N3RtvOvBNqnqcqvbFJS315q7KAY6r6TMzjYMlEdPQFbk/5BWPScDRwBpVXanelA2vhnCeQcAUVf0eQFUrJtzsAkwXkSXA/wHVTicvIklAa1Wd44rG4S12VKFi8rxFQLcgp0gBtgS8ngxc5CbguwFvCg7w5kAaLd5Ke7PxpurpCsQBL7h4J+MtXlRhgaqucdtL8Focj4jI6apaEFBvM960HsZYEjGNVlXz/ZSy/+9FvHuWKo55CnhaVY8BfhlQv7YqWhhlBB+zLAp8D/UmSZyBtwDRVcDrAfFeHpA8u6pqNnAHkI/XkkjHWwGwwq6A834LDMBLJn+T/ZcOjndxGGNJxDRKy4FUETnKvb4mYN9avCVQEW8t+FRXPhO4SkTauX1tXXkSsNFtB65PXYi3VOl+3H/0PwSMd1wPzKlcrxrfcmAL5UW8bqkvAlpI04HbAmYpPj4g3jzXLXU93pKyBxCRTsBuVX0VeAz3mTg98SYSNMaSiGnwKo+JPKyqxcBI4H03sB64bO6bQFvXDXQz3h9tVHUp8Bdgjoh8DVRMuf0nYLKIfAJ8H3Ced4HLKgbWK8U0HG+M4hu8RZYeCPWLUdVdwCoR6R5QtgjYAfwnoOqDeF1X34hIlnsN8CwwXEQ+w0sGuwjuGGCB+xzuBR4CEJFkvC7CcE3lb+o5m8XXmCgjIpcBA1T1D+51J7xxj6NdC8PP974D2KGqL/n5PiZ6WEvEmCijqm/jdbvhbmD8HO+qKl8TiLMd72IAYwBriRhjjDkE1hIxxhhTa5ZEjDHG1JolEWOMMbVmScQYY0ytWRIxxhhTa/8PjrbCd+1msLEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot mean income in each age group\n", "grouped = gss.groupby('educ')\n", "mean_income_by_educ = grouped['realinc'].mean()\n", "plt.plot(mean_income_by_educ, 'o', alpha=0.5)\n", "\n", "# Plot the predictions\n", "pred = results.predict(df)\n", "plt.plot(df['educ'], pred, label='Age 30')\n", "\n", "# Label axes\n", "plt.xlabel('Education (years)')\n", "plt.ylabel('Income (1986 $)')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Logistic regression" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Predicting a binary variable" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: 0.588510\n", " Iterations 6\n", "Intercept -1.685223\n", "C(sex)[T.2] -0.384611\n", "age -0.034756\n", "age2 0.000192\n", "educ 0.221860\n", "educ2 -0.004163\n", "dtype: float64\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hUZfbA8e+ZFJKQAoQAIfQmoXdcQASpgoAdESzrsljWsrZd61r3t+rqomsFlRV7wYZYEBQEUZFQpdcAIRACgSRAes7vjzuJkzBJJpNpCe/nee6TKbe8w6733LedV1QVwzAMwyjP5u8CGIZhGIHJBAjDMAzDKRMgDMMwDKdMgDAMwzCcMgHCMAzDcCrY3wXwlMaNG2ubNm38XQzDMIxaZfXq1UdUNc7Zd3UmQLRp04akpCR/F8MwDKNWEZG9FX1nmpgMwzAMp0yAMAzDMJwyAcIwDMNwqs70QRiGYRQUFJCSkkJubq6/ixJwwsLCaNGiBSEhIS4fYwKEYRh1RkpKClFRUbRp0wYR8XdxAoaqcvToUVJSUmjbtq3Lx5kmJsMw6ozc3FxiY2NNcChHRIiNja12zcoECMMw6hQTHJxz59/FqwFCRMaKyDYR2Ski91Sy36UioiLSz/6+jYjkiMg6+/aK1wpZXATfPgDHKhwKbBiGcUbyWoAQkSDgReB8oAswRUS6ONkvCrgVWFnuq12q2su+3eCtcpKxB9a8Ca+PgtR1XruMYRhnBhHhqquuKn1fWFhIXFwcF1xwQaXHLV26tMp9fM2bNYgBwE5V3a2q+cD7wCQn+z0GPAX4Z9hB4w5w3UKwhcD/xsGOxX4phmEYdUP9+vXZuHEjOTk5ACxatIiEhAQ/l8o93gwQCcB+h/cp9s9KiUhvoKWqLnByfFsRWSsiP4jIOc4uICIzRCRJRJLS09PdL2mTRJi+GGLbwbuXWzUKwzAMN51//vl8+eWXALz33ntMmTKl9Ltff/2VQYMG0bt3bwYNGsS2bdtOO/7kyZNcd9119O/fn969e/P555/7rOyOvDnM1VmPSOn6piJiA2YC1zrZ7yDQSlWPikhf4DMR6aqqWWVOpjobmA3Qr1+/mq2dGh0Pf/waPrwa5t8CmQdg2D1gOrwMo3b6+h449Jtnz9msO5z/RJW7XXHFFTz66KNccMEFbNiwgeuuu47ly5cD0LlzZ5YtW0ZwcDCLFy/mvvvu4+OPPy5z/D//+U/OO+885syZw/HjxxkwYAAjR46kfv36nv09VfBmgEgBWjq8bwGkOryPAroBS+29682A+SIyUVWTgDwAVV0tIruAToB3s/HVi4IrP4QvboMfnoDMFJjwLARVPrFky8FMvtmYxoHjOSQ0CGdst6Ykxsd4taiGYQSuHj16kJyczHvvvce4cePKfJeZmck111zDjh07EBEKCgpOO/7bb79l/vz5PP3004A1fHffvn0kJib6pPwlvBkgVgEdRaQtcAC4Ariy5EtVzQQal7wXkaXAXaqaJCJxQIaqFolIO6AjsNuLZf1dUAhMehFiWlpBIusAXD4Xwpzf8LcczGT2sj3EhIcQHxNGZk4Bs5ftYcbQtiZIGIY/ufCk700TJ07krrvuYunSpRw9erT08wcffJDhw4fz6aefkpyczLBhw047VlX5+OOPOeuss3xY4tN5rQ9CVQuBm4GFwBbgQ1XdJCKPisjEKg4fCmwQkfXAPOAGVc3wVllPIwLD77UCRfJyeH0MHN/ndNdvNqYREx5CTHgINpHS199sTCvdZ8vBTGYu2s5dH61n5qLtbDmY6atfYhiGn1x33XX84x//oHv37mU+z8zMLO20fuONN5weO2bMGJ5//nlUrZbztWvXerWsFfHqPAhV/UpVO6lqe1X9p/2zf6jqfCf7DrM3LaGqH6tqV1Xtqap9VPULb5azQr2nwbSPISsVXh0BB1aftsuB4zlEhZWtiEWFBXPguDWCoaSGkZlTUKaGYYKEYdRtLVq04Lbbbjvt87/97W/ce++9DB48mKKiIqfHPvjggxQUFNCjRw+6devGgw8+6O3iOiUlEaq269evn3ptwaDDW+Hdy+BEOlzyKiROKP1q5qLtZOYUEBP+ez9FyfvbR3Wq8nvDMDxny5YtPm+nr02c/fuIyGpV7edsf5NqwxVNOsP076BpF/jgKvjpBbAH1rHdmpKZU0BmTgHFqqWvx3ZrClRdwzAMwwhUJkC4KrIJXLPAqj18ez98eQcUFZAYH8OMoW2JCQ/hYGYuMeEhZTqoExqEk51bWOZU2bmFJDQI98evMAzDcJlJ910doRFw2Vz47mFY8Rxk7IbL3iAxvmGFI5bGdmvK7GV7AKvmkJ1bSGZOAZP7t/BhwQ3DMKrP1CCqy2aDUY/aRzitgNdGwtFdFe5eVQ3DMAwjUJkahLt6T4NG7eD9qfDqeXD5m9DuXKe7JsbHmIBgGEatY2oQVah0DkPrQfDn7yGqGbx9MST9z38FNQzD8DCXAoSIJIjIIBEZWrJ5u2CBwKU5DI3awp++hXbDYMFfrfwvRYUVndIwjDouKCiIXr16lW7Jycleu9Ybb7zBzTff7LXzV9nEJCJPApOBzUDJrA4FlnmtVAHCcZY0UPr3m41pZZuMwmJgygfWwkMrX4b0rXDpHIho5I9iG4bhR+Hh4axbVzfWlnGlBnEhcJaqjlPVCfatqlQZdUK15jAEBVu5XyY+D8k/Wv0Sh7f4qKSGYQSyoqIi7r77bvr370+PHj2YNWsWYC0SdO6553L55ZfTqVMn7rnnHt555x0GDBhA9+7d2bXLGgDzxRdfMHDgQHr37s3IkSNJS0s77Rrp6elccskl9O/fn/79+7NixYoal9uVTurdQAj27KpnkoQG4afNgq5yDkOfq6HxWfDBNGuE08WvQudxFe9vGIZXPPnrk2zN2OrRc3Zu1Jm/D/h7pfvk5OTQq1cvANq2bcunn37K66+/TkxMDKtWrSIvL4/BgwczevRoANavX8+WLVto1KgR7dq1Y/r06fz6668899xzPP/88zz77LMMGTKEX375BRHhtdde46mnnuKZZ54pc93bbruN22+/nSFDhrBv3z7GjBnDli01e0h1JUCcAtaJyHc4BAlVvbVGVw5Qjqm7Q4OEtKw8aBRRvTkMrQbCjCXWCKf3r4Tz7odz7jJrSxjGGcBZE9O3337Lhg0bmDdvHmAl7NuxYwehoaH079+f+Ph4ANq3b18aOLp3786SJUsASElJYfLkyRw8eJD8/Hzatm172nUXL17M5s2bS99nZWWRnZ1NVFSU27/FlQAx377VeeVTd2fnFlKsSkFhEQczrZrD5P4tXBuyGtMCrvvGWnzo+8fh0Ea48CUI9e2CH4ZxpqrqSd+XVJXnn3+eMWPGlPl86dKl1KtXr/S9zWYrfW+z2SgstAa83HLLLdxxxx1MnDiRpUuX8vDDD592jeLiYn7++WfCwz2XpaHKAKGqc0UkFGvBHoBtqnr6Chd1gLNO6dax9U9LrOfyAkEh4VYTU9NusPhhOLIDrnjbmj9hGMYZY8yYMbz88sucd955hISEsH379mqtU+2YInzu3LlO9xk9ejQvvPACd999NwDr1q0rbepyV5Wd1CIyDNgBvAi8BGyvq8NcXemUrnb6bhEY8leYNs9afGj2MNixyIu/wjCMQDN9+nS6dOlCnz596NatG9dff31p7cAVDz/8MJdddhnnnHMOjRs3drrPf//7X5KSkujRowddunThlVdeqXG5q0z3LSKrgStVdZv9fSfgPVXtW+Ore5An0n27kpq7Rum7M/ZY2WDTNsLw++GcO63UHYZheIRJ9105b6T7DikJDgCquh1rVFOdU1Xqbqhh+u6SSXXdL4Ulj1sjnXKzPP0zDMMwPMKVTuokEXkdeMv+fipw+tJqdUBJYj3H/oXyndJuDX11FBph9Usk9IWF91vzJa54hy2FzVzr1zAMw/ARVwLEjcBfgFsBwZpB/ZI3C+VPVSXW80j6bhE4+0Zo1h0+upbiWcNY1eQuMpuNLdOvYbK+Gkb1qSpihpSfxp3VQ6tsYlLVPFX9j6perKoXqepMVXVp0pyIjBWRbSKyU0TuqWS/S0VERaSfw2f32o/bJiJjKjrW1zyavrvNELh+GQfDO3J16qNMOvgswVpYOpLqm42nz5Y0DKNiYWFhHD161K2bYV2mqhw9epSwsLBqHVdhDUJEPlTVy0XkN6zcS+Uv2KOyE4tIENbIp1FACrBKROar6uZy+0Vh1U5WOnzWBbgC6Ao0BxaLSCdVdb7Ct495NH13dHOeazmTy4+9Sr+D79EsezMLOv+L4rCmZllSw6imFi1akJKSQnp6ur+LEnDCwsJo0aJ6C5VV1sR0m/3vBW6WZwCwU1V3A4jI+8AkrKR/jh4DngLucvhsEvC+vaayR0R22s/3s5tlCWjxDaNYEHYLh6J7MmrnY0xbN415bR7G1nSIv4tmGLVKSEiI01nGhnsqbGJS1YP2lzep6l7HDbjJhXMnAPsd3qfYPyslIr2Blqq6oLrH2o+fISJJIpJUm58YSkZPJdUfyjs93iArOJapO+9gWu67UBwQlSbDMM5ArgxzHeXks/NdOM5ZL1FpU5WI2ICZwJ3VPbb0A9XZqtpPVfvFxcW5UKTA5NivsTm/KR/2eoPUVhOIWz2THc+MZNaXP1U8Ec8wDMNLKuuDuBGrptBORDY4fBUFuJJHNgVo6fC+BZBa7jzdgKX2EQfNgPkiMtGFY+scx36NLQczeeb4HQxr1Y1xKc8wde1U3kx7AM6/1IxqMgzDZyqcSS0iMUBD4F+A4wikbFXNqPLEIsHAdmAEcABYhTUje1MF+y8F7lLVJBHpCryL1e/QHPgO6FhZJ7UnZlIHCsfZ2rGndjF+6700yknm+ybX8FuHG0jJzDdzJQzD8IjKZlJXWINQ1UwgE5hiP0kTIAyIFJFIVd1X2UVVtVBEbgYWAkHAHFXdJCKPAkmqWmGGWPt+H2J1aBcCfwmUEUyeUlnCvwPHc4iPsYajHY1oz7s95zJo6/8x4vAbtDu1gYWdH+dQTrCZK2EYhle5kotpAvAfrCf5w0BrYIuqdvV+8VxXm2oQjmnFHSfbldzsneV7WrrtMMNzF/OXUy9TEBTOwg4PsT58gGs5oAzDMCpQ01xMjwNnA9tVtS1Wk1HN17Krw7YczGTmou3c9dF6Zi7afloHs2NacZvIaRPjnOWEOnaqgLR2F/Fuz7mcDInloi1/ZcLB/3Iow3ReG4bhHa4EiAJVPQrYRMSmqkuAmiUZr8NcSQdeVcI/Z7O1z+kQS73gYDIi2vFezzdYGz+Zvgff5859N1nrTBiGYXiYK7mYjotIJFYOpndE5DBWv4DhhLNFh0o+L+krcCXhX/nZ2iWBByAqLJTP429jdVBPrj3yb5g1FM5/CnpPM8uaGobhMa7UICZhrUt9O/ANsAuY4M1C1WaupAN3Ja14ec5qFfEDLmROt7fZEdoZ5t9M6utTeOnrpAqbtgzDMKrDlQBxB5CgqoWqOldV/wtc4uVy1VoJDcLJzi1bwXJWO3An4V9ifAy3j+rE05f1ZGy3pizekk5qcUO+6PEi8xv/mSYp3zJl9RUMYGPVK90ZhmFUwZVRTIeBI1hDTZfYP1ujqn18UD6XBcoopqpGKHlK+ZFOP+8+StPszdyb8wzxRamsbn4lXzf5M5H165tRToZhVKimo5gOAGOBJ0Tk7pJzeqpwdY1H04FXonxT1oncQg7V78xNUc+yodnF9Et9h+u3T6fo0EaPXtcwjDOHK53UqOo+ETkXeFlEPgJcXD7tzOTRdOAVKN/RHRkWTFZOATHhUXzf7h72NBzCyB2P8tfd18PPB2HgjWb9a8MwqsWVO0YSgKrmquofgaVAqDcLZVStfEd3s6h6nMwrpFl0PYpVWRc+kCfb/o+cVufCwvvgzYlwbK+/i20YRi1SZR9EbREofRC+VD5dR6em9dmedrLs+0MnSEiex6S0Fwm2QdDYf0Gfq81wWMMwgMr7ICpL1lejFeV87UwMEJUp31luy9zPBXsep+OptdBhFEx8HqLj/V1MwzD8zK1kfdR8RTnDj8pP2KNBK94963nOyfyc85JfgJfOhnFPQ/dLTW3CMAynqlxRrvxqcg6ryhkBzNmEvcjwUL4KnwA3roDGneCT6fDBNMhO81MpDcMIZBUGCBHJFpEsJ1u2iGT5spBG9VU6YS+2PVz3DYx6FHYsghcHwLr3oI70RxmG4RmV1SCiVDXayRalqtG+LKS3Ldq7iNzCXH8Xw6OqTOdhC4LBt1m1ibiz4LMb4N3LIfOAfwtuGEbAcHkUk8OCQYA1N8JbhXKHu53UuzN3c+FnF9IltgvPDX+OpvUrzodU21S2KFEZxUXw66vw3SNgC4bRj0Gfa1zum3D5OoZhBBy3RjE5HDwReIY6vGDQ9/u+597l9xIREsGzw5+lZ1xPD5eulsjYDfNvheTl0HYoXPCs1RxVCV+lFjEMwztqmmrjMer4gkHntTqPt8e9TVhQGH/85o98tvMzfxfJPxq1g6vnwwUzIXUdvDwIlv8HigoqPKSqxY8Mw6i9vLpgkIiMFZFtIrJTRO5x8v0NIvKbiKwTkR9FpIv98zYikmP/fJ2IvFKtX+WGjg078t749+jTpA8PrniQJ399ksLiM3DZC5sN+l0Hf1kJHUdZzU6zh0HKaqe7u5Le3DCM2smVAFF+waDncGHBIBEJAl4Ezge6AFNKAoCDd1W1u6r2Ap7CWvu6xC5V7WXfbnDlx9RUg7AGvDLqFaYmTuXtLW9z0+KbyMw7Q9NlRzeHyW/D5HfgVAa8NgK+vgfyssvs5kp6c8MwaidXFwzKofoLBg0AdqrqblXNB963n6uUqjoOl62PkxnbvhZsC+aeAffw6KBHWZW2iisWXMG2jG3+Lpb/JF5g1Sb6T4eVr8CLZ8OWL0qHxLqz+JFhGLVDlQFCVU+qapHjgkH2JqeqJAD7Hd6n2D8rQ0T+IiK7sGoQtzp81VZE1orIDyJyjgvX86iLOl7E/8b8j/yifK76+iq+3vO1r4sQOMKiYfzTcN1CCIuxJte9OxmOJfssvblhGL7nyiimbE5/ss/EyvJ6p6ruruC4y4Axqjrd/v4qYICq3lLB/lfa979GROoBkap6VET6Ap8BXcvVOBCRGcAMgFatWvXdu9fzE7yP5BzhzqV3subwGq7ucjW3972dYJtLWdLrpqICWDkLlvwfaDEMvQsG3QrBJsGvYdRGNR3F9B/gbqyn/xbAXcCrWE1Gcyo5LgVo6fC+BZBayf7vAxcCqGpeSS1FVVdjNWudtiyaqs5W1X6q2i8uLs6Fn1J9jcMb89ro15jSeQpvbn6TGYtmcDTHlQpUHRUUAoNuhptXWZ3Y3z8GrwyGPcv8XTLDMDzMlRrESlUdWO6zX1T1bBFZr6pOJw2ISDCwHWtY7AFgFXClqm5y2Kejqu6wv54APKSq/UQkDshQ1SIRaQcsB7qrakZF5fRFNtf5u+bz6M+P0qBeA2YOm0n3uO5evZ4/VHvS245F8OWdcHwvdL0YRj8OMae1JBqGEaBqWoMoFpHLRcRm3y53+K7C6KKqhcDNwEJgC/Chqm4SkUftk+8AbhaRTSKyDrgDuMb++VBgg4isB+YBN1QWHHxlYvuJvHn+mwRJENd8cw0fbP2AurKeBvw+6S0zp4D4mDAycwqYvWwPWw5mnrbfzEXbueuj9cxMbs3WSxbBuffAtq/ghX6w/BkozPPTrzAMw1NcqUG0A54D/oAVEH7BGtF0AOirqj96u5Cu8OV6EMdzj3Pvj/fy44EfGdd2HA/94SEiQiJ8cm1vmrloe5llTIHS97ePslr4Kp05HXYMFt4PWxdYk+7GPgmdRvvr5xiG4YIa1SDsw1QnqGpjVY2zv96pqjmBEhx8rUFYA14c8SK39L6Fb5K/YcqXU9h93GlffcAp8/S/aHuZ2oErk94qnTndsA1c8Q5M+wQkCN69zBrtdHSXr36eYRgeVGWAEJFOIvKdiGy0v+8hIg94v2iBzSY2ZvSYwexRszmed5wrvryCr3Z/5e9iVaqqJiRXJr25NHO6wwi48ScY9Rgk/wgvDrRqFjnHvffjAkxlgdgwagtX+iBeBe4FCgBUdQNwhTcLVZsMjB/IRxM+IrFRIn9f/nce/+Vx8ooCs/29qrxJrkx6c3nmdHAoDL4VblkDPa+An1+E5/vAqtegqG6nMHG1L8cwAp0rASJCVX8t91nd/i+8mppENOG1Ma/xx65/5INtHzDtq2kkZyb7u1inqerp35VJb9WeOR3VFCa9ANcvg7hEa8TTK0Ng53de+53+ZhIYGnWFKwHiiIi0xz5iSUQuBQ56tVS1UIgthDv63cGLI17k0MlDTF4wmQW7F/i7WGW48vSfGB/D7aM68fRlPbl9VKfThri6PXM6vgdcu8DK71SYA29fDG9fCmmbPfb7AoVJYGjUFa6OYpoNDAKOAXuAaaqa7PXSVYMvRzFV5dDJQ/x92d9Zc3gNF3a4kHsH3BsQo5wCZu2GwjxrNvbyp63kf72mwvD7rASBdYAro8EMI1DUaMEgh5PUB2yqml3lzn4QSAECoLC4kJfXv8yrG16lbUxbnj73aTo27OjvYgXW6m+nMmDZ0/DrbGsluz/8xVoGNax2r2gbMIHYMFzgVoAQkTsqO6mq/qey730t0AJEiV8O/sI9y+7hRMEJ7u53N5efdTni4lKeZ4xjyfDdY7BxHkQ0hnP/Dn2v9Up+J18FyIAKxIZRCXcDxEOVnVRVH/FA2TwmUAMEWAn/HvjxAVakrmBYy2E8MugRGoU18nexAs+B1fDtP2Dvj9CgFQy7D3pcDrYgj5zePNkbxuk80sQU6AI5QAAUazHvbHmHmatnElMvhn8O+SeDmg/yd7ECj6o1wum7R+DQBmvk03kPQOfxUMOal+kbMIzT1TQXk+EBNrFxVZereG/8e8SExnD9ouv596p/k1+U7++iBRYR6DgSZvwAl/4Pigvgg6nw2kjY/UONTm1GFxlG9ZgA4WNnNTqL9y94n8lnTebNzW8y9aup7Druv1QUATvj12aDbhfDTSthwn8h+yC8ORHeuACSV7h1SrM8qmFUjwkQfhAWHMYDZz/A8+c9T9rJNC7/4nLe2vwWxVrs03LUihm/QcHQ9xprRvbYJ+HIdnhjHMydCPtWVutUZnlUw6geV+ZBOBvNlAmsVtV1XimVGwK9D6IiR3KO8PBPD/NDyg8MaDaAxwc/TnxkvE+uXSvb5PNPQdIcWPEsnEyH9iOsORQtnDahnsaMLjKMsmrUSS0i7wL9gC/sH43HWvynM/CRqj7lwbK6rbYGCABV5dOdn/Lkr09iExv3DLiHie0nen047F0frSc+Jgybw3WKVTmYmcvTlzldBypw5J+08jr9+CzkZFiBYujd0PoP/i6ZYdQqNe2kjgX6qOqdqnonVrCIw1rU51qPlfIMJiJc3PFi5k2cR6eGnXhgxQPcvvR2MnK9u0ZSrW6TD61vTar7628w8mE4uB7+Nxb+Nx52L7VGQxmGUSOuBIhWgONQmwKgtarmAIGZtrSWahnVkjlj5nBH3ztYlrKMiz6/iG+Tv/Xa9Vxtkw/YjmyAepEw5HYrUIz5FxzdCW9OgtdHw/ZvTaAwjBpwpYnpQeAi4HP7RxOA+cAzwGxVnerVErqoNjcxObPj2A4eWPEAm49uZkybMdw38D6vTK6rqk2+1k0uK8iFdW9bTU+Z+6FZdxj8V+hyodXhfYYzfTBGeTWeKCci/YDBgAA/qmrA3YnrWoAAKCgu4I2Nb/DS+peIDo3m/oH3M7qNb5fwrJUd2QCF+fDbh1agOLoDGrS21qfoNRVCakETmhfUumBv+IQnAkQQ0BQofQRT1X0uHDcWaz3rIOA1VX2i3Pc3AH8BioATwAxV3Wz/7l7gT/bvblXVhZVdqy4GiBI7ju3gwRUPsunoJka3Hs39Z9/vs1Qd7nRkB9RTanExbPsKfpwJB5KgfhwMvB76T4fwhh65hKd+r7f/3WptsDe8qkad1CJyC5AGLAIWAF/a/1Z1XBDwInA+0AWYIiJdyu32rqp2V9VewFPAf+zHdsFata4rMBZ4yX6+M1LHhh15e9zb3NbnNpbsX8KFn13Igt0L8EWalOp2ZAfc3AqbDRIvgOmL4dovIb4XfP84/KcrfPU3yNhTo9N76vf64t/NzCQ3qsuVTurbgLNUtauq9rDf0Hu4cNwAYKeq7lbVfOB9YJLjDqqa5fC2PvZFiez7va+qeaq6B9hpP98ZK9gWzPTu0/nwgg9pGdWSe5ffy43f3UjqiVSvXre6k8sCdjU1EWgzBKbNgxt+hC4TrfkUz/eBD6ZZk+7cCLie+r2unKemgwVq9ag1wy9cCRD7sSbGVVeC/dgSKfbPyhCRv4jILqwaxK3VPHaGiCSJSFJ6erobRQzwETpOdGjYgTfPf5N7BtzDmrQ1XPj5hbyz5R2Kiou8cr3qriBXK55Sm3WHi16xRj4N/ivsWQ5zRlv5njZ+AkUFLp/KU7+3qvN4ooZhZpIb1eVKgNgNLBWRe0XkjpLNheOczfI67RFNVV9U1fbA34EHqnnsbFXtp6r94uLiXChSWQHXHOKiIFsQUxOn8tmkz+jbtC9P/PoEV399NTuO7fDK9apahtRRrXpKjY6HkQ/BHZth3NNw6ijM+yM82wOW/RtOHqnyFJ76vVWdxxM1FbeXizXOWK4EiH1Y/Q+hQJTDVpUUoKXD+xZAZe0h7wMXunmsWwK2OcRFzSOb89KIl3jinCfYn72fy7+4nGdXP0tOof+e1mvlU2pofRjwZ7hlNUx5H+LOsvdTdIFPb4TUijPKeOr3VnUeT9VUqhPsDcNr60GISDCwHRgBHMBKz3Glqm5y2Kejqu6wv54APKSq/USkK/AuVr9Dc+A7oKOqVtiO4s4oplqdaqKcY7nHeCbpGT7f9TkJkQncN/A+hrYY6peyBNQoJncd3mothbr+fSg4CS0HQr8/QZdJEBJWZldfjGIyI5AMb3F3RblnVfWvIvIFzpt3Jrpw4XHAs1jDXOeo6j9F5FEgSVXni8hzwEis2dnHgJtLAoiI3A9cBxQCf1XVryu7ljsBoi7+R7fq0Coe/+VxdvEphw8AACAASURBVGfuZlTrUfyt/99oVr+Zv4vlE14JTDnHYd07Vt6njN0QEWvNpej3R2jUzjMFd4GZw2B4i7sBoq+qrhaRc519r6o1W73Fw9wJEHX1P7qCogLe2PQGszbMIkiCuKX3LVzR+QqCbbVrJnF1bvhe/9+yuBj2LLVGPm39CrTIShDY/0/QcYxPZmkHcs0skMtmVM7tiXL2uQdzVXWatwrnKe5OlKvL/8fen72ff678JysOrKBTw07cN/A++jbt6+9iuaS6N3yf1gazUmHNm7B6LmSnQmQz6HUl9LnKp7WKQFFXH7TOFDVN970QmGCfyxCw6vJM6ppQVRbvW8xTq57i0MlDjG83njv73klcRPVHfXlT+UB9JDuXkOAgl2/4fulPKiqEHQthzVvWXy2GNudAn2sgccJpfRWeFigPN3WxqfZMUlmAcKVenAysEJH5wMmSD1X1P54pnuFNIsKo1qMY3Hwwr298nf9t/B9L9i3hpl43cWXilYTYQqo+iRc43txCg4S0rDxaNoooHW68fOdRBndoBPxevspG7SQ0CD/tJuX14bVBwdB5vLVlpVp9FWvegk+mQ1gD6H4Z9JoCzfuw5VCWR2/mjk/tjkO0/fHUfuB4DvExZYNhwM19MdziSoBItW82XBveagSgiJAIbul9C5PaT+LJVU/ydNLTfLLjE/4+4O8Maj7I69evLCAs255Odm4hzWLqYRNrqHHDiBA2p2bT5Kzfb/CV3fDHdmvK7GVW2gzHZo7J/Vt4/bcBEN3cWrBoyJ2QvMxqglr7Fqx6lbyGHdkZch7FTcYRHxPvkZu54xBtoPTvNxvTfB4g/BKcDZ+och6Eqj6iqo9g5Ul6xuG9UQu1im7FiyNe5IXzXiC/KJ/rF13PLd/dQnJmsteuWX5C4qbULPYcOUlBURE2EQqKlMh6Qew8XFpBJTE+imOnXJ9fEDCTwGw2aDcMLp0Dd26DC57laGEYEw7P4vaNF3Hx5tsYcOJ7GtcrqtF8m0CasV4r574YLnGlD6Ib8BZQkj70CHC143yGQGD6IKovvyift7e8zewNs8krzGNK4hRu6HkD0aHRHr1O+TbqRZvTCLZBeGgwZ7eL5efdR8nLLyS/WBndxRqSm5lTQH5hEXFRYX5vY6+puz5aT5fQNLqmf03i4S+Jzk8j3xbOb1Hn0PeC662AUs1RUIHW7h8o/SFG9dW0k/on4H5VXWJ/Pwz4P1X1frtENZgA4b4jOUd4Ye0LfLLjExrUa8Bfev2FSzpd4rFhseU7kMsHhPTsXFbuziAyLJihneICfhRMdW+GZW7mWkyLrLW0O/gVXY8vIawo20pB3vViq8+iRT8rsaALZTAjhwxPqGmAWK+qPav6zN9MgKi5rRlbefLXJ0lKS6JDgw7c3vd2zkk4B3HhhlWZ8k+7zgLC3qMnaR4TRl6RevQJ1NNPtu7cmCs65vrBzemc/au1sNG2b6AoD2JaQdcLoetF0Lx3pcHCPLUbnlDTAPEpsAarmQlgGtBPVS+s+CjfMwHCM1SV7/Z9x8zVM9mXvY8BzQZwZ7876RJbfikP1zm7QboTEKp7Q3T3KdsbKS+qLHtupjUBb9OnsOt7KC6Ahm2sQNH1ImjWw6WahWFUV00DREPgEWAIVpbVH4BHVPWYpwtaEyZAeFZBUQEfbv+QWetncSzvGOPbjefW3rfSPLK5W+er6dOuOzd7d27mVV3HJ/Mtco7B1i+t1OO7l1qzthu2seZWJE6EhH5WZ7hheEBN50G0UdVbq97NqEtCgkKYmjiVie0nMmfjHN7a/BbfJn/L1MSp/Knbn2gQ1qBa50uMj6lR84c7wzrdGZ9f1XV8MqQzvCH0nmZtJ4/C1gXW9ssr8NPz1sztxAusgNF6MAT5Zy6LUfe58hjyHxHZKiKP2bOsGmeQqNAobutzGwsuWsD5bc9n7qa5nP/J+cxaP4tTBad8Vg53hnW6s1ZDVdfx+ZDO+rHQ9xqY+hH8bRdc8jq0HADr3oU3J8FT7WHen+C3eVZiQcPwIJfSfYtIM+ByYDIQDXygqo97uWzVYpqYfGPHsR08v/Z5luxfQqOwRszoMYPLOl1GaFCoV6/rjeYiV6+zJ/0Eh7LzaNUogoQG4XRqWp/taSf92zmcfwp2L4FtX8H2hXAyHWzB0HoQdDofOo2B2PYuncp0dp/ZatQHUe5E3YG/AZNV1bt3hGoyAcK31qev57k1z7Hq0Cqa12/OTb1uYny78V7LGOuNDmdXrrP3yEnW7j9On1YNaBVbPzCHkxYXw4HVVrDY9jWkb7E+b9QeOo6GjiOh9RCnuaHMcFmjpp3UiVg1h0uBo1grv32sqoc9XdCaMAHC91SVn1N/5rm1z7H56GbaRLfhhp43MLbNWIJsQR6/nq+edB2vsy/jFPHR9WjTOLL0+4BPRJexB3Ysgh3fQvJyKMyFkAhoOxQ6jIT251lZZ0UCbsKd4Xs1DRC/AO8BH6mqx5f99BQTIPynZGjsS+tfYsexHbSLaceNPW9kdJvR2KR2j7ap9asOFuRA8o9WsNjxLRxLtj6PaQXth/Pm4XZkxw8iP/T3QQe16vcZNeaxJqZAZgKE/xVrMYv3LualdS+xK3MXHRp04KZeNzGi1YhaGyjq1BO2qrUq3u4lsGsJ7FkGeVkoQlpkZ/bH9Gd/TD+2hHQlIjK69v0+wy01rUF0BP4FdAFKGzFVNaBWRjEBInAUFRfx7d5veWndSyRnJdOhQQdm9JjB6NajvdL05E11uo2+qJDkDcvYumI+iTlraHFqE0FaSBFB5DXrTUSn4db6Fi36Q2iEv0treElNA8SPwEPATGAC8Ef7cQ95uqA1YQJE4CkqLuKb5G94dcOr7MrcRZvoNkzvPp1x7cb5bR0Kd9T1UT4lvy89I4O+sp1zQ7fQOH0lpK61FkGyhVhpP1oPsuZdtBoIYXXn95/pahogVqtqXxH5TVW72z9brqrnuHDhscBzQBDwmqo+Ue77O4DpQCGQDlynqnvt3xUBv9l33aeqEyu7lgkQgatYi/lu33fM3jCbrRlbSYhM4Lpu13Fhhwu9PjzWqIHcTNj3C+z9ydpS10BxISDQrBu0PBtaDoSW/aFBa5MKpJaqaYBYAZwDzAO+Bw4AT6jqWVUcFwRsB0YBKcAqYIqqbnbYZziwUlVPiciNwDBVnWz/7oSqRjo5tVMmQAQ+VWVZyjJmbZjFb0d+Iy48jqu6XMVlnS4jMtTl/6kNf8k/BQeSfg8YB1ZD/gnru8im1gS+FgOsJqn4nqZZqpaoaYDoD2wBGgCPYU2U+7eq/lLFcX8AHlbVMfb39wKo6r8q2L838IKqDra/NwGijlJVfjn4C69vfJ2VB1cSFRLF5M6TmZo4lcbhjV0+T11v+gl4xUVweDPsXwn7f7X+loySkiBo0gUS+kBCX2uL61ztdS8M73MrQIjIW6p6lYjcpqrPuXHRS4Gxqjrd/v4qYKCq3lzB/i8Ah0pmaItIIbAOq/npCVX9zMkxM4AZAK1ateq7d+/e6hbT8LNNRzYxZ+McFu1dRIgthEkdJnFt12tpFd2q0uPqdOdxDfg9aJ44bNUsDqyGA2usv7n2FCDB4VbTVHzP37e4RAg2zYz+5G6A2AycD8wHhmFlci2lqhlVXPQyYEy5ADFAVW9xsu804GbgXFXNs3/WXFVTRaQdVtPWCFXdVdH1TA2idtubtZc3Nr3B5zs/p7C4kOEth3NN12vo3aS30/Uo6tTw0wr4Kr25p8tRRsnQ2gOrIXUdHFwPhzZAXpb1vS0EmnSGpt2t4NG0q/W6fqzb5TWqx90AcStwI9AOq9/B8b9SrWqYq6tNTCIyEngeKzg4nZ0tIm8AC1R1XkXXMwGibjiSc4T3tr7HB9s+IDMvk26x3bi669WMbD2yzMinWj+BrQq+Sm/ujXJUqbgYju2xgkXJlrYJTjr85x/ZzAoYcZ2hSaJV04jrBPWi3LumUaGa9kG8rKo3unHRYKxO6hFYAWYVcKXjWtb2fod5WE1ROxw+bwicUtU8EWkM/AxMcuzgLs8EiLolpzCHL3Z9wVub3yI5K5lm9ZsxtfNULup4ETH1Yup8DcKd3+csaKZl5bA+JZOuzWPcWpjJp6lGThyGtI1WsDhk/3tku7XSXomYllbQiDsLYjtA444Q2xEim5hRVG6q0XoQ7gQH+3GFInIzsBBrmOscVd0kIo8CSao6H/g3EAl8ZG9GKBnOmgjMEpFirJTkT1QWHIy6Jzw4nMvPupxLO13K8pTlzN08l2dWP8NL619ifLvxnN16It+stW4Ijk+2k/u38HPJPcOdtSzKr1WRnp3Lqj3HiAwLJj4mjMycAmYv2+Py8qjxMWGs23eczFP5RIYF0zgyzKVyuC2yCUSeZ+WKKlFcZHV8H95iJSFM3waHt/6eY6pEvRho3MFKUNionX1ra/2NiDXBw00m1YZRa2zL2Ma7W9/ly91fkleUR9eGfWjKSIJyu9KyYaRbHbJ+79StgCfSmy/bnk52biF/aN+o9OZe1TnKX/fn3UfJsr8/u12sS+fwieJiyEqBIzvg6E6rpnFkh5WoMHM/4HBfC42CRm2suRoNWkODVmW3sGh//YqA4G4fRL2SDuPawASIM8fx3ON8vONjPtj2AQdPHiS+fjyXdbqMizpeVO1hsv4aCVVVYPJEevPNqVn0bBlNk6jfF0iqqp+mfDNVenYua/Yep6C4mHHd42vHaLHCPDi21+rnyNhtbceS4fg+ayu/0FW9aIhOgOjmEJMA0S2s11HxENXM2sIb1dllXt0NEGtUtU/JcFevltADTIA48xQWF/LD/h94b+t7rDy0kmBbMCNajWDyWZPp17Sf09FPjvzVj+Hqzb+mtRt3fp8rCyYFSi3LLapw6igc3/t7wMg8AFn2LfNA2c7yErZgazJgVDOo3wTqN4b6cdYWaX8fEWstFxveqFZNEnS3DyJURK4BBonIxeW/VNVPPFVAw3BHsC2YEa1HMKL1CPZk7uGj7R/x+c7PWZi8kLYxbbm80+VMaD+BmHqeW7PaE1xdX7um63iP7daU2cv2AK730zg7JijIxiMTu9TeoOBIxH5zb2xN3nOmMA+yUuFEGmQfsrYThyA7DbIPQmaKlafqZDpokfNzBIf9HizCG1i5q+pFW81Z9aKt92HREBoJofUdtkhr7Y6QCAiuByHh4McEl5XVIIYAU7GWGp1f7mtV1eu8XLZqMTUIAyC3MJeFyQv5cNuHbDiygVBbKCNbj+SSjpfQr1m/MmnH/VWD8OUQXXdqIYHaLxNwioutSYAnj1jB4tRRyDkGORlwKsP+2r7lZkFepv1vlpUE0VW2ECvghIRBUD1rNnpQqPV5kH1r2hUmVHs+M+BmDUJVfwR+FJEkVX3drSsbho+FBYcxqcMkJnWYxNaMrXyy4xMW7F7AV3u+omVUSy7qcBGTOkyiSUQTt56wPaH8aCOA7NxCEhqEV3KUe9yphdS05nLGsNkgopG1xVXjgULVymGVm2X1h+SfgPyT1paXbf0tzLUWeyrMg8IcKMi1/hYVWFtxwe+vi/KtYOEFrsyDCAVuAIbaP/oBeEVVC7xSIjeZGoRRkdzCXBbvW8wnOz5h1aFVBEkQgxMGM6n9JJqF9Oa7zcd8+rRs0oQYgaSmE+VeA0KAufaPrgKKSlJoBAoTIAxX7Mvaxyc7PuGLXV9wOOcw0aHRnN/2fCa1n0S3xt2q7Nj2FNOM4xnm37Hmahog1qtqz6o+8zcTIIzqKCouYuXBlXy+63O+2/cdeUV5tI1py8T2ExnXdhzNI5v7u4hnnEDJPXWmqWmAWANcVpIoz548b56q9vF4SWvABAjDXdn52Xyb/C2f7/qctYfXAtCnSR/GtxvP6NajaRDWwM8lrPsCJffUmahGqTaAu4ElIrIbK2Ffa6xlRw2jTogKjeKSTpdwSadLSMlO4as9X/Hl7i957JfH+Nev/2JI8yGMbzeeoS2GEhFSe8a31yauDv115K9hymcSV3IxfSciHYGzsALE1to0w9owqqNFVAtm9JjBn7v/ma0ZW/ly95d8vedrlqYsJTw4nKEthjK2zViGJAwhLDis6hMaLvFE7inw3miwM5VLyzvZA8IGL5fFMAKGiJAYm0hibCK3972dNYfX8M2eb1i8bzELkxcSERzBsJbDGNtmLIMSBlEvqJ6/i1yruXOz99cw5TOJSdZnGNVQWFzIqkOrWJi8kMX7FpOZl0lEcARDWwxlZOuRnJNwjmmGcoMnck+ZUUzuqVEndW1hAoThawXFBaw8uJLFexezZP8SMnIzqBdUj0HNBzGq9SiGthhaYZoP43TmZu8fNR3F9DEwB/hatTrzw33LBAjDn4qKi1hzeA3f7fuOxXsXk3YqjSAJok/TPgxvOZxhLYfRMqqlv4tpGKepaYAYiTVq6WzgI+ANVd3q8VLWkAkQRqAo1mI2HdnEkv1LWLJ/CTuP7wSgQ4MODG85nHNbnku32G4E+TEJm2GU8EgTk4jEAFOA+4H9wKvA24GScsMECCNQ7c/az9KUpSzdv5TVaasp0iIa1mvI4ITBnJNwDoMTBpumKMNvahwgRCQWmIaVZiMVeAcYAnRX1WGeK6r7TIAwaoPMvEx+Sv2JZSnL+PHAjxzPO45NbPSM68k5CecwKGEQiY0Sy2SdNQxvqmkT0ydAZ+AtrOalgw7fJVV0Yvv3Y4HnsNakfk1Vnyj3/R3AdKAQSAeuU9W99u+uAR6w7/q4qs6lEiZAGLVNUXERG49uZFnKMpanLGdLxhYAGoU14uz4sxnUfBCDmg8iLiLOzyU16rKaBohxqvpVuc+qXI5URIKA7cAoIAVYBUxR1c0O+wwHVqrqKRG5ERimqpNFpBGQBPTDWlx2NdBXVY9VdD0TIIza7kjOEX5O/ZkVqSv4OfVnMnIzAOjUsBMD4wdydvzZ9G3al/oh9f1cUqMujbiqcS6m8nmXnH3m5Lg/AA+r6hj7+3sBVPVfFezfG3hBVQeLyBSsYHG9/btZwFJVfa+i65kAYdQlxVrMtoxtrEhdwS+pv7D28Fryi/MJlmC6Ne7GwPiBDIwfSM+4noQGhfq7uGeUupYk0K1cTCLSDEgAwu0375I8yNGAKzOBErA6s0ukAAMr2f9PwNeVHJvgwjUNo06wia10Jvf07tPJLcxlXfo6Vh5cycqDK3n1t1eZtWEWobZQejbpSb+m/ejXtB894nqYFCAOvPGk707eqNqqslQbY4BrgRbAfxw+zwbuc+HczhLrO62uiMg0rOakc6tzrIjMAGYAtGrVyoUiGUbtFBYcxtnxZ3N2/NmAlYE26VASSWlJrDq0ilkbZvGyvkyILYTujbvTt2lfejfpTc8mPYkOjfZoWWrLMqaOT/rxMWFk5hQwe9meGj/pn0lJAitbcnQuMFdELlHVj904dwrgODOoBdYIqDLs8yzuB8516NdIAYaVO3apkzLOBmaD1cTkRhkNo1aKCo1ieKvhDG81HLACxtrDa1l1aBVJh5KYs3EORVqEIHRs2JHeTXrTu0lvejXpRfP6zd1eGMmdm25Fx4xMjGN72kmvBQ1PPemXD271goTs3MIzIklgZU1M01T1baCNfbRRGar6HyeHOVoFdBSRtsAB4ArgynLX6A3MAsaq6mGHrxYC/yciDe3vRwP3VvVjDONMFRUaxdAWQxnawloZ+FTBKX478htrDq9h3eF1fLHrCz7Y9gEAjcMb06NxD3o26UmPxj3o2rgr4cGu3dzcuek6OybjRB7Pf7+Ls9vFevTp3pEnnvSdBbfUzFxsItAoos4nCaysialkqESkOydW1UIRuRnrZh8EzFHVTSLyKJCkqvOBf9vP/5H9iWafqk5U1QwReQwryAA8qqoZ7pTDMM5EESERpR3ZYCUZ3HFsB+vT17M+fT0b0jfw/f7vAQiSIDo17ETXxl3pFtuNbo270b5Be4Jtp98e3LnpOjvmUFYuBUXFXm3H90Q6cGfBrXVsffILi4gJDymtVUzu36LO9T9A5U1Ms+x/H3H35PbhsV+V++wfDq9HVnLsHKwcUIZh1FCwLbi00/uKzlcAcCz3GL8d+Y11h9fx25HfWJi8kHnb5wEQFhRG50ad6da4m3Vco0TaxrR166br7JiMkwXE1i87+qp8oKlpv4Un0oFXFBAPZhZWumpdXRkGW1kT038rO1BVb/V8cQzD8JWGYQ3LNEupKvuy97HxyEY2HtnIpqObmLd9HrlFuQDUC6pHy/rtyDkZT5O8djQPb0dwUXNO5gZVetN1dqMOtslpN17HQOOJDubE+BhmDG1b5kZd3Sd9dwKitzrH/aGyJqbVPiuFYRh+JyK0jm5N6+jWjG83HrBmeydnJbP56Ga2ZmxlS8YWUk/9yoFT37P2FIAQH5HAq9sS6XS4E2c1PIsODTuQlR3Jt5vSS2/M5TukbxnRnsVb0snMKXD6dO+pDubE+Jga3ZTdqYXUpWGwVY1iMgzjDBZkC6J9g/a0b9CeCe0nAFZNI+VECtsztrP9mLVty9jGor2LSo+zEUpMUAtiQ1uTeqw5a9LjmXH22dzWtntpnql2cZEVPt0HylBSd2ohgVJ2T6isielZVf2riHyBkzkIqjrRqyUzDCMgiQgto1rSMqolI1qPKP38VMEpth/bzis//URabjKn9AAH8taSU7wEgJuWQ/jP4bSJbkObmDa0jWlLl47tGBfdhlbRLcuMpAqk9aarWwsJpLLXVGVNTG/Z/z7ti4IYhlG7RYRE0KtJLyILhI6NRlhDQYG84myO5u9jX/Ye+nQoYE/WHjakb+DrPV+XOb5pRFNaR7emVXQrwiKbsCM9lCZFCTSLaE5OXpDPhpIGQud4oHA13XcoVkZXBbapar63C1ZdJheTYQSGmYu2n/YEXfLeceRPTmEOe7P2kpyZzN6svezL3kdyVjL7svZxPO94mXOG2xrQKrolnRq1oUVUC5pHNqd5/eY0j2xO0/pNCbGF4AmeyrNUm0YxuZWLyeHg8cArwC6sFBhtReR6Vf268iMNwzgTufoEHR4cTudGnencqPNp58jMy2Rf1j5STqSQkp1S+jcpLYkFuxegDq3eNrERFx5H88jmNKvfjGb1m9E0ommZ143CGrm0xkagdI4HiioDBPAMMFxVdwKISHvgS35PrGcYhlHKE8NLY+rF0D2uO93jup/2XUFRAYdOHiL1ZCqpJ1J//3sild/Sf2Px3sUUFJdd6DLYFkxceJy1RcTROLwxTSKaEBceR2x4LLFhscSGx7L/WBYJDaLKHFtbO5g9wZUAcbgkONjtBg5XtLNhGIY3n6BDgkJoGd2SltEtnX5frMUcyz3GoVOHSDuZxqGTh0g7lcaRnCOkn0pnb9ZeVh1aRVZ+lvPz59Yn3BZDWFA0YbYobMVRRIfGMHfTOhrUa0BMvRhi6sUQHRpd+reuplyvbBTTxfaXm0TkK+BDrD6Iy/g9BYZhGEZAsYnNqhWEx9I1tmuF++UV5XEk5whHc45aW+5RtqYfYMWevagtmyJOcLwgjdyiHewvOMGGpMIKzxUeHE5USBSRoZFEhkYSFRJF/ZD6RIVGERkSSURIBBHBEUSERBAeHF76Piw4jLCgMOoF1yM8KJx6wfUICwojNCjUaaoTX6usBBMcXqfxeyrudKDh6bsbhmHUHvWC6pEQmUBCpMNSM51gSzuHDubGVgdz52bRnCw4ybHcY2TmZ5KVl1Xmb2ZeJicKTpCdn82J/BNk5Wdx4MQBThac5ETBCXIKq99EZRMbobZQQoJCCLWFlgaNYFswQRJU5m/Hhh156A8PefBfx1LZRLk/evxqhmEYXuSJ0UMVNY+V1A5a4rxpqzLFWkxuYS4nC05yqvAUpwpOcarwFLmFueQW5ZJXmEduUW7p+/yifPKL8ikoLrBeF1vvC4sLKdIiioqLKCwuJDM3l7SsHDKOZTLzxHaPj5ZyZcnRMKzV3roCpdMDVfU6j5XCA8wwV8M4s7k7RLU2DUl15KkhuZUNc6163Jc1Ya4Z1gpzP2At3pPt8tUNwzB8wHGIqk2k9PU3G9MqPKbkJpuZU1Amsd6Wg5k+LLl73Pm91eVKL0gHVb1MRCap6lwReRdrjQfDMAyfqepJ350cSLU5sZ4vcj65UoMoGVB8XES6ATFAG4+VwDAMowquPOknNAgnO7fsSKOqciAdOJ5DVFjZ5+TaMu/Bnd9bXa4EiNn2pT8fBOYDm4EnPVYCwzCMKrjSnDK2W1MycwrIzCmgWLX09dhuTSs8ry9ust7izu+trioDhKq+pqrHVPUHVW2nqk1KVpszDMPwBVee9EtmcMeEh3AwM5eY8JAqO2x9cZP1Fnd+b3W5kospFngYGIw1UW458JiqHvVYKQzDMCrhagrt6s7g9kRaEH/yds4nVzqp3weWAZfY308FPgAqXE+6hIiMBZ4DgoDXVPWJct8PBZ4FegBXqOo8h++KgN/sb/eZ9ScM48zlzRTadSWxnje40gfRSFUfU9U99u1xoEFVB4lIEPAicD7QBZgiIl3K7bYPuBZ418kpclS1l30zwcEwzmC+aE4xTudKDWKJiFyBlYsJ4FKsbK5VGQDsVNXdACLyPjAJq5MbAFVNtn9XXI0yG4ZxBjJP+r5XWbK+bKw+BwHuAN62f2UDTgBVJf5IAPY7vE8BBlajbGEikgQUAk+o6mdOyjgDmAHQqlWrapzaMAzDt2rjjO0Km5hUNUpVo+1/baoabN9sqhrtwrnF2WmrUbZW9unfVwLP2tehKF/G2araT1X7xcXFVePUhmEYvlNbZ2y7lE9WRCYCQ+1vl6rqAhcOS4EyWa1aAKmuFkxVU+1/d4vIUqA31qp2hmEYtUptnbFdZSe1iDwB3IbVd7AZuM3+WVVWAR1FpK19TesrsCbaVUlEGopIPfvrxlhDbDdXfpRhGEZgqq0ztl2pQYwDeqlqMYCIzAXWAvdUdpCqForIzVh5m4KAOaq6SUQeBZJUIWdleAAACjJJREFUdb6I9Ac+xVpfYoKIPKKqXYFEYJa989qG1QdhAoRhGLWSq/M4Ao2rSxY1ADLsr12uD6nqV8BX5T77h8PrVVhNT+WP+wk4fTFawzCMWsib8zi8yZUA8S9grYgswep4Hgrc69VSGYZheIk/RhPV1hnblS4YJCKC9YRfCPTHChArVfWQb4rnOrNgkGEYVfHUIjt1SWULBlVag1BVFZHPVLUvLnYwG4ZhBKraOprIX1xpYvpFRPrb+wsMwzBqLV8ssuNNvm4ecyUX03CsILFLRDaIyG8issFrJTIMw/CS2rz+gz8m27lSgzjfa1c3DMPwodo0mqh8beFIdq7Pm8cqrEGISJiI/BW4GxgLHFDVvSWbV0pjGIbhRbUlK6yz2sLynUfJKyxb+/F281hlNYi5WOtRL+f3lN23ea0khmEYPlAbssI660xvGBHC5tRsmpz1e3OYt5vHKuuD6KKq0+zLi14KnOO1UhiGYRilnKXmSIyP4tgp3y6PWlmAKCh5oaqFlexnGIZheJCzzvSwkGCGdIj1afNYZU1MPUUky/5agHD7e8GaIuFKym/DMAzDBY6d0qFBQlpWHjSK8OuEvgoDhKoG+awUhmEYZzDHGd7xMWFk5xZSrEpBYREHMwv9lprD1WR9hmEYhpc465RuHVufmPAQbh/VyW/lcmWinGEYhuFFgbpehAkQhmEYfhaoM7xNgDAMw/Czsd2alg5b9dUQVleYAGEYhuFngTrD23RSG4ZhBIBAnOFtAoRhGIYf+GNlu+ryahOTiIwVkW0islNE7nHy/VARWSMihSJyabnvrhGRHfbtGm+W0zAMw5f8kbrbHV4LECISBLzI74n+pohIl3K77QOuBd4td2wj4CFgIDAAeEhEGnqrrIZhGL7kOO/BJlL6+puNaf4uWhnebGIaAOxU1d0AIvI+MAnYXLKDqibbvysud+wYYJGqZti/X4SVcvw9L5bXMAzDLdVtLqotK9t5s4kpAdjv8D7F/pnHjhWRGSKSJCJJ6enpbhfUMAzDXe40FwXqvIfyvBkgxMln6sljVXW2qvZT1X5xcXHVKpxhGIYnuNNcFKjzHsrzZoBIAVo6vG8BpPrgWMMwDJ9xJ01GoM57KM+bfRCrgP9v715jpKzuOI5/fwVvoBXBSyhQwYTghVaE1YI0xGrrLQZvmLD2hS+0fdOmYExMTWMTNW1i0rQ2TdpovCQ2jW29VM2mkRJEo74AF2WVdaVqpLoqgqIY26a1+u+Lc0bH7bPuyszwnKG/TzKZmTMzyy/Pc2b+POeZOWeupDnAa8BK4NJxvnYN8NOmE9NnAte0P6KZWWtmTDmI3f/84OOJ9mB8w0Ul/u5hpI4dQeRFhr5P+rAfAv4YEYOSrpe0HEDSyZKGgUuAmyUN5tfuAm4gFZkngesbJ6zNzErSLcNFe0IR4z0tULaenp7o7++vO4aZ/R/qhh+9jUbSpojoqXrMv6Q2M2tRNwwX7QlP1mdmZpVcIMzMrJILhJmZVXKBMDOzSi4QZmZWaZ/5mqukncDfOvhPHA681cG/307O2hnO2hndkrVbcsLny3p0RFTOVbTPFIhOk9Q/2neFS+OsneGsndEtWbslJ7Qvq4eYzMyskguEmZlVcoEYv1vqDvA5OGtnOGtndEvWbskJbcrqcxBmZlbJRxBmZlbJBcLMzCq5QFSQNEvSeklDkgYlrcrtUyWtlfRCvj5srL+1F7IeKGmjpIGc9brcPkfShpz1D5L2rzsrgKQJkp6W1Jfvl5pzm6RnJW2W1J/bitv/AJKmSLpH0vO5zy4pMaukeXl7Ni7vSVpdYlYASVfm99QWSXfl91qp/XVVzjkoaXVua3m7ukBU+w9wVUQcBywGvifpeOCHwLqImAusy/fr9i/g9Ig4EVgAnC1pMXAj8Iuc9R3g8hozNltFWkCqodScAN+IiAVN3ycvcf8D/BJ4KCKOBU4kbd/iskbE1rw9FwCLgH8Af6LArJJmAD8AeiJiPjCBtCpmcf1V0nzgO8AppP1/nqS5tGO7RoQvY1yAB4BvAVuB6bltOrC17mwjck4CngK+RvoV5cTcvgRYU0C+mbmjng70ASoxZ86yDTh8RFtx+x/4IvAy+QsnJWcdke9M4IlSswIzgFeBqaR1c/qAs0rsr6QVOW9tun8tcHU7tquPIMYgaTZwErABOCoi3gDI10fWl+wTedhmM7ADWAu8BLwbadlXgGFSh6/bTaSO+1G+P40ycwIE8BdJmyR9N7eVuP+PAXYCd+Shu1slTabMrM1WAnfl28VljYjXgJ8BrwBvALuBTZTZX7cAyyRNkzQJOBeYRRu2qwvEZ5B0MHAvsDoi3qs7z2gi4sNIh+0zSYeZx1U9be+m+jRJ5wE7ImJTc3PFU0v53vXSiFgInEMaYlxWd6BRTAQWAr+JiJOAv1PAEM1nyeP2y4G7684ymjxefz4wB/gSMJnUF0aqvb9GxBBp6Gst8BAwQBomb5kLxCgk7UcqDr+LiPty85uSpufHp5P+x16MiHgXeIR03mSKpMaSsjOB1+vKlS0FlkvaBvyeNMx0E+XlBCAiXs/XO0jj5KdQ5v4fBoYjYkO+fw+pYJSYteEc4KmIeDPfLzHrN4GXI2JnRHwA3AecSrn99baIWBgRy4BdwAu0Ybu6QFSQJOA2YCgift700IPAZfn2ZaRzE7WSdISkKfn2QaSOPQSsB1bkp9WeNSKuiYiZETGbNLzwcER8m8JyAkiaLOmQxm3SePkWCtz/EbEdeFXSvNx0BvAcBWZt0ssnw0tQZtZXgMWSJuXPg8Z2La6/Akg6Ml9/GbiItH1b3651n2Ap8QJ8nXTo+AywOV/OJY2ZryNV53XA1AKyfhV4OmfdAvw4tx8DbAReJB3KH1B31qbMpwF9pebMmQbyZRD4UW4vbv/nXAuA/twH7gcOKzjrJOBt4NCmtlKzXgc8n99XvwUOKLG/5qyPkQrYAHBGu7arp9owM7NKHmIyM7NKLhBmZlbJBcLMzCq5QJiZWSUXCDMzq+QCYdYGki6UFJKOrTuLWbu4QJi1Ry/wOOlHgGb7BBcIsxblObuWkqZ+XpnbviDp13l+/j5Jf5a0Ij+2SNKjeSLANY3pEMxK4wJh1roLSOsx/BXYJWkhabqD2cBXgCtIU0M35vj6FbAiIhYBtwM/qSO02Vgmjv0UMxtDL2niQUgTEfYC+wF3R8RHwHZJ6/Pj84D5wNo0xQ8TSNNJmxXHBcKsBZKmkWamnS8pSB/4QZoBtvIlwGBELNlLEc32mIeYzFqzArgzIo6OiNkRMYu0wttbwMX5XMRRpAkKIa3ydYSkj4ecJJ1QR3CzsbhAmLWml/89WriXtMjMMGkm0JtJKxLujoh/k4rKjZIGSDMFn7r34pqNn2dzNesQSQdHxPt5GGojaZW67XXnMhsvn4Mw65y+vJjT/sANLg7WbXwEYWZmlXwOwszMKrlAmJlZJRcIMzOr5AJhZmaVXCDMzKzSfwFMmMsChxTx1QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Recode grass\n", "gss['grass'].replace(2, 0, inplace=True)\n", "\n", "# Run logistic regression\n", "results = smf.logit('grass ~ age + age2 + educ + educ2 + C(sex)', data=gss).fit()\n", "print(results.params)\n", "\n", "# Make a DataFrame with a range of ages\n", "df = pd.DataFrame()\n", "df['age'] = np.linspace(18, 89)\n", "df['age2'] = df['age'] ** 2\n", "\n", "# Set the education level to 12\n", "df['educ'] = 12\n", "df['educ2'] = df['educ'] ** 2\n", "\n", "# Generate predictions for men and women\n", "df['sex'] = 1\n", "pred1 = results.predict(df)\n", "\n", "df['sex'] = 2\n", "pred2 = results.predict(df)\n", "\n", "grouped = gss.groupby('age')\n", "favor_by_age = grouped['grass'].mean()\n", "plt.plot(favor_by_age, 'o', alpha=0.5)\n", "\n", "plt.plot(df['age'], pred1, label='Male')\n", "plt.plot(df['age'], pred2, label='Female')\n", "\n", "plt.xlabel('Age')\n", "plt.ylabel('Probability of favoring legalization')\n", "plt.legend()\n", "plt.savefig('../images/brfss-logreg.png')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }