MalawiFlood_modelling_regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from tqdm import tqdm\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import LinearRegression, Lasso, Ridge\n",
    "from sklearn.metrics import mean_squared_error as mse\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "import lightgbm as lgb\n",
    "from keras.models import Sequential, Model\n",
    "from keras.layers import Dense, Embedding, LSTM, Dropout, GRU, Conv1D, MaxPooling1D, Input, concatenate\n",
    "from keras.datasets import imdb\n",
    "from keras.callbacks import TensorBoard, EarlyStopping, ModelCheckpoint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Read in data files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "train = pd.read_pickle('features.pkl')\n",
    "predict = pd.read_csv('Train.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Set up training and testing sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_cols = ['sum','mean','elevation','poly_rainfall','poly_elevation',\n",
    "                'elev_diff','rainfall_diff','water_count','water_dist']\n",
    "features = train[feature_cols]\n",
    "target = train['target']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train-test split (80-20)\n",
    "x_train, x_val, y_train, y_val = train_test_split(features, target, test_size=0.2, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Define loss function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rmse(y_true,y_pred):\n",
    "    return np.sqrt(mse(y_true,y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4a. Prediction using Lasso regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 51/51 [00:00<00:00, 107.93it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE for lasso: 0.2183069118585445, alpha = 0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Lasso regression\n",
    "best_lasso_model = None\n",
    "best_lasso_score = float('inf')\n",
    "best_alpha = 0\n",
    "for a in tqdm(range(0,51)):\n",
    "    lasso_model = Lasso(alpha=a/10).fit(x_train, y_train)\n",
    "    val_pred_lasso = lasso_model.predict(x_val)\n",
    "    val_score_lasso = rmse(val_pred_lasso,y_val)\n",
    "    if val_score_lasso < best_lasso_score:\n",
    "        best_alpha = a\n",
    "        best_lasso_score = val_score_lasso\n",
    "        best_lasso_model = lasso_model\n",
    "print(f'RMSE for lasso: {best_lasso_score}, alpha = {best_alpha}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4b. Prediction using Ridge regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 50/50 [00:00<00:00, 175.59it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE for ridge: 0.21830698469763316, alpha = 1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Ridge regression\n",
    "best_ridge_model = None\n",
    "best_ridge_score = float('inf')\n",
    "best_alpha = 1\n",
    "for a in tqdm(range(1,51)):\n",
    "    ridge_model = Ridge(alpha=a/10).fit(x_train,y_train)\n",
    "    val_pred_ridge = ridge_model.predict(x_val)\n",
    "    val_score_ridge = rmse(val_pred_ridge,y_val)\n",
    "    if val_score_ridge < best_ridge_score:\n",
    "        best_alpha = a\n",
    "        best_ridge_score = val_score_ridge\n",
    "        best_ridge_model = ridge_model\n",
    "print(f'RMSE for ridge: {best_ridge_score}, alpha = {best_alpha}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4c. Prediction using Random Forests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE for Random Forest: 0.12262997851269453\n"
     ]
    }
   ],
   "source": [
    "# RandomForest \n",
    "rf_model = RandomForestRegressor(n_estimators = 500,\n",
    "                                min_samples_split = 2,\n",
    "                                min_samples_leaf = 1,\n",
    "                                max_samples = 0.8\n",
    "                                )\n",
    "\n",
    "# param_grid = [\n",
    "#     { 'min_samples_split':[2,4,8,16],\n",
    "#      'max_features': [2, 4, 6, 8, 9],\n",
    "#      'max_samples': [0.6,0.7,0.8],\n",
    "#      'min_impurity_decrease':[0.01,0.02,0.05]\n",
    "#     }\n",
    "# ]\n",
    "# grid_search = GridSearchCV(rf_model, param_grid, cv=5,scoring='neg_root_mean_squared_error', n_jobs=-1)\n",
    "# grid_search.fit(x_train, y_train)\n",
    "# rf_model = grid_search.best_estimator_\n",
    "rf_model.fit(x_train,y_train)\n",
    "val_pred_rf = rf_model.predict(x_val)\n",
    "val_score_rf = rmse(val_pred_rf,y_val)\n",
    "print(f'RMSE for Random Forest: {val_score_rf}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4d. Prediction using LGBM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training until validation scores don't improve for 30 rounds\n",
      "[100]\ttraining's rmse: 0.102463\tvalid_1's rmse: 0.118202\n",
      "[200]\ttraining's rmse: 0.0965204\tvalid_1's rmse: 0.116839\n",
      "Early stopping, best iteration is:\n",
      "[237]\ttraining's rmse: 0.0951186\tvalid_1's rmse: 0.116717\n"
     ]
    }
   ],
   "source": [
    "# LGBM\n",
    "lgb_params = {\n",
    "        'boosting_type': 'gbdt',         \n",
    "        'objective': 'regression',       \n",
    "        'metric': ['rmse'],             \n",
    "        'subsample': 0.5,                \n",
    "        'subsample_freq': 1,\n",
    "        'learning_rate': 0.05,           \n",
    "        'num_leaves': 2**8,            \n",
    "        'min_data_in_leaf': 2**4,      \n",
    "        'feature_fraction': 0.5,\n",
    "        'n_estimators': 5000,            \n",
    "        'early_stopping_rounds': 30,     \n",
    "        'verbose': -1,\n",
    "            } \n",
    "train_set = lgb.Dataset(x_train, y_train)\n",
    "val_set = lgb.Dataset(x_val, y_val)\n",
    "lgb_model = lgb.train(lgb_params, train_set, num_boost_round = 2000, valid_sets = [train_set, val_set], verbose_eval = 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4e. Prediction using deep neural network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [],
   "source": [
    "# training configurations\n",
    "from keras import backend as K\n",
    "\n",
    "def keras_rmse(y_true, y_pred):\n",
    "        return K.sqrt(K.mean(K.square(y_pred - y_true))) \n",
    "    \n",
    "num_epochs = 100\n",
    "batch_size = 4\n",
    "num_nodes = 256\n",
    "num_layers = 3\n",
    "dropout = 0.2\n",
    "loss_fn = keras_rmse\n",
    "optimizer = 'adagrad'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 13172 samples, validate on 3294 samples\n",
      "Epoch 1/100\n",
      "13172/13172 [==============================] - 11s 822us/step - loss: 0.1111 - accuracy: 0.8409 - val_loss: 0.1003 - val_accuracy: 0.8391\n",
      "\n",
      "Epoch 00001: val_accuracy improved from -inf to 0.83910, saving model to mlp_best_model.h5\n",
      "Epoch 2/100\n",
      "13172/13172 [==============================] - 10s 772us/step - loss: 0.1044 - accuracy: 0.8439 - val_loss: 0.0969 - val_accuracy: 0.8427\n",
      "\n",
      "Epoch 00002: val_accuracy improved from 0.83910 to 0.84274, saving model to mlp_best_model.h5\n",
      "Epoch 3/100\n",
      "13172/13172 [==============================] - 10s 775us/step - loss: 0.1047 - accuracy: 0.8440 - val_loss: 0.0959 - val_accuracy: 0.8437\n",
      "\n",
      "Epoch 00003: val_accuracy improved from 0.84274 to 0.84366, saving model to mlp_best_model.h5\n",
      "Epoch 4/100\n",
      "13172/13172 [==============================] - 10s 788us/step - loss: 0.1017 - accuracy: 0.8447 - val_loss: 0.0953 - val_accuracy: 0.8430\n",
      "\n",
      "Epoch 00004: val_accuracy did not improve from 0.84366\n",
      "Epoch 5/100\n",
      "13172/13172 [==============================] - 11s 799us/step - loss: 0.1029 - accuracy: 0.8447 - val_loss: 0.0950 - val_accuracy: 0.8437\n",
      "\n",
      "Epoch 00005: val_accuracy did not improve from 0.84366\n",
      "Epoch 6/100\n",
      "13172/13172 [==============================] - 11s 836us/step - loss: 0.1022 - accuracy: 0.8450 - val_loss: 0.0964 - val_accuracy: 0.8427\n",
      "\n",
      "Epoch 00006: val_accuracy did not improve from 0.84366\n",
      "Epoch 7/100\n",
      "13172/13172 [==============================] - 10s 787us/step - loss: 0.1018 - accuracy: 0.8451 - val_loss: 0.0964 - val_accuracy: 0.8434\n",
      "\n",
      "Epoch 00007: val_accuracy did not improve from 0.84366\n",
      "Epoch 00007: early stopping\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<keras.callbacks.callbacks.History at 0x7fc890db6ac8>"
      ]
     },
     "execution_count": 135,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# MLP\n",
    "mlp_model = Sequential()\n",
    "mlp_model.add(Dense(num_nodes, input_dim=x_train.shape[1], activation='sigmoid'))\n",
    "mlp_model.add(Dropout(dropout))\n",
    "for i in range(num_layers-1):\n",
    "    mlp_model.add(Dense(num_nodes, activation='sigmoid'))\n",
    "    mlp_model.add(Dropout(dropout))\n",
    "mlp_model.add(Dense(1, activation='sigmoid'))\n",
    "mlp_model.compile(loss=loss_fn, optimizer=optimizer, metrics=['accuracy'])\n",
    "es = EarlyStopping(monitor='val_loss', mode='min', verbose=1, patience=2)\n",
    "mc = ModelCheckpoint('mlp_best_model.h5', monitor='val_accuracy', mode='max', verbose=1, save_best_only=True)\n",
    "mlp_model.fit(x_train, y_train,\n",
    "                batch_size=batch_size,\n",
    "                epochs=num_epochs,\n",
    "                validation_data=(x_val, y_val),\n",
    "                callbacks=[es,mc],\n",
    "                verbose=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4f. Prediction using a stacking ensemble"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stacking ensemble\n",
    "\n",
    "# creating stacked dataset\n",
    "lasso_pred = best_lasso_model.predict(x_val)\n",
    "ridge_pred = best_ridge_model.predict(x_val)\n",
    "rf_pred = rf_model.predict(x_val)\n",
    "lgb_pred = lgb_model.predict(x_val)\n",
    "mlp_pred = mlp_model.predict(x_val)\n",
    "mlp_pred = mlp_pred.reshape(mlp_pred.shape[0])\n",
    "preds = [lasso_pred, ridge_pred, rf_pred, lgb_pred, mlp_pred]\n",
    "stacked_X = None\n",
    "for pred in preds:\n",
    "    if stacked_X is None:\n",
    "        stacked_X = pred\n",
    "    else:\n",
    "        stacked_X = np.dstack((stacked_X, pred))\n",
    "\n",
    "stacked_X = stacked_X.reshape(stacked_X.shape[1],stacked_X.shape[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train-test split\n",
    "cut = stacked_X.shape[0]//5\n",
    "stacked_x_val = stacked_X[:cut]\n",
    "stacked_x_train = stacked_X[cut:]\n",
    "stacked_y_val = y_val[:cut]\n",
    "stacked_y_train = y_val[cut:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 2636 samples, validate on 658 samples\n",
      "Epoch 1/100\n",
      "2636/2636 [==============================] - 2s 863us/step - loss: 0.1400 - accuracy: 0.8338 - val_loss: 0.1367 - val_accuracy: 0.8237\n",
      "Epoch 2/100\n",
      "2636/2636 [==============================] - 1s 560us/step - loss: 0.1380 - accuracy: 0.8338 - val_loss: 0.1365 - val_accuracy: 0.8237\n",
      "Epoch 3/100\n",
      "2636/2636 [==============================] - 2s 610us/step - loss: 0.1394 - accuracy: 0.8338 - val_loss: 0.1357 - val_accuracy: 0.8237\n",
      "Epoch 4/100\n",
      "2636/2636 [==============================] - 1s 546us/step - loss: 0.0893 - accuracy: 0.8467 - val_loss: 0.0676 - val_accuracy: 0.8450\n",
      "Epoch 5/100\n",
      "2636/2636 [==============================] - 2s 611us/step - loss: 0.0709 - accuracy: 0.8517 - val_loss: 0.0657 - val_accuracy: 0.8450\n",
      "Epoch 6/100\n",
      "2636/2636 [==============================] - 2s 591us/step - loss: 0.0712 - accuracy: 0.8524 - val_loss: 0.0655 - val_accuracy: 0.8450\n",
      "Epoch 7/100\n",
      "2636/2636 [==============================] - 1s 563us/step - loss: 0.0705 - accuracy: 0.8517 - val_loss: 0.0689 - val_accuracy: 0.8450\n",
      "Epoch 8/100\n",
      "2636/2636 [==============================] - 1s 562us/step - loss: 0.0722 - accuracy: 0.8513 - val_loss: 0.0650 - val_accuracy: 0.8450\n",
      "Epoch 9/100\n",
      "2636/2636 [==============================] - 2s 689us/step - loss: 0.0694 - accuracy: 0.8524 - val_loss: 0.0670 - val_accuracy: 0.8450\n",
      "Epoch 10/100\n",
      "2636/2636 [==============================] - 2s 589us/step - loss: 0.0717 - accuracy: 0.8494 - val_loss: 0.0660 - val_accuracy: 0.8450\n",
      "Epoch 11/100\n",
      "2636/2636 [==============================] - 2s 604us/step - loss: 0.0703 - accuracy: 0.8520 - val_loss: 0.0647 - val_accuracy: 0.8450\n",
      "Epoch 12/100\n",
      "2636/2636 [==============================] - 1s 569us/step - loss: 0.0681 - accuracy: 0.8520 - val_loss: 0.0659 - val_accuracy: 0.8450\n",
      "Epoch 13/100\n",
      "2636/2636 [==============================] - 2s 569us/step - loss: 0.0700 - accuracy: 0.8513 - val_loss: 0.0709 - val_accuracy: 0.8450\n",
      "Epoch 14/100\n",
      "2636/2636 [==============================] - 2s 622us/step - loss: 0.0680 - accuracy: 0.8517 - val_loss: 0.0742 - val_accuracy: 0.8450\n",
      "Epoch 00014: early stopping\n"
     ]
    }
   ],
   "source": [
    "# use MLP for stacking ensemble\n",
    "model = Sequential()\n",
    "model.add(Dense(num_nodes, input_dim=stacked_x_train.shape[1], activation='sigmoid'))\n",
    "model.add(Dropout(dropout))\n",
    "model.add(Dense(1, activation='sigmoid'))\n",
    "model.compile(loss=keras_rmse, optimizer='adam', metrics=['accuracy'])\n",
    "es = EarlyStopping(monitor='val_loss', mode='min', verbose=1, patience=3)\n",
    "mc = ModelCheckpoint('ensemble_best_model.h5', monitor='val_accuracy', mode='max', verbose=1, save_best_only=True)\n",
    "history=model.fit(stacked_x_train, stacked_y_train,\n",
    "                batch_size=batch_size,\n",
    "                epochs=100,\n",
    "                validation_data=(stacked_x_val, stacked_y_val),\n",
    "               callbacks=[es,mc],\n",
    "                 verbose=1)\n",
    "# score, accuracy = model.evaluate(stacked_x_test, stacked_y_test,batch_size=batch_size)\n",
    "# # print(f'Test score: {score}')\n",
    "# print(f'Final test accuracy: {accuracy}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare features for prediction\n",
    "predict = pd.read_csv('Train.csv')\n",
    "predict = predict.reset_index()\n",
    "predict.columns = [\"square_idx\" if x == \"index\" else x for x in predict.columns]\n",
    "cols_2019 = [col for col in predict.columns if '2019' in col] + ['square_idx']\n",
    "predict = predict[cols_2019]\n",
    "fixed = train[['id','square_idx','poly_idx','poly_elevation','elevation','elev_diff','water_count','water_dist']]\n",
    "precip = [i for i in range(1,18)]\n",
    "predict.columns = precip + ['square_idx']\n",
    "predict['sum'] = predict[precip].apply(lambda x: x.sum(),axis=1)\n",
    "predict['mean'] = predict[precip].apply(lambda x: x.mean(),axis=1)\n",
    "predict = predict.merge(fixed,on='square_idx')\n",
    "poly_rainfall = predict.groupby('poly_idx')['sum'].mean().to_frame().reset_index()\n",
    "poly_rainfall = poly_rainfall[poly_rainfall.poly_idx!=-1]\n",
    "poly_rainfall.columns = ['poly_idx','poly_rainfall']\n",
    "predict = predict.merge(poly_rainfall,on='poly_idx',how='left')\n",
    "predict.loc[predict.poly_rainfall.isna(),'poly_rainfall'] = predict['sum']\n",
    "predict['rainfall_diff'] = predict['sum'] - predict['poly_rainfall']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [],
   "source": [
    "# individual model prediction\n",
    "pred_features = predict[feature_cols]\n",
    "lasso_pred = best_lasso_model.predict(pred_features)\n",
    "ridge_pred = best_ridge_model.predict(pred_features)\n",
    "rf_pred = rf_model.predict(pred_features)\n",
    "lgb_pred = lgb_model.predict(pred_features)\n",
    "mlp_pred = mlp_model.predict(pred_features)\n",
    "mlp_pred = mlp_pred.reshape(mlp_pred.shape[0])\n",
    "preds = [lasso_pred, ridge_pred, rf_pred, lgb_pred, mlp_pred]\n",
    "stacked_X = None\n",
    "for pred in preds:\n",
    "    if stacked_X is None:\n",
    "        stacked_X = pred\n",
    "    else:\n",
    "        stacked_X = np.dstack((stacked_X, pred))\n",
    "\n",
    "stacked_X = stacked_X.reshape(stacked_X.shape[1],stacked_X.shape[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {},
   "outputs": [],
   "source": [
    "# final prediction\n",
    "predict['target_2019'] = model.predict(stacked_X)\n",
    "submission = predict[['id','target_2019']]\n",
    "submission.columns = ['Square_ID','target_2019']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 223,
   "metadata": {},
   "outputs": [],
   "source": [
    "submission.to_csv('submission.csv',index=False)"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}