Practical-10a-Weather-AirQuality.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><h1>7SSG2059 Geocomputation 2016/17</h1></center>\n",
    "\n",
    "<h1><center>Practical 10a: Analysis of Relationships in Weather and Air Quality Data</h1></center>\n",
    "\n",
    "<p><center><i>James Millington, 27 November 2016</i></center>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Practical 10 is split into two notebooks - one examining relationships in Heathrow Weather and Air Quality data (10a) and one examining relationships in NS-SeC and house prices data (10b). The two notebooks are self-contained and can be used independently. You should decide which sets of data you most likely want to use for your final report, and work through the corresponding notebook during supervised practical time. This will give you the basics of analyses that you can then build on for your final report. Of course, you are welcome to work through both notebooks, although you are unlikely to be able to complete both during class time. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helper Functions\n",
    "\n",
    "Before getting to the data and code in this notebook, you should first run the code in the next three code blocks. These code blocks:\n",
    "\n",
    "1. import packages required for functionality in the remiander of the notebook and set `matplotlib` font parameters \n",
    "2. define a function to help interpret OLS regression output\n",
    "3. define a function to plot a histogram to file\n",
    "\n",
    "Take a quick look at the code when running these blocks, but don't spend too much time as we will return to look at the function definitions more closely later in the notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#import packages required for functionalit below and set matplotlib font parameters \n",
    "import os\n",
    "import pandas as pd\n",
    "import seaborn as sb      \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt    #see http://matplotlib.org/users/pyplot_tutorial.html\n",
    "import statsmodels.api as sm       #see http://statsmodels.sourceforge.net/stable/  \n",
    "\n",
    "#set matplotlib font params\n",
    "plt.rcParams['axes.titlesize'] = 20\n",
    "plt.rcParams['axes.labelsize'] = 20\n",
    "plt.rcParams['xtick.labelsize'] = 14\n",
    "plt.rcParams['ytick.labelsize'] = 14\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#define function to help interpret OLS regression output\n",
    "def mod_diagnostics(model, data):\n",
    "    \n",
    "    \"\"\"\n",
    "    Output to file model diagnostics for an OLS model\n",
    "    \n",
    "    Input:\n",
    "        model - statsmodels.regression.linear_model.OLS object\n",
    "        data  - pandas.DataFrame containing data for model\n",
    "        \n",
    "    Output:\n",
    "        XX-XX-OLS_SampleXX_Summary.txt contains the model summary output\n",
    "        XX-XX-OLS_SampleXX_ResidHist.png is histogram of the residuals\n",
    "        XX-XX-OLS_SampleXX_StdResid.png is a plot of standardised residuals against fitted values\n",
    "        \n",
    "        if model is univariate: XX-XX_OLS_SampleXX_Regression.png is a scatter plot with regression line\n",
    "        \n",
    "    Requires:\n",
    "        statsmodels.api\n",
    "        pandas\n",
    "        numpy\n",
    "        matplotlib.pyplot\n",
    "    \"\"\"\n",
    "    \n",
    "    fitted = model.fit()\n",
    "    dep = model.endog_names\n",
    "    indep_names = \"\"\n",
    "    \n",
    "    #create a string containing list of indep names for output files\n",
    "    for name in model.exog_names[1:]:            #we don't want 0 element as that is the intercept\n",
    "        indep_names += \"{0}_\".format(name)\n",
    "\n",
    "\n",
    "    #Want to include name of DataFrame in the output filename but currently DataFrame does not have a name attribute\n",
    "    #So for now use nobs from fitted  (Dan potential solution: pass data in a dictionary and access the label)\n",
    "    samplesize = str(int(fitted.nobs))\n",
    "    \n",
    "    f1 = open(\"{0}-{1}OLS_Sample{2}_Summary.txt\".format(dep, indep_names, samplesize), \"w\")\n",
    "    f1.write(fitted.summary().as_text())\n",
    "    f1.close()\n",
    "\n",
    "    #calculate standardized residuals ourselves\n",
    "    fitted_sr = (fitted.resid / np.std(fitted.resid)) \n",
    "\n",
    "    #Histogram of residuals\n",
    "    ax = plt.hist(fitted.resid)\n",
    "    plt.xlabel('Residuals')\n",
    "    plt.savefig('{0}-{1}OLS-Sample{2}_ResidHist.png'.format(dep, indep_names, samplesize), bbox_inches='tight')\n",
    "    plt.close()\n",
    "\n",
    "    #standardized residuals vs fitted values\n",
    "    ax = plt.plot(fitted.fittedvalues, fitted_sr, 'bo')\n",
    "    plt.axhline(linestyle = 'dashed', c = 'black')\n",
    "    plt.xlabel('Fitted Values')\n",
    "    plt.ylabel('Standardized Residuals')                \n",
    "    plt.savefig('{0}-{1}OLS-Sample{2}_StdResid.png'.format(dep, indep_names, samplesize), bbox_inches='tight')\n",
    "    plt.close()\n",
    "  \n",
    "    \n",
    "    if(len(model.exog_names) == 2):  #univariate model (with intercept)\n",
    "            \n",
    "        indep = model.exog_names[1]\n",
    "        \n",
    "        #scatter plot with regression line \n",
    "        ax = plt.plot(data[indep], data[dep], 'bo')\n",
    "        x = np.arange(data[indep].min(), data[indep].max(), 0.1)    #list of values to plot the regression line using\n",
    "        plt.plot(x, fitted.params[1]*x + fitted.params[0], '-', c = 'black')  #plot a line using the standard equation with parms from the model\n",
    "        \n",
    "        plt.xlabel(indep)\n",
    "        plt.ylabel(dep)                \n",
    "        plt.savefig('{0}-{1}OLS_Sample{2}_Regression.png'.format(dep, indep, samplesize), bbox_inches='tight')\n",
    "        plt.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#define function to plot histogram to file\n",
    "def plot_hist(series):\n",
    "    \n",
    "    \"\"\"\n",
    "    Output to file a simple histogram\n",
    "    \n",
    "    Input:\n",
    "        series - pandas.Series containing data (may also be able to take a numpy array)\n",
    "        \n",
    "    Output:\n",
    "        XX-SampleXX-Hist.png - the histogram image\n",
    "        \n",
    "    Requires:\n",
    "        pandas\n",
    "        matplotlib.pyplot\n",
    "    \"\"\"\n",
    "\n",
    "    out_name = \"{0}-Sample{1}-Hist.png\".format(series.name, len(series))\n",
    "    plt.hist(series.dropna())\n",
    "    plt.xlabel(series.name)\n",
    "    plt.ylabel('Count')\n",
    "    plt.savefig(out_name, bbox_inches='tight')           #save the figure\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heathrow Weather and Air Quality Data\n",
    "\n",
    "The additional data we'll use with the Heathrow Weather data are air quality data have been downloaded from the Air Quality England [website](http://www.airqualityengland.co.uk/) (AQE 2016) for the [Hounslow Hatton Cross site](http://www.airqualityengland.co.uk/site/latest?site_id=HS7) (site HS7). This site was chosen as it is near Heathrow Airport. \n",
    "\n",
    "Air pollution is an important aspect of the ongoing argument about the construction of the third runway at Heathrow (e.g. GLA 2012). In particular, although Nitrogen Dioxide (NO2) concentrations around Heathrow, are lower than in the centre of London, they are still often above recommended levels (e.g. Heathrow 2012). By looking at relationships between weather and air quality we may begin to better understand the drivers of pollution.\n",
    "\n",
    "In Practical 9 we used code to clean the air quality data and join it to the weather data to create a single time-series of data. This code is copied in the next code block - you will need to run this code if you have not already done so to create the `HeathrowAQWeather2016.pkl` file. However, if you have already created `HeathrowAQWeather2016.pkl` you can skip that code and simply load the data into memory (the subsequent code block)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##ONLY run this code block if you have NOT already created HeathrowAQWeather2016.pkl in Practical 9\n",
    "\n",
    "#read weather data\n",
    "metDF = pd.read_pickle(\"CleanedHeathrowData2016.pkl\")\n",
    "\n",
    "#read the aq data\n",
    "aqDF = pd.read_csv('AirQuality_HattonCross_2016.csv', header=0, skiprows=3, usecols=range(0,10), skipfooter = 1, engine = 'python')\n",
    "\n",
    "#rename columns and drop those not needed\n",
    "aqCN = ['Date', 'Time', 'PM10', 'PM10_su', 'NO', 'NO_su', 'NO2', 'NO2_su', 'NOx', 'NOx_su']\n",
    "aqDF.columns = aqCN\n",
    "del aqDF['PM10_su'] \n",
    "del aqDF['NO_su']   \n",
    "del aqDF['NO2_su'] \n",
    "del aqDF['NOx_su'] \n",
    "\n",
    "#fix issue with AQ data time formatting and set TimeDateIndex\n",
    "aqDF.Time.replace(to_replace = '24:00:00', value = '00:00:00', inplace = True) \n",
    "aqDF[\"DT\"] = pd.to_datetime(aqDF.Date + aqDF.Time, format='%d/%m/%Y%H:%M:%S')  #create a new series containing a datetime object\n",
    "aqDF.index = aqDF[\"DT\"]     \n",
    "\n",
    "import datetime as dt\n",
    "oneday = dt.timedelta(days=1)                           #create a timedelta object of days = 1\n",
    "aqDF.loc[aqDF['Time'] == '00:00:00','DT'] += oneday     #update the DT cell for rows with Time == 00:00 by adding oneday\n",
    "del aqDF['DT'] \n",
    "aqDF = aqDF.drop_duplicates()\n",
    "\n",
    "aqfirstDate = aqDF.index[0]\n",
    "aqlastDate = aqDF.index[len(aqDF.Date) - 1]\n",
    "aqDF = aqDF.reindex(index=pd.date_range(start = aqfirstDate, end = aqlastDate, freq = '1H'), fill_value = None)   #http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.reindex.html\n",
    "\n",
    "#select only column we want for the final joined dataframe\n",
    "metDF = metDF.select_dtypes(include=['float64'])\n",
    "del metDF['WindGust'] \n",
    "del metDF['LocID']\n",
    "\n",
    "del aqDF['Date'] \n",
    "del aqDF['Time']\n",
    "\n",
    "#join!\n",
    "aqmetDF = aqDF.join(metDF)\n",
    "\n",
    "#reset pressure values that don't make sense\n",
    "aqmetDF.loc[aqmetDF.Pressure == 0] = None\n",
    "\n",
    "#write the data to file\n",
    "aqmetDF.to_pickle(\"HeathrowAQWeather2016.pkl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## run this code if you HAVE already created HeathrowAQWeather2016.pkl\n",
    "#import pandas as pd                                   #already imported above but would be needed otherwise\n",
    "aqmetDF = pd.read_pickle(\"HeathrowAQWeather2016.pkl\")  #assumes file is saved in the same folder as this notebook file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: the previous code block assumes `HeathrowAQWeather2016.pkl` is saved in the same folder as this notebook file, but it is possible to read (and write) data to other folders by specifying the 'path' we want to use. The following code block shows one way to do this (as James discussed in Week 9 lecture). **ONLY** run the next code block if you want to read data from a location other than the folder in which this notebook files is saved - it is more for your information for future use. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##ONLY run this code if you want to read data from a location other than the folder in which this notebook files is saved\n",
    "#import os             #already imported above but would be needed otherwise\n",
    "\n",
    "#set the path to the directory where we want to read and save from\n",
    "path = os.path.join(os.path.expanduser(\"~\"),\"Google Drive\",\"Teaching\",\"2016-17\",\"Undergrad\",\"Geocomp\",\"Week10\",\"Practical\")\n",
    "os.chdir(path)\n",
    "\n",
    "#the following line would now read the pkl file from my Week 9 folder (specified above)\n",
    "#nsvalDF.to_pickle(\"HeathrowAQWeather2016.pkl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exploratory Analysis and Plotting\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step in our analysis should be looking at the shape of the distributions of variables and any possible relationships between variables. Let's use a pairplot to do this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sb.pairplot(aqmetDF.dropna(axis = 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Describe the distributions and relationships of the **air quality** variables (edit this text block):\n",
    "\n",
    "**A: **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We should also check for missing data in our time series. We'll use a loop to automate the plotting of a time series for each of the air quality variables: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for name in ['PM10', 'NO', 'NO2', 'NOx']:\n",
    "    \n",
    "    fig = plt.figure(name)   #create a new figure (closes any existing)\n",
    "    fig = aqmetDF[name].plot()  \n",
    "    plt.xlabel('Date')\n",
    "    plt.ylabel(name)\n",
    "    plt.legend(loc = 0)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the time series plots you can see that there is a fair chunk of missing data for June (and there also other smaller chunks of missing data in other places). For your final report you may want to think about what time period you want to analyse, and you may also want to do some interpolation. For the remainder of this practical we will use data for September only (without any interpolation). \n",
    "\n",
    "To subset for September:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "aqmetDF_sept = aqmetDF['2016-09']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "Repeat the exploatory plots we just did for the entire data set (i.e. pairplot and time series plots) for the `aqmetDF_sept` dataframe just created:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "For the September time series you just plotted, what differences or similarities can you see between the plots? (edit this text block to answer)\n",
    "\n",
    "**A: **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correlation and Relationships\n",
    "\n",
    "Now that we have all the weather and air quality data together in one DataFrame we can start to look at correlations and relationships. Let’s look directly at some correlation matrices for the variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "corrmat = aqmetDF_sept.corr()\n",
    "print \"Pearson correlation coefficient matrix:\", '\\n', corrmat, '\\n'\n",
    "\n",
    "corspmat = aqmetDF_sept.corr(method = \"spearman\")\n",
    "print \"Spearman rank correlation coefficient matrix:\", '\\n', corspmat, '\\n'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you recall, the degree of linearity in the relationships between variables (e.g. look back to your pairplot) should determine which type of correlation we should use. We’ll focus here on NO2 as that has been highlighted as a particular issue for air quality around Heathrow (but feel free to examine other variables for your final report). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "For the weather variable with the strongest correlation to NO2, create jointplots to examine their relationships (replace `???` with the names of the appropriate variable):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import spearmanr\n",
    "sb.jointplot(???, \"NO2\", data=aqmetDF_sept, stat_func=spearmanr)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Can you tell if the relationship is linear? Is it a Positive or Negative relationship? How would you describe the distribution of the variables? (edit this text block to answer)\n",
    "\n",
    "**A: **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As all these variables have been measured over time we can also look at their relationships through time. For now, we'll focus on the weather variable with the strongest relationship to NO2. Check you understand how the code below produces a two-axis time series of NO2 with Wind Speed (consult the comments):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = sb.axes_style()                                                #save default sb style settings\n",
    "sb.set_style(\"darkgrid\", {'axes.grid': False})                     #turn off the grid\n",
    "\n",
    "fig, ax1 = plt.subplots()                                          #set up the figure and left axis\n",
    "ax1.set_ylim(0,25)                                                 #modify y-axis limits   (to create room for legends later) \n",
    "\n",
    "name1 = \"WindSpeed\"\n",
    "ax1 = aqmetDF_sept[name1].plot(style='-r', label = name1)          #plot \n",
    "ax1.set_ylabel(name1, color='r')                                   #set first y-axis label\n",
    "\n",
    "name2 = \"NO2\"\n",
    "ax2 = ax1.twinx()                                                  #create a twin of the left axis for right\n",
    "ax2.set_ylim(0,100)    \n",
    "ax2 = aqmetDF_sept[name2].plot(style='-b', label = name2)          #plot\n",
    "ax2.set_ylabel(name2, color='b')                                   #change second axis label\n",
    "\n",
    "ax2.legend(loc = 1)                                                #add legend \n",
    "ax1.legend(loc = 2)                                                #add legend \n",
    "ax1.set_xlabel(\"Date\", color = 'black')\n",
    "\n",
    "#uncomment next two lines to save the figure to an image file (e.g. for use in reports)\n",
    "#plt.savefig('{0}{1}.png'.format(name1, name2), bbox_inches='tight') #save the figure\n",
    "#plt.close()                                                        #close plot\n",
    "                \n",
    "sb.set_style(a)                                                    #revert style back to defaults\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the plot just produced you should be able to see how NO2 is generally greater when Wind Speed is lower (indicating the negative relationship found through the correlation above). However, the pattern could be clearer – there’s a lot of variation from hour to hour which makes the relationship more difficult to see. Maybe the use of running means would help here… "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Create (and save to file) a two-axis plot with the 12-hour running mean of WindSpeed on the left axis and the 12-hour running mean of NO2 on the right axis. Copy the code from the last code block and edit appropriately to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 12-hour running mean plot has made it much easier to see a relationship between the two variables. Consequently, it's likely that using such a running mean will help to strengthen the correlation between the variables. To check this, we first need to calculate the running mean for all our variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "aqmetDF_sept_12mw = aqmetDF_sept.rolling(window=12, center = True).mean() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check you understand what the dataframe just created contains (e.g. by using code to find out about it). Now we can check the correlation between the 12-hour running means of NO2 and Wind Speed. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Create a `jointplot` to visualise the relationship between the 12-hour running means of Wind Speed and NO2 (use the `aqmetDF_sept_12mw` dataframe just created):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Create correlation matrices (pearson and spearman) for the 12-hour running means of all the variables (use the `aqmetDF_sept_12mw` dataframe just created):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Compare the correlations between NO2 and Windspeed for the raw data and the 12-hour running mean. How have the correlations changed? Are they stronger or weaker? What is your physical interpretation for any differences? (edit this code block to answer)\n",
    "\n",
    "**A: **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple Linear Regression\n",
    "\n",
    "We seem to have identified a reasonably good relationship between Wind Speed and NO2. Maybe we could use regression to try to predict what level of NO2 we would expect at Heathrow over a 12 hour period from the mean wind speed over that period. Regression would allow us to do this.\n",
    "\n",
    "To fit a regression we can use the OLS method in the statsmodels package. The `statsmodels.api` was imported with alias `sm` at the top of this notebook, so we can use it to fit the regression between our variables as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#create OLS object\n",
    "NO2_WS_RM_mod = sm.OLS.from_formula(\"NO2 ~ WindSpeed\", aqmetDF_sept_12mw, missing = 'drop')  #use the missing argument with value drop to tell python to ignore missing data\n",
    "#fit the regression\n",
    "fitted_NO2_WS_RM_mod = NO2_WS_RM_mod.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how there are two steps to fitting the regression model. First, we create a OLS model object by specifying the 'formula' and the data to use - see that 'formula' does not use the `=` symbol and instead relates variables using `~`. \n",
    "\n",
    "The second line above then actually fits the regression model (using the `fit` method) and puts this in a 'fitted model' object. We can then get a summary of the regression model using the `summary` method with the fitted model object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print \"Summary of NO2_WS_RM_mod\", '\\n', fitted_NO2_WS_RM_mod.summary()           #output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "Using your model output, answer the following questions. You will find it useful to refer to your lecture notes (including Week 9) and [Johnson (2014)](http://connor-johnson.com/2014/02/18/linear-regression-with-python/) provides a nice overview of the model output and what it all means:\n",
    "\n",
    "Qa)\tWhat is the model r2? Does this make sense given the value of Pearson’s r you got earlier? \n",
    "\n",
    "**A: **\n",
    "\n",
    "Qb)\tWhat is the value of the model intercept? \n",
    "\n",
    "**A: **\n",
    "\n",
    "Qc)\tWhat is the parameter value for the slope of the regression? Is it statistically significantly different from zero (at 95% confidence)? \n",
    "\n",
    "**A: **\n",
    "\n",
    "Qd)\tReplace ???: For every 1 m/s of Wind Speed increase, we would expect ??? µg/m3 of NO2 ??? [more/less]\n",
    "\n",
    "**A: **\n",
    "\n",
    "Qe)\tAccording to the output are the residuals normally distributed (look at `Prob(Omnibus`) and `Prob(JB)` and compare to lecture slides)? \n",
    "\n",
    "**A: **\n",
    "\n",
    "Qf)\tAccording to the Durbin-Watson test is there structure (auto-correlation) in the residuals? \n",
    "\n",
    "**A: **\n",
    "\n",
    "Qg)\tGiven your answers to Qe) and Qf) do you think the assumptions of linear regression have been violated? \n",
    "\n",
    "**A: **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check further if the assumptions of linear regression have been violated, we can look at plots of the residuals. The `mod_diagnostics` function defined above provides these for us..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### More Model Diagnostics\n",
    "\n",
    "Go and look at the `mod_diagnostics` function at the top of this notebook now to check you understand it. Identify where it does the following:\n",
    "- Fits the model passed to it \n",
    "- Writes the fit model summary to file on disk\n",
    "- Saves the histogram of the residuals to file using a string format\n",
    "- Creates a scatter plot of standardized residuals\n",
    "- Only creates a scatter plot with a regression line when the number of independent variables is equal to 1 (and think about why we don’t try to do this when we have more than one independent variable)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the `mod_diagnostics` function to create diagnostic plots for the `NO2_WS_RM_mod` model. Remember two things about `mod_diagnostics` function:\n",
    "1. it takes a OLS model object as an argument, NOT a fitted model object\n",
    "2. it writes its output to file (so you'll need to check your hard disk for output) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mod_diagnostics(NO2_WS_RM_mod, aqmetDF_sept_12mw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Go and look at the output just created by `mod_diagnostics` - if you can't find it on file check above if you have set the working directory somewehere else. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the histogram of the residuals you should be able to see that residuals are positively skewed. But does this matter? Not really, as we have a large sample size here. \n",
    "\n",
    "What is more problematic is the structure in the residuals, which you should be able to see in the standardised residuals plot. So we need to do a little more work to get a valid regression. \n",
    "\n",
    "There are a couple of things we might do:\n",
    "1.\tWe should think about how we may have introduced structure into the input variables ourselves (which can be the root of structure in the residuals)\n",
    "2.\tWe should think about whether we can transform the distribution of the input variables to make them more normal (which often helps to make residuals more normal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sampling from the Running Mean\n",
    "\n",
    "The structure we can see in the residuals of the last model is also present in the regression plot itself (as hopefully you can see for yourself). We may have introduced this structure ourselves by taking the running mean. Think about it: the value of every hour in our running mean data is weighted by nearby hours. These values are no longer independent. \n",
    "\n",
    "One way around might be to sample the running mean data to only use independent values. Given we used a 12-hour running mean, using data for every 12th hour should do this... (think about that). \n",
    "\n",
    "The code below creates a new dataframe containing only values for 6am and 6pm (by building on [this](http://stackoverflow.com/a/10567298) SO answer). Check you generally understand what is going on here and what the new dataframe contains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hour = aqmetDF_sept_12mw.index.hour\n",
    "selector = ((hour == 6) | (hour == 18))   \n",
    "sixes = aqmetDF_sept_12mw[selector]\n",
    "sixes.info()\n",
    "print sixes.head()\n",
    "\n",
    "corrmat_sixes = sixes.corr()\n",
    "print \"\\nPearson correlation coefficient matrix:\", '\\n', corrmat_sixes, '\\n'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Answer the following questions (edit this code block):\n",
    "\n",
    "What has happened to the correlation between NO2 and Wind Speed for the sampled data?\n",
    "\n",
    "**A: **\n",
    "\n",
    "What has happened to the number of values in this dataframe compared to the previous one?\n",
    "\n",
    "**A: **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Create a joint plot from the sampled dataframe just created to check that the relationship between the NO2 and Wind Speed is still roughly linear:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, looks good, let's try fitting a regression to these data to see what we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "NO2_WS_six_mod = sm.OLS.from_formula(\"NO2 ~ WindSpeed\", sixes)  #use the missing argument with value drop to tell python to ignore missing data\n",
    "mod_diagnostics(NO2_WS_six_mod, sixes)\n",
    "print NO2_WS_six_mod.fit().summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Compare the results for this regression on the sampled data to the previous model:\n",
    "\n",
    "Q: (How) have model parameters changed? \n",
    "\n",
    "**A: **\n",
    "\n",
    "Q: How does the model fit comapre to the previous model? _Hint: compare the confidence intervals_ \n",
    "\n",
    "**A: **\n",
    "\n",
    "Q: Have we successfully removed structure from the residuals? \n",
    "\n",
    "**A: **\n",
    "\n",
    "Q: Are residuals normally distributed? \n",
    "\n",
    "**A: **\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transform input variables\n",
    "\n",
    "There is some ambiguity about whether the residuals in our new regression are normal. As our sample size is much smaller than previously we should maybe take this quite seriously and try our second option from above (transformation) to further improve the situation. \n",
    "\n",
    "You can see from the last jointplot we made that the distribution of the NO2 data is heavily positively skewed - if we take a log transform of these data and use them in our model it may help with the issue of the normality of residuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sixes['logNO2'] = np.log(sixes['NO2'])     #transform the data by taking natural log\n",
    "\n",
    "#output some plots to quickly compare \n",
    "fig1 = sixes.logNO2.hist()\n",
    "plt.title('Log NO2')\n",
    "\n",
    "fig2 = plt.figure()\n",
    "fig2 = sixes.NO2.hist()\n",
    "plt.title('NO2')\n",
    "\n",
    "#plot_hist(sixes['logNO2'])      #note that the plot_hist function defined above plots a histogram in an image file on disk\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The log of NO2 seems to be less skewed than the original NO2 data so maybe that will help us to meet the assumptions of linear regression (given a low sample size for the sampled data)\n",
    "\n",
    "#### Task\n",
    "Fit a regression model for the logNO2 data for the sample running mean data, using Wind Speed as a predictor. Print the output in this notebook and to file (e.g. using the `mod_diagnostics` function):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Check the output of the model summary and plots of residuals \n",
    "\n",
    "Qa)\tLooking at the plots of residuals do you think this model meets the assumptions of linear regression better?\n",
    "\n",
    "**A: **\n",
    "\n",
    "Qb)\tHow much difference is there between the explanatory power (r2) of this new model and the previous one that use all the running mean data points?\n",
    "\n",
    "**A: **\n",
    "\n",
    "Qc)\tUsing r2 how much of the variation in mean 12-hourly NO2 can we explain by mean 12-hourly wind speed? \n",
    "\n",
    "**A: **\n",
    "\n",
    "Qd) Replace ?? to answer: For every 1 m/s of Wind Speed increase, we would expect a ???% NO2 ??? [increase/decrease]\n",
    "\n",
    "_Hint: remember the dependent variable is log(NO2) and read [this](http://stats.stackexchange.com/a/18639) CV answer and see Table 2 of Lin et.al (2013)_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiple Linear Regression\n",
    "\n",
    "Finally, how about if we wanted to include multiple variables in our model to predict NO2? For example, what if we wanted include the second-most correlated weather variable with NO2 in our model? \n",
    "\n",
    "First we need to check for colinearity between the predictor variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "In the last correlation matrix created here, check if the second-most correlated variable with NO2 meets the general assumption about colinearity (with Wind Speed) highlighted in previous lectures.\n",
    "\n",
    "Q: If we were to include Wind Speed and the second-most correlated variable with NO2 in a multiple linear regression model to predict log(NO2), would we be violating the assumption of independence between the predictors? Why?\n",
    "\n",
    "**A: **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming there's no colinearity, we fit and summarise the multiple regression model in a very similar way to our simple linear regression. We just change the equation to include an extra variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "logNO2_WS_Pr_sixes_mod = sm.OLS.from_formula(\"logNO2 ~ WindSpeed + Pressure\", data = sixes)  #use the missing argument with value drop to tell python to ignore missing data\n",
    "mod_diagnostics(logNO2_WS_Pr_sixes_mod, sixes)\n",
    "print logNO2_WS_Pr_sixes_mod.fit().summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Task\n",
    "\n",
    "Qa) Do you think the model meets the assumptions of linear regression?\n",
    "\n",
    "**A: **\n",
    "\n",
    "Qb) What percentage of variance in logNO2 is explained by this model?\n",
    "\n",
    "**A: **\n",
    "\n",
    "Qc) What percentage of variance in logNO2 is explained by Pressure?\n",
    "\n",
    "**A: **\n",
    "\n",
    "Qd) Replace ???: If Wind Speed is ???, for every 1 mb of Pressure increase we would expect a ??? % NO2 ??? [increase/decrease]\n",
    "\n",
    "*Hint: see [this page](http://sites.stat.psu.edu/~ajw13/stat200/mos/12_multregr/12_multregr_print.html) for help*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final Project\n",
    "\n",
    "Think about what the results above tell us about the key weather variables that influence air quality.\n",
    "\n",
    "Start thinking about your final project and what data you might analyse for it. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### References\n",
    "- AQE (2016) _Air Quality England_ [Online] Available at: http://www.airqualityengland.co.uk/ \n",
    "- GLA (2012) Air and noise pollution around a growing Heathrow Airport [Online] Available from: http://www.london.gov.uk/mayor-assembly/london-assembly/publications/tackling-air-and-noise-pollution-around-heathrow \n",
    "- Heathrow (2012) Heathrow Air Quality [Online] Available from: http://www.heathrow.com/file_source/Company/Static/PDF/Communityandenvironment/air-quality-strategy_LHR.pdf\n",
    "- Johnson (2014) Linear Regression with Python [Online] Available at: http://connor-johnson.com/2014/02/18/linear-regression-with-python/\n",
    "- Lin et al. (2013) Too Big to Fail: Large Samples and the p-Value Problem _Information Systems Research_ 24 906–917 DOI: [10.1287/isre.2013.0480](http://dx.doi.org/10.1287/isre.2013.0480)\n",
    "- Lumley et al. (2002) _Annu. Rev. Public Health_ 23:151–69 DOI: [10.1146/annurev.publhealth.23.100901.140546](http://dx.doi.org/10.1146/annurev.publhealth.23.100901.140546)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}