Practical-09-Correlation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "geopyter": {
     "Contributors": "James Millington (james.millington@kcl.ac.uk)",
     "git": {
      "active_branch": "master",
      "author.name": "Jon Reades",
      "authored_date": "2017-08-17 19:06:58",
      "committed_date": "2017-08-17 19:06:58",
      "committer.name": "Jon Reades",
      "sha": "5e3b396ae18a982d693c4bfd86c721c7e1e21051"
     }
    }
   },
   "source": [
    "<center><h1>7SSG2059 Geocomputation</h1></center>\n",
    "\n",
    "<center><h1>Practical 9: Correlation and Regression</h1></center>\n",
    "\n",
    "<center><h2>Relationships and Explanations</h2></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Correlation and Regression\n",
    "\n",
    "Correlation and regression are useful tools to understand relationships between variables in our data, and begin to explain possible causes. This week, we will look at some possible ways that might use these tools to analyse the data for your final report. \n",
    "\n",
    "Specifically we will look at:\n",
    "1. correlation and regression for ALL LSOAs\n",
    "2. correlation and regression for LSOAs grouped by borough\n",
    "3. mapping correlations by borough"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup\n",
    "\n",
    "As usual we will be using pandas and geopandas for data analysis, with `matplotlib` and `seaborn` for visualisation. Let's load those now: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt    #Plotting library used by seaborn, see http://matplotlib.org/users/pyplot_tutorial.html\n",
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import matplotlib as mpl\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "import statsmodels.api as sm "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And ignore any warnings for now..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings \n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Later we will import further packages for regression and mapping. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, in our setup, let's load the data into a pandas `DataFrame`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\n",
    "    'https://github.com/kingsgeocomp/geocomputation/blob/master/data/LondonLSOAData.csv.gz?raw=true',\n",
    "    compression='gzip', low_memory=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Covariance and Correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correlation & Covariance Matrices\n",
    "\n",
    "Let's start by seeing how we can calculate covariance using pandas. Remember, covariance is like an unstandardised version of correlation. Handily, pandas has [a method](http://pandas.pydata.org/pandas-docs/stable/computation.html#covariance) to calculate a matrix (like a table) of covariance values between `series` in a `DataFrame`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "covmat = df.cov()\n",
    "covmat.iloc[1:10,1:10] # Notice what integer indexing ('*i*nteger *loc*ation') gives us?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a lot of numbers... what do they all mean? \n",
    "\n",
    "Usually, I'd suggest you pause and have a think about what these numbers tell us. But as discussed in lecture, **correlation coefficients** are often more useful for comparing and understanding relationships than **covariance values**. The important thing here is to understand the _structure_ of the covariance matrix produced above, before we move on to correlation.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 1.1\n",
    "\n",
    "To ensure you understand the structure of the covariance matrix just created, retrieve values (and format to 2 decimal places) in the code cells below to provide values for:\n",
    "\n",
    "- The covariance between `HHOLDRES` and `Owned` is ... (1784.69)\n",
    "- The covariance between `PM10max` and `SocialRented` is ... (124.98)\n",
    "\n",
    "(_Ask_ if you are not sure how the matrix is structured)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to the covariance method, pandas has [a method](http://pandas.pydata.org/pandas-docs/stable/computation.html#correlation) for calculating the correlation between all `series` in a `DataFrame`. If we don't specify what particular correlation we want, the `corr` method calculates Pearson's _r_ correlation coefficient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corrmat = df.corr()\n",
    "corrmat.iloc[1:10,1:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix produced has the same structure as for the covariance matrix. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 1.2 \n",
    "\n",
    "Let's identify the Pearson correlation coefficient (i.e. value) for the same pairs of variables as we did for covariance (edit the cells below again, providing values to three decimal places): \n",
    "\n",
    "- The Pearson correlation between `HHOLDRES` and `Owned` is ... (0.051)\n",
    "- The Pearson correlation between `PM10max` and `SocialRented` is ... (0.144)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TASK:** Compare the Pearson _r_ values to the covariance values. Check you understand why the values are different and why correlation values are often more useful. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correlation Heatplot\n",
    "\n",
    "Even though the 'standardised' correlation values are a bit easier to read than the covariance values, it would still be useful to think about how we can visualise these numbers for quick reference. The seaborn `heatmap` [plot](http://seaborn.pydata.org/generated/seaborn.heatmap.html) is useful in this circumstance:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1,figsize=(13,11.5))\n",
    "sns.heatmap(corrmat, cmap='Spectral')\n",
    "plt.title(\"Pearson Correlation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 2.1\n",
    "\n",
    "Compare the plot just created to the Pearson correlation matrix to check you understand how the heatmap plot represents the matrix. For example, why is there a diagonal line of high values? Why does the bottom right corner have the pattern it does?\n",
    "\n",
    "> Your answer here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 2.2 \n",
    "\n",
    "Why do you think we used a 'Spectral' colormap here?\n",
    "\n",
    "> Your answer here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember you could save this plot to file (in the same directory where your notebook is saved) for use in your reports by doing something like:\n",
    "```python\n",
    "plt.savefig('Heatmap-Pearson.png', bbox_inches='tight', dpi=200)\n",
    "plt.close()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pairplot\n",
    "\n",
    "Another plot that is often useful is the seaborn `pairplot` which produces scatter plots for all pairs of 'series' in a `DataFrame`. Read more [here](http://seaborn.pydata.org/tutorial/distributions.html#visualizing-pairwise-relationships-in-a-dataset). \n",
    "\n",
    "Our LSOA dataframe is quite large and you would find, if you tried to run the `pairplot` function for the entire dataframe, that it would take a very long time to calculate and the plot produced would be huge. So, to use the pairplot function here we will subset the dataframe to fewer `series` and run the pairplot function on that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's remind ourselves of the series in our dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's produce the pairplot for all series between `GreenspaceArea` and `PrivateRented`. So first, create the subset DataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sub = df.loc[???]  # subset all rows, and columns from GSA to PR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now run the pairplot for the subset DataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = sns.pairplot(sub)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It may have taken a few seconds to run that function - it would have been much longer had we tried to do it for the entire DataFrame! (which is why I suggested we take a subset of the DataFrame)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 3.1\n",
    "\n",
    "Take a look at the pairplot produced and check you can see how it is a set of scatter plots for each pair of variables with a histogram for each individual variable. Pretty nice.\n",
    "\n",
    "> Any comments or observations here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 3.2\n",
    "\n",
    "Let's see if these observations hold for some of the air quality variables, but to help you learn some Python syntax I want you to use a 'list comprehension' to select only the NO variables from the data frame while _also_ excluding the NO standard deviation variables. The basic format for a list comprehension is:\n",
    "```python\n",
    "[x for x in <list>]\n",
    "```\n",
    "The output of a list comprehension is _also_ a list. This is basically a really clever way to write 'do something to every element in the list without needing a for loop'. So, for instance:\n",
    "```python\n",
    "[2*x for x in <list>]\n",
    "```\n",
    "Would double every element of the list and return _that_ as a new list.\n",
    "\n",
    "The real magic is that we can _also_ filter data:\n",
    "```python\n",
    "[x for x in <list> if x % 2==0]\n",
    "```\n",
    "would return only the even elements of a list.\n",
    "\n",
    "With those as clues, think about what the `NO...` columns have in column and _then_ think about how you would remove the NO columns that are measuring the standard deviation. This can all be done in one line, but you might want to start out by trying to do it in two list comprehensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NOx = df.loc[???]\n",
    "fig = sns.pairplot(NOx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 3.3 \n",
    "\n",
    "Looking at the pairplot for NO variables check you can see which variables have normal-like distributions and which are non-normal, and which variables have strong relationships and which are weak. Think about possible reasons for any differences. Also look closely at the scatter plots themselves; do you notice anything that might be influencing your impressions?\n",
    "\n",
    "> Your answer here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jointplot\n",
    "\n",
    "Above we have seen how to calculate correlation matrices and plots for _all_ series (variables) in a DataFrame. But to zoom in on individual relationships, such as between `HHOLDRES` and `Owned` or between `PM10max` and `SocialRented`, let's create `jointplot`s for these pairs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.jointplot(x=\"HHOLDRES\", y=\"Owned\", data=df, marker='.', color='red', joint_kws={'s':2})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.jointplot(x=\"PM10max\", y=\"SocialRented\", data=df, marker='.', joint_kws={'s':2})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the jointplots above, do you think a Pearson correlation is appropriate? Let's check the Spearman rank correlation coefficients for our data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corspmat = df.???(method=???)\n",
    "print(\"Spearman rank correlation coefficient matrix:\")\n",
    "corspmat.iloc[1:10,1:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 4.1 \n",
    "\n",
    "Identify the Spearman correlation coefficient (i.e. value) for the same pairs of variables as we did for Pearson's r above (edit the cells below again, providing values to three decimal places): \n",
    "\n",
    "- The Spearman correlation between `HHOLDRES` and `Owned` is ...\n",
    "- The Spearman correlation between `PM10max` and `SocialRented` is ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"{:0.3f}\".format(corspmat.loc['HHOLDRES'][['Owned']].values[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"{:0.3f}\".format(corspmat.loc['PM10max'][['SocialRented']].values[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 4.2 \n",
    "\n",
    "Compare the Spearman rank correlation coefficients you have just entered above, to the corresponding Pearson correlation coefficients. For each pair of variables, which correlation do you think is most appropriate?\n",
    "\n",
    "> Your answer here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 4.3 \n",
    "\n",
    "For the pairplots above, enter code in the cells below to write (i.e. save) images to file on your hard disk as `.png` files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression\n",
    "\n",
    "Now that we have seen how you might use correlation to examine the data, let's move on to look at regression by picking up on some of the ideas we examined in in the notebook on _transformations_.\n",
    "\n",
    "We'll look at the possible relationship between pollution and the presence of major roads: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1,figsize=(8,8))\n",
    "g = sns.regplot(x=df.RoadsArea, y=df.NOxmax)\n",
    "fig.suptitle('NOx (Max) against Roads Area\"', fontsize=18, color=\"k\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like our _RoadsArea_ data are very heavily skewed: there are quite a lot of very low values in the data, but there also seems to a be at least one major _NOxmax_ outlier that is _so_ different in scale that we should consider whether to even keep it in the analysis. This is what we would call a 'leverage point' in a regression model: it alters the entire regression! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dealing with Leverage\n",
    "\n",
    "Let's identify which LSOA has such a high _NOxmax_ value (using a new bit of code, you can investigate what it's doing later if you like):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The outlier LSOAs are: \" + \",\".join(df[df.NOxmax > 2000]['LSOA11NM'].values))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 5.1\n",
    "\n",
    "Where is Hillingdon? What mught be happening there that means it has such high NOx values?\n",
    "\n",
    "> Your answer here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try stripping that LSOA out and make the plot again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first remove the LSOA with high NOx\n",
    "dfn = df[???]\n",
    "\n",
    "# then plot\n",
    "fig, ax = plt.subplots(1,1,figsize=(8,8))\n",
    "g = sns.regplot(x=dfn.???, y=dfn.???)\n",
    "fig.suptitle('NOx (Max) against Roads Area\"', fontsize=18, color=\"k\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distribution of these variables is _still_ quite skewed; let's see how the Spearman rank correlation (which is more robust to skewed distributions and Pearson's r) would look for these variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import spearmanr  #import spearmanr function from scipy.stats \n",
    "from scipy.stats import pearsonr  #import pearson function from scipy.stats \n",
    "\n",
    "??? = spearmanr(dfn.RoadsArea.values, dfn.NOxmax.values)\n",
    "??? = pearsonr(dfn.RoadsArea.values, dfn.NOxmax.values)\n",
    "\n",
    "print(\"Spearman Correlation results: {0:0.5f} correlation with p-value {0:0.5f}\".format(???, ???))\n",
    "print(\"Pearson Correlation results: {0:0.5f} correlation with p-value {0:0.5f}\".format(???, ???))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in the code above we imported the `spearmanr` and `pearsonr` from the `scipy.stats` [package](https://docs.scipy.org/doc/scipy/reference/stats.html). This is what `pandas` is doing in the background (because, again, it's working with numpy arrays and scipy is designed to work with this as well).\n",
    "\n",
    "Note how the Spearman rank correlation coefficient is greater than the Pearson correlation coefficient, but the p-value is as well indicating a less significant result. Think about how ranking the data (which is what Spearman implicitly does) might improve the correlation, compared to looking at the absolute values. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to transform the data so that they look more like a normal distribution.Let's try a log transform:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logRA = pd.Series(np.log(???))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's examine these trandformed data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(logRA.describe().round(3))\n",
    "sns.distplot(logRA, color='green')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now this _is_ interesting: the output of the graph shows what seems to be two quite different things going on in our data! We've obviously got the LSOAs that contain _no_ major roads, but then we've got something else that is _much_ closer to 'normal' (though obviously not properly normal as there is clear evidence of negative skew). Technically, this is closer to _log-normal_ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we go on to see if we can use regression to estimate the importance of roads for NOx we'll see how transforms are useful for regression. \n",
    "\n",
    "First, though let's remove the LSOAs from our data that have no major roads and plot the distribution of the remaining data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hrds = dfn.loc[dfn.RoadsArea > 0].copy() # hrds == has (major) roads\n",
    "sns.distplot(hrds.RoadsArea)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is this distribution normal? How would you describe it?\n",
    "\n",
    "How many LSOAs do we have left? \n",
    "\n",
    "#### Task 6.1 \n",
    "\n",
    "Add a line of code here to check how many LSOAs are in the `hrds` dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performing a Regression in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thinking back to what Lumley _et al._ (2002, p.166) said about how, \n",
    "\n",
    ">\"…linear regression [does] not require any assumption of Normal distribution in sufficiently large samples. Previous simulations studies show that “sufficiently large” is often under 100, and even for our extremely non-Normal medical cost data it is less than 500.\"\n",
    "\n",
    "Hopefully, it's clear that in the `hrds` data we maybe don't need to worry about the fact that the `RoadsArea` data (where roads are present) has a log-normal distribution. So let's just fit a regression with our un-transformed data and see what we get.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In python, linear regression can be performed using functions available in the `statsmodels` [package](http://www.statsmodels.org/stable/), and specifically using the [OLS function](http://www.statsmodels.org/devel/examples/notebooks/generated/ols.html) from the `statsmodels.api`. So let's import `statsmodels.api` first:     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the `OLS` function from `statsmodels` to fit a regression requires we create an `OLS` object first, then use the `fit` method on that object. To create the `OLS` object we can use the `from_formula` method to pass the equation of the model we want to fit (as well as indicating what the data are that we are using):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NOxmax_roads_mod = sm.OLS.from_formula(\"NOxmax ~ RoadsArea\", data = hrds) \n",
    "NOxmax_roads_mod_fit = NOxmax_roads_mod.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happened...? Well, it looks like nothing happened but we have indeed now fit a regression model! \n",
    "\n",
    "To check this we we should look at a summary of the model (by using the `summary` method on the `fit model` object we just created) to understand what the model can tell us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NOxmax_roads_mod_fit.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 5.2\n",
    "\n",
    "Interpret the results of your regression by answering the following questions:\n",
    "- What is r-squared value?\n",
    "- What is p-value for the `RoadsArea` variable?\n",
    "- What is the effect size of the `RoadsArea` variable?  (you might need to look at LSOA metadata to check the units of the variables we are modelling!)\n",
    "\n",
    "> Your answers here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But now we need to check for problems in our residuals (look back to your lecture notes about this). \n",
    "\n",
    "First, let's plot a histogram of the residuals, using the `.resid` method applied to the fit OLS model object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.hist(???.???)\n",
    "plt.xlabel('Residuals')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 5.3\n",
    "\n",
    "What do you think? Are the residuals normally distributed? But, given the number of values we used to fit the model does this question matter?\n",
    "\n",
    "> Your answer here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's no built-in method for calculating standardised residuals, so we do that first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate standardized residuals ourselves\n",
    "NOxmax_roads_mod_sr = (NOxmax_roads_mod_fit.resid / np.std(NOxmax_roads_mod_fit.resid)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plot these in a scatter plot against the predicted values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.plot(NOxmax_roads_mod_fit.fittedvalues, NOxmax_roads_mod_sr, 'r.')\n",
    "plt.axhline(linestyle = 'dashed', c = 'black')\n",
    "plt.xlabel('Predicted y Values')\n",
    "plt.ylabel('Standardized Residuals (z)')                "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These certainly don't _look_ like constant variance in the standardised residuals.\n",
    "\n",
    "Thinking back to what Lumley _et al._ (2002, p.166) said about how, \n",
    "\n",
    "> \"Linear regression does assume that the variance of the outcome variable is approximately constant\"\n",
    "\n",
    "So despite having such a large data set, we _do_ need to sort this out. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So let's try fitting the model with _transformed_ data. We'll try a log transform, so first we need to calculate new values for our variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hrds['logRoadsArea'] = np.log(hrds['RoadsArea'])\n",
    "hrds['logNOxmax']    = np.log(hrds['NOxmax'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's just quickly check that we created the new `Series` correctly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hrds.iloc[1:5,-5:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like we're okay, so we'll move on to now fit the regression with our new log transformed variables: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logNOxmax_logRoads_mod     = sm.OLS.???('logNOxmax ~ logRoadsArea', data=hrds)\n",
    "logNOxmax_logRoads_mod_fit = logNOxmax_logRoads_mod.fit() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before looking at the summary, let's look at the residuals to see if we have overcome the issue we had with the un-transformed data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 5.4 \n",
    "\n",
    "Add code in the two code blocks below to plot:\n",
    "- a histogram of residuals\n",
    "- a scatterplot of standardised residuals vs predicted values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should see that the histogram of the residuals looks quite normal, and the variance of the standardised resiudals is reasonably constant. So using the log transformed variables seems to have helped! \n",
    "\n",
    "So now let's look at the summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logNOxmax_logRoads_mod_fit.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 5.5 \n",
    "\n",
    "Interpret the results of your regression by answering the following questions:\n",
    "- What is the r-squared value?\n",
    "- What is p-value for the `logRoadsArea` variable?\n",
    "- What is the effect size of the `logRoadsArea` variable?  (remember we now have log values! So see Table 2 of Lin et al. in lecture slides) \n",
    "- What is the confidence interval?\n",
    "\n",
    "> Your answers here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should have noted that the $r^{2}$ value is not very large, indicating not much of the variation in NOx is explained by variation in RoadsArea (and this is also shown by the small effect size). \n",
    "\n",
    "A quick look at a scatter plot of the two variables shows why the _r2_ is so poor: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = sns.regplot(x=hrds.logRoadsArea, y=hrds.logNOxmax)  \n",
    "plt.xlabel('log(Roads Area)')\n",
    "plt.ylabel('log(NOx max)')    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's worth noting that regression plots (while very useful) _can_ be hard to read accurately because the number of data points in the main cluster is masked by the size of the points themselves. For instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1,figsize=(14,8))\n",
    "h = sns.regplot(x=hrds.logRoadsArea, y=hrds.logNOxmax, marker='.', ax=ax)  \n",
    "plt.xlabel('log(Roads Area)')\n",
    "plt.ylabel('log(NOx max)')    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Residuals Geographically\n",
    "\n",
    "Here is something that is easy to do, but all too rarely done by researchers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(logNOxmax_logRoads_mod_fit.resid)\n",
    "plt.xlabel('Residuals')\n",
    "\n",
    "print(logNOxmax_logRoads_mod_fit.resid.sample(n=15)) # Just so that you can see this is a data series!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Join the residuals back on to the HRDS data frame and \n",
    "# then join that to the geometry of LSOAs so that we can\n",
    "# plot them.\n",
    "hrds['Residuals'] = logNOxmax_logRoads_mod_fit.resid\n",
    "lsoas = gpd.read_file('https://github.com/kingsgeocomp/geocomputation/raw/master/data/src/LSOAs.gpkg', driver='GPKG')\n",
    "geom = lsoas.geometry\n",
    "p = pd.merge(lsoas, hrds, left_on='lsoa11cd', right_on='LSOA11CD', how='inner')\n",
    "gdf = gpd.GeoDataFrame(p, geometry=geom)\n",
    "\n",
    "# And now plot it.\n",
    "fig, ax = plt.subplots(1,1,figsize=(15,12))\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "gdf.plot(column='Residuals', ax=ax, cmap=mpl.cm.get_cmap('viridis'), linewidth=0, legend=True, cax=cax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we're looking for is pattern in the residuals because that implies that our model is _systematically_ under- or over-predicting errors. We'd normally interpret this in the context of missing variable(s) from the model. There is a _slight_ suggestion of some patterning to the data but, if so, it's very mild. So our model isn't missing something _obvious_ but, that said, it's also the case that we are still doing a pretty bad job of 'explaining' the data with this model (ca. 6%)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression by Borough\n",
    "\n",
    "One reason there may be a poor relationship between roads and NOx for our entire data set is that there may be variation in that relationship across London. For example, would we expect the same influence of roads on NOx in central London compared to the suburbs?\n",
    "\n",
    "So let's have a look at how we could calculate correlations or fit regressions for data **at the borough-level**. To do that we'll need to use some of the techniques we learned for grouping data in Week 7. In particular, we use the `groupby` pandas  method again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grouped = hrds.groupby(???)\n",
    "print(grouped.groups.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now create DataFrames for individual boroughs by using `get_group` methods on the `groupby` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Hackney = grouped.get_group('Hackney')\n",
    "sns.jointplot(x='NOxmax', y='RoadsArea', data=Hackney) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this plot we can see there might indeed be a stronger relationship for the borough of Hackney. But also we notice that there are far fewer data points (as now we are looking at only the LSOAs in one borough, not across the whole of London).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 6.1 \n",
    "\n",
    "Given that it looks like there may be a stronger relationship for the borough of Hackney, fit a linear regression using the `Hackney` DataFrame for `NOxmax ~ RoadsArea`. Refer to and reuse code from above if you need to:\n",
    "- first create a model object\n",
    "- then use the `.fit()` method to actually fit the regression\n",
    "- finally, use the `.summary()` method to interpret the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should be able to see from your output that the relationship for Hackney has a much larger r-squared value. \n",
    "\n",
    "But the model was fit using only 104 data points (compared to >3000 for the whole of London) so now we should be a little more worried about the distribution of the residuals. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 6.2 \n",
    "\n",
    "Plot a histogram to check the distribution of residuals for your Hackney regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should see from your histogram that the distribution is not very normal, and actually looks more like a log-normal distribution. From last `jointplot` for the Hackney data, it looks like _NOxmax_ is the problem (in the sense that this variable is farther from a normal distribution than _RoadsArea_), so let's try a regression model with a log transform of this variable.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 6.3 \n",
    "Fit a regression model for `logNOxmax ~ RoadsArea` for the Hackney data and print the summary to allow interpretation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 6.4 \n",
    "Plot a histogram to check the distribution of the residuals you have just fit: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From your histogram you should be able to see that the residuals for this model are much more normal, and we can probably be happy with this (given we still have ~100 data points). But note that in other cases we might use the log transform for _both_ variables (i.e. both _RoadsArea_ and _NOxmax_).\n",
    "\n",
    "Sticking with this model, we now need to check the standardised residuals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 6.5 \n",
    "Calculate standardised residuals for your latest model and plot them against the predicted values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variance of the standardised residuals looks pretty good so let's go back and interpret the summary of the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 6.6 \n",
    "\n",
    "Interpret the results of your regression by answering the following questions:\n",
    "- What is the r-squared value?\n",
    "- What is p-value for the `RoadsArea` variable?\n",
    "- What is the effect size of the `RoadsArea` variable?  (remember one of the variables is a log transform! So see Table 2 of Lin et al. in lecture slides) \n",
    "- What is the confidence interval?\n",
    "\n",
    "> Your answers here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Looping through boroughs\n",
    "\n",
    "Now we have seen how we can investigate correlations and relationships for a single borough, let's think about how we could automate this for all boroughs using grouping and a loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import pearsonr\n",
    "from scipy.stats import spearmanr\n",
    "\n",
    "boroughs = hrds.groupby('LAD11NM')\n",
    "bnames   = grouped.groups.keys()\n",
    "\n",
    "r  = []\n",
    "sp = []\n",
    "n  = []\n",
    "\n",
    "for name in bnames:\n",
    "\n",
    "    borough = boroughs.get_group(name)\n",
    "    \n",
    "    y = np.log(borough['NOxmax'])\n",
    "    X = np.log(borough['RoadsArea'])\n",
    "    \n",
    "    pr  = pearsonr(X,y)\n",
    "    spn = spearmanr(X,y)\n",
    "    \n",
    "    r.append(pr[0])\n",
    "    sp.append(spn[0])\n",
    "    n.append(len(y))\n",
    "\n",
    "    \n",
    "rSummary  = pd.Series(r, index=bnames)\n",
    "spSummary = pd.Series(sp, index=bnames)\n",
    "nSummary  = pd.Series(n, index=bnames)\n",
    "    \n",
    "mySummary = pd.concat([rSummary, spSummary, nSummary], axis=1)\n",
    "mySummary = mySummary.rename(columns={0: \"Pearson r\", 1: 'Spearman rho', 2: 'n'})\n",
    "\n",
    "mySummary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 7.1 \n",
    "\n",
    "Check you understand the code above by adding comments to it to explain its function, and by answering the following questions:\n",
    "- What package do the `pearsonr` and `spearmanr` functions come from?\n",
    "- Why do we need to use `pr[0]` and not just `pr` to access the Pearson correlation coefficient?\n",
    "- Why do we need `axis = ` for the pandas `append` method?\n",
    "\n",
    "> Your answers here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 7.2 \n",
    "\n",
    "It might even be worth seeing if there are any boroughs for which $r$ and $\\rho$ are wildly different as this might suggest something worth investigating further... so how would you compare the two correlation values in a single figure?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mapping Borough-level Relationships\n",
    "\n",
    "Given that we have now calculated borough-level correlations for _all_ boroughs, maybe it would be nice to visualise this using a map. Let's see how we can combine the results from this week so far with code for mapping from previous weeks (look back to the previous notebook if you need a reminder). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we have spatial information in `boros` and we have the results from all the regressions for each borough in our `mySummary` DataFrame (but without spatial information). So to plot the regression results spatially, we need to `merge` the two DataFrames together into a single DataFrame. \n",
    "\n",
    "But what is there a common column we can use?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(mySummary.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boros = gpd.read_file('https://github.com/kingsgeocomp/geocomputation/raw/master/data/src/Boroughs.gpkg', driver='GPKG')\n",
    "print(boros.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Doesn't look like there's a common column... So what are we going to do?\n",
    "\n",
    "The solution here, lies in the fact that the `mySummary` has an `index` (see documentation [here](http://pandas.pydata.org/pandas-docs/stable/indexing.html), and discussion [here](https://stackoverflow.com/questions/27238066/what-is-the-point-of-indexing-in-pandas) ) that matches one of the columns in `df`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(mySummary.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where did this come from? \n",
    "\n",
    "Well, you may have noticed above that when we made the `mySummary` DataFrame in the loop above, we set the `index` for each `Series` as the value of `bnames`. In turn, `bnames` came from the `keys` of the `GroupBy` object, which in turn was set using the _LAD11NM1_ columns of `df`. That's the long explanation of saying that the `index` of `mySummary` uses the same values as _LAD11NM_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now in our [merge](https://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.merge.html), we can use the `right_index` argument, rather than specfying a column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cordf = pd.merge(boros, mySummary, left_on='NAME', right_index=True, how='inner')\n",
    "print(cordf.head())\n",
    "print(cordf.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your `cordf` DataFrame should have 61 columns - check you can see the columns of `mySummary` on the far right of the `cordf` DataFrame. \n",
    "\n",
    "Also check you understand why the first five rows of the `r`, `sp` and `n` columns all have the same value..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Right, so now we should be able to plot the correlations (the `r` column) as a map... but first we need to convert our `cordf` Pandas DataFrame to a GeoPandas DataFrame, setting the geometry appropriately (look back to Week 8 if you need a reminder of what's going on here): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcordf = gpd.GeoDataFrame(cordf)\n",
    "gcordf = gcordf.set_geometry('geometry')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Now_ our plotting of the map should work..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1,figsize=(15,12))\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "gcordf.plot(column='Pearson r', legend=True, ax=ax, cax=cax)  #include ax argument!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you've managed to get it to all come together, your map should look something like the one below, though here you can see that it's actually the result of a _1:Many_ (_1:n_) join since each LSOA in a given borough has been given the same value whereas our example now only contains boroughs, not LSOAs:\n",
    "\n",
    "![PySAL Logo](https://kingsgeocomputation.files.wordpress.com/2017/11/cordf_r_map.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well done! \n",
    "\n",
    "Have a look at your map and think about what it shows (you may want to add a legend to understand what values the colours refer to). Which borough has the strongest correlation between `NOxmax` and `RoadsArea`? Which has the weakest? Is this the same for the Spearman correlations? etc... "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary\n",
    "\n",
    "In this practical we have looked at several different ways to calculate (correlation matrix) and visualise (heatplot, pairplot) the correlations between many variables in a dataset. We then saw how we could fit regressions in python (using `statsmodels` functions) for particular pairs of variables; this included thinking about the assumptions of regression and whether transforms of the data were needed to meet those assumptions. Then towards the end we looked at how we could combine grouping of the data to calculate and map borough-level correlations. \n",
    "\n",
    "Hopefully you will find many of these techniques useful for analysing the data for your final report. For example, here we focused on on particular pair of variables and one particular borough, but you could look at others for your report.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional Exercises\n",
    "\n",
    "Finally, here are some exercises to help reinforce what you have learned above... and throughout the module. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 8.1\n",
    "\n",
    "Build on the (looping) code for calculating correlations for all boroughs, but instead of calculating correlations you should fit regressions (for specified variables in a df) for ALL boroughs and summarise the results in a table. You may find this [SO question and answer](https://stackoverflow.com/questions/24088439/how-to-apply-ols-from-statsmodels-to-groupby) useful as a guide. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task 8.2\n",
    "\n",
    "You may have noticed above that we repeated quite a lot of code, but with slight variations in variable and object names, when fitting and analysing the regression models. In circumstances like that (when you are repeating code), it can he useful to write yourself a 'helper function' to speed up your analysis. \n",
    "\n",
    "For this exercise, write a helper function to: \n",
    "- read a statsmodels.api OLS model and the data it uses\n",
    "- output (to a .png file) a histogram of the residuals and a plot of standardised residuals against fitted values\n",
    "\n",
    "You should be able to use much of the code from above, but may also need to use some string formatting functions.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Credits!\n",
    "\n",
    "#### Contributors:\n",
    "The following individuals have contributed to these teaching materials: James Millington (james.millington@kcl.ac.uk), Jon Reades (jonathan.reades@kcl.ac.uk)\n",
    "\n",
    "#### License\n",
    "These teaching materials are licensed under a mix of the MIT and CC-BY licenses...\n",
    "\n",
    "#### Acknowledgements:\n",
    "Supported by the [Royal Geographical Society](https://www.rgs.org/HomePage.htm) (with the Institute of British Geographers) with a Ray Y Gildea Jr Award."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "geopyter": {
   "Contributors": [
    "James Millington (james.millington@kcl.ac.uk)"
   ],
   "git": {
    "active_branch": "master",
    "author.name": "Jon Reades",
    "authored_date": "2017-08-17 19:06:58",
    "committed_date": "2017-08-17 19:06:58",
    "committer.name": "Jon Reades",
    "sha": "5e3b396ae18a982d693c4bfd86c721c7e1e21051"
   },
   "libs": {
    "datetime": "?",
    "matplotlib": "1.5.1",
    "os": "?",
    "pandas": "0.20.3",
    "scipy": "0.19.0",
    "seaborn": "0.7.0"
   }
  },
  "kernelspec": {
   "display_name": "GSA2019",
   "language": "python",
   "name": "gsa2019"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}