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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in Python, module can be imported by a command similar to 'import numpy as np'. \n",
    "# it is a good practice to maintain a section at the beginning of the notebook to import all necessary modules.\n",
    "# for new module, could use pip to install it. \n",
    "# for example 'pip install numpy'\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd \n",
    "import math\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
    "from statsmodels.graphics.regressionplots import *\n",
    "from sklearn import datasets, linear_model\n",
    "from patsy import dmatrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# since in Python, there is no default MASS module and Boston dataset, I will read in the Boston dataset from CSV. The data is in the ./data folder.\n",
    "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": [
    "# use the commands we learned in the previous chapeter to exame the data.\n",
    "list(Boston) # or Boston.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Boston.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Boston.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to have similar formula notation as R, use the following import. \n",
    "# One thing to note is that the ' ' for the formula part in fitting step and the '.fit()' at the end.\n",
    "# import statsmodels.formula.api as smf, we would use smf to call the model. Of course, there are other ways to run linear regression in pythin, such as sklearn.\n",
    "lm = smf.ols ('medv~lstat', data = Boston).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lm.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use dir() to get a list of all the attributes an object has\n",
    "dir(lm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can try a few \n",
    "print(lm.params)\n",
    "print(lm.conf_int())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# provide prediction for 3 observations\n",
    "lm.predict(pd.DataFrame({'lstat':[5, 10, 15]}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the fitted line, we only take two extreme points to make the plot\n",
    "X_new = pd.DataFrame({'lstat': [Boston.lstat.min(), Boston.lstat.max()]})\n",
    "preds = lm.predict(X_new)\n",
    "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the statistics of the linear regression mostly stored in lm.get_influence(), for example, the cookdistances, leverage.\n",
    "dir(lm.get_influence())\n",
    "# for example, the following identifies the observation with the largest leverage \n",
    "np.argmax(lm.get_influence().hat_matrix_diag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from statsmodels.graphics.regressionplots import * just as a reference\n",
    "plot_leverage_resid2(lm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# as mentioned above. For machine learning models, sklearn is the most common used module, but sklearn is a little bit less on statistics.\n",
    "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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# we can still use smg.ols to run multiple linear regression.\n",
    "lm = smf.ols ('medv~lstat+age', data = Boston).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lm.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.\n",
    "all_columns = \"+\".join(Boston.columns.difference([\"medv\"]))\n",
    "my_formula = \"medv~\" + all_columns\n",
    "lm = smf.ols(my_formula, data=Boston).fit()\n",
    "print(lm.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "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.\n",
    "# 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",
    "    \n",
    "print(vif_coeff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.6.4 Interaction Terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we use * to add interaction terms\n",
    "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# there is anova function built in already in statsmodels.  \n",
    "# if you know what to do, use the key words to google it and likely you will find a very good answer. \n",
    "# here we compare the models with one order of stat and two orders of stats. \n",
    "# by looking at the p value that will reject the null hypothesis that the coefficent of lstat**2 equals 0.\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 \n",
    "\n",
    "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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# let us create dummy variables using get_dummies, then exclude the first dummy column\n",
    "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# then the model buliding will be the same with all numerrical variables.\n",
    "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": [
    "# the interpretation of the coefficients are holding everything fixed, Medium shelve location is associated with an average\n",
    "# increase of sale around 2.0046. \n",
    "lm_carseats_dummy.summary() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compapre the summary of two models, one with explicit encoding of dummy varible, while the other used the built-in function.\n",
    "lm_carseats_wo_dummy = smf.ols('Sales ~ Income + Advertising + Price + Age + C(ShelveLoc)', \n",
    "                            data = Carseats).fit()\n",
    "lm_carseats_wo_dummy.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.6.7 Writing Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let us write a simple function to print current time. \n",
    "# yhe key word in Python for user defined function is 'def'. \n",
    "# pay attention to the ':'. The difference betwwen R (others) and Python is that Python \n",
    "# forces you to obey its indentation rules. For example, the following function won't work because of the extra space in front of 'print'.\n",
    "def print_current_time_wrong():\n",
    "    from datetime import datetime # this is very bad practice !!! \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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of Chapter 3."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "interpreter": {
   "hash": "4548a0e672c5b3a287feee7b2962606840aa548749d1830ef724408652b0c250"
  },
  "kernelspec": {
   "display_name": "Python 2.7.16 64-bit ('base': conda)",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}