Tutorial13 Model Stacking.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><img src=\"http://sydney.edu.au/images/content/about/logo-mono.jpg\"></center>\n",
    "\n",
    "<center><h1>Statistical Learning and Data Mining (QBUS6810)</h1></center>\n",
    "<center><h2>Tutorial 13: Model Stacking</h2></center>\n",
    "<br>\n",
    "\n",
    "In this notebook we build on the analysis of Tutorial 12 to consider model stacking. \n",
    "\n",
    "<a class=\"buttom\" href=\">#Data:-California-Housing\">Data: California Housing</a> <br>\n",
    "<a class=\"buttom\" href=\"#Previous-results\">Previous results</a> <br>\n",
    "<a class=\"buttom\" href=\"#Adding-new-models\">Generalised additive modelling</a> <br>\n",
    "<a class=\"buttom\" href=\"#Model-stacking\">Model stacking</a> <br>\n",
    "<a class=\"buttom\" href=\"#Model-evaluation\">Model Evaluation</a> <br>\n",
    "\n",
    "This notebook relies on the following libraries and settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Packages\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Methods\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.linear_model import LassoCV, RidgeCV\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "from sklearn.ensemble import BaggingRegressor, RandomForestRegressor, GradientBoostingRegressor\n",
    "from statlearning import generalised_additive_regression\n",
    "from statlearning import rmse_jack, r2_jack, mae"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Data processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>MedInc</th>\n",
       "      <th>HouseAge</th>\n",
       "      <th>AveRooms</th>\n",
       "      <th>AveBedrms</th>\n",
       "      <th>Population</th>\n",
       "      <th>AveOccup</th>\n",
       "      <th>Latitude</th>\n",
       "      <th>Longitude</th>\n",
       "      <th>MedianHouseValue</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>8.3252</td>\n",
       "      <td>41.0</td>\n",
       "      <td>6.984127</td>\n",
       "      <td>1.023810</td>\n",
       "      <td>322.0</td>\n",
       "      <td>2.555556</td>\n",
       "      <td>37.88</td>\n",
       "      <td>-122.23</td>\n",
       "      <td>4.526</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>8.3014</td>\n",
       "      <td>21.0</td>\n",
       "      <td>6.238137</td>\n",
       "      <td>0.971880</td>\n",
       "      <td>2401.0</td>\n",
       "      <td>2.109842</td>\n",
       "      <td>37.86</td>\n",
       "      <td>-122.22</td>\n",
       "      <td>3.585</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>7.2574</td>\n",
       "      <td>52.0</td>\n",
       "      <td>8.288136</td>\n",
       "      <td>1.073446</td>\n",
       "      <td>496.0</td>\n",
       "      <td>2.802260</td>\n",
       "      <td>37.85</td>\n",
       "      <td>-122.24</td>\n",
       "      <td>3.521</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>5.6431</td>\n",
       "      <td>52.0</td>\n",
       "      <td>5.817352</td>\n",
       "      <td>1.073059</td>\n",
       "      <td>558.0</td>\n",
       "      <td>2.547945</td>\n",
       "      <td>37.85</td>\n",
       "      <td>-122.25</td>\n",
       "      <td>3.413</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3.8462</td>\n",
       "      <td>52.0</td>\n",
       "      <td>6.281853</td>\n",
       "      <td>1.081081</td>\n",
       "      <td>565.0</td>\n",
       "      <td>2.181467</td>\n",
       "      <td>37.85</td>\n",
       "      <td>-122.25</td>\n",
       "      <td>3.422</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   MedInc  HouseAge  AveRooms  AveBedrms  Population  AveOccup  Latitude  \\\n",
       "0  8.3252      41.0  6.984127   1.023810       322.0  2.555556     37.88   \n",
       "1  8.3014      21.0  6.238137   0.971880      2401.0  2.109842     37.86   \n",
       "2  7.2574      52.0  8.288136   1.073446       496.0  2.802260     37.85   \n",
       "3  5.6431      52.0  5.817352   1.073059       558.0  2.547945     37.85   \n",
       "4  3.8462      52.0  6.281853   1.081081       565.0  2.181467     37.85   \n",
       "\n",
       "   Longitude  MedianHouseValue  \n",
       "0    -122.23             4.526  \n",
       "1    -122.22             3.585  \n",
       "2    -122.24             3.521  \n",
       "3    -122.25             3.413  \n",
       "4    -122.25             3.422  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv('cal_housing.csv')\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data=data[data['AveRooms']<data['AveRooms'].quantile(.99)]\n",
    "data=data[data['Population']<data['Population'].quantile(.99)]\n",
    "data=data[data['AveOccup']<data['AveOccup'].quantile(.99)]\n",
    "\n",
    "response = data.columns[-1]\n",
    "predictors= list(data.columns[:-1])\n",
    "    \n",
    "train = data.sample(frac=0.2, random_state=1)\n",
    "test = data[data.index.isin(train.index)==False].copy()\n",
    "\n",
    "y_train = np.log(train[response])\n",
    "y_test = np.log(test[response])\n",
    "X_train = train[predictors]\n",
    "X_test = test[predictors]\n",
    "\n",
    "X_train_tree = X_train.copy()\n",
    "X_test_tree = X_test.copy()\n",
    "\n",
    "mu=X_train.mean()\n",
    "sigma=X_train.std()\n",
    "\n",
    "X_train=(X_train-mu)/sigma\n",
    "X_test=(X_test-mu)/sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Previous results\n",
    "\n",
    "This selection replicates the results of the previous tutorial, where we found that gradient boosting was the most accurate method for the test data among a range of candidates including a linear regression, a generalised additive model based on regression splines, and a random forest. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# OLS\n",
    "ols = LinearRegression()\n",
    "ols.fit(X_train, y_train)\n",
    "\n",
    "# Lasso\n",
    "lasso = LassoCV(cv=5)\n",
    "lasso.fit(X_train, y_train)\n",
    "\n",
    "# Ridge\n",
    "alphas = np.exp(np.linspace(-10,20,500)) \n",
    "ridge = RidgeCV(alphas=alphas, cv=5)\n",
    "ridge.fit(X_train, y_train)\n",
    "\n",
    "# Decision Tree\n",
    "tree = DecisionTreeRegressor(max_depth=9, min_samples_leaf=10)\n",
    "tree.fit(X_train_tree, y_train)\n",
    "\n",
    "# Bagged trees\n",
    "bag = BaggingRegressor(n_estimators=1000, random_state=1)\n",
    "bag.fit(X_train_tree, y_train)\n",
    "\n",
    "# Random forest\n",
    "rf = RandomForestRegressor(n_estimators=5000, max_features = 3, min_samples_leaf= 1)\n",
    "rf.fit(X_train_tree, y_train)\n",
    "\n",
    "# Gradient boosting\n",
    "gb = GradientBoostingRegressor(learning_rate= 0.05, max_depth = 5, n_estimators= 750)\n",
    "gb.fit(X_train_tree, y_train)\n",
    "\n",
    "# Regression splines \n",
    "nonlinear=['MedInc', 'HouseAge', 'AveOccup', 'AveRooms', 'Population', 'AveBedrms', 'Latitude', 'Longitude']\n",
    "dfs=[8,8,5,3,2,6,18,18]\n",
    "spline_dfs=dict(zip(nonlinear,dfs))\n",
    "\n",
    "gam = generalised_additive_regression()\n",
    "gam.fit(X_train, y_train, spline_dfs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Test RMSE</th>\n",
       "      <th>SE</th>\n",
       "      <th>Test R2</th>\n",
       "      <th>SE</th>\n",
       "      <th>Test MAE</th>\n",
       "      <th>SE</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>OLS</th>\n",
       "      <td>0.324</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.676</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.244</td>\n",
       "      <td>0.002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ridge</th>\n",
       "      <td>0.324</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.676</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.244</td>\n",
       "      <td>0.002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso</th>\n",
       "      <td>0.324</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.676</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.244</td>\n",
       "      <td>0.002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GAM</th>\n",
       "      <td>0.267</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.779</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.196</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decision tree</th>\n",
       "      <td>0.313</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.697</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.225</td>\n",
       "      <td>0.002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bagged trees</th>\n",
       "      <td>0.255</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.799</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.181</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Random forest</th>\n",
       "      <td>0.249</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.808</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.179</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gradient boosting</th>\n",
       "      <td>0.232</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.833</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.164</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   Test RMSE     SE  Test R2     SE  Test MAE     SE\n",
       "OLS                    0.324  0.003    0.676  0.005     0.244  0.002\n",
       "Ridge                  0.324  0.003    0.676  0.005     0.244  0.002\n",
       "Lasso                  0.324  0.003    0.676  0.005     0.244  0.002\n",
       "GAM                    0.267  0.003    0.779  0.004     0.196  0.001\n",
       "Decision tree          0.313  0.003    0.697  0.006     0.225  0.002\n",
       "Bagged trees           0.255  0.003    0.799  0.004     0.181  0.001\n",
       "Random forest          0.249  0.003    0.808  0.004     0.179  0.001\n",
       "Gradient boosting      0.232  0.003    0.833  0.004     0.164  0.001"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "columns=['Test RMSE', 'SE', 'Test R2', 'SE', 'Test MAE', 'SE']\n",
    "rows=['OLS', 'Ridge', 'Lasso', 'GAM', 'Decision tree', 'Bagged trees', 'Random forest', 'Gradient boosting']\n",
    "results=pd.DataFrame(0.0, columns=columns, index=rows) \n",
    "\n",
    "methods=[ols, ridge, lasso, gam, tree, bag, rf, gb]\n",
    "\n",
    "for i, method in enumerate(methods):\n",
    "    if i < 4:\n",
    "        y_pred=method.predict(X_test)\n",
    "    else:\n",
    "        y_pred=method.predict(X_test_tree)\n",
    "    \n",
    "    results.iloc[i,0], results.iloc[i,1] = rmse_jack(y_test, y_pred)\n",
    "    results.iloc[i,2], results.iloc[i,3] = r2_jack(y_test, y_pred)\n",
    "    results.iloc[i,4], results.iloc[i,5] = mae(y_test, y_pred)\n",
    "\n",
    "results.round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Generalised additive modelling \n",
    "\n",
    "In order to obtain accuracy gains from model stacking, we diversify our range of methods to consider a more tailored generalised additive model (GAM). Our GAM will combine individual regression splines for most with a local regression for the spatial effect.\n",
    "\n",
    "The model has the form: \n",
    "\n",
    "\\begin{equation}\n",
    "Y=\\beta_0+\\left(\\sum_{j=1}^{p-2}\\,f_j(X_j)\\right)+g\\,(\\text{latitude},\\text{longitude})+\\varepsilon.\n",
    "\\end{equation}\n",
    "\n",
    "We will estimate $f_j(\\cdot)$ with natural splines and $g(\\cdot,)$ with a kernel weighted average (local constant model). The flexibility of a nonparametric regression is well suited for estimating the spatial effects.    \n",
    "\n",
    "To illustrate the use of local regression as building block, the next cell estimates a local linear regression with three predictors:  median income, latitude, and longitude. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 41min 8s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "from statsmodels.nonparametric.kernel_regression import KernelReg\n",
    "kernel = KernelReg(y_train, X_train.iloc[:,[0,6,7]], var_type=3*'c', reg_type='ll')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While a local regression can be very flexible, it is has the disadvantages of being (i) subject to the curse of dimensionality, which we mitigate by considering only a few predictors (ii) very computationally costly. \n",
    "\n",
    "We estimate the generalised additive regression by backfitting (see Chapter 9 of ESL), which we explicly demonstrate below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 2h 48min 50s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# Degrees of freedom for the regression splines\n",
    "nonlinear=['MedInc', 'HouseAge', 'AveOccup', 'AveRooms', 'Population', 'AveBedrms']\n",
    "dfs=[8,8,5,3,2,6]\n",
    "spline_dfs=dict(zip(nonlinear,dfs))\n",
    "\n",
    "# Backfitting algorithm\n",
    "resid1 = np.copy(y_train)  # initialising the residual \n",
    "for i in range(10):\n",
    "    gam2 = generalised_additive_regression()\n",
    "    gam2.fit(X_train.iloc[:,:-2], resid1, spline_dfs)\n",
    "    resid2 = y_train - gam2.predict(X_train.iloc[:,:-2])\n",
    "    \n",
    "    gam1 = KernelReg(resid2, X_train.iloc[:,[6,7]], var_type=2*'c', reg_type='lc')\n",
    "    resid1 = y_train-gam1.fit(X_train.iloc[:,[6,7]])[0]\n",
    "    \n",
    "    \n",
    "gam2 = generalised_additive_regression()\n",
    "gam2.fit(X_train.iloc[:,:-2], resid1, spline_dfs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Model Stacking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The <TT>statlearning</TT> module provides functionality for model stacking based on linear and local regressions as meta models. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statlearning import stack_design_matrix, linear_stack, local_stack"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To fit a model stack, we should ideally use cross validation predictions as features, or a validation set. It would be extremely time consuming to conduct cross validation for the GAM, so that we use fitted values for illustrative purposes only.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 39.9 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "X_train_stack= stack_design_matrix([gb], [X_train_tree], y_train, cv=10)\n",
    "X2 = gam1.fit(X_train.iloc[:,[6,7]])[0] + gam2.predict(X_train)\n",
    "X_train_stack = np.hstack((X_train_stack, X2.reshape((-1,1))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We fit the model stack as follows. The <TT>normalise=True</TT> option leads to a weighted average. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.028,  0.972])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack1 = linear_stack()\n",
    "stack1.fit(X_train_stack, y_train, normalise=True)\n",
    "stack1.beta.round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stack places almost all the weight on the GAM, since we \"unfairly\" used fitted values rather than cross validation predictions for it. We therefore change the stack to be a simple average for to compute the model evaluation results below.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "stack1.beta=np.array([0.5,0.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For illustrative purposes, we can estimate a local regression as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 7min 4s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "stack2 = local_stack()\n",
    "stack2.fit(X_train_stack, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Model evaluation\n",
    "\n",
    "The results show that our generalised additive model (labelled Splines + Local) outperforms boosting, and that the average of the two models outpeforms both of them by a meaningful margin. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Test RMSE</th>\n",
       "      <th>SE</th>\n",
       "      <th>Test R2</th>\n",
       "      <th>SE</th>\n",
       "      <th>Test MAE</th>\n",
       "      <th>SE</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>OLS</th>\n",
       "      <td>0.314</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.696</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.235</td>\n",
       "      <td>0.002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Regression Splines</th>\n",
       "      <td>0.267</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.780</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.195</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Local Regression</th>\n",
       "      <td>0.276</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.764</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.201</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decision tree</th>\n",
       "      <td>0.313</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.697</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.225</td>\n",
       "      <td>0.002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bagged trees</th>\n",
       "      <td>0.255</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.799</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.181</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Random forest</th>\n",
       "      <td>0.249</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.808</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.179</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gradient boosting</th>\n",
       "      <td>0.232</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.833</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.163</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Splines + Local</th>\n",
       "      <td>0.227</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.841</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.156</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear Stack</th>\n",
       "      <td>0.214</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.858</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.148</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Local stack</th>\n",
       "      <td>0.228</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.839</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.155</td>\n",
       "      <td>0.001</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                    Test RMSE     SE  Test R2     SE  Test MAE     SE\n",
       "OLS                     0.314  0.003    0.696  0.005     0.235  0.002\n",
       "Regression Splines      0.267  0.003    0.780  0.004     0.195  0.001\n",
       "Local Regression        0.276  0.003    0.764  0.005     0.201  0.001\n",
       "Decision tree           0.313  0.003    0.697  0.006     0.225  0.002\n",
       "Bagged trees            0.255  0.003    0.799  0.004     0.181  0.001\n",
       "Random forest           0.249  0.003    0.808  0.004     0.179  0.001\n",
       "Gradient boosting       0.232  0.003    0.833  0.004     0.163  0.001\n",
       "Splines + Local         0.227  0.003    0.841  0.004     0.156  0.001\n",
       "Linear Stack            0.214  0.003    0.858  0.004     0.148  0.001\n",
       "Local stack             0.228  0.003    0.839  0.004     0.155  0.001"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "columns=['Test RMSE', 'SE', 'Test R2', 'SE', 'Test MAE', 'SE']\n",
    "rows=['OLS', 'Regression Splines', 'Local Regression', 'Decision tree', 'Bagged trees', \n",
    "      'Random forest', 'Gradient boosting', 'Splines + Local', 'Linear Stack', 'Local stack']\n",
    "results=pd.DataFrame(0.0, columns=columns, index=rows) \n",
    "\n",
    "methods=[ols, gam, kernel, tree, bag, rf, gb]\n",
    "\n",
    "y_preds = np.zeros((len(y_test),len(rows)))\n",
    "for i, method in enumerate(methods):\n",
    "    if i < 2:\n",
    "        y_preds[:,i] = method.predict(X_test)\n",
    "    elif i==2:\n",
    "        y_preds[:,i] = method.fit(X_test.iloc[:,[0,6,7]])[0]\n",
    "    else:\n",
    "        y_preds[:,i] = method.predict(X_test_tree)\n",
    "        \n",
    "y_preds[:,7] = gam1.fit(X_test.iloc[:,[6,7]])[0] + gam2.predict(X_test)\n",
    "\n",
    "# This block truncates implausible predictions\n",
    "y_preds[y_preds>1.609]=1.609\n",
    "y_preds[y_preds<-2]=-2\n",
    "\n",
    "# In case of numerical instability: \n",
    "y_preds[np.isnan(y_preds[:,7]), 7]=y_preds[np.isnan(y_preds[:,7]),6]\n",
    "\n",
    "X_test_stack = y_preds[:,[6,7]]\n",
    "y_preds[:,-2] = stack1.predict(X_test_stack)\n",
    "y_preds[:,-1] = stack2.predict(X_test_stack)\n",
    "\n",
    "for i in range(len(rows)):\n",
    "    results.iloc[i,0], results.iloc[i,1] = rmse_jack(y_test, y_preds[:,i])\n",
    "    results.iloc[i,2], results.iloc[i,3] = r2_jack(y_test, y_preds[:,i])\n",
    "    results.iloc[i,4], results.iloc[i,5] = mae(y_test, y_preds[:,i])\n",
    "\n",
    "results.round(3)"
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}