Chapter_7_sec_7.8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7.8 Lab: Non-linear Modeling "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "import pandas as pd \n",
    "\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "import statsmodels.api as sm\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "In this lab, we will use Wage data. Let us read in the CSV data ans look at a sample of this data.\n",
    "\"\"\"\n",
    "Wage = pd.read_csv('data/Wage.csv', header=0, na_values='NA')\n",
    "print Wage.shape\n",
    "Wage.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.8.1 Polynomial Regression and Step Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will examine how to fit a polynomial regression model on the wage dataset. As all the techniques, we have multiple ways to do this. Here I will use sklearn as we alreadly used statsmodel.api before in Chapter 3.  If you are looking for more built-in functions around p-value, significance, confidence intervie, etc., I would recommend to use statsmodel.api. \n",
    "\n",
    "But scikit-learn does not have built error estimates for doing inference. But this problem forces us to think about a more general method to find Confidence Interview (key word: Bootstrap) \n",
    "\n",
    "Numpy also has a nice function to do ploynomial regression: https://www.ritchieng.com/machine-learning-polynomial-regression/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_deg = 4\n",
    "X = Wage.age\n",
    "y = Wage.wage\n",
    "X = X.reshape(X.shape[0], 1)\n",
    "y = y.reshape(y.shape[0], 1)\n",
    "\n",
    "polynomial_features= PolynomialFeatures(degree=n_deg)\n",
    "X_poly = polynomial_features.fit_transform(X)\n",
    "\n",
    "\n",
    "reg = LinearRegression()\n",
    "reg.fit(X_poly, y)\n",
    "\n",
    "# get coefficients and compare with the numbers as the end of page 288.\n",
    "print reg.intercept_, reg.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now create a grid of values for age at which we want predictionsm and the call the generic predict() function "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate a sequence of age values spanning the range\n",
    "age_grid = np.arange(Wage.age.min(), Wage.age.max()).reshape(-1,1)\n",
    "\n",
    "# generate test data use PolynomialFeatures and fit_transform\n",
    "X_test = PolynomialFeatures(degree=n_deg).fit_transform(age_grid)\n",
    "\n",
    "# predict the value of the generated ages\n",
    "y_pred = reg.predict(X_test)\n",
    "\n",
    "# creating plots\n",
    "plt.plot(age_grid, y_pred, color='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Decide on the polynomial to use. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the book, the authors did this by using hypothesis testing. ANOVA using F-test was explanied. In order\n",
    "to use the ANOVA function, $M_1$ and $M_2$ must be nested model: the predictors in $M_1$ must be a subset of the predictors in $M_2$. statsmodel.api has a nice built-in function to do that. \n",
    "\n",
    "As an alternative to using hypothesis tests and ANOVA, we could choose the polynomial degree using cross-validation, as discussed in before. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X1 = PolynomialFeatures(1).fit_transform(X)\n",
    "X2 = PolynomialFeatures(2).fit_transform(X)\n",
    "X3 = PolynomialFeatures(3).fit_transform(X)\n",
    "X4 = PolynomialFeatures(4).fit_transform(X)\n",
    "X5 = PolynomialFeatures(5).fit_transform(X)\n",
    "fit1 = sm.GLS(y, X1).fit()\n",
    "fit2 = sm.GLS(y, X2).fit()\n",
    "fit3 = sm.GLS(y, X3).fit()\n",
    "fit4 = sm.GLS(y, X4).fit()\n",
    "fit5 = sm.GLS(y, X5).fit()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "print(sm.stats.anova_lm(fit1, fit2, fit3, fit4, fit5, typ=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The row of the above take shows the fit1 to the quadratic model fit2 is $2.36*10^{-32}$, indicating that a quadratic model is significant informative to a linear model. Similarly, the cubic model is significnat informative to a quadratic model ($p = 1.68 * 10^{-2}$).Hence, either a cubic or a quartic polynomial appear to provide a reasonable fit to the data, but lower- or higher-order models are not justified."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the book, the authors also discussed logistic regression and the polynomial terms. In python, sm.GLM function provided some functions similar to glm() in R."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logistic_model = sm.GLM ((y>250), X4, family=sm.families.Binomial())\n",
    "logistic_fit = logistic_model.fit()\n",
    "print(logistic_fit.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to fit a step function, we use the cut() function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "age_cut, bins = pd.cut(Wage.age, bins=4, retbins=True, right=True)\n",
    "age_cut.value_counts(sort=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here cut() automatically picked the cutpoints at 33.5, 49, and 64.5 years of age. We could also have specified our own cutpoints directly using the breaks option  (set bins into a sequence of scalars, e.g. [0, 10, 20, 40, 100]). Note in the following code, I manually added a constant column and dropped the lowest value bin (17.938, 33.5] dummy variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "age_cut_dummies = pd.get_dummies(age_cut)\n",
    "age_cut_dummies = sm.add_constant(age_cut_dummies)\n",
    "fit_age_cut = sm.GLM(Wage.wage, age_cut_dummies.drop(age_cut_dummies.columns[1], axis=1)).fit()\n",
    "print(fit_age_cut.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.8.2 Splines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to fit regression splines in python, we use the spatsy library. In Section 7.4, we saw that regression splines can be fit by constructing an appropriate matrix of basis functions. The bs() function generates the entire matrix of bs() basis functions for splines with the specified set of knots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from patsy import dmatrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have prespecified knots at ages 25, 40, and 60. This produces a spline with six basis functions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "age_grid = np.arange(Wage.age.min(), Wage.age.max()).reshape(-1,1)\n",
    "spline_basis1 = dmatrix(\"bs(Wage.age, knots=(25,40,60), degree=3, include_intercept=False)\", {\"Wage.age\": Wage.age}, return_type='dataframe')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spline_fit1 = sm.GLM(Wage.wage, spline_basis1).fit()\n",
    "spline_fit1.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Another approach is to fix the degree of freedom and let the code to automatically choose the knots.\n",
    "spline_basis2 = dmatrix(\"bs(Wage.age, df=6, include_intercept=False)\",\n",
    "                        {\"Wage.age\": Wage.age}, return_type='dataframe')\n",
    "spline_fit2 = sm.GLM(Wage.wage, spline_basis2).fit()\n",
    "spline_fit2.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Package patsy also has a nice function to do natural spline, cr()\n",
    "spline_basis3 = dmatrix(\"cr(Wage.age, df=4)\", {\"Wage.age\": Wage.age}, return_type='dataframe')\n",
    "spline_fit3 = sm.GLM(Wage.wage, spline_basis3).fit()\n",
    "spline_fit3.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Finally, let us makesome predictions\n",
    "pred1 = spline_fit1.predict(dmatrix(\"bs(age_grid, knots=(25,40,60), include_intercept=False)\",{\"age_grid\": age_grid}, return_type='dataframe'))\n",
    "pred2 = spline_fit2.predict(dmatrix(\"bs(age_grid, df=6, include_intercept=False)\",{\"age_grid\": age_grid}, return_type='dataframe'))\n",
    "pred3 = spline_fit3.predict(dmatrix(\"cr(age_grid, df=4)\", {\"age_grid\": age_grid}, return_type='dataframe'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the splines and error bands\n",
    "plt.scatter(Wage.age, Wage.wage, facecolor='None', edgecolor='k', alpha=0.1)\n",
    "plt.plot(age_grid, pred1, color='r', label='Cubic spine with knots at [25, 40, 60]')\n",
    "plt.plot(age_grid, pred2, color='g', label='Cubic spine with df=6')\n",
    "plt.plot(age_grid, pred3, color='b', label='Natural spline df=4')\n",
    "plt.legend()\n",
    "plt.xlim(15,85)\n",
    "plt.ylim(0,350)\n",
    "plt.xlabel('age')\n",
    "plt.ylabel('wage')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.8.2 GAMs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we construct the basis of all the perdictors. \n",
    "age_basis = dmatrix(\"cr(Wage.age, df=5)\", {\"Wage.age\": Wage.age}, return_type='dataframe')\n",
    "year_basis = dmatrix(\"cr(Wage.year, df=4)\", {\"Wage.year\": Wage.year}, return_type='dataframe').drop (['Intercept'], axis = 1)\n",
    "education_dummies = pd.get_dummies(Wage.education).drop([education_dummies.columns[0]], axis = 1)\n",
    "\n",
    "x_all = pd.concat([age_basis, year_basis, education_dummies], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gam1_fit = sm.OLS(Wage.wage, x_all).fit()\n",
    "gam1_fit.summary()"
   ]
  }
 ],
 "metadata": {
  "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
}