Regression.ipynb

{
 "cells": [
  {
   "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": "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",
    "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diabetes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make this a little more readable, we can print just the description"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": "code",
   "execution_count": null,
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for feature in df.columns.get_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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using column of data frame as target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "progression_male.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "progression_female.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "stats.ttest_ind(progression_male, progression_female, equal_var=True)"
   ]
  },
  {
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}