6800_Day4-5.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "6800_Day4-5",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/heckelei/PromotionskollegModule6800_2020/blob/master/6800_Day4-5.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "71QCGinKTAAf",
        "colab_type": "text"
      },
      "source": [
        "# Day 4-5: Code used during lecture and lab assignment\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MRwTIIF1s764",
        "colab_type": "text"
      },
      "source": [
        "## Instructions\n",
        "\n",
        "The notebook combines 'code used during lecture' with the corresponding lab assignment (see further down)\n",
        "Please add answers/discussion/comments to the notebook as comments or text box. Do not create another file in addition.\n",
        "When you are done with your assignment, save the notebook in drive and add your last name to the name of the file\n",
        "Upload your final notebook to https://uni-bonn.sciebo.de/s/mTpqLLBN9Wu71Ku latest September 30th, 2020. The password for access is the same as before\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UBjm-4D-Tg-L",
        "colab_type": "text"
      },
      "source": [
        "## Code used during lecture\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "f3ez9AQ5HaUN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "import scipy as sc\n",
        "from sklearn.linear_model import LinearRegression, Lasso\n",
        "from sklearn import linear_model"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DO8K6AbHaROt",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Set the numpy random seed\n",
        "np.random.seed(100)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VPoRqPZ8eM0k",
        "colab_type": "text"
      },
      "source": [
        "First generate simulated data for known DGP"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iwzJukupeFFb",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# code for creating simulated data\n",
        "N = 10000\n",
        "K = 3\n",
        "sig = 0.3\n",
        "mean = np.random.randint(0,40,K)\n",
        "# TODO repalce with automatic option to allow for different K\n",
        "cov = np.array([[30,0,0],[0,5,0],[0,0,15]])\n",
        "mean"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vWcJTDE3HkAH",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# create characteristics/explantory variables\n",
        "X = np.random.multivariate_normal(mean, cov,N)\n",
        "X.shape"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JSg6hRH6HnG1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Check that cov is correct\n",
        "np.cov(X,rowvar=False)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tsLE30gThHQ7",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Add a constant to X\n",
        "Xc = np.concatenate([np.ones([N,1]),X],axis=1)\n",
        "# Get beta values\n",
        "beta = pd.DataFrame(np.random.randint(0, 20, (Xc.shape[1],1)),index=['const']+[f\"beta{i}\" for i in range(1,Xc.shape[1])],columns=['beta_true'])\n",
        "beta\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KR3RLFaeg3mH",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# create random white noise error\n",
        "errStd = 50\n",
        "err = np.random.normal(0, errStd, (N,1))\n",
        "\n",
        "# Create dependent variable y=Xb+e\n",
        "Y = np.matmul(Xc,beta.values)+err\n",
        "Y"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1SUBJF0NgoDc",
        "colab_type": "text"
      },
      "source": [
        "Run LASSO on simulated data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eR3bOlMqgxez",
        "colab_type": "text"
      },
      "source": [
        "Compare to the true data generating process\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Y3NL86L-62l5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Estimate the model with Lasso \n",
        "# Check how alpha affects the estimation\n",
        "res = beta.transpose().copy()\n",
        "lstBetaCols = list(res.columns)\n",
        "\n",
        "lstAlpha = np.logspace(-5, 1, num = 10, base = 2)\n",
        "# Loop over a range of alpha value\n",
        "for alpha in lstAlpha:\n",
        "  # Estimate Lasso\n",
        "  modLasso = Lasso(normalize=True, fit_intercept=True, alpha=alpha)\n",
        "  modLasso.fit(X, Y[:,0])\n",
        "\n",
        "  # Add estimated coef to result dataframe\n",
        "  # Note, I rounded alpha in the variable name just to make it more readable\n",
        "  res.loc[f'beta_hat_alpha_{round(alpha,3)}','alpha'] = alpha\n",
        "  res.loc[f'beta_hat_alpha_{round(alpha,3)}','const'] = modLasso.intercept_\n",
        "  res.loc[f'beta_hat_alpha_{round(alpha,3)}',lstBetaCols[1:]] = modLasso.coef_.transpose()\n",
        "\n",
        "res"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6aGJLIeLWBJO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Plot the estimated coeficients across different levels of alpha\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "fig, ax = plt.subplots();\n",
        "for strBeta in ['beta1','beta2','beta3']:\n",
        "  ax.plot(res['alpha'],res[strBeta],label=strBeta)\n",
        "\n",
        "ax.plot([10e5,0.031250],[0,0],color='black')\n",
        "ax.set_xscale('log')\n",
        "ax.set_xlabel('alpha (log scale)')\n",
        "ax.invert_xaxis()\n",
        "ax.set_xlim(1,0.031250)\n",
        "ax.spines[\"top\"].set_visible(False)\n",
        "ax.spines[\"right\"].set_visible(False)\n",
        "ax.legend();"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yuc5bo9ha0cI",
        "colab_type": "text"
      },
      "source": [
        "Now add a additional variable that is correlated to one of the existing variables but does not have an effect on the dependent variable "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gw4K37xHLUGg",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "beta.values.squeeze()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wQYemQhShcJr",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Lasso Simulation example\n",
        "# 1) Create a variable (X4) that is linearly dependent with X3\n",
        "#    The strength of the relationship is governed by the parameter \"paramX4\"\n",
        "#    A white noise error term is added. X4 = paramX4 * X3 + err\n",
        "# 2) The strength of the relationship (paramX4) is varied in a loop \n",
        "# 3) In each step X4 is created and the Lasso is estimated. Alpha in the \n",
        "#    Lasso estimation is automatically determined using LassoCV (best alpha\n",
        "#    selected by cross-validation). \n",
        "# 4) Results are saved given: \n",
        "#     - The correlation coef between X4 and X3\n",
        "#     - The estimated coef from Lasso\n",
        "#   => This allows to see the highest coerrelation coef for which Lasso can \n",
        "#     recover the true model\n",
        "\n",
        "\n",
        "from sklearn.linear_model import LassoCV\n",
        "\n",
        "res = pd.DataFrame(columns=['corrcoef','beta1','beta2','beta3','beta4'])\n",
        "res.loc['true_model',['beta1','beta2','beta3']] = beta.iloc[1:].values.squeeze()\n",
        "for errStdX4 in [0,1,2,3,4,5,7.5,10,15,20,25,30,45,50,55,100,200]:\n",
        "  paramX4 = 2 # This parameter specifies the strength of the relationship between \n",
        "              # X3 and X4\n",
        "  X4 = paramX4*X[:,[2]]+np.random.normal(0, errStdX4, (N,1))\n",
        "\n",
        "  XPlus4 = np.concatenate([X,X4],axis=1)\n",
        "\n",
        "  res.loc[f\"errStdX4_{errStdX4}\",'corrcoef'] = np.corrcoef(XPlus4[:,2],XPlus4[:,3])[0,1]\n",
        "  # print(XPlus4[1,:])\n",
        "\n",
        "  modLassoCV = LassoCV(cv=5,normalize=True,fit_intercept=True)\n",
        "  modLassoCV.fit(XPlus4,  Y[:,0])\n",
        "  coef = modLassoCV.coef_\n",
        "  res.loc[f\"errStdX4_{errStdX4}\",'alpha'] = modLassoCV.alpha_\n",
        "  res.loc[f\"errStdX4_{errStdX4}\",['beta1','beta2','beta3','beta4']] = coef\n",
        "res"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IoA_xEY6rLjw",
        "colab_type": "text"
      },
      "source": [
        "# Lab 4a "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "faiHfDS04JDD",
        "colab_type": "text"
      },
      "source": [
        "The first part of today's lab will have you see if you can recover the true DGP using a LASSO"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lQ3JrmJXWWe1",
        "colab_type": "text"
      },
      "source": [
        "Task:\n",
        "\n",
        "- Below you will find a DGP similar to the one above\n",
        "- Try to determine the optimal alpha value in order to recover the \"best\" specification\n",
        "- Check if this \"best\" specification is close to the true DGP"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "S70WWzurdjGT",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# re-run the DGP code to get a different set of coefficients and data"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FBsQBpNDH0XV",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sklearn.datasets import make_spd_matrix"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TPNXANYLuhEI",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Call the function to generate data\n",
        "\n",
        "# Set the random seed such that everybody gets the same dataset\n",
        "np.random.seed(11)\n",
        "# Get beta values\n",
        "K = 5\n",
        "errStd = 3\n",
        "beta = pd.DataFrame(np.random.randint(-5, 5, (K+1,1)),index=['const']+[f\"beta{i}\" for i in range(1,K+1)],columns=['beta_true'])\n",
        "print(beta)\n",
        "\n",
        "# Sklearn function to generate a random symmetric, positive-definite covariance matrix\n",
        "# https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_spd_matrix.html\n",
        "Xcov = make_spd_matrix(K)\n",
        "# Generate random mean for X variables \n",
        "Xmean = np.random.randint(-20,20,K)\n",
        "\n",
        "# create characteristics/explantory variables\n",
        "X = np.random.multivariate_normal(Xmean, Xcov,N)\n",
        "\n",
        "# Add a constant to X\n",
        "Xc = np.concatenate([np.ones([N,1]),X],axis=1)\n",
        "\n",
        "# create random white noise error\n",
        "err = np.random.normal(0, errStd, (N,1))\n",
        "\n",
        "# Create dependent variable y=Xb+e\n",
        "Y = np.matmul(Xc,beta.values)+err\n",
        "\n",
        "\n",
        "# errStdXnoiseA = 20\n",
        "# errStdXnoiseB = 30\n",
        "\n",
        "errStdXnoiseA = 2\n",
        "errStdXnoiseB = 2\n",
        "            \n",
        "# XnoiseA = 2*X[:,[1]]+2*X[:,[2]]+np.random.normal(0, errStdXnoiseA, (N,1))\n",
        "# XnoiseB = 2*X[:,[2]]+1*X[:,[4]]+np.random.normal(0, errStdXnoiseB, (N,1))\n",
        "XnoiseA = X[:,[1]]*X[:,[3]]*np.random.normal(0, errStdXnoiseB, (N,1))+np.random.normal(0, 1, (N,1))\n",
        "XnoiseB = X[:,[2]]*X[:,[4]]*np.random.normal(0, errStdXnoiseB, (N,1))+np.random.normal(0, 2, (N,1))\n",
        "\n",
        "XNoise = np.concatenate([X,XnoiseA, XnoiseB],axis=1)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jBYrPexPMXFH",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sklearn.model_selection import train_test_split\n",
        "#  This randomly splits the data into 80% train and 20% test data\n",
        "X_train, X_test, Y_train, Y_test = train_test_split(XNoise, Y, test_size = 0.2)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CxjvSoJecW_F",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# ================================\n",
        "#  Task\n",
        "# ================================\n",
        "\n",
        "# Test several alpha values and check how this effects the R2 in the test set\n",
        "\n",
        "# Prepare dataframe to hold results (no need to change anything here)\n",
        "lstBeta = [f\"beta{k}\" for k in range(1,K+1)]\n",
        "res = pd.DataFrame(columns=lstBeta+['betaNoiseA','betaNoiseB'])\n",
        "res.loc['true_model',lstBeta] = beta.iloc[1:].values.squeeze()\n",
        "\n",
        "# This is the place to set alpha values that sould be tested\n",
        "lstAlpha = [...]\n",
        "\n",
        "for alpha in  lstAlpha:\n",
        "\n",
        "  # Specify the model you want to use\n",
        "  modLasso = ...\n",
        "  # Fit the model to data\n",
        "  modLasso. ...\n",
        "  \n",
        "  # Get the coeficients of the model\n",
        "  coef = \n",
        "\n",
        "  # Save results in list (no need to change anything here)\n",
        "  res.loc[f\"alpha_{alpha}\",lstBeta+['betaNoiseA','betaNoiseB']] = coef\n",
        "  res.loc[f\"alpha_{alpha}\",'alpha'] = alpha\n",
        "  res.loc[f\"alpha_{alpha}\",'R2_train'] = modLasso.score(X_train, Y_train[:,0])\n",
        "  res.loc[f\"alpha_{alpha}\",'R2_test'] = modLasso.score(X_test, Y_test[:,0])\n",
        "\n",
        "res\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HyTMoYPFd0BG",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# ================================\n",
        "#  Optional: Alternative approach using LassoCV \n",
        "# ================================\n",
        "# One alternative to setting the alphas manually is to use \n",
        "# CV this is implemented in LassoCV\n",
        "\n",
        "\n",
        "# TODO Try to implement a LassoCV estimation\n",
        "\n",
        "\n",
        "# Prepare a dataframe to hold results \n",
        "# (no need to change something here, as long as you call your model \"model\")\n",
        "print('Best Alpha: ', model.alpha_)\n",
        "lstBeta = [f\"beta{k}\" for k in range(1,K+1)]\n",
        "res = pd.DataFrame(columns=lstBeta+['betaNoiseA', 'betaNoiseB'])\n",
        "res.loc['true_model',lstBeta] = beta.iloc[1:].values.squeeze()\n",
        "res.loc['beta_hat',:] = model.coef_\n",
        "res\n",
        "# You can also check the example here:\n",
        "#  https://scikit-learn.org/stable/auto_examples/linear_model/plot_lasso_model_selection.html\n",
        "# And try to prepare similar plots"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JZm0tXJ8DLDM",
        "colab_type": "text"
      },
      "source": [
        "# Lab 4b"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "czQzN5tNDZua",
        "colab_type": "text"
      },
      "source": [
        "The objective of the second part of today's lab is to determine the effects of protected areas on forest cover. To do this, we take the following strategy:\n",
        "\n",
        "1) We estimate a model predicting forest cover where we use all of the observations that are not in a protected area. \n",
        "\n",
        "2) We use this model to predict forest cover for *all* observations. This gives as a prediction about the expected forest cover given the characteristics of the cell (such as slope, elevation, etc, just not protected status).\n",
        "\n",
        "3) We compare the predicted forest cover with the actual forest cover. The difference between the two is the estimated effect of forest protected areas. \n",
        "\n",
        "=> Most of the steps for the estimation should be familiar from previous labs."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "w9dV8MX8FRjV",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn.preprocessing import MinMaxScaler\n",
        "from sklearn.metrics import mean_squared_error, r2_score\n",
        "# import seaborn for visualization\n",
        "import seaborn as sns"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "urcpRQs9QFxO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Download data for amazon AWS store (this should only take 10s independent \n",
        "# of your internet connection)\n",
        "!wget http://www.ilr.uni-bonn.de/agpo/courses/ml/brazil_all_data_v2.gz\n",
        "\n",
        "# load data into dataframe\n",
        "df = pd.read_parquet('brazil_all_data_v2.gz')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yxNVCmreJT9M",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Define target (dependent) variable (% forest cover for 2018)\n",
        "strY = 'perc_treecover'\n",
        "\n",
        "# Define a list of features names (explantory variables)\n",
        "lstX = [\n",
        "  # 'wdpa_2017', => Exclude protected areas \n",
        "  'population_2015',\n",
        "  'chirps_2017',\n",
        "  'maize',\n",
        "  'soy',\n",
        "  'sugarcane',\n",
        "  'perm_water',\n",
        "  'travel_min',\n",
        "  'cropland',\n",
        "  'mean_elev',\n",
        "  'sd_elev',\n",
        "  'near_road',\n",
        "  # 'wdpa_2017_lag_1st_order',\n",
        "  'chirps_2017_lag_1st_order',\n",
        "  'population_2015_lag_1st_order',\n",
        "  'maize_lag_1st_order',\n",
        "  'soy_lag_1st_order',\n",
        "  'sugarcane_lag_1st_order',\n",
        "  'perc_treecover_lag_1st_order',\n",
        "  'perm_water_lag_1st_order',\n",
        "  'travel_min_lag_1st_order',\n",
        "  'cropland_lag_1st_order',\n",
        "  'mean_elev_lag_1st_order',\n",
        "  'sd_elev_lag_1st_order',\n",
        "  'near_road_lag_1st_order',\n",
        " ]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "koS5lAM3I9m7",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Split that sample in those observations that ...\n",
        "#... are protected areas\n",
        "df_PA = df.loc[~(df['wdpa_2017']==0),:]\n",
        "\n",
        "#... are NOT protected areas\n",
        "df_NoPA = df.loc[df['wdpa_2017']==0,:]\n",
        "\n",
        "print('Number of observations without protected area', df_NoPA.shape[0])\n",
        "print('Number of observations with protected area', df_PA.shape[0])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JgJR5As0FDgu",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Select the target variable \n",
        "# ... for the observations without protected areas \n",
        "Y_all_NoPA = df_NoPA[strY]\n",
        "# ... for the observations with protected areas \n",
        "Y_PA = df_PA[strY]\n",
        "\n",
        "# Get the features\n",
        "# ... for the observations without protected areas \n",
        "X_all_NoPA =  df_NoPA.loc[:,lstX]\n",
        "# ... for the observations wit protected areas \n",
        "X_PA_raw =  df_PA.loc[:,lstX]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9kH38Z7Ufk-I",
        "colab_type": "text"
      },
      "source": [
        "Choose a ML model and features to predict deforestation"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ePR5RYQ9ftsL",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Split the data into train and test using only the the observations \n",
        "# without protected areas\n",
        "X_train_raw, X_test_raw, Y_train, Y_test = train_test_split(X_all_NoPA, Y_all_NoPA, test_size = 0.2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BcujfVW9E_kb",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "# Scale data to 0-1 range using sklearn MinMaxScalar object \n",
        "# (see: https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.MinMaxScaler.html) \n",
        "scaler = MinMaxScaler()\n",
        "# Use only the train data to fit the MinMaxScalar \n",
        "scaler.fit(X_train_raw)\n",
        "\n",
        "# Apply the MinMax transformation to the train and test data \n",
        "X_train = scaler.transform(X_train_raw)\n",
        "X_test = scaler.transform(X_test_raw)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "a-IGittsEV46",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Define a function to print model stats\n",
        "def printModStats(mod,X_train, Y_train,X_test, Y_test, showPlot=True):\n",
        "  # Inspect the model performancefor those observations that are\n",
        "  # not protected areas\n",
        "  print('Score in train', mod.score(X_train, Y_train))\n",
        "  print('Score in test', mod.score(X_test, Y_test))\n",
        "\n",
        "  # Get predicted values for the test set\n",
        "  Y_test_had_Tree = mod.predict(X_test)\n",
        "\n",
        "  # Calculate MSE in test set\n",
        "  mse_ols_sklearn  = mean_squared_error(Y_test,Y_test_had_Tree)\n",
        "  print('\\nMean squared error: ',mse_ols_sklearn)\n",
        "  # The coefficient of determination: 1 is perfect prediction\n",
        "  R2_ols_sklearn = r2_score(Y_test,Y_test_had_Tree)\n",
        "  print('Coefficient of determination: ',R2_ols_sklearn)\n",
        "\n",
        "  if showPlot:\n",
        "    # plot Y vs Y-hat\n",
        "    h = sns.jointplot(Y_test_had_Tree, Y_test, kind=\"hex\")\n",
        "    h.set_axis_labels('Y predicted', 'Y true');\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8HKK5AhGf6eI",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Run XGBoost. In contrast to day 2 we now use XGBRegressor as our task is \n",
        "# a regression task and not a classification task\n",
        "import xgboost as xgb\n",
        "model_xgb = xgb.XGBRegressor(objective=\"reg:squarederror\", random_state=42)\n",
        "model_xgb.fit(X_train, Y_train)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dqIcCCsKEmX_",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "printModStats(model_xgb,X_train, Y_train,X_test, Y_test, showPlot=True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dSLkW7g224nj",
        "colab_type": "text"
      },
      "source": [
        "Now we use the trained model to make predictions for those areas without protected areas"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LW40wwuk6joa",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# First we need to scale the data for those observations without \n",
        "# protected areas\n",
        "X_PA = scaler.transform(X_PA_raw)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "u_wtagGgJwmN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Lets define a function for our conterfactual that can be reused below \n",
        "def calculateConterfactual(mod,X_PA):\n",
        "  # Then use the train model to make predictions for protected areas\n",
        "  Y_PA_had = mod.predict(X_PA)\n",
        "\n",
        "  # Now we can compare the mean forest cover in protected area that we observer\n",
        "  print(f'Predicted mean forest cover in protected Areas {np.mean(Y_PA_had):0.2f} %' )\n",
        "  print(f'Actual mean forest cover in protected Areas {np.mean(Y_PA):0.2f} %', )\n",
        "  print(f'Mean difference between predicted and observed forest cover: {np.mean(Y_PA-Y_PA_had):0.2f} pp')\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Djf7N-zj7M37",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Call function for our xgb model \n",
        "calculateConterfactual(model_xgb,X_PA)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9pxCQ8okFsQU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# ================================\n",
        "#  Task\n",
        "# ================================\n",
        "\n",
        "# Until now we have always used the default hyperparameters for XGBoost\n",
        "# Lets see if we can improve model performance by optimizing some of those.\n",
        "# Usually this is a rather complex and time consuming task with many \n",
        "# different ways to approach this. Here we only try to optimize two paramters\n",
        "# manually. Specifically we want to optimize \"max_depth\" (Typical values: 3-10) \n",
        "# and \"min_child_weight\" (Typical values: 1-5).\n",
        "\n",
        "# Hint: The caculation is quite time consuming. Coordinate in your team \n",
        "#       to try as many combination as possible.\n",
        "\n",
        "# In case you are interested you can find a complete guide for doing \n",
        "# hyperparameter optimization for XGB here: \n",
        "# https://www.analyticsvidhya.com/blog/2016/03/complete-guide-parameter-tuning-xgboost-with-codes-python/\n",
        "# This is NOT required for this task here, just to give you an idea of how\n",
        "# many paramter can be optimized (Note, this referes to an older sklearn version, \n",
        "# therefore some of the code might need to be adjusted) \n",
        "\n",
        "# Specify model with \"max_depth\" and \"min_child_weight\" set explicitly\n",
        "model_xgb_opt = ...\n",
        "\n",
        "# Fit model to train data \n",
        "model_xgb_opt....\n",
        "\n",
        "# Print stats\n",
        "printModStats(model_xgb_opt,X_train, Y_train,X_test, Y_test, showPlot=True)\n",
        "\n",
        "# Do conterfactual \n",
        "calculateConterfactual(model_xgb_opt,X_PA)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wsNLu-jzHKOk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# ================================\n",
        "#  Optional Task\n",
        "# ================================\n",
        "\n",
        "# Try alternative ML models and check how they perform compared to XGBoost. \n",
        "# If you are using sklearn you often only have to replace one line in oder to \n",
        "# change your model. LightGBM for example is an alterantive to XGBoost that \n",
        "# gained popularity.  \n",
        "# https://lightgbm.readthedocs.io/en/latest/pythonapi/lightgbm.LGBMRegressor.html#lightgbm.LGBMRegressor.  "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y0v6H7tP8Oq8",
        "colab_type": "text"
      },
      "source": [
        "## Alternative approach \n",
        "\n",
        "- Train a model with \"protected areas\" included \n",
        "- Select only the protected areas (in the test set)\n",
        "- make a prediction for those areas for forest cover\n",
        "- make a prediction for those areas for forest cover where we set the explanatory variable for \"protected areas\" equal to zero i.e. we make a prediction for forest cover under the assumption that those areas where not protected"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lyIRMR3X4xCk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Define binary variable for deforestration in 2018\n",
        "\n",
        "Y_all = df['perc_treecover']\n",
        "# Define a list of features names (explantory variables)\n",
        "lstX = [\n",
        "   'wdpa_2017', \n",
        "  'population_2015',\n",
        "  'chirps_2017',\n",
        "  'maize',\n",
        "  'soy',\n",
        "  'sugarcane',\n",
        "  'perm_water',\n",
        "  'travel_min',\n",
        "  'cropland',\n",
        "  'mean_elev',\n",
        "  'sd_elev',\n",
        "  'near_road',\n",
        "  'wdpa_2017_lag_1st_order',\n",
        "  'chirps_2017_lag_1st_order',\n",
        "  'population_2015_lag_1st_order',\n",
        "  'maize_lag_1st_order',\n",
        "  'soy_lag_1st_order',\n",
        "  'sugarcane_lag_1st_order',\n",
        "  'perc_treecover_lag_1st_order',\n",
        "  'perm_water_lag_1st_order',\n",
        "  'travel_min_lag_1st_order',\n",
        "  'cropland_lag_1st_order',\n",
        "  'mean_elev_lag_1st_order',\n",
        "  'sd_elev_lag_1st_order',\n",
        "  'near_road_lag_1st_order',\n",
        " ]\n",
        "\n",
        "# Get the explanatory Variables\n",
        "X_all =  df.loc[:,lstX]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "X1KM91dB5pwy",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "\n",
        "# Split the data into train and test using only the the observations \n",
        "# with protected areas\n",
        "\n",
        "# Use only the train data to fit the MinMaxScalar \n",
        "\n",
        "\n",
        "# Apply the MinMax transformation to the train and test data \n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "42uzaA5-5xMQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# train...\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZDyt913g8F5-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Get only those observations which are protected areas\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "B6wk-utf9AwZ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# do prediction for these as if they were not protected areas\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OAeCcDPG8iJZ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "chjqTO7h81qp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# And compare with the case they are protected\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6kMcZPeZ_s5i",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "# Now we can compare the mean forest cover in protected area that we observer\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9JNWwRkK9lRV",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# => Discuss pro and cons of the two approaches"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}