Chapter_3_sec_6.1_6.7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.6 Lab: Linear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.6.1 Libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In Python, module can be imported by a command similar to 'import numpy as np'. It is a good practice to maintain a section at the beginning of the notebook to import all necessary modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd \n",
    "import math\n",
    "\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.graphics.regressionplots import *\n",
    "from sklearn import datasets, linear_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Since in Python, there is no default MASS module and Boston dataset, I will read in the Boston dataset from CSV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Boston = pd.read_csv('data/Boston.csv', header=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.6.2 Simple Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(Boston)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Boston.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Boston.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To have similar formula notation as R, use the following import. One thing to note is that the ' ' for the formula part in fitting step and the '.fit()' at the end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm = smf.ols ('medv~lstat', data = Boston).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print lm.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use dir() to get a list of all the attributes an object has"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir(lm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print lm.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print lm.conf_int()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm.predict(pd.DataFrame({'lstat':[5, 10, 15]}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_new = pd.DataFrame({'lstat': [Boston.lstat.min(), Boston.lstat.max()]})\n",
    "print X_new\n",
    "preds = lm.predict(X_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Boston.plot(kind='scatter', x='lstat', y='medv')\n",
    "plt.plot(X_new, preds, c='red', linewidth=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ((ax1, ax2), (ax3, ax4))= plt.subplots(2, 2)\n",
    "ax1.plot(Boston.lstat, lm.predict(),'ro')\n",
    "ax2.plot(lm.predict(), lm.resid, 'go')\n",
    "ax3.plot(lm.predict(), lm.resid_pearson, 'bo')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### the statistics of the linear regression mostly stored in lm.get_influence(), for example, the cookdistances, leverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir(lm.get_influence())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The following identifies the observation with the largest leverage "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.argmax(lm.get_influence().hat_matrix_diag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.graphics.regressionplots import *\n",
    "plot_leverage_resid2(lm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For machine learning models, sklearn is the most common used module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = pd.DataFrame(Boston.lstat)\n",
    "y = Boston.medv\n",
    "print x.shape\n",
    "\n",
    "model = linear_model.LinearRegression()\n",
    "model.fit(x, y)\n",
    "print model.intercept_\n",
    "print model.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.6.3 Multiple Linear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To do multiple linear regression, we can use those two modules again. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lm = smf.ols ('medv~lstat+age', data = Boston).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print lm.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### If we want to use all the variable. We can use the following trick to manually construct the list. In Python, most of time, you have to manully construct the variable list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_columns = \"+\".join(Boston.columns.difference([\"medv\"]))\n",
    "my_formula = \"medv~\" + all_columns\n",
    "lm = smf.ols(my_formula, data=Boston).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unlike R, Python is not fully up speeded to all the statistics. If you want to have the VIF of the variables in LM, you have to code a little bit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from patsy import dmatrices\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y, X = dmatrices(my_formula, data=Boston, return_type='dataframe')\n",
    "vif_coeff = {}\n",
    "for i in range(X.shape[1]):\n",
    "    vif_coeff[X.columns[i]] = variance_inflation_factor(np.array(X.dropna()),i)\n",
    "print vif_coeff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.6.4 Interaction Terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm = smf.ols('medv~lstat * age', data=Boston).fit()\n",
    "print lm.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.6.5 Non-linear Transformations of the Predictors "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm_order1 = smf.ols('medv~ lstat', data=Boston).fit()\n",
    "lm_order2 = smf.ols('medv~ lstat+ I(lstat ** 2.0)', data=Boston).fit()\n",
    "print lm_order2.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ((ax1, ax2), (ax3, ax4))= plt.subplots(2, 2)\n",
    "ax1.plot(Boston.lstat, lm_order1.predict(),'ro')\n",
    "ax3.plot(lm_order1.predict(), lm_order1.resid, 'go')\n",
    "ax4.plot(lm_order1.predict(), lm_order1.resid_pearson, 'bo')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### if we added in the second order, we can see the residues are more random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ((ax1, ax2), (ax3, ax4))= plt.subplots(2, 2)\n",
    "ax1.plot(Boston.lstat,  lm_order2.predict(),'ro')\n",
    "ax2.plot(Boston.lstat ** 2.0,  lm_order2.predict(),'ro')\n",
    "ax3.plot(lm_order2.predict(), lm_order2.resid, 'go')\n",
    "ax4.plot(lm_order2.predict(), lm_order2.resid_pearson, 'bo')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### There is anova function built in already in statsmodels.  If you know what to do, use the key words to google it and likely you will find a very good answer. Here we compare the models with one order of stat and two orders of stats. By looking at the p value that will reject the null hypothesis that the coefficent of lstat**2 equals 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "table = sm.stats.anova_lm(lm_order1, lm_order2)\n",
    "print table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm_log = smf.ols('medv~ np.log(rm)', data=Boston).fit()\n",
    "lm_log.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.6.6 Qualitative Predictors "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### I prepared the Carseats file from .Rdata. And it is saved under the data folder.  Let us load them in and explore this dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Carseats = pd.read_csv('data/Carseats.csv', header=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(Carseats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Carseats.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Carseats.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm_carseats = smf.ols('Sales ~ Income + Advertising + Price + Age', data = Carseats).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm_carseats.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### let us create dummy variables using get_dummies, then exclude the first dummy column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ShelveLoc_dummies = pd.get_dummies(Carseats.ShelveLoc, prefix='ShelveLoc').iloc[:,1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Carseats_dummy = pd.concat([Carseats, ShelveLoc_dummies], axis=1)\n",
    "Carseats_dummy.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Then the model buliding will be the same with all numerrical variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm_carseats_dummy = smf.ols('Sales ~ Income + Advertising + Price + Age + ShelveLoc_Good + ShelveLoc_Medium', \n",
    "                            data = Carseats_dummy).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm_carseats_dummy.summary() \n",
    "# The interpretation of the coefficients are holding everything fixed, Medium shelve location is associated with an average\n",
    "# increase of sale around 2.0046. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### statsmodels has the built-in support for categorical variables, just add C() to indicate that the variable will be categorical. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lm_carseats_wo_dummy = smf.ols('Sales ~ Income + Advertising + Price + Age + C(ShelveLoc)', \n",
    "                            data = Carseats).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm_carseats_wo_dummy.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compapre the summary of above two models, one with explicit encoding of dummy varible, while the other used the built-in function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.6.7 Writing Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let us write a simple function to print current time. The key word in Python for user defined function is 'def'. Pay attention to the ':'. The difference betwwen R (others) and Python is that Python forces you to obey its indentation rules. For example, the following function won't work because of the extra space in front of 'print'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_current_time_wrong():\n",
    "    from datetime import datetime\n",
    "     print str(datetime.now())    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def print_current_time():\n",
    "    from datetime import datetime\n",
    "    print str(datetime.now()) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_current_time()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}