ordinary_least_squares_in_python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ordinary Least Squares in Python\n",
    "\n",
    "Linear Regression, also called Ordinary Least-Squares (OLS) Regression, is probably the most commonly used technique in Statistical Learning. It is also the oldest, dating back to the eighteenth century and the work of [Carl Friedrich Gauss](http://en.wikipedia.org/Carl_Freidrich_Gauss) and [Adrien-Marie Legendre](http://en.wikipedia.org/Adrien-Marie_Legendre). It is also one of the easier and more intuitive techniques to understand, and it provides a good basis for learning more advanced concepts and techniques. This posting explains how to perform linear regression using the [statsmodels](http://statsmodels.sourceforge.net/) Python package, we will discuss the single variable case and defer multiple regression to a future post.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statsmodels\n",
    "\n",
    "Statsmodel is a Python library designed for more statistically-oriented approaches to data analysis, with an emphasis on econometric analyses. It integrates well with the [pandas](http://pandas.pydata.org) and [numpy](http://numpy.org) libraries we covered in a previous post. It also has built in support for many of the statistical tests to check the quality of the fit and a dedicated set of plotting functions to visualize and diagnose the fit. [Scikit-learn](http://scikit-learn.org/stable/) also has support for linear regression, including many forms of regularized regression lacking in statsmodels, but it lacks the rich set of statistical tests and diagnostics that have been developed for linear models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Regression and Ordinary Least Squares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear regression one of the simplest and most commonly used modeling techniques. It makes very strong assumptions about the relationship between the predictor variables (the X) and the response (the Y). It assumes that this relationship takes the form:\n",
    "\n",
    "$y = \\beta_0 + \\beta_1 * x$\n",
    "\n",
    "Ordinary Least Squares is the simplest and most common estimator in which the two $\\beta$s are chosen to minimize the square of the distance between the predicted values and the actual values. Even though this model is quite rigid and often does not reflect the true relationship, this still remains a popular approach for several reasons. For one, it is computationally cheap to calculate the coefficients. It is also easier to interpret than more sophisticated models, and in situations where the goal is understanding a simple model in detail, rather than estimating the response well, they can provide insight into what the model captures. Finally, in situations where there is a lot of noise, it may be hard to find the true functional form, so a constrained model can perform quite well compared to a complex model which is more affected by noise.\n",
    "\n",
    "The resulting model is represented as follows:\n",
    "\n",
    "$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_1 * x$\n",
    "\n",
    "Here the hats on the variables represent the fact that they are estimated from the data we have available. The $\\beta$s are termed the parameters of the model or the coefficients. $\\beta_0$ is called the constant term or the intercept."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"skip-to-python\"></div>\n",
    "##Ordinary Least Squares Using Statsmodels\n",
    "\n",
    "The statsmodels package provides several different classes that provide different options for linear regression. Getting started with linear regression is quite straightforward with the OLS module. \n",
    "\n",
    "To start with we load the [Longley dataset](http://stat.ethz.ch/R-manual/R-patched/library/datasets/html/longley.html) of US macroeconomic data from the [Rdatasets](http://vincentarelbundock.github.io/Rdatasets/) website."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# load numpy and pandas for data manipulation\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# load statsmodels as alias ``sm``\n",
    "import statsmodels.api as sm\n",
    "\n",
    "# load the longley dataset into a pandas data frame - first column (year) used as row labels\n",
    "df = pd.read_csv('longley.csv', index_col=0)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.isna().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the variable Total Derived Employment (``'Employed'``) as our response ``y`` and Gross National Product (``'GNP'``) as our predictor ``X``. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We take the single response variable and store it separately. We also add a constant term so that we fit the intercept of our linear model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = df.Employed  # response/Traget\n",
    "X = df.GNP  # predictor/input\n",
    "X = sm.add_constant(X)  # Adds a constant term to the predictor\n",
    "X.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we perform the regression of the predictor on the response, using the ``sm.OLS`` class and and its initialization ``OLS(y, X)`` method. This method takes as an input two array-like objects: ``X`` and ``y``. In general, ``X`` will either be a numpy array or a pandas data frame with shape ``(n, p)`` where ``n`` is the number of data points and ``p`` is the number of predictors. ``y`` is either a one-dimensional numpy array or a pandas series of length ``n``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "est = sm.OLS(y,X)\n",
    "est"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then need to fit the model by calling the OLS object's ``fit()`` method. Ignore the warning about the kurtosis test if it appears, we have only 16 examples in our dataset and the test of the kurtosis is valid only if there are more than 20 examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "est = est.fit()\n",
    "est.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After visualizing the relationship we will explain the summary. First, we need the coefficients of the fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make sure that graphics appear inline in the iPython notebook\n",
    "%pylab inline\n",
    "\n",
    "# We pick 100 hundred points equally spaced from the min to the max\n",
    "X_prime = np.linspace(X.GNP.min(), X.GNP.max(), 100)[:, np.newaxis]\n",
    "X_prime = sm.add_constant(X_prime)  # add constant as we did before\n",
    "X_prime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we calculate the predicted values\n",
    "y_hat = est.predict(X_prime)\n",
    "\n",
    "y_hat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,9))\n",
    "plt.scatter(X.GNP, y, alpha=0.3)  # Plot the raw data\n",
    "plt.xlabel(\"Gross National Product\")\n",
    "plt.ylabel(\"Total Employment\")\n",
    "plt.plot(X_prime[:, 1], y_hat, 'r', alpha=0.9)  # Add the regression line, colored in red"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Statsmodels also provides a formulaic interface that will be familiar to users of R. Note that this requires the use of a different api to statsmodels, and the class is now called ``ols`` rather than ``OLS``. The argument ``formula`` allows you to specify the response and the predictors using the column names of the input data frame ``data``. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# import formula api as alias smf\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "# formula: response ~ predictors\n",
    "est = smf.ols(formula='Employed ~ GNP', data=df).fit()\n",
    "est.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This summary provides quite a lot of information about the fit. The pats of the table we think are the most important are bolded in the description below.\n",
    "\n",
    "The left part of the first table provides basic information about the model fit:\n",
    "\n",
    "<table>\n",
    "<thead>\n",
    "<tr>\n",
    "<th>Element</th>\n",
    "<th>Description</th>\n",
    "</tr>\n",
    "</thead>\n",
    "<tbody>\n",
    "<tr>\n",
    "<td>Dep. Variable</td>\n",
    "<td>Which variable is the response in the model</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Model</td>\n",
    "<td>What model you are using in the fit</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Method</td>\n",
    "<td>How the parameters of the model were calculated</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>No. Observations</td>\n",
    "<td>The number of observations (examples)</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>DF Residuals</td>\n",
    "<td>Degrees of freedom of the residuals. Number of observations - number of parameters</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>DF Model</td>\n",
    "<td>Number of parameters in the model (not including the constant term if present)</td>\n",
    "</tr>\n",
    "</tbody>\n",
    "</table>\n",
    "\n",
    "The right part of the first table shows the goodness of fit\n",
    "\n",
    "<table>\n",
    "<thead>\n",
    "<tr>\n",
    "<th>Element</th>\n",
    "<th>Description</th>\n",
    "</tr>\n",
    "</thead>\n",
    "<tbody>\n",
    "<tr>\n",
    "<td>**R-squared**</td>\n",
    "<td>The [coefficient of determination](http://en.wikipedia.org/wiki/Coefficient_of_determination). A statistical measure of how well the regression line approximates the real data points</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>**Adj. R-squared**</td>\n",
    "<td>The above value adjusted based on the number of observations and the degrees-of-freedom of the residuals</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>F-statistic</td>\n",
    "<td>A measure how significant the fit is. The mean squared error of the model divided by the mean squared error of the residuals</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Prob (F-statistic)</td>\n",
    "<td>The probability that you would get the above statistic, given the null hypothesis that they are unrelated</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Log-likelihood</td>\n",
    "<td>The log of the likelihood function.</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>AIC</td>\n",
    "<td>The [Akaike Information Criterion](http://en.wikipedia.org/wiki/Akaike_information_criterion). Adjusts the log-likelihood based on the number of observations and the complexity of the model.</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>BIC</td>\n",
    "<td>The [Bayesian Information Criterion](http://en.wikipedia.org/wiki/Bayesian_information_criterion). Similar to the AIC, but has a higher penalty for models with more parameters.</td>\n",
    "</tr>\n",
    "</tbody>\n",
    "</table>\n",
    "\n",
    "The second table reports for each of the coefficients\n",
    "\n",
    "<table>\n",
    "<thead>\n",
    "<tr>\n",
    "<th></th>\n",
    "<th>Description</th>\n",
    "</tr>\n",
    "</thead>\n",
    "<tbody>\n",
    "<tr>\n",
    "<td></td>\n",
    "<td>The name of the term in the model</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>coef</td>\n",
    "<td>The estimated value of the coefficient</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>std err</td>\n",
    "<td>The basic standard error of the estimate of the coefficient. More sophisticated errors are also available.</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>t</td>\n",
    "<td>The t-statistic value. This is a measure of how statistically significant the coefficient is.</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><strong>P &gt; |t|</strong></td>\n",
    "<td>P-value that the null-hypothesis that the coefficient = 0 is true. If it is less than the confidence level, often 0.05, it indicates that there is a statistically significant relationship between the term and the response.</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>**[95.0% Conf. Interval]**</td>\n",
    "<td>The lower and upper values of the 95% confidence interval</td>\n",
    "</tr>\n",
    "</tbody>\n",
    "</table>\n",
    "\n",
    "Finally, there are several statistical tests to assess the distribution of the residuals \n",
    "\n",
    "<table>\n",
    "<thead>\n",
    "<tr>\n",
    "<th>Element</th>\n",
    "<th>Description</th>\n",
    "</tr>\n",
    "</thead>\n",
    "<tbody>\n",
    "<tr>\n",
    "<td>Skewness</td>\n",
    "<td>A measure of the symmetry of the data about the mean. Normally-distributed errors should be symmetrically distributed about the mean (equal amounts above and below the line).</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Kurtosis</td>\n",
    "<td>A measure of the shape of the distribution. Compares the amount of data close to the mean with those far away from the mean (in the tails).</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Omnibus</td>\n",
    "<td>D'Angostino's test. It provides a combined statistical test for the presence of skewness and kurtosis.</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Prob(Omnibus)</td>\n",
    "<td>The above statistic turned into a probability</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Jarque-Bera</td>\n",
    "<td>A different test of the skewness and kurtosis</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Prob (JB)</td>\n",
    "<td>The above statistic turned into a probability</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Durbin-Watson</td>\n",
    "<td>A test for the presence of autocorrelation (that the errors are not independent.) Often important in time-series analysis</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Cond. No</td>\n",
    "<td>A test for multicollinearity (if in a fit with multiple parameters, the parameters are related with each other).</td>\n",
    "</tr>\n",
    "</tbody>\n",
    "</table>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a final note, if you don't want to include a constant term in your model, you can exclude it using the minus operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Fit the no-intercept model\n",
    "est_no_int = smf.ols(formula='Employed ~ GNP - 1', data=df).fit()\n",
    "\n",
    "# We pick 100 hundred points equally spaced from the min to the max\n",
    "X_prime_1 = pd.DataFrame({'GNP': np.linspace(X.GNP.min(), X.GNP.max(), 100)})\n",
    "X_prime_1 = sm.add_constant(X_prime_1)  # add constant as we did before\n",
    "\n",
    "y_hat_int = est.predict(X_prime_1)\n",
    "y_hat_no_int = est_no_int.predict(X_prime_1)\n",
    "\n",
    "fig = plt.figure(figsize=(8,4))\n",
    "splt = plt.subplot(121)\n",
    "\n",
    "splt.scatter(X.GNP, y, alpha=0.3)  # Plot the raw data\n",
    "plt.ylim(30, 100)  # Set the y-axis to be the same\n",
    "plt.xlabel(\"Gross National Product\")\n",
    "plt.ylabel(\"Total Employment\")\n",
    "plt.title(\"With intercept\")\n",
    "splt.plot(X_prime[:, 1], y_hat_int, 'r', alpha=0.9)  # Add the regression line, colored in red\n",
    "\n",
    "splt = plt.subplot(122)\n",
    "splt.scatter(X.GNP, y, alpha=0.3)  # Plot the raw data\n",
    "plt.xlabel(\"Gross National Product\")\n",
    "plt.title(\"Without intercept\")\n",
    "splt.plot(X_prime[:, 1], y_hat_no_int, 'r', alpha=0.9)  # Add the regression line, colored in red"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But notice that this may not be the best idea... :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlation and Causation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly there is a relationship or correlation between GNP and total employment. So does that mean a change in GNP cause a change in total employment? Or does a change in total employment cause a change in GNP? This is a subject we will explore in the next post."
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}