1survey.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook 1: Survey data manipulation\n",
    "Replication code for:\n",
    "- Figure S1\n",
    "- Figure S3\n",
    "- Figure S6\n",
    "- Figure S7\n",
    "- Figure S14\n",
    "- Figure S15\n",
    "- Table S7\n",
    "- Table S8\n",
    "- Table S9\n",
    "- Table S10\n",
    "- Table S11\n",
    "- Table S16"
   ]
  },
  {
   "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 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",
    "from scipy.stats import percentileofscore\n",
    "import matplotlib.ticker as mtick\n",
    "from multiprocessing import Pool\n",
    "from sklearn.model_selection import train_test_split, cross_validate, KFold\n",
    "from sklearn.linear_model import LinearRegression, Lasso, LassoCV, RidgeCV, ElasticNet, Ridge\n",
    "from sklearn.linear_model import lars_path\n",
    "from sklearn.metrics import r2_score\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.linear_model.base import _rescale_data\n",
    "from sklearn.model_selection import GridSearchCV, cross_val_score, cross_val_predict\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "from IPython.display import display, clear_output\n",
    "\n",
    "from helpers import *\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure S7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pmt_feature_selection(df, outname):\n",
    "    # Normalize continuous variables\n",
    "    continuous_columns = [col for col in components if col[:4] == 'cont']\n",
    "    for col in continuous_columns:\n",
    "        df[col] = (df[col] - min(df[col]))/(max(df[col]) - min(df[col]))\n",
    "\n",
    "    # Convert categorical variables to dummies\n",
    "    categorical_columns = [col for col in components if col[:3] == 'cat']\n",
    "    df = pd.get_dummies(df, columns=categorical_columns, drop_first=True, prefix_sep='__')\n",
    "\n",
    "    # Train/test split\n",
    "    train, test = train_test_split(df, test_size=.25, random_state=6, shuffle=True)\n",
    "\n",
    "    xcols = train.drop(['consumption', 'weight'], axis=1).columns\n",
    "    x = train.drop(['consumption', 'weight'], axis=1).values\n",
    "    y = train['consumption'].values\n",
    "    weights = train['weight'].values\n",
    "    x_test = test.drop(['consumption', 'weight'], axis=1).values\n",
    "    y_test = test['consumption'].values\n",
    "    weights_test = test['weight'].values\n",
    "    x = pd.DataFrame(x, columns = xcols)\n",
    "    \n",
    "    # Define models\n",
    "    fit_params={'sample_weight':weights}\n",
    "    rf = GridSearchCV(RandomForestRegressor(random_state=9), \n",
    "                       param_grid={'max_depth':range(1, 8, 2), 'n_estimators':[50, 100]}, \n",
    "                       cv=3, n_jobs=5)\n",
    "    lr = LinearRegression()\n",
    "    \n",
    "    # Stepwise forward selection of predictors\n",
    "    for (model, name) in [(lr, 'LR'), (lr, 'RF')]:\n",
    "        used_features = []\n",
    "        unused_features = components\n",
    "        train_scores, test_scores = [], []\n",
    "        for i in range(30):\n",
    "            potential_model_test_scores, potential_model_train_scores = [], []\n",
    "            for feature in unused_features:\n",
    "                clear_output(wait=True)\n",
    "                print(('%i ' + feature) % i)\n",
    "                x_selected = x[inmultiple(used_features, x.columns) + inmultiple([feature], x.columns)]\n",
    "                results = cv(model, x_selected, y, weights)\n",
    "                potential_model_train_scores.append(results[0])\n",
    "                potential_model_test_scores.append(results[1])\n",
    "            best_idx = np.argmax(potential_model_test_scores)\n",
    "            best_feature = unused_features[best_idx]\n",
    "            used_features.append(best_feature)\n",
    "            train_scores.append(potential_model_train_scores[best_idx])\n",
    "            test_scores.append(potential_model_test_scores[best_idx])\n",
    "            print(train_scores[-1], test_scores[-1])\n",
    "            unused_features = list(set(unused_features) - set([best_feature]))\n",
    "        pd.DataFrame([used_features, train_scores, test_scores]).T\\\n",
    "            .to_csv('outputs/survey/' + outname + '_' + name + '.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Full PMT\n",
    "survey = pd.read_csv('data/survey2018.csv')\n",
    "components = [col for col in survey.columns if col[:3] == 'cat' or col[:3] == 'bin' or col[:4] == 'cont']\n",
    "df = survey[components + ['consumption', 'weight']].copy()\n",
    "pmt_feature_selection(df, 'pmt_feature_selection')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rural-only PMT\n",
    "survey = pd.read_csv('data/survey2018.csv')\n",
    "survey = survey[survey['milieu'] == 'rural']\n",
    "components = [col for col in survey.columns if col[:3] == 'cat' or col[:3] == 'bin' or col[:4] == 'cont']\n",
    "df = survey[components + ['consumption', 'weight']].copy()\n",
    "pmt_feature_selection(df, 'pmt_feature_selection_rural')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(style='white', font_scale=2.3)\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(20, 7))\n",
    "\n",
    "# Linear Regression\n",
    "features = pd.read_csv('outputs/survey/pmt_feature_selection_LR.csv')\n",
    "features.columns = ['feature', 'train', 'test']\n",
    "ax[0].plot(range(len(features['train'])), features['train'], color='mediumseagreen', label='Train')\n",
    "ax[0].plot(range(len(features['test'])), features['test'], color='indianred', label='Test')\n",
    "ax[0].scatter(range(len(features['train'])), features['train'], color='mediumseagreen')\n",
    "ax[0].scatter(range(len(features['test'])), features['test'], color='indianred')\n",
    "ax[0].axvline(11, color='grey', dashes=[2, 2])\n",
    "ax[0].set_ylabel('r2 Score')\n",
    "ax[0].set_xlabel('Number of Features')\n",
    "simpleaxis(ax[0])\n",
    "\n",
    "# Random Forest\n",
    "features = pd.read_csv('outputs/survey/pmt_feature_selection_RF.csv')\n",
    "features.columns = ['feature', 'train', 'test']\n",
    "ax[1].plot(range(len(features['train'])), features['train'], color='mediumseagreen', label='Train')\n",
    "ax[1].plot(range(len(features['test'])), features['test'], color='indianred', label='Test')\n",
    "ax[1].scatter(range(len(features['train'])), features['train'], color='mediumseagreen')\n",
    "ax[1].scatter(range(len(features['test'])), features['test'], color='indianred')\n",
    "ax[1].axvline(11, color='grey', dashes=[2, 2])\n",
    "ax[1].set_ylabel('r2 Score')\n",
    "ax[1].set_xlabel('Number of Features')\n",
    "simpleaxis(ax[1])\n",
    "\n",
    "ax[0].set_title('Ridge Regression')\n",
    "ax[1].set_title('Random Forest')\n",
    "ax[0].set_xlim(0, 30)\n",
    "ax[1].set_xlim(0, 30)\n",
    "ax[0].set_ylim(.2, .8)\n",
    "ax[1].set_ylim(.2, .9)\n",
    "ax[1].legend(loc='lower right')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table S8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "survey = pd.read_csv('data/survey2018.csv')\n",
    "components = [col for col in survey.columns if col[:3] == 'cat' or col[:3] == 'bin' or col[:4] == 'cont']\n",
    "df = survey[components + ['consumption', 'weight']].copy()\n",
    "\n",
    "# Normalize continuous variables\n",
    "continuous_columns = [col for col in components if col[:4] == 'cont']\n",
    "for col in continuous_columns:\n",
    "    df[col] = (df[col] - min(df[col]))/(max(df[col]) - min(df[col]))\n",
    "\n",
    "# Convert categorical variables to dummies\n",
    "categorical_columns = [col for col in components if col[:3] == 'cat']\n",
    "df = pd.get_dummies(df, columns=categorical_columns, drop_first=True, prefix_sep='__')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get selected variables, process x and y\n",
    "number_vars = 12\n",
    "selected_vars = pd.read_csv('outputs/survey/pmt_feature_selection_LR.csv')\n",
    "selected_vars.columns = ['feature', 'train', 'test']\n",
    "selected_vars = selected_vars['feature'][:number_vars]\n",
    "x = df[inmultiple(selected_vars, df.columns)]\n",
    "y = np.log(df['consumption'])\n",
    "\n",
    "# Define model\n",
    "model = RidgeCV()\n",
    "\n",
    "# Get r2 over cross validation, produce OOS predictions over cross validation\n",
    "r2 = np.mean(cross_val_score(model, x, y, cv=KFold(n_splits=10, shuffle=True, random_state=12), \n",
    "                             fit_params={'sample_weight':df['weight']}))\n",
    "print('R2 score: %.2f' % r2)\n",
    "df['pmt'] = cross_val_predict(model, x, y, cv=KFold(n_splits=10, shuffle=True, random_state=12), \n",
    "                                     fit_params={'sample_weight':df['weight']})\n",
    "\n",
    "# Fit model on all data and report coefficients\n",
    "model.fit(x, y, sample_weight=df['weight'])\n",
    "coefficients = pd.DataFrame([list(x.columns), model.coef_]).T\n",
    "coefficients.columns = ['Feature', 'Coefficient']\n",
    "coefficients = coefficients.sort_values('Coefficient', ascending=False)\n",
    "coefficients['Coefficient'] = coefficients['Coefficient'].apply(lambda x: '%.2f' % x)\n",
    "coefficients.to_csv('outputs/survey/pmt_coefficients.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table S9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "survey = pd.read_csv('data/survey2018.csv')\n",
    "survey = survey[survey['milieu'] == 'rural']\n",
    "components = [col for col in survey.columns if col[:3] == 'cat' or col[:3] == 'bin' or col[:4] == 'cont']\n",
    "df = survey[components + ['consumption', 'weight']].copy()\n",
    "\n",
    "# Normalize continuous variables\n",
    "continuous_columns = [col for col in components if col[:4] == 'cont']\n",
    "for col in continuous_columns:\n",
    "    df[col] = (df[col] - min(df[col]))/(max(df[col]) - min(df[col]))\n",
    "\n",
    "# Convert categorical variables to dummies\n",
    "categorical_columns = [col for col in components if col[:3] == 'cat']\n",
    "df = pd.get_dummies(df, columns=categorical_columns, drop_first=True, prefix_sep='__')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get selected variables, process x and y\n",
    "number_vars = 12\n",
    "selected_vars = pd.read_csv('outputs/survey/pmt_feature_selection_rural_LR.csv')\n",
    "selected_vars.columns = ['feature', 'train', 'test']\n",
    "selected_vars = selected_vars['feature'][:number_vars]\n",
    "x = df[inmultiple(selected_vars, df.columns)]\n",
    "y = np.log(df['consumption'])\n",
    "\n",
    "# Define model\n",
    "model = RidgeCV()\n",
    "\n",
    "# Get r2 over cross validation, produce OOS predictions over cross validation\n",
    "r2 = np.mean(cross_val_score(model, x, y, cv=KFold(n_splits=10, shuffle=True, random_state=12), \n",
    "                             fit_params={'sample_weight':df['weight']}))\n",
    "print('R2 score: %.2f' % r2)\n",
    "df['pmt'] = cross_val_predict(model, x, y, cv=KFold(n_splits=10, shuffle=True, random_state=12), \n",
    "                                     fit_params={'sample_weight':df['weight']})\n",
    "\n",
    "# Fit model on all data and report coefficients\n",
    "model.fit(x, y, sample_weight=df['weight'])\n",
    "coefficients = pd.DataFrame([list(x.columns), model.coef_]).T\n",
    "coefficients.columns = ['Feature', 'Coefficient']\n",
    "coefficients = coefficients.sort_values('Coefficient', ascending=False)\n",
    "coefficients['Coefficient'] = coefficients['Coefficient'].apply(lambda x: '%.2f' % x)\n",
    "coefficients.to_csv('outputs/survey/pmt_coefficients_rural.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table S7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "survey = pd.read_csv('data/survey2018.csv')\n",
    "components = [col for col in survey.columns if col[:3] == 'bin']\n",
    "df = survey[components].copy()\n",
    "\n",
    "pca = PCA(n_components=1)\n",
    "scaler = MinMaxScaler()\n",
    "assets = scaler.fit_transform(df)\n",
    "df['assetindex'] = pca.fit_transform(assets)\n",
    "print('PCA variance explained: %.2f' % (100*pca.explained_variance_ratio_[0]) + '%')\n",
    "\n",
    "basis_vector = pd.DataFrame([df.columns, pca.components_[0]]).T\n",
    "basis_vector.columns = ['Asset', 'Magnitude']\n",
    "basis_vector = basis_vector.sort_values('Magnitude', ascending=False)\n",
    "basis_vector.to_csv('outputs/survey/asset_index_basis_vector.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table S10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "survey = pd.read_csv('data/survey2018.csv')\n",
    "means = survey.groupby('occupation_poverty', as_index=False).agg('mean')\n",
    "counts = survey.groupby('occupation_poverty', as_index=False).agg('count')\\\n",
    "    .rename({'consumption':'count'}, axis=1)\n",
    "sums = survey.groupby('occupation_poverty', as_index=False).agg('sum')\\\n",
    "    .rename({'weight':'total_weight'}, axis=1)\n",
    "occupations = means[['occupation_poverty', 'consumption']]\\\n",
    "    .merge(counts[['occupation_poverty', 'count']])\\\n",
    "    .merge(sums[['occupation_poverty', 'total_weight']])\n",
    "occupations = occupations.rename({'occupation_poverty':'occupation_category'}, axis=1)\n",
    "occupations['Proportion'] = 100*occupations['total_weight']/occupations['total_weight'].sum()\n",
    "occupations = occupations.drop('total_weight', axis=1).round(2)\n",
    "occupations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure S3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "survey = pd.read_csv('data/survey2018.csv')\n",
    "fig, ax = plt.subplots(1, figsize=(12, 7))\n",
    "sns.kdeplot(survey[survey['formal_occupation'] == 0]['consumption'], shade=True, label='Informal Occupation')\n",
    "sns.kdeplot(survey[survey['formal_occupation'] == 1]['consumption'], shade=False, color='indianred', \n",
    "           dashes=[2, 2], label='Formal Occupation')\n",
    "ax.set_xlabel('Consumption')\n",
    "ax.set_ylabel('Density')\n",
    "simpleaxis(ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table S16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "survey = pd.read_csv('data/survey2018.csv')\\\n",
    "    .rename({'prefecture':'survey_prefecture', 'canton':'survey_canton'}, axis=1)\n",
    "home_locations = pd.read_csv('data/inferred_home_locations2018.csv')\\\n",
    "    .rename({'region':'phone_region', 'prefecture':'phone_prefecture', 'canton':'phone_canton'}, axis=1)\n",
    "df = survey.merge(home_locations, on='phone_number', how='inner')\n",
    "df = pd.DataFrame(np.repeat(df.values, df['weight'], axis=0), columns=df.columns)\n",
    "\n",
    "for spatial in ['prefecture', 'canton']:\n",
    "    for combo in [('survey', 'phone'), ('survey', 'voter'), ('phone', 'voter')]:\n",
    "        percent_matching = len(df[df[combo[0] + '_' + spatial] == df[combo[1] + '_' + spatial]])/len(df)\n",
    "        print(spatial, combo[0] + '<-->' + combo[1], ('%.2f' % percent_matching) + '%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table S11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check phone matching\n",
    "survey = pd.read_csv('data/survey2018.csv')\n",
    "survey['phone'] = survey['phone_number'].apply(lambda x: 0 if pd.isnull(x) else 1)\n",
    "phone = pd.read_csv('outputs/ml/consumption/LGBM/oos_predictions.csv')\n",
    "survey['matched'] = survey['phone_number'].apply(lambda x: 1 if x in set(phone['phone_number']) else 0)\n",
    "\n",
    "# Transform variables for aggregation\n",
    "survey['rural'] = survey['milieu'].apply(lambda x: 1 if x == 'rural' else 0)\n",
    "survey['female'] = survey['gender'].apply(lambda x: 1 if x == 'F' else 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = ['consumption', 'pmt', 'formal_occupation', 'rural', 'female', 'age']\n",
    "\n",
    "overall = pd.DataFrame(survey[cols].agg('mean'))\n",
    "overall.columns = ['Overall']\n",
    "\n",
    "by_phone = survey.groupby('phone')[cols].agg('mean').T.round(2)\n",
    "by_phone.columns = ['Phone Number', 'No Phone Number']\n",
    "\n",
    "have_phone = survey[survey['phone'] == 1].copy()\n",
    "by_match = have_phone.groupby('matched')[cols].agg('mean').T.round(2)\n",
    "by_match.columns = ['Matching', 'Not Matching']\n",
    "\n",
    "table = pd.DataFrame()\n",
    "table['Overall'] = overall['Overall']\n",
    "table['Phone Number'] = by_phone['Phone Number']\n",
    "table['No Phone Number'] = by_phone['No Phone Number']\n",
    "table['Matching'] = by_match['Matching']\n",
    "table['Not Matching'] = by_match['Not Matching']\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure S1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check phone matching\n",
    "survey = pd.read_csv('data/survey2018.csv')\n",
    "df = pd.DataFrame(np.repeat(survey.values, survey['weight'], axis=0), columns=survey.columns)\n",
    "fsec = df.groupby('fsec', as_index=False).agg('count')[['fsec', 'uid']]\\\n",
    "    .rename({'uid':'count'}, axis=1)\n",
    "fsec['count'] = 100*fsec['count']/fsec['count'].sum()\n",
    "\n",
    "sns.set(font_scale=1.5, style='white')\n",
    "fig, ax = plt.subplots(1, figsize=(10, 5))\n",
    "barlist = ax.bar(fsec['fsec'], fsec['count'])\n",
    "\n",
    "colors = ['darkgreen', 'forestgreen', 'limegreen', 'lightgreen', 'darkkhaki', 'wheat', 'navajowhite', \n",
    "          'burlywood']\n",
    "for b in range(len(barlist)):\n",
    "    barlist[b].set_color(colors[b])\n",
    "\n",
    "for i in range(len(fsec)):\n",
    "    ax.annotate(fsec['count'].values[i].round(2), (fsec['fsec'].values[i]-0.5, fsec['count'].values[i] + 1))\n",
    "    \n",
    "ax.yaxis.set_major_formatter(mtick.PercentFormatter())\n",
    "ax.set_ylim(0, 20)\n",
    "ax.set_xlabel('Days')\n",
    "\n",
    "simpleaxis(ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure S6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Panel A\n",
    "survey_indiv = pd.read_csv('data/survey_indiv2018.csv')\n",
    "\n",
    "sns.set(style='white', font_scale=1.5)\n",
    "fig, ax = plt.subplots(1, figsize=(10, 5))\n",
    "\n",
    "ages = sorted(survey_indiv['age'].unique())\n",
    "\n",
    "counts = []\n",
    "for age in ages:\n",
    "    counts.append(len(survey_indiv[survey_indiv['age'] == age]))\n",
    "ax.fill_between(ages, 0, np.array(counts)/max(counts), color='wheat')\n",
    "\n",
    "for gender, color in [('M', 'blue'), ('F', 'indianred')]:\n",
    "    subset = survey_indiv[survey_indiv['gender'] == gender].copy()\n",
    "    means, stds = [], []\n",
    "    for age in ages:\n",
    "        means.append(subset[subset['age'] == age]['own_phone'].mean())\n",
    "        stds.append(subset[subset['age'] == age]['own_phone'].mean())\n",
    "    ax.scatter(ages, means, color=color)\n",
    "    for i in range(len(means)):\n",
    "        ax.axvline(ages[i], means[i] - stds[i], means[i] + stds[i], color=color)\n",
    "\n",
    "ax.set_title('Individual Phone Ownership')\n",
    "ax.set_xlabel('Age')\n",
    "ax.set_ylabel('Proportion of Individuals with 1+ Phones')\n",
    "simpleaxis(ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Panel B\n",
    "survey = pd.read_csv('data/survey2018.csv')\n",
    "\n",
    "means, stds, counts = [], [], []\n",
    "ages = sorted(survey['age'].unique())\n",
    "for age in ages:\n",
    "    means.append(survey[survey['age'] == age]['own_phone'].mean())\n",
    "    stds.append(survey[survey['age'] == age]['own_phone'].mean())\n",
    "    counts.append(len(survey[survey['age'] == age]))\n",
    "\n",
    "sns.set(style='white', font_scale=1.5)\n",
    "fig, ax = plt.subplots(1, figsize=(10, 5))\n",
    "ax.fill_between(ages, 0, np.array(counts)/max(counts), color='wheat')\n",
    "ax.scatter(ages, means, color='black')\n",
    "for i in range(len(means)):\n",
    "    ax.axvline(ages[i], means[i] - stds[i], means[i] + stds[i], color='black')\n",
    "    \n",
    "ax.set_title('Household Head Phone Ownership')\n",
    "ax.set_xlabel('Age')\n",
    "ax.set_ylabel('Proportion of HH with 1+ Phones')\n",
    "simpleaxis(ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Panel C\n",
    "prefectures = gpd.read_file('data/shapefiles/prefectures.geojson')\n",
    "survey_indiv = pd.read_csv('data/survey_indiv2018.csv')\n",
    "ownership = survey_indiv.groupby('prefecture', as_index=False).agg('mean')[['prefecture', 'own_phone']]\n",
    "ownership = prefectures.merge(ownership, on='prefecture', how='inner')\n",
    "ownership['own_phone'] = ownership['own_phone']*100\n",
    "ownership['own_phone'] = pd.cut(ownership['own_phone'], 10)\n",
    "\n",
    "sns.reset_orig()\n",
    "fig, ax = plt.subplots(1, figsize=(10, 12))\n",
    "ownership.plot(ax=ax, column='own_phone', cmap='Reds', legend=True, edgecolor='lightgrey', linewidth=1)\n",
    "ax.axis('off')\n",
    "ax.set_title('Individual Phone Ownership by Prefecture', fontsize='xx-large')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure S14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "survey = pd.read_csv('data/survey2020.csv')\n",
    "survey = survey[survey['responded'] == 1]\n",
    "\n",
    "sns.set(style='white', font_scale=1.5)\n",
    "fig, ax = plt.subplots(1, 3, figsize=(20, 5))\n",
    "\n",
    "sns.kdeplot(survey['draw_probability'], ax=ax[0], shade=True, label='', color='indianred')\n",
    "ax[0].set_xlabel('Draw Probability')\n",
    "ax[0].set_ylabel('Density')\n",
    "ax[0].set_title('Draw Probabilities')\n",
    "simpleaxis(ax[0])\n",
    "\n",
    "sns.kdeplot(survey['response_probability'], ax=ax[1], shade=True, label='', color='orange')\n",
    "ax[1].set_xlabel('Response Probability')\n",
    "ax[1].set_ylabel('Density')\n",
    "ax[1].set_title('Response Probabilities')\n",
    "simpleaxis(ax[1])\n",
    "\n",
    "sns.kdeplot(survey['weight'], ax=ax[2], shade=True, label='', color='mediumseagreen')\n",
    "ax[2].set_xlabel('Survey Weight')\n",
    "ax[2].set_ylabel('Density')\n",
    "ax[2].set_title('Survey Weights')\n",
    "simpleaxis(ax[2])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure S15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(10, 6))\n",
    "\n",
    "survey = pd.read_csv('data/survey2020.csv')\n",
    "survey['response_bin'] = pd.qcut(survey['response_probability'], 20)\n",
    "survey = survey.groupby('response_bin', as_index=False).agg('mean')\n",
    "\n",
    "ax.scatter(survey['response_probability'], survey['responded'], s=100)\n",
    "grid = [0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1]\n",
    "ax.plot(grid, grid, color='lightgrey', dashes=[2, 2])\n",
    "\n",
    "ax.set_xlabel('True Response Probability (20 Bins)')\n",
    "ax.set_ylabel('Predicted Response Probability (20 Bins)')\n",
    "simpleaxis(ax)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "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
}