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",
    "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",
    "from patsy import dmatrix\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in this lab, we will use Wage data. Let us read in the CSV data ans look at a sample of this data.\n",
    "Wage = pd.read_csv('data/Wage.csv', header=0, na_values='NA')\n",
    "print(Wage.shape)\n",
    "print(Wage.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.8.1 Polynomial Regression and Step Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "We will examine how to fit a polynomial regression model on the wage dataset. As all the techniques, \n",
    "we have multiple ways to do this. Here I will use sklearn as we alreadly used statsmodel.api before in Chapter 3. \n",
    "If you are looking for more built-in functions around p-value, significance, confidence intervie, etc., \n",
    "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 \n",
    "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/\n",
    "\"\"\"\n",
    "\n",
    "n_deg = 4\n",
    "X = Wage.age\n",
    "y = Wage.wage\n",
    "X = X.values.reshape(X.shape[0], 1)\n",
    "y = y.values.reshape(y.shape[0], 1)\n",
    "\n",
    "polynomial_features= PolynomialFeatures(degree=n_deg)\n",
    "X_poly = polynomial_features.fit_transform(X)\n",
    "\n",
    "reg = LinearRegression()\n",
    "reg.fit(X_poly, y)\n",
    "\n",
    "# get coefficients and compare with the numbers \n",
    "print(reg.intercept_)\n",
    "print(reg.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we now create a grid of values for age at which we want predictionsm and the call the generic predict() function \n",
    "# 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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Next we need to decide the order of the polynomial.\n",
    "In the book, the authors did this by using hypothesis testing.  ANOVA using F-test was explanied. \n",
    "In order to use the ANOVA function, two models $M_1$ and $M_2$ must be nested model: \n",
    "the predictors in $M_1$ must be a subset of the predictors in $M_2$. \n",
    "statsmodel.api has a nice built-in function to do that. \n",
    "\n",
    "As an alternative, we could choose the polynomial degree using cross-validation, as discussed in before. \n",
    "Actually, the cross-validation approach is more commonly used in practice. \n",
    "\"\"\"\n",
    "\n",
    "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sm.stats.anova_lm(fit1, fit2, fit3, fit4, fit5, type=1))\n",
    "\n",
    "\"\"\"\n",
    "The above table, we fit five different models and sequentially compare the simpler model to the more complex model.\n",
    "The summary above shows the quadratic model fit2 is significantly better than fit1 at p value of $2.36*10^{-32}$.\n",
    "Similarly, the cubic model is significnatly better than the quadratic model ($p = 1.68 * 10^{-3}$).\n",
    "The p-value comparing the cubic and degree-4 polynomials, fit3 and fit4, is approximately 0.05 \n",
    "while the degree-5 polynomial fit5 seems unnecessary because its p-value is 0.37. \n",
    "Hence, either a cubic or a quartic polynomial appear to provide a reasonable fit to the data, \n",
    "but lower- or higher-order models are not justified.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in the book, the authors also discussed logistic regression and the polynomial terms. \n",
    "# in python, sm.GLM function provided some functions similar to glm() in R.\n",
    "logistic_model = sm.GLM ((y>250), X4, family=sm.families.Binomial())\n",
    "logistic_fit = logistic_model.fit()\n",
    "print(logistic_fit.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in python, we could use the pd.cut() function to fit a step function.\n",
    "age_cut, bins = pd.cut(Wage.age, bins=4, retbins=True, right=True)\n",
    "age_cut.value_counts(sort=False)\n",
    "\n",
    "\"\"\" \n",
    "Here cut() automatically picked the cutpoints at 33.5, 49, and 64.5 years of age. \n",
    "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]). \n",
    "Note in the following code, I manually added a constant column and dropped the lowest value bin (17.938, 33.5] dummy variable.\n",
    "\"\"\""
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in order to fit regression splines in python, we use the spatsy library. \n",
    "# from patsy import dmatrix\n",
    "\n",
    "\"\"\" \n",
    "In the content of section 7.4, we saw that regression splines can be fit by constructing an appropriate matrix of basis functions. \n",
    "The bs() function generates the entire matrix of bs() basis functions for splines with the specified set of knots. \n",
    "By default, cubic splines are produced. Here we have prespecified knots at ages 25, 40, and 60. \n",
    "This produces a spline with six basis functions. \n",
    "\"\"\"\n",
    "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": [
    "# now we can fit the model using the spline basis functions\n",
    "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 using 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 make some 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.3 GAMs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we now fit a GAM to predict wage using natural spline functions of year and age, treating education as a qualitative (i.e. categorical) predictor.\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)\n",
    "education_dummies = education_dummies.drop([education_dummies.columns[0]], axis = 1)\n",
    "\n",
    "# we concatenate all the predictors\n",
    "x_all = pd.concat([age_basis, year_basis, education_dummies], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# fit a model and print the summary\n",
    "gam1_fit = sm.OLS(Wage.wage, x_all).fit()\n",
    "print(gam1_fit.summary())\n",
    "\n",
    "\"\"\" \n",
    "We could apply similar analysis procedure to this analysis, \n",
    "such as ANOVA, construction of a classification model and visually inspecting the model performance.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of Chapter 7"
   ]
  }
 ],
 "metadata": {
  "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
}