2018/materials/boulder/day5-causal-inference/Day 5 - Case Study.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "pd.options.display.max_columns = 100\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.lines as mlines\n",
    "import seaborn as sb\n",
    "\n",
    "import json\n",
    "\n",
    "import networkx as nx\n",
    "\n",
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "state_codings = {1: \"Alabama\", 2: \"Alaska\", 4: \"Arizona\", 5: \"Arkansas\", 6: \"California\", \n",
    " 8: \"Colorado\", 9: \"Connecticut\", 10: \"Delaware\", 11: \"District Of Columbia\",\n",
    " 12: \"Florida\", 13: \"Georgia\", 15: \"Hawaii\", 16: \"Idaho\", 17: \"Illinois\", \n",
    " 18: \"Indiana\", 19: \"Iowa\", 20: \"Kansas\", 21: \"Kentucky\", 22: \"Louisiana\", \n",
    " 23: \"Maine\", 24: \"Maryland\", 25: \"Massachusetts\", 26: \"Michigan\", \n",
    " 27: \"Minnesota\", 28: \"Mississippi\", 29: \"Missouri\", 30: \"Montana\", \n",
    " 31: \"Nebraska\", 32: \"Nevada\", 33: \"New Hampshire\", 34: \"New Jersey\", \n",
    " 35: \"New Mexico\", 36: \"New York\", 37: \"North Carolina\", 38: \"North Dakota\", \n",
    " 39: \"Ohio\", 40: \"Oklahoma\", 41: \"Oregon\", 42: \"Pennsylvania\", \n",
    " 43: \"Puerto Rico\", 44: \"Rhode Island\", 45: \"South Carolina\", \n",
    " 46: \"South Dakota\", 47: \"Tennessee\", 48: \"Texas\", 49: \"Utah\", \n",
    " 50: \"Vermont\", 52: \"Virgin Islands\", 51: \"Virginia\", 53: \"Washington\", \n",
    " 54: \"West Virginia\", 55: \"Wisconsin\", 56: \"Wyoming\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_df = pd.read_csv('accident_counts.csv',encoding='utf8')\n",
    "counts_df['WEEKDAY'] = pd.to_datetime(counts_df[['YEAR','MONTH','DAY']]).apply(lambda x:x.weekday())\n",
    "print('{0:,} rows'.format(len(counts_df)))\n",
    "counts_df.head()\n",
    "\n",
    "population_estimates = pd.read_csv('population_estimates.csv',encoding='utf8',index_col=0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploratory data analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.factorplot(x='HOUR',y='ST_CASE',hue='WEEKDAY',data=counts_df,\n",
    "              aspect=2,order=range(24),palette='nipy_spectral',dodge=.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.factorplot(x='MONTH',y='ST_CASE',hue='WEEKDAY',data=counts_df,\n",
    "              aspect=2,order=range(1,13),palette='nipy_spectral',dodge=.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changes in fatal accidents following legalization?\n",
    "\n",
    "One interesting source of exogenous variation Colorado and Washington's legalization of cannabis in 2014. If cannabis usage increased following legalization and this translated into more impaired driving, then there should be an increase in the number of fatal auto accidents in these states after 2014. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "annual_state_counts_df = counts_df.groupby(['STATE','YEAR']).agg({'ST_CASE':np.sum,'DRUNK_DR':np.sum,'FATALS':np.sum}).reset_index()\n",
    "annual_state_counts_df = pd.merge(annual_state_counts_df,population_estimates,\n",
    "                                  left_on=['STATE','YEAR'],right_on=['State','Year']\n",
    "                                 )\n",
    "annual_state_counts_df = annual_state_counts_df[['STATE','YEAR','ST_CASE','DRUNK_DR','FATALS','Population']]\n",
    "annual_state_counts_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize accident statistics for Colorado, Washington, New Mexico (similar to Colorado), and Oregon (similar to Washington)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_cols = ['ST_CASE','DRUNK_DR','FATALS']\n",
    "annual_co_counts = annual_state_counts_df[(annual_state_counts_df['STATE'] == \"Colorado\") & (annual_state_counts_df['YEAR'] > 2010)].set_index('YEAR')[_cols]\n",
    "annual_wa_counts = annual_state_counts_df[(annual_state_counts_df['STATE'] == \"Washington\") & (annual_state_counts_df['YEAR'] > 2010)].set_index('YEAR')[_cols]\n",
    "annual_nm_counts = annual_state_counts_df[(annual_state_counts_df['STATE'] == \"New Mexico\") & (annual_state_counts_df['YEAR'] > 2010)].set_index('YEAR')[_cols]\n",
    "annual_or_counts = annual_state_counts_df[(annual_state_counts_df['STATE'] == \"Oregon\") & (annual_state_counts_df['YEAR'] > 2010)].set_index('YEAR')[_cols]\n",
    "\n",
    "# Make the figures\n",
    "f,axs = plt.subplots(3,1,figsize=(10,6),sharex=True)\n",
    "\n",
    "# Plot the cases\n",
    "annual_co_counts.plot.line(y='ST_CASE',c='blue',ax=axs[0],legend=False,lw=3)\n",
    "annual_wa_counts.plot.line(y='ST_CASE',c='green',ax=axs[0],legend=False,lw=3)\n",
    "annual_nm_counts.plot.line(y='ST_CASE',c='red',ls='--',ax=axs[0],legend=False,lw=3)\n",
    "annual_or_counts.plot.line(y='ST_CASE',c='orange',ls='--',ax=axs[0],legend=False,lw=3)\n",
    "axs[0].set_ylabel('Fatal Incidents')\n",
    "\n",
    "# Plot the drunk driving cases\n",
    "annual_co_counts.plot.line(y='DRUNK_DR',c='blue',ax=axs[1],legend=False,lw=3)\n",
    "annual_wa_counts.plot.line(y='DRUNK_DR',c='green',ax=axs[1],legend=False,lw=3)\n",
    "annual_nm_counts.plot.line(y='DRUNK_DR',c='red',ls='--',ax=axs[1],legend=False,lw=3)\n",
    "annual_or_counts.plot.line(y='DRUNK_DR',c='orange',ls='--',ax=axs[1],legend=False,lw=3)\n",
    "axs[1].set_ylabel('Drunk Driving')\n",
    "\n",
    "# Plot the fatalities\n",
    "annual_co_counts.plot.line(y='FATALS',c='blue',ax=axs[2],legend=False,lw=3)\n",
    "annual_wa_counts.plot.line(y='FATALS',c='green',ax=axs[2],legend=False,lw=3)\n",
    "annual_nm_counts.plot.line(y='FATALS',c='red',ls='--',ax=axs[2],legend=False,lw=3)\n",
    "annual_or_counts.plot.line(y='FATALS',c='orange',ls='--',ax=axs[2],legend=False,lw=3)\n",
    "axs[2].set_ylabel('Total Fatalities')\n",
    "\n",
    "# Plot 2014 legalization\n",
    "for ax in axs:\n",
    "    ax.axvline(x=2014,c='r')\n",
    "\n",
    "# Stuff for legend\n",
    "b = mlines.Line2D([],[],color='blue',label='Colorado',linewidth=3)\n",
    "g = mlines.Line2D([],[],color='green',label='Washington',linewidth=3)\n",
    "r = mlines.Line2D([],[],color='red',linestyle='--',label='New Mexico',linewidth=3)\n",
    "o = mlines.Line2D([],[],color='orange',linestyle='--',label='Oregon',linewidth=3)\n",
    "axs[2].legend(loc='lower center',ncol=4,handles=[b,g,r,o],fontsize=12,bbox_to_anchor=(.5,-.75))\n",
    "\n",
    "f.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "annual_state_counts_df['Treated'] = np.where(annual_state_counts_df['STATE'].isin(['Colorado','Washington']),1,0)\n",
    "annual_state_counts_df['Time'] = np.where(annual_state_counts_df['YEAR'] >= 2014,1,0)\n",
    "annual_state_counts_df = annual_state_counts_df[annual_state_counts_df['YEAR'] >= 2010]\n",
    "\n",
    "annual_state_counts_df.query('STATE == \"Colorado\"')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll specify a difference-in-differences design with \"Treated\" states (who legalized) and \"Time\" (when they legalized), while controlling for differences in population. The `Treated:Time` interaction is the Difference-in-Differences estimate, which is not statistically significant. This suggests legalization did not increase the risk of fatal auto accidents in these states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "m_cases = smf.ols(formula = 'ST_CASE ~ Treated*Time + Population', \n",
    "                  data = annual_state_counts_df).fit()\n",
    "print(m_cases.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_cases = smf.ols(formula = 'FATALS ~ Treated*Time + Population', \n",
    "                  data = annual_state_counts_df).fit()\n",
    "print(m_cases.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_cases = smf.ols(formula = 'DRUNK_DR ~ Treated*Time + Population', \n",
    "                  data = annual_state_counts_df).fit()\n",
    "print(m_cases.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changes in fatal accidents on 4/20?\n",
    "\n",
    "There's another exogenous source of variation in this car crash data we can exploit: the unofficial cannabis enthusiast holiday of April 20. If consumption increases on this day compared to a week before or after (April 13 and 27), does this explain changes in fatal auto accidents?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "population_estimates.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select only data after 2004 in the month of April\n",
    "april_df = counts_df.query('MONTH == 4 & YEAR > 2004').set_index(['STATE','HOUR','MONTH','DAY','YEAR'])\n",
    "\n",
    "# Re-index the data to fill in missing dates\n",
    "ix = pd.MultiIndex.from_product([state_codings.values(),range(0,24),[4],range(1,31),range(2005,2017)],\n",
    "                                names = ['STATE','HOUR','MONTH','DAY','YEAR'])\n",
    "\n",
    "april_df = april_df.reindex(ix).fillna(0)\n",
    "april_df.reset_index(inplace=True)\n",
    "\n",
    "# Add in population data\n",
    "april_df = pd.merge(april_df,population_estimates,\n",
    "                    left_on=['STATE','YEAR'],right_on=['State','Year'])\n",
    "april_df = april_df[[i for i in april_df.columns if i not in ['Year','State']]]\n",
    "\n",
    "# Inspect\n",
    "april_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate whether day is a Friday, Saturday, or Sunday\n",
    "april_df['Weekday'] = pd.to_datetime(april_df[['YEAR','MONTH','DAY']]).apply(lambda x:x.weekday())\n",
    "april_df['Weekend'] = np.where(april_df['Weekday'] >= 4,1,0)\n",
    "\n",
    "# Treated days are on April 20\n",
    "april_df['Fourtwenty'] = np.where(april_df['DAY'] == 20,1,0)\n",
    "april_df['Legal'] = np.where((april_df['STATE'].isin(['Colorado','Washington'])) & (april_df['YEAR'] >= 2014),1,0)\n",
    "\n",
    "# Examine data for a week before and after April 20\n",
    "april_df = april_df[april_df['DAY'].isin([13,20,27])]\n",
    "\n",
    "# Inspect Colorado data\n",
    "april_df.query('STATE == \"Colorado\"').sort_values(['YEAR','DAY'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect the main effect of cases on April 20 are compared to the week before and after."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.factorplot(x='DAY',y='ST_CASE',data=april_df,kind='bar',palette=['grey','green','grey'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate the models. The `Fourtwenty` and `Legal` dummy variables (and their interaction) capture whether fatal accidents increased on April 20 compared to the week beforehand and afterwards, while controlling for state legality, year, whether it's a weekend, and state population. We do not observe a statistically signidicant increase in the number of incidents, alcohol-involved incidents, and total fatalities on April 20."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_cases_420 = smf.ols(formula = 'ST_CASE ~ Fourtwenty*Legal + YEAR + Weekend + Population', \n",
    "                  data = april_df).fit()\n",
    "print(m_cases_420.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_drunk_420 = smf.ols(formula = 'DRUNK_DR ~ Fourtwenty*Legal + YEAR + Weekend + Population', \n",
    "                  data = april_df).fit()\n",
    "print(m_drunk_420.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_fatal_420 = smf.ols(formula = 'FATALS ~ Fourtwenty*Legal + YEAR + Weekend + Population', \n",
    "                  data = april_df).fit()\n",
    "print(m_fatal_420.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case study: Wikipedia pageview dynamics\n",
    "\n",
    "On November 8, 2016, California voters passed Proposition 64 legalizing recreational use of cannabis. On January 1, 2018, recreational sales began. The following two files capture the daily pageview data for the article [Cannabis in California](https://en.wikipedia.org/wiki/Cannabis_in_California) as well as the daily pageview for all the other pages it links to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ca2016 = pd.read_csv('wikipv_ca_2016.csv',encoding='utf8',parse_dates=['timestamp']).set_index('timestamp')\n",
    "ca2018 = pd.read_csv('wikipv_ca_2018.csv',encoding='utf8',parse_dates=['timestamp']).set_index('timestamp')\n",
    "ca2016.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are two plots of the pageview dynamics for the seed \"Cannabis in California\" article."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f,axs = plt.subplots(2,1,figsize=(10,5))\n",
    "ca2016['Cannabis in California'].plot(ax=axs[0],color='red',lw=3)\n",
    "ca2018['Cannabis in California'].plot(ax=axs[1],color='blue',lw=3)\n",
    "axs[0].axvline(pd.Timestamp('2016-11-08'),c='k',ls='--')\n",
    "axs[1].axvline(pd.Timestamp('2018-01-01'),c='k',ls='--')\n",
    "\n",
    "for ax in axs:\n",
    "    ax.set_ylabel('Pageviews')\n",
    "    \n",
    "f.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a visualization of the local hyperlink network around \"Cannabis in California.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = nx.read_gexf('wikilinks_cannabis_in_california.gexf')\n",
    "\n",
    "f,ax = plt.subplots(1,1,figsize=(10,10))\n",
    "\n",
    "g_pos = nx.layout.kamada_kawai_layout(g)\n",
    "\n",
    "nx.draw(G = g,\n",
    "        ax = ax,\n",
    "        pos = g_pos,\n",
    "        with_labels = True,\n",
    "        node_size = [dc*(len(g) - 1)*10 for dc in nx.degree_centrality(g).values()],\n",
    "        font_size = 7.5,\n",
    "        font_weight = 'bold',\n",
    "        node_color = 'tomato',\n",
    "        edge_color = 'grey'\n",
    "       )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What kinds of causal arguments could you make from these pageview data and the hyperlink networks?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Appendix 1: Cleaning NHTSA FARS Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"accidents.csv\" is a ~450MB file after concatenating the raw annual data from [NHTSA FARS](https://www.nhtsa.gov/research-data/fatality-analysis-reporting-system-fars)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_accident_df = pd.read_csv('accidents.csv',encoding='utf8',index_col=0)\n",
    "all_accident_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [state population estimates for 2010-2017](https://www.census.gov/data/tables/2017/demo/popest/state-total.html) from the U.S. Census Buureau."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "population_estimates = pd.read_csv('census_pop_estimates.csv')\n",
    "\n",
    "_cols = [i for i in population_estimates.columns if \"POPESTIMATE\" in i] + ['NAME']\n",
    "population_estimates_stacked = population_estimates[_cols].set_index('NAME').unstack().reset_index()\n",
    "\n",
    "population_estimates_stacked.rename(columns={'level_0':'Year','NAME':'State',0:'Population'},inplace=True)\n",
    "population_estimates_stacked['Year'] = population_estimates_stacked['Year'].str.slice(-4).astype(int)\n",
    "population_estimates_stacked = population_estimates_stacked[population_estimates_stacked['State'].isin(state_codings.values())]\n",
    "population_estimates_stacked.dropna(subset=['Population'],inplace=True)\n",
    "population_estimates_stacked.to_csv('population_estimates.csv',encoding='utf8')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Groupby-aggregate the data by state, month, day, and year counting the number of cases, alcohol-involved deaths, and total fatalities. Save the data as \"accident_counts.csv\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gb_vars = ['STATE','HOUR','DAY','MONTH','YEAR']\n",
    "agg_dict = {'ST_CASE':len,'DRUNK_DR':np.sum,'FATALS':np.sum}\n",
    "counts_df = all_accident_df.groupby(gb_vars).agg(agg_dict).reset_index()\n",
    "\n",
    "counts_df['STATE'] = counts_df['STATE'].map(state_codings)\n",
    "counts_df = counts_df.query('YEAR > 1999')\n",
    "counts_df.to_csv('accident_counts.csv',encoding='utf8',index=False)\n",
    "counts_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Appendix 2: Get page outlinks\n",
    "\n",
    "Load libraries and define two functions:\n",
    "\n",
    "* `get_page_outlinks` - Get all of the outlinks from the current version of the page.\n",
    "* `get_pageviews` - Get all of the pageviews for an article over a time range "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "import requests, json\n",
    "from bs4 import BeautifulSoup\n",
    "from urllib.parse import urlparse, quote, unquote\n",
    "import networkx as nx\n",
    "\n",
    "def get_page_outlinks(page_title,lang='en',redirects=1):\n",
    "    \"\"\"Takes a page title and returns a list of wiki-links on the page. The \n",
    "    list may contain duplicates and the position in the list is approximately \n",
    "    where the links occurred.\n",
    "    \n",
    "    page_title - a string with the title of the page on Wikipedia\n",
    "    lang - a string (typically two letter ISO 639-1 code) for the language \n",
    "        edition, defaults to \"en\"\n",
    "    redirects - 1 or 0 for whether to follow page redirects, defaults to 1\n",
    "    \n",
    "    Returns:\n",
    "    outlinks_per_lang - a dictionary keyed by language returning a dictionary \n",
    "        keyed by page title returning a list of outlinks\n",
    "    \"\"\"\n",
    "    \n",
    "    # Replace spaces with underscores\n",
    "    page_title = page_title.replace(' ','_')\n",
    "    \n",
    "    bad_titles = ['Special:','Wikipedia:','Help:','Template:','Category:','International Standard','Portal:','s:','File:','Digital object identifier','(page does not exist)']\n",
    "    \n",
    "    # Get the response from the API for a query\n",
    "    # After passing a page title, the API returns the HTML markup of the current article version within a JSON payload\n",
    "    req = requests.get('https://{2}.wikipedia.org/w/api.php?action=parse&format=json&page={0}&redirects={1}&prop=text&disableeditsection=1&disabletoc=1'.format(page_title,redirects,lang))\n",
    "    \n",
    "    # Read the response into JSON to parse and extract the HTML\n",
    "    json_string = json.loads(req.text)\n",
    "    \n",
    "    # Initialize an empty list to store the links\n",
    "    outlinks_list = [] \n",
    "    \n",
    "    if 'parse' in json_string.keys():\n",
    "        page_html = json_string['parse']['text']['*']\n",
    "\n",
    "        # Parse the HTML into Beautiful Soup\n",
    "        soup = BeautifulSoup(page_html,'lxml')\n",
    "        \n",
    "        # Remove sections at end\n",
    "        bad_sections = ['See_also','Notes','References','Bibliography','External_links']\n",
    "        sections = soup.find_all('h2')\n",
    "        for section in sections:\n",
    "            if section.span['id'] in bad_sections:\n",
    "                \n",
    "                # Clean out the divs\n",
    "                div_siblings = section.find_next_siblings('div')\n",
    "                for sibling in div_siblings:\n",
    "                    sibling.clear()\n",
    "                    \n",
    "                # Clean out the ULs\n",
    "                ul_siblings = section.find_next_siblings('ul')\n",
    "                for sibling in ul_siblings:\n",
    "                    sibling.clear()\n",
    "        \n",
    "        # Delete tags associated with templates\n",
    "        for tag in soup.find_all('tr'):\n",
    "            tag.replace_with('')\n",
    "\n",
    "        # For each paragraph tag, extract the titles within the links\n",
    "        for para in soup.find_all('p'):\n",
    "            for link in para.find_all('a'):\n",
    "                if link.has_attr('title'):\n",
    "                    title = link['title']\n",
    "                    # Ignore links that aren't interesting or are redlinks\n",
    "                    if all(bad not in title for bad in bad_titles) and 'redlink' not in link['href']:\n",
    "                        outlinks_list.append(title)\n",
    "\n",
    "        # For each unordered list, extract the titles within the child links\n",
    "        for unordered_list in soup.find_all('ul'):\n",
    "            for item in unordered_list.find_all('li'):\n",
    "                for link in item.find_all('a'):\n",
    "                    if link.has_attr('title'):\n",
    "                        title = link['title']\n",
    "                        # Ignore links that aren't interesting or are redlinks\n",
    "                        if all(bad not in title for bad in bad_titles) and 'redlink' not in link['href']:\n",
    "                            outlinks_list.append(title)\n",
    "\n",
    "    return outlinks_list\n",
    "\n",
    "def get_pageviews(page_title,lang='en',date_from='20150701',date_to=str(datetime.today().date()).replace('-','')):\n",
    "    \"\"\"Takes Wikipedia page title and returns a all the various pageview records\n",
    "    \n",
    "    page_title - a string with the title of the page on Wikipedia\n",
    "    lang - a string (typically two letter ISO 639-1 code) for the language edition,\n",
    "        defaults to \"en\"\n",
    "        datefrom - a date string in a YYYYMMDD format, defaults to 20150701\n",
    "        dateto - a date string in a YYYYMMDD format, defaults to today\n",
    "        \n",
    "    Returns:\n",
    "    revision_list - a DataFrame indexed by date and multi-columned by agent and access type\n",
    "    \"\"\"\n",
    "    quoted_page_title = quote(page_title, safe='')\n",
    "    \n",
    "    df_list = []\n",
    "    #for access in ['all-access','desktop','mobile-app','mobile-web']:\n",
    "        #for agent in ['all-agents','user','spider','bot']:\n",
    "    s = \"https://wikimedia.org/api/rest_v1/metrics/pageviews/per-article/{1}.wikipedia.org/{2}/{3}/{0}/daily/{4}/{5}\".format(quoted_page_title,lang,'all-access','user',date_from,date_to)\n",
    "    json_response = requests.get(s).json()\n",
    "    if 'items' in json_response:\n",
    "        df = pd.DataFrame(json_response['items'])\n",
    "        df_list.append(df)\n",
    "\n",
    "        concat_df = pd.concat(df_list)\n",
    "        concat_df['timestamp'] = pd.to_datetime(concat_df['timestamp'],format='%Y%m%d%H')\n",
    "        concat_df = concat_df.set_index(['timestamp','agent','access'])['views'].unstack([1,2]).sort_index(axis=1)\n",
    "        concat_df[('page','page')] = page_title\n",
    "        return concat_df\n",
    "    else:\n",
    "        print(\"Error on {0}\".format(page_title))\n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get all of the links from \"Cannabis in California\", and add the seed article itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ca_links = get_page_outlinks('Cannabis in California') + ['Cannabis in California']\n",
    "\n",
    "link_d = {}\n",
    "\n",
    "for l in list(set(ca_links)):\n",
    "    link_d[l] = get_page_outlinks(l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a network object of these hyperlinks among articles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = nx.DiGraph()\n",
    "seed_edges = [('Cannabis in California',l) for l in link_d['Cannabis in California']]\n",
    "#g.add_edges_from(seed_edges)\n",
    "for page,links in link_d.items():\n",
    "    for link in links:\n",
    "        if link in link_d['Cannabis in California']:\n",
    "            g.add_edge(page,link)\n",
    "\n",
    "print(\"There are {0:,} nodes and {1:,} edges in the network.\".format(g.number_of_nodes(),g.number_of_edges()))\n",
    "nx.write_gexf(g,'wikilinks_cannabis_in_california.gexf')            \n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the pageviews for the articles linking from \"Cannabis in California\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pvs_2016 = {}\n",
    "pvs_2018 = {}\n",
    "\n",
    "pvs_2016['Cannabis in California'] = get_pageviews('Cannabis in California',\n",
    "                                                   date_from='20160801',date_to='20170201')\n",
    "pvs_2018['Cannabis in California'] = get_pageviews('Cannabis in California',\n",
    "                                                   date_from='20171001',date_to='20180301')\n",
    "\n",
    "for page in list(set(ca_links)):\n",
    "    pvs_2016[page] = get_pageviews(page,date_from='20160801',date_to='20170201')\n",
    "    pvs_2018[page] = get_pageviews(page,date_from='20171001',date_to='20180301')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function `cleanup_pageviews` to make a rectangular DataFrame with dates as index, pages as columns, and pageviews as values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cleanup_pageviews(pv_dict):\n",
    "    _df = pd.concat(pv_dict.values())\n",
    "    _df.reset_index(inplace=True)\n",
    "    _df.columns = _df.columns.droplevel(0)\n",
    "    _df.columns = ['timestamp','pageviews','page']\n",
    "    _df = _df.set_index(['timestamp','page']).unstack('page')['pageviews']\n",
    "    return _df\n",
    "\n",
    "cleanup_pageviews(pvs_2016).to_csv('wikipv_ca_2016.csv',encoding='utf8')\n",
    "cleanup_pageviews(pvs_2018).to_csv('wikipv_ca_2018.csv',encoding='utf8')\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}