3ml.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook 3: Machine Learning with Mobile Phone Data\n",
    "Replication code for:\n",
    "- Figure S10\n",
    "- Figure S11\n",
    "- Table S4\n",
    "- Table S12\n",
    "- Table S19"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "import os\n",
    "import time\n",
    "import json\n",
    "import shutil\n",
    "import dask.dataframe as ddf\n",
    "import random\n",
    "from joblib import dump, load\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import geopandas as gpd\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier\n",
    "from sklearn.metrics import r2_score, roc_auc_score, mean_squared_error\n",
    "from sklearn.linear_model import LinearRegression, LassoCV, RidgeCV, Lasso\n",
    "from sklearn.model_selection import cross_validate, KFold, GridSearchCV, cross_val_predict, cross_val_score\n",
    "from sklearn.base import BaseEstimator, TransformerMixin, clone\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.feature_selection import VarianceThreshold\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.impute import SimpleImputer\n",
    "from lightgbm import LGBMRegressor\n",
    "from scipy.stats import spearmanr\n",
    "from joblib import dump, load\n",
    "\n",
    "from helpers import *\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cross validation initialization\n",
    "inner_strat = KFold(n_splits=3, shuffle=True, random_state=15)\n",
    "outer_strat = KFold(n_splits=5, shuffle=True, random_state=17)\n",
    "\n",
    "# LGBM Model\n",
    "lgbm_grid = {'winsorizer__limits':[(0., 1.), (0.005, .995)],\n",
    "             'model__min_data_in_leaf':[10, 20, 50], \n",
    "             'model__num_leaves':[5, 10, 20],\n",
    "             'model__learning_rate':[0.05, 0.075, 0.1],\n",
    "             'model__n_estimators':[20, 50, 100]}\n",
    "lgbm_cv = Pipeline([('winsorizer', Winsorizer()), \n",
    "                     ('model', LGBMRegressor(n_jobs=-1, random_state=2))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ml(feats_fname, survey_fname, outcome, toplevel_dir, rural_restriction=False, grid=None):\n",
    "    \n",
    "    # Load features\n",
    "    feats = pd.read_csv(feats_fname)\n",
    "    print('CDR observations: %i (%i unique)' % (len(feats), len(feats['phone_number'].unique())))\n",
    "    \n",
    "    # Load labels\n",
    "    survey = pd.read_csv(survey_fname)\n",
    "    if rural_restriction:\n",
    "        survey = survey[survey['milieu'] == 'rural']\n",
    "    if outcome in ['pmt', 'consumption', 'cons']:\n",
    "        survey[outcome] = np.log(survey[outcome])\n",
    "    print('Survey observations: %i (%i unique)' % (len(survey), len(survey['phone_number'].unique())))\n",
    "    \n",
    "    # Merge features and labels\n",
    "    survey = survey.dropna(subset=[outcome])\n",
    "    print('Survey observations with outcome: %i (%i unique)' % \n",
    "          (len(survey), len(survey['phone_number'].unique())))\n",
    "    df = survey[['phone_number', outcome, 'weight']].merge(feats, on='phone_number')\n",
    "    print('Final number of matched observations: %i (%i unique)' % (len(df), len(df['phone_number'].unique())))\n",
    "    \n",
    "    # Get X and Y \n",
    "    x = df.drop(['phone_number', 'weight', outcome], axis=1)\n",
    "    y = df[outcome]   \n",
    "    \n",
    "    modelname, model, grid = 'LGBM', lgbm_cv, lgbm_grid\n",
    "    print('Running tuned ' + modelname + '....')\n",
    "    \n",
    "    # Make folder structure\n",
    "    if not os.path.isdir(toplevel_dir):\n",
    "        os.mkdir(toplevel_dir)\n",
    "    outfolder = toplevel_dir + modelname + '/'\n",
    "    if os.path.isdir(outfolder):\n",
    "        shutil.rmtree(outfolder)\n",
    "    os.mkdir(outfolder)\n",
    "\n",
    "    # Train model over all data, tuning hyperparameters with 3 fold CV\n",
    "    model = GridSearchCV(estimator=model,\n",
    "                         param_grid=grid,\n",
    "                         cv=inner_strat,\n",
    "                         verbose=0,\n",
    "                         scoring='r2',\n",
    "                         return_train_score=True,\n",
    "                         refit=True,\n",
    "                         n_jobs=-1)\n",
    "    model.fit(x, y)\n",
    "\n",
    "    # Write the model itself\n",
    "    dump(model, outfolder + 'model')\n",
    "\n",
    "    # Record tuning results\n",
    "    resultsdf = pd.DataFrame(model.cv_results_)\n",
    "    resultsdf.to_csv(outfolder + 'tuning.csv', index=False)\n",
    "    best_row = resultsdf.iloc[resultsdf['mean_test_score'].argmax()]\n",
    "    print('Best r2 --- Train: %.2f (%.2f), Test: %.2f (%.2f)' % (best_row['mean_train_score'], \n",
    "                                                                 best_row['std_train_score'],\n",
    "                                                                 best_row['mean_test_score'],\n",
    "                                                                 best_row['std_test_score']))\n",
    "    print('Best Hyperparameters: ' + json.dumps(resultsdf.iloc[resultsdf['mean_test_score'].argmax()]['params']))\n",
    "\n",
    "    # Feature importances\n",
    "    try:\n",
    "        imports = model.best_estimator_.named_steps['model'].feature_importances_\n",
    "    except:\n",
    "        imports = model.best_estimator_.named_steps['model'].coef_\n",
    "    imports = pd.DataFrame([x.columns, imports]).T\n",
    "    imports.columns = ['Feature', 'Importance']\n",
    "    imports = imports.sort_values('Importance', ascending=False)\n",
    "    imports.to_csv(outfolder + 'importances.csv')\n",
    "    \n",
    "    # Out of sample predictions\n",
    "    oos_preds = cross_val_predict(model.best_estimator_, x, y, cv=outer_strat)\n",
    "    oos_results = pd.DataFrame([list(df['phone_number']), list(y), oos_preds]).T\n",
    "    oos_results.columns = ['phone_number', 'true', 'predicted']\n",
    "    oos_results.to_csv(outfolder + 'oos_predictions.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Run machine learning model for each outcome of interest\n",
    "for outcome in ['consumption', 'pmt', 'assetindex']:\n",
    "    ml('data/features2018.csv', 'data/survey2018.csv', outcome, 'outputs/ml/' + outcome + '/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use \"old\" machine learning model to predict on \"new\" data\n",
    "old_model = load('outputs/ml/pmt2018/LGBM/model')\n",
    "new_data = pd.read_csv('data/features2018.csv')\n",
    "predictions = old_model.predict(new_data.drop('phone_number', axis=1))\n",
    "predictions = pd.DataFrame([list(new_data['phone_number']), predictions]).T\n",
    "predictions.columns = ['phone_number', 'prediction']\n",
    "predictions.to_csv('outputs/ml/consumption/temporality/old_model_new_data.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table S12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = []\n",
    "outcomes = ['consumption', 'pmt', 'assetindex']\n",
    "outcome_names = ['Consumption', 'PMT', 'Asset Index']\n",
    "\n",
    "for outcome in outcomes:\n",
    "    results = pd.read_csv('outputs/ml/' + outcome + '/LGBM/oos_predictions.csv')\n",
    "    singlefeature = pd.read_csv('data/single_feature2018.csv')\n",
    "    results = results.merge(singlefeature, on='phone_number', how='inner')\n",
    "    table.append([np.corrcoef(results['true'], results['predicted'])[0][1],\n",
    "                  np.corrcoef(results['true'], results['single_feature'])[0][1]])\n",
    "    \n",
    "table = pd.DataFrame(table).T\n",
    "table.columns = outcome_names\n",
    "table = table.round(2)\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table S4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outcome = 'consumption'\n",
    "importances = pd.read_csv('outputs/ml/' + outcome + '/LGBM/importances.csv')\\\n",
    "    [['Feature', 'Importance']]\\\n",
    "    .sort_values('Importance', ascending=False)\n",
    "importances = importances[:10]\n",
    "importances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure S10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(font_scale=2, style='white')\n",
    "outcome = 'consumption'\n",
    "poverty_line = 7.5\n",
    "display_features = ['feature0', 'feature1', 'feature2', 'feature3']\n",
    "\n",
    "results = pd.read_csv('outputs/ml/' + outcome + '/LGBM/oos_predictions.csv')\n",
    "features = pd.read_csv('data/features2018.csv')[['phone_number'] + display_features]\n",
    "results = results.merge(features, on='phone_number', how='inner')\n",
    "\n",
    "fig, ax = plt.subplots(2, 2, figsize=(20, 8))\n",
    "ax = ax.flatten()\n",
    "for f, feat in enumerate(display_features):\n",
    "    sns.kdeplot(results[results['true'] < poverty_line][feat], ax=ax[f], shade=True, label='Poor')\n",
    "    sns.kdeplot(results[results['true'] >= poverty_line][feat], ax=ax[f], shade=True, label='Nonpoor')\n",
    "    ax[f].set_title(feat.title())\n",
    "    simpleaxis(ax[f])\n",
    "ax[0].set_ylabel('Density')\n",
    "ax[2].set_ylabel('Density')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure S11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outcome = 'consumption'\n",
    "results = pd.read_csv('outputs/ml/' + outcome + '/LGBM/oos_predictions.csv')\n",
    "homes = pd.read_csv('data/inferred_home_locations2018.csv')\n",
    "results = results.merge(homes, on='phone_number', how='inner')\n",
    "\n",
    "prefectures = gpd.read_file('data/shapefiles/prefectures.geojson')\n",
    "cantons = gpd.read_file('data/shapefiles/cantons.geojson')\n",
    "\n",
    "by_prefecture = results.groupby('prefecture', as_index=False).agg('mean')[['prefecture', 'predicted']]\n",
    "counts_by_prefecture = results.groupby('prefecture', as_index=False).agg('count')[['prefecture', 'predicted']]\\\n",
    "    .rename({'predicted':'count'}, axis=1)\n",
    "by_prefecture = by_prefecture.merge(counts_by_prefecture, on='prefecture', how='inner')\n",
    "by_prefecture = prefectures.merge(by_prefecture, on='prefecture', how='inner')\n",
    "\n",
    "by_canton = results.groupby('canton', as_index=False).agg('mean')[['canton', 'predicted']]\n",
    "counts_by_canton = results.groupby('canton', as_index=False).agg('count')[['canton', 'predicted']]\\\n",
    "    .rename({'predicted':'count'}, axis=1)\n",
    "by_canton = by_canton.merge(counts_by_canton, on='canton', how='inner')\n",
    "by_canton = cantons.merge(by_canton, on='canton', how='inner')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Maps\n",
    "sns.reset_orig()\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10, 10))\n",
    "\n",
    "prefectures.plot(ax=ax[0], color='lightgrey')\n",
    "by_prefecture.plot(ax=ax[0], column='predicted', legend=True, legend_kwds={'shrink':0.5})\n",
    "\n",
    "cantons.plot(ax=ax[1], color='lightgrey')\n",
    "by_canton.plot(ax=ax[1], column='predicted', legend=True, legend_kwds={'shrink':0.5})\n",
    "\n",
    "ax[0].axis('off')\n",
    "ax[1].axis('off')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Validation against survey data\n",
    "sns.set(font_scale=1.5, style='white')\n",
    "survey = pd.read_csv('data/survey2018.csv')\n",
    "survey_by_prefecture = survey.groupby('prefecture', as_index=False).agg('mean')[['prefecture', outcome]]\n",
    "survey_by_canton = survey.groupby('canton', as_index=False).agg('mean')[['canton', outcome]]\n",
    "\n",
    "survey_by_prefecture = by_prefecture.merge(survey_by_prefecture, on='prefecture', how='inner')\n",
    "survey_by_prefecture['consumption'] = np.log(survey_by_prefecture['consumption'])\n",
    "survey_by_canton = by_canton.merge(survey_by_canton, on='canton', how='inner')\n",
    "survey_by_canton['consumption'] = np.log(survey_by_canton['consumption'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(20, 6))\n",
    "\n",
    "ax[0].scatter(survey_by_prefecture[outcome], survey_by_prefecture['predicted'], \n",
    "              s=survey_by_prefecture['count']*150, alpha=.5, c='darkorange', edgecolor='black')\n",
    "ax[1].scatter(survey_by_canton[outcome], survey_by_canton['predicted'], \n",
    "              s=survey_by_canton['count']*200, alpha=.5, c='dodgerblue', edgecolor='black')\n",
    "\n",
    "p = np.corrcoef(survey_by_prefecture[outcome], survey_by_prefecture['predicted'])[0][1]\n",
    "s = spearmanr(survey_by_prefecture[outcome], survey_by_prefecture['predicted'])[0]\n",
    "ax[0].annotate('Pearson: %.2f\\nSpearman: %.2f' % (p, s), xy=[.6, .2], xycoords='axes fraction', \n",
    "                 bbox=dict(facecolor='none', edgecolor='black', boxstyle='round,pad=1'))\n",
    "\n",
    "p = np.corrcoef(survey_by_canton[outcome], survey_by_canton['predicted'])[0][1]\n",
    "s = spearmanr(survey_by_canton[outcome], survey_by_canton['predicted'])[0]\n",
    "ax[1].annotate('Pearson: %.2f\\nSpearman: %.2f' % (p, s), xy=[.1, .8], xycoords='axes fraction', \n",
    "                 bbox=dict(facecolor='none', edgecolor='black', boxstyle='round,pad=1'))\n",
    "\n",
    "simpleaxis(ax[0])\n",
    "simpleaxis(ax[1])\n",
    "\n",
    "for a in range(len(ax)):\n",
    "    ax[a].set_title('Validation Against Survey Data')\n",
    "    ax[a].set_xlabel('Average Consumption from Survey')\n",
    "    ax[a].set_ylabel('Average Phone-Inferred Consumption')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table S19"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lgbm_cv = Pipeline([('winsorizer', Winsorizer()), \n",
    "                     ('model', LGBMRegressor(n_jobs=-1, random_state=2))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def classifier(feats_fname, survey_fname, outcome, toplevel_dir, rural_restriction=False, grid=None):\n",
    "    \n",
    "    # Load features\n",
    "    feats = pd.read_csv(feats_fname)\n",
    "    print('CDR observations: %i (%i unique)' % (len(feats), len(feats['phone_number'].unique())))\n",
    "    \n",
    "    # Load labels\n",
    "    survey = pd.read_csv(survey_fname)\n",
    "    if rural_restriction:\n",
    "        survey = survey[survey['milieu'] == 'rural']\n",
    "    if outcome in ['pmt', 'consumption', 'cons']:\n",
    "        survey[outcome] = np.log(survey[outcome])\n",
    "    print('Survey observations: %i (%i unique)' % (len(survey), len(survey['phone_number'].unique())))\n",
    "    \n",
    "    # Merge features and labels\n",
    "    survey = survey.dropna(subset=[outcome])\n",
    "    print('Survey observations with outcome: %i (%i unique)' % \n",
    "          (len(survey), len(survey['phone_number'].unique())))\n",
    "    df = survey[['phone_number', outcome, 'weight']].merge(feats, on='phone_number')\n",
    "    print('Final number of matched observations: %i (%i unique)' % (len(df), len(df['phone_number'].unique())))\n",
    "    \n",
    "    # Get X and Y \n",
    "    x = df.drop(['phone_number', 'weight', outcome], axis=1)\n",
    "    y = df[outcome]   \n",
    "    \n",
    "    modelname, model, grid = 'LGBM', lgbm_cv, lgbm_grid\n",
    "    print('Running tuned ' + modelname + '....')\n",
    "    \n",
    "    # Make folder structure\n",
    "    if not os.path.isdir(toplevel_dir):\n",
    "        os.mkdir(toplevel_dir)\n",
    "    outfolder = toplevel_dir + modelname + '/'\n",
    "    if os.path.isdir(outfolder):\n",
    "        shutil.rmtree(outfolder)\n",
    "    os.mkdir(outfolder)\n",
    "\n",
    "    # Train model over all data, tuning hyperparameters with 3 fold CV\n",
    "    model = GridSearchCV(estimator=model,\n",
    "                         param_grid=grid,\n",
    "                         cv=inner_strat,\n",
    "                         verbose=0,\n",
    "                         scoring='roc_auc',\n",
    "                         return_train_score=True,\n",
    "                         refit=True,\n",
    "                         n_jobs=-1)\n",
    "    model.fit(x, y)\n",
    "\n",
    "    # Write the model itself\n",
    "    dump(model, outfolder + 'model')\n",
    "\n",
    "    # Record tuning results\n",
    "    resultsdf = pd.DataFrame(model.cv_results_)\n",
    "    resultsdf.to_csv(outfolder + 'tuning.csv', index=False)\n",
    "    best_row = resultsdf.iloc[resultsdf['mean_test_score'].argmax()]\n",
    "    print('Best r2 --- Train: %.2f (%.2f), Test: %.2f (%.2f)' % (best_row['mean_train_score'], \n",
    "                                                                 best_row['std_train_score'],\n",
    "                                                                 best_row['mean_test_score'],\n",
    "                                                                 best_row['std_test_score']))\n",
    "    print('Best Hyperparameters: ' + json.dumps(resultsdf.iloc[resultsdf['mean_test_score'].argmax()]['params']))\n",
    "\n",
    "    # Feature importances\n",
    "    try:\n",
    "        imports = model.best_estimator_.named_steps['model'].feature_importances_\n",
    "    except:\n",
    "        imports = model.best_estimator_.named_steps['model'].coef_\n",
    "    imports = pd.DataFrame([x.columns, imports]).T\n",
    "    imports.columns = ['Feature', 'Importance']\n",
    "    imports = imports.sort_values('Importance', ascending=False)\n",
    "    imports.to_csv(outfolder + 'importances.csv')\n",
    "    \n",
    "    # Out of sample predictions\n",
    "    oos_preds = cross_val_predict(model.best_estimator_, x, y, cv=outer_strat)\n",
    "    oos_results = pd.DataFrame([list(df['phone_number']), list(y), oos_preds]).T\n",
    "    oos_results.columns = ['phone_number', 'true', 'predicted']\n",
    "    oos_results.to_csv(outfolder + 'oos_predictions.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run machine learning model for each outcome of interest\n",
    "for outcome in ['responded']:\n",
    "    classifier('data/features2020.csv', 'data/survey2020.csv', outcome, 'outputs/ml/' + outcome + '/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outcome = 'responded'\n",
    "importances = pd.read_csv('outputs/ml/' + outcome + '2020/LGBM/importances.csv')\\\n",
    "    [['Feature', 'Importance']]\\\n",
    "    .sort_values('Importance', ascending=False)\n",
    "importances = importances[:10]\n",
    "importances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}