Loneliness Process Data - Wales.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import zipfile as zp\n",
    "import pandas as pd\n",
    "#from pypac import PACSession as Session #or use requests below if non-ONS\n",
    "from requests import Session\n",
    "from io import BytesIO\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating Summary Prescription Data for Wales\n",
    "\n",
    "This code has been created to work with 2016, 2017 and 2018 data, earlier data may have format inconsistencies that need to be dealt with."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set File Location of Prescribing Data\n",
    "\n",
    "These files are the \"General Practice Prescribing Data Extract\" zip files that come from: http://www.primarycareservices.wales.nhs.uk/general-practice-prescribing-data-extrac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path to zip files\n",
    "path = r\"[CHANGE THIS PATH]\\Wales\\\\\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in drug data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get drug data (NB some drugs duplicated for illnesses)\n",
    "drug_data = pd.read_csv(path + r\"..\\drug_list.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Process Prescribing Data\n",
    "\n",
    "This code iterates over the monthly prescribing data, ultimately producing an aggregate table.\n",
    "\n",
    "Note, if you want to find prescribing for non-loneliness related diseases, all you need to do is provide a different set of drug_data and edit the code_loneliness function.\n",
    "\n",
    "**Filename format:** gpdata[year][month].zip (e.g. gpdata201801.zip)\n",
    "\n",
    "**Columns:** HB, Locality, PracticeID, BNFCode, BNFName, Items, NIC, ActCost, Quantity, DDD, ADQ, Period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# names of data files\n",
    "fn_addr = \"Address\"\n",
    "fn_prescribe = \"GPData\"\n",
    "\n",
    "# column names\n",
    "col_bnfname = \"BNFName\"\n",
    "col_items = \"Items\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to find loneliness related prescribing\n",
    "def code_loneliness(x):\n",
    "    out = {}\n",
    "    # coding by illness categories\n",
    "    for illness in drug_data['illness'].unique():\n",
    "        out[illness] = x[col_bnfname].str.contains(\"|\".join(drug_data[drug_data['illness'] == illness]['medication']),\n",
    "                                                 case=False, \n",
    "                                                 regex=True).astype('int16')\n",
    "    # Make dataframe\n",
    "    out = pd.DataFrame(out)\n",
    "    # Add loneliness related disease binary - avoids double counting some drugs.\n",
    "    out['loneliness'] = x[col_bnfname].str.contains(\"|\".join(drug_data['medication'].unique()),\n",
    "                                                  case = False, \n",
    "                                                  regex = True).astype('int16')\n",
    "    # Return dataframe multiplied by counts of items.\n",
    "    return out.multiply(x[col_items], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make dictionary for aggregation\n",
    "agg_cols = {col : 'sum' for col in drug_data['illness'].unique()}\n",
    "agg_cols[col_items] = 'sum'\n",
    "agg_cols['loneliness'] = 'sum'\n",
    "for key in ['Date','HB','pcstrip','Street','Area','Posttown','Postcode']:\n",
    "    agg_cols[key] = 'first'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gpdata201801.zip\n",
      "gpdata201802.zip\n",
      "gpdata201803.zip\n",
      "gpdata201804.zip\n",
      "gpdata201805.zip\n",
      "gpdata201806.zip\n",
      "gpdata201807.zip\n",
      "gpdata201808.zip\n",
      "gpdata201809.zip\n",
      "gpdata201810.zip\n",
      "gpdata201811.zip\n",
      "gpdata201812.zip\n"
     ]
    }
   ],
   "source": [
    "monthly_data = []\n",
    "\n",
    "for file in os.listdir(path + \"Prescriptions\"):\n",
    "    with zp.ZipFile(path + \"Prescriptions\\\\\" + file) as zipf:\n",
    "        zip_names = zipf.namelist()\n",
    "        \n",
    "        # Deal with Address Files\n",
    "        addr_name = next((filename for filename in zip_names if fn_addr in filename), None)\n",
    "        # Open address file in pandas, set header.\n",
    "        addr = pd.read_csv(zipf.open(addr_name), \n",
    "                           #header=1, \n",
    "                           #names = [\"Date\", \"PracCode\", \"PracName\",\"CenterName\",\n",
    "                           #         \"Street\", \"Town\", \"Town2\", \"Postcode\"], \n",
    "                           usecols = range(8))\n",
    "        \n",
    "        # Deal with prescription info\n",
    "        prescribe_name = next((filename for filename in zip_names if fn_prescribe in filename), None)\n",
    "        # Open prescribing files in pandas.\n",
    "        prescribe = pd.read_csv(zipf.open(prescribe_name))\n",
    "        prescribe.columns = prescribe.columns.str.strip()\n",
    "        # Rename 'period' column to 'date'\n",
    "        prescribe.rename(columns = {'Period': 'Date'}, inplace = True) \n",
    "        # Get counts of prescribing dataframe for loneliness related diseases\n",
    "        loneliness_prescribing = code_loneliness(prescribe[[col_bnfname, col_items]])\n",
    "        # merge dataframes\n",
    "        prescribe = prescribe.merge(loneliness_prescribing, left_index=True, right_index=True)\n",
    "        del loneliness_prescribing\n",
    "\n",
    "        # merge in address information\n",
    "        prescribe = prescribe.merge(addr, left_on = 'PracticeID', right_on = 'PracticeId')\n",
    "        del addr\n",
    "        \n",
    "        # Create uniform postcode field\n",
    "        prescribe['pcstrip'] = prescribe['Postcode'].str.replace(\"\\s\",\"\")\n",
    "\n",
    "        # get a summary - grouping by PracCode\n",
    "        summary = prescribe.groupby('PracticeID', as_index=False).agg(agg_cols)\n",
    "        del prescribe\n",
    "\n",
    "        monthly_data.append(summary)\n",
    "        print(file)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# concatenate all the monthly data together.\n",
    "data = pd.concat(monthly_data, ignore_index = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save aggregated data\n",
    "data.to_csv(path + \"processed_data.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Add Postcode Location\n",
    "\n",
    "Postcode location is pulled in from the latest ONS NSPL (National Statistics Postcode Lookup). There is a guide to fields here: http://geoportal.statistics.gov.uk/datasets/0a404beab6f544be8fb72d0c2b12d62b\n",
    "\n",
    "NB - If using an ONS laptop `pip install pypac`, if not comment pypac import above and use requests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(path + \"processed_data.csv\", index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in postcode lookup data\n",
    "# This is the persistent link to the latest ONS NSPL\n",
    "postcode_url = \"http://geoportal.statistics.gov.uk/datasets/055c2d8135ca4297a85d624bb68aefdb_0.csv\"\n",
    "\n",
    "with Session() as session:\n",
    "    response = session.get(postcode_url, verify = False)\n",
    "\n",
    "field_dtypes = {'objectid': 'int32', 'pcd':'str', 'pcd2': 'str', 'pcds':'str', 'dointr':'str','doterm':'str',\n",
    "                'usertype':'int8','oseast1m': 'float', 'osnorth1m': 'float', 'osgrdind':'int8', 'lat':'float', \n",
    "                'long':'float', 'X':'float', 'Y':'float', 'imd': 'int8',\n",
    "                'oa11':'str', 'cty': 'str', 'ced':'str', 'laua': 'str', 'ward': 'str', 'hlthau':'str',\n",
    "                'ctry': 'str','pcon': 'str','eer': 'str','teclec': 'str','ttwa': 'str','pct': 'str','nuts': 'str',\n",
    "                'park': 'str','lsoa11': 'str','msoa11': 'str','wz11': 'str','ccg': 'str','bua11': 'str',\n",
    "                'buasd11': 'str','ru11ind': 'str','oac11': 'str','lep1': 'str','lep2': 'str','pfa': 'str',\n",
    "                'ced': 'str','nhser': 'str','rgn': 'str','calncv': 'str','stp': 'str'}\n",
    "\n",
    "pc = pd.read_csv(BytesIO(response.content), dtype = field_dtypes)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create pcstrip for matching\n",
    "pc['pcstrip'] = pc['pcd'].str.replace(\"\\s\",\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NB - here I'm joining 2011 LSOA, 2011 MSOA, Rural-Urban Indicator, Region (formally GOR), Local Authority Area (effectively district), and IMD score (NB separate for each country). However, you can add any of the geography codes available in the NSPL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_temp = data.merge(pc[['pcstrip','oseast1m','osnrth1m','lsoa11','msoa11','ru11ind','rgn','laua','imd']], \n",
    "                       how = 'left',\n",
    "                       on = 'pcstrip')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CF34LG     12\n",
       "LD3OAW     12\n",
       "OO000OO     1\n",
       "Name: pcstrip, dtype: int64"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check for missing postcodes\n",
    "data_temp[data_temp['oseast1m'].isnull()]['pcstrip'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clean Missing Postcodes - appear to be typos.\n",
    "new_pcs = {'CF34LG': 'CF33LG',    # Rumney Primary Care Centre\n",
    "           'LD3OAW': 'LD30AW',    # Haygarth Doctors\n",
    "           'OO000OO': 'CF433HB'}  # \"THE SURGERY, PENRHYS, FERNDALE RHONDDA\"(Tylorstown Surgery, in real life)\n",
    "\n",
    "data['pcstrip'] = data['pcstrip'].map(new_pcs).fillna(data['pcstrip'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merge data\n",
    "data = data.merge(pc[['pcstrip','oseast1m','osnrth1m','lsoa11','msoa11','ru11ind','rgn','laua','imd']], \n",
    "                  how = 'left', \n",
    "                  on = 'pcstrip')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save aggregated data\n",
    "data.to_csv(path + \"processed_data_with_postcodes.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check Postcodes \n",
    "\n",
    "Some Practice Codes have more than one postcode associated with them. Possible reasons for this are:\n",
    "* Practices move to a new location.\n",
    "* Practices are assigned a new postcode but don't physically move.\n",
    "* Practice postcodes are wrongly entered at a particular wave and subsequently fixed.\n",
    "\n",
    "There are 10 practices codes which have more than 1 postcode assigned to them, this is about 1.7% of unique practices.\n",
    "\n",
    "All 10 practices have 2 postcodes associated with them.\n",
    "\n",
    "2 of these practices fall within the same LSOA, 8 fall within 2 different LSOAs.\n",
    "\n",
    "**We'll ignore this for now** - this will require some more advanced cleaning - useful to be aware of though."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PracticeID                 object\n",
       "depression                  int64\n",
       "alzheimers                  int64\n",
       "blood pressure              int64\n",
       "hypertension                int64\n",
       "diabeties                   int64\n",
       "cardiovascular disease      int64\n",
       "insomnia                    int64\n",
       "addiction                   int64\n",
       "social anxiety              int64\n",
       "Items                       int64\n",
       "loneliness                  int64\n",
       "Date                        int64\n",
       "HB                         object\n",
       "pcstrip                    object\n",
       "Street                     object\n",
       "Area                       object\n",
       "Posttown                   object\n",
       "Postcode                   object\n",
       "oseast1m                  float64\n",
       "osnrth1m                  float64\n",
       "lsoa11                     object\n",
       "msoa11                     object\n",
       "ru11ind                    object\n",
       "rgn                        object\n",
       "laua                       object\n",
       "imd                          int8\n",
       "dtype: object"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check for 1 postcode per Practice Code\n",
    "pc_prac_counts = data.groupby('PracticeID')['pcstrip'].unique().map(len)\n",
    "# 10 practices have more than 1 postcode associated with it.\n",
    "pc_prac_counts[pc_prac_counts > 1].count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.017889087656529516"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Practice codes with multiple associated postcodes account forc. 1.7% of the data\n",
    "pc_prac_counts[pc_prac_counts > 1].count()/len(pc_prac_counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2    10\n",
       "Name: pcstrip, dtype: int64"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pc_prac_counts[pc_prac_counts > 1].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2    8\n",
       "1    2\n",
       "Name: lsoa11, dtype: int64"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Count of these Practices falling within the same LSOA\n",
    "(data[data['PracticeID'].isin(pc_prac_counts[pc_prac_counts > 1].index)]\n",
    " .groupby('PracticeID')['lsoa11']\n",
    " .unique()\n",
    " .map(len)\n",
    " .value_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Subsetting the Data\n",
    "\n",
    "## Use only General Practice surgeries\n",
    "\n",
    "Use the 'GP Practice Analysis' data from: http://www.primarycareservices.wales.nhs.uk/gp-practice-analysis to get GP surgery codes and number of registered patients.\n",
    "\n",
    "Note that Wales publishes only one file per year, so we'll just have to use the end-of-year patient counts rather than month-by-month counts.\n",
    "\n",
    "**Filename format:** PracticeItems[month][year].xlsx (e.g. PracticeItemsDecember2018.xlsx)\n",
    "\n",
    "**Columns:** LHBName, PracticeID, CurrentRegisteredGPsInPractice, SeniorPartner, Street, Area, PostTown, County, Postcode, Dispensing, Number of Prescribing Patients, Number of Dispensing Patients, Total Number of Patients (Including Temporary Residents), Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, Grand Total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(path + \"processed_data_with_postcodes.csv\", index_col = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>PracticeID</th>\n",
       "      <th>Number_of_Patients</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>W93007</td>\n",
       "      <td>4953</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>W93009</td>\n",
       "      <td>4531</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>W93032</td>\n",
       "      <td>8570</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>W93035</td>\n",
       "      <td>7093</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>W93063</td>\n",
       "      <td>5589</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  PracticeID  Number_of_Patients\n",
       "0     W93007                4953\n",
       "1     W93009                4531\n",
       "2     W93032                8570\n",
       "3     W93035                7093\n",
       "4     W93063                5589"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get GP files\n",
    "gp_path = path + r\"GP data\\\\\"\n",
    "file = \"PracticeItemsDecember2018.xlsx\"\n",
    "\n",
    "gp_data = pd.read_excel(gp_path + file)\n",
    "gp_data.rename(columns = {'Total Number of Patients (Including Temporary Residents)': 'Number_of_Patients'}, inplace = True)\n",
    "\n",
    "gp_data = gp_data[['PracticeID', 'Number_of_Patients']]\n",
    "\n",
    "gp_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6009, 27)"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get the unique codes for GP surgeries and subset the prescribing data according to these codes.\n",
    "gp_ids = gp_data['PracticeID'].unique()\n",
    "data = data[data['PracticeID'].isin(gp_ids)].copy()\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>PracticeID</th>\n",
       "      <th>depression</th>\n",
       "      <th>alzheimers</th>\n",
       "      <th>blood pressure</th>\n",
       "      <th>hypertension</th>\n",
       "      <th>diabeties</th>\n",
       "      <th>cardiovascular disease</th>\n",
       "      <th>insomnia</th>\n",
       "      <th>addiction</th>\n",
       "      <th>social anxiety</th>\n",
       "      <th>...</th>\n",
       "      <th>Postcode</th>\n",
       "      <th>oseast1m</th>\n",
       "      <th>osnrth1m</th>\n",
       "      <th>lsoa11</th>\n",
       "      <th>msoa11</th>\n",
       "      <th>ru11ind</th>\n",
       "      <th>rgn</th>\n",
       "      <th>laua</th>\n",
       "      <th>imd</th>\n",
       "      <th>Number_of_Patients</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>W00005</td>\n",
       "      <td>1542</td>\n",
       "      <td>17</td>\n",
       "      <td>1791</td>\n",
       "      <td>1807</td>\n",
       "      <td>329</td>\n",
       "      <td>1448</td>\n",
       "      <td>99</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>...</td>\n",
       "      <td>CF83 4AZ</td>\n",
       "      <td>312151.0</td>\n",
       "      <td>189346.0</td>\n",
       "      <td>W01001331</td>\n",
       "      <td>W02000307</td>\n",
       "      <td>D1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000018</td>\n",
       "      <td>9</td>\n",
       "      <td>8031.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>W00007</td>\n",
       "      <td>1006</td>\n",
       "      <td>15</td>\n",
       "      <td>949</td>\n",
       "      <td>955</td>\n",
       "      <td>186</td>\n",
       "      <td>730</td>\n",
       "      <td>76</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "      <td>...</td>\n",
       "      <td>LL11 3NS</td>\n",
       "      <td>328624.0</td>\n",
       "      <td>351132.0</td>\n",
       "      <td>W01000364</td>\n",
       "      <td>W02000083</td>\n",
       "      <td>D1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000006</td>\n",
       "      <td>104</td>\n",
       "      <td>7216.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>W00019</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>CH5 1XS</td>\n",
       "      <td>330760.0</td>\n",
       "      <td>367615.0</td>\n",
       "      <td>W01000335</td>\n",
       "      <td>W02000066</td>\n",
       "      <td>C1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000005</td>\n",
       "      <td>29</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>W00021</td>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>NP44 7XX</td>\n",
       "      <td>329326.0</td>\n",
       "      <td>193774.0</td>\n",
       "      <td>W01001504</td>\n",
       "      <td>W02000335</td>\n",
       "      <td>C1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000020</td>\n",
       "      <td>3</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>W00022</td>\n",
       "      <td>16</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>NP44 7XX</td>\n",
       "      <td>329326.0</td>\n",
       "      <td>193774.0</td>\n",
       "      <td>W01001504</td>\n",
       "      <td>W02000335</td>\n",
       "      <td>C1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000020</td>\n",
       "      <td>3</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 28 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "  PracticeID  depression  alzheimers  blood pressure  hypertension  diabeties  \\\n",
       "0     W00005        1542          17            1791          1807        329   \n",
       "1     W00007        1006          15             949           955        186   \n",
       "2     W00019           1           0               3             3          0   \n",
       "3     W00021           8           0               3             3          1   \n",
       "4     W00022          16           0               2             2          2   \n",
       "\n",
       "   cardiovascular disease  insomnia  addiction  social anxiety  ...  Postcode  \\\n",
       "0                    1448        99          0               3  ...  CF83 4AZ   \n",
       "1                     730        76          0              11  ...  LL11 3NS   \n",
       "2                       1         1          0               0  ...   CH5 1XS   \n",
       "3                       0         4          0               0  ...  NP44 7XX   \n",
       "4                       2         2          0               0  ...  NP44 7XX   \n",
       "\n",
       "   oseast1m  osnrth1m     lsoa11     msoa11 ru11ind        rgn       laua  \\\n",
       "0  312151.0  189346.0  W01001331  W02000307      D1  W99999999  W06000018   \n",
       "1  328624.0  351132.0  W01000364  W02000083      D1  W99999999  W06000006   \n",
       "2  330760.0  367615.0  W01000335  W02000066      C1  W99999999  W06000005   \n",
       "3  329326.0  193774.0  W01001504  W02000335      C1  W99999999  W06000020   \n",
       "4  329326.0  193774.0  W01001504  W02000335      C1  W99999999  W06000020   \n",
       "\n",
       "   imd  Number_of_Patients  \n",
       "0    9              8031.0  \n",
       "1  104              7216.0  \n",
       "2   29                 0.0  \n",
       "3    3                 0.0  \n",
       "4    3                 0.0  \n",
       "\n",
       "[5 rows x 28 columns]"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Merge on the basis of PracticeID - some patient counts are missing (zero)\n",
    "data = data.merge(gp_data, how = 'left', on = ['PracticeID'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save aggregated data\n",
    "data.to_csv(path + \"processed_data_with_postcodes_GPs.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate Statistics from Prescribing Counts\n",
    "\n",
    "## Percentages At Postcode Level\n",
    "\n",
    "Aggregate observations to postcodes and compute percentages for 'depression', 'alzheimers', 'blood pressure', 'hypertension', 'diabetes', 'cardiovascular disease', 'insomnia', 'addiction', 'social anxiety', and 'loneliness'.\n",
    "\n",
    "## Outlier Removal\n",
    "\n",
    "Should we remove some GPs on the basis that they have very low/high values which might indicate they are not accessible to the general population, and instead represent specialist services?\n",
    "\n",
    "Currently, we won't do this, but it's an advanced task to look into.\n",
    "\n",
    "## Standardise Percentages\n",
    "\n",
    "Should we standardise within time points, or across them? Or standardise with GP practices or across them?\n",
    "\n",
    "Can't standardize within GPs, as can't then compare between GPs.\n",
    "\n",
    "Can't standardise across GPs within years, as can't then compare between years.\n",
    "\n",
    "Can't standardise across GPs across years, as can't then disambiguate changes to rank order over time.\n",
    "\n",
    "<u> First Step </u>\n",
    "\n",
    "Take the average percentage of disease groups within postcodes annually - this is then an annual summary measure of prescribing by postcode. Aggregation entire depends on desired time frame for analysis.\n",
    "\n",
    "NB, this is a mean of percentages - could also calculate an overall percentage by summing monthly counts by year and dividing through by monthly sum of items.\n",
    "\n",
    "<u> Second Step </u>\n",
    "\n",
    "z-score standardise for earliest year observed across GPs. Store mean and standard deviation.\n",
    "\n",
    "z-score standardise subsequent years wrt baseline mean and standard deviation.\n",
    "\n",
    "OR\n",
    "\n",
    "Use min-max normalisation by year (decile normalisation?). This standardises the different percentages to the same range - although in theory they're comparable anyway...\n",
    "\n",
    "## Aggregation and Percentages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = pd.read_csv(path + \"processed_data_with_postcodes_GPs.csv\", index_col = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>PracticeID</th>\n",
       "      <th>depression</th>\n",
       "      <th>alzheimers</th>\n",
       "      <th>blood pressure</th>\n",
       "      <th>hypertension</th>\n",
       "      <th>diabeties</th>\n",
       "      <th>cardiovascular disease</th>\n",
       "      <th>insomnia</th>\n",
       "      <th>addiction</th>\n",
       "      <th>social anxiety</th>\n",
       "      <th>...</th>\n",
       "      <th>Postcode</th>\n",
       "      <th>oseast1m</th>\n",
       "      <th>osnrth1m</th>\n",
       "      <th>lsoa11</th>\n",
       "      <th>msoa11</th>\n",
       "      <th>ru11ind</th>\n",
       "      <th>rgn</th>\n",
       "      <th>laua</th>\n",
       "      <th>imd</th>\n",
       "      <th>Number_of_Patients</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>W00007</td>\n",
       "      <td>1006</td>\n",
       "      <td>15</td>\n",
       "      <td>949</td>\n",
       "      <td>955</td>\n",
       "      <td>186</td>\n",
       "      <td>730</td>\n",
       "      <td>76</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "      <td>...</td>\n",
       "      <td>LL11 3NS</td>\n",
       "      <td>328624.0</td>\n",
       "      <td>351132.0</td>\n",
       "      <td>W01000364</td>\n",
       "      <td>W02000083</td>\n",
       "      <td>D1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000006</td>\n",
       "      <td>104</td>\n",
       "      <td>7216.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>W00019</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>CH5 1XS</td>\n",
       "      <td>330760.0</td>\n",
       "      <td>367615.0</td>\n",
       "      <td>W01000335</td>\n",
       "      <td>W02000066</td>\n",
       "      <td>C1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000005</td>\n",
       "      <td>29</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>W00021</td>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>NP44 7XX</td>\n",
       "      <td>329326.0</td>\n",
       "      <td>193774.0</td>\n",
       "      <td>W01001504</td>\n",
       "      <td>W02000335</td>\n",
       "      <td>C1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000020</td>\n",
       "      <td>3</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>W00022</td>\n",
       "      <td>16</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>NP44 7XX</td>\n",
       "      <td>329326.0</td>\n",
       "      <td>193774.0</td>\n",
       "      <td>W01001504</td>\n",
       "      <td>W02000335</td>\n",
       "      <td>C1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000020</td>\n",
       "      <td>3</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>W00023</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>NP44 7XX</td>\n",
       "      <td>329326.0</td>\n",
       "      <td>193774.0</td>\n",
       "      <td>W01001504</td>\n",
       "      <td>W02000335</td>\n",
       "      <td>C1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000020</td>\n",
       "      <td>3</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>W00024</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>NP44 7XX</td>\n",
       "      <td>329326.0</td>\n",
       "      <td>193774.0</td>\n",
       "      <td>W01001504</td>\n",
       "      <td>W02000335</td>\n",
       "      <td>C1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000020</td>\n",
       "      <td>3</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>W00030</td>\n",
       "      <td>67</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>15</td>\n",
       "      <td>62</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>SA72 6SY</td>\n",
       "      <td>195752.0</td>\n",
       "      <td>203509.0</td>\n",
       "      <td>W01000603</td>\n",
       "      <td>W02000138</td>\n",
       "      <td>D1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000009</td>\n",
       "      <td>-65</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>W00031</td>\n",
       "      <td>10</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>LL18 5UJ</td>\n",
       "      <td>300309.0</td>\n",
       "      <td>376019.0</td>\n",
       "      <td>W01000191</td>\n",
       "      <td>W02000050</td>\n",
       "      <td>D1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000004</td>\n",
       "      <td>-126</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>W00032</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>LL18 5UJ</td>\n",
       "      <td>300309.0</td>\n",
       "      <td>376019.0</td>\n",
       "      <td>W01000191</td>\n",
       "      <td>W02000050</td>\n",
       "      <td>D1</td>\n",
       "      <td>W99999999</td>\n",
       "      <td>W06000004</td>\n",
       "      <td>-126</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>9 rows × 28 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "  PracticeID  depression  alzheimers  blood pressure  hypertension  diabeties  \\\n",
       "1     W00007        1006          15             949           955        186   \n",
       "2     W00019           1           0               3             3          0   \n",
       "3     W00021           8           0               3             3          1   \n",
       "4     W00022          16           0               2             2          2   \n",
       "5     W00023           0           0               0             0          0   \n",
       "6     W00024           0           0               0             0          0   \n",
       "7     W00030          67           0               0             0          0   \n",
       "8     W00031          10           0               1             1          2   \n",
       "9     W00032           0           0               0             0          0   \n",
       "\n",
       "   cardiovascular disease  insomnia  addiction  social anxiety  ...  Postcode  \\\n",
       "1                     730        76          0              11  ...  LL11 3NS   \n",
       "2                       1         1          0               0  ...   CH5 1XS   \n",
       "3                       0         4          0               0  ...  NP44 7XX   \n",
       "4                       2         2          0               0  ...  NP44 7XX   \n",
       "5                       0         1          0               0  ...  NP44 7XX   \n",
       "6                       0         0          0               0  ...  NP44 7XX   \n",
       "7                       0        15         62               0  ...  SA72 6SY   \n",
       "8                       0         1          0               0  ...  LL18 5UJ   \n",
       "9                       0         0          0               0  ...  LL18 5UJ   \n",
       "\n",
       "   oseast1m  osnrth1m     lsoa11     msoa11 ru11ind        rgn       laua  \\\n",
       "1  328624.0  351132.0  W01000364  W02000083      D1  W99999999  W06000006   \n",
       "2  330760.0  367615.0  W01000335  W02000066      C1  W99999999  W06000005   \n",
       "3  329326.0  193774.0  W01001504  W02000335      C1  W99999999  W06000020   \n",
       "4  329326.0  193774.0  W01001504  W02000335      C1  W99999999  W06000020   \n",
       "5  329326.0  193774.0  W01001504  W02000335      C1  W99999999  W06000020   \n",
       "6  329326.0  193774.0  W01001504  W02000335      C1  W99999999  W06000020   \n",
       "7  195752.0  203509.0  W01000603  W02000138      D1  W99999999  W06000009   \n",
       "8  300309.0  376019.0  W01000191  W02000050      D1  W99999999  W06000004   \n",
       "9  300309.0  376019.0  W01000191  W02000050      D1  W99999999  W06000004   \n",
       "\n",
       "   imd  Number_of_Patients  \n",
       "1  104              7216.0  \n",
       "2   29                 0.0  \n",
       "3    3                 0.0  \n",
       "4    3                 0.0  \n",
       "5    3                 0.0  \n",
       "6    3                 0.0  \n",
       "7  -65                 0.0  \n",
       "8 -126                 0.0  \n",
       "9 -126                 0.0  \n",
       "\n",
       "[9 rows x 28 columns]"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[1:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make dictionary for aggregation\n",
    "# counts to sum\n",
    "agg_cols = {col : 'sum' for col in drug_data['illness'].unique()}\n",
    "agg_cols[col_items] = 'sum'\n",
    "agg_cols['loneliness'] = 'sum'\n",
    "agg_cols['Number_of_Patients'] = 'sum'\n",
    "\n",
    "# Other data to preserve\n",
    "for key in ['HB','Street','Area','Posttown','Postcode',\n",
    "            'oseast1m', 'osnrth1m', 'lsoa11', 'msoa11','ru11ind', 'rgn', 'laua', 'imd']:\n",
    "    agg_cols[key] = 'first'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do aggregation and produce counts by postcode by date.\n",
    "data = data.groupby(['pcstrip','Date'], as_index=False).agg(agg_cols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate percentages\n",
    "perc_cols = drug_data['illness'].unique()\n",
    "target_cols = perc_cols + '_perc'\n",
    "\n",
    "# Percentages for discrete illness groups\n",
    "data[target_cols] = data[perc_cols].divide(data[col_items], axis=0) * 100\n",
    "\n",
    "# Overall percentage for loneliness realted disease prescribing\n",
    "data['loneliness_perc'] = data['loneliness'].divide(data[col_items], axis=0) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Standardisation (z-scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "pcstrip                         object\n",
       "Date                             int64\n",
       "depression                       int64\n",
       "alzheimers                       int64\n",
       "blood pressure                   int64\n",
       "hypertension                     int64\n",
       "diabeties                        int64\n",
       "cardiovascular disease           int64\n",
       "insomnia                         int64\n",
       "addiction                        int64\n",
       "social anxiety                   int64\n",
       "Items                            int64\n",
       "loneliness                       int64\n",
       "Number_of_Patients             float64\n",
       "HB                              object\n",
       "Street                          object\n",
       "Area                            object\n",
       "Posttown                        object\n",
       "Postcode                        object\n",
       "oseast1m                       float64\n",
       "osnrth1m                       float64\n",
       "lsoa11                          object\n",
       "msoa11                          object\n",
       "ru11ind                         object\n",
       "rgn                             object\n",
       "laua                            object\n",
       "imd                              int64\n",
       "depression_perc                float64\n",
       "alzheimers_perc                float64\n",
       "blood pressure_perc            float64\n",
       "hypertension_perc              float64\n",
       "diabeties_perc                 float64\n",
       "cardiovascular disease_perc    float64\n",
       "insomnia_perc                  float64\n",
       "addiction_perc                 float64\n",
       "social anxiety_perc            float64\n",
       "loneliness_perc                float64\n",
       "dtype: object"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Firstly aggregate percentages by postcodes by year.\n",
    "#data['Year'] = data['Date'].dt.year\n",
    "data['Year'] = 2018\n",
    "\n",
    "# Aggregation\n",
    "cols = {'Number_of_Patients': 'mean', 'HB': 'first', 'oseast1m': 'first', 'osnrth1m': 'first',\n",
    "        'lsoa11': 'first', 'msoa11': 'first', 'ru11ind': 'first', 'rgn': 'first', 'laua':'first', 'imd': 'first',\n",
    "        'depression_perc': 'mean', 'alzheimers_perc': 'mean', 'blood pressure_perc': 'mean', 'hypertension_perc': 'mean',\n",
    "        'diabeties_perc': 'mean', 'cardiovascular disease_perc': 'mean', 'insomnia_perc': 'mean', 'addiction_perc': 'mean',\n",
    "        'social anxiety_perc': 'mean', 'loneliness_perc': 'mean'}\n",
    "\n",
    "data = data.groupby(['pcstrip','Year'], as_index=False).agg(cols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The mean value returns a value broadly in the centre of the distribution of respective disease classes.\n",
    "# Therefore we'll go with an un-truncated arithmetic mean.\n",
    "# Can always revisit this assumption later.\n",
    "\n",
    "per_cols = ['depression_perc', 'alzheimers_perc', 'blood pressure_perc', 'hypertension_perc', \n",
    "            'diabeties_perc', 'cardiovascular disease_perc', 'insomnia_perc', 'addiction_perc',\n",
    "            'social anxiety_perc', 'loneliness_perc']\n",
    "\n",
    "# Get mean and std for baseline (2016)\n",
    "mean_std = data[data['Year'] == 2018][per_cols].agg(['mean','std'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make new column names.\n",
    "std_cols = [col[:-4] + 'zscore' for col in per_cols]\n",
    "\n",
    "zscores = []    \n",
    "# z-score standardise for each year by baseline mean and std \n",
    "for year in [2016,2017,2018]:\n",
    "    zscores.append((data.loc[data['Year'] == year, per_cols] - mean_std.loc['mean', per_cols]) / mean_std.loc['std', per_cols])\n",
    "\n",
    "zscores = pd.concat(zscores).sort_index()\n",
    "data[std_cols] = zscores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzIAAAFpCAYAAABH4x0wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAF0dJREFUeJzt3X+snXd9H/D3p7gZPwyEEjDIybigpqys0Qa2EIits5uqpAQR/oANNFjapYv6Y4wNpuKOP0CbpsJYC1RCraqkXbrSGpr+AMW0rGX2ulVNWgcYAdKMNEUQlkFRSVozNpr1uz/ucbhxbJ9z7XPvcz7Xr5d05POc+73H7+d57vO1336e57jGGAEAAOjkm6YOAAAAsFmKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANDOru38zS655JKxtrY2d9xXv/rVPO5xj9v6QOdJzuXqkLNDxmRzOW+//fYvjzGessWRls58Mo0OOTtkTHZeTnPJ6rMOq8E6zLfwfDLG2LbHvn37xiKOHj260LipyblcHXJ2yDjG5nImOT62cR5Y1sN8Mo0OOTtkHGPn5TSXrD7rsBqsw3yLzicuLQMAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2lFkAACAdhQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAALiBrh45k7dCRqWPAeVNkAACAdhQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2lFkAACAdhQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaGehIlNV/6KqPlVVn6yqX66qR1fVM6vqtqr6TFW9r6ou2uqwAAAAyQJFpqr2JvlnSfaPMb4jyaOSvCrJ25O8c4xxeZKvJLluK4MCAACctOilZbuSPKaqdiV5bJL7knxXkptnX78pycuXHw8AAOCR5haZMcYXkvz7JJ/LeoF5IMntSe4fYzw4G3Zvkr1bFRIAAGCjGmOcfUDVk5L8apJ/kOT+JL8yW37LGONbZ2MuS/KhMcYVp/n+65NcnyR79uzZd/jw4bmhTpw4kd27d29uTSYg53J1yNkhY7K5nAcPHrx9jLF/iyMthflkeh1ydsiY7Lyc5pLVd3Id7vjCA0mSK/Y+ceJEm7eT9kNnW70OC88nY4yzPpK8MsmNG5b/UZKfTvLlJLtmr70wyYfnvde+ffvGIo4ePbrQuKnJuVwdcnbIOMbmciY5PuYcu6v4MJ9Mo0PODhnH2Hk5zSWr7+Q6PONNt4xnvOmWacOco520Hzrb6nVYdD5Z5B6ZzyV5QVU9tqoqyZVJPp3kaJJXzMZcm+QDC9csAACA87DIPTK3Zf2m/o8muWP2PT+b5E1J3lBVdyd5cpIbtzAnAADAQ3YtMmiM8ZYkbznl5XuSPH/piQAAAOZY9OOXAQAAVoYiAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAHABuOMLD2Tt0JGpY8DSKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SxUZKrq4qq6uar+qKrurKoXVtW3VNVvV9VnZr8+aavDAgAAJIufkXl3kt8aY/yNJH8ryZ1JDiX5yBjj8iQfmS0DAABsublFpqqekOQ7k9yYJGOMr48x7k9yTZKbZsNuSvLyrQoJAACw0SJnZJ6V5E+T/HxVfayqbqiqxyXZM8a4L0lmvz51C3MCAAA8pMYYZx9QtT/JrUleNMa4rareneTPk7xujHHxhnFfGWM84j6Zqro+yfVJsmfPnn2HDx+eG+rEiRPZvXv3plZkCnIuV4ecHTImm8t58ODB28cY+7c40lKYT6bXIWeHjMnOy2kuWX1f+rMH8sWvfWP5ir1PnC7MOdoJ+8E6zLfwfDLGOOsjydOSfHbD8t9NciTJXUmePnvt6Unumvde+/btG4s4evToQuOmJudydcjZIeMYm8uZ5PiYc+yu4sN8Mo0OOTtkHGPn5TSXrL6f+sXfGM940y0PPTraCfvBOsy36Hwy99KyMcb/SvL5qnr27KUrk3w6yQeTXDt77dokH1i0ZQEAAJyPXQuOe12S91bVRUnuSfL9Wb+/5v1VdV2SzyV55dZEBAAAeLiFiswY4+NJTned2pXLjQMAADDfov+PDAAAwMpQZAAAgHYUGQCAC9DaoSNZO3Rk6hhwzhQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2lFkAACAdhQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHZ2TR0AAICts3boSJLkjVdMHASWzBkZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2lFkAACAdhQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYA4AK2duhI1g4dmToGbJoiAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANDOwkWmqh5VVR+rqltmy8+sqtuq6jNV9b6qumjrYgIAAHzDZs7IvD7JnRuW357knWOMy5N8Jcl1ywwGAABwJgsVmaq6NMnVSW6YLVeS70py82zITUlevhUBAQAATrXoGZl3JfnRJH81W35ykvvHGA/Olu9NsnfJ2QAAAE6rxhhnH1D10iQvGWP8cFUdSPIvk3x/kt8fY3zrbMxlST40xrjiNN9/fZLrk2TPnj37Dh8+PDfUiRMnsnv37k2uyvaTc7k65OyQMdlczoMHD94+xti/xZGWwnwyvQ45O2RMdl5Oc8nquuMLDyRJ9jwm+eLXzjzuir1P3KZE567zfjjJOsy36HyySJH58SSvTfJgkkcneUKSX0/y4iRPG2M8WFUvTPLWMcaLz/Ze+/fvH8ePH58b/tixYzlw4MDccVOTc7k65OyQMdlczqpq85ePjcwn0+iQs0PGZOflNJesrrVDR5Ikb7ziwfzEHbvOOO6zb7t6uyKds8774STrMN+i88ncS8vGGD82xrh0jLGW5FVJ/vMY4x8mOZrkFbNh1yb5wHnkBQAAWNj5/D8yb0ryhqq6O+v3zNy4nEgAAABnd+bzi6cxxjiW5Njs+T1Jnr/8SAAAAGd3PmdkAAAAJqHIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAAALWDt0JGuHjkwdA5hRZAAAgHYUGQAAoB1FBgAAaEeRAQDAPUC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO3smjoAAADLt3boyNQRYEs5IwMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANDO3CJTVZdV1dGqurOqPlVVr5+9/i1V9dtV9ZnZr0/a+rgAAACLnZF5MMkbxxjfnuQFSX6kqp6T5FCSj4wxLk/ykdkyAADAlptbZMYY940xPjp7/hdJ7kyyN8k1SW6aDbspycu3KiQAAMBGm7pHpqrWkjw3yW1J9owx7kvWy06Spy47HAAAwOnUGGOxgVW7k/yXJP92jPFrVXX/GOPiDV//yhjjEffJVNX1Sa5Pkj179uw7fPjw3N/rxIkT2b1794KrMB05l6tDzg4Zk83lPHjw4O1jjP1bHGkpzCfT65CzQ8Zk5+W8EOaSP3ng/yVJrtj7xC3Ntyx3fOGBhy3veUzyxa/N/75VXr8ux83ZWIf5Fp1PFioyVfXNSW5J8uExxk/OXrsryYExxn1V9fQkx8YYzz7b++zfv38cP3587u937NixHDhwYO64qcm5XB1ydsiYbC5nVbX5y8dG5pNpdMjZIWOy83JeCHPJ9/3WV5Mkn33b1VsdaynWDh152PIbr3gwP3HHrrnft8rr1+W4ORvrMN+i88kin1pWSW5McufJEjPzwSTXzp5fm+QD5xIUAABgs+bX8uRFSV6b5I6q+vjstX+V5G1J3l9V1yX5XJJXbk1EAACAh5tbZMYY/y1JneHLVy43DgAAwHyb+tQyAACAVaDIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAA8JC1Q0eydujI1DFgLkUGAABoR5EBAADaUWQAAIB2dk0dAACA5XBvCxcSZ2QAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2lFkAACAdhQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2lFkAADOwdqhI1k7dGTqGHDBUmQAAIB2FBkAAKAdRQYAAGhn19QBAAA4P+7V4ULkjAwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAADwCGuHjvhYZ1aaIgMAALSjyAAAAO3smjoAAEAnp15udXL5s2+7eoo4cMFyRgYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2vHxywAAS+BjmOGRNn5c+bKPDWdkAACAdhQZAACgHUUGAABoxz0yAABb4NR7ZpZ5r8B23o/j3h9WlTMyAABAO4oMAADQznldWlZVVyV5d5JHJblhjPG2paQCAGhq4yVkp1ve+Nqpl2vNe33e77WVtvLSuHN57zO9x7zls/0e895jM07dN4vu682892a/d6ddJnjOZ2Sq6lFJ3pPke5M8J8mrq+o5ywoGAABwJudzadnzk9w9xrhnjPH1JIeTXLOcWAAAAGd2PkVmb5LPb1i+d/YaAADAlqoxxrl9Y9Urk7x4jPEDs+XXJnn+GON1p4y7Psn1s8VnJ7lrgbe/JMmXzynY9pJzuTrk7JAx2VzOZ4wxnrKVYZbFfLISOuTskDHZeTnNJavPOqwG6zDfQvPJ+RSZFyZ56xjjxbPlH0uSMcaPn9MbPvy9j48x9p/v+2w1OZerQ84OGZM+ObdDl20h5/J0yJjI2c1O2A7WYTVYh+U5n0vL/jDJ5VX1zKq6KMmrknxwObEAAADO7Jw/fnmM8WBV/dMkH876xy//3BjjU0tLBgAAcAbn9f/IjDE+lORDS8qy0c9uwXtuBTmXq0PODhmTPjm3Q5dtIefydMiYyNnNTtgO1mE1WIclOed7ZAAAAKZyPvfIAAAATGLSIlNVV1XVXVV1d1UdOs3X/1pVvW/29duqam37Uy6U8w1V9emq+kRVfaSqnrGKOTeMe0VVjara9k+bWCRjVf392fb8VFX90nZnnGWYt8//elUdraqPzfb7SybI+HNV9aWq+uQZvl5V9VOzdfhEVT1vuzNupw7ziblkucwnS81oPllAVb2jqv5otg1+vaounjrTIhY9pldVVV02O0bunB3Lr58607mqqkfNjvVbps5yLqrq4qq6eXYc3Fnrn2I8nTHGJI+sf0DAHyd5VpKLkvz3JM85ZcwPJ/mZ2fNXJXnfiuY8mOSxs+c/tKo5Z+Men+R3k9yaZP+qZUxyeZKPJXnSbPmpq7gts35t6A/Nnj8nyWcnyPmdSZ6X5JNn+PpLkvxmkkrygiS3bXfGFdtnk84n5pJJtqf5ZPGc5pPFttP3JNk1e/72JG+fOtMCmRc6plf5keTpSZ43e/74JP+j2zpsWJc3JPmlJLdMneUc89+U5Admzy9KcvGUeaY8I/P8JHePMe4ZY3w9yeEk15wy5pqsb7AkuTnJlVVV25gxWSDnGOPoGON/zxZvTXLpNmdMFtueSfJvkvy7JP9nO8PNLJLxnyR5zxjjK0kyxvjSNmdMFss5kjxh9vyJSf7nNuZbDzDG7yb5s7MMuSbJL4x1tya5uKqevj3ptl2H+cRcslzmkyUynyxmjPGfxhgPzhanOkY3a9FjemWNMe4bY3x09vwvktyZZO+0qTavqi5NcnWSG6bOci6q6glZ/0ePG5NkjPH1Mcb9U2aassjsTfL5Dcv35pE/lA+NmU0cDyR58rakO02GmdPl3Oi6rP+r1Xabm7OqnpvksjHGVKczF9mW35bk26rq96rq1qq6atvSfcMiOd+a5DVVdW/WP7nvddsTbVM2+7PbWYf5xFyyXOaT7XUhzSeL+seZ5hjdrB2172aXBT83yW3TJjkn70ryo0n+auog5+hZSf40yc/PLo+7oaoeN2WgKYvM6f4l9NSPUFtkzFZbOENVvSbJ/iTv2NJEp3fWnFX1TUnemeSN25bokRbZlruyfjnIgSSvTnLDBNcgL5Lz1Un+wxjj0qxfcvEfZ9t4lazC8bNdOswn5pLlMp9sr6mPn21TVb9TVZ88zeOaDWPenOTBJO+dLunCdsy+q6rdSX41yT8fY/z51Hk2o6pemuRLY4zbp85yHnZl/RLUnx5jPDfJV5NMes/Vef0/Mufp3iSXbVi+NI88nX5yzL1VtSvrp9zPdup7KyySM1X13UnenOTvjTH+7zZl22hezscn+Y4kx2ZX0zwtyQer6mVjjOMrkvHkmFvHGH+Z5E+q6q6s/0XkD7cn4kMZ5uW8LslVSTLG+P2qenSSS5JMcenKmSz0s7tDdJhPzCXLZT7ZXhfMfDLG+O6zfb2qrk3y0iRXjtmNAituR+y7qvrmrJeY944xfm3qPOfgRUleNvswj0cneUJV/eIY4zUT59qMe5PcO8Y4eTbs5kxcZKa8WWhXknuSPDPfuPnsb54y5kfy8Jtz37+iOZ+b9RvpLl/l7XnK+GPZ/pv9F9mWVyW5afb8kqyfDn/yCub8zSTfN3v+7Vn/Q6Em2O9rOfPNuVfn4Tfn/sF251uxfTbpfGIumWR7mk82l9V8Mn8bXZXk00meMnWWTWTe1DG9io/Zz90vJHnX1FmWtD4H0vdm//+a5Nmz529N8o5J80y8MV6S9U+e+OMkb5699q+TvGz2/NFJfiXJ3Un+IMmzVjTn7yT5YpKPzx4fXMWcp4yd6i8f87ZlJfnJ2R8UdyR51Spuy6x/stDvzf5A+HiS75kg4y8nuS/JX2b9X0muS/KDSX5ww7Z8z2wd7phif6/YPpt8PjGXbPv2NJ8sntF8sth2ujvrhfjkMfozU2daMPcjfgY7PZL8naxfDveJDdv+JVPnOo/1OZC+ReZvJzk+2xe/kdmnQk71qFkoAACANlbtZkIAAIC5FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2vn/bUUhBCY5VroAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot zscores for loneliness\n",
    "f, [ax1, ax2, ax3] = plt.subplots(1,3, figsize = (14,6), sharey = True)\n",
    "\n",
    "# Note that there appears to be increasing variation in lonelines prescribing over time.\n",
    "# These means are comparable as standardised using 2016 mean and std.\n",
    "# data[data['Year'] == 2016]['loneliness_zscore'].hist(bins=100, ax = ax1)\n",
    "# data[data['Year'] == 2017]['loneliness_zscore'].hist(bins=100, ax = ax2)\n",
    "data[data['Year'] == 2018]['loneliness_zscore'].hist(bins=100, ax = ax3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save aggregated data\n",
    "data.to_csv(path + \"processed_data_with_postcodes_GPs_stats.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create Loneliness Variable\n",
    "\n",
    "The actual loneliness variable we work with is the sum of the standardised scores of: depression, alzheimers, hypertension, insomnia, addiction and social anxiety, for each year of interest.\n",
    "\n",
    "This means that the loneliness variable is actually an equally weighted index of the above domains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sum function ignores NAs\n",
    "data['loneills'] =  data[['depression_zscore', 'alzheimers_zscore', 'hypertension_zscore', 'insomnia_zscore',\n",
    "                          'addiction_zscore','social anxiety_zscore']].sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzIAAAFpCAYAAABH4x0wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGkBJREFUeJzt3X+MpPddH/D3hxwpaS51EkwW16RZUENEFIvAnaIgpHJH+OHGEgkVVESCGpH2gBaEhFtlC1JxS1GPihCpEuovJcVFhYNSaCxfWjDpXVMQCdyByTm4aUI4gR3LaUTicjQFDN/+sXPx+rzrmdmdndnP+PWSRjs/np15f59n93v7vmfmeWqMEQAAgE4+Y9UBAAAA5qXIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0cW+aL3XzzzWNzc3Pqcn/0R3+U5z//+Ycf6IDkXKwOOTtkTObLefny5Y+PMT7nkCMtnPlkNTrk7JAxWb+c6z6XHFSX7T2PdRxTsp7j6jammeeTMcbSLidOnBizuHDhwkzLrZqci9UhZ4eMY8yXM8mlscR5YFEX88lqdMjZIeMY65dz3eeSg+qyveexjmMaYz3H1W1Ms84n3loGAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAMxgc+t8NrfOrzoGMKHIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0M7UIlNVn1VVv1ZVv1VVH6iqfzy5//Or6n1V9aGq+umqeu7hxwUAAJhtj8wfJ/nKMcYXJ3l1ktur6rVJfjjJ28YYL0/yiSRvPryYAAAAT5paZMa2a5Obnzm5jCRfmeRnJ/ffk+SNh5IQAADgBjN9RqaqnlNVDyT5WJL7k/xOkk+OMZ6YLPJwklsPJyIAAMBT1Rhj9oWrXpjk55P8oyT/bozxVyf3vzTJu8YYt+3yPWeSnEmSjY2NE+fOnZv6OteuXcvx48dnzrUqci5Wh5wdMibz5Tx9+vTlMcbJQ460EOaT1euQs0PGZP1yrvtckiRXHnk8SXLbrTfN/Zpdtvc81nFMyXqOq9uYZp5PxhhzXZL8QJJ/kOTjSY5N7vuyJL8w7XtPnDgxZnHhwoWZlls1ORerQ84OGceYL2eSS2POeeAoXMwnq9EhZ4eMY6xfznWfS8YY42VvuW+87C33zbz8Tl229zzWcUxjrOe4uo1p1vlklqOWfc5kT0yq6nlJvirJQ0kuJPmGyWJ3JnnnfF0LAABgf47NsMwtSe6pqudk+zM1PzPGuK+qfjvJuar6p0l+M8nbDzEnAADAp00tMmOM9yf5kl3u/0iS1xxGKAAAgGcy01HLAAAAjhJFBgAAaEeRAQAA2lFkAACAdhQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQBgoTa3zmdz6/yqY7DmFBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2lFkAACAdhQZAACgHUUGAIBDtbl1Pptb51cdgzWjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0c2zVAQAA6On6SS6vnr3jKbdhGeyRAQAA2lFkAACAdhQZAACgHUUGAABoR5EBAADamVpkquqlVXWhqh6qqg9U1fdM7r+7qh6pqgcml9cfflwAAIDZDr/8RJK7xhi/UVUvSHK5qu6fPPa2McaPHF48AACAp5taZMYYjyZ5dHL9D6vqoSS3HnYwAACAvdQYY/aFqzaTvCfJq5J8b5JvTfJ/klzK9l6bT+zyPWeSnEmSjY2NE+fOnZv6OteuXcvx48dnzrUqci5Wh5wdMibz5Tx9+vTlMcbJQ460EOaT1euQs0PGZP1yrvtckiRXHnk8SXLbrTfN/Zpdtvc8rl27lt99/M+SPLlOrq+j6268fz/rbtnWdVt1GtPM88kYY6ZLkuNJLif5G5PbG0mek+3P2fxQkndMe44TJ06MWVy4cGGm5VZNzsXqkLNDxjHmy5nk0phxHjhKF/PJanTI2SHjGOuXc93nkjHGeNlb7hsve8t9My+/U5ftPY8LFy48bZ1cv73X/R2s67bqZNb5ZKajllXVZyb5T0n+wxjj5yYF6LExxp+NMf48yb9N8pr5uhYAAMD+zHLUskry9iQPjTF+dMf9t+xY7OuTPLj4eAAAAE83y1HLvjzJtyS5UlUPTO77viRvqqpXJxlJrib59kNJCAAAcINZjlr2y0lql4fetfg4AAAA0830GRkAAICjRJEBAADaUWQAAIB2FBkAAKCdWY5aBgDAs9Tm1vlPX7969o4VJoGnskcGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgCAmWxunX/KCTJhlRQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgCAuWxunc+VRx5fdQye5RQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHaOrToAAADraXPr/KojsMbskQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2plaZKrqpVV1oaoeqqoPVNX3TO5/cVXdX1Ufmnx90eHHBQAAmG2PzBNJ7hpjfFGS1yb5e1X1yiRbSd49xnh5kndPbgMAABy6qUVmjPHoGOM3Jtf/MMlDSW5N8oYk90wWuyfJGw8rJAAAwE5zfUamqjaTfEmS9yXZGGM8mmyXnSQvWXQ4AACA3dQYY7YFq44n+e9JfmiM8XNV9ckxxgt3PP6JMcbTPidTVWeSnEmSjY2NE+fOnZv6WteuXcvx48dnHMLqyLlYHXJ2yJjMl/P06dOXxxgnDznSQphPVq9Dzg4Zk/XLue5zSZJceeTxJMltt94092t22d67uT7uG208L3nsU/M9137W3bJ13lZ76TamWeeTmYpMVX1mkvuS/MIY40cn930wyakxxqNVdUuSi2OMVzzT85w8eXJcunRp6utdvHgxp06dmrrcqsm5WB1ydsiYzJezqtr88bGT+WQ1OuTskDFZv5zrPpckyebW+STJ1bN3zP06Xbb3bq6P+0Z33fZE3nrl2FzPtZ91t2ydt9Veuo1p1vlklqOWVZK3J3noeomZuDfJnZPrdyZ5536CAgAAzGuWGv3lSb4lyZWqemBy3/clOZvkZ6rqzUl+L8k3Hk5EAACAp5paZMYYv5yk9nj4dYuNAwAAMN1cRy0DAAA4ChQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2lFkAACAdhQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQDg0za3zmdz6/yqY8BUigwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO8dWHQAAgGen3U68efXsHStIQkf2yAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7UwtMlX1jqr6WFU9uOO+u6vqkap6YHJ5/eHGBAAAeNIse2R+PMntu9z/tjHGqyeXdy02FgAAwN6mFpkxxnuS/MESsgAAAMzkIJ+R+a6qev/krWcvWlgiAACAKWqMMX2hqs0k940xXjW5vZHk40lGkh9McssY49v2+N4zSc4kycbGxolz585Nfb1r167l+PHjs41gheRcrA45O2RM5st5+vTpy2OMk4ccaSHMJ6vXIWeHjMn65Vz3uSRJrjzyeJLktltvmvs1u2zv5MlxTrPxvOSxT+3vNa6vw91eaz/rd5E6batZdRvTrPPJvorMrI/d6OTJk+PSpUtTX+/ixYs5derU1OVWTc7F6pCzQ8ZkvpxV1eaPj53MJ6vRIWeHjMn65Vz3uSRJNrfOJ0munr1j7tfpsr2TJ8c5zV23PZG3Xjm2r9e4vg53e639rN9F6rStZtVtTLPOJ/t6a1lV3bLj5tcneXCvZQEAABZtao2uqp9KcirJzVX1cJIfSHKqql6d7beWXU3y7YeYEQAA4CmmFpkxxpt2ufvth5AFAABgJgc5ahkAAMBKKDIAAEA7igwAANCOIgMAALSzv4N/AwDAPs16rhp4JvbIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7x1YdAACA1dvcOr/qCLu6nuvq2TtWnISjxh4ZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHacEBMA4FnsqJ8IE/ZijwwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0M7UIlNV76iqj1XVgzvue3FV3V9VH5p8fdHhxgQAAHjSLHtkfjzJ7Tfct5Xk3WOMlyd59+Q2AADAUkwtMmOM9yT5gxvufkOSeybX70nyxgXnAgAA2NN+PyOzMcZ4NEkmX1+yuEgAAADPrMYY0xeq2kxy3xjjVZPbnxxjvHDH458YY+z6OZmqOpPkTJJsbGycOHfu3NTXu3btWo4fPz5L/pWSc7E65OyQMZkv5+nTpy+PMU4ecqSFMJ+sXoecHTIm65dz3eeSJLnyyONJkttuvWnu1zzK2/v6uOa18bzksU8tOEyeXL835trPet+Po7yt9qvbmGadT/ZbZD6Y5NQY49GquiXJxTHGK6Y9z8mTJ8elS5emvt7Fixdz6tSpqcutmpyL1SFnh4zJfDmrqs0fHzuZT1ajQ84OGZP1y7nuc0mSbG6dT5JcPXvH3K9zlLf39XHN667bnshbrxxbcJon1++Nufaz3vfjKG+r/eo2plnnk/2+tezeJHdOrt+Z5J37fB4AAIC5zXL45Z9K8qtJXlFVD1fVm5OcTfLVVfWhJF89uQ0AALAUU/cHjjHetMdDr1twFgAAgJns961lAAAAK6PIAAAA7SgyAABAO4oMAADQzuIP/g0AwJG33/PHwFFhjwwAANCOIgMAALSjyAAAAO0oMgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAO06ICQDAkTHriTqvL3f17B2HGYcjzB4ZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHYUGQAAoB1FBgAAaEeRAQAA2lFkAACAdhQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGjn2KoDAABw+Da3zidJrp69Y8VJDub6OPa6v/v4mJ09MgAAQDuKDAAA0I4iAwAAtKPIAAAA7SgyAABAOwc6allVXU3yh0n+LMkTY4yTiwgFAADwTBZx+OXTY4yPL+B5AAAAZuKtZQAAQDsHLTIjyS9W1eWqOrOIQAAAANPUGGP/31z1l8cYH62qlyS5P8l3jzHec8MyZ5KcSZKNjY0T586dm/q8165dy/Hjx/eda1nkXKwOOTtkTObLefr06ctdPt9mPlm9Djk7ZEzWL+e6zyVJcuWRx5Mkt91609yvOW09HuS5Z3Xja1y/vV8bz0se+9SBYx3YjeM56Drs8rs5j25jmnU+OVCRecoTVd2d5NoY40f2WubkyZPj0qVLU5/r4sWLOXXq1EJyHSY5F6tDzg4Zk/lyVlWbPz52Mp+sRoecHTIm65dz3eeSJNncOp8kuXr2jrlfZ9p6PMhzz+rG17h+e7/uuu2JvPXKIj5ufTA3jueg67DL7+Y8uo1p1vlk328tq6rnV9ULrl9P8jVJHtzv8wEAAMzqIDV6I8nPV9X15/nJMcZ/XUgqAACAZ7DvIjPG+EiSL15gFgAAgJk4/DIAANCOIgMAALSjyAAAAO0oMgAAQDuKDADAIdjcOn/gc7UAe1NkAACAdhQZAACgHUUGAABoR5EBAADaUWQAAIB2FBkAAKAdRQYAAGhHkQEAANpRZAAAgHaOrToAAEBnm1vnkyRXz96x4iTPTtfX/163E9tmXdkjAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0o8gAAADtKDIAAEA7igwAANCOE2ICAMxhtxMuTlv+rtueyLdunf/0iRn3eo5ZT645y0kf583Jk5zktAd7ZAAAgHYUGQAAoB1FBgAAaEeRAQAA2lFkAACAdhQZAACgHUUGAABoR5EBAADacUJMAIAFmOUElId5kspZn/vZeKLMG8e814lJZz0R6Twnypz3JKfPtNyiTtS5c9x7nUh1kScDPawTjNojAwAAtKPIAAAA7SgyAABAO4oMAADQjiIDAAC0c6AiU1W3V9UHq+rDVbW1qFAAAADPZN9Fpqqek+THkvz1JK9M8qaqeuWiggEAAOzlIHtkXpPkw2OMj4wx/iTJuSRvWEwsAACAvR2kyNya5Pd33H54ch8AAMChqjHG/r6x6huTfO0Y429Pbn9LkteMMb77huXOJDkzufmKJB+c4elvTvLxfQVbLjkXq0PODhmT+XK+bIzxOYcZZlHMJ0dCh5wdMibrl3Pd55KD6rK957GOY0rWc1zdxjTTfHKQIvNlSe4eY3zt5PY/TJIxxj/b1xM+9bkvjTFOHvR5Dpuci9UhZ4eMSZ+cy9BlXci5OB0yJnI+26zjelzHMSXrOa51HFNysLeW/XqSl1fV51fVc5N8U5J7FxMLAABgb8f2+41jjCeq6ruS/EKS5yR5xxjjAwtLBgAAsId9F5kkGWO8K8m7FpRlp39zCM95GORcrA45O2RM+uRchi7rQs7F6ZAxkfPZZh3X4zqOKVnPca3jmPb/GRkAAIBVOchnZAAAAFZipUWmqm6vqg9W1YeramuXx/9CVf305PH3VdXm8lPOlPN7q+q3q+r9VfXuqnrZUcy5Y7lvqKpRVUs/esUsGavqb07W5weq6ieXnXGSYdo2/ytVdaGqfnOy3V+/gozvqKqPVdWDezxeVfUvJmN4f1V96bIzLpP5ZHkZdyy3srlk8vrmk8VlNJ8sQVXdXVWPVNUDk8vSt/WizDpPdFJVV6vqymTbXFp1nv3a7fe5ql5cVfdX1YcmX1+0yowLM8ZYySXbBwj4nSRfkOS5SX4ryStvWObvJvlXk+vflOSnj2jO00n+4uT6dx7VnJPlXpDkPUnem+TkUcuY5OVJfjPJiya3X3IU12W232v6nZPrr0xydQU5/1qSL03y4B6Pvz7Jf0lSSV6b5H3LznjEtpn5ZEEZJ8utbC6ZY12aT2bPaT5Zznq+O8nfX3WOBYxjpnmi2yXJ1SQ3rzrHAsbxtN/nJP88ydbk+laSH151zkVcVrlH5jVJPjzG+MgY40+SnEvyhhuWeUOSeybXfzbJ66qqlpgxmSHnGOPCGOP/Tm6+N8nnLTljMtv6TJIfzPYP8/9bZriJWTL+nSQ/Nsb4RJKMMT625IzJbDlHkr80uX5Tko8uMd92gDHek+QPnmGRNyT592Pbe5O8sKpuWU66pTOfLDHjxCrnksR8slDmE+Y06zzBCuzx+7zz38B7krxxqaEOySqLzK1Jfn/H7Ycn9+26zBjjiSSPJ/nspaTbJcPEbjl3enO2/9dq2abmrKovSfLSMcZ9ywy2wyzr8guTfGFV/UpVvbeqbl9auifNkvPuJN9cVQ9n+8h9372caHOZ92e3M/PJ4nSYSxLzybI9m+aTw/Zdk7fnvaPx23vW9edhJPnFqrpcVWdWHWbBNsYYjybJ5OtLVpxnIQ50+OUD2u1/Qm88hNosyxy2mTNU1TcnOZnkKw410e6eMWdVfUaStyX51mUF2sUs6/JYtt8Ocirb/xP9P6rqVWOMTx5ytp1myfmmJD8+xnhrVX1Zkp+Y5Pzzw483s6Pw+7Ms5pPF6TCXJOaTZTsKvz8tVNUvJfncXR76/iT/Mtt7M8fk61uTfNvy0i3Muv48fPkY46NV9ZIk91fV/5zs3eCIWmWReTjJS3fc/rw8fXf69WUerqpj2d7l/ky7vg/DLDlTVV+V7UnqK8YYf7ykbDtNy/mCJK9KcnHybprPTXJvVX3dGGNZH2ibdZu/d4zxp0l+t6o+mO0/RH59ORE/nWFazjcnuT1Jxhi/WlWfleTmJKt468peZvrZXRPmk8XpMJck5pNlezbNJwcyxviqWZarqn+bZJV7NQ9iLX8exhgfnXz9WFX9fLbfQrcuReaxqrpljPHo5G2hR2l+2bdVvrXs15O8vKo+v6qem+0P3957wzL3Jrlzcv0bkvy3MfmU0hJNzTl5m8W/TvJ1K3oPdjIl5xjj8THGzWOMzTHGZrbfe7/sPzxm2eb/Odsfdk5V3Zztt4Z8ZIkZk9ly/l6S1yVJVX1Rks9K8r+XmnK6e5P8rcnRhl6b5PHru5XXkPlkSRmPyFwyNeeE+WRxnk3zyaG54XNFX59k16PENTDLz3UrVfX8qnrB9etJviZ9t89udv4beGeSd64wy+Ks8kgD2T4Kyv/K9pEvvn9y3z/J9j+KyfZk/h+TfDjJryX5giOa85eSPJbkgcnl3qOY84ZlL2Y1Rxqati4ryY8m+e0kV5J801Fcl9k+stCvZPtILQ8k+ZoVZPypJI8m+dNs/+/Ym5N8R5Lv2LEuf2wyhiur2N5HbJuZTxaU8YZlVzKXzLguzSezZzSfLGc9/8Rk/b0/239Y3rLqTAcYy9N+rjtfsn0Ett+aXD7QeUx7/D5/dpJ3J/nQ5OuLV51zEZeaDBgAAKCNlZ4QEwAAYD8UGQAAoB1FBgAAaEeRAQAA2lFkAACAdhQZAACgHUUGAABoR5EBAADa+f86uPWPnpWjjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot zscores for loneills\n",
    "f, [ax1, ax2, ax3] = plt.subplots(1,3, figsize = (14,6), sharey = True)\n",
    "\n",
    "# Note that there appears to be increasing variation in lonelines prescribing over time.\n",
    "data[data['Year'] == 2016]['loneills'].hist(bins=100, ax = ax1)\n",
    "data[data['Year'] == 2017]['loneills'].hist(bins=100, ax = ax2)\n",
    "data[data['Year'] == 2018]['loneills'].hist(bins=100, ax = ax3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save aggregated data\n",
    "data.to_csv(path + \"final_data.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}