04-Scoring Neighbourhoods.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scoring Neighbourhoods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The scoring methodology is based on 4 variables (household income, house prices, % of \"favourable\" occupations, % of graduates) which are combined using PCA to score each neighbourhood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below code loads in relevant function libraries and the data to build the SES scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Needed on a Mac\n",
    "import matplotlib as mpl\n",
    "mpl.use('TkAgg')\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# For reproducibility\n",
    "import random\n",
    "import numpy as np\n",
    "r_state = 42\n",
    "random.seed(r_state) \n",
    "np.random.seed(r_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import pysal as ps\n",
    "import seaborn as sns\n",
    "\n",
    "from sklearn import decomposition  \n",
    "from sklearn.preprocessing import scale  \n",
    "from sklearn import preprocessing \n",
    "from sklearn import linear_model\n",
    "from sklearn import model_selection\n",
    "#from sklearn import cross_validation\n",
    "\n",
    "from scipy.stats import boxcox\n",
    "from scipy.stats import spearmanr\n",
    "from scipy.stats import pearsonr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lkp = os.path.join('data','lkp')\n",
    "src = os.path.join('data','src')\n",
    "\n",
    "analytical = os.path.join('data','analytical')\n",
    "canonical  = os.path.join('data','canonical')\n",
    "converted  = os.path.join(canonical,'converted')\n",
    "greenspace = os.path.join(canonical,'greenspace')\n",
    "dwelling   = os.path.join(canonical,'dwellings')\n",
    "travel     = os.path.join(canonical,'travel')\n",
    "household  = os.path.join(canonical,'households')\n",
    "housing    = os.path.join(canonical,'housing')\n",
    "work       = os.path.join(canonical,'work')\n",
    "scores     = os.path.join(canonical,'scores')\n",
    "\n",
    "for d in [analytical,canonical,converted,greenspace,dwelling,travel,household,housing,work,scores]:\n",
    "    if not os.path.exists(d):\n",
    "        os.makedirs(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_checks(df, selected_cols=None, prefix='Test'):\n",
    "    sns.set(rc={\"figure.figsize\": (12, 3)})\n",
    "    if not selected_cols:\n",
    "        selected_cols = df.columns\n",
    "    for d in selected_cols:\n",
    "        print(\"Working on \" + d)\n",
    "        fig = plt.figure(d)\n",
    "        sns.distplot(df[d], color='green', hist=True, rug=True, norm_hist=False)\n",
    "        fig = plt.gcf() # *G*et the *C*urrent *F*igure environment so that the next command works\n",
    "        plt.savefig(\"{0}-{1}-Check.pdf\".format(prefix, d.replace(':',' - ')), bbox_inches=\"tight\")\n",
    "        plt.close()\n",
    "    print(\"Done.\")\n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Scoring Data for 2001 and 2011"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build the Scoring Data Set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df11 = pd.DataFrame()\n",
    "df01 = pd.DataFrame()\n",
    "for d in ['Income.csv','Values.csv','Occupations.csv','Qualifications.csv']:\n",
    "    if 'Values' in d:\n",
    "        loc = housing\n",
    "    else:\n",
    "        loc = work\n",
    "    tmp11_df = pd.read_csv(os.path.join(loc,d.replace('.csv','-2011.csv')))\n",
    "    tmp01_df = pd.read_csv(os.path.join(loc,d.replace('.csv','-2001.csv')))\n",
    "    \n",
    "    # Needed for GeoConvert-ed data, will have no effect on other dfs\n",
    "    tmp01_df.rename(columns=\n",
    "                  {'2011 census lower super output areas and data zones  (2001 codes used in scotland)':'lsoacd'},\n",
    "                 inplace=True)\n",
    "    \n",
    "    if df11.shape[0] == 0:\n",
    "        df11 = tmp11_df\n",
    "        df01 = tmp01_df\n",
    "    else: \n",
    "        df11 = pd.merge(df11, tmp11_df, how='outer', left_on='lsoacd', right_on='lsoacd')\n",
    "        df01 = pd.merge(df01, tmp01_df, how='outer', left_on='lsoacd', right_on='lsoacd')\n",
    "\n",
    "print(\"Shape of 2001 data frame: \" + str(df01.shape))\n",
    "print(\"Shape of 2011 data frame: \" + str(df11.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rename = {\n",
    "    'Mean':'Mean_inc',\n",
    "    'Median_x':'Median_inc',\n",
    "    'Median_y':'Median_hp',\n",
    "    'Total_x':'Total (Occupations)',\n",
    "    'Total_y':'Total (Qualifications)'\n",
    "}\n",
    "df01.rename(columns=rename, inplace=True)\n",
    "df11.rename(columns=rename, inplace=True)\n",
    "print(\"Columns renamed to remove ambiguity.\")\n",
    "\n",
    "#  Set the index of dataframe to LSOA\n",
    "df01.set_index('lsoacd', inplace=True)\n",
    "df11.set_index('lsoacd', inplace=True)\n",
    "print(\"Datasets indexed to LSOA\")\n",
    "\n",
    "df11.sample(3, random_state=r_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transform Property Prices and Income\n",
    "\n",
    "These are so heavily skewed that some kind of transformation may be required to prevent high prices or incomes dominating the scoring metric. This is something of a judgement call, so for completeness we've included three possible ways of working with the data but have selected only one of them for our analysis.\n",
    "\n",
    "<span style=\"color:red;font-weight:bolder;size:14pt;\">You can run _all_ of the following options, but you can only select _one_ for further work in the PCA Processing section.</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Option 1: No Transform\n",
    "\n",
    "In this option we pass through household income and property price 'as is'. These will consequently weigh more heavily in the final score and highlight slightly different features across London."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df_transformed = pd.DataFrame({\n",
    "    'hpu_01' : df01['Median Property Price'],\n",
    "    'hpu_11' : df11['Median Property Price'],\n",
    "    'hhu_01' : df01['Median Income'],\n",
    "    'hhu_11' : df11['Median Income']\n",
    "}, index=df01.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Option 2: Box-Cox Transform\n",
    "\n",
    "In this option with use Box-Cox transforms on these two variables so that they are pseudo-normal. Note that we need to use the same transform on both years so as to ensure that the results are comparable across 2001/2011. Since skills and occupation remain skewed (albeit much less so) this will tend to highlight changes in 'human factor' issues and deemphasise financial changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# We want to use the same transform so that change is detectable\n",
    "hpb01, lmd01a = boxcox(df01['Median Property Price'])  #  Process 2001 data\n",
    "print(\"2001 property price transform lambda: \" + str(lmd01a))\n",
    "hpb11 = boxcox(df11['Median Property Price'], lmbda=lmd01a)\n",
    "print(\"Property prices transformed using same Box-Cox lambda.\")\n",
    "\n",
    "print(\" \")\n",
    "\n",
    "hhb01, lmd01b = boxcox(df01['Median Income'])  #  Process 2001 data\n",
    "print(\"2001 income transform lambda: \" + str(lmd01b))\n",
    "hhb11 = boxcox(df11['Median Income'], lmbda=lmd01b)\n",
    "print(\"Household income transformed using same Box-Cox lambda.\")\n",
    "\n",
    "df_transformed2 = pd.DataFrame({\n",
    "    'hpb_01': hpb01,\n",
    "    'hpb_11': hpb11,\n",
    "    'hhb_01': hhb01,\n",
    "    'hhb_11': hhb11\n",
    "}, index=df01.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df_transformed = pd.merge(df_transformed, df_transformed2, how='inner', left_index=True, right_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Option 3: In-Between Transform\n",
    "\n",
    "Transforming for normality has quite a strong impact on the score, so possibly it would be better to select an intermeidate transformation that has less of an impact. Such as a fairly common log-transform on prices, and $x^{2/3}$ on incomes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# We want to use the same transform so that change is detectable\n",
    "hpl01 = np.log(df01['Median Property Price'])  #  Process 2001 data\n",
    "hpl11 = np.log(df11['Median Property Price'])\n",
    "print(\"Property prices transformed using natural log.\")\n",
    "\n",
    "print(\" \")\n",
    "\n",
    "hhl01 = np.power(df01['Median Income'], 2.0/3.0)  #  Process 2001 \n",
    "hhl11 = np.power(df11['Median Income'], 2.0/3.0)\n",
    "print(\"Household income transformed using same exponent.\")\n",
    "\n",
    "df_transformed3 = pd.DataFrame({\n",
    "    'hpl_01': hpl01,\n",
    "    'hpl_11': hpl11,\n",
    "    'hhl_01': hhl01,\n",
    "    'hhl_11': hhl11\n",
    "}, index=df01.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df_transformed = pd.merge(df_transformed, df_transformed3, how='inner', left_index=True, right_index=True)\n",
    "print(\"Final shape: \" + str(df_transformed.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Slow. Uncomment to explore further.\n",
    "#plot_checks(df_transformed, prefix='Transforms')\n",
    "print(\"Done printing out results of transforms.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df_transformed.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Occupation Share\n",
    "\n",
    "Each LSOA is assessed on the share of occupations in terms of their likely contribution to neighbourhood change. Note that shares are not transformed as they are much less likely to be as heavily skewed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Process occupational data\n",
    "def process_occ_data(df):\n",
    "    #  Columns of interest\n",
    "    occ = ['Managerial','Professional','Technical','Administrative','Skilled','Personal Service','Customer Service','Operators','Elementary']\n",
    "    \n",
    "    # Integrate results into Occupations datasets -- \n",
    "    # right now we don't replicate Jordan's approach of\n",
    "    # grouping them into 'knowledge worker' and 'other'. \n",
    "    # We've added calculation of the HHI since it's a nice \n",
    "    # measure of diversity and might help us to unpick\n",
    "    # different types of change.\n",
    "    occ_data = pd.DataFrame()\n",
    "    for c in occ:\n",
    "        # Only if we want percentages for each group\n",
    "        #occ_data[c+'_pct']   = (df.loc[:,c] / df.loc[:,'Total_occ'])*100\n",
    "        occ_data[c+'_share'] = (df.loc[:,c] / df.loc[:,'Total (Occupations)'])**2\n",
    "    \n",
    "    # Calculate the HHI to get at dominance by a group/trade\n",
    "    #occ_data['hhi_occ'] = occ_data[[s for s in occ_data.columns if '_share' in s]].sum(axis=1)\n",
    "    \n",
    "    # Drop the share columns\n",
    "    occ_data.drop([s for s in occ_data.columns if '_share' in s], axis=1, inplace=True)\n",
    "    \n",
    "    # Add the 'knowledge worker' share -- this is columns 0-2 of the data frame\n",
    "    occ_data['kw_pct'] = (df.loc[:,occ[0:3]].sum(axis=1) / df.loc[:,'Total (Occupations)'])*100\n",
    "    \n",
    "    return occ_data\n",
    "\n",
    "occ01 = process_occ_data(df01)  #  Processed 2001 occupation data\n",
    "occ11 = process_occ_data(df11)  #  Processed 2011 occupation data\n",
    "\n",
    "print(\"Allows you to examine how occupations data will be processed.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "occ01.sample(3, random_state=r_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#g = sns.jointplot(y=occ01.hhi_occ, x=occ01.kw_pct, kind='kde', stat_func=spearmanr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Qualification Share\n",
    "\n",
    "Each LSOA is scored on the % of graduates (or equivalent), which corresponds to \"level 4-5\" in the Census methodology. Note that shares are not transformed at this stage as they're much less likely to be heavily skewed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Process qualifications data\n",
    "def process_quals_data(df):\n",
    "    # Columns of interest\n",
    "    quals  = ['No qualifications','Level 1','Level 2','Level 3','Level 4 and above','Other']\n",
    "    squals = ['Students In employment','Students Unemployed','Students inactive'] # Not currently used\n",
    "    \n",
    "    #  Integrate results into Qualifications datasets  -- \n",
    "    # right now we don't replicate Jordan's approach of\n",
    "    # grouping them into 'highly educated' and 'other'. \n",
    "    # We've added calculation of the HHI since it's a nice \n",
    "    # measure of diversity and might help us to unpick\n",
    "    # different types of change.\n",
    "    quals_data = pd.DataFrame()\n",
    "    for c in quals:\n",
    "        #quals_data[c+'_pct']   = (df.loc[:,c] / df.loc[:,'Total_qual'])*100\n",
    "        quals_data[c+'_share'] = (df.loc[:,c] / df.loc[:,'Total (Qualifications)'])**2\n",
    "    \n",
    "    # Calculate the HHI to get at dominance by a group/trade\n",
    "    #quals_data['hhi_quals'] = quals_data[[s for s in quals_data.columns if '_share' in s]].sum(axis=1)\n",
    "    \n",
    "    # Drop the share columns\n",
    "    quals_data.drop([s for s in quals_data.columns if '_share' in s], axis=1, inplace=True)\n",
    "    \n",
    "    # The 'highly educated' share -- this is columns 0-2 of the data frame\n",
    "    quals_data['he_pct'] = (df.loc[:,quals[4]] / df.loc[:,'Total (Qualifications)'])*100\n",
    "    \n",
    "    return quals_data\n",
    "\n",
    "qual01 = process_quals_data(df01)  #  Qualification data 2001\n",
    "qual11 = process_quals_data(df11)  #  Qualification data 2011 \n",
    "\n",
    "print(\"Allows you to examine how qualifications data will be processed.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "qual01.sample(3, random_state=r_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#g = sns.jointplot(y=qual01.hhi_quals, x=qual01.he_pct, kind='kde', stat_func=spearmanr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking on Distributions of Other Variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "quals_01 = process_quals_data(df01)\n",
    "quals_11 = process_quals_data(df11)\n",
    "occ_01   = process_occ_data(df01)\n",
    "occ_11   = process_occ_data(df11)\n",
    "\n",
    "df_test = pd.concat([quals_01, quals_11, occ_01, occ_11], join='inner', axis=1)\n",
    "df_test.columns = ['Qualifications 2001','Qualifications 2011','Occupations 2001','Occupations 2011']\n",
    "df_test.sample(3, random_state=r_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Slow. Uncomment to explore futher.\n",
    "#plot_checks(df_test, prefix='Transforms')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA Processing\n",
    "\n",
    "PCA was used to combine data on the four variables into a single score.  First, Principal Components Analysis finds the single component which maximises the variance in the data.  Data for the scoring variables from both the 2001 and 2011 datasets were combined (vertically) into a single dataset which was used to find the PCA transformation.\n",
    "\n",
    "<span style=\"color:red;font-weight:bolder\">You will have a choice to make here: which transform (if any) do you want to use on the skewed data? Note that your choice here should be the same as the one used in notebook 6; there is nothing, however, preventing your running this next section once using each type of transform so that you output all three scoring files and can more quickly experiment with different options in notebook 6.</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Don't forget to set the house price and income scores to the \n",
    "# transform/non-transform that you want to use!\n",
    "to_use = 'Untransformed' # Choices: ['Untransformed','Box-Cox','Log']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#  Indicators\n",
    "quals_score_01    = process_quals_data(df01)  #  Qualifications\n",
    "occ_score_01      = process_occ_data(df01)    #  Occupation\n",
    "quals_score_11    = process_quals_data(df11)  #  Qualifications\n",
    "occ_score_11      = process_occ_data(df11)    #  Occupation\n",
    "\n",
    "house_pr_score_01 = None  #  House Prices\n",
    "hh_inc_score_01   = None  #  Household Income\n",
    "house_pr_score_11 = None  #  House Prices\n",
    "hh_inc_score_11   = None  #  Household income\n",
    "\n",
    "if to_use == 'Untransformed':\n",
    "    house_pr_score_01 = df_transformed['hpu_01']  #  House Prices\n",
    "    hh_inc_score_01   = df_transformed['hhu_01']  #  Household Income\n",
    "    house_pr_score_11 = df_transformed['hpu_11']  #  House Prices\n",
    "    hh_inc_score_11   = df_transformed['hhu_11']  #  Household Income\n",
    "elif to_use == 'Box-Cox':\n",
    "    house_pr_score_01 = df_transformed['hpb_01']  #  House Prices\n",
    "    hh_inc_score_01   = df_transformed['hhb_01']  #  Household Income\n",
    "    house_pr_score_11 = df_transformed['hpb_11']  #  House Prices\n",
    "    hh_inc_score_11   = df_transformed['hhb_11'] \n",
    "elif to_use == 'Log':\n",
    "    house_pr_score_01 = df_transformed['hpl_01']  #  House Prices\n",
    "    hh_inc_score_01   = df_transformed['hhl_01']  #  Household Income\n",
    "    house_pr_score_11 = df_transformed['hpl_11']  #  House Prices\n",
    "    hh_inc_score_11   = df_transformed['hhl_11'] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Safety Checks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can't have NaN/Non-Finite values for PCA processing, so although we could fill in missing values by taking a weighted mean of the surrounding LSOAs wherever values are missing, the more effective (and less problematic) way is simply to drop them. Note that this means the data sets are _not_ the same size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "checks = {\n",
    "    \"Qualifications 2001\":quals_score_01,\n",
    "    \"Qualifications 2011\":quals_score_11,\n",
    "    \"Occupations 2001\":occ_score_01,\n",
    "    \"Occupations 2011\":occ_score_11,\n",
    "    \"House Prices 2001\":house_pr_score_01,\n",
    "    \"House Prices 2011\":house_pr_score_11,\n",
    "    \"Incomes 2001\":hh_inc_score_01,\n",
    "    \"Incomes 2011\":hh_inc_score_11\n",
    "}\n",
    "\n",
    "for k, v in checks.items():\n",
    "    if (np.isnan(v.values).any()):\n",
    "        print(\"Have null values in data set: \" + k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Create dataset of indicator data - 2001\n",
    "res_01 = pd.concat([house_pr_score_01,quals_score_01,occ_score_01,hh_inc_score_01], axis=1)\n",
    "res_11 = pd.concat([house_pr_score_11,quals_score_11,occ_score_11,hh_inc_score_11], axis=1)\n",
    "\n",
    "if to_use is 'Untransformed':\n",
    "    res_01.columns = ['House Prices 2001','Percentage with Level 4+ Qualifications 2001',\n",
    "                      'Percentage of Knowledge Workers 2001','Household Income 2001']\n",
    "    res_11.columns = ['House Prices 2011','Percentage with Level 4+ Qualifications 2011',\n",
    "                      'Percentage of Knowledge Workers 2011','Household Income 2011']\n",
    "else:\n",
    "    res_01.columns = ['House Prices 2001 (' + to_use + ' Transformed)','Percentage with Level 4+ Qualifications 2001',\n",
    "                      'Percentage of Knowledge Workers 2001','Household Income 2001 (' + to_use + ' Transformed)']\n",
    "    res_11.columns = ['House Prices 2011 (' + to_use + ' Transformed)','Percentage with Level 4+ Qualifications 2011',\n",
    "                      'Percentage of Knowledge Workers 2011','Household Income 2011 (' + to_use + ' Transformed)']\n",
    "\n",
    "# Create dataset of indicator data\n",
    "X_01 = res_01.values\n",
    "X_11 = res_11.values\n",
    "\n",
    "#  Join 2001 and 2011 datasets and sanity-check\n",
    "SES_inds = np.concatenate((X_01, X_11), axis=0)\n",
    "\n",
    "print(\"Any infinite values? \" + str(~np.isfinite(SES_inds).any()))\n",
    "print(\"Any NaN values? \" + str(np.isnan(SES_inds).any()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Median removal and Unit scaling\n",
    "scaler = preprocessing.RobustScaler()\n",
    "scaler.fit(SES_inds)\n",
    "SES_inds = scaler.transform(SES_inds)\n",
    "\n",
    "print(\"Data scaled and transformed.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This next section just gives us a sense of loadings and how much each component counts if we were performing full PCA:\n",
    "\n",
    "Explained Variance:\n",
    "* _Untransformed_: `[0.78793941, 0.151054, 0.04878766, 0.01221892]`\n",
    "* _Box Cox_: `[0.78728497, 0.15370062, 0.03948151, 0.01953289]`\n",
    "* _Log_: `[0.79813204, 0.14576833, 0.03743246, 0.01866718]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pca_full = decomposition.PCA()                           # Use all Principal Components\n",
    "pca_full.fit(SES_inds)                                   # Train model on data\n",
    "SES_full_T = pd.DataFrame(pca_full.transform(SES_inds))  # Transform data using model\n",
    "\n",
    "print(\"The amount of explained variance of the SES score using each component is...\")\n",
    "print(pca_full.explained_variance_ratio_)\n",
    "\n",
    "# Adapted from https://stackoverflow.com/questions/22984335/recovering-features-names-of-explained-variance-ratio-in-pca-with-sklearn\n",
    "i = np.identity(SES_inds.shape[1])  # identity matrix\n",
    "\n",
    "coef = pca_full.transform(i)\n",
    "\n",
    "loadings = pd.DataFrame(coef, index=res_01.columns)\n",
    "loadings.to_csv(os.path.join(scores,to_use + '-Loadings-2011.csv.gz'), compression='gzip', index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the transform we'll actually use. Notice that we limit `n_components` to 1 since we can only have a single score in our prediction code. So a PCA model is fitted and the transformation for the first component is used to assign each LSOA a score, before the 2001 and 2011 results are separated.\n",
    "\n",
    "I make it that the explained variances for each approach are:\n",
    "* _Untransformed_: 0.78794\n",
    "* _Box Cox_: 0.78728\n",
    "* _Natural log_: 0.79813"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Fitting PCA Model to derive SES score\n",
    "pca = decomposition.PCA(n_components=1)             # Only need 1st Principal Component\n",
    "pca.fit(SES_inds)                                   #  Train model on data\n",
    "SES_inds_T = pd.DataFrame(pca.transform(SES_inds))  #  Transform data using model\n",
    "\n",
    "print(\"The amount of explained variance of the SES score is: {0:6.5f}\".format(pca.explained_variance_ratio_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Split transformed data into 2001 and 2011 datasets\n",
    "#  Note the way we do this to deal with missing data (if any)\n",
    "scores_01 = SES_inds_T.loc[0:len(X_01)-1,0]\n",
    "scores_11 = SES_inds_T.loc[len(X_01):,0]\n",
    "\n",
    "# Create dfs from the two sets of scores\n",
    "res_01 = res_01.assign(scores=pd.Series(scores_01).values)\n",
    "res_11 = res_11.assign(scores=pd.Series(scores_11).values)\n",
    "\n",
    "#res.columns = ['LSOANM','PRICE-01','QUALS-01','OCC-01','INCOME-01','PRICE-11',\n",
    "#               'QUALS-11','OCC-11','INCOME-11','SES_01','SES_11']\n",
    "\n",
    "# Join them together so we've got a single df for 2001 and 2011\n",
    "res = res_01.merge(res_11, how='outer', suffixes=('_01','_11'), left_index=True, right_index=True)\n",
    "\n",
    "# Rename columns for consistency with Jordan's code\n",
    "res.rename(columns={'scores_01':'SES_01', 'scores_11':'SES_11'}, inplace=True)\n",
    "\n",
    "# Sanity check\n",
    "res.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below code computes other metrics for the LSOAs including:  SES ascent, SES percentile ascent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#  Compute rank of LSOA in 2001 (so low rank = 'low status')\n",
    "res['RANK_01'] = res.SES_01.rank(ascending=False)\n",
    "\n",
    "#  Compute rank of LSOA in 2011 (so low rank = 'low status')\n",
    "res['RANK_11'] = res.SES_11.rank(ascending=False)\n",
    "\n",
    "#  Compute amount by which LSOA has ascended (so +ve = status improvement; -ve = status decline)\n",
    "res.loc[:,'SES_ASC'] = res.loc[:,'SES_11'] - res.loc[:,'SES_01']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import re \n",
    "#  Calculate LSOA percentile score in 01\n",
    "res.loc[:,'SES_PR_01'] = res.RANK_01.rank(ascending=False, pct=True) * 100\n",
    "\n",
    "#  Calculate LSOA percentile score in 11\n",
    "res.loc[:,'SES_PR_11'] = res.RANK_11.rank(ascending=False, pct=True) * 100\n",
    "\n",
    "#  Calculate percentile change (so +ve = 'moved up' in the world; -ve = 'moved down')\n",
    "res.loc[:,'SES_PR_ASC'] = res.loc[:,'SES_PR_11'] - res.loc[:,'SES_PR_01']\n",
    "\n",
    "inp = res.loc[:,[x for x in res.columns if 'SES' not in x and 'RANK' not in x]]\n",
    "\n",
    "# Tidy up the naming\n",
    "inp.rename(columns=lambda x: re.sub('_11',' 2011',re.sub('_01',' 2001',x)), inplace=True)\n",
    "inp.rename(columns=lambda x: re.sub('kw_pct','Knowledge Worker Percentage',x), inplace=True)\n",
    "inp.rename(columns=lambda x: re.sub('he_pct','Highly-Educated Percentage',x), inplace=True)\n",
    "inp.rename(columns=lambda x: re.sub('hp','Property Prices (Transformed)',x), inplace=True)\n",
    "inp.rename(columns=lambda x: re.sub('hh','Household Income (Transformed)',x), inplace=True)\n",
    "\n",
    "# Save to file (note that we are also saving some info about the input variables as we use these as well)\n",
    "res[\n",
    "    ['RANK_01','RANK_11','SES_01','SES_11','SES_ASC','SES_PR_01','SES_PR_11','SES_PR_ASC']\n",
    "].to_csv(os.path.join(analytical,to_use + '-Scores.csv.gz'), compression='gzip', index=True) \n",
    "inp[\n",
    "    [x for x in inp.columns if '2001' in x]\n",
    "].to_csv(os.path.join(scores,to_use + '-Inputs-2001.csv.gz'), compression='gzip', index=True)\n",
    "inp[\n",
    "    [x for x in inp.columns if '2011' in x]\n",
    "].to_csv(os.path.join(scores,to_use + '-Inputs-2011.csv.gz'), compression='gzip', index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diagnostics\n",
    "\n",
    "If you are comfortable with the output of the code above you do not need to run the blocks below as these are simply sanity checks to help you (me) envision the output effectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Sanity check\n",
    "res[['SES_01','SES_11','RANK_01','RANK_11','SES_PR_01','SES_PR_11','SES_PR_ASC']].sample(5, random_state=r_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The lowest-ranked (highest status) LSOAs\n",
    "res.loc[res['RANK_01'] < 5,:].sort_values('RANK_01')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# The highest-ranked (lowest status) LSOAs\n",
    "res.loc[res['RANK_01'] > (res.RANK_01.max()-5),:].sort_values('RANK_01')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Biggest falls in percentile status\n",
    "res.sort_values('SES_PR_ASC').head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Biggest gains in percentile status\n",
    "res.sort_values('SES_PR_ASC', ascending=False).head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "g = sns.jointplot(x='SES_01', y='SES_11', data=res, kind='scatter', s=3, color='k', height=7, ratio=5, space=0, linewidth=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "g = sns.jointplot(x='SES_PR_01', y='SES_PR_11', data=res, kind='scatter', s=3, color='k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "g = sns.jointplot(x='RANK_01', y='RANK_11', data=res, kind='scatter', s=3, color='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automated Mapping\n",
    "\n",
    "With a little tweaking (as a result of changes I made and didn't propagate) you could automatically map the results in Python directly; however, I found that the overall results were rather better in QGIS and so haven't updated [notebook 10](10-Mapping Scores.ipynb)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "ML Gentrification 3",
   "language": "python",
   "name": "mlgent"
  },
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}