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 2017/18</h1></center>\n",
    "\n",
    "<center><h1>Practical 9: Correlation and Regression</h1></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. 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 `panda`s 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 sb\n",
    "import matplotlib.pyplot as plt    #Plotting library used by seaborn, see http://matplotlib.org/users/pyplot_tutorial.html\n",
    "%matplotlib inline  \n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also make use of `numpy` for some functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "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) # The 'low memory' option means pandas doesn't guess data types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Covariance and Correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 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",
    "\n",
    "print \"Covariance matrix:\", '\\n', covmat, '\\n'"
   ]
  },
  {
   "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:** To ensure you understand the structure of the covariance matrix just created, add values (to 2 decimal places) in the code cells below to provide values for:\n",
    "\n",
    "- The covariance between `HHOLDRES` and `Owned` is ...\n",
    "- The covariance between `PM10max` and `SocialRented` is ...\n",
    "\n",
    "(_Ask_ if you are not sure how the matrix is structured)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The covariance between HHOLDRES and Owned is "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The covariance between PM10max and SocialRented is "
   ]
  },
  {
   "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",
    "print \"Pearson correlation coefficient matrix:\", '\\n', corrmat, '\\n'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix produced has the same structure as for the covariance matrix. \n",
    "\n",
    "**Task:** 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 ...\n",
    "- The Pearson correlation between `PM10max` and `SocialRented` is ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The Pearson correlation between HHOLDRES and Owned is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The Pearson correlation between PM10max and SocialRented is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** Which pair of variables has the strongest relationship? \n",
    "\n",
    "_Your answer here_ [Edit this cell]"
   ]
  },
  {
   "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": [
    "sb.heatmap(corrmat)\n",
    "plt.title(\"Pearson Correlation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** Compare the plot just created to the Pearson correlation matrix to check you understand how the heatmap plot represents the matrix. "
   ]
  },
  {
   "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.heatmap(corrmat)\n",
    "plt.title(\"Pearson Correlation\")\n",
    "plt.savefig('Heatmap-Pearson.png', bbox_inches='tight')\n",
    "plt.close()"
   ]
  },
  {
   "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[:,'GreenspaceArea':'PrivateRented']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now run the pairplot for the subset DataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = sb.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!\n",
    "\n",
    "**Task:** 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. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can note two things from the pairplot above:\n",
    "1. Several of the variable are not normally distributed\n",
    "2. The relationships between these variables are not particularly clear \n",
    "\n",
    "Let's see if these observations hold for some of the air quality variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NO = df.loc[:,'NOxmean':]\n",
    "fig = sb.pairplot(NO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** 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 string relationships and which are weak. Think about possible reasons for any differences. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#remove!\n",
    "_Answer: Generally the strongest relationships are between the same measures of NOx and NO2 (e.g. mean vs mean and sd vs sd)\n",
    "Mean and Min values are most normal-like, sd and max are generally non-normal_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jointplot\n",
    "\n",
    "Above we have seen how to calculate correlation matrices and plots for _all_ series (variuables) in a DataFrame. But what if we wanted to focus on on specific pairs of variables? To do this we can use the seaborn `jointplot`.\n",
    "\n",
    "For example, above we looked at relationships between `HHOLDRES` and `Owned`, and between `PM10max` and `SocialRented`. Let's create `jointplot`s for these pairs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.jointplot(x=\"HHOLDRES\", y=\"Owned\", data=df) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.jointplot(x=\"PM10max\", y=\"SocialRented\", data=df) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that by default the Pearson correlation coefficient is presented; you can prevent this using `stat_func` argument (`stat_func = None`) or to specify a different function (we'll see an example of this below). "
   ]
  },
  {
   "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.corr(method = \"spearman\")\n",
    "print \"Spearman rank correlation coefficient matrix:\", '\\n', corspmat, '\\n'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** 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": [
    "#The Spearman correlation between HHOLDRES and Owned is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The Spearman correlation between PM10max and SocialRented is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** 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_ [Edit this cell]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer:\n",
    "HHOLDRES and Owned: weak relationship so distribution/linearity not much of an issue and Peason is fine\n",
    "PM10max and SocialRented: non-normal distribution of variables means spearman more appropriate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** For the pairplots above, write in the cell below to write images to file on your hard disk as `.png` files. "
   ]
  },
  {
   "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 analysis you did in Week 7. If you look back to Week 7, you'll remember that you looked at the possible relationship between pollution and the presence of major roads: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sb.jointplot(x=df.RoadsArea, y=df.NOxmax, size=6)\n",
    "g.fig.suptitle('NOx (Max) against Roads Area\"', fontsize=18,color=\"r\",alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From Week 7:\n",
    "\n",
    ">It looks like our RoadsArea data is very heavily skewed and that there are a lot of very low values in the data. There also seems to a be at least one major outlier that is _so_ different in scale that I'd want to know whether to even keep it in the analysis -- it looks like what would be called a 'leverage point' in a regression model: something that is so 'out of whack' that it alters the entire regression! \n",
    "\n",
    ">Let's try stripping that out and running this code again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The outlier LSOAs are: \" + \",\".join(df[df.NOxmax > 4000]['LSOA11NM'].values))\n",
    "df = df[df.NOxmax < 4000]\n",
    "\n",
    "g = sb.jointplot(x=df.RoadsArea, y=df.NOxmax, size=6)\n",
    "g.fig.suptitle('NOx (Max) against Roads Area\"', fontsize=18, color=\"r\",alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Week 7, you then went on to look at how a transformation might be useful to address the skew of the data. Before we do that again this week to see why it's useful for regression, let's see how the Spearman rank correlation would look for these variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import spearmanr\n",
    "sb.jointplot(x=\"RoadsArea\", y=\"NOxmax\", data=df, stat_func=spearmanr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in the code above we imported the `spearmanr` funtction from the `scipy.stats` [package](https://docs.scipy.org/doc/scipy/reference/stats.html), then we passed that to the `jointplot`.\n",
    "\n",
    "Also, note how the Spearman rank correlation coefficient is greater than the Pearson correlation coefficient (in the previous plot). Think about how ranking the data might improve the correlation, compared to looking at the absolute values themselves. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Right, back to Week 7...\n",
    "\n",
    ">To create a new series in the data frame containing the natural log of the original value it’s a similar process to what we've done before, but since pandas doesn't provide a log-transform operator (i.e. you can’t call `df.Owned.log()` ) we need to use the `numpy` package:\n",
    "```python\n",
    "series = pd.Series(np.log(df.RoadsArea))\n",
    "```\n",
    ">Try performing the transformation and then `describe()` the results in the coding area below. Is it more clear to you now why a log-transform is a non-linear transformation?\n",
    "\n",
    "Here's the code for you:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option('display.float_format', lambda x: '{:,.4f}'.format(x))\n",
    "\n",
    "series = pd.Series(np.log(df.RoadsArea + 1))\n",
    "print(series.describe())\n",
    "sb.distplot(series, color='green')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Week 7:\n",
    ">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": [
    "From this point on we'll leave Week 7 behind, but 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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hrds = df.loc[df.RoadsArea > 0] # hrds == has (major) roads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many LSOAs do we have left? \n",
    "\n",
    "**Task:** 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": [
    "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` variable 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 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 (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": [
    "print(NOxmax_roads_mod_fit.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** 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!) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add your answers here _Click to edit this cell_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But now 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 fir OLS model object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.hist(NOxmax_roads_mod_fit.resid)\n",
    "plt.xlabel('Residuals')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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",
    "More important (given we fit the model with >3000 data points), is to check the distribution of the standardised residuals (gain, look back to your lecture notes). "
   ]
  },
  {
   "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, 'bo')\n",
    "plt.axhline(linestyle = 'dashed', c = 'black')\n",
    "plt.xlabel('Predicted y Values')\n",
    "plt.ylabel('Standardized Residuals (z)')                "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing this plot to what we saw in lecture, it looks like we don't have 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": [
    "You likely got a warning there, so let's just quickly check that we created the new `series` correctly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hrds.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like we're okay, so we'll ignore the warning and 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.from_formula('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.\n",
    "\n",
    "**Task:** 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": [
    "#a histogram of residuals\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#a scatterplot of standardised residuals vs predicted values\n",
    "       "
   ]
  },
  {
   "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": [
    "print logNOxmax_logRoads_mod_fit.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** 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? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add your answers here _Click to edit cell_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should have noted that the r2 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 = sb.regplot(x=hrds.logRoadsArea, y=hrds.logNOxmax)  \n",
    "plt.xlabel('log(RoadsArea)')\n",
    "plt.ylabel('log(NOx max)')    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression by Brorough\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 boroughs level. To do that we'll need to use some of the techniques we learned for grouping data in Week 6. In particular, we use the `groupby` pandas  method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grouped = hrds.groupby('LAD11NM')\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",
    "sb.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:** 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:** 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 not normally distributed), so let's try a regression model with a log transform of this variable.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** 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": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** 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...\n",
    "\n",
    "Sticking with this model, we now need to check standardised residuals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** 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 goog so let's go back and interpret the summary of the model\n",
    "\n",
    "**Task:** 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?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add your answers here _Click to edit cell_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Looping through borough\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: 'r', 1: 'sp', 2: 'n'})\n",
    "\n",
    "print mySummary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task:** 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? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add answers here _(double-click to edit)_"
   ]
  },
  {
   "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 through map. Let's see how we can combine the results so far with code for mapping from Week 8 (look back to the previous notebook if you need a reminder). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, import the relevant packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pysal as ps\n",
    "import geopandas as gpd\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can load the spatial data. Remember these data were donwloaded using last week using the `download_gsa_geodata` function we defined - you may need to copy those data to an appropriate location to match the code below, or change the code so that the path points to the location of the data on your computer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shp_path = os.path.join('shapes','lsoas','Lower_Layer_Super_Output_Areas_December_2011_Generalised_Clipped__Boundaries_in_England_and_Wales.shp')\n",
    "print(\"Loading data from: \" + shp_path)\n",
    "\n",
    "gdf = gpd.read_file(shp_path)\n",
    "gdf.set_index('lsoa11cd', drop=True, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check what this looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print gdf.sample(3)\n",
    "print df.sample(3)\n",
    "print df.index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Currently, the data does not include a series (column) that provides useful borough-level labels. But we have those in the original DataFrame we loaded at the start of the notebook, so let's merge the two:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df.set_index('LSOA11CD', drop=True, inplace=True)  #may need to set the index for df if you get an error (if uncomment this line)\n",
    "\n",
    "gdf2 = pd.merge(gdf, df, left_index=True, right_index=True)\n",
    "gdf2.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `gdf2` DataFrame should now have 63 columns, including the series `LAD11NM`. If the following code gives `False` you might need to check the `merge` in the last code cell worked properly:   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('LAD11NM' in gdf2.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming the `LAD11NM` column exists in the `gdf2` DataFrame, we now need to do _another_ `merge` to join the correlation data for each borough (held in the mySummary DataFrame we created in the loop). We do this using the `LAD11NM` column we just included in the merge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cordf = pd.merge(gdf2, mySummary, left_on = 'LAD11NM', right_index = True)   #causes geometry column to become geometry_x!! \n",
    "cordf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your `cordf` DataFrame should have 66 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 [NB, you will likely get an error here]:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cordf.plot(column= 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did you get an error? Something like:\n",
    "\n",
    "`AttributeError: No geometry data set yet (expected in column 'geometry'.`\n",
    "\n",
    "See if you can work out why we get this error by looking at the column names of the cordf DataFame and comparing them to the error...\n",
    "\n",
    "<br></br>\n",
    "<br></br>\n",
    "<br></br>\n",
    "<br></br>\n",
    "<br></br>\n",
    "<br></br>\n",
    "<br></br>\n",
    "<br></br>\n",
    "<br></br>\n",
    "<br></br>\n",
    "\n",
    "You may have spotted that somewhere in all that merging of DataFrames, the column containing the geometry information of the polygons has been renamed. Originally it was called `geometry`, but it seems to have ended up now as `geometry_x`. Let's use the `rename` method to get it back to its orginal state: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cordf = cordf.rename(columns={'geometry_x': 'geometry'})\n",
    "cordf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Now_ our plotting of the map should work..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cordf.plot(column= 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you've finally managed to get it to all come together, your map should look like this:\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": [
    "## Exercises\n",
    "\n",
    "Finally, here are some exercises to help you to reinforce, and extent upon, what you have learned above and throughout the module. We'll provide suggested answers to these next week. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1**\n",
    "\n",
    "Build on the (looping) code for loopingcalculating correlations for all boroughs to do similar, but instead of calculating correlations, fit regressions (for specified variables in a DF) for ALL boroughs, summarising 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 to guide you. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Excercise 2**\n",
    "\n",
    "You may have noticed above that we repeated quite a lot of code, but with slight variations in 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}