String_Processing/Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Source Information\n",
    "---\n",
    "**Created by**:\n",
    "    \n",
    "**Updated by**: October 01, 2024 by Gloria Seo\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goal\n",
    "This tutorial outlines the steps to perform regression analysis using the diabetes dataset in Python. By using scikit-learn, you can create, fit, and evaluate a linear regression model, allowing you to understand the relationships between various factors and disease progression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we'll explore regression models using Python's scikit learn (sklearn) package and the built data set. Please keep in mind that although regression is considered one of the simplest or most basic machine learning techniques, a thorough understanding of the assumptions and limitations is essential for a correct interpretation of the results. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start by loading the linear_model and diabetes data set from sklearn. Note that we're only loading the components that we need for this exercise since the entire sklearn package is extremely large."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Required Modules for the Jupyter Notebook\n",
    "\n",
    "Before running the notebook, make sure to load the following modules.\n",
    "\n",
    "**Module: linear_model, mean_squared_error, r2_score, numpy, matplotlib, pandas, sklearn, scipy** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import linear_model\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "from sklearn.datasets import load_diabetes\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the Diabetes Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, load the diabetes dataset using the load_diabetes() function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "diabetes = load_diabetes()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the diabetes data set. We are interested in how the disease progression depends on factors such as age, sex, BMI (body mass index) and blood pressure. Note that these factors have been mean-centered and scaled by the standard deviation.\n",
    "\n",
    "The disease progression is the *dependent* variable and age, sex, BMI etc. are the *independent* variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'data': array([[ 0.03807591,  0.05068012,  0.06169621, ..., -0.00259226,\n",
       "          0.01990842, -0.01764613],\n",
       "        [-0.00188202, -0.04464164, -0.05147406, ..., -0.03949338,\n",
       "         -0.06832974, -0.09220405],\n",
       "        [ 0.08529891,  0.05068012,  0.04445121, ..., -0.00259226,\n",
       "          0.00286377, -0.02593034],\n",
       "        ...,\n",
       "        [ 0.04170844,  0.05068012, -0.01590626, ..., -0.01107952,\n",
       "         -0.04687948,  0.01549073],\n",
       "        [-0.04547248, -0.04464164,  0.03906215, ...,  0.02655962,\n",
       "          0.04452837, -0.02593034],\n",
       "        [-0.04547248, -0.04464164, -0.0730303 , ..., -0.03949338,\n",
       "         -0.00421986,  0.00306441]]),\n",
       " 'target': array([151.,  75., 141., 206., 135.,  97., 138.,  63., 110., 310., 101.,\n",
       "         69., 179., 185., 118., 171., 166., 144.,  97., 168.,  68.,  49.,\n",
       "         68., 245., 184., 202., 137.,  85., 131., 283., 129.,  59., 341.,\n",
       "         87.,  65., 102., 265., 276., 252.,  90., 100.,  55.,  61.,  92.,\n",
       "        259.,  53., 190., 142.,  75., 142., 155., 225.,  59., 104., 182.,\n",
       "        128.,  52.,  37., 170., 170.,  61., 144.,  52., 128.,  71., 163.,\n",
       "        150.,  97., 160., 178.,  48., 270., 202., 111.,  85.,  42., 170.,\n",
       "        200., 252., 113., 143.,  51.,  52., 210.,  65., 141.,  55., 134.,\n",
       "         42., 111.,  98., 164.,  48.,  96.,  90., 162., 150., 279.,  92.,\n",
       "         83., 128., 102., 302., 198.,  95.,  53., 134., 144., 232.,  81.,\n",
       "        104.,  59., 246., 297., 258., 229., 275., 281., 179., 200., 200.,\n",
       "        173., 180.,  84., 121., 161.,  99., 109., 115., 268., 274., 158.,\n",
       "        107.,  83., 103., 272.,  85., 280., 336., 281., 118., 317., 235.,\n",
       "         60., 174., 259., 178., 128.,  96., 126., 288.,  88., 292.,  71.,\n",
       "        197., 186.,  25.,  84.,  96., 195.,  53., 217., 172., 131., 214.,\n",
       "         59.,  70., 220., 268., 152.,  47.,  74., 295., 101., 151., 127.,\n",
       "        237., 225.,  81., 151., 107.,  64., 138., 185., 265., 101., 137.,\n",
       "        143., 141.,  79., 292., 178.,  91., 116.,  86., 122.,  72., 129.,\n",
       "        142.,  90., 158.,  39., 196., 222., 277.,  99., 196., 202., 155.,\n",
       "         77., 191.,  70.,  73.,  49.,  65., 263., 248., 296., 214., 185.,\n",
       "         78.,  93., 252., 150.,  77., 208.,  77., 108., 160.,  53., 220.,\n",
       "        154., 259.,  90., 246., 124.,  67.,  72., 257., 262., 275., 177.,\n",
       "         71.,  47., 187., 125.,  78.,  51., 258., 215., 303., 243.,  91.,\n",
       "        150., 310., 153., 346.,  63.,  89.,  50.,  39., 103., 308., 116.,\n",
       "        145.,  74.,  45., 115., 264.,  87., 202., 127., 182., 241.,  66.,\n",
       "         94., 283.,  64., 102., 200., 265.,  94., 230., 181., 156., 233.,\n",
       "         60., 219.,  80.,  68., 332., 248.,  84., 200.,  55.,  85.,  89.,\n",
       "         31., 129.,  83., 275.,  65., 198., 236., 253., 124.,  44., 172.,\n",
       "        114., 142., 109., 180., 144., 163., 147.,  97., 220., 190., 109.,\n",
       "        191., 122., 230., 242., 248., 249., 192., 131., 237.,  78., 135.,\n",
       "        244., 199., 270., 164.,  72.,  96., 306.,  91., 214.,  95., 216.,\n",
       "        263., 178., 113., 200., 139., 139.,  88., 148.,  88., 243.,  71.,\n",
       "         77., 109., 272.,  60.,  54., 221.,  90., 311., 281., 182., 321.,\n",
       "         58., 262., 206., 233., 242., 123., 167.,  63., 197.,  71., 168.,\n",
       "        140., 217., 121., 235., 245.,  40.,  52., 104., 132.,  88.,  69.,\n",
       "        219.,  72., 201., 110.,  51., 277.,  63., 118.,  69., 273., 258.,\n",
       "         43., 198., 242., 232., 175.,  93., 168., 275., 293., 281.,  72.,\n",
       "        140., 189., 181., 209., 136., 261., 113., 131., 174., 257.,  55.,\n",
       "         84.,  42., 146., 212., 233.,  91., 111., 152., 120.,  67., 310.,\n",
       "         94., 183.,  66., 173.,  72.,  49.,  64.,  48., 178., 104., 132.,\n",
       "        220.,  57.]),\n",
       " 'frame': None,\n",
       " 'DESCR': '.. _diabetes_dataset:\\n\\nDiabetes dataset\\n----------------\\n\\nTen baseline variables, age, sex, body mass index, average blood\\npressure, and six blood serum measurements were obtained for each of n =\\n442 diabetes patients, as well as the response of interest, a\\nquantitative measure of disease progression one year after baseline.\\n\\n**Data Set Characteristics:**\\n\\n  :Number of Instances: 442\\n\\n  :Number of Attributes: First 10 columns are numeric predictive values\\n\\n  :Target: Column 11 is a quantitative measure of disease progression one year after baseline\\n\\n  :Attribute Information:\\n      - age     age in years\\n      - sex\\n      - bmi     body mass index\\n      - bp      average blood pressure\\n      - s1      tc, T-Cells (a type of white blood cells)\\n      - s2      ldl, low-density lipoproteins\\n      - s3      hdl, high-density lipoproteins\\n      - s4      tch, thyroid stimulating hormone\\n      - s5      ltg, lamotrigine\\n      - s6      glu, blood sugar level\\n\\nNote: Each of these 10 feature variables have been mean centered and scaled by the standard deviation times `n_samples` (i.e. the sum of squares of each column totals 1).\\n\\nSource URL:\\nhttps://www4.stat.ncsu.edu/~boos/var.select/diabetes.html\\n\\nFor more information see:\\nBradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani (2004) \"Least Angle Regression,\" Annals of Statistics (with discussion), 407-499.\\n(https://web.stanford.edu/~hastie/Papers/LARS/LeastAngle_2002.pdf)',\n",
       " 'feature_names': ['age',\n",
       "  'sex',\n",
       "  'bmi',\n",
       "  'bp',\n",
       "  's1',\n",
       "  's2',\n",
       "  's3',\n",
       "  's4',\n",
       "  's5',\n",
       "  's6'],\n",
       " 'data_filename': '/cm/shared/apps/spack/0.17.3/cpu/b/opt/spack/linux-rocky8-zen/gcc-8.5.0/anaconda3-2021.05-q4munrgvh7qp4o7r3nzcdkbuph4z7375/lib/python3.8/site-packages/sklearn/datasets/data/diabetes_data.csv.gz',\n",
       " 'target_filename': '/cm/shared/apps/spack/0.17.3/cpu/b/opt/spack/linux-rocky8-zen/gcc-8.5.0/anaconda3-2021.05-q4munrgvh7qp4o7r3nzcdkbuph4z7375/lib/python3.8/site-packages/sklearn/datasets/data/diabetes_target.csv.gz'}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diabetes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make this a little more readable, we can print just the description"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploring the Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can examine the dataset to understand its structure and contents:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".. _diabetes_dataset:\n",
      "\n",
      "Diabetes dataset\n",
      "----------------\n",
      "\n",
      "Ten baseline variables, age, sex, body mass index, average blood\n",
      "pressure, and six blood serum measurements were obtained for each of n =\n",
      "442 diabetes patients, as well as the response of interest, a\n",
      "quantitative measure of disease progression one year after baseline.\n",
      "\n",
      "**Data Set Characteristics:**\n",
      "\n",
      "  :Number of Instances: 442\n",
      "\n",
      "  :Number of Attributes: First 10 columns are numeric predictive values\n",
      "\n",
      "  :Target: Column 11 is a quantitative measure of disease progression one year after baseline\n",
      "\n",
      "  :Attribute Information:\n",
      "      - age     age in years\n",
      "      - sex\n",
      "      - bmi     body mass index\n",
      "      - bp      average blood pressure\n",
      "      - s1      tc, T-Cells (a type of white blood cells)\n",
      "      - s2      ldl, low-density lipoproteins\n",
      "      - s3      hdl, high-density lipoproteins\n",
      "      - s4      tch, thyroid stimulating hormone\n",
      "      - s5      ltg, lamotrigine\n",
      "      - s6      glu, blood sugar level\n",
      "\n",
      "Note: Each of these 10 feature variables have been mean centered and scaled by the standard deviation times `n_samples` (i.e. the sum of squares of each column totals 1).\n",
      "\n",
      "Source URL:\n",
      "https://www4.stat.ncsu.edu/~boos/var.select/diabetes.html\n",
      "\n",
      "For more information see:\n",
      "Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani (2004) \"Least Angle Regression,\" Annals of Statistics (with discussion), 407-499.\n",
      "(https://web.stanford.edu/~hastie/Papers/LARS/LeastAngle_2002.pdf)\n"
     ]
    }
   ],
   "source": [
    "print(diabetes.DESCR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear regression fitting function (linear.fit) expects a 2D arrays for the data (# samples x # features). We'll start off by working with a single feature at a time. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing Data for Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to prepare the data for regression analysis by defining the independent (features) and dependent (target) variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract column corresponding to BMI from data set and convert to (n x 1) arrays\n",
    "bmi              = diabetes.data[:, np.newaxis, 2]\n",
    "disease_progress = diabetes.target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we'll fit the model and use the model to calculate the expected disease progression. We also calculate the $R^2$ coefficient, which is the percentage of the change in the dependent variable that can be attributed to the change in the independent variable.\n",
    "\n",
    "$R^2 = 1 - \\frac{\\Sigma (y - ypred)^2}{\\Sigma (y - ymean)^2}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a Regression Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can create a linear regression model using LinearRegression from **sklearn**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept:  152.1334841628967\n",
      "Coefficient:  949.4352603839491\n",
      "Variance score (R2): 0.34\n"
     ]
    }
   ],
   "source": [
    "# Create and fit the model\n",
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(bmi, disease_progress)\n",
    "\n",
    "# Apply the model (predict the disease progression from BMI using linear model)\n",
    "disease_progress_pred = regr.predict(bmi)\n",
    "\n",
    "print('Intercept: ', regr.intercept_)\n",
    "print('Coefficient: ', regr.coef_[0])\n",
    "print('Variance score (R2): %.2f' % r2_score(disease_progress, disease_progress_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing Results\n",
    "Finally, you can visualize the results using **matplotlib**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3aklEQVR4nO2de5RcR33nPzWtGVkjybbUlo0ca2achCxrE8APDOaVgELADlkDBgKRCSchKPSBJYesk7VWCTEnmQVEvMmydsgqwQHSkzgksGsg4eEYkzUJMZZsyZKx5acky2+MHdvI2Eaq/aNvj+703HrcW/fZ/fucU2e6b9+6t+r29O9b9fvVQ2mtEQRBEASAsaoLIAiCINQHEQVBEARhHhEFQRAEYR4RBUEQBGEeEQVBEARhniVVFyCE4447Ts/MzFRdDEEQhEaxffv272mt1yR91mhRmJmZYdu2bVUXQxAEoVEopfaZPhP3kSAIgjCPiIIgCIIwj4iCIAiCME9hoqCUOkop9R2l1E6l1M1KqQ9Hxy9WSt2rlNoRpXNjeTYppe5QSu1RSr2uqLIJgiAIyRQZaH4aeI3W+kml1DjwLaXUV6LP/lhr/Ufxk5VSpwBvB04FTgT+SSn1U1rrQwWWURAEQYhRWE9B93gyejseJdvqe+cBV2itn9Za3w3cAZxVVPkEQRDqwtzcHDMzM4yNjTEzM8Pc3FxlZSk0pqCUaimldgAPAVdpra+LPnq/UuompdTlSqlV0bEfA+6JZT8QHRu85kal1Dal1LaHH364yOILgiAUztzcHBs3bmTfvn1ordm3bx8bN26sTBgKFQWt9SGt9YuAk4CzlFLPBz4J/ATwIuB+4JLodJV0iYRrbtVan6m1PnPNmsS5F4IgCI1h8+bNHDx4cMGxgwcPsnnz5krKU8roI631Y8A3gddrrR+MxOIw8OcccREdANbFsp0E3FdG+QRBGE7q5JYxsX///lTHi6bI0UdrlFLHRq+XAT8H3KqUWhs77U3A7uj1F4G3K6WWKqVOBp4LfKeo8gmCMNzUzS1jYmpqKtXxoimyp7AWuEYpdRNwPb2YwpeBLUqpXdHxVwMfBNBa3wx8Dvgu8FXgfTLySBCErNTNLWNidnaWycnJBccmJyeZnZ2tpkBa68amM844QwuCMLx0u109PT2tlVJ6enpad7td77xKKU0vLrkgKaUKLHE2QuqZBWCbNthVpRu8R/OZZ56pZUE8QRhO+u6feGt/cnKSrVu3smHDBmf+mZkZ9u1bvO7b9PQ0e/fuzbOojUMptV1rfWbSZ7LMhSAItSTU/VM7t0xDEFEQhBGl7iNzQkflbNiwga1btzI9PY1Siunpae9exigj7iNBGEFCXTNlIO6f4hD3kSAIC2jCyBxx/1SDiIIgjCB1mzCVhLh/qkHcR4IwgohrZrQR95EgCAsQ14xgQkRBEEYQcc0IJkQUBGFE2bBhA3v37uXw4cPs3bu3EkGo+7DYUaTIndcEQRCMDA6L7S9YB0iPpUKkpyAIFSAt5GYMix1FpKcgCCUjLeQeTRgWO4pIT0EQSkZayD3qto+A0ENEQRBKRlrIPWRYbDZuuAGmp+GNbyzm+iIKQqOosy/et2xpWsh1rG9eZZJhsem44QZQCs44A/bvhyuvhI98pIAbmTZaaEKSTXZGi263qycnJxdsmDI5OVn4hiR5l8333DrWt45lGna2b9caktOnP53tmlg22ancsIckEYXRYnp6OnEnrenp6aqLlrpsPjtt1bG+dSyTL2XvbhbKDTeYxQC0/t3fzX5tEQVhKKjz9opFlK1O9e0b1KTy1OU7sNGkHo5LDN73Pq0PHw67h00UJKYgNIY6j1Ypomx1qW9/CG3SAnp96vAd2GjCiK8dO3oxg9NPT/78fe+Dw4fh0kt75xWFiILQGMoerZImoFpE2WZnZ5mYmFhwbGJiovTROUkGNU4TRgzVecTXzp09I3/aacmfdzrliME8pi5EaAKOAr4D7ARuBj4cHV8NXAXcHv1dFcuzCbgD2AO8znUPcR+NHmX5hbO4G/IuW7fb1ePj4wvKMD4+XrrLw+TGIool1NEFM0gdYyE7dtjdRL/xG1ofOlTMvakipgAoYEX0ehy4DngpsAW4KDp+EfCx6PUpkYAsBU4G7gRatnuIKAhFUQcjUocy1KkcIdQpprBzZ3Vi0McmCoW5j6J7Pxm9HY+SBs4DPhMd/wzwxuj1ecAVWuuntdZ30+sxnFVU+YThJ2Q8fR3cDXUoAwzHJLM6zInYtavn/nnhC5M/f8974NAh+LM/g7EqHfsmtcgjAS1gB/AkR3oEjw2c82j091LggtjxTwFvsV1fegqCidCWoal13Gq1ShvSWKcWetOGc9aJXbvsPYP3vKf4nsEgVD0kFTgWuAZ4vkUULksQhfMTrrUR2AZsm5qaKuqZCQ0n1KAmicpgKtr9kJfLQwx6NbjE4N3vLl8M+lQuCr0y8PvAhfSCyGujY2uBPdHrTcCm2PlfA862XVN6CvUmqzHKw4jlMcY/Xo5Wq1Vaqz1+33a7rdvtduZnUSdf+qiwe3d9xaBPJaIArAGOjV4vA64F3gB8nIWB5i3R61NZGGi+Cwk0N5asxigvI5a366WsiWR5G/E6uaCGHZcY/OqvVi8GfaoShRcANwI3AbuBD0XH28DV9IakXg2sjuXZTG/U0R7gHNc9RBTqS1ZjlJcRa6pxbaqYjTI339wcMehTC/dREUlEob5kNUZ5GrE8felluWHyNuLtdlt6CgXhEoN3vat+YtBHREEonap7CkVQRsA2z/onTX4D9MTEhMQUAvjud91i8KMfVV1KOyIKQulUHVNoKnnW3yQw7Xa7gJIPP7fcYheDd76z/mLQR0RBqARXy9r0eZVDKOswfDOvMkg8IR+GSQz6iCgIuZKH0aq6R5BUh6rLlIam7sdQFnn8j956q10MNmxonhj0EVEQcqOuQ0bTYKpDFUHZLMaryTu3lUFovYdZDPqIKAi5kZcxr9K1YaqDKRVVpqzGK813UIQ7rA4uNhtZ/0f37LGLwTve0Xwx6COiIORGXsY8VFxM7h8fY2VbCrrMnkLWZ+D7HRQlCHXvfaT9H73tNrsY/NIvaf3ssyVXomBEFITcyKunEGJckvJOTEwsGn5pup5tVE6ZBi+N8Uq75Ibp+XY6nQVusna7nap+VfdS8izjKIpBHxEFITfybClmNRpp3D+DhqDb7SbGDvp1KNOQ+RqvLIvzpXlGaTbu8RWyTqez6NyyehSu/9Hbb7eLwVvfOrxi0EdEQciVqn3Kadw/cWNlMq5pW8tZGXxunU7HKbDdbtfYM7At412Ui8xHyLrdrvH+ZY18SvofFTE4goiCMFRk7Snk6frKa8RQp9MxXsvVQ7DFcYoKpvv0FG33LiPuMcgdd9jF4PzzR0cM+ogoCENF2phC3/CEGkTTvfMeMeTKE+/hpCmna6/lNM/AZsh971N00FrEwIyIgjB0+I4+8vHHpzGIRY8Y8skTF0KbAU1yV01MTCy6TpqYgg+mZ6SU8upRhLqY7rzTLgZvepPWzzwTWMmGI6IgjCyu1vbExESqETlZh+QW0VPIYkAHA+1FxFNMvZROp7PgvLznqtx1l4iBLyIKwsjiam23Wi09Njbm3QrP2rrN4irx6eVUsY6RzX0Ud9X1A+SmWEFePQWXGJx3nojBICIKQi0pI8iYNuDqMk6h8yuyBKhtdSh7xVNb/dM+m9CYwt13ixhkRURBqB1lzYz1aW2nbYVXMSTXtC5TWaLgEqfp6elMLf8sz9IlBr/4iyIGLkQUhNrhMiB5+r5dBi1NTyEv0hrDrLOf8xAtXzdW0bvm7d2rtVJhYlD1HJu6ECQKwBrgvwFbgcv7yZWvjCSi0CziP0ibcel2i9kxrNvtegtCkbuTZekl+bbCi3Bv+Qa884oRLK7DlIZDRjH4hV/Q+umns1y3fus2lUWoKPwr8DHgbcD5/eTKV0YSUWgOvm4cm3HJo/VucsPEg81Fz3DO6mbxMWhFBMJdwfqsMQX381mn4VmjGJx7rp8YhD6bYSRUFHa4zqkqiSg0B5/Wpo8RCh1tU4fWYlY3i4/ro4ghsy6RThp9FOKe6fUMzGJwzjnpxKCP7ER3hFBR+EPgXNd5VaRhEoUifZ0h186rXC5D7+uuyKNVV7VfucgWq+natnWStLYbzLKEdP9+rcfHk4UAtD7qqG/qH/4w+/Wlp3CEUFF4AjgM/DB6/QTwuEe+dcA1wC3AzcBvRscvBu4FdkTp3FieTcAdwB7gda57DIsoFPmjC/UxF72JfNIPsqiYQl0o+/seTFncTkUKqUsM4Ct62bJjCwmYS0whgyhkTcBa4PTo9UrgNuCUSBQuTDj/FGAnsBQ4GbgTaNnuMSyiUEXL0efaeZYryxj2pNFHVbfy86KsnqHP3gv9PGUbzHvu0XrpUlvP4J81LM31+QzL/08owaIA/Cfgj6L0Bp88Cde4EnitRRQ2AZti778GnG275rCIQpG+zpBr512uNLNdTflDDVcRwzXrbGSqHMpqwiUGr32tDnITCW6CRAH4KHA18GtRugr4qCvfwDVmgP3A0ZEo7AVuoje8dVV0zqXABbE8nwLeknCtjcA2YNvU1FTxT68ERqGn0Cfrcg95xBlMrpWso42Srtc3wkUa1TTGu05+9AMHtF62zCwG69dr/dRTpRdrJAkVhZuAsdj7FnCTK1/s/BXAduDN0fsTomuMAbNEcx6AyxJEwTr0dVh6CqMQU+iT1kj5TpwKuberXlnH75c1Q7vIpSTyQMSgfuQhCqtj71f7igIwTs8N9FuGz2eA3dHrkXUfaT38o4/6pHVJuQyvq9XrM2HOdp2Q8ftFtMhNz8O1t0IVLq4DB7RevtwsBq9+tYhBFvIZ9hsmCu8A9gGfBj4D3A283SOfAj4L/MnA8bWx1x8Erohen8rCQPNdjEiguYlk/cdM21PwnThlKmOadY+ShMlWXh/BynsMvO1edYln3HuvXQzGxr6pL7/8r6suZiPJq+cXJAo6MuT0gs3nAc/xzPOKqNA3ERt+CvwVsCs6/sUBkdhMb9TRHuAc1z1EFKqhTJeUK5aQxZ+eplWfdvy+r9hlxTSaqIh7peXee7VeudIsBnCNhqW1KGtTyStGlEkUgOdFf09PSqZ8ZSYRhWoI/cdM08sIEaA0m9ebrumqazwIPni/Inz3ZfZKfLnvPh8xOCqXstZ9tFfR5DUiMKsobI3+XpOQvmHKV2YSUaiGspcLyNtVNVj+pF3B4vdOEqVOp7OoTGUYrDxGYeXFffdpfcwxZjF45Su1Xrfup3Irax2C5lVTaU+hCUlEoRrqNMyxT9IciHa7vWhPYpOg+Qarp6d7ex1XZZy63W4p+yzbuP9+rVetsovBwYNHylvFrPhhpRYxBeCtwMro9e8CXwBOc+UrI4koVEPSMhRlGqWk8ph8++Pj47rdbs8bdJf7xaeFn9U45dWTMM32zopvuVxi8PKXa/2DH2S/vosye6h1dlPVYfTRTdHfVwDX0gs2X+fKV0YSUaiGpNZqlesSuYx93Fj7BJ9dLa8sxqmurg+fcj3wgNbttlkMXvayZDHIm7J6CnX9rvIkVBRujP5+BPjl+LGqk4hCNeT547S1enxbRK6ActxY+w5TtdUlS/198lTROrWVyyUGZ59djhj0KctYj4KbKlQUvgz8b3pDRY+lN49gpytfGUlEoRry6sbbfuRpDECankL/vn3j6yMkacqd9ZlV1TpNLtcaDQ8axeAlL9H6yScLLZaRMoRzFPZdCBWFSeDNwHOj92uBn3flKyOJKFRDXi0p23XS3MPW+s9rm8uke6YxTq77VNU6XXjf+opBmUhPwS0KPwEsjV7/LPAB4FhXvjKSiEI15NWqtbXI0rbWkkYf9Y21y0XlO+S0yGdWVeu02+3qZcumNNxvFIMXv3g0xKCPxBTcorADWAL8JD0X0h8D/+jKV0YSUaiOPLrxefUUXOV0/cDLGnJqe2ZVtE4fekjrE05IFoJRFIM4dR59lAehonBD9Pd3gP8cvb7Rla+MJKLQbPKKKdjIMygc733kRdkzorXuicHatWYxOOMMrZ94IvfbCjUiVBSuo7co3m7g5OjYble+MpKIQvWEtKhc4+1DW2vdbtdo3G1umayL8KUtb5Lw5bkfw2KXWlu3Wg+IGAjBonAK8AngHdH7k4GLXPnKSCIK1VLmwnh5lC2PnoIpf5b6FOkyWlietoZ7jGJw2mkiBqNGkCj08rMM+A8+55aZRBSqxWbUXK3mol00tuv77PSWdmOfLAY+j+Cy6Tn3yuMWg8cf976VMESE9hR+kd5S1ndH718EfNGVr4wkolAtNjeLq9Ucsk9CaNl8ezIu4fK5n83Ahw6Hjbub4s/tk5/8nIZ9RjGAG3Vv5RphVAkVhe3AMcSCy8AuV74yUlNEoa4jGULLZTJqpjX/48Yui4smj7IVMXop6/2yuJzsvZjVHmKwIvjZCs0nVBSui/7eGDvmvUdzkakJolDXMc95lKvT6TgNu6nVnMVFE1o/yLZ4nI94Zn2e+UyCW6XhbosY7NSwslb/f0K1hIrCp4BfprdT2nOB/wX8mStfGakJolDF+POyyhXSU9A6vYvGRpJxHRzdVLRRLH8JhhMtQqD1unXf1+vWnbrgO6lTT1WojlBRmARmgeuj9IfAUa58ZaQmiEJd11FxbTPpY9xCYgpxQnsttvx5i3LVrsBefdZaxeCkkx7V//7vpRZLaBiZRQFoAf9kO6fK1ARRqENPIcmQmcrVbreDF6IbvIdPKzXE2NrKkacoV+0KvO8+sxCA1uPjt+mtW/+2lLIIzSa0p/BF4BjXeVWkJohC1YbEtrZP0vEkd4tJxFx1K6vuNsOfpyhXJfAuMQCtH3us0CIIQ0aoKHwO2E8vtvCJfnLlKyM1QRS0rtbl4NOaj5cr60J0Va7n46pjXsJUtivw/vvdYvDAA4XcWhhyQkXhXUnJla+M1BRRqJK0hixPQ550naR757GchavHEnr9PALivuUQMRCKJkgUsiZgHXANcAtwM/Cb0fHVwFXA7dHfVbE8m4A76E2We53rHqMuCj5GxjZCyHdo5cTExIJ9jn0nf5kEKW5E07TkbfU1fVaE4GTpdfjU84EH3GJw//2pil8aVQfghXSE9hR20RuOGk/X0ltCu23JtxY4PXq9EriN3jpKW4jWTgIuAj4WvT4F2ElvZ7eT6S3T3bKVbZRFwdeY2oyaj/FdsWJFJkNoa1l3Oh3neYOt77wmeqV1HdnqAeZ5D4NG0har8RGD++7zLnLpVB03E9ITKgpb6O3P/NNRmgX+O/BfgS+58seucyXw2qgXsFYfEY49+kgvYVPs/K8BZ9uuOcqikMbN0+12vecODObzae0n4VrGot+atJ2Ttb4+efIYdmsygD4T83ppTaPFoE9VAXghO6Gi8C+mY3gudwHM0AtWHw08NvDZo9HfS4ELYsc/Bbwl4VobgW3AtqmpqQIfW71JGyvIc20eVz5X3rgxHRsbS/ys796yDSt1lSPLPArfFr7NALrr7haDe++1f/91oq5zcQQzoaKwE3hJ7P1ZwM7o9Y0e+VfQWz/pzdF7kyhcliAK59uuLT0F/9ZZltaczai6WoG2XoZv8mltZ+kpmHpNSXM0xsfH9cTEhLMccQNorvdxQyUGfaSn0DxsojCGm18H/kIpdbdSai/wF8CvK6WW03MrGVFKjQOfB+a01l+IDj+olFobfb4WeCg6foBecLrPScB9HuUbSWZnZ5mcnFxwbHJyktnZ2VzOB5iamko8rpSy5gPYsGED733ve1FKWc8z0Wq1OHjwoPUcV/lNdT506FDi+Y888siiez777LOsXLmS6elplFK0Wq3EvPFntfi5HUfPTj5sLOuBAz1ZOPFEmJubY2ZmhrGxMWZmZpibmzPmqwNZ/reEGmNSi8FEb6XUY1Ocr4DPAn8ycPzjLAw0b4len8rCQPNdSKDZStoRH/Hz2+22c0RRkm9cKbUgUOyi0+kYW+YYWuc+PYRWq+Usx+DaR/2gsI9ra7DOtmdijim0nT2De+5xP/MmBG1l9FGzINB9dAzwP4j8+MAleMxwBl4R/VPfBOyI0rlAG7ia3pDUq4HVsTyb6Y062gOc47rHqIlCXj+8vIaBZrlP0j3TLMPhayxtdTR95jub2/VMHn7YLgRJYtBHXDFCGYSKwueBDwM/HqXfB77gyldGGiVRyLMFWZbhsRl2l8D4juDJGkMxraoa8oy/9z23GBx11HOt15OgrVAGoaKww+dYFWmURCHtEFRbS9bHTZIHWQycycWVtsxZjWuWnpGPGMA6r+uZeivtdjuojIIQJ1QUvg28Ivb+5cC3XfnKSKMkCr5GzmfJh6xzD9KStkdiK3sZo63S8sgjbjHYty/dNV2i0NSYg1AvQkXhhfQCwHujdCPwAle+MtIwioKpFehr5FznmT7v76OQd11MK7SmraMp6B3/3OfeLpeVTwvcRwz27s32zFziLzEHIQ8yiwK9/RQ+Hr0+Gjjadn7ZadhEwWZEfXcQcxkVmyumCOKjj1qtll6/fr3RWLvKHu8xJG1YnyQMvm4WHxH5/vfdYgBTut1uZ3bvuIx+GTEHcU8NP6E9hW+4zqkqDZso2FrxSe6ENAvguXoKtgBsVvyXe+jdP6/eUBZs1/QTg+T8Pj0U1zOL5y+6pyDuqdEgVBQuobfRzjuBN/eTK18ZadhEIc0MYJMR73Q6i2bgTkxMLIgppB2qmdUgmAyYSfh87x/aWo73OGxzKOCYIDHIarRNwjw476IIoy3uqdEgVBT+MiFd7spXRho2UcjDiI6Pjy8ydOPj44smV4XELXzJQ+Sy9IZs+PVe3GJw113pekKh7h3TvUw9xqzIkNjRIEgU6pyGTRRswdQkA5hGRHwMZt4Gwbd8aYPcIQbSXqajnWJw552LyxIXMdsEONMM65BnmXcLXnoKo0GQKNCbsPYlegu3PERvCeyTXfnKSMMmClonu4PSBmazGPZuN9vy2q5r+raks1w7jSvFPkvaLQYnnvjKoIB1klsPFvfiBstr22e6iBa8xBRGg1BR+Dd68YQlUboAuM6Vr4w0jKKQRFp3T1rDbjPePgbB5vIxGe88RMe3VWuun1sMLrnk/+YytNX2XfmUt6y5JabyC8NFqCgsEgDg31z5ykhNE4W8f2ymmMJgi9RlxEwGq7/onK3Mvi3LIoKkvu6uxfVb6RSD22/Pt/eUZla26fvwGYYrCD6EisJH6a1mOgNMA78D/B69vZZXu/IXmZokCkV1y01r+PiIj6sF2y+jrcxpWuv9c/uGNlQYfe99xJi6xeC2246U1+b6SppJbnvmaXoKrn0spAUvhBIqCndb0l2u/EWmJolC3QJ4Pv5+n1ayLb/tXnkJos911607xVsM+rjEMv4MfGZut9ttvWTJkkXXSYop1O1/RRg+gkShzqlJolC3oX5pewimMpuEo9VqOe8VMvO3j62F/vjjdiEArffsSb6uaytPn97S4DXGx8f1ihUrFtQ/aahwUl5xFQl5IqJQA+rW+nO5KGyupfiKnT49Bd9RUnkZvieecIvBli1XWq9hi7P4TqZL833bem7iKhLyRkShBrhcHWWP+PARqW63q8fHxxedE58h7XMdV68kL5H0EYNbb/W7VhqXV5r6mXqGtt6GCIKQNyIKNcG2fEHZY8OT7jkxMbFoi07XbmSuZTNM7pC0RtOGjxjccku2Z+QbsA8dQurquQlCngSJAr29li8APhS9nwLOcuUrI9VRFLK0+KtyLcXL2m63F/UKfOMKphFQvhPXstb5ySfdYvDd7xbw4BJIM+kwCVtvQ5aYEPImVBQ+CVwG3BK9XwVc78pXRqqbKGRt8duMZFnYfOg2A57FyJmSb+/IRwxuvjnvJ5SeNA0E2wx16SkIeRMqCjdEf2+MHdvpyldGqpsoZG3x+4zgKRqXsU5ryNMEX10C08dHDHbvLuoJ+ZM1PtTpdGTUkVAKoaJwHb3NdvrisCYuEFWmuolC1klHVfcUXK3UtMs0aJ2up2Ba/6fPD37QDDHQOjw+JEtMCGUQKgob6O2ncACYBfYAb3XlKyPVTRR8x6tnnRWchjTGxWbA44Fn0zmDO6P1YxRJy20sX77cu64+YrBrV/r6FklZ8aG61FdoJkGi0MvP84D3Ae8H/qNnnsvpraq6O3bsYuBeYEeUzo19tgm4IxKd1/nco26ikHUUSt6jj9JeL62rJ0k4kkYpjY+PLxrN5DOJ7+BBfzHo19e2sVDe2AxyWdtlykqmQgihPYWfAJZGr38W+ABwrEe+VwGnJ4jChQnnngLsBJYCJwN3Ai3XPeomClpnX/I4z5Zf2tZqlqCwbxq8p61sPmKwc+fi8puGzcYn2eWFyyCX0VOoarSaMDyEisIOektm/yS9lvwfA//oyhflnfEUhU3Aptj7rwFnu65fR1EYpIofcNrYRtbhoz4pSfwG77Vs2epMYtDHdv+8cX2fSRP+XDGTtNRtyRShedhEYQw3h7XWP6K3N/P/1Fp/EFjrkc/E+5VSNymlLldKrYqO/RhwT+ycA9GxRSilNiqltimltj388MMBxSiH2dlZJicnFxybnJxkdna2sHtOTU0lHldKsW/fPrTW7Nu3j40bNzI3N8eGDRvYunUr09PTKKWYnp6m3W4XUpb4veAoQPPUU48Y8+/Y0ZOFF7wgl+KkYm5ujpmZGcbGxpiZmWFubo79+/cnnhs/rpRa8Nng+1BM36/puCCkwqQW/URv9NE7gN1EO64Ra/078s6wsKdwAr2RTGP0gtaXR8cvAy6Infcp4HzX9ZvQU9A6m2soxJ2UxwzbPHoPJj/3U0+53UQ33uj/TGyb+HQ6He/n5qr/5OSkc4Z32p5h1v8NiSkIIRDoPjoF+ATwjuj9ycBFrnw6QRRMn9Fw91HeI0Hy+NFnjW2YrmEaNWRKSfsP+4gBnDb/DH1mXPfPS1qjadAop3l+pmfWbret300a153PktumcqedGCcjlYQ4QaIQkgZFAVgbe/1B4Iro9aksDDTfRUMCzUW02oqIQ9iEwcdQuILRfWOYdC0fMTjqqLMXXC9pBznbM+l2u85z03wvNr+9zcianlPSsGRTryPPCWxl9SpEeJpFkCgAzwX+HvhuZKzvwmNzHeBvgPuBZ+nFCN4N/BWwC7iJ3tyHuEhspjfqaA9wjuv6uiaiUIQB9w0kpm0tZpmZ7CqTTVR++EO3GGzfHj76qY/v+SEi6Ppe07ju0qSs/09lDHQQd1bzCBWFbwHrI0M+TW8E0Ydd+cpIVYhCHm4ZFz4/5Cw/xCwzk9OUqY+PGGzbduT8EGMZX1ratk5TWhEMMXTx55ymTK56ZqGMkUplCI+QL6GisD36uyt27FpXvjJS2aKQRwDXB581cLL8EF2i4IovuIykjxhcf/3ia4cazn6dO51Opnyu55VXsH8wJcUnTP9PWXepK8NgyxDZ5hEqCv9Cb7TQF+jNaH4TsMeVr4xUtij4ujny9gErpRaNokn7Q/QxUlmN5NNPu8XgO98xXy9EEAbr3Ol0vEWmSKPlqtfgnhO2JbeT4itpei1Fu3akp9A8QkXhxcAK4CTgL4HPAy915SsjlS0KNt/w4HIOWfH9gaX9IfoaqTRkFQOt/VvScYPpGg6adI+s7rJQssRgBsvtqrfPtZKuV0SQWWIKzSJIFOZPhOW+55aV6tRTyMvApAky57W+UVpD8elP/7VTDOCsxPL49g5MeZOEJGn4a8izyoM8W8+uQLVPXcoQhiKvL+RLaE/hbHojj/ZH718I/KkrXxmpipiC6YeZlyvCZUwGx+779lDSrE5qwqdnAGcZr+87IS6pLnExSTKSLsMYdyu1Wq3Mk9p8scWf0hpNHxF1xZKKFkURhWYRKgrXAetYuMmO14zmolMVo4/SujDSYvsBZ/1xr1+/PrHMY2NjXj/eZ57xEYOXOMUyq3HzFZNWq2Wc6FWFe8MmZGn3WHDV39YoKdrnL+6j5hEsCtHfG2PHRnbntSpbXVl+3K5JXbZW3TPPaN1qucTgpd5G3uUGSQqo2+qdlJK+i6oDoXncPyQ+UvTooKqfr5CeUFH4e+BlwA3ABHAh0UzkqlNVk9eq6ipn+XHbgpQmQ+onBgtnICcNr+wfd4may6hn2doz9LnlSZ73z9IoKdpoV/18hfSEisJxwBzwIL1Nc7pA25WvjFSHGc1lkvbH7bP0Q/wazzyj9cSEXQw+9KGvWt1bSSJkc3/51MdH2GzGyPe5mcQ+tBGQt1FOW56ie7fSU2geQaJQ5zRqopD2x+3vdmlp+KFVDK69dmE5TEYpTaDcx6h3u4t3VYNePGRsbMzb2Luem2nCYNK8gSyLE5bpc0/6fors3UpMoXmE9hS2AEcD48DVwPeILXNdZRo1UdA6XSvR7XZpafiBtxj4kMaV4NPCNJ3Td0v5GiPbc+t2zduEmibCpW0Fl+VyrDqoLqOPmkGoKOyI/r4J+AywmhEONDcJc0+hpeEJqxiccMJbMrUw07gSfAyYS2Rcxt6n7GkC2TaRqwPiyhF8CBWFm6O/fw68Pno9sqKQR4uoulZjS8PjVjGAV86fPzExYdzDwP+e9jxJ8wfizydrSz1NOWw9qrx6CmUhQV/Bh1BR+ChwK3AjPRfSGqJhqlWnOiyIV7V/2SUw3W5XT02drOFRhxi8yruV7GOQfUQv6Vn47KPg87zStJhN5/aHyDbJXy49BcGHIFHo5WcV0aY3wCTwHJ98Rae6LHOR5geX54/WJTA/+pHWxxxjEwKt4We8xSDvVmcat02r1UrVswrdAS0+Z6JJ/nIJ+go+ZBIF4DXR3zcnJVO+MlNdFsRLYyTz6t53u12ja2Nq6mS9erVdDFqtn1uwREaaYZ99wxhqKNPOP/B5Ji63k2lmcZMMv4thqotQDFlF4cPR379MSJeb8pWZhrWn4OMSSh7vP6bhe46ewauN9xy8ZlJMQSml169fn0trNG1PYfAZxJ9RkpvHJQhZvr+6kFR/EQLBl0yi0IQ0jDEFn3ssNqZjGh50iMFrEo1knKSgb9L4fVsPIvR52pJtLSHbkNIid8mrAp/nJi4jwUbWnsJv2ZIpX5lpGEcf+fQkjhhAtxhcdVXY9p5p3EpZl21wuXzAvIxGmjINSxDWt4fVtHoJ5ZFVFH4/Sn8N3A5cEqXbgL8w5SszhYhCXf2uPjGHqakZDfdZxeDrXz9yTVvvI976Dk3Lly8PWp7atYlRljL5CF9dvntffHtuTesBCeWRSRTmT4CvAytj71cCX3XlKyNlFYU6Gwdba/bQIa2f8xyzEIDWxx9/QaJ/eXBdona77eWHD01phMFU935A3McIur7TshsDRdxPegpCKKGicCuwNPZ+KXCrK18ZKaso2Axv2evGDJIkWMuWLderVtmXo4DXGo2DaQ0fH0Ob1W3TT4MB4rR17xt2lyHs17FOvb+iGh8SUxBCCRWFzcBO4GJ67qQdwCaPfJfTW1V1d+zYauAqeu6oq4BVsc82AXcAe4DXua6vA0TBZgzz3DTdhUlsjhjBMd1q2d1Exx//K97GOa1BV0otEsV2u71oRJIr+dTZ55nktZNZWZiELO2ciyRk9JEQQpAo9PJzOvCbUTrNM8+ronxxUdgCXBS9vgj4WPT6lEh4lgInA3cSTZazpbx7CmkMZ2jX3NaKPHRI66kpW69A6698pXedtGP90yZXz8n1zOI9hdCWc13jQCZ8vhtp0QtVECwKWRMwMyAKe4C10eu1wB59pJewKXbe14CzXdfPO6aQxlimCeLF3R/xQOzi6yrdat1jFYN/+IfF1/RJpmWmXcHd0AD18uXL58+vYh2hKoVEfP9CXamTKDw28Pmj0d9LiS3HDXwKeIvhmhuBbcC2qampTA8kKeia1sj6/pD9xuIrDXd6iYH/NReLwqAbbHx8XK9YsSLxfNtw1NA4Q1IqgqoHFPh+T3mMEmpaL0qoliaIwmUJonC+6/pZegqu4Zl5xxTsQuMWgy9/Od01W62WdUG5fmu93W4vOm/QR1+0ayp+3yKMmG1AQVn4uNqKdEUKQhJ1EoXK3UcuQ5H36COzYb3NKgZf+lKWa/ol16S0+CisMkRh8L55UbdlpIuKqdRB/IRmUSdR+DgLA81botensjDQfBcFBZpthqKcMeW3WsXgyiuzXDP/ZBrGaktZRjgl3TcvYTAJX97GMs3/Tdb/MZug1E38hPpTiSgAfwPcDzwLHADeDbTpbel5e/R3dez8zfRGHe0BzvG5R549BVtQNYRut6uXLZvUcItVDFauvGDBqqUuw1L0pDPoxSE6nU7m2cRZU7vdDnrm/WeUNHR2YmIi195IWa4bW29AegpCWioRhTJSnjGFIlqVhw9rfcopZiHopfO8W8yD8wVsxto02sgmgklpyZIlznkJefQQBlOoUbWJf56UZZBdPVyJKQhpEFEYIKkLb/PTp+3q+4jBmjUbnYYxbliyDKM1Gay+IUmzAqrtHmmu4TvSa7DuaV0uZblUbHXP0yhniYUJggkRBQ98/PSu1tfhw1r/9E+7xSCNEXWNXLGNaDHNAu6vRxTqFupfyzfGMfj8ut2u9dr9c7K0gotuwfsuvZGXce50Oon3SLvooCBoLaLgha+fPsmo+IjBBz7wz5kmyPkaW5Mxsi0bESII8Xv5bnCzfv36RfNDli9fbn3OWY17kS6VNDGdvERI4gZCnogoeBI3oK4WrNY9MXjb2+xi8Hd/1zu3qBFD8eGjaVwHeZZnOlp7J0tsISlmETfeIW6golwqaZ5dXu4qGWEk5ImIQgZsLbPDh7X+7d+2i8HnPte7Tpbx/ml6CC5DZzKMviuk+pY5ZDRUu902Gu86tpDTuv/ywPc5SGxB8EFEIQNJLoJlyyb1G96w2yoG8Lb5H6NpWOTgjzrph2wyAoMrbNqMgM2FYhMqH6McT66RTj4imOZ7qHpkjS2AX1Q5fZ5DHZ+VUE9EFDJyxHgqffTRf2oVg4mJdy76MZr85T4/WFOQOC4kJiPgCv728/tsZO/yn9uW1LAZzDSt6bq1fl3Pvahyup5DHXtVQj0RUcjI4cNab9pk6xVo/a1vZfPP+xiNeIs+qRVqcu+43Bv9lrlPy31w8cAVK1YsmGRnczHFezWdTidoMlkdhaFO5dFa4g6CPyIKKfERg2uvPXJ+lrH+fXyMS95B6n7L0dWy9HFHpBmnPygw/TkLLsQt4of0FARfRBQ8OXxY682b7WKwdOl67267KfVn1Xa73UXul6SWc56rlQ4OV7UZW1u9Wq2W1U2V58zhuhu7uvQaRDwFX0QUHPiIAbzCaIzSjFuPG32T62VsbMwr6Jw2jY2NeS2d0b+3zzXXr18fbIhcRjUPt0hRhrsoQxyycF4dBEqoNyIKBg4f1vr3fs9fDGzGyGeOAyx0qfi27NOuVmpLJrIustdqtYIMkY9RdQmU655FtqCL6MVIi18oGhGFBD70IbsYfPOb2X/wtny+4mHL5wrw+opS/JohLqoQfJ6xj2DZjGaR7qcigrt1d5cJzUdEYQBbEPmaa46cl7XFZhuymKU1buqZZDHgtl3msqRWq5XpO+jja1S7XfMezy6jWeSonCIMuIwiEopGRGGA5z3PLgZx0rpG4jGAvhELjQ2YgrZZewvT0+Y1+NOm9evXZ/oO+tjKMfi8fYfa+t4jj5Z3Ea4e6SkIRSOiMMAVV2j9/Of3av+Nb2S6RCIuAxHioklaDTNraz/UXZTFUJmGo7pmffuOhrKVpWgffd7BXYkpCEUjolASJqMVugxE35AnGYWsMQqfnkK73Xa6bHwXpksy/EuWLNGdTsc5K7pv7NevX+8lHqYyNGlUTtPKKzQLEYWCCXEN5bm4mk8ZfGMK/fNM6/j7lslVLh/B7Aui6Vn1R0AJguCHiEIB2JagSGPk0+x+5mqVJxn6iYkJ497P8dbo8uXL5w10f2Ja/5wlS5Zkap33yWN3N5uwSADWjvQ6hEFEFHImj5E7cUPmKww+rfI0BsC1tpKtB5SmdR4S1O6XI2RBvVFG4hNCEiIKOZPHyB3bEtgrVqzI3Cr3IY0IpRkeaRIk1/DZwXjD4GqwtmduirVkJa9WdV1a5zKSSUiidqIA7AV2ATv6hQNWA1cBt0d/V7muU9bS2Vk2qMnSGva5dyiuGMGgwfUdMupqkdrua3JxDS7BkSQeee5RnFeruk6tc5nzICRRV1E4buDYFuCi6PVFwMdc1yl7kx2fDWoGW7u+qYyWW9oJb30D7dpPwfZM4vWyTUAbrH/SiKVWq2WMj+RBXq3qOrXO61QWoT40RRT2AGuj12uBPa7rVLUdZ9aYgi1P1pabb48ibZkHV1PNOoJosF6+LVfbXhFFiUJereo6tc7r1GsR6kMdReFu4AZgO7AxOvbYwDmPuq5TpCi4fthpRg3FxSTPlluaH3yaOIgpiGzL47OLm60cg+elFa48GMaegtb1iW8I9aGOonBi9Pd4YCfwKl9RADYC24BtU1NThTwwrd0/7DRGdtCFMugWGR8fz/RDTWN8fAXMZmhdeX32KHYJWdo5H3ka2mGMKQhCErUThQUFgIuBC6mZ+yjPJSsGt7f02VjHp3Vn680M5retkzQ42ik+TyGOz1pLvq4s0yiltG65vF0ywzb6SBCSqJUoAMuBlbHX/wq8Hvg4CwPNW1zXqmr0kdbZewq+AVmflqbpWu12e1H+8fHxRWLkmrU8KAyudYpCW+1pnmle9xSEUaRuovDj9FxGO4Gbgc3R8TZwNb0hqVcDq13XqnpG86DhTZr567uncbzFa2qRJ43QSRIPW4t+cOVWrbVxRFDSstiDC9vZRMv03ExCa+t9JQ1JFZeMIGSjVqKQZ6p67aMkA+dyG7h6CrZho74TxlyurUFjajs3Td19nleWPaHjz0ZcMoIQjohCjUgateS7PLSvq8THDRO/VpqeQgg+Rl8CtIJQPCIKNSHJ6A3OyrW18n2No0/AdnDtpaRz8pwtbKvbYCBeegOCUCwiCjXBJ8hsCx6nwTW0c7DX0el05nsMptFHoaQZQisIQnGIKNQE35Zyni6UOrlk8pyjIQhCdmyiMIZQGlNTU87jGzZsYOvWrUxPT6OUYnp6mq1bt7Jhw4ZM98z7eqEopazvBUGoFtUTjWZy5pln6m3btlVdDG/m5ubYuHEjBw8enD82OTlZqZEuk5mZGfbt27fo+PT0NHv37i2/QIIwoiiltmutz0z6THoKJVK3VnvZ7N+/P9VxQRDKZ0nVBRg1NmzYMDIiEGdubo6xsTEOHTq06DOTW00QhPKRnoJQOH23WZIgTE5OMjs7W0GpBEFIQkRBKJzNmzcviKP0abVaI+U+E4QmIKIgFI4pZnD48GERBEGoGSIKQuH4DMUVBKEeiCgIhTM7O8vk5OSCYxJLEIR6IqIgFM6oD8UVhCYhk9cEQRBGDJm8JgiCIHghoiAIgiDMI6IgCIIgzCOiIAiCIMwjoiAIgiDM0+jRR0qph4H4WszHAd+rqDhlMyp1lXoOH6NS1zrXc1prvSbpg0aLwiBKqW2mYVbDxqjUVeo5fIxKXZtaT3EfCYIgCPOIKAiCIAjzDJsobK26ACUyKnWVeg4fo1LXRtZzqGIKgiAIQhjD1lMQBEEQAhBREARBEOZpnCgopVYrpa5SSt0e/V1lOO9ypdRDSqndWfJXTYp6vl4ptUcpdYdS6qLY8YuVUvcqpXZE6dzySu/GVO7Y50op9Yno85uUUqf75q0bgXXdq5TaFX2HtV4S2KOez1NKfVsp9bRS6sI0eetEYD3r/31qrRuVgC3ARdHri4CPGc57FXA6sDtL/qqTTzmBFnAn8OPABLATOCX67GLgwqrrYaibsdyxc84FvgIo4KXAdb5565RC6hp9thc4rup65FTP44EXA7Px/80mfach9WzK99m4ngJwHvCZ6PVngDcmnaS1/n/A97PmrwE+5TwLuENrfZfW+hngiihf3fEp93nAZ3WPfwOOVUqt9cxbJ0Lq2iSc9dRaP6S1vh54Nm3eGhFSz0bQRFE4QWt9P0D09/iS85eFTzl/DLgn9v5AdKzP+yN3xOU1c5O5ym07xydvnQipK4AGvq6U2q6U2lhYKcMJ+V6a9J2GlrX23+eSqguQhFLqn4DnJHy0ueyyFEkO9VQJx/pjjD8J/EH0/g+AS4BfS1vGgrCV23WOT946EVJXgJdrre9TSh0PXKWUujXqBdeNkO+lSd9paFlr/33WUhS01j9n+kwp9aBSaq3W+v6oi/1QysuH5s+NHOp5AFgXe38ScF907Qdj1/pz4Mv5lDoXjOX2OGfCI2+dCKkrWuv+34eUUv+HnvuiVkYkwqeeReQtm6CyNuH7bKL76IvAu6LX7wKuLDl/WfiU83rguUqpk5VSE8Dbo3wM+KTfBOxOyF8VxnLH+CLwK9HInJcC/x650Xzy1onMdVVKLVdKrQRQSi0Hfp56fY9xQr6XJn2nmcvamO+z6kh32gS0gauB26O/q6PjJwL/GDvvb4D76QV7DgDvtuWvW0pRz3OB2+iNiNgcO/5XwC7gJnr/tGurrtNA/RaVG3gv8N7otQIuiz7fBZzpqnNdU9a60hvhsjNKN9e9rh71fE70W3wceCx6fXTTvtOs9WzK9ynLXAiCIAjzNNF9JAiCIBSEiIIgCIIwj4iCIAiCMI+IgiAIgjCPiIIgCIIwj4iCIGREKXUoWu1yp1LqBqXUy6LjM0oprZT6g9i5xymlnlVKXRq9v3hwBU1BqAMiCoKQnae01i/SWr8Q2AR8JPbZXcAbYu/fSm9suiDUGhEFQciHo4FHY++fAm5RSp0Zvf8l4HOll0oQUlLLtY8EoSEsU0rtAI4C1gKvGfj8CuDtSqkHgEP01sg5sdQSCkJKRBQEITtPaa1fBKCUOhv4rFLq+bHPv0pvhdoHgb8tv3iCkB5xHwlCDmitvw0cB6yJHXsG2A78F+DzFRVNEFIhPQVByAGl1PPobdX4CDAZ++gS4J+11o8olbQUvyDUCxEFQchOP6YAvZVO36W1PhQ3/lrrm5FRR0KDkFVSBUEQhHkkpiAIgiDMI6IgCIIgzCOiIAiCIMwjoiAIgiDMI6IgCIIgzCOiIAiCIMwjoiAIgiDM8/8BcjkcR60ULo4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(bmi, disease_progress,  color='black')\n",
    "plt.plot(bmi, disease_progress_pred, color='blue', linewidth=3)\n",
    "plt.xlabel('BMI')\n",
    "plt.ylabel('disease progression')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we used all the the available data to train the model. Often in machine learning applications we want to set aside some of the data to test the model. This allows us to determine if our model has predictive value.\n",
    "\n",
    "Let's go back and divide our data into training and testing sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "bmi_train = bmi[:-20] # All but last 20 elements\n",
    "bmi_test  = bmi[-20:] # Last 20 elements\n",
    "\n",
    "disease_progress_train = disease_progress[:-20] # All but last 20 elements\n",
    "disease_progress_test  = disease_progress[-20:] # Last 20 elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept:  152.91886182616167\n",
      "Coefficient:  938.2378612512633\n",
      "Variance score (R2): 0.47\n"
     ]
    }
   ],
   "source": [
    "# Create and fit the model\n",
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(bmi_train, disease_progress_train)\n",
    "\n",
    "# Apply the model\n",
    "disease_progress_pred = regr.predict(bmi_test)\n",
    "\n",
    "print('Intercept: ', regr.intercept_)\n",
    "print('Coefficient: ', regr.coef_[0])\n",
    "print('Variance score (R2): %.2f' % r2_score(disease_progress_test, disease_progress_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjcElEQVR4nO3de5xcdX3/8dc7XAJbQIgEuYTdBUQqtBJgQQvWH16qFLRRKv6wCwWkXfBH5VKqgqsV1H3gHa1VJCqKskKD0B8RaC1SReHHxQQChosSzG4SiBBALiUQSPbz++OcnRmS2Z0zu3Pm+n4+HvOYM+d8z8xnZ3fnM9/L+X4VEZiZmQHMaHQAZmbWPJwUzMyswEnBzMwKnBTMzKzAScHMzAo2b3QA07HjjjtGb29vo8MwM2spixcvfjwiZpc71tJJobe3l0WLFjU6DDOzliJpdKJjbj4yM7MCJwUzMytwUjAzswInBTMzK3BSMDOzAicFM7MWMjw8TG9vLzNmzKC3t5fh4eGaPn9LD0k1M+skw8PDDAwMsHbtWgBGR0cZGBgAoL+/vyav4ZqCmVmLGBwcLCSEcWvXrmVwcLBmr+GkYGbWIlasWFHV/qlwUjAzaxHd3d1V7Z8KJwUzsxYxNDREV1fXy/Z1dXUxNDRUs9dwUjAzaxH9/f3Mnz+fnp4eJNHT08P8+fNr1skMoFZeo7mvry88IZ6ZWXUkLY6IvnLHcqspSNpK0h2S7pZ0r6Tz0/2zJN0g6cH0foeSc86VtEzSbyS9I6/YzMysvDybj9YBb4mI/YG5wBGS3gCcA9wYEXsDN6aPkbQvcCywH3AE8A1Jm+UYn5mZbSS3pBCJ/0kfbpHeApgHXJruvxR4d7o9D7giItZFxHJgGXBIXvGZmdmmcu1olrSZpCXAY8ANEXE78KqIWA2Q3u+UFt8NWFly+qp038bPOSBpkaRFa9asyTN8M7OOk2tSiIgNETEXmAMcIulPJimuck9R5jnnR0RfRPTNnl12NTkzM5uiugxJjYingJ+T9BU8KmkXgPT+sbTYKmD3ktPmAI/UIz4zM0vkOfpotqTt0+2tgbcBDwALgRPSYicA16TbC4FjJc2UtAewN3BHXvGZmdmm8pwldRfg0nQE0QxgQURcK+lWYIGkk4EVwDEAEXGvpAXAfcB64LSI2JBjfGZmthFfvGZm1mEacvGamZm1HicFMzMrcFIwM7MCJwUzMytwUjAzswInBTMzK3BSMDOzAicFMzMrcFIwM7MCJwUzMytwUjAzswInBTMzK3BSMDOzAicFMzMrcFIwM7MCJwUzMytwUjAzswInBTMzK3BSMDOzAicFMzMrcFIwM7MCJwUzMytwUjAzswInBTMzK3BSMDOzAicFMzMrcFIwM7MCJwUzMyvILSlI2l3SzyTdL+leSWek+8+T9LCkJentyJJzzpW0TNJvJL0jr9jMzKy8zXN87vXA2RFxp6RtgcWSbkiPXRgRXywtLGlf4FhgP2BX4KeSXhMRG3KM0czMSuRWU4iI1RFxZ7r9LHA/sNskp8wDroiIdRGxHFgGHJJXfGZmtqm69ClI6gUOAG5Pd/2DpHskXSJph3TfbsDKktNWMXkSMTOzGss9KUjaBrgKODMingEuAvYC5gKrgS+NFy1zepR5vgFJiyQtWrNmTT5Bm5l1qFyTgqQtSBLCcERcDRARj0bEhogYA75FsYloFbB7yelzgEc2fs6ImB8RfRHRN3v27DzDNzPrOHmOPhLwHeD+iPhyyf5dSoq9B1iabi8EjpU0U9IewN7AHXnFZ2Zmm8pz9NFhwPHAryUtSfd9DHi/pLkkTUMjwCkAEXGvpAXAfSQjl07zyCMzs/rKLSlExM2U7ye4fpJzhoChvGIyM7PJ+YpmMzMrqFhTkDQb+Hugt7R8RHwgv7DMzKwRsjQfXQP8Evgp4DZ+M7M2liUpdEXER3OPxMzMKlq6FK64At76Vnjzm2v//Fn6FK4tnbTOzMyyGx4epre3lxkzZtDb28vw8HDVz7F2LZx0Ekjwp38KQ0NJUvjVr2ofb5aawhnAxyS9CLyU7ouI2K724ZiZtY/h4WEGBgZYu3YtAKOjowwMDADQ399f8fzLL4e/+ZvyxyKSW61VrClExLYRMSMitkq3t3VCMDOrbHBwsJAQxq1du5bBwcEJz3noIdhvv6RWMFFCADj/fDgkhylDM12nIOmvgDelD38eEdfWPhQzs/ayYsWKTPtffBE+8hH46lcrP+ePfwzvfGctoiuvYk1B0mdJmpDuS29npPvMzGwS3d3dk+6/9tqkRjBz5uQJ4cwzYd26pLkoz4QA2WoKRwJz0wnskHQpcBdwTp6BmZm1uqGhoZf1KQBstdVebLHFrajcfA8l9t0XrrkGXv3qnIPcSNYrmrcv2X5FDnGYmbWd/v5+5s+fT3f3nsClQPDCC8tYtmziGZ6Hh5Mawb331j8hQLaawgXAXZJ+RjKX0ZuAc3ONysysDVx6KZx4Yj8w+UijE0+Er38durrqEtakKiaFiLhc0s+Bg0mSwkcj4vd5B2Zm1opWroQJuhJeZs4cuP765LqDZjJh85GkP07vDwR2IVkEZyWwa7rPzMyAsTF473uTTuNKCeGb30zKr1zZfAkBJu9T+Mf0/ktlbl/MOS4zawG1uFq3lS1cmCSCzTaDq66avOzvf5/0FZxyChU7mRtpwuajiBhI73OYXcPMWt10r9ZtVY8/DjvtlO1q4quugqOPzj+mWspyncIxkrZNtz8u6WpJB+Qfmpk1s6lcrduqIuCDH0y+4c+ePXlCmDcPNmxIyrRaQoBso48+ERFXSnoj8A6SpqNvAq/PNTIza2pZr9ZtZTfdBIcfnq3s8uXQ25tnNPWR5TqF8TUUjgIuiohrgC3zC8nMWkGlq3Vb1bPPws47J7WCSgnh298uTkzXDgkBsiWFhyVdDLwPuF7SzIznmVkbGxoaomujgfVdXV0MDbXmMuvnnZckgu22g0cfnbjcoYcmcxVFwMkn1y28usnSfPQ+4AjgixHxlKRdgA/nG5aZNbvxzuTBwUFWrFhBd3c3Q0NDLdXJfOedcNBB2couXZrMXtruFBW60CXtBayKiHWSDgdeB3w/Ip7KPboK+vr6YtGiRY0Ow8xayAsvwMEHJx/ylXzuc8nspe1G0uKI6Ct3LEsz0FXABkmvBr4D7AH8sIbxmZnl7mtfS5qHtt568oSw997w3HNJ81A7JoRKsjQfjUXEeklHA1+JiK9JuivvwMzMpuu3v4V99slW9rbb4PUeU5mppvCSpPcDfwuML66zRX4hmZlN3fr1yfrFUuWE8OEPF0cPOSEkstQUTgJOBYYiYrmkPYDL8g3LzKw6w8Nw3HGVy227LaxYAdtvn3tILSnLLKn3Sfoo0J0+Xg545TUza7hVq2D33bOV/clP4O1vzzeedpBlmot3AUuA/0wfz5W0MOe4zMzKGhtLmoakygnhpJOS8hFOCFll6VM4DzgEeAogIpaQjEAyM6ubT3+6OCNpJatXJ4ngkkuae0bSZpQlKayPiKc32ldxfkBJu0v6maT7Jd0r6Yx0/yxJN0h6ML3foeSccyUtk/QbSe+o7kcxs3YzMlKsFfzzP09edsGCYqfxzjvXJby2lKWjeamkvwE2k7Q3cDrw/zKctx44OyLuTGdZXSzpBuBE4MaI+Kykc4BzgI9K2hc4FtgP2BX4qaTXRMSGCZ7fzNrUdtslcxBlsX59ttqDZZOlpvAhkg/qdSQXrT0NnFnppIhYHRF3ptvPAvcDuwHzSFawJr1/d7o9D7giItalndnLSJqtzKwDfOtbxVpBpYRw663FWoETQm1NWlOQtBmwMCLeBkx5knRJvcABwO3AqyJiNSSJQ9JOabHdgNtKTluV7tv4uQaAAWj92RjNOt3jjyfrE2Tx3vfClVfmG49VqCmkTTdrJb1iqi8gaRuSqTLOjIhnJitaLoQyMc2PiL6I6Jud9a/JzJrK3LnFBWsqGZ9ywgmhPrL0KbwA/DrtD3hufGdEnF7pRElbkCSE4Yi4Ot39qKRd0lrCLsBj6f5VQOkAsznAIxniM7MWsHBhsipZFtdeC0cdlW88Vl6WpHBdequKJJFMoHd/RHy55NBC4ASSC+BOAK4p2f9DSV8m6WjeG7ij2tc1s+bx3HOwzTbZyh58MNzh//iGy3JF86WVykzgMOB4klrGknTfx0iSwQJJJwMrgGPS17lX0gLgPpKRS6d55JFZa5o3L6kZZPHEEzBrVr7xWHYVk4KkX7Np2/7TwCLgMxHxRLnzIuJmyvcTALx1gnOGgNZctsmsw91yC7zxjdnKfve7cOKJuYZjU5Sl+eg/SNZpHl9D4ViSD/unge8B78olMjNreuvWwVZbZSs7a1ZSK7DmliUpHBYRh5U8/rWkWyLiMEkZ5iQ0s3bT1weLF2cru2JF9knrrPGyXLy2jaTCTOOSDgHGu47W5xKVmTWdm28uXlxWKSFccEHx4jInhNaSpabwd8Al6fUGAp4BTpb0R8AFeQZnZo01NlbdFcPjM5ha66pYU4iIX0XEnwJzgbkR8bp033MRsSD3CNvQ8PAwvb29zJgxg97eXoaHhxsdktnLvP/92Wck/eUvi7UCJ4TWl2X00SuATwJvSh/fBHyqzMyplsHw8DADAwOsXbsWgNHRUQYGBgDo7+9vZGjW4R54AF772mxlDzoIFi3KNx5rDEVMPgu2pKuApRQnsTse2D8ijs45tor6+vpiUYv9Zfb29jI6OrrJ/p6eHkZGRuofkHW8ar7dv/ACzJyZXyxWH5IWR0RfuWNZOpr3iohPRsTv0tv5wJ61DbFzrFixoqr9Znn42MeKncaV/OhHxeYhJ4T2lyUpPC+pcEmKpMOA5/MLqb1NNLOrZ3y1LKbTH7V6dTERXFBhiEhXVzER/PVfTzNoaylZksKpwNcljUgaAf4VOCXXqNrY0NAQXV1dL9vX1dXF0JAv5LbJjfdHjY6OEhGF/qhKiWE8Eey6a+XXeOqpJBE891zFotamJk0K6XoKx0XE/sDrgNdFxAERcU9domtD/f39zJ8/n56eHiTR09PD/Pnz3clsFQ0ODhYGKIxbu3Ytg4ObLnVy0UXZm4e++tVireAVU54k39pFlo7m/46It9Qpnqq0Ykez2VTNmDGDcv+vkhgbG+PZZ5NlLLOq8K9vbWy6Hc13SVoo6XhJR4/fahyjmVUwUb+T9CRStoSwcmWxVtDKfK1PfrIkhVnAE8BbSCa/exfwzjyDMrNNvbw/6l0kkxcHY2PbT3re2WcXE8GcOTkHWQdT7VuxbCo2HzUzNx9ZJ6lmRlJo/drARHytz/RNq/lI0p6SfixpjaTHJF0jaY/ah2lm5Wy2WdJhnCUhLF3aHs1Dk/G1PvnK0nz0Q2ABsAvJMplXAlfkGZRZp/vJT4qjh8bGJi/7nvcUE8F++9UnvkbytT75ypIUFBE/iIj16e0yNl2JzcymaXxCOQmOOKJy+Q0bknOuvjr/2JqJr/XJV5ak8DNJ50jqldQj6SPAdZJmSfLKqmbTdOCBSSKYkeG/8aqrirWCLOXbka/1yVeW6xSWT3I4IqJh8yC5o9la1ZIlcMAB2cu3cx+B1d9kHc0Vp86OCHcqm9VINTOSPvdcMgeRWT11aAXUrH5e//rsU05ceGGxecgJwRohy3KcZlalhx6CV786e3k3D1mzcFIwq6FqmocefRR22im/WMymIsvFa5J0nKR/Th93Szok/9DMWsPAQPbmoaOOKjYPOSFYM8pSU/gGMEYy99GngGeBq4CDc4zLrKk9+SS88pXZy7t5yFpFlo7m10fEacALABHxB2DLXKMya1LjNYIsCeGuu9p/yglrP1mSwkvpYjsBIGk2Sc3BWoynG56ac87J3jy0667FRDB3bu6hmdVcluajfwH+HdhJ0hDwXuDjuUZlNTc+3fD4yl3j0w0DvhK0jGpnJB0bq66T2axZVawpRMQw8BHgAmA18O6IuLLSeZIuSWdVXVqy7zxJD0takt6OLDl2rqRlkn4j6R1T+3FsItUs5djJxmsEWRLCddcVawVOCNYusow+2gtYHhFfB5YCfyFp+wzP/T2g3LReF0bE3PR2ffoa+wLHAvul53wjbbKyGvF0wxP73veyNw9BMREceWTlsmatJkufwlXABkmvBr4N7EEynfakIuIXwJMZ45gHXBER6yJiObAM8LDXGvJ0wy9XOiPpSSdVLv/ii+40ts6QJSmMRcR64GjgqxFxFsnaClP1D5LuSZuXdkj37QasLCmzKt23CUkDkhZJWrRmzZpphNFZPN1wYjwRZJlh9FOfKiaCLbbIPzazZpB19NH7gb8Frk33TfVf5CJgL2AuSf/El9L95SruZb+TRcT8iOiLiL7Zs2dPMYzO08nTDd9009Sahz7xiXzjMmtGWUYfnQScCgxFxPJ0Kc7LpvJiEfHo+Lakb1FMMquA3UuKzgEemcpr2MT6+/s7IgmMq6bz98knYYcdKpcza3dZRh/dFxGnR8Tl6ePlEfHZqbyYpNJmp/eQdFwDLASOlTQzTTp7A3dM5TWss3V3Z68VHHNMsVbghGCWqFhTkLQ3yXDUfYHCQL1Ki+tIuhw4HNhR0irgk8DhkuaSNA2NAKekz3WvpAXAfcB64LSI2FD9j2Od6Le/hX32yV7encVmE8vSfPRdkg/0C4E3kzQnVfweFhHvL7P7O5OUHwI6q9fTpqWa5qEHH6xuKmuzTpWlo3nriLiRZOnO0Yg4j2RyPLO6O/ro7M1De+xRbB5yQjDLJktN4QVJM4AHJf0D8DDgSX+tbp5+GrbfPnt5Nw+ZTV2WmsKZQBdwOnAQcBxwQo4xmQHFGkGWhPCLX/jiMrNaqFhTiIhfAUiKiMhw7afZ1H3ta3D66dnLOwmY1VaW0Ud/RtJBvA3QLWl/4JSI+D95B2ed4aWXYMsqVujwjKRm+cnSfPQV4B3AEwARcTfwphxjsg4x3jyUJSHMnv33SDPo6enlhz/0OhBmecmSFIiIlRvt8jUENiXXXlvdlBOXXTZMV9cfsWbNt4mIwjoQXiDILB9ZRh+tlHQoEJK2JOlwvj/fsKydRGSbgG7cunXF2kNv78TrQHTSlB1m9ZLlX/VU4DSSWUtXkUxmd1qOMVmb2Hnn7DOSfvnLxdFDpc1JXgfCrL6yjD56HPBXMsvkrrvgwAOzl680eqi7u5vR0dGy+82s9rKsvPZ5SdtJ2kLSjZIel3RcPYKz1jHeT5AlITz5ZPZrCrwOhFl9ZWk+entEPAO8k6T56DXAh3ONylrCX/5l9k7jU06Z2oyknbwOhFkjZOloHl9Q50jg8oh4Uh4k3rFWrYLdd69cblwtLi7rtHUgzBopS03hx5IeAPqAGyXNBl7INyxrNuM1giwJ4aGHPOWEWavKssjOOcCfAX0R8RLwHDAv78Cs8Y466r7MzUOHHlpMBHtOutKGmTWzCZuPJL0lIv5b0tEl+0qLXJ1nYNYYzz4L2203/mjfiuVdGzBrL5P1Kfwv4L+Bd5U5FjgptJVquoluvRXe8Ib8YjGzxpkwKUTEJ9N7z4zapi6+GE49NWvpdcBWSOINbxjLMSoza6TJmo/+cbITI+LLtQ/H8rZhA2yeZcxZwcurEL5ozKy9TdbRvG166wM+SDLNxW4k015Ubmy2pjJvXtJElCUhLFhQnIiulC8aM2t/EyaFiDg/Is4HdgQOjIizI+JsktXX5tQrQJu6e+4pDiVduLBy+fHRQ8cc44vG8jI8PExvby8zZsygt7fXs71a01FUGD6SXqOwf0SsSx/PBO6OiD+uQ3yT6uvri0WLFjU6jKZS7Yykzz8PW22VXzxWNDw8zMDAwMtmfe3q6nKytbqTtDgi+sody/Lx8QPgDknnSfokcDtwaS0DtOk766zsM5L++MfFWoETQv0MDk48DbhZs8gyS+qQpP8A/jzddVJE3JVvWJbFypWQtd939mx47LF847HJeRpwawWZxqFExJ3AnTnHYhlVc03BM8/AttvmF4tl52nArRVU0fpsjfSVr2SfkfTii4vNQ04IzcPTgFsrcFJoYk89VUwEZ51Vufx4IhgYyD20umqXETse0WWtoOLoo2bWrqOPZs6EF1/MVnb16mTZy3blETtmtTfd0UdWBwsWFGsFlRLC4GCxVtDOCQE8Yses3qqa8KAaki4hWa3tsYj4k3TfLODfgF5gBHhfRPwhPXYucDKwATg9In6SV2zNYt266oaEtnClbso8YsesvvKsKXwPOGKjfecAN0bE3sCN6WMk7QscC+yXnvMNSZvlGFtDnXpqUiPIkhAeeKCzF6yZaGSOR+yY5SO3pBARvwCe3Gj3PIoXvl0KvLtk/xURsS4ilgPLgEPyiq0Rli4tNg9dfPHkZfv7i4lgn33qE1+z8ogds/rKrfloAq+KiNUAEbFa0k7p/t2A20rKrUr3bULSADAAzf9tsdopJ8bGqrsGoROMdyYPDg6yYsUKuru7GRoaciezWU6apaO53Edh2QaTiJgfEX0R0Td79uycw5qa66/PPuXELbcUawXNmhAaPSS0v7+fkZERxsbGGBkZcUIwy1G9awqPStolrSXsAoxPvLAKKF0Sfg7wSJ1jm5Ynn4Rdd006jys580y48MLcQ6qJjYeEjo6OMpBeCOEPZ7P2U++awkLghHT7BOCakv3HSpopaQ9gb+COOsc2JWeemXzDf+UrKyeEl15KagSNTgjVfPP3kFCzzpLnkNTLgcOBHSWtAj4JfBZYIOlkYAVwDEBE3CtpAXAfsB44LSI25BXbdN1yC7zxjdnKLlsGe+2VbzzVqPabv4eEmnUWX9Gc0f/8D7z2tbBqVeWyF11UzdrH9dXb21t2Uraenh5GRkamXd7Mmp+vaJ6GoaGkeWjbbSdPCAcdBC+8kDQPNWtCgOq/+XtIqFlncVIoo3QZy49/fPKyd9+dJIJFi5I5i5pdtReDeRI3s87ipJBatw4OPDBJBPvvP3nZz3ymOIz0da+rT3y1MpVv/h4SatY5Oj4pfPObxSkn7ppkPbndd4dnn00SQSsPvPE3fzObTEd2NL/0EuyxBzz8cOWyN98Mhx02heDMzJqUO5o38oEPTJ4Qzjij2DzkhGBmnaTeVzQ3hSVLNt231VZJopg1q+7hmJk1jY6sKVxyCRx3XLJ93XVJjeD5550QzMw6sk/BzKyTuU/BzMwycVIwM7MCJwUzMytwUshRoxenMTOrVkcOSa0HL05jZq3INYWceHEaM2tFTgo58eI0ZtaKnBRyUu0U1a3EfSWb8ntibSMiWvZ20EEHRbO67LLLoqurK4DCraurKy677LJGhzYt7fpzTYffE2s1wKKY4HO14R/s07k1c1KISD4senp6QlL09PS0xYdET0/Pyz78xm89PT0Vz23H9yNieu+JWSNMlhQ8zYVVZcaMGZT7m5HE2NjYhOdtPBoLksV92mEth6m+J2aN4mkurGam2lfSzqOx2rn/yDqPk4JVZSrLeUJ7j8aa6nti1oycFKwqU13Os52/TXuJU2srE3U2tMJtuh3N7drx2Yw8QseseTBJR3PH1hTGOz5HR0eJiMI0FB5fng9/mzZrDR07+qi3t5fR0dFN9vf09DAyMjLNyMzMmpdHH5XRzh2fZmZT1bFJoZ07Ps3MpqohSUHSiKRfS1oiaVG6b5akGyQ9mN7vkGcMHkZoZrapRtYU3hwRc0vatc4BboyIvYEb08e5ccenVcMT3lmnaKbmo3nApen2pcC783qh8X/w448/HoAf/OAHjIyMOCFYWR6pZp2kIaOPJC0H/kAyXv3iiJgv6amI2L6kzB8iYpMmJEkDwABAd3f3QeVGEE2mnefgsXx4pJq1m8lGHzUqKewaEY9I2gm4AfgQsDBLUig1lSGp/ge3annCO2s3TTckNSIeSe8fA/4dOAR4VNIuAOn9Y3m8toeiWrU8Us06Sd2TgqQ/krTt+DbwdmApsBA4IS12AnBNHq/vf3CrlkeqWSdpRE3hVcDNku4G7gCui4j/BD4L/IWkB4G/SB/XnP/BrVoeqWadpCOnuRgeHmZwcJAVK1bQ3d3N0NCQ/8HNrGM0XUdzrXjlNTOz6jVdR7OZmTUnJwUzMytwUjAzswInBTMzK3BSMDOzgpYefSRpDbDxnBU7Ao83IJzJOKbsmjGuZowJmjMux5RdI+PqiYjZ5Q60dFIoR9KiiYZaNYpjyq4Z42rGmKA543JM2TVrXG4+MjOzAicFMzMraMekML/RAZThmLJrxriaMSZozrgcU3ZNGVfb9SmYmdnUtWNNwczMpshJwczMClouKUiaJekGSQ+m92WX7JR0hKTfSFom6ZyS/XMl3SZpiaRFkg5phrjSYx9Kj90r6fPNEFN6/J8khaQdGx2TpC9IekDSPZL+XdL204yn0s8uSf+SHr9H0oFZz613TJJ2l/QzSfenf0NnNDqmkuObSbpL0rW1imm6cUnaXtKP0r+n+yX9WRPEdFb6u1sq6XJJW9UipqpEREvdgM8D56Tb5wCfK1NmM+AhYE9gS+BuYN/02H8Bf5luHwn8vEniejPwU2Bm+ninRseUHt8d+AnJRYI7NjomkpX6Nk+3P1fu/CpimfRnL/kb+Q9AwBuA27Oe24CYdgEOTLe3BX7b6JhKjv8j8EPg2unGU6u4gEuBv0u3twS2b/DvbzdgObB1+ngBcGKt3q+st5arKQDzSH6ZpPfvLlPmEGBZRPwuIl4ErkjPAwhgu3T7FcAjTRLXB4HPRsQ6KKxf3eiYAC4EPkLyvtXCtGKKiP+KiPVpuduAOdOIpdLPPh7v9yNxG7C9kjXEs5xb15giYnVE3AkQEc8C95N80DQsJgBJc4CjgG/XIJaaxCVpO+BNwHcAIuLFiHiqkTGlxzYHtpa0OdBF7T6fMmvFpPCqiFgNkN7vVKbMbsDKkserKP5znAl8QdJK4IvAuU0S12uAP5d0u6SbJB3c6Jgk/RXwcETcXYNYahLTRj5A8o1rqrK8zkRlssZYz5gKJPUCBwC3N0FMXyH5YjFWg1hqFdeewBrgu2mz1reVrBnfsJgi4mGSz6QVwGrg6Yj4rxrEVJXN6/2CWUj6KbBzmUODWZ+izL7xb7ofBM6KiKskvY/km8LbmiCuzYEdSKqTBwMLJO0ZaT2y3jFJ6kqf4+0Znyf3mDZ6jUFgPTBcXXTVvc4kZbKcOxXTiSk5KG0DXAWcGRHPNDImSe8EHouIxZIOr0EsNYmL5P/tQOBDEXG7pK+SNGd+olExpf1r84A9gKeAKyUdFxGXTTOmqjRlUoiICT+kJT06XlVOq1zlmllWkbSFj5tDsRp2AjDeAXclVVRpc45rFXB1mgTukDRGMmHWmgbFtBfJH+fdksb33ynpkIj4fYNiGn+OE4B3Am+tlDQrmPR1KpTZMsO59Y4JSVuQJIThiLi6BvFMN6b3An8l6UhgK2A7SZdFxHENjiuAVRExXpP6EUlSaGRMbwOWR8QaAElXA4cCdU0Kde3AqMUN+AIv76j8fJkymwO/I/lQG+/s2S89dj9weLr9VmBxk8R1KvCpdPs1JNVLNTKmjcqNUJuO5um+T0cA9wGzaxBLxZ+dpC28tFPwjmretzrHJOD7wFdq8Tddi5g2KnM4te1onlZcwC+BfdLt84AvNPj393rgXpK+BJH0uX2olr/LTD9DvV+wBm/6K4EbgQfT+1np/l2B60vKHUky+uIhYLBk/xuBxekv63bgoCaJa0uSbwRLgTuBtzQ6po2ea4TaJIXpvk/LSBLmkvT2zWnGs8nrkCToU9NtAV9Pj/8a6KvmfatnTOnfdgD3lLw/RzYypo2e43BqmBRq8PubCyxK36//C+zQBDGdDzxA8jnwA9LRiPW8eZoLMzMraMXRR2ZmlhMnBTMzK3BSMDOzAicFMzMrcFIwM7MCJwWzKZK0Qclsu3dLulPSoen+XiWzyn66pOyOkl6S9K/p4/Mk/VOjYjebiJOC2dQ9HxFzI2J/kjm0Lig59juSq67HHUNyYZJZU3NSMKuN7YA/lDx+HrhfUl/6+H+TTIVs1tSacu4jsxaxtaQlJHP67AK8ZaPjVwDHSvo9sIFkfptd6xqhWZWcFMym7vmImAuQrtr1fUl/UnL8P4FPA48C/1b/8Myq5+YjsxqIiFtJZrWdXbLvRZJ5ts4mmbnUrOm5pmBWA5L+mGQpxidIZrkc9yXgpoh4Ip2C3KypOSmYTd14nwIkM1+eEBEbSj/8I+JePOrIWohnSTUzswL3KZiZWYGTgpmZFTgpmJlZgZOCmZkVOCmYmVmBk4KZmRU4KZiZWcH/B4SN6AshAJB8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(bmi_test, disease_progress_test,  color='black')\n",
    "plt.plot(bmi_test, disease_progress_pred, color='blue', linewidth=3)\n",
    "plt.xlabel('BMI')\n",
    "plt.ylabel('disease progression')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, we looked at a simple linear regression model that depends on a single variable. We generalize this to a Multiple Linear Regression model of the form\n",
    "\n",
    "$y_i = b + a_1 x1_i + a_2 x2_i + ...$\n",
    "\n",
    "In the example below, we use all ten baseline variables to build our model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept:  152.1334841628965\n",
      "Coefficient:  [ -10.01219782 -239.81908937  519.83978679  324.39042769 -792.18416163\n",
      "  476.74583782  101.04457032  177.06417623  751.27932109   67.62538639]\n",
      "Variance score (R2): 0.52\n"
     ]
    }
   ],
   "source": [
    "# Create and fit the model\n",
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(diabetes.data, diabetes.target)\n",
    "\n",
    "# Apply the model (predict the disease progression from BMI using linear model)\n",
    "disease_progress_pred = regr.predict(diabetes.data)\n",
    "\n",
    "print('Intercept: ', regr.intercept_)\n",
    "print('Coefficient: ', regr.coef_)\n",
    "print('Variance score (R2): %.2f' % r2_score(diabetes.target, disease_progress_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combining sklearn with pandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, we've been doing things the hard way - selecting columns from our 2D array of data, saving as new arrays, making sure that we retain multiple dimensions even if we only want one column.\n",
    "\n",
    "A much easier approach is to use pandas and work directly with data frames. This has the additional advantage that you can the dataframe methods to operate on the data. In the following examples, we train using all the independent variables, two variables and each variable individually"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import pandas as pd\n",
    "df = pd.DataFrame(diabetes.data, columns=diabetes.feature_names)\n",
    "disease_progress = pd.DataFrame(diabetes.target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>age</th>\n",
       "      <th>sex</th>\n",
       "      <th>bmi</th>\n",
       "      <th>bp</th>\n",
       "      <th>s1</th>\n",
       "      <th>s2</th>\n",
       "      <th>s3</th>\n",
       "      <th>s4</th>\n",
       "      <th>s5</th>\n",
       "      <th>s6</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.038076</td>\n",
       "      <td>0.050680</td>\n",
       "      <td>0.061696</td>\n",
       "      <td>0.021872</td>\n",
       "      <td>-0.044223</td>\n",
       "      <td>-0.034821</td>\n",
       "      <td>-0.043401</td>\n",
       "      <td>-0.002592</td>\n",
       "      <td>0.019908</td>\n",
       "      <td>-0.017646</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-0.001882</td>\n",
       "      <td>-0.044642</td>\n",
       "      <td>-0.051474</td>\n",
       "      <td>-0.026328</td>\n",
       "      <td>-0.008449</td>\n",
       "      <td>-0.019163</td>\n",
       "      <td>0.074412</td>\n",
       "      <td>-0.039493</td>\n",
       "      <td>-0.068330</td>\n",
       "      <td>-0.092204</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.085299</td>\n",
       "      <td>0.050680</td>\n",
       "      <td>0.044451</td>\n",
       "      <td>-0.005671</td>\n",
       "      <td>-0.045599</td>\n",
       "      <td>-0.034194</td>\n",
       "      <td>-0.032356</td>\n",
       "      <td>-0.002592</td>\n",
       "      <td>0.002864</td>\n",
       "      <td>-0.025930</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        age       sex       bmi        bp        s1        s2        s3  \\\n",
       "0  0.038076  0.050680  0.061696  0.021872 -0.044223 -0.034821 -0.043401   \n",
       "1 -0.001882 -0.044642 -0.051474 -0.026328 -0.008449 -0.019163  0.074412   \n",
       "2  0.085299  0.050680  0.044451 -0.005671 -0.045599 -0.034194 -0.032356   \n",
       "\n",
       "         s4        s5        s6  \n",
       "0 -0.002592  0.019908 -0.017646  \n",
       "1 -0.039493 -0.068330 -0.092204  \n",
       "2 -0.002592  0.002864 -0.025930  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>151.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>75.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>141.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       0\n",
       "0  151.0\n",
       "1   75.0\n",
       "2  141.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "disease_progress.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Build the multiple linear regression model using all independent variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [ -10.01219782 -239.81908937  519.83978679  324.39042769 -792.18416163\n",
      "  476.74583782  101.04457032  177.06417623  751.27932109   67.62538639]\n",
      "Variance score (R2): 0.52\n"
     ]
    }
   ],
   "source": [
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(df, disease_progress)\n",
    "disease_progress_pred = regr.predict(df)\n",
    "\n",
    "print('Intercept: ', regr.intercept_)\n",
    "print('Coefficient: ', regr.coef_[0])\n",
    "print('Variance score (R2): %.2f' % r2_score(disease_progress, disease_progress_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Build the multiple linear regression model using age, sex and bmi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [133.01372901 924.81645876]\n",
      "Variance score (R2): 0.35\n"
     ]
    }
   ],
   "source": [
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(df[[\"age\", \"bmi\"]], disease_progress)\n",
    "disease_progress_pred = regr.predict(df[[\"age\", \"bmi\"]])\n",
    "\n",
    "print('Intercept: ', regr.intercept_)\n",
    "print('Coefficient: ', regr.coef_[0])\n",
    "print('Variance score (R2): %.2f' % r2_score(disease_progress, disease_progress_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Build the ordinary linear regression model using bmi - set aside data for testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept:  [152.91886183]\n",
      "Coefficient:  [938.23786125]\n",
      "Variance score (R2): 0.47\n"
     ]
    }
   ],
   "source": [
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(df[[\"bmi\"]][:-20], disease_progress[:-20])\n",
    "disease_progress_pred = regr.predict(df[[\"bmi\"]][-20:])\n",
    "\n",
    "print('Intercept: ', regr.intercept_)\n",
    "print('Coefficient: ', regr.coef_[0])\n",
    "print('Variance score (R2): %.2f' % r2_score(disease_progress[-20:], disease_progress_pred[-20:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Perform ordinary linear regression, looping over all features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Feature:  age\n",
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [304.18307453]\n",
      "Variance score (R2): 0.04\n",
      "\n",
      "Feature:  sex\n",
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [69.71535568]\n",
      "Variance score (R2): 0.00\n",
      "\n",
      "Feature:  bmi\n",
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [949.43526038]\n",
      "Variance score (R2): 0.34\n",
      "\n",
      "Feature:  bp\n",
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [714.7416437]\n",
      "Variance score (R2): 0.19\n",
      "\n",
      "Feature:  s1\n",
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [343.25445189]\n",
      "Variance score (R2): 0.04\n",
      "\n",
      "Feature:  s2\n",
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [281.78459335]\n",
      "Variance score (R2): 0.03\n",
      "\n",
      "Feature:  s3\n",
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [-639.14527932]\n",
      "Variance score (R2): 0.16\n",
      "\n",
      "Feature:  s4\n",
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [696.88303009]\n",
      "Variance score (R2): 0.19\n",
      "\n",
      "Feature:  s5\n",
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [916.13872282]\n",
      "Variance score (R2): 0.32\n",
      "\n",
      "Feature:  s6\n",
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [619.22282068]\n",
      "Variance score (R2): 0.15\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for feature in df.columns.values:\n",
    "    print(\"Feature: \", feature)\n",
    "    regr = linear_model.LinearRegression()\n",
    "    regr.fit(df[[feature]], disease_progress)\n",
    "    disease_progress_pred = regr.predict(df[[feature]])\n",
    "\n",
    "    print('Intercept: ', regr.intercept_)\n",
    "    print('Coefficient: ', regr.coef_[0])\n",
    "    print('Variance score (R2): %.2f' % r2_score(disease_progress, disease_progress_pred))\n",
    "    print()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using column of data frame as target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import pandas as pd\n",
    "df = pd.DataFrame(diabetes.data, columns=diabetes.feature_names)\n",
    "df['Disease progression'] = diabetes.target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>age</th>\n",
       "      <th>sex</th>\n",
       "      <th>bmi</th>\n",
       "      <th>bp</th>\n",
       "      <th>s1</th>\n",
       "      <th>s2</th>\n",
       "      <th>s3</th>\n",
       "      <th>s4</th>\n",
       "      <th>s5</th>\n",
       "      <th>s6</th>\n",
       "      <th>Disease progression</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.038076</td>\n",
       "      <td>0.050680</td>\n",
       "      <td>0.061696</td>\n",
       "      <td>0.021872</td>\n",
       "      <td>-0.044223</td>\n",
       "      <td>-0.034821</td>\n",
       "      <td>-0.043401</td>\n",
       "      <td>-0.002592</td>\n",
       "      <td>0.019908</td>\n",
       "      <td>-0.017646</td>\n",
       "      <td>151.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-0.001882</td>\n",
       "      <td>-0.044642</td>\n",
       "      <td>-0.051474</td>\n",
       "      <td>-0.026328</td>\n",
       "      <td>-0.008449</td>\n",
       "      <td>-0.019163</td>\n",
       "      <td>0.074412</td>\n",
       "      <td>-0.039493</td>\n",
       "      <td>-0.068330</td>\n",
       "      <td>-0.092204</td>\n",
       "      <td>75.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.085299</td>\n",
       "      <td>0.050680</td>\n",
       "      <td>0.044451</td>\n",
       "      <td>-0.005671</td>\n",
       "      <td>-0.045599</td>\n",
       "      <td>-0.034194</td>\n",
       "      <td>-0.032356</td>\n",
       "      <td>-0.002592</td>\n",
       "      <td>0.002864</td>\n",
       "      <td>-0.025930</td>\n",
       "      <td>141.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.089063</td>\n",
       "      <td>-0.044642</td>\n",
       "      <td>-0.011595</td>\n",
       "      <td>-0.036656</td>\n",
       "      <td>0.012191</td>\n",
       "      <td>0.024991</td>\n",
       "      <td>-0.036038</td>\n",
       "      <td>0.034309</td>\n",
       "      <td>0.022692</td>\n",
       "      <td>-0.009362</td>\n",
       "      <td>206.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.005383</td>\n",
       "      <td>-0.044642</td>\n",
       "      <td>-0.036385</td>\n",
       "      <td>0.021872</td>\n",
       "      <td>0.003935</td>\n",
       "      <td>0.015596</td>\n",
       "      <td>0.008142</td>\n",
       "      <td>-0.002592</td>\n",
       "      <td>-0.031991</td>\n",
       "      <td>-0.046641</td>\n",
       "      <td>135.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        age       sex       bmi        bp        s1        s2        s3  \\\n",
       "0  0.038076  0.050680  0.061696  0.021872 -0.044223 -0.034821 -0.043401   \n",
       "1 -0.001882 -0.044642 -0.051474 -0.026328 -0.008449 -0.019163  0.074412   \n",
       "2  0.085299  0.050680  0.044451 -0.005671 -0.045599 -0.034194 -0.032356   \n",
       "3 -0.089063 -0.044642 -0.011595 -0.036656  0.012191  0.024991 -0.036038   \n",
       "4  0.005383 -0.044642 -0.036385  0.021872  0.003935  0.015596  0.008142   \n",
       "\n",
       "         s4        s5        s6  Disease progression  \n",
       "0 -0.002592  0.019908 -0.017646                151.0  \n",
       "1 -0.039493 -0.068330 -0.092204                 75.0  \n",
       "2 -0.002592  0.002864 -0.025930                141.0  \n",
       "3  0.034309  0.022692 -0.009362                206.0  \n",
       "4 -0.002592 -0.031991 -0.046641                135.0  "
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept:  [152.13348416]\n",
      "Coefficient:  [133.01372901 924.81645876]\n",
      "Variance score (R2): 0.35\n"
     ]
    }
   ],
   "source": [
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(df[['age', 'bmi']], df[['Disease progression']])\n",
    "disease_progress_pred = regr.predict(df[['age', 'bmi']])\n",
    "\n",
    "print('Intercept: ', regr.intercept_)\n",
    "print('Coefficient: ', regr.coef_[0])\n",
    "print('Variance score (R2): %.2f' % r2_score(df[['Disease progression']], disease_progress_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with categorical variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keep in mind that sex is a categorical value and cannot be used for linear regression. You might recall when we looped over the features that the $R^2$ value for dependence of disease progression on sex was exactly zero. We'll need to take a different approach.\n",
    "\n",
    "This was easy to overlook since sex appears as a numerical value in the data set. Male and female were likely assigned the values $\\pm 1$, which then became 0.050680 and -0.044642 after normalizing the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "progression_male = np.array(df['Disease progression'][df['sex'] > 0])\n",
    "progression_female = np.array(df['Disease progression'][df['sex'] < 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "78.26358288372121"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "progression_male.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "75.7441070979943"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "progression_female.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_indResult(statistic=0.9041147550244715, pvalue=0.3664292946519826)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#from scipy import stats\n",
    "stats.ttest_ind(progression_male, progression_female, equal_var=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Submit Ticket\n",
    "If you find anything that needs to be changed, edited, or if you would like to provide feedback or contribute to the notebook, please submit a ticket by contacting us at:\n",
    "\n",
    "Email: consult@sdsc.edu\n",
    "\n",
    "We appreciate your input and will review your suggestions promptly!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}