notebooks/xgboost_predictor.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.model_selection import cross_validate\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.metrics.scorer import make_scorer\n",
    "import xgboost as xgb\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import sys\n",
    "sys.path.append('..')\n",
    "from utils import preprocess, missing, evaluate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Globals and load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "411 indicators included\n"
     ]
    }
   ],
   "source": [
    "target = 'SI.POV.DDAY'\n",
    "predict_year=2010\n",
    "#percent of input Indicators to use (set to 100 for full set of input features)\n",
    "percent = 50\n",
    "\n",
    "#Load the data from disk\n",
    "input_dir = '.\\\\..\\\\data\\\\'\n",
    "data_input = \"cleaned_data.pkl\"\n",
    "data = pd.read_pickle(input_dir + data_input)\n",
    "\n",
    "#Possible subset of data choosen to reduce calulation time\n",
    "#For percetages less than 100% we try to choose a subset that represents the spread of variables\n",
    "\n",
    "if percent == 100:\n",
    "    pass\n",
    "else: \n",
    "    num_indicators_original = data.shape[1]\n",
    "    step = int(100/percent)\n",
    "    data_new = data.iloc[:,::step].copy()\n",
    "    #Add the target column if not already included\n",
    "    if target not in data_new.columns:\n",
    "        data_new[target] = data[target]\n",
    "    data = data_new\n",
    "    \n",
    "print(data.shape[1], \"indicators included\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Window data and preprocess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 3min 26s\n"
     ]
    }
   ],
   "source": [
    "%time data_regressors, data_targets = \\\n",
    "        preprocess.window_data(data, lag=3,num_windows=10, step=1, predict_year=2010, \\\n",
    "                         target=target, impute_type='interpolation')\n",
    "\n",
    "#Break up into training and testing data.\n",
    "\n",
    "idx = pd.IndexSlice\n",
    "data_train_regressors = data_regressors.loc[idx[:,2:10],:]\n",
    "data_train_targets = data_targets.loc[idx[:,2:10],:]\n",
    "data_test_regressors = data_regressors.loc[idx[:,1],:]\n",
    "data_test_targets= data_targets.loc[idx[:,1],:]\n",
    "\n",
    "#For Training, only consider windows that don't have a missing target as they offer nothing to training\n",
    "#Therefore, remove those observations from both the training regressors and targets datasets.\n",
    "data_train_regressors_subset = data_train_regressors[~np.isnan(list(data_train_targets.values.flatten()))]\n",
    "data_train_targets_subset = data_train_targets[~np.isnan(list(data_train_targets.values.flatten()))]\n",
    "\n",
    "#For testing, also remove windows with no target variable as it is impossible to measure preformance.\n",
    "data_test_regressors_subset = data_test_regressors[~np.isnan(list(data_test_targets.values.flatten()))]\n",
    "data_test_targets_subset = data_test_targets[~np.isnan(list(data_test_targets.values.flatten()))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = data_train_regressors_subset.values\n",
    "y_train = data_train_targets_subset.values.ravel()\n",
    "X_test = data_test_regressors_subset.values\n",
    "y_test = data_test_targets_subset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### XGBoost Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Out-of-the-box using the Scikit-learn interface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE of XGBoost out-of-the-box is: 5.154362890679558\n"
     ]
    }
   ],
   "source": [
    "XGB = xgb.XGBRegressor(random_state=42 ,objective='reg:squarederror', subsample=0.9)\n",
    "XGB.fit( X_train,y_train)\n",
    "#Make predictions\n",
    "predictions = XGB.predict(X_test) \n",
    "\n",
    "mse= mean_squared_error(y_test, predictions)\n",
    "print(\"RMSE of XGBoost out-of-the-box is:\", np.sqrt(mse))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tuning of the Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_folds = 5\n",
    "\n",
    "scorer = make_scorer(mean_squared_error ,greater_is_better=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 1. Tune the number of estimators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal number of estimators: 100\n"
     ]
    }
   ],
   "source": [
    "model = xgb.XGBRegressor(random_state=42,\n",
    "                         objective='reg:squarederror',\n",
    "                         max_depth=5, \n",
    "                         min_child_weight = 1, \n",
    "                         gamma = 0, \n",
    "                         subsample=0.9, \n",
    "                         colsample_bytree = 0.8, \n",
    "                         scale_pos_weight = 1)\n",
    "\n",
    "param = model.get_xgb_params()\n",
    "data_matrix = xgb.DMatrix(X_train, label=y_train)\n",
    "cvresult = xgb.cv(param, data_matrix, num_boost_round=model.get_params()['n_estimators'], nfold=cv_folds,\n",
    "            metrics='rmse', early_stopping_rounds=50)\n",
    "#Set the optimised number of estimators\n",
    "model.set_params(n_estimators=cvresult.shape[0])\n",
    "print(\"Optimal number of estimators:\", cvresult.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE of xgboost after tuning (step 1) is: 5.189437244183729\n"
     ]
    }
   ],
   "source": [
    "model.fit( X_train,y_train)\n",
    "#Make predictions\n",
    "predictions = model.predict(X_test) \n",
    "\n",
    "mse= mean_squared_error(y_test, predictions)\n",
    "print(\"RMSE of xgboost after tuning (step 1) is:\", np.sqrt(mse))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 2. Tune max_depth and min_child_weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE of xgboost after tuning (step 2) is: 5.494248333090531\n"
     ]
    }
   ],
   "source": [
    "params = {\n",
    " 'max_depth':range(3,10,2),\n",
    " 'min_child_weight':range(1,6,2)\n",
    "}\n",
    "\n",
    "grid_model = GridSearchCV(model, param_grid = params, scoring=scorer,\n",
    "                        n_jobs=4,iid=False, cv=5)\n",
    "grid_model.fit(X_train,y_train)\n",
    "\n",
    "#Score the best model using the test data\n",
    "model = grid_model.best_estimator_\n",
    "model.fit( X_train,y_train)\n",
    "#Make predictions\n",
    "predictions = model.predict(X_test) \n",
    "\n",
    "mse= mean_squared_error(y_test, predictions)\n",
    "print(\"RMSE of xgboost after tuning (step 2) is:\", np.sqrt(mse))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result of our test shows that the performance of the model after tuning was actually worse than before. The model is generalising very poorly. This may be a reflection on some kind of bias introduced in creating our training and data subsets. It would be worth looking at how I decided to discard any countries early on that did not have target values for the target year, 2010. It may have made more sense to window the data and split into training and test subsets and then, after, discard any observations that did not have a target value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'max_depth': 5, 'min_child_weight': 5}"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_model.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 3. Tune Gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE of xgboost after tuning (step 3) is: 5.141798455470761\n"
     ]
    }
   ],
   "source": [
    "params = {\n",
    " 'gamma':[i/10.0 for i in range(0,5)]\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "grid_model = GridSearchCV(model, param_grid = params, scoring=scorer,\n",
    "                        n_jobs=4,iid=False, cv=5)\n",
    "grid_model.fit(X_train,y_train)\n",
    "\n",
    "#Score the best model using the test data\n",
    "model = grid_model.best_estimator_\n",
    "model.fit( X_train,y_train)\n",
    "#Make predictions\n",
    "predictions = model.predict(X_test) \n",
    "\n",
    "mse= mean_squared_error(y_test, predictions)\n",
    "print(\"RMSE of xgboost after tuning (step 3) is:\", np.sqrt(mse))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'gamma': 0.0}"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_model.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 4. Tune Regularization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE of xgboost after tuning (step 4) is: 5.141798730456352\n"
     ]
    }
   ],
   "source": [
    "params = {\n",
    " 'reg_alpha':[1e-5, 1e-2, 0.1, 1]\n",
    "}\n",
    "\n",
    "grid_model = GridSearchCV(model, param_grid = params, scoring=scorer,\n",
    "                        n_jobs=4,iid=False, cv=5 ,return_train_score=True)\n",
    "grid_model.fit(X_train,y_train)\n",
    "\n",
    "#Score the best model using the test data\n",
    "model = grid_model.best_estimator_\n",
    "model.fit( X_train,y_train)\n",
    "#Make predictions\n",
    "predictions = model.predict(X_test) \n",
    "\n",
    "mse= mean_squared_error(y_test, predictions)\n",
    "print(\"RMSE of xgboost after tuning (step 4) is:\", np.sqrt(mse))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'reg_alpha': 1e-05}"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_model.best_params_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE of xgboost after tuning (step 5) is: 5.141798730456352\n"
     ]
    }
   ],
   "source": [
    "params = {\n",
    " 'reg_lambda':[1e-5, 1e-2, 0.1, 1]\n",
    "}\n",
    "\n",
    "grid_model = GridSearchCV(model, param_grid = params, scoring=scorer,\n",
    "                        n_jobs=4,iid=False, cv=5 ,return_train_score=True)\n",
    "grid_model.fit(X_train,y_train)\n",
    "\n",
    "#Score the best model using the test data\n",
    "model = grid_model.best_estimator_\n",
    "model.fit( X_train,y_train)\n",
    "#Make predictions\n",
    "predictions = model.predict(X_test) \n",
    "\n",
    "mse= mean_squared_error(y_test, predictions)\n",
    "print(\"RMSE of xgboost after tuning (step 5) is:\", np.sqrt(mse))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### XGBoost Model with no Imputation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 2.88 s\n"
     ]
    }
   ],
   "source": [
    "%time data_regressors, data_targets = \\\n",
    "        preprocess.window_data(data, lag=3,num_windows=10, step=1, predict_year=2010, \\\n",
    "                         target=target)\n",
    "\n",
    "#Break up into training and testing data.\n",
    "\n",
    "idx = pd.IndexSlice\n",
    "data_train_regressors = data_regressors.loc[idx[:,2:10],:]\n",
    "data_train_targets = data_targets.loc[idx[:,2:10],:]\n",
    "data_test_regressors = data_regressors.loc[idx[:,1],:]\n",
    "data_test_targets= data_targets.loc[idx[:,1],:]\n",
    "\n",
    "#For Training, only consider windows that don't have a missing target as they offer nothing to training\n",
    "#Therefore, remove those observations from both the training regressors and targets datasets.\n",
    "data_train_regressors_subset = data_train_regressors[~np.isnan(list(data_train_targets.values.flatten()))]\n",
    "data_train_targets_subset = data_train_targets[~np.isnan(list(data_train_targets.values.flatten()))]\n",
    "\n",
    "#For testing, also remove windows with no target variable as it is impossible to measure preformance.\n",
    "data_test_regressors_subset = data_test_regressors[~np.isnan(list(data_test_targets.values.flatten()))]\n",
    "data_test_targets_subset = data_test_targets[~np.isnan(list(data_test_targets.values.flatten()))]\n",
    "\n",
    "X_train_miss = data_train_regressors_subset.values\n",
    "y_train_miss = data_train_targets_subset.values.ravel()\n",
    "X_test_miss = data_test_regressors_subset.values\n",
    "y_test_miss = data_test_targets_subset\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Out-of-the-box xgboost on data without imputation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hand-tuning of the XGBoost model. I pick some important paramneters and play around until I get a good result. I'm sure there is more accuracy that can be obtained from this model by gridsearching but I think this is enough to illustrate that using a XGBoost (or perhaps any tree-based predictive algo) without any imputation done on the input data gives by far the best results of any of the models tried."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE of XGBoost out-of-the-box is: 3.0772104636217508\n"
     ]
    }
   ],
   "source": [
    "XGB = xgb.XGBRegressor(n_estimators=200,  objective='reg:squarederror',max_depth=7, subsample=0.87, reg_lambda=0.2)\n",
    "XGB.fit( X_train_miss,y_train_miss)\n",
    "#Make predictions\n",
    "predictions = XGB.predict(X_test_miss) \n",
    "\n",
    "mse= mean_squared_error(y_test_miss, predictions)\n",
    "print(\"RMSE of XGBoost out-of-the-box is:\", np.sqrt(mse))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}