08-Neighbourhood Prediction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neighbourhood Prediction\n",
    "\n",
    "**You should almost certainly run the [Script](08-Neighbourhood Prediction.py) instead since I cannot guarantee that the Jupyter server will not timeout after a period of seeming inactivity and cause potential data loss.**\n",
    "\n",
    "However, this shows in a slightly more accessible form the same content as appears in the script so you are welcome to use this for exploratory purposes provided that you understand the likely impact of attempting to run the full GridSearch that is at the heart of this analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import re\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "import sklearn\n",
    "print('Your scikit-learn version is {}.'.format(sklearn.__version__))\n",
    "print('Please check it is at least 0.18.0.')\n",
    "\n",
    "from sklearn.preprocessing import scale\n",
    "from sklearn import linear_model\n",
    "from sklearn import tree\n",
    "from sklearn import preprocessing\n",
    "from sklearn import feature_selection\n",
    "from sklearn import model_selection\n",
    "from sklearn import metrics  \n",
    "from sklearn import ensemble\n",
    "\n",
    "from sklearn.externals.six import StringIO\n",
    "#from sklearn.model_selection import GridSearchCV\n",
    "#from sklearn.feature_selection import SelectKBest \n",
    "#from sklearn.feature_selection import f_regression\n",
    "\n",
    "from timeit import default_timer as timer\n",
    "import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "analytical = os.path.join('data','analytical')\n",
    "output     = os.path.join(os.path.expanduser('~'),'Documents','Dropbox','ESRC Gentrification','data','analytical')\n",
    "\n",
    "def load_status_scores(dtype):\n",
    "    status = pd.read_csv(os.path.join(analytical,dtype+'-Scores.csv.gz'), index_col=0)  # SES scores\n",
    "    \n",
    "    # Scores\n",
    "    status.drop(['RANK_01','RANK_11'], axis=1, inplace=True)\n",
    "    status.rename(columns={\n",
    "        'SES_01':'SES 2001',\n",
    "        'SES_11':'SES 2011',\n",
    "        'SES_ASC':'SES Ascent 2001-2011',\n",
    "        'SES_PR_01':'SES 2001 Percentile', # 99 = High-status\n",
    "        'SES_PR_11':'SES 2011 Percentile', # 99 = High-status\n",
    "        'SES_PR_ASC':'SES Percentile Ascent 2001-2011'\n",
    "    }, inplace=True)\n",
    "    return status\n",
    "\n",
    "def load_predictors(dtype):\n",
    "    \n",
    "    return status\n",
    "\n",
    "def classifier_report(clf, y_true, y_hat):\n",
    "    \n",
    "    txt = ''\n",
    "    \n",
    "    # If the task is regression evaluate using regression metrics, \n",
    "    # otherwise evaluate using classification metrics\n",
    "    txt += \"R2:        {0:8.5f}\".format(metrics.r2_score(y_true, y_hat)) + \"\\n\" #  R2 - Coefficient of determination\n",
    "    txt += \"MSE:       {0:8.5f}\".format(metrics.mean_squared_error(y_true, y_hat)) + \"\\n\"  #  Mean squared error regression loss\n",
    "    txt += \"MAE:       {0:8.5f}\".format(metrics.mean_absolute_error(y_true, y_hat)) + \"\\n\"  #  Mean absolute error regression loss\n",
    "    txt += \"Expl. Var: {0:8.5f}\".format(metrics.explained_variance_score(y_true, y_hat)) + \"\\n\"  # Explained variance regression score function\n",
    "    txt += \"\\n\"\n",
    "    \n",
    "    #print(metrics.accuracy_score(y_true, y_pred))  #  Accuracy Score\n",
    "    #print(metrics.classification_report(y_true, y_pred, target_names=[\"Unascended\",\"Ascended\"]))  #  Classification Report\n",
    "    #print(metrics.confusion_matrix(y_true, y_pred))  #  Confusion Matrix\n",
    "    #print()\n",
    "    return txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Information About Variable Utility (Optional)\n",
    "\n",
    "The code below evaluates the significance of each variable using the F-regression function in Scikit-Learn, and then sorts the results in ascending order.  The results are merged with data from the variable database.\n",
    "\n",
    "We don't actually use the output of this next step to perform feature selection as the Random Forest will take care of that for us. This is simply a way of understanding the relative utility of different variables to linear-type models (of which the RF is _not_ one)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "to_use = 'Untransformed'\n",
    "\n",
    "SES = load_status_scores(to_use)  # SES scores in 2011\n",
    "\n",
    "#  Read the transformed data\n",
    "d01_trs2 = pd.read_csv(os.path.join(analytical,to_use+'-2001-Data-Transformed_and_Scaled.csv.gz'), index_col=0)\n",
    "d11_trs2 = pd.read_csv(os.path.join(analytical,to_use+'-2011-Data-Transformed_and_Scaled.csv.gz'), index_col=0)\n",
    "\n",
    "# Data about variables used later in process\n",
    "vardb = pd.read_csv(os.path.join('data','variables.csv'), index_col=False)\n",
    "vardb.drop('Description', axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s1 = set(vardb.Predictor.values)\n",
    "s2 = set(d01_trs2.columns.values)\n",
    "if s2.difference(s1):\n",
    "    print(s2.difference(s1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kb = feature_selection.SelectKBest(feature_selection.f_regression, k='all')             #  Evaluate f-regression to evaluate all variables\n",
    "kb.fit(d01_trs2, SES.loc[:,'SES Ascent 2001-2011']) #  Pass variable data from 2001 to find correlation with SES 11\n",
    "\n",
    "# Check this!!!\n",
    "print(\"Max f-test value: \" + str(np.max(kb.scores_)))\n",
    "f_test = kb.scores_\n",
    "f_test /= np.max(f_test) # Normalise by maximum value (http://scikit-learn.org/stable/auto_examples/feature_selection/plot_f_test_vs_mi.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Format results and write to file\n",
    "results = pd.DataFrame(data=f_test, index=d01_trs2.columns)\n",
    "results.reset_index(inplace=True)\n",
    "results.columns = ['Predictor','Score']\n",
    "results.to_csv(os.path.join(analytical,to_use+'-Variable Results.csv'), index=False)\n",
    "\n",
    "# Formatted results\n",
    "fresults = vardb.loc[:,['Predictor','Title','Category','Group']].merge(results, on='Predictor', how='left')\n",
    "fresults[['Title','Category','Score','Group']].sort_values(by='Score', ascending=False)\n",
    "fresults.to_csv(os.path.join(analytical,to_use+'-Variable Importance.csv'), index=False)\n",
    "\n",
    "# Results are...\n",
    "fresults.sort_values(by='Score', ascending=False).head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Developing the Model\n",
    "\n",
    "The code below is concerned with building the best predictive model for the period of 2001-2011."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Can override to_use here\n",
    "to_use = 'Untransformed'\n",
    "\n",
    "SES = load_status_scores(to_use)  # SES scores in 2011\n",
    "\n",
    "#  Read the transformed data\n",
    "d01_trs2 = pd.read_csv(os.path.join(analytical,to_use+'-2001-Data-Transformed_and_Scaled.csv.gz'), index_col=0)\n",
    "d11_trs2 = pd.read_csv(os.path.join(analytical,to_use+'-2011-Data-Transformed_and_Scaled.csv.gz'), index_col=0)\n",
    "\n",
    "# Data about variables used later in process\n",
    "vardb = pd.read_csv(os.path.join('data','variables.csv'), index_col=False)\n",
    "vardb.drop('Description', axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To evaluate the models most reliably a portion of the dataset must be kept as holdout to evaluate the classifier on independently.  The code below splits the data into training and test sets using a test size of 20%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_full = d01_trs2\n",
    "\n",
    "X_train, X_test, y_train, y_test = model_selection.train_test_split(\n",
    "    d01_trs2, SES['SES Ascent 2001-2011'], test_size=0.2, random_state=r_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Start Logging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "log = open(os.path.join(to_use+'-Fit.txt'),'w')\n",
    "print(\"Data Transform: \" + to_use, file=log)\n",
    "print(\"\", file=log)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple Linear Regression\n",
    "\n",
    "To begin with modelling was attempted using 1R (i.e. 1 rule) modelling.  This code below uses Stochastic Loss Gradient to build a simple linear regression estimator using each variable separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Evaluate each predictor using simple linear regression\n",
    "preds_ls = list(d01_trs2.columns)  #  List of predictors\n",
    "df = pd.DataFrame(X_train, columns=d01_trs2.columns)\n",
    "results_dict = dict()  #  Store results here\n",
    "\n",
    "#  Loop over each predictor and evaluate it's performance\n",
    "for p in preds_ls:\n",
    "    #print(\"Evaluating: \" + p)\n",
    "    #X_sing = df[p].as_matrix()  #  Get predictor data from dataset   \n",
    "    clf = linear_model.SGDRegressor(loss='squared_loss', penalty=None, random_state=r_state, max_iter=1000, tol=1e-3) #  Build Stochastic Gradient Descent estimator\n",
    "    clf.fit(X_train[[p]],y_train)\n",
    "    y_pred = clf.predict(X_test[[p]])\n",
    "    sc  = metrics.r2_score(y_test, y_pred, multioutput='variance_weighted')\n",
    "    mse = metrics.mean_squared_error(y_test, y_pred)  #  Mean squared error regression loss\n",
    "    mae = metrics.mean_absolute_error(y_test, y_pred)  #  Mean absolute error regression loss\n",
    "    var = metrics.explained_variance_score(y_test, y_pred)  # Explained variance regression score function\n",
    "    \n",
    "    results_dict[p] = [sc, mse, mae, var] \n",
    "\n",
    "results = pd.DataFrame.from_dict(results_dict, orient='index').sort_values(by=0, ascending=False)\n",
    "results.reset_index(inplace=True)\n",
    "results.columns = ['Predictor','R2-Score','MSE','MAE','Explained Variance']\n",
    "df = vardb.loc[:,['Predictor','Category']].merge(results, on='Predictor', how='left' )\n",
    "df.sort_values(by=['R2-Score'], ascending=False, inplace=True)\n",
    "df.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# And for logging purposes\n",
    "p   = df['Predictor'].iloc[0]\n",
    "print(p)\n",
    "clf = linear_model.SGDRegressor(loss='squared_loss', penalty=None, random_state=r_state, max_iter=1000, tol=1e-3)\n",
    "clf.fit(X_train[[p]],y_train)\n",
    "y_pred = clf.predict(X_test[[p]])\n",
    "\n",
    "print(\"Singular Regression results:\")\n",
    "print(classifier_report(clf, y_test, y_pred), file=log)\n",
    "print(classifier_report(clf, y_test, y_pred)) # clf, y_test, y_hat\n",
    "print(\"\", file=log)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remember that these are the results on the training data set using a fixed random seed and not the full data set, so changes to the seed/split will certainly change the results!**\n",
    "\n",
    "For the Untransformed data set we get the following results:\n",
    "\n",
    "| Predictor | Category | R2-Score | MSE | MAE | Explained Variance |\n",
    "| --------- | -------- | -------- | ---- |---- | ------------- |\n",
    "| House Prices | Scoring Metric | 0.54164 | 0.28576 | 0.33353 | 0.54405 |\n",
    "\n",
    "For the Box-Cox transformed data set we get the following:\n",
    "\n",
    "| Predictor | Category | R2-Score | MSE | MAE | Explained Variance |\n",
    "| --------- | -------- | -------- | ---- |---- | ------------- |\n",
    "| House Prices (Box-Cox Transformed) | Scoring Metric | 0.186940 | 0.075678 | 0.199408 | 0.198013 |\n",
    "\n",
    "And for the Log-transformed data we get: \n",
    "\n",
    "| Predictor | Category | R2-Score | MSE | MAE | Explained Variance |\n",
    "| --------- | -------- | -------- | ---- |---- | ------------- |\n",
    "| G. Wholesale and retail | Industry of Employment | 0.087389 | 0.091154 | 0.219888 | 0.098506 |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiple Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "clf = linear_model.LinearRegression(fit_intercept=True, copy_X=True)\n",
    "clf.fit(X_train,y_train)\n",
    "y_pred = clf.predict(X_test)\n",
    "\n",
    "print(\"Multiple Regression results:\", file=log)\n",
    "print(classifier_report(clf, y_test, y_pred), file=log)\n",
    "print(classifier_report(clf, y_test, y_pred))\n",
    "print(\"\", file=log)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For Multiple Linear Regression I get:\n",
    "\n",
    "| Model | R2-Score | MSE | MAE | Explained Variance |\n",
    "| ----- | -------- | ---- |---- | ------------- |\n",
    "| Singular Regression | 0.54164 | 0.28576 | 0.33353 | 0.54405 |\n",
    "| Multiple-Regression | 0.63932 | 0.22486 | 0.30493 | 0.64028 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predicting the 'Future' (2001 > 2011)\n",
    "\n",
    "The code below trains the model on training sets and then predicts the entire results of 2011."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Baseline without Hyperparameter Tuning\n",
    "\n",
    "Accepting only the default parameters for the models so that we have a baseline before tuning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extemely Random Trees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "clf = ensemble.ExtraTreesRegressor(n_jobs=-1, random_state=r_state, n_estimators=100)  \n",
    "clf.fit(X_train,y_train)\n",
    "y_pred = clf.predict(X_test)\n",
    "\n",
    "print(\"Default Extra Trees results:\", file=log)\n",
    "print(classifier_report(clf, y_test, y_pred), file=log)\n",
    "print(classifier_report(clf, y_test, y_pred))\n",
    "print(\"\", file=log)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Single best-performing option from Notebook 7\n",
    "clf = ensemble.ExtraTreesRegressor(n_estimators=180, n_jobs=-1, random_state=r_state)  \n",
    "clf.fit(X_train,y_train)\n",
    "y_pred = clf.predict(X_test)\n",
    "\n",
    "print(\"Individually Tuned Extra Trees results:\")\n",
    "print(classifier_report(clf, y_test, y_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "log.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the two Extra-Tress Regressors I get:\n",
    "\n",
    "| Model | R2-Score | MSE | MAE | Explained Variance |\n",
    "| ----- | -------- | ---- |---- | ------------- |\n",
    "| Singular Regression | 0.54164 | 0.28576 | 0.33353 | 0.54405 |\n",
    "| Multiple-Regression | 0.63932 | 0.22486 | 0.30493 | 0.64028 |\n",
    "| Un-Tuned RF | 0.66722 | 0.20747 | 0.28518 | 0.67187 |\n",
    "| Single-Tuned RF| 0.69477 | 0.19029 | 0.26110 | 0.69934 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GridSearchCV (Run from Script)\n",
    "\n",
    "Using a grid search to tune the hyperparameters. **This should actually be done using the included script of the same name since Jupyter may timeout while running the GridSearch (or, at least, it did for me).**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extremely Random Trees\n",
    "\n",
    "Similar to RandomForests we can tuen the following hyperparameters:\n",
    "1. `n_estimators` (number of trees; difficult to overfit so large _n_ probably simplest starting point).\n",
    "2. `max_depth` (maximum depth of trees; can encourage overfitting since more depth == more complexity).\n",
    "3. `max_features` (maximum number of features to consider at each split; allows more complex models so may lead to overfitting).\n",
    "4. `min_samples_leaf` (also helps to control depth and reduce overfitting by preventing splits that hold outliers).\n",
    "\n",
    "Note that the permutations can pile up rather quickly when we incorporate additional parameters such as feature determination and bootstrapping with cross-validation. On a 2.9GHz Core i5 Mac with 16GB of RAM this is working out at about 16s per fold, so 100 candidates with 7 folds == 700 fits == 10,500s == 175m == 3hrs.\n",
    "\n",
    "**Broadly: 100 fits ~= 45 minutes on a laptop.**\n",
    "\n",
    "<span style=\"color:red;weight:bold\">I ended up moving the code below to a script since I was experiencing timeouts on Jupyter.</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Use a grid over parameters of interest -- search grid\n",
    "# partly extracted from testing with notebook 7 and party\n",
    "# from playing with grid ranges here (since results produced\n",
    "# by manipulating one parameter separately from the others \n",
    "# don't always replicate well as the single tuned parameter \n",
    "# for the ensemble as a whole). In other words, just because\n",
    "# max_depth==10 was the best result from manipulating _only_\n",
    "# tree depth doesn't mean that it will be the best when you\n",
    "# start manipulating all the main hyperparameters together.\n",
    "param_grid = {\n",
    "    \"n_estimators\"      : [int(x) for x in np.arange(start=160, stop=211, step=20)] +  \n",
    "                           [int(x) for x in np.arange(start=1300, stop=1501, step=100)] +\n",
    "                           [int(x) for x in np.arange(start=1800, stop=2001, step=100)],\n",
    "    \"max_depth\"         : [None], # [int(x) for x in np.arange(start=10, stop=141, step=90)]+[None],\n",
    "    \"min_samples_leaf\"  : [1,2,4], #\n",
    "    \"max_features\"      : [None] # [0.7, 0.85, None], # For regression normally n_features (i.e. auto)\n",
    "}\n",
    "\n",
    "print(\"Estimators: \" + str(param_grid['n_estimators']))\n",
    "print(\"Depth: \" + str(param_grid['max_depth']))\n",
    "print(\"Minimum Samples Leaf: \" + str(param_grid['min_samples_leaf']))\n",
    "print(\"Maximum Features: \" + str(param_grid['max_features']))\n",
    "print(\"Number of permutations: \" + str(len(param_grid['n_estimators']) * len(param_grid['max_depth']) * len(param_grid['max_features']) * len(param_grid['min_samples_leaf'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "clf = ensemble.ExtraTreesRegressor(n_jobs=-1, random_state=r_state) # Can be 'mae' or 'mse' -- should presumably match scoring below\n",
    "start = timer()\n",
    "# There is some disagreement about whether cross-validation or bootstrapping \n",
    "# is needed for ExtraTrees (or even RandomForests) regressors:\n",
    "# https://stats.stackexchange.com/questions/279163/cross-validation-in-extratreesregressor\n",
    "scoring = {'mae':'neg_mean_absolute_error', 'mse':'neg_mean_squared_error'} #, 'r2':'r2'}\n",
    "cv = model_selection.GridSearchCV(estimator=clf, param_grid=param_grid, cv=4, n_jobs=6, verbose=0, scoring='neg_mean_squared_error')\n",
    "cv.fit(X_train, y_train)\n",
    "duration = timer() - start\n",
    "print(\"Execution complete in: {0:15.1f}s\".format(duration) + \" (\" + str(datetime.timedelta(seconds=duration)) + \")\")\n",
    "print(\"Best score: \" + str(cv.best_score_))\n",
    "print(\"Done.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"Best score: \" + str(cv.best_score_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "log = open(os.path.join(output,to_use+'-Fit.txt'),'a')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"Params: \", file=log)\n",
    "print(param_grid, file=log)\n",
    "print(\"Best Cross-Validation score: \" + str(cv.best_score_), file=log)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "best_clf = cv.best_estimator_ # Extract the best estimator from the GridSearch\n",
    "best_clf.fit(X_train, y_train)\n",
    "y_pred  = best_clf.predict(X_test)\n",
    "\n",
    "print(\"Best parameters from Cross-Validation: \" + str(cv.best_params_), file=log)\n",
    "print(\"Best parameters from Cross-Validation: \" + str(cv.best_params_))\n",
    "print(\"\", file=log)\n",
    "\n",
    "print(\"Cross-check against full spec of model: \", file=log)\n",
    "print(best_clf.get_params, file=log)\n",
    "print(best_clf.get_params)\n",
    "print(\"\", file=log)\n",
    "\n",
    "print(\"Tuned Extra Trees result:\", file=log)\n",
    "print(classifier_report(best_clf, y_test, y_pred), file=log)\n",
    "print(classifier_report(best_clf, y_test, y_pred))\n",
    "print(\"\", file=log)\n",
    "\n",
    "# Create a data frame of feature importance so that we\n",
    "# can inspect later...\n",
    "fi = pd.DataFrame.from_dict({'feature':X_test.columns.values, 'importance':best_clf.feature_importances_})\n",
    "fi.sort_values(by='importance', ascending=False, inplace=True)\n",
    "fi.to_csv(os.path.join(analytical,to_use+'-Feature Importance.csv.gz'), compression='gzip', index=False)\n",
    "\n",
    "print(\"Feature Importances (5 Biggest):\", file=log)\n",
    "print(fi.head(5), file=log)\n",
    "print(fi.head(5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "log.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outputting Final Results \n",
    "\n",
    "Best performing model from testing across grid:\n",
    "```\n",
    "Cross-check against full spec of model: \n",
    "<bound method BaseEstimator.get_params of ExtraTreesRegressor(bootstrap=False, criterion='mse', max_depth=None,\n",
    "          max_features=0.85, max_leaf_nodes=None,\n",
    "          min_impurity_decrease=0.0, min_impurity_split=None,\n",
    "          min_samples_leaf=2, min_samples_split=2,\n",
    "          min_weight_fraction_leaf=0.0, n_estimators=1400, n_jobs=-1,\n",
    "          oob_score=False, random_state=42, verbose=0, warm_start=False)>\n",
    "\n",
    "Tuned Extra Trees result:\n",
    "R2:         0.69899\n",
    "MSE:        0.18766\n",
    "MAE:        0.25969\n",
    "Expl. Var:  0.70261\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Can override to_use here for other transformations\n",
    "to_use = 'Untransformed'\n",
    "\n",
    "SES = load_status_scores(to_use)  # SES scores in 2011\n",
    "\n",
    "#  Read the transformed data\n",
    "d01_trs2 = pd.read_csv(os.path.join(analytical,to_use+'-2001-Data-Transformed_and_Scaled.csv.gz'), index_col=0)\n",
    "d11_trs2 = pd.read_csv(os.path.join(analytical,to_use+'-2011-Data-Transformed_and_Scaled.csv.gz'), index_col=0)\n",
    "\n",
    "# Data about variables used later in process\n",
    "vardb = pd.read_csv(os.path.join('data','variables.csv'), index_col=False)\n",
    "vardb.drop('Description', axis=1, inplace=True)\n",
    "\n",
    "X_full = d01_trs2\n",
    "\n",
    "X_train, X_test, y_train, y_test = model_selection.train_test_split(\n",
    "    d01_trs2, SES['SES Ascent 2001-2011'], test_size=0.2, random_state=r_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "best_clf = ensemble.ExtraTreesRegressor(bootstrap=False, criterion='mse', max_depth=None,\n",
    "          max_features=0.85, max_leaf_nodes=None,\n",
    "          min_impurity_decrease=0.0, min_impurity_split=None,\n",
    "          min_samples_leaf=2, min_samples_split=2,\n",
    "          min_weight_fraction_leaf=0.0, n_estimators=1400, n_jobs=-1,\n",
    "          oob_score=False, random_state=42, verbose=0, warm_start=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "best_clf.fit(X_train, y_train)\n",
    "y_pred  = best_clf.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"Cross-check against full spec of model: \")\n",
    "print(best_clf.get_params)\n",
    "print(\"\")\n",
    "\n",
    "print(\"Tuned Extra Trees result:\")\n",
    "print(classifier_report(best_clf, y_test, y_pred))\n",
    "print(\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the two Extra-Tress Regressors I get:\n",
    "\n",
    "| Model | R2-Score | MSE | MAE | Explained Variance |\n",
    "| ----- | -------- | ---- |---- | ------------- |\n",
    "| Singular Regression | 0.54164 | 0.28576 | 0.33353 | 0.54405 |\n",
    "| Multiple-Regression | 0.63932 | 0.22486 | 0.30493 | 0.64028 |\n",
    "| Un-Tuned RF | 0.66722 | 0.20747 | 0.28518 | 0.67187 |\n",
    "| Single-Tuned RF| 0.69477 | 0.19029 | 0.26110 | 0.69934 |\n",
    "| Fully-Tuned RF | 0.69739 | 0.18866 | 0.26012 | 0.70101 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Output Graph Representation of 1 Tree\n",
    "\n",
    "This is used in the article to illustrate how a decision tree within the Random Forest works to split the data so as to make predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.tree import export_graphviz\n",
    "\n",
    "t = best_clf.estimators_[0]\n",
    "feature_names = X_test.columns.values\n",
    "export_graphviz(t, out_file=os.path.join(analytical,to_use + \"-tree.dot\"), filled=True, rounded=True, feature_names=feature_names)\n",
    "os.system('dot -Tpng ' + os.path.join(analytical,to_use + \"-tree.dot\") + ' -o ' + os.path.join(analytical,to_use + \"-tree.png\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Save Feature Importances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create a data frame of feature importance so that we\n",
    "# can inspect later...\n",
    "fi = pd.DataFrame.from_dict({'feature':X_test.columns.values, 'importance':best_clf.feature_importances_})\n",
    "\n",
    "fi = vardb.loc[:,['Predictor','Category']].merge(\n",
    "    pd.DataFrame.from_dict({'feature':X_test.columns.values, 'importance':best_clf.feature_importances_}), \n",
    "    left_on='Predictor', right_on='feature', how='left' )\n",
    "fi.drop(['feature'], axis=1, inplace=True)\n",
    "fi.sort_values(by='importance', ascending=False, inplace=True)\n",
    "fi.to_csv(os.path.join(analytical,to_use+'-Feature_Importance.csv.gz'), compression='gzip', index=False)\n",
    "\n",
    "print(\"Feature Importances (5 Biggest):\")\n",
    "print(fi.head(5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SES = load_status_scores(to_use) \n",
    "y_pr = best_clf.predict(X_full)\n",
    "\n",
    "predicted11 = pd.DataFrame(\n",
    "    {'lsoacd':        pd.Series(d01_trs2.index), \n",
    "     'SES Ascent 2001-2011 (Predicted)': pd.Series(y_pr)})  #  Combine with list of areas\n",
    "predicted11.set_index('lsoacd', inplace=True)\n",
    "predicted11.sample(3, random_state=r_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "predicted11 = predicted11.merge(SES, left_index=True, right_index=True, how='inner')\n",
    "\n",
    "predicted11['SES 2011 (Predicted)'] = predicted11.loc[:,'SES 2001'] \\\n",
    "                                      + predicted11.loc[:,'SES Ascent 2001-2011 (Predicted)']\n",
    "predicted11['Score Divergence'] = predicted11.loc[:,'SES 2011 (Predicted)'] \\\n",
    "                                   - predicted11.loc[:,'SES 2011']\n",
    "predicted11['Ascent Divergence'] = predicted11.loc[:,'SES Ascent 2001-2011 (Predicted)'] \\\n",
    "                                   - predicted11.loc[:,'SES Ascent 2001-2011']\n",
    "\n",
    "predicted11.sort_index(axis=1, inplace=True)\n",
    "\n",
    "predicted11.to_csv(os.path.join(analytical,to_use+'-Predicted Ascent 2001-2011.csv.gz'), compression='gzip', index=True)\n",
    "\n",
    "# Sanity check\n",
    "print(\"Results data frame has \" + str(predicted11.shape[0]) + \" rows.\")\n",
    "predicted11.sample(5, random_state=r_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure('SES Divergence')\n",
    "sns.distplot(predicted11['Score Divergence'], kde=True)      \n",
    "fig = plt.gcf() # *G*et the *C*urrent *F*igure environment so that the next command works\n",
    "plt.savefig(\"{0}-{1}.pdf\".format(to_use, 'SES Ascent 2001-2011-Divergence'), bbox_inches=\"tight\")\n",
    "plt.close()\n",
    "print(\"Done.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure('SES Divergence Scatter')\n",
    "sns.jointplot(x='SES Ascent 2001-2011', y='SES Ascent 2001-2011 (Predicted)', data=predicted11, kind='scatter')     \n",
    "fig = plt.gcf() # *G*et the *C*urrent *F*igure environment so that the next command works\n",
    "plt.savefig(\"{0}-{1}.pdf\".format(to_use, 'SES Ascent 2001-2011-Divergence (Scatter)'), bbox_inches=\"tight\")\n",
    "plt.close()\n",
    "print(\"Done.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predicting the _Future_ (2011 > 2021)\n",
    "\n",
    "The code below this is used to make predictions of 2021.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Make future predictions\n",
    "y_pred_21 = best_clf.predict(d11_trs2)  #  Make predictions using data from 2011\n",
    "predicted21 = pd.DataFrame({\n",
    "        'lsoacd': pd.Series(d11_trs2.index),\n",
    "        'SES Ascent 2011-2021 (Predicted)': pd.Series(y_pred_21)})\n",
    "predicted21.set_index('lsoacd', inplace=True)\n",
    "\n",
    "predicted21.to_csv(os.path.join(analytical,to_use+'-Predicted Ascent 2011-2021.csv.gz'), compression='gzip', index=True)  #  Write results to csv\n",
    "predicted21.sample(3, random_state=r_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pdf = predicted11.merge(predicted21, left_index=True, right_index=True, how='left')  #  Integrate SES 2021 predictions into SES score data\n",
    "pdf['SES 2021 (Predicted)'] = pdf.loc[:,'SES 2011'] + pdf.loc[:,'SES Ascent 2011-2021 (Predicted)']  # Compute SES score in 2021\n",
    "\n",
    "#  Compute rank in 2012\n",
    "pdf['SES 2021 Percentile'] = pdf.loc[:,'SES 2021 (Predicted)'].rank(ascending=True, pct=True)*100\n",
    "\n",
    "#  Compute change in LSOA ranking from 2011 to 2021\n",
    "pdf['SES Percentile Ascent 2011-2021'] = pdf.loc[:,'SES 2021 Percentile'] \\\n",
    "                                          - pdf.loc[:,'SES 2011 Percentile']\n",
    "\n",
    "pdf[['SES 2001','SES 2011','SES 2021 (Predicted)','SES Ascent 2001-2011','SES Ascent 2011-2021 (Predicted)']].sample(3, random_state=r_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pdf.columns.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Write results to file\n",
    "cols = ['SES 2001','SES 2011 (Predicted)','SES 2011','SES 2021 (Predicted)',\n",
    "        'SES 2001 Percentile','SES 2011 Percentile','SES 2021 Percentile',\n",
    "        'SES Ascent 2001-2011','SES Ascent 2001-2011 (Predicted)','SES Ascent 2011-2021 (Predicted)',\n",
    "        'SES Percentile Ascent 2001-2011','SES Percentile Ascent 2011-2021',\n",
    "        'Score Divergence','Ascent Divergence']\n",
    "\n",
    "pdf = pdf[cols]\n",
    "fcols = ['Score Divergence','Ascent Divergence']\n",
    "pdf.loc[:, fcols] = pdf[fcols].astype(float).applymap('{0:.15f}'.format)\n",
    "pdf.to_csv(os.path.join(analytical,to_use+'-Predictions.csv.gz'), compression='gzip', index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}