Skip to content

Instantly share code, notes, and snippets.

@drcjar

drcjar/ipfjes-analysis.ipynb

Created July 4, 2018 12:50
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save drcjar/cc6e106c92c17cc3a67085b3bdcd9fe5 to your computer and use it in GitHub Desktop.
Save drcjar/cc6e106c92c17cc3a67085b3bdcd9fe5 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
ipfjes-analysis.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import psycopg2 as pg\n",
"import pandas.io.sql as psql\n",
"import datetime as dt\n",
"import numpy as np\n",
"import scipy.stats as stats\n",
"from ukpostcodeutils import validation\n",
"import geopy.distance\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sbn\n",
"import collections\n",
"#from collections import OrderedDict\n",
"%matplotlib inline\n",
"pd.options.display.float_format = '{:,.2f}'.format"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def dict_from_db():\n",
" \"\"\"\n",
" Returns a dict of dataframes from the database tables\n",
" which are populated by the study application\n",
" \"\"\"\n",
" conn = pg.connect(\"dbname='carlplaying' user='drcjar' password='drcjar_is_not_fake'\")\n",
"\n",
" cursor = conn.cursor()\n",
" cursor.execute(\"select relname from pg_class where relkind='r' and relname !~ '^(pg_|sql_)';\")\n",
" e = cursor.fetchall()\n",
"\n",
" table_name_list = [i[0] for i in e]\n",
"\n",
" select_template = 'SELECT * FROM {table_name}'\n",
" frames_dict = {}\n",
" for tname in table_name_list:\n",
" query = select_template.format(table_name = tname)\n",
" frames_dict[tname] = pd.read_sql(query, conn)\n",
" \n",
" return frames_dict\n",
"\n",
"def notepisodes():\n",
" \"\"\"\n",
" Returns a list of episode_ids that are not in fact episodes for analysis purposes\n",
" e.g patients who have not yet been interviewed\n",
" This list is made in the quality check notebook\n",
" \"\"\"\n",
" notepisodes = []\n",
"\n",
" with open(\"notepisodes.txt\", \"r\") as f:\n",
" for line in f:\n",
" notepisodes.append(int(line.strip()))\n",
" \n",
" return notepisodes\n",
"\n",
"def study_df(df):\n",
" \"\"\"\n",
" Makes the main study dataframe from our occupational dataframe\n",
" \"\"\"\n",
" df = df[['patient_id', 'soc_job_ft', 'soc_job_fk_id', 'job_tasks', 'start_year', 'end_year']].copy() # e.g for debugging\n",
" \n",
" df = df[~df.patient_id.isin(notepisodes())] # discard not episodes\n",
" \n",
" # add participant type\n",
" pt_lookup = ddf['ipfjes_studyparticipantdetails'][['episode_id', 'participant_type']].copy()\n",
" pt_lookup.index = pt_lookup['episode_id']\n",
" pt_lookup = pt_lookup['participant_type'].to_dict()\n",
" df['pt'] = df['patient_id'].map(lambda x: pt_lookup.get(x)).copy()\n",
" \n",
" # add job titles\n",
" df = df.replace('', np.nan) # stick a NaN in missing soc_job_ft and replace with a job title from our lookup\n",
" job_lookup = ddf['ipfjes_socjob'][['id','name']]\n",
" job_lookup.index = job_lookup['id']\n",
" job_lookup = job_lookup['name'].to_dict()\n",
" df['jt'] = df['soc_job_fk_id'].map(lambda x: job_lookup.get(x))\n",
" df.soc_job_ft.fillna(df.jt, inplace=True)\n",
" del df['jt']\n",
" del df['soc_job_fk_id']\n",
" \n",
" # add soc 90 code\n",
" soc90_lookup = ddf['ipfjes_soccode'][['soc90', 'title']]\n",
" soc90_lookup.index = soc90_lookup['title']\n",
" soc90_lookup = soc90_lookup['soc90'].to_dict()\n",
" df['soc90'] = df.soc_job_ft.map(lambda x: soc90_lookup.get(x))\n",
" \n",
" # add soc 2000 code\n",
" soc2000_lookup = ddf['ipfjes_soccode'][['soc2000', 'title']]\n",
" soc2000_lookup.index = soc2000_lookup['title']\n",
" soc2000_lookup = soc2000_lookup['soc2000'].to_dict()\n",
" df['soc2000'] = df.soc_job_ft.map(lambda x: soc2000_lookup.get(x))\n",
"\n",
" # add socioeconomic code\n",
" ssec_lookup = pd.read_csv('nssec/nssecderivationtablessoc2000_tcm77-264822.csv') \n",
" # nb we're using sheet fm2kac even though we don't know what it means. we should know.\n",
" ssec_lookup = ssec_lookup[['SOC2000', 'ssec']] \n",
" # simplified socioeconomic classification using soc2000\n",
" ssec_lookup.index = ssec_lookup['SOC2000']\n",
" ssec_lookup = ssec_lookup['ssec'].to_dict()\n",
" df['ssec'] = df.soc2000.map(lambda x: ssec_lookup.get(x))\n",
"\n",
" # add dob\n",
" dob_lookup = ddf['ipfjes_demographics'][['id', 'date_of_birth']]\n",
" dob_lookup.index = dob_lookup['id']\n",
" dob_lookup = dob_lookup['date_of_birth'].to_dict()\n",
" df['dob'] = df.patient_id.map(lambda x: dob_lookup.get(x))\n",
" \n",
" # add age in years and agegroup\n",
" df['years'] = ((dt.datetime.now() - pd.to_datetime(df.dob)).dt.days / 365).astype(int)\n",
" df['agegroup'] = pd.cut(df.years, range(25,100,5), \n",
" labels=[\n",
" '25 to 29',\n",
" '30 to 34',\n",
" '35 to 39',\n",
" '40 to 44',\n",
" '45 to 49',\n",
" '50 to 54',\n",
" '55 to 59',\n",
" '60 to 64',\n",
" '65 to 69',\n",
" '70 to 74',\n",
" '75 to 79',\n",
" '80 to 84',\n",
" '85 to 90',\n",
" '90 to 94'],\n",
" right=False) \n",
" \n",
" # add ethnicity\n",
" eth_lookup = ddf['ipfjes_demographics'][['id', 'ethnicity_ft']]\n",
" eth_lookup.index = eth_lookup['id']\n",
" eth_lookup = eth_lookup['ethnicity_ft'].to_dict()\n",
" df['ethnicity'] = df.patient_id.map(lambda x: eth_lookup.get(x))\n",
" \n",
" # add smoking\n",
" es_lookup = ddf['ipfjes_smokinghistory'][['episode_id', 'ever_smoked']]\n",
" es_lookup.index = es_lookup['episode_id']\n",
" es_lookup = es_lookup['ever_smoked'].to_dict()\n",
" df['es'] = df.patient_id.map(lambda x: es_lookup.get(x))\n",
" \n",
" cs_lookup = ddf['ipfjes_smokinghistory'][['episode_id', 'current_smoker']]\n",
" cs_lookup.index = cs_lookup['episode_id']\n",
" cs_lookup = cs_lookup['current_smoker'].to_dict()\n",
" df['cs'] = df.patient_id.map(lambda x: cs_lookup.get(x))\n",
" \n",
" # add participant_id\n",
" ct_lookup = ddf['ipfjes_demographics'][['patient_id', 'hospital_number']].copy()\n",
" ct_lookup.index = ct_lookup['patient_id']\n",
" ct_lookup = ct_lookup['hospital_number'].to_dict()\n",
" df['participant_id'] = df['patient_id'].map(lambda x: ct_lookup.get(x)).copy()\n",
" \n",
" # add centre_id\n",
" # note that we encode centre id in first two digits of hospital number\n",
" # note patient_id == episode_id in the demographics table\n",
" df['centre'] = df['participant_id'].str[:2].astype(int) \n",
" \n",
" # add gp long/lat\n",
" centroids = pd.read_csv('ONS_Postcode_Directory_Latest_Centroids.csv', usecols=[0,1,3,4,5,12,13,44])\n",
" centroids['pcd'] = centroids['pcd'].str.strip()\n",
" centroids['pcd'] = centroids['pcd'].str.replace(\" \",\"\")\n",
" centroids['pcd'] = centroids['pcd'].str.replace(\" \",\"\")\n",
" \n",
" gp_lookup = ddf['ipfjes_demographics'][ddf['ipfjes_demographics']['contact_details'].notna()] \n",
" # quality issues with gp postcode addressed in quality notebook\n",
" gp_lookup = gp_lookup[['patient_id', 'contact_details']]\n",
" gp_lookup.index = gp_lookup.patient_id\n",
" gp_lookup['contact_details'] = gp_lookup['contact_details'].str[-8:]\n",
" gp_lookup['contact_details'] = gp_lookup['contact_details'].str.strip()\n",
" gp_lookup['contact_details'] = gp_lookup['contact_details'].str.replace(\" \",\"\")\n",
" gp_lookup = gp_lookup[gp_lookup['contact_details'].map(lambda x: validation.is_valid_postcode(str(x)))]\n",
" gp_lookup.columns = ['patient_id', 'pcd']\n",
" gp_lookup_centroids = pd.merge(gp_lookup, centroids, on='pcd')\n",
" gp_lookup_centroids['gp_coords'] = list(zip(gp_lookup_centroids.Y, gp_lookup_centroids.X))\n",
" \n",
" gp_lookup_centroids = gp_lookup_centroids[['patient_id', 'gp_coords']].copy()\n",
" gp_lookup_centroids.index = gp_lookup_centroids['patient_id']\n",
" gp_lookup_centroids = gp_lookup_centroids['gp_coords'].to_dict()\n",
" df['gp_coords'] = df['patient_id'].map(lambda x: gp_lookup_centroids.get(x)).copy()\n",
" \n",
" trustpcds = pd.read_csv('trustpcd.csv')\n",
" trustpcds['pcd'] = trustpcds['pcd'].str.strip()\n",
" trustpcds['pcd'] = trustpcds['pcd'].str.replace(\" \",\"\")\n",
" trustpcds_centroids = pd.merge(trustpcds, centroids, on='pcd')\n",
" trustpcds_centroids['centre_coords'] = list(zip(trustpcds_centroids.Y, trustpcds_centroids.X))\n",
" \n",
" tp_lookup_centroids = trustpcds_centroids[[\"Centre ID\", \"centre_coords\"]].copy()\n",
" tp_lookup_centroids.index = tp_lookup_centroids['Centre ID']\n",
" tp_lookup_centroids = tp_lookup_centroids['centre_coords'].to_dict()\n",
" df['centre_coords'] = df['centre'].map(lambda x: tp_lookup_centroids.get(x)).copy()\n",
" \n",
" df['distfromcentre'] = df.apply(dist_from_centre, axis=1)\n",
" # pack years to be added here\n",
" \n",
" df = cumrisk(df) # add cum risk\n",
" return df\n",
" \n",
"def asbexposed():\n",
" \"\"\"\n",
" Returns a df of at asbestos exposed SOC2000 codes from David Coggon\n",
" Definition of asbestos-exposed job groups: \n",
" those groups with significantly elevated PMR for \n",
" pleural cancer during the whole period of 1979-80, 1982-2010\n",
" \"\"\"\n",
" df = pd.read_csv('coggon_soc2000_meso_pmr.csv')\n",
" df.soc2000 = df.soc2000.astype(str)\n",
" return df\n",
" \n",
"def asbexposed2():\n",
" \"\"\"\n",
" Returns a df of at asbestos exposed SOC2000 codes from Julian Peto\n",
" Definition of asbestos-exposed job groups: \n",
" those groups with significantly elevated PMR for \n",
" pleural cancer during the whole period of 1979-80, 1982-2010\n",
" \"\"\"\n",
" df = pd.read_excel('./data/occ80007.xls', skiprows=2) \n",
" df = df[(df.Lower > 100) & (df.Deaths > 30) & (df.PMR > 200)]\n",
" df['soc1990'] = df['SOC Occupation Code'].astype(str)\n",
" return df\n",
"\n",
"def clean_list(messylist):\n",
"#cleans a list containing ranges by expanding ranges e.g 1-4 becomes 1,2,3,4. \n",
"#like the ones we have in table 2.3.2 www.hse.gov.uk/research/rrpdf/rr696.pdf\n",
"\n",
" nums = []\n",
" expandnums = []\n",
"\n",
" messylist = messylist.split(\",\")\n",
" messylist = [s.strip() for s in messylist]\n",
"\n",
" for s in messylist:\n",
" if '-' in s:\n",
" s = s.split(\"-\")\n",
" s = [int(i) for i in s]\n",
" expandnums.append(range(s[0], s[1]+1))\n",
" else:\n",
" nums.append(int(s))\n",
" \n",
" expandnums = [item for sublist in expandnums for item in sublist]\n",
" cleanlist = nums + expandnums\n",
" cleanlist.sort()\n",
" \n",
" return cleanlist\n",
"\n",
"def makejobcatlookup():\n",
" #makes a jobgrouplookup that maps soc codes to asbestos exposure group\n",
" #from table 2.3.2 www.hse.gov.uk/research/rrpdf/rr696.pdf\n",
" #and ***edited to remove duplicates so it doesn't break the reverse lookup***\n",
" #uses soc 1990\n",
" #need to add nec to office and low risk\n",
" \n",
" #manually added by me 5:'551'\n",
" #?579 woodworking n.e.c\n",
" #?814\n",
" \n",
" # (2, clean_list(\"\"\"570, 920,\n",
" # 532, 913,\n",
" # 521,\n",
" # 507,\n",
" # 111, 500-506, 509, 885, 886, 895, 921, 923, 924, 913,\n",
" # 530, 579\"\"\")),\n",
" \n",
" jobclassification = ((5, clean_list(\"\"\"\n",
" 101, 120-127, 130-132, 139, 154, 155, 160, 169, 170, 172-\n",
" 179, 190, 191, 220-224, 230-235, 239-242, 250-253, 260,\n",
" 261, 270, 271, 290-293, 320, 340-347, 361-363, 370, 371,\n",
" 380, 381, 383-387, 390, 392, 399-401, 410-412, 420, 421,\n",
" 430, 450-452, 459-462, 490, 551, 556, 559, 560, 569, 580-582,\n",
" 592, 594, 595, 598, 610, 619-622, 630, 640-644, 650-652,\n",
" 659-661, 670, 671, 673, 690, 691, 699, 700-703, 710, 719-\n",
" 722, 730, 732, 790-792, 900-904, 950-956, 958,\n",
" 134, 394, 140, 152, 800, 550, 551, 555, 557, 957, 199,\n",
" 364, 360, 100, 209, 331, 350, 393,\n",
" 491, 201, 593, 932, 102, 614, 313, 202, 613, \n",
" 530, 303, 563, 312, 562, 103, 591, 802\"\"\")),# I added last 2 rows (beginning 491 and 530)\n",
" \n",
" (4, clean_list(\"\"\"540,\n",
" 310,\n",
" 210-219,\n",
" 440-441,\n",
" 150-151, 600-601,\n",
" 731, 870-875,\n",
" 113, 153, 171, 304, 348, 396, 531, 542-544, 553, 561, 569,\n",
" 571, 590, 596, 597, 599, 611, 612, 615, 619, 631, 642, 672,\n",
" 699, 733, 801, 809, 811, 820, 822, 824, 825, 829, 910, 911-\n",
" 913, 919, 923, 924, 931, 933, 934, 940, 941, 955, 956, 958,\n",
" 990, 999,\n",
" 812, 395, 141, 142, 572, 112, 823, 552, 821, 814\"\"\")),\n",
" (3, clean_list(\"\"\"516, 913,\n",
" 881-884, 922,\n",
" 200, 300-302, 309,\n",
" 110, 260, 262, 311,\n",
" 510-515, 517-519,\n",
" 520, 522-529,\n",
" 535-537,\n",
" 830-844,\n",
" 850-869,\n",
" 887-892, 894, 897-899,\n",
" 391, 394, 554, 530\"\"\")), # I added 530 - blacksmith\n",
" (2.3, clean_list(\"\"\"111, 500-506, 509, 885, 886, 895, 923, 924\"\"\")),\n",
" (2.2, clean_list(\"\"\"532, 521, 507, 913\"\"\")), \n",
" (2.1, clean_list(\"\"\"570, 920\"\"\")),\n",
" (1, clean_list(\"\"\"533, 534,\n",
" 541,\n",
" 893, 896, 921, 929, 990,\n",
" 880, 332, 903, 169, 173, 174, 239, 385, 463, 620, 621,\n",
" 630, 900, 930, 952, 953\"\"\")))\n",
"\n",
"\n",
" jobclassification = collections.OrderedDict(jobclassification) \n",
"\n",
" #invert the dictionary to make our look up\n",
" jobcatlookup = {value: key for key in jobclassification for value in jobclassification[key]}\n",
" return jobcatlookup\n",
"\n",
"def oddsratio(df):\n",
" \"\"\"\n",
" Takes a df and adds a column 'exposed' based on whether or not\n",
" soc2000 code is in list of Coggon exposed socs.\n",
" Then makes a summary dataframe where a patient having any job\n",
" that is exposed counts as being exposed and uses this \n",
" dataframe to calculate the odds ratio.\n",
" \"\"\"\n",
" df['exposed'] = df.soc2000.isin(ae.soc2000)\n",
" df.index = df.patient_id\n",
" exposedany = df.groupby(level='patient_id').exposed.any().to_frame()\n",
" casestatus = df.drop_duplicates(subset='patient_id').pt.to_frame()\n",
" summarydf = exposedany.merge(casestatus, right_index=True, left_index=True)\n",
" table = summarydf.groupby('pt').exposed.value_counts().values\n",
" table = ([[table[1], table[3]], [table[0], table[2]]])\n",
" oddsratio, pvalue = stats.fisher_exact(table)\n",
" print(\"OddsR: \", oddsratio, \"p-Value:\", pvalue)\n",
" \n",
"def highest_ssec(df):\n",
" \"\"\"\n",
" Takes a df makes a summary dataframe of highest\n",
" socioeconomic, read lowest, code for each participant. \n",
" \"\"\"\n",
" df.index = df.patient_id\n",
" highssec = df.groupby(level='patient_id').ssec.min().to_frame()\n",
" casestatus = df.drop_duplicates(subset='patient_id').pt.to_frame()\n",
" summarydf = highssec.merge(casestatus, right_index=True, left_index=True)\n",
" print('\\n# Social class (based on lowest social class code job)\\n')\n",
" print (summarydf.groupby('pt').ssec.describe())\n",
" \n",
"def cumrisk(df):\n",
" \"\"\"\n",
" Takes a df and calculates av cummulative risk \n",
" exposed y/n * years for cases and controls\n",
" using Coggon exposure codes.\n",
" \"\"\"\n",
" df.end_year = df.end_year.fillna(0)\n",
" df.start_year = df.start_year.fillna(0)\n",
" df['exposed'] = df.soc2000.isin(ae.soc2000)\n",
" df['duration'] = df.end_year.astype(int) - df.start_year.astype(int) \n",
" df = df[df['duration'] > 0] # temp hack\n",
" df['risk'] = df['duration'] * df['exposed']\n",
" print('\\n# Average cummulative risk\\n')\n",
" print(df.groupby('pt').risk.sum() / df.groupby('pt').patient_id.nunique())\n",
" return df\n",
"\n",
"def dist_from_centre(x): \n",
" \"\"\"\n",
" Calculate distance from patients GP practice to their\n",
" recruitment centre.\n",
" \"\"\"\n",
" return geopy.distance.vincenty(x.gp_coords, x.centre_coords).km\n",
"\n",
"def ever_vs_never(df):\n",
" \"\"\"\n",
" OR for ever vs never exposed\n",
" \"\"\"\n",
" he = df.groupby(['pt','participant_id']).jobcat.min().reset_index() # highest exposed\n",
" he = he.groupby(['pt','jobcat']).participant_id.count().reset_index()\n",
" he.index = he.jobcat\n",
" \n",
" a = he[he['pt'] == 'case'][['participant_id']]\n",
" b = he[he['pt'] == 'control'][['participant_id']]\n",
" c = pd.concat([a,b], axis=1)\n",
" c.columns = ['cases', 'controls']\n",
" \n",
" ever = c.iloc[0:6][['cases', 'controls']].sum()\n",
" ever.name = 'ever'\n",
" never = c.iloc[6:][['cases', 'controls']].sum()\n",
" never.name = 'never'\n",
" evn = pd.concat([ever, never], axis=1).transpose()\n",
" evn['odds ratio'] = (evn['cases'].astype(float) / evn['controls'].astype(float)) \\\n",
" / (evn.loc['never']['cases'].astype(float) / evn.loc['never']['controls'].astype(float))\n",
" \n",
" oddsratio, pvalue = stats.fisher_exact([evn['cases'].values, evn['controls'].values])\n",
" evn['p-value'] = pvalue\n",
" print(evn)\n",
" return evn\n",
"\n",
"def bycat(df):\n",
" \"\"\"\n",
" OR by cat of exposure\n",
" \"\"\"\n",
" he = df.groupby(['pt','participant_id']).jobcat.min().reset_index() # highest exposed\n",
" he = he.groupby(['pt','jobcat']).participant_id.count().reset_index()\n",
" he.index = he.jobcat\n",
" \n",
" a = he[he['pt'] == 'case'][['participant_id']]\n",
" b = he[he['pt'] == 'control'][['participant_id']]\n",
" c = pd.concat([a,b], axis=1)\n",
" c.columns = ['cases', 'controls']\n",
" c.index = jobgroups.values()\n",
" Construction = c.iloc[1] + c.iloc[2] + c.iloc[3]\n",
" Construction.name = 'Construction'\n",
" c = c.append(Construction)\n",
" c['odds ratio'] = (c['cases'].astype(float) / c['controls'].astype(float)) \\\n",
" / (c.iloc[6]['cases'].astype(float) / c.iloc[6]['controls'].astype(float))\n",
" \n",
" odds = []\n",
" p = []\n",
" for cat in range(len(c)):\n",
" t = c.iloc[[cat,6]]\n",
" oddsratio, pvalue = stats.fisher_exact([t['cases'].values, t['controls'].values])\n",
" odds.append(oddsratio)\n",
" p.append(pvalue)\n",
" \n",
" c['p-value'] = p\n",
" \n",
" print(c.iloc[[0,7,4,5,6]])\n",
" print('\\n')\n",
" print(c.iloc[[0,1,2,3,4,5,6]])\n",
" \n",
" return c \n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/drcjar/anaconda3/envs/ipfjes/lib/python3.5/site-packages/IPython/core/interactiveshell.py:2850: DtypeWarning: Columns (12,13) have mixed types. Specify dtype option on import or set low_memory=False.\n",
" if self.run_code(code, result):\n",
"/home/drcjar/anaconda3/envs/ipfjes/lib/python3.5/site-packages/geopy/point.py:74: UserWarning: Point coordinates should be finite. NaN or inf has been passed as a coordinate. This will cause a ValueError exception in the future versions of geopy.\n",
" UserWarning)\n",
"/home/drcjar/anaconda3/envs/ipfjes/lib/python3.5/site-packages/ipykernel_launcher.py:334: SettingWithCopyWarning: \n",
"A value is trying to be set on a copy of a slice from a DataFrame.\n",
"Try using .loc[row_indexer,col_indexer] = value instead\n",
"\n",
"See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"# Average cummulative risk\n",
"\n",
"pt\n",
"case 9.40\n",
"control 9.17\n",
"dtype: float64\n"
]
}
],
"source": [
"ddf = dict_from_db() # dict of dataframes\n",
"ne = notepisodes() # not episodes\n",
"ae = asbexposed() # asbestos exposed soc2000 codes dataframe\n",
"ae2 = asbexposed2() # asbestos exposed soc1990 codes dataframe\n",
"df = study_df(ddf['ipfjes_occupationalhistory']) # main study dataframe\n",
"\n",
"# (2,'Construction'), \n",
"# this kinda asb exposed 3 - it more closely follows peto \n",
"jobgroups = ((1,'Non-construction high risk occupations'), \n",
" (2.1,'Carpenters'), \n",
" (2.2,'Plumber, electrician, painter or decorator'),(2.3, 'Other construction'),\n",
" (3,'Medium risk industrial'), (4,'Low risk industrial'),\n",
" (5,'Office and other low risk'))\n",
"jobgroups = collections.OrderedDict(jobgroups)\n",
"\n",
"jobcatlookup = makejobcatlookup()\n",
"df['jobcat'] = df['soc90'].astype(int).map(lambda x: jobcatlookup.get(x))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def resultsf(df):\n",
" print('# Number of cases and controls\\n')\n",
" print (df.groupby('pt').patient_id.nunique()) # group by participant type and count number of unique patient ids\n",
"\n",
" print('\\n# Age breakdown\\n')\n",
" print (df.drop_duplicates(subset='patient_id').groupby('pt').years.describe())\n",
" # df.boxplot(by='pt', column='years')\n",
"\n",
" print('\\n# Ethnicity\\n')\n",
" print (df.drop_duplicates(subset='patient_id').groupby('pt').ethnicity.value_counts())\n",
" # drop duplicate patient ids (they exist because our main table has a row for each job), count ethnicities\n",
"\n",
" print('\\n# Ever smoked\\n')\n",
" print (df.drop_duplicates(subset='patient_id').groupby('pt').es.value_counts())\n",
" # drop duplicate patient ids (they exist because our main table has a row for each job), count ever smokers\n",
"\n",
" print('\\n# Ever smoked(%)\\n')\n",
" print(df[df['es'] == 'Yes'].drop_duplicates(subset='patient_id').groupby('pt').es.value_counts() / df.groupby('pt').patient_id.nunique() * 100)\n",
"\n",
" print('\\n# Current smoker\\n')\n",
" print (df.drop_duplicates(subset='patient_id').groupby('pt').cs.value_counts())\n",
" # drop duplicate patient ids (they exist because our main table has a row for each job), count ever smokers\n",
"\n",
" print('\\n# Current smoker(%)\\n')\n",
" print(df[df['cs'] == 'Yes'].drop_duplicates(subset='patient_id').groupby('pt').cs.value_counts() / df.groupby('pt').patient_id.nunique() * 100)\n",
"\n",
" print('\\n# Social class (based on all jobs)\\n')\n",
" print (df.groupby('pt').ssec.describe())\n",
" # this is calculated based on job soc; currently is looking at all jobs which probably isnt desirable\n",
"\n",
" highest_ssec(df)\n",
" # same as Social class above but based on lowest social code job (which is how is normally done)\n",
"\n",
" print('\\n# Asbestos exposed job(%) of all jobs (according to Coggon)\\n')\n",
" print(df[df.soc2000.isin(ae.soc2000)].groupby('pt').patient_id.count() / df.groupby('pt').patient_id.count() * 100)\n",
" # this is calculated based on job soc; currently is looking at all jobs which probably isnt desirable\n",
"\n",
" print('\\n# Asbestos exposed job(%) of all jobs(according to Peto)\\n')\n",
" print(df[df.soc90.isin(ae2.soc1990)].groupby('pt').patient_id.count() / df.groupby('pt').patient_id.count() * 100)\n",
" # this is calculated based on job soc; currently is looking at all jobs which probably isnt desirable\n",
"\n",
" print('\\n# Patients with at least 1 asbestos exposed job(%) (according to Coggon)\\n')\n",
" print(df[df.soc2000.isin(ae.soc2000)].drop_duplicates(subset='patient_id').groupby('pt').patient_id.count() / df.groupby('pt').patient_id.nunique() * 100)\n",
"\n",
" print('\\n# Patients with at least 1 asbestos exposed job(%) (according to Peto)\\n')\n",
" print(df[df.soc90.isin(ae2.soc1990)].drop_duplicates(subset='patient_id').groupby('pt').patient_id.count() / df.groupby('pt').patient_id.nunique() * 100)\n",
"\n",
" cumrisk(df)\n",
"\n",
" # print('\\n# OR for IPF if at least 1 asbestos exposed job (according to Coggon)\\n')\n",
" # oddsratio(df)\n",
"\n",
" print('\\n# Distance from centre (by centre) in km\\n')\n",
" df['distfromcentre'] = df.apply(dist_from_centre, axis=1)\n",
" print(df[~df['gp_coords'].isnull()].drop_duplicates(subset='patient_id').groupby(['centre','pt']).distfromcentre.describe())\n",
"\n",
" print('\\n# Distance from centre (study av) in km\\n')\n",
" print(df[~df['gp_coords'].isnull()].drop_duplicates(subset='patient_id').groupby(['pt']).distfromcentre.mean())\n",
" # figures don't look right"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"# Number of cases and controls\n",
"\n",
"pt\n",
"case 310\n",
"control 129\n",
"Name: patient_id, dtype: int64\n",
"\n",
"# Age breakdown\n",
"\n",
" count mean std min 25% 50% 75% max\n",
"pt \n",
"case 310.00 74.63 7.95 50.00 70.00 75.00 80.00 93.00\n",
"control 129.00 71.72 8.24 40.00 66.00 72.00 78.00 93.00\n",
"\n",
"# Ethnicity\n",
"\n",
"pt ethnicity \n",
"case White 301\n",
" Asian / Asian British 7\n",
" Black / African/ Caribbean/ Black British 1\n",
" Mixed / Multiple ethnic groups 1\n",
"control White 127\n",
" Asian / Asian British 1\n",
" Other ethnic group 1\n",
"Name: ethnicity, dtype: int64\n",
"\n",
"# Ever smoked\n",
"\n",
"pt es \n",
"case Yes 232\n",
" No 78\n",
"control Yes 92\n",
" No 37\n",
"Name: es, dtype: int64\n",
"\n",
"# Ever smoked(%)\n",
"\n",
"pt es \n",
"case Yes 74.84\n",
"control Yes 71.32\n",
"dtype: float64\n",
"\n",
"# Current smoker\n",
"\n",
"pt cs \n",
"case No 226\n",
" Yes 6\n",
"control No 82\n",
" Yes 10\n",
"Name: cs, dtype: int64\n",
"\n",
"# Current smoker(%)\n",
"\n",
"pt cs \n",
"case Yes 1.94\n",
"control Yes 7.75\n",
"dtype: float64\n",
"\n",
"# Social class (based on all jobs)\n",
"\n",
" count mean std min 25% 50% 75% max\n",
"pt \n",
"case 1,392.00 4.65 2.02 1.10 3.00 5.00 6.00 7.00\n",
"control 528.00 4.59 2.04 1.10 3.00 5.00 6.25 7.00\n",
"\n",
"# Social class (based on lowest social class code job)\n",
"\n",
" count mean std min 25% 50% 75% max\n",
"pt \n",
"case 310.00 3.16 1.78 1.10 2.00 3.00 5.00 7.00\n",
"control 129.00 3.03 1.64 1.10 1.20 3.00 4.00 7.00\n",
"\n",
"# Asbestos exposed job(%) of all jobs (according to Coggon)\n",
"\n",
"pt\n",
"case 20.83\n",
"control 18.75\n",
"Name: patient_id, dtype: float64\n",
"\n",
"# Asbestos exposed job(%) of all jobs(according to Peto)\n",
"\n",
"pt\n",
"case 10.78\n",
"control 9.09\n",
"Name: patient_id, dtype: float64\n",
"\n",
"# Patients with at least 1 asbestos exposed job(%) (according to Coggon)\n",
"\n",
"pt\n",
"case 45.48\n",
"control 42.64\n",
"Name: patient_id, dtype: float64\n",
"\n",
"# Patients with at least 1 asbestos exposed job(%) (according to Peto)\n",
"\n",
"pt\n",
"case 24.52\n",
"control 24.81\n",
"Name: patient_id, dtype: float64\n",
"\n",
"# Average cummulative risk\n",
"\n",
"pt\n",
"case 9.40\n",
"control 9.17\n",
"dtype: float64\n",
"\n",
"# Distance from centre (by centre) in km\n",
"\n",
" count mean std min 25% 50% 75% max\n",
"centre pt \n",
"1 case 22.00 17.73 18.22 2.32 6.96 17.02 20.35 93.31\n",
" control 12.00 6.36 2.90 0.82 5.10 5.35 8.13 12.18\n",
"2 case 11.00 10.08 4.14 4.45 7.16 8.32 13.79 15.51\n",
" control 1.00 8.37 nan 8.37 8.37 8.37 8.37 8.37\n",
"3 case 33.00 27.30 38.81 2.57 5.57 16.00 31.05 210.19\n",
" control 14.00 11.40 12.09 2.39 4.94 9.31 13.01 50.79\n",
"4 case 23.00 29.32 21.31 2.35 7.65 30.49 44.42 77.25\n",
" control 20.00 20.79 51.67 0.99 5.83 7.71 10.67 238.30\n",
"5 case 30.00 38.12 30.25 3.40 18.93 27.39 45.70 106.99\n",
" control 9.00 2.81 1.81 0.79 1.50 2.27 4.07 5.56\n",
"6 case 11.00 76.08 25.28 27.03 53.76 88.60 94.81 101.71\n",
"7 case 7.00 45.83 46.51 2.73 20.38 30.40 53.84 139.22\n",
" control 6.00 13.84 5.84 3.21 12.75 14.82 17.85 19.25\n",
"8 case 29.00 23.81 16.11 4.31 11.62 18.20 34.42 58.76\n",
" control 14.00 6.38 4.01 2.46 3.82 4.54 8.24 16.74\n",
"9 case 28.00 24.08 20.70 2.34 6.58 15.04 38.85 78.19\n",
" control 16.00 12.19 10.95 1.88 5.32 8.44 17.14 43.85\n",
"10 case 6.00 13.82 15.66 4.83 5.99 7.54 11.01 45.38\n",
" control 6.00 30.80 67.94 0.85 1.95 2.76 6.18 169.41\n",
"11 case 7.00 28.10 37.66 1.90 2.82 4.23 48.51 87.91\n",
" control 5.00 24.73 22.98 3.15 3.57 20.56 41.86 54.53\n",
"13 case 7.00 7.29 6.81 1.93 2.61 5.51 8.82 20.73\n",
"14 case 18.00 30.70 21.69 5.61 12.03 24.78 45.14 84.95\n",
" control 5.00 9.93 6.57 1.69 4.59 11.37 15.69 16.32\n",
"15 case 7.00 15.97 13.37 1.43 5.95 8.94 28.33 32.88\n",
"16 case 14.00 17.77 7.51 2.62 12.69 19.00 22.12 31.27\n",
"17 case 3.00 4.34 3.03 0.84 3.45 6.06 6.09 6.12\n",
"18 case 3.00 15.01 20.68 1.93 3.09 4.25 21.55 38.85\n",
"19 case 3.00 14.09 7.47 8.63 9.83 11.04 16.82 22.61\n",
" control 4.00 15.52 9.50 2.41 11.59 18.58 22.51 22.52\n",
"\n",
"# Distance from centre (study av) in km\n",
"\n",
"pt\n",
"case 27.05\n",
"control 13.15\n",
"Name: distfromcentre, dtype: float64\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/drcjar/anaconda3/envs/ipfjes/lib/python3.5/site-packages/geopy/point.py:74: UserWarning: Point coordinates should be finite. NaN or inf has been passed as a coordinate. This will cause a ValueError exception in the future versions of geopy.\n",
" UserWarning)\n"
]
}
],
"source": [
"resultsf(df) "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"# Number of cases and controls\n",
"\n",
"pt\n",
"case 80\n",
"control 77\n",
"Name: patient_id, dtype: int64\n",
"\n",
"# Age breakdown\n",
"\n",
" count mean std min 25% 50% 75% max\n",
"pt \n",
"case 80.00 76.64 7.84 56.00 71.00 77.00 82.25 93.00\n",
"control 77.00 72.87 7.97 56.00 67.00 73.00 78.00 93.00\n",
"\n",
"# Ethnicity\n",
"\n",
"pt ethnicity \n",
"case White 80\n",
"control White 76\n",
" Other ethnic group 1\n",
"Name: ethnicity, dtype: int64\n",
"\n",
"# Ever smoked\n",
"\n",
"pt es \n",
"case Yes 63\n",
" No 17\n",
"control Yes 56\n",
" No 21\n",
"Name: es, dtype: int64\n",
"\n",
"# Ever smoked(%)\n",
"\n",
"pt es \n",
"case Yes 78.75\n",
"control Yes 72.73\n",
"dtype: float64\n",
"\n",
"# Current smoker\n",
"\n",
"pt cs \n",
"case No 62\n",
" Yes 1\n",
"control No 48\n",
" Yes 8\n",
"Name: cs, dtype: int64\n",
"\n",
"# Current smoker(%)\n",
"\n",
"pt cs \n",
"case Yes 1.25\n",
"control Yes 10.39\n",
"dtype: float64\n",
"\n",
"# Social class (based on all jobs)\n",
"\n",
" count mean std min 25% 50% 75% max\n",
"pt \n",
"case 428.00 4.92 1.94 1.10 3.00 6.00 7.00 7.00\n",
"control 326.00 4.78 1.96 1.10 3.00 5.00 7.00 7.00\n",
"\n",
"# Social class (based on lowest social class code job)\n",
"\n",
" count mean std min 25% 50% 75% max\n",
"pt \n",
"case 80.00 3.29 1.82 1.10 2.00 3.00 5.00 7.00\n",
"control 77.00 3.26 1.56 1.10 2.00 3.00 5.00 7.00\n",
"\n",
"# Asbestos exposed job(%) of all jobs (according to Coggon)\n",
"\n",
"pt\n",
"case 25.70\n",
"control 17.48\n",
"Name: patient_id, dtype: float64\n",
"\n",
"# Asbestos exposed job(%) of all jobs(according to Peto)\n",
"\n",
"pt\n",
"case 14.49\n",
"control 10.12\n",
"Name: patient_id, dtype: float64\n",
"\n",
"# Patients with at least 1 asbestos exposed job(%) (according to Coggon)\n",
"\n",
"pt\n",
"case 53.75\n",
"control 41.56\n",
"Name: patient_id, dtype: float64\n",
"\n",
"# Patients with at least 1 asbestos exposed job(%) (according to Peto)\n",
"\n",
"pt\n",
"case 31.25\n",
"control 29.87\n",
"Name: patient_id, dtype: float64\n",
"\n",
"# Average cummulative risk\n",
"\n",
"pt\n",
"case 11.06\n",
"control 8.43\n",
"dtype: float64\n",
"\n",
"# Distance from centre (by centre) in km\n",
"\n",
" count mean std min 25% 50% 75% max\n",
"centre pt \n",
"1 case 7.00 5.28 1.83 2.32 4.48 4.59 6.78 7.51\n",
" control 11.00 5.83 2.36 0.82 5.10 5.35 7.13 9.21\n",
"2 case 6.00 6.67 1.36 4.45 6.09 7.16 7.20 8.32\n",
" control 1.00 8.37 nan 8.37 8.37 8.37 8.37 8.37\n",
"3 case 12.00 5.03 2.26 2.57 3.33 4.53 5.92 9.04\n",
" control 8.00 5.57 2.85 2.39 2.89 5.42 7.44 9.81\n",
"4 case 8.00 6.23 2.52 2.35 4.91 7.11 7.57 9.57\n",
" control 15.00 6.18 2.75 0.99 4.93 5.91 8.18 9.75\n",
"5 case 2.00 4.67 1.79 3.40 4.04 4.67 5.30 5.93\n",
" control 9.00 2.81 1.81 0.79 1.50 2.27 4.07 5.56\n",
"7 case 1.00 2.73 nan 2.73 2.73 2.73 2.73 2.73\n",
" control 1.00 3.21 nan 3.21 3.21 3.21 3.21 3.21\n",
"8 case 6.00 5.86 1.63 4.31 4.61 5.44 6.75 8.44\n",
" control 11.00 4.67 1.84 2.46 3.52 4.37 5.34 8.82\n",
"9 case 10.00 5.72 2.21 2.34 4.43 5.80 6.58 9.01\n",
" control 11.00 6.21 2.45 1.88 4.36 6.08 8.44 9.55\n",
"10 case 4.00 6.42 1.56 4.83 5.54 6.19 7.07 8.50\n",
" control 5.00 3.07 2.44 0.85 1.84 2.27 3.25 7.16\n",
"11 case 4.00 2.94 1.14 1.90 2.03 2.82 3.73 4.23\n",
" control 2.00 3.36 0.29 3.15 3.26 3.36 3.46 3.57\n",
"13 case 5.00 3.72 1.88 1.93 2.18 3.04 5.51 5.94\n",
"14 case 3.00 7.01 1.52 5.61 6.20 6.78 7.71 8.63\n",
" control 2.00 3.14 2.05 1.69 2.41 3.14 3.86 4.59\n",
"15 case 4.00 5.57 3.78 1.43 2.83 5.95 8.69 8.94\n",
"16 case 2.00 5.17 3.60 2.62 3.89 5.17 6.44 7.71\n",
"17 case 3.00 4.34 3.03 0.84 3.45 6.06 6.09 6.12\n",
"18 case 2.00 3.09 1.64 1.93 2.51 3.09 3.67 4.25\n",
"19 case 1.00 8.63 nan 8.63 8.63 8.63 8.63 8.63\n",
" control 1.00 2.41 nan 2.41 2.41 2.41 2.41 2.41\n",
"\n",
"# Distance from centre (study av) in km\n",
"\n",
"pt\n",
"case 5.36\n",
"control 5.05\n",
"Name: distfromcentre, dtype: float64\n"
]
}
],
"source": [
"resultsf(df[df['distfromcentre'] < 10].copy()) # sensitivity analysis"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x7fb9a33b4898>"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztnX2cHGWV77+nu2cmQxLyBmaJkxjcgAosCTKyxCC7AqugLLrXEHVlYe9ll3v3qotvlxdXWdzPvZ8lsiu64t2VRa6w11VZUBOju6IJ3BhBcKITQkBxQCQvGMKYt4mTeek594+qnvRLVXdVd1dXddf5fj7zme7qp+s59dRMP/38znnOEVXFMAzDSC+ZuA0wDMMw4sUmAsMwjJRjE4FhGEbKsYnAMAwj5dhEYBiGkXJsIjAMw0g5NhEYhmGkHJsIDMMwUo5NBIZhGCknF7cBQTjhhBN06dKlcZthGIbRVmzduvUlVT2xVru2mAiWLl3KwMBA3GYYhmG0FSLyyyDtTBoyDMNIOTYRGIZhpBybCAzDMFKOTQSGYRgpxyYCwzCMlGMTgVEXwyNjbNt5gOGRsbhNMQyjQdoifNRIFusGd3P9/Y/TlckwMTXFJ99xJpeteHncZhmGUSe2IjBCMTwyxvX3P87RiSkOj01ydGKK6+5/3FYGhtHG2ERghGLX/lG6MqV/Nl2ZDLv2j8ZkkWEYjWITgRGKvnm9TExNlRybmJqib15vTBYZhtEoNhEYoVgwq4dPvuNMZnRlmN2TY0ZXhk++40wWzOqJ2zTDMOrEnMVGaC5b8XJWLTuBXftH6ZvXa5OAYbQ5NhEYdbFgVo9NAIbRIZg0ZBiGkXJsIjAMw0g5NhEYhmGknEh9BCLyHHAYyAOTqtovIvOBrwJLgeeANaq6P0o7DMMwDH9asSJ4o6quUNV+9/kNwEZVPQXY6D43DMMwYiIOaehtwN3u47uBt8dgg2EYhuES9USgwAMislVErnGPLVTVFwDc3y+L2AbDMAyjClHvI1ilqntE5GXAd0Xkp0Hf6E4c1wAsWbIkEuOGR8ZsU5RhGKkn0olAVfe4v18Uka8D5wB7ReQkVX1BRE4CXvR57x3AHQD9/f3abNsslbJhGIZDZNKQiMwUkdmFx8CbgCeA9cBVbrOrgHVR2eCHpVI2DMM4RpQrgoXA10Wk0M+/qup/iMiPgHtF5GrgeeDyCG3wpJBK+SjHsmgWUimbRGQYRtqIbCJQ1WeB5R7Hh4ELo+o3CJZK2TAM4xip3FlsqZQNwzCOkdrso5ZK2TAMwyG1EwFYKmXDMAxIqTRktDfDI2Ns23nAorwMo0mkekVgtB+2/8Mwmo+tCIy2wfZ/GEY02ERgtA2F/R/FFPZ/GIZRPzYRGG2D7f8wjGiwicBoG2z/h2FEgzmLjbbC9n8YRvOxicBoO2z/h2E0F5OGEk7YmPlOirHvpGsxjCRjK4IEEzZmvpNi7DvpWgwj6diKIKGEjZnvpBj7TroWw2gHbCJIKGFj5jspxr6TrsUw2gGbCBJK2Jj5Toqx76RrMYx2wCaChBI2Zr6TYuw76VoMox0Q1abXhW86/f39OjAwELcZsTA8MhYqZj5s+yTTSddiGHEgIltVtb9WO4saSjhhY+Y7Kca+k67FMJKMSUOGYRgpxyYCwzCMlGMTgWEYRsqxiaCMJKQ1SIINhmGkB3MWF5GEtAZJsMEwjHRhKwKXJKQ1SIINhmGkD5sIXJKQ1iAJNhiGkT5sInBJQlqDJNhgGEb6sInAZcGsHtac3VdybE1/X0s3NFlqBcMw4sCcxS7DI2Pcu3VXybF7B3Zx7YWntvSD2EoxGobRaiJfEYhIVkR+IiIb3Ocni8ijIvJzEfmqiHRHbUMQkqTPL5jVw/LFc20SMAyjJbRCGroWeKro+VrgNlU9BdgPXN0CG6oyPDLGwdEJxvPx6PO2b8AwjDiJVBoSkT7grcD/Aj4kIgJcAPyx2+Ru4GbgH6O0oxrFcfv5qSm6ssKMXHY6hj/qb+W2b8AwjLiJ2kfwaeA6YLb7fAFwQFUn3ee7gNg+9Yrj9o/irAZ6cvC595zF6YvmRD4JePV/3f2Ps2rZCSYLGYbRMiKThkTkUuBFVd1afNijqWdBBBG5RkQGRGRg3759kdjo5RfozmaZ09vt+0EcRsap1TZqv0QQW02WMgwjyhXBKuAyEXkLMAM4HmeFMFdEcu6qoA/Y4/VmVb0DuAOcwjRRGBg2bj+MjBOkbZT7BoL0b7KUYRgQ4YpAVW9U1T5VXQq8C9ikqu8BHgRWu82uAtZFZUMtwsTth0n/ELRtVPsGgvRv6SwMwygQxz6C64GviMj/BH4CfCEGG6YJGrdfkHEKWj4ck3HK3xOmbRT7BoL0H8ZGwzA6m5ZMBKr6EPCQ+/hZ4JxW9BuUaiURC3VzZ3ZnA8s4YSWfZpdkDNJ/1OksrN6wYbQPlmKiCusGd7Nq7SauuPNRLr19C2v6+wLJOHGnigjSf5Q2Fo/bqrWbWD+4u+FzGoYRHaIaiR+2qfT39+vAwEBL+xweGWPV2k0cnTj2rXlGV4YN7zuPI+P5QN904/5WHKT/ZtvoN24/uP4CWxkYRosRka2q2l+rneUa8sFPQz8ynmf54rmBztFsyScsQfpvto3mezCM9qOjpaFGYuQtJXR92LgZRvvRsRNBozp13Dp/u2LjZhjtR0f6CJqpU8et87crNm6GET+p9hE0U6eOW+dvV2zcDKN96EhpKAqd2nLyGIbRqXTkiqCgU19Xlken3m+olpPHMIxOpiMnAmhe6gZLFW0YRqfTsRMBNEentrh4wzA6nY6eCJpBLX+DX3TM8MgYO/Yc4tDoBMf3dnH6ouMbWpVYBE51bIwMo35sIqhBNX+Dn+9g3eBuPnzvIJNF80cuA59asyK0b8H8E7WxMTKMxujIfQRRUP6Ns1ouord+9vuMTVaOa09OePiGCwN/Y7W8PbWxMTIMf4LuI+jI8NEoWDCrh+WL51bk8y+mK5NhcOcBsuI9rFkJV4Yy6lKWfrRTqGxcY2QYnYRJQ3Xi5ztYsXgueZ3yfE9ew+1liCNvT7vJLJbbyDAax1YEdeKXU2fZwtncuno5ubKRzWXg1tXLQ8kVrc7b047lKy23kWE0jq0IGsBvr0LheDOihqIoZelHu4bKtnKMDKMTSfVE0IyQQ7+9Cgtm9XD+qSc2aqJnH1GFSrazzGK5jQyjflI7EbSbFl4gSrubnZrDMIz2IJXho+0actgqu21zlmF0BhY+WoV2DTlsld3lobKGYXQ2qZwI6tXC446vb2cN3zCM5JLKiaCekMNGS182AwuVNAwjClLpIygQVAtPmk/BNHzDMIKQ6lKVQQkacpi0+HoLlTQMo5mkUhoKS9zafBy+ibj9IUY82H1PJ6leEQQlzvj6OPY7tOseC6Mx7L6nl1T7CMLSam0+Dt9E0vwhRmuw+96ZxL6PQERmiMhjIrJNRHaIyCfc4yeLyKMi8nMR+aqIdEdlQ7NpdXx9HPsd2nWPhdEYdt/TTZQ+gjHgAlVdDqwALhaRc4G1wG2qegqwH7g6Qhvamjh8E3H7Q4x4sPuebiKbCNRhxH3a5f4ocAFwn3v8buDtUdnQ7sSxb8D2KqQTu+/pJlIfgYhkga3AMuBzwK3AD1V1mfv6YuDfVfWMaudJio8gLuLYN2B7FdKJ3ffOIhH7CFQ1D6wQkbnA14HXeDXzeq+IXANcA7BkyZLIbGwH4tg3YHsV0knhvhfCSG1CSActCR9V1QMi8hBwLjBXRHKqOgn0AXt83nMHcAc4K4JW2GkYhoWRppEoo4ZOdFcCiEgvcBHwFPAgsNptdhWwLiobDMMIRzuWKzUaJ8qooZOAB0XkceBHwHdVdQNwPfAhERkCFgBfiNAGwzBCYGGk6SSQNCQiV6vqF8qO3aKqN/i9R1UfB87yOP4scE5YQw3DiB4LI00nQVcEq0XkPYUnIvK/geYU5DUMIzFYGGk6Ceos/k/AehGZAi4Bfq2q/z06swzDiIvLVrycVctOsDDSFFF1IhCR+UVP/wz4BvAD4G9EZL6q/jpK4wzDiAcLH04XtVYEW3Hi/KXo91vdHwVeGal1HUaSNuskyZakYWNjpI2qE4GqntwqQzqdJMVmJ8mWpGFjY6SRQM5iEblcRGa7jz8mIl8TkYqIIMObJMVmJ8mWpGFjY6SVoFFDH1fVwyJyHvBmnGRx/xSdWZ1FkmKzk2RL0rCxMdJK0Ikg7/5+K/CPqroOaJs6AnGTpNhsL1vG86W2pLVcYZLuk2G0kqATwW4R+TywBvi2iPSEeG/qSVJsdsGWrqxMH8tPTfGDoZcARyNftXYTV9z5KKvWbmL94O6W2xgXSbpPhtFKAqWhFpHjgIuB7ar6cxE5CfgdVX0gagOhc9JQJyUaZXhkjNffspGxyWP3fkZXhg3vO49Lb9+S+nKFSblPhtEoTUlDLSLHq+ohYAbwkHtsPk71sfb/ZG4xSYnN3rV/lO5slrHJyeljXZkMgzsP0JXJcJSpkuO79o8mwu5WkZT7ZBitotY+gn8VkT8EXgKew9lHUMD2ESSYat9q/bTwFYvnhtLI4/zmbN/a6yfI2Nn4pota+wguBRCRQVV9bWtMMhqlVix8QQu/rqzNsoWzPY97fRDEGW9vsf71E2TsbHzTR1Afwe3A3ar6o+hNqqRTfAStYHhkjFVrNwXS+f2+9dX6Nhimj2YTZ9/tTpCxs/HtLIL6CIJG/lwAPCIiz4jI4yKy3a0zYEREvSGcYWLhF8zqYfniuRX/4H7H6+mj2dTTd1rDYcsJMna2lyKdBM0+ekmkVhglNLI0b0UsfJzx9mH7NpnjGEHGzvZSpJNAKwJV/aXXT9TGpZFG0xy0IhY+znj7MH1byohSgoyd7aVIJy0pXm8Ep7A0bySEsxX55OPMWR+072aMZacRZOysHkH6sIkgYfgtzWd2Z9m280Dgf8xWxMLHGW8fpO9OkDmiCOMMMnZB762FmXYGNhEkDK/QzjX9fVx6+xbTuUPiFybbLh9YSfdvJN0+IziBwkfjJo3ho4VvWjO7s5b2oUHa8Vtr0sM4k26f4dDs8FGjxRRCOI+M5y2cr0FqhcMmkaSHcSbdPiMcNhE0wNDew9w3sJOhvYcj66PZOrfF1LcHSfdvJN0+Ixw2EdTJTd/YzkW3beYj9z3ORbdt5qZ12yPpp5nhfGlOMd1uJD2MM+n2GeEwH0EdDO09zEW3ba44/r0Pns+yhbMj6bNRnds03fYk6f6NpNuXdpqShtrwZnDnAd/jUU0EjYZqWkx9e5L0lNhJt88IRkdLQ7X08GqvV3ttxeK5nufzO54Ewmq6jfgSzA9hGO1Fx64IasU4V3u91nuXLZzNlSuXcM8jz08fu3LlkshWA80gTEx9I/HhFltuGO1HR/oIaunh1V4HAmvpQ3sPM7jzACsWz030JFBMlCmmzQ9hGMki9n0EIrJYRB4UkadEZIeIXOseny8i3xWRn7u/5zW771oxztVeDxMfvWzhbFb3L26bSQCiTTFtseWG0Z5E6SOYBD6sqq8BzgXeKyKnATcAG1X1FGCj+7yp1NLDq73eN6+X8Xy+5LWxyTwzu7ORat9J0dWD+BL8bPUaO4stN4zkE5mPQFVfAF5wHx8WkaeAlwNvA37fbXY38BBwfTP7rqWHV3t93eBuJvKlctmUKhd/ZjMiwoxctunad5J09VpjV83WLUMvMVU0dLkMFltuGG1AS3wEIrIU2AycATyvqnOLXtuvqlXloXr3EQTRw4tfHx4Z4/W3bGRssvaYNEv7Tqqu7jV2YX0rPbkMD99g/gHDiIvE7CMQkVnA/cAHVPWQiAR93zXANQBLliypq+9aMc7lr+/aP0pWMkDe9z0FujIZduw5xJzeroY20zQrvr/ZG3u8xq6arYXHxa91Z22fgmG0A5FOBCLShTMJfElVv+Ye3isiJ6nqCyJyEvCi13tV9Q7gDnBWBFHaWaBvXi95nardEBidmOTP7xmgO9uYnNOMnC2tkpZq2Wq5ZwyjPYkyakiALwBPqeqnil5aD1zlPr4KWBeVDWFZMKuHW1cvJ1c2Kj1ZIZeBrqwwuydHT04QEcYmGy+B2GjOllaWY6xmq+WeMYz2JcoVwSrgT4DtIjLoHvsocAtwr4hcDTwPXB6hDaEplOnbsecQoCya08uR8fz0N9td+0c5ODrOe7/0Eybyk9PvayRdQyOlAVudOqKarVbi0DDakyijhrYAfg6BC6PqtxksmNXD+aee6Pva8MhY02WQenO2xJEOuJqtlnvGMNqPjs41FBVJkkGSZIthGO1JR6aYaBVJSsGbJFsMw0gGiQkf7WSSJIMkyRbDMNoLk4YoTZnQSOrqZtlgNB8bX8PwJ/UrguIY/NGJyappJKKK109SiolOxMbXMKqT6hVBeQz+5BRM5NUzHj+qeP1W7gNIIza+hlGbVE8EXmmTiwmaurrZNljq5uZh42sYtUm1NOQVg19M0NTVYSmO8IljH0CasPE1jNqkekVQHoNfnEaiPB6/WfH66wZ3s2rtJq6481FWrd3ED4Zesn0AEWL7LAyjNraPgNJv6ECo1NVh+6mWxtn2AUSH7bMw0ojtIwhBeQx+mNTVYaiWF6ha+UijcWyfhWH4k2ppqJmUx6l7xa23g15t8faGkT5sRdAEyuPU15zdx71bd1XErS+Y1cOas/u454fPT793TX9fYr6pWry9YaQTWxE0iFec+j0/fN4zbn14ZIx7t+4qef+9A7sS8e3b4u0NI73YRNAgtfYiwDE/QJJj2r1sy2aEHXsOmVRkGB2OSUMNUmsvArRHOUev6zgylufqu3/km3LDMIzOwFYEDVIcp35ct/dwfuiiUxNfznHBrB4+fulpFcf9Um4YhtE52IqgCRRKNN798HP8w6ahitfnz+yuaJvEmPYzFs1hZneWI+N5z9ejLIFpGEZ82IqgBkN7D3PfwE6G9h6umq56waweLlu+yPMcKxbPLXm+YFZPxb6BsGGbQcJVw9I3r5d8lQ2GSZGxwmIhsYZRHVsRVOGmb2wvCfUUYFZPzjdd9bKFs7ly5RLueeTYe65cuYRlC2dX7Sds2GZFuGp/H/cOVIarhqUgXV1XJS13u60GLCTWMGpjKSZ8GNp7mItu2xyobSFNROFDcmjvYQZ3HmDF4rk1J4FqaSf8UlyUt69lT1jCpNxIMmHH1jA6DUsx0SCDOw8EbluunS9bOLvmBFCgWtoJrw8rr/a17AlLmJQbSSbs2BpGWrGJwIdyXb8aXtp50CRnYdNOhA1XrUWzk7ENj4yxY89BQDh90fGxfuDWGltLRFeJjUk6sYnABy+9v5hsRjiuy1s7D6NLl+vytbR4r/ZePoIg/8TN1s/XDe7mI/+2jYm8IzfmMvCpNSti0+Srja35DiqxMUkv5iOowcAvhnnXPz/K5FTpOHVn4c6rXsfpi+ZURP/Uo0uH/SZW3r6e9zdTPx8eGeP1t2xibLL0G3hPTnj4hgtj/XbpNVbmOyjFxqQzMR9Bk+jKZenOZZgsi63PZbLM6e2u+CepV5cOmybZS8cP8/5m6+e79o+SzUjF8azEr8mXj435DiqxMUk3NhF4UKxzL5ozg/xU5aopr946vJcuPTaZZ2Z3NlT/Ueu0M7uzjOWbl+6ib15vqHGKk3ZIB95qbEzSjW0oK2Pd4G7O/duNXHnXj7jyrse4+DObeefr+ujKHvu2m8vArauXe35IF6eR6HHfk8kIl96+hfWDuwP1X1zKMsh7wrJucDeX3r4FcWXBGV2ZhtNdLJjVw62rzww8TnGS5FQfcWFjkm7MR1BENZ37W+9/A3sOjhI0GmZo72He8tktjE8G11xbodN69dGdFb79l28IHPJa6/xJiRqqhUXIVGJj0lnE7iMQkbuAS4EXVfUM99h84KvAUuA5YI2q7o/KhrBU07mPjOc5/9SXBT7XkfE8PdlMyURQS3NthU7r1UdPzj+/UFgWzOoJNU5xYuUrK7ExSSdRSkNfBC4uO3YDsFFVTwE2us9bTiH3zMAvhrlvYCcDvxhm284DzOzOeurck1N5+ub1hsrvE1RzLT5HK3Tamd1ZxiZLP/Sb3Yfl9gmGjZORFCJbEajqZhFZWnb4bcDvu4/vBh4Cro/KBi8KsdITk1Pkiz7zcxnIZTOcs3Qe3x8aLnmPInxm49Mlsfp+5SgLBNkf4BW3HWZPQb3XnskI5JWerCAZiaQPi0Wvjo2TkSQi9RG4E8GGImnogKrOLXp9v6rOq3WeZvkIguTpqRc/Ld9Pc63mD4Dm5/fx9A3kMnz7/ec1xTfg14fFoldi42S0iqA+gsRGDYnINSIyICID+/bta8o5g5SVrJfykpOFZT9QkXLaz5Zif0DfvF527R9tmmzg1V9PNuPpG6hXsoiiFGcnyidJLllqpJNW7yPYKyInqeoLInIS8KJfQ1W9A7gDnBVBMzoPkqenXop19iDL/mr+gChkg6D+h0b6braPo1PlE4vZN5JGq1cE64Gr3MdXAeta2XlxrHS2MjgIcHII9eRkeg/AjK4MPbkM5cFEGYGeXGXM9fDIGNff/zhHJ6aqlnj0i9sGAr2/kWv3ixMPansjfQSlUVuSjMXsG0kjyvDRL+M4hk8QkV3AXwO3APeKyNXA88DlUfXvR3GpyCf3HOTmbz5Zsm/guK4sn3vPWczp7Z4u23hwdIL3funHHB6bnG43szvH597zWub0dpVo+WFCQL3KVm7beSCyENJaZTKbEb7arFKcnZ7yIMklS430EWXU0Lt9Xrowqj6DUoiV7pvXyyc2PFny2ng+z+mL5gDHHLZ+S/lFc2ZwZDzP/iPjNdsGXfZHLRtUixNvVt/NiEVPgnwS9eaqIONkG7yMVpD6ncXrB3fzoXsHKSwKurLCu89ZXJHWGShN/eyGjwIcnZgqCcUsb+unbftp4OsHdwd6f1TjEVffSbIlCf6JJNhgtDdBo4ZSPxH4pZUopjysc2Z3lktv3+IZhho0BLRWCGGc3wST9C00DluSEN6ZBBuM9if2FBPtwq79o3RnM1UngoI2XQgD9dLx/dpW67eaBh7nVv8kpRmIw5Yk+CeSYIORHlI7ETjJ0Q5xaHSc8Xz1kNKCNl1IqHZodNL3PcVtq32TbUQDL9gOWlEYp7iNV/9Dew8zuPMAKxbPZd7Mbs/iNgUneRKK18exIkiCfyIJNhjpIZUTwbrB3Xy4yC+QEcc3MCOX5ehkfrrUYoE1/X1sGXqppAyjwHTK5eL2hba1tN2wJSr9bO/KCn9/+fKS8/tpyzd9Yzv3/PBY6c2MONFPfj6PvCoizrjEoVHHpZHXe286zQYjPaTOR+D4BDYyNlleelL41JoVfOjen1C+2bY7K4hIhXzUlQEEJora9+QEKG1bTdsN843Xz/aeXIaHbzjmW/DSlv/vfzmH1Z//YdXz16KVGnUSNPIk+EqSYIPRvpiPwIdd+0fJSgYoKz2ZzXB0Ik8uk2U8X/qaiOC1/yyTyZAVYaKofVYylDeupu2G0cD9bM9mZPr8ftry5p+/FKiParRSo06CRp4EX0kSbDA6n8TmGoqKvnm95LVS389PKTO6MoxPVubeUVW81k2qWnGuvE5VpLKulX66GbYXzu+nLZ9/ygmB+/Gj+DqGR8bY/PSLbH56X6BrCJPCG2pr5GHHrxNzFkWBjVM6Sd2KwCmpuJxrvzJY8uE+PjnFtV8ZJO/xif+ucxbT/4r5fOCrgxR/xv/uyfO5vH9xxf6Cf33smA7fla1M81yv9l2wvXzfw62rj53fT1vuP3kBV65cwj2P+PgI+vu4d6C2j2DBrB7WDe4u8ZfkMvCpNSt8r6H8emul8K52HYX+w4yfxeMHw8YpvaTORwDB9g4U05PL8KWrvTX2733w/OnoG6/9BT054eEbLpz+oG6G9h1n1FC1cp7F11lsR63U32F8KGHHLwm+hnbAxqkzMR9BFYLsHSgmmxFfjX1w5wFW9y/23V/Qnc2W6NrN0L6dcpAn1mzjdb5lC2eX1B8obuP3nuJj1cp5el2D1/WWE8aHEnb8kuBraAdsnNJNR/sI/PTOmd1OmGhQJvJTvGxWt+9rBZ3cS9cez5eWuZzZXemM9osPD6PXhtV2h/Ye5r6BnQztPVxxHj/tf3hkjIOjE57lPPPqfQ1epTHLqRYfX35dYePr++b1+o533Hp4lP2HPXfYcY177NJAK8e4Y1cEfnpneRw+OFp5RkAVTx/BRF752PonK44LcOPXnwCO6eRrzu4ridWfnKKkzOXoxCRaFFaUy+AZHx5Grw2r7ZbvJ7hy5RL+5m2/U1X7L+4jPzVFRijxl7zzdYt9r6G8NGbBH1ErPt7vusLE128ZeqnEzsJ4B9nrESVR6vH1nDvMvgXzJURPq8e4I30Efnrnhvedx1s/+33POPxvvf889hw8yp/e9VgVEaM63VkAYdxrNvGheA9ALfv9SmGG0XaH9h7mots2Vxy/77+ey3u+8Jin9v+t97/BN7eSX5+1SmPWio9vRi4mr3MU7nX59XTKHolGz93ofTEap5lj3PalKhvBrxTgd3b8Cq8CZSKwZWgfz+4bqdgDEAYhg0i4E3RnK0sUPvLMMOXxquWlDAvSziPPDHte6449Bz2XlYNu+cxyvvPkXs9Lz0qGQdf3UY1y+8KUxvRi1/5RcmW+iOI+FszqqZnPaceeQ2TKrqo76309QUpFNmupHmWpyqjLYFqZzeiJY4w7Uhry0jt/M5Hn1gee9mx/dGKKm7/5VMP9TuSnyPmVPvN7T5kOWy7beLUrb1PuvD06mefP7xmgO1uZGmLF4rmedtz98HOeK5m8TrFi8dyaJT7Lr6PRUpxP7D7IyFjppDE6MRk41866wd1cd9/jFSuciSnv66mVx6eZS/Uo8wg1cu5GS6wazSGOMe7IFUF5KcCeXMbTwdlsslnhry87nVzRqHZlhStXLpm2JZdxjnmVKBzae9hzEujJHWvn1SY/pdNlM3tyGVSVsUn1LPE4b2a3Z9SP1ySQy8Ctq5ezbOHs6fGc6ehfFXz80tMqIpDqLcU5PDLG32zYUdFH0NVWocyll8z1yXecWXI9QUpFNrtspt/YNENaqffcjZZYNVmoecQxxh25IoDSUoA/eX4/N3+z0tl1zI/sAAAQXklEQVQbhByQyQbT/buyGc5YNIdHP3pRRZz/tReeOq29gndGTz/Z5sZLXj39zaxam7OWzPMsq1kcBrhr/yjHdWVLXi9nRi7DDZe8mj9cvmjavsJ4PvjTF7n5mztKvq3P7M5yhlvVrZh6S3H6pdLoygYLZ/QKhTyuO8s/XfFazj/1Zb62hTlfo6GVUZaqrOfcjZZYNZpLq8e4YyeCYs5YdHz9bxYnlUQQJvNTHBwdp29eb0Wcf3k8vNeN9ZNtzlt2QqA2BSdstWXlzO4sYzXSbiOUTALFNr/x1S/jY+ueKDk+GTD1NgRb9gZJpVENrz6mVKdLkNaieNPdsoWz616q1xqPKPMIhT231zWO5ad8V4CWAyl6WjnGHSkNgaN3rlq7iSvufJQr7nqMVy2cWdd5JhXKg2Uy4oRc5spGbzyvXPMvW1m1dhPrB3eH7mvZwtlcuXJJybErVy4p2QBWq82CWT2sObuv5PU1/X3TqRkuvX0L4k5sM7oyzOjKlEhXtZahhWVr8bVPqRMiWxjvwvUX34PCsSDL3kIqjXLOWTov0D9GkD68bAPH/3LRbZv5yH2Pc9Ftm7lp3fa6lup+508qxdc4o8u5uaLKpbdvSbztRuOkJny0mRRCPvcfGeeSz2yumCigsZC68m+kYdpUC50tD5nszgrf/ss3BArnLCZIio5a6bhr9ecX5vq9D57vOyZednr1ETZVd6HPoGPUziGWQ3sP85Z/+H6JFNouthuVpDrFRJC0Bo1QHPI5oyvHhIfe3oiGXJ4GIkwbP6130EOb78llp8M5wyxDg6ToqJWOu1Z/fr6QwZ0HAk8Efn2ETdVd6DPoGLVzuoYj43l6clnG894+JqMz6UhpyEvvbCZj+SkmJvMcHB2vSF9QYDw/xc5fH/FN1VAej158bHhkjA3b9vDFHzxbkQLCj41P/orr79vGL/Yd9tSzvUImj05MlmjATnqJfWx++sWqETFO2obq4xs0HbcXwyNjzOjy1qb9fCS1zhckTYVfqu6lC44L1V+U4X9B9zLUu+chrO1+qUrCYOkq4qcjpSGA9YO7S7bLv+4V8/j+0HBT7ZrRlWEir56hqcKxPWF+qRoKsdoK08dGJyYp/6JdSAHhx5tue4in9x6Zfn7S8d3sH52siAdfP7i7JIV1sW0KNUtgFvBK05HLOOk0CmkkCqGi5SkLasXeF4/PbybyJWNbaxxqna98LLxsu+LOH7Kl6O9EgB7XJxBm34Df+Rsh6F6GRvc8BLXdL1VJFNdk1EdQaahjJwI4phF7pYduNX6pGnpyGUAr0l6U46eNb3zyV1x9z9aK459ecyYnnzi7RM/20/b9UmMETX9Rcq6iNBKF9mF8D14pIW685NXTUVFhCJumotq11aOTN7PMZFC/Q7P8E63y4bSrL6VdSHWKiQKFNARHxvM1UyREjV+qhmxG3Jj56vhp5g88udfz+CPP/roiBYNfCmm/1BiFEpjFeG1/L6Y8jUSQVBDVzt2dzXDWknmhJwG/81VLU1Ht2urZ4h/m2msRNO1As9IT1LK9mg8nKJauIjl09ERQTTdvNZNTeU+dfjI/VTNNM1Rq4wVddeUr53u2X/nK+WzYtpsv/uAX0/pt37xeTxlLmfLcK+EVt1/L/9KIFu6dxnuKg6MTdenHYctdVru2iSknpr5eLbtRHTyodt+ofyKo5u/nqwnjw7F0FcmhY6UhP918YnLKM9V01GQEPv3OFcAx3fzI+GRJiuSseKfBLtdey3XVeb05Xjg0Pv36Scd3s/fweMm5C+dYP7ibD967rUJ773/F/IoSmH4+gmIN+ehkHlWltyvXFI232ef207tr+Q7gWMlOyUig8pp+NEsHD6rd1+ufCKv537Rue0np03p8BFH4UoxjpNpH4Keb3/jmU/nb73gnnqtFV6ZyY1lYCpo7wCPPvMT7vjwY6H3FuqufrnrLH53BI8/+mpWvnM9192/3TIlRKKv5+ls2lvgkCrosULMEZoFiDRm8U2bUi1OK8yB/fs+Ap51h+whb7rK8ZKeXjymoLc3WwcPsZQhzT+rV/IPsealFM30pRimJ3kcgIhcDnwGywJ2qekszz++nm3/vZ/vqPqejoTc+ae7aP8ryxXNDOa6LY+f9YtRPPnE2b3/tYrbtPOBr6+DOA5yycDbd2Sxjk5Vx4ssXz61ZArNAkJQZ9bJgVg9zert97QzbV7mtteL8y9sHyY/kR7P3FATdyxA2PUG9+zaC7HmphaWriJ+W+whEJAt8DrgEOA14t4ic1sw+3nTaQs/jbzvzpAbO2kChgiIK36DDaKnFbWvpqn3zen1zI61YPLdtdNkkpWpuxJZ2Ge9maP5G+xKHs/gcYEhVn1XVceArwNua2cGFp/1WRW6hVy2cyRWvP7kiT085Xh/3V65cwt9dfiZdNWoNFKeXfsOyBSWvZTPCrauP5afxyxlUK9dQrbw3C2b18HeXL6c8OKhwnnZJI5ykVM2N2NIu4x0kz5XRubTcRyAiq4GLVfXP3Od/Avyuqr7P7z317iPY+OSveODJvbzptIVceNpvTR8f2nuYLUMvccKsHl79W7MZeO7XbN9ziN9ZdDx/cLrT7pFnhnlp5CjnLTuxRJ/fsecgIJy+6Hj2Hxlny9A+Tpg1g5W/7XzwF2udxf2s/O0FvvHY5RprEN01SEnBR555iZdGxj1j8NtFl43SzrDnbsSWdhnvZmj+RnJIrLNYRC4H3lw2EZyjqu8va3cNcA3AkiVLzv7lL3/ZUjsNwzDanSRvKNsFLC563gfsKW+kqneoar+q9p94YjAHpmEYhhGeOCaCHwGniMjJItINvAtYH4MdhmEYBjGEj6rqpIi8D/gOTvjoXapaWaDWMAzDaAmx7CNQ1W8D346jb8MwDKOUjs41ZBiGYdSmLVJMiMg+IEzY0AmAd7mpZJBk+5JsG5h9jWL2NUaS7fOy7RWqWjPapi0mgrCIyECQkKm4SLJ9SbYNzL5GMfsaI8n2NWKbSUOGYRgpxyYCwzCMlNOpE8EdcRtQgyTbl2TbwOxrFLOvMZJsX922daSPwDAMwwhOp64IDMMwjIB01EQgIheLyM9EZEhEbkiAPXeJyIsi8kTRsfki8l0R+bn7e16M9i0WkQdF5CkR2SEi1ybJRhGZISKPicg2175PuMdPFpFHXfu+6qYqiQURyYrIT0RkQwJte05EtovIoIgMuMcScW9dW+aKyH0i8lP3b3BlUuwTkVe541b4OSQiH0iKfa6NH3T/L54QkS+7/y91/f11zETQioI3dfBF4OKyYzcAG1X1FGCj+zwuJoEPq+prgHOB97pjlhQbx4ALVHU5sAK4WETOBdYCt7n27Qeujsk+gGuBp4qeJ8k2gDeq6oqisMKk3FtwqhT+h6q+GliOM46JsE9Vf+aO2wrgbOA3wNeTYp+IvBz4S6BfVc/ASdfzLur9+1PVjvgBVgLfKXp+I3BjAuxaCjxR9PxnwEnu45OAn8VtY5Ft64A/SKKNwHHAj4Hfxdk0k/O67y22qQ/nw+ACYANOXaNE2Ob2/xxwQtmxRNxb4HjgF7h+yqTZV2bTm4AfJMk+4OXATmA+TqqgDcCb6/3765gVAccGpsAu91jSWKiqLwC4v18Wsz0AiMhS4CzgURJkoyu9DAIvAt8FngEOqGqhmHGc9/nTwHUwXZB4AcmxDZzC1Q+IyFa3vgck596+EtgH/B9XWrtTRGYmyL5i3gV82X2cCPtUdTfwd8DzwAvAQWArdf79ddJE4FVH0kKiAiAis4D7gQ+o6qG47SlGVfPqLM/7cMqcvsarWWutAhG5FHhRVbcWH/ZoGuff4CpVfS2OXPpeETk/RlvKyQGvBf5RVc8CjhCvTOWJq7FfBvxb3LYU4/om3gacDCwCZuLc53IC/f110kQQqOBNAtgrIicBuL9fjNMYEenCmQS+pKpfcw8nykYAVT0APITjy5grIoXMuXHd51XAZSLyHE7d7QtwVghJsA0AVd3j/n4RR98+h+Tc213ALlV91H1+H87EkBT7ClwC/FhV97rPk2LfRcAvVHWfqk4AXwNeT51/f500EbRLwZv1wFXu46twdPlYEBEBvgA8paqfKnopETaKyIkiMtd93Ivzx/8U8CCwOk77VPVGVe1T1aU4f2ubVPU9SbANQERmisjswmMcnfsJEnJvVfVXwE4ReZV76ELgSRJiXxHv5pgsBMmx73ngXBE5zv0/LoxffX9/cTtimuxAeQvwNI6O/FcJsOfLOPrdBM43oKtxdOSNwM/d3/NjtO88nKXj48Cg+/OWpNgInAn8xLXvCeAm9/grgceAIZwle0/M9/n3gQ1Jss21Y5v7s6Pw/5CUe+vasgIYcO/vN4B5CbPvOGAYmFN0LEn2fQL4qfu/8S9AT71/f7az2DAMI+V0kjRkGIZh1IFNBIZhGCnHJgLDMIyUYxOBYRhGyrGJwDAMI+XkajcxjOQiIjcDIzi5azar6vd82r0deFpVn3SfvxpnI5gCq1X1mdZYHBwRWQEsUtVvx22L0dnYisDoCFT1Jr9JwOXtOFlpi5+vU9WziicBcUjK/8UKnH0dFRTtHjWMhrF9BEbbISJ/BVyJk2RwH06yrTNwNnXdJyK34OSHmQQewNl+vwEnMddB4J+BjwN5nA2I/xn4d5xdmStxJonXAx/FyR/0LVW93u17BCfd+UU4aX4/CnwSWIKTq2m9mxJ9LU42SAX+WVU/KyJnA58CZuFkifxTVX1BRB7CSfb3RmAuzsbDR3E2BfUCu4G/xcmztAgno+1LwJ8At+BsaOsBPqeqn298hI3UEdeuOPuxn3p+cHLDb8fZ9Xk8zoflR3BqP6zGScv7M459yZnr/v4ijgRUOM/NwEfcx0txMoie6z5fhLOF/0Qc+XQT8Hb3NQUucR9/HWei6cLJpz/oHv8LnPxNhXTA8902DwMnusfeCdzlPn4I+Hv38VuA77mP/xS4vczmrUCv+/wa4GPu4x6cXbonx32P7Kf9fmx5abQbbwC+rqq/ARCR8nxSh4CjwJ0i8i2clUAQfqmqP3Qfvw54SFX3uX18CTgfJw3COPAfbrvtwJiqTojIdpwJBZzVwj+pmw5YVX8tImfgrFq+66SGIYuTfqRAIeHf1qLzeLFeVUfdx28CzhSRQm6ZOcApOHn+DSMwNhEY7YivnqmqkyJyDk4SrncB78PJDFqLI0WPvdJJF5hQ1UL/UzhV1FDVqSLdXjxsFGCHqq70Oe+Y+ztP9f/Lcjvfr6rfqdLeMGqSFKeYYQRlM/BHItLrZtf8w+IX3doKc9SJtPkAjsMV4DAwO2AfjwK/JyInuHr/u4H/F8LGB4D/VpgYRKQgV50oIivdY10icnqN89Sy+TvAX7ipxBGRU91Mo4YRCpsIjLZCVX8MfBUnU+r9wPfLmswGNojI4zgf3h90j38F+B9uNazfrtHHCzilTh/Eyd75Y1UNk274Thwfw+Misg34Y1Udx/FhrHWPDeI4pKvxIHCaWzz9nT79PAn8WESeAD6PrfKNOrCoIcMwjJRjKwLDMIyUYxOBYRhGyrGJwDAMI+XYRGAYhpFybCIwDMNIOTYRGIZhpBybCAzDMFKOTQSGYRgp5/8DWQ0aDfHjWMIAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fb9ae225748>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df[(~df['gp_coords'].isnull()) & (df.risk > 0) & (df.distfromcentre < 80)].plot(x='distfromcentre', y='risk', kind='scatter')"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x7fb9ae1bb668>"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X2UG/V97/H3V9oHG9v4Ye34YK+NSQ1NAsWbZkNwDNyUcFpCXOBejJM2iek9pDRtSNM8FNOn3D6ccy+GNjQpuWkIyQ20hIRCEnNI2iZguK4hIbHp2oTQkm0K7BquwWYNtrH3QfrdP2ZkS9qRNNJqpNHM53XOnpVGo5nfzGj3p/n+vvMdc84hIiLplWl3A0REpL3UEYiIpJw6AhGRlFNHICKScuoIRERSTh2BiEjKqSMQEUk5dQQiIimnjkBEJOW62t2AMBYvXuxWrVrV7maIiHSUXbt27XfOLak1X0d0BKtWrWLnzp3tboaISEcxs2fDzKfQkIhIyqkjEBFJOXUEIiIpp45ARCTl1BGIiKScOoIEOnB4nN0jBzlweLzdTRGRDtAR6aMS3tahvWy+dw/dmQyT+Tw3XnE2lw4sb3ezRCTGdEaQIAcOj7P53j0cm8xzaHyKY5N5rrt3j84MRKQqdQQJMjp2lO5M6SHtzmQYHTvaphaJSCdQR5Ag/QtnM5nPl0ybzOfpXzi7TS0SkU6gjiBB+ub2cuMVZzOrO8O83i5mdWe48Yqz6Zvb2+6miUiMabA4YS4dWM661YsZHTtK/8LZ6gREpCZ1BAnUN7dXHYCIhKbQkIhIyqkjEBFJOXUEIiIpF+kYgZk9AxwCcsCUc27QzBYBXwdWAc8AG51zY1G2Q0REKmvFGcEvOecGnHOD/vPrgQedc6cDD/rPRUSkTdoRGroMuN1/fDtweRvaICIivqg7Agd818x2mdk1/rSlzrkXAPzfr4u4DSIiUkXU1xGsc849b2avA75nZv8W9o1+x3ENwMqVK6Nqn4hI6kV6RuCce97//SLwTeAcYJ+ZnQLg/36xwntvdc4NOucGlyxZEmUzRURSLbKOwMzmmNm8wmPgl4EfA/cBV/mzXQVsjaoNIiJSW5ShoaXAN82ssJ6vOuf+ycx+BNxtZlcDzwFXRtgGERGpIbKOwDn3M2BNwPQDwDujWq+IiNRHVxaLiKScOgIRkZRTRyAiknLqCKQuBw6Ps3vkIAcOj7e7KSLSJLoxjYS2dWgvm+/dQ3cmw2Q+z41XnM2lA8vb3SwRmSGdEUgoBw6Ps/nePRybzHNofIpjk3muu3ePzgxEEkAdgYQyOnaU7kzpx6U7k2F07GibWiQizaKOQELpXzibyXy+ZNpkPk//wtltapGINIs6Agmlb24vN15xNrO6M8zr7WJWd4Ybrzibvrm97W6aiMyQBosltEsHlrNu9WJGx47Sv3C2OgGRhFBHIHXpm9urDkAkYRQaipG45ejHrT0iEg2dEcRE3HL049YeEYmOzghiIG45+nFrj4hESx1BDMQtRz9u7RGRaKkjiIG45ejHrT0iEi11BDEQtxz9uLVHRKJlzrl2t6GmwcFBt3PnznY3I3IHDo/HKkc/bu0RkfqY2S7n3GCt+ZQ1FCNxy9GPW3tEJBoKDYmIpJw6AhGRlFNHICKScqnvCKIuo6AyDSISd6keLI66jILKNIhIJ0jtGUHUZRRUpkFEOkVqO4KoyyioTIOIdIrUdgRRl1FQmQYR6RSp7QiiLqOgMg0i0ilSX2Ii6jIKKtMgIu0SmxITZpYFdgJ7nXPrzew04GvAIuBx4APOuYmo21FJ1GUUVKZBROKuFaGhjwJPFT3fAtzsnDsdGAOubkEbRESkgkg7AjPrB94N3OY/N+BC4B5/ltuBy6Nsg4iIVBf1GcFfA9cBhfSZPuCgc27Kfz4K6AorEZE2iqwjMLP1wIvOuV3FkwNmDRytNrNrzGynme186aWXImljNe0sDVFt3SpZISLNFuVg8TrgUjO7BJgFnIx3hrDAzLr8s4J+4PmgNzvnbgVuBS9rKMJ2TtPO0hDV1q2SFSIShcjOCJxzf+Cc63fOrQLeC2xzzr0PeAjY4M92FbA1qjY0op2lIaqtWyUrRCQq7bigbDPwcTMbxhsz+FIb2lBRO0tDVFu3SlaISFRaUn3UOfcw8LD/+GfAOa1YbyPaWRqi1rpn2i5d3CYiQVJbYqKSdpaGqLbumbZr69Be1m3Zxvtve4x1W7Zx39DeiLdGRDpF6ktMVNLOb8/V1t1Iuw4cHmfdlm0cmzxxRjGrO8Mjmy/UmYFIgsWmxESnamdpiGrrbqRdhfGFY5zoCArjC+oIREShoRRQSWwRqUYdQQqoJLaIVKPQUEpcOrCcdasXK2tIRKZRR5AiKoktIkEUGgqgej4ikiY6Iyijej4ikjY6Iyiiej4ikkbqCIqono+IpJE6giL15tunfSwh7dsvkhQaIyhSyLe/rmyMICjTJu1jCWnffpEkUa2hALXq+aS9dk/at1+kU4StNaTQUIC+ub2sWbGg4j+1uIwltCs0E5ftF5HmUGioAXGo3dPO0Ewctl9EmkdnBA1od+2edqe5tnv7RaS5dEbQoHbW7olDWWnVLhJJjlR2BM266cxMavfMpA1xCc2odpFIMqSuI4hD2uNM21BPmquISC2pSh+NQ9pjM9ugm9GLSDVKHw0Qh7THZrahVpqriEgYqeoI4hBbj0MbRESKpaojiEPaYxzaICJSLFVjBAVxiK3HoQ0ikmxhxwhSlzUE8Uh7jEMbREQgZaGhJFIp6NbQfm4+7dP4SOUZQVLE4ZqINNB+bj7t03jRGUGHane9obTQfm4+7dP4iawjMLNZZvZDM9ttZk+a2Z/5008zs8fM7Kdm9nUz64mqDUkWh2si0kD7ufm0T+MnyjOCceBC59waYAC42MzOBbYANzvnTgfGgKsjbENi6XqE1tB+bj7t0/iJrCNwnsP+027/xwEXAvf4028HLo+qDUmm6xFaQ/u5+bRP4yfS6wjMLAvsAlYDnwNuAn7gnFvtv74C+Efn3FnVltPqW1V2El2P0Braz82nfRq9WFxH4JzLAQNmtgD4JvDGoNmC3mtm1wDXAKxcuTKyNnY6XY/QGoX9XEh51D+vmdNnNz5akj7qnDtoZg8D5wILzKzLOTcF9APPV3jPrcCt4J0RtKKdItUo5VGSKsqsoSX+mQBmNhu4CHgKeAjY4M92FbA1qjaINItSHiXJoswaOgV4yMz2AD8Cvuecux/YDHzczIaBPuBLEbZBpCmU8ihJFio0ZGZXO+e+VDbtBufc9ZXe45zbA7w5YPrPgHPqbahIOynlUZIs7BnBBjN7X+GJmf1vYEk0TRKJH6U8SpKFHSz+b8B9ZpYH3gW87Jz7neiaJRI/lw4sZ93qxUp5lMSp2hGY2aKipx8EvgU8Avy5mS1yzr0cZeNE4kYpj5JEtc4IduHl+VvR73f7Pw54faStk+OSdPFNkrZFJAmqdgTOudNa1RCpLEn560naFpGkCDVYbGZXmtk8//Efm9k3zGxaRpA0X5Ly15O0LSJJEjZr6E+cc4fM7DzgV/CKxf1tdM2SgiTlrydpW0SSJGxHkPN/vxv4vHNuK6D7CLRAUP76RM7LX++0W/0pF18knsJ2BHvN7AvARuA7ZtZbx3tlBgr5691ZOz4tl8/zmQeeZt2Wbbz/tsdYt2Ub9w3tbWMrw1Euvkg8hSpDbWYnARcDTzjnfmpmpwC/4Jz7btQNBJWhPnB4nLff8CDjU5WP1azuDI9svrAj/qkqa0ikNZpShtrMTnbOvQrMAh72py3Cu/tYev8zt9jo2FF6slnGp6YqzlOItXfCP1bl4ovES63rCL5qZr8K7AeewbuOoEDXEbRIUGy9XHmsPcpv3Un/Rl9t+5K+7XGj/d0ata4jWA9gZkPOuV9sTZOkXCG2fl1R/v3GwX7u3jlako9f+EOJMlc/6dcBVNu+pG973Gh/t07YMYJbgNudcz+KvknTpX2MoKD821HQt6UDh8dZt2UbxyZPnEE0a/wgymXHQbXtAxK97XGT9M9aqzT7VpUXAh8ys2eBI/glJ5xzZ8+gjVKn8th6UKy9kKt/jBN/QM0aPwiz7E4+la+2fYXHUexXmS7Kz7FMF7YjeFekrZCmiTJXv9ayO/1Uvtb26RqI1tE1J60V6loA59yzQT9RN07qF2WufrVlJ6F8RLXt0zUQraX93VqhxgjaTWME9Wt11tDukYO8/7bHODR+IsV1Xm8Xf//Bt7FmxYKmrj9qyhqKD+3vmWn2GIE0QSs/1FHm6gctO46n8o3u72r7rtpr+qfVfLrmpDXUEbRIp8fPawlKcW3nqXyr93fSj68km0JDLZCmVLg4fCtu9f5O0/GVzhI2NKTCcS2QpvLLfXN7WbNiQVv/AbZ6f6fp+EoyqSOoYXjfIe7ZOcLwvkMNLyOO8fMka/X+1vGVTqeOoIpPfesJLrp5O5+8Zw8X3bydT219oqHlKBWutVq9v3V8pdNpjKCC4X2HuOjm7dOmP/CxC1i9dF5Dy4xD/DxNWr2/dXwlbpQ+OkNDIwcrTm+0I1AqXGu1en/r+EqnSnRoKOhWjpVu71g+faDCRVCVpneCTru1pYi0RmLPCILyuh0E5npXygHftHYld3z/uePL3LR2ZcNnA+2mPHcRqSSRYwRBed29XQYY41Olud73X3se62/ZUTEHfHjfIYZGDjKwYkHHdgLKcxdJp7ZfR2BmK8zsITN7ysyeNLOP+tMXmdn3zOyn/u+FzV53UF531jJkM1YyrTuTYWjkYNUc8NVL57FhcEXHdgKgPHcRqS7KMYIp4BPOuTcC5wIfNrM3AdcDDzrnTgce9J83VVBed87lyeVLz34m83kGVixgIpebNn1OT3ZG8fQ4xeP7F84O3EbluYsIRDhG4Jx7AXjBf3zIzJ4ClgOXAe/wZ7sdeBjY3Mx1V6p7A0yb9uQLrzKZK+0g3nrqQtbfsqPheHrc4vE7hvdT3Ad2ZVCeu4gc15IxAjNbBWwHzgKec84tKHptzDlXNTzU6HUElW7lWJgG8PYbHmR8qvo+qCeeHrd4fPB4SYZHr9f4gEjStX2MoKghc4F7gd9zzr1ax/uuMbOdZrbzpZdeamjdQXVviqeNjh0la7V3QQbjyefDNT1u8fig9vRkNT4gIidE2hGYWTdeJ3Cnc+4b/uR9ZnaK//opwItB73XO3eqcG3TODS5ZsiSS9vUvnE3O5WvO99pkjt+8Yyf3De0Ntcw41Z2JW3tEJH6izBoy4EvAU865Txe9dB9wlf/4KmBrVG2opW9uLzdtWENX0V7ozhqb1q70001PGJ8Kd+vFuNWdiVt7RCR+IhsjMLPzgH8BngAKX0n/EHgMuBtYCTwHXOmce7nasqKuNXTg8Lgf+nGcuWw+fXN72f70i3zo7x/ntYkT2Tb13HoxbnVn4tYeEYle22sNOed2AFbh5XdGtd5G9M3t5YIzSsNPZy6bT95NTzcNG1KJW92ZuLVHpFn0JWfmEltiYqbidutFEZkubqnanUodQRWXDixn3erF+rYhEkMHDo+z+d49HJvMc8yPPl937x7WrV6sv9U6qSOoQSEVkXgqpEYf40RWXCFVW3+z9Ul0GepqapWAmGmJiDiVmGi1NG+7tI5So5snlWcEteKKM407pjlumeZtl9bSOF7zJLIMdTW1SkDMtERE3EpMtFKat13aR1lDlcWmxETc1CoBMdMSEXErMdFKad52aZ+gUjJSn9R1BLXiijOJOx44PM4rRyeZyKUzbqmYrUhnSl1HUKvkQqMlGbYO7WXdlm18+M7HyeXzdGctdSUdVM5CpDOlboygoFZcsZ64Y6VbY35x0+DxkhVpopitSDy0vcRE3NW6PqCe6weC8pl7slnmz+5J5T9CXXsh0llSFxpqVLXc+HbdClL5+iLSDKk9I6hHrdz4dtwKUvn6ItIsOiOoobieyaHxKY5Nlt6XoPB68X2Ps5kM61YvblubRETqoY6ghqDc+IzB7Y8+w/C+Q225FaTy9UWkmRQaqiEoN/61iTyf3TbMZ7cNs3Fwectz5/sXzubYVOmYxLGpnPL1RaQhOiOooTg3/qSe6bvr7p17+fhFZ7Q8d7487bcT0oBFJJ50RhBC4b4Etz/6DJ/dNjzt9UVzenhk84Uty50fHTvK7O4uDo1PHZ82u7tL5XdFpCE6I6igPDWzb24vl65ZFjjvgF/npLjeSZjUzsI8w/sO1ZUGGudSDkppFek8OiMIUCk1c/XSeWxau5I7vv/c8Xk3rV3J6qXzQr0/aB6AY5N5erOGZSxUGmhcy+8qpVWkM6W2xEQlYUopD+87xNDIQQZWLJjWCYR5f9A8leat1da4lHJQCWqR+FGJiQaFuf3d6qXzpnUA9bw/aJ5K81YTp1IOum2gSOdSR1AmTPy92jfxMO8PmqfSvOUaPQs4cHicJ59/BTDOXHZy0/85V9vuOJ25NCoJ2yBSiTqCMn1ze9k42F8yDrBxsP/4H3+tOHiY+H3xPDB9jKDSP5pGY/Bbh/byyX/Yffzq564MfHrjQFPj95W2e8fw/o4fN9DYhySdxgjKVIt1A6Hj4GG+QRbmmdOT5chErua8jcTgDxwe5+03bGN8qvTbem+X8ej172z6t9vi7Ybw+yuuNPYhnUxjBA2qFusuPA4TBw8Tv59pqeswMfjRsaNkMzZtetaiid8Xb9PukYMdP26gsQ9JA3UElGYB1Yrxl782PpVjTk+24rKbFVue05NlvIFbYPYvnE0uP/2sL+eiv+4gztc7hJWEbRCpJfUXlH3qW09w0c3b+eQ9e7jo5u185sGnK95usbjcRG/W+5adyRjrb9nBfUN7py27cPvK99/2GOu2bAucJ4ytQ3tZf8sOzA/jzerOhC5l0Te3l5s2nE139sRZQVcGbtqwJvJvtEm4dWUStkGkllSPEQzvO8RFN2+fNv2Bj13Awjk9Fb/JD+87xCV/s4OJqfquFWgkthy0nJ6s8Z3fPb9iCmul5USZNVRr3Z2ecZOEbZD0afsYgZl9GVgPvOicO8uftgj4OrAKeAbY6Jwbi6oNtQyNHKw4fcPgiop/8EcmcvRmMyUdQZhrBRqJLQctp7fLG1yuR9/cXi4443V1vadZ4nS9Q6OSsA0ilUQZGvoKcHHZtOuBB51zpwMP+s9b6sDhcbY//SL3736eqVxwLv/AigVV6wA1eq1AI7HlZsaok1gHKInbJNJqkZ0ROOe2m9mqssmXAe/wH98OPAxsjqoN5crz6QEMKA6OZQzu+MEz3L1zFJd3jOccs7q9/rKQP17vtQIzqQfUrOUkMRc+idsk0g6RjhH4HcH9RaGhg865BUWvjznnFtZaTjPGCCrl09cjaBwg7LUCM40tz2Q5ScyFT+I2iTRb28cIZsrMrgGuAVi5cuWMl1cpn74e5TH+Zl8rUO9ywnYOMxmviOsgqfL7pVPE9W+oWKs7gn1mdopz7gUzOwV4sdKMzrlbgVvBOyOY6Yor5dPXI0754/WERRodZ4hz6EX5/dIJ4vw3VKzV1xHcB1zlP74K2NqqFQfl0wfpzhpBJw69XfHJHz9weJzN9+7h2GSeQ+NTHJvMc929eyoOmDaSC1/vOlpN+f0Sd3H/GyoWZfroXXgDw4vNbBT4H8ANwN1mdjXwHHBlVOsPcunAchac1MNv/d1OjhbFluf0Zvn9X/55Xr9kDmB8+M7HS24DeVJ3lr/9wFu44IwlrWxuRY2ERQq32wx7itoJoZd6t0mklTrhb6ggyqyhX6vw0jujWmcYZy47mfIA0VQuz6+uWUbf3F4OHB4PCDnkePXoBAcOj8fiADYaFqlnvCLq0Euz4qbVtqkTYrOSXDP9G2rl5zd1JSYKIYWuoi3PO3hkeH/J68VlJCZzcO1dQ7ztfz7QcJmIZmpFWCTKdTSr9Ea71yFSzUz+hlr9+U1liYmgVNKg21Fe8tkdTJRddBZV+eZGtOIbQ7PX0Yq0T6WWSpzU+zfUzM9vx6ePRml07Cg92UxJR1AeuzsykaMra5RXcoiqfHMjWlH2oNnraEXctJNis5J89f4NtePzm8qOIKik8/hUnpGXjxyP3428fITJgBIUzSrfXM+3hKAb2AAtKSLX7DOCVqR9KrVUOlk7Pr+p6wgKeb2Fks7dWWMy55jI5bn2riEyBmZW8ZqD97y1cjG6etsQJre4MG9xuYupXB7HiTZGcevJetsZVrNKZrR7HSJRacfnN1VjBEGxt3rNNNZcT/yvnvY2e+wi6jh7J45viLRSMz6/GiMIEBR7q9dMY3X1xP/qaW+zxy6ijlN24viGSCu18vObqvTR/oWzmcjVV8e/3GQ+z5yebGDp4zAlkeuJ/wXNW0n52MXwvkPcs3OE4X2HAucvbmtQu6u1s9Z2JrE0dBK3SaQgVWcEO4b3E6bcUMYomc+Aub1dTObzbHxLP+tv2TEtbh42nl5P/K943lpjBMW3nvzUt57gjh88d3w5m9au5M8v+4Xjz4vbenRyCjNjVle2pN2V2rljeH/V7eyU2ir1SOI2iRRLzRhBPfH2nqwxUXTPgt6uDF/cNMiy+bNYf8uOaXHz+689L3B6tXh6VFlD1W6/uXrpvJr7oVqpbaDquEES8/eTuE2SHhojKFNPvD1jpber6clmmD+7myMTucC4+dDIwbrj6fXE/yrNG3TryWq331y9dF7N/VCt1PbuGtuZxPz9JG6TSLnUjBF44wPh4u3l4aPxqRwjL7/GnJ7stLj5RC7HrO4s41NTJdOL4/7NiC8HxfzLpx04PM6s7mzg+wdWePcDqjXuUD5eUdz2WuMbQWMw1cZUmqEZ+7baMsKO6aR9DCHt219NJ+yb1JwR7BjeX3Kz+WrKO4yJnOPau/6Vrgyc+/o+dgwfOP7aVB4+cfdQyRXIXRmOx/2bEV8OivnjKJl2/uo+fvTsGN2Z6X37+av7WL10HjB9jCJojKDwTTeo7dXGN8rHYLoyVBxTaYZm7NtaywgzppP2MYS0b381nbJvUjFG4NUWepDxqdJt7c4YkzO8WU2Q3q4Mj15/IVA9ph5GpZh/PYLWWR77Lx+vqBYbDzt/b1cGcCX7vVnx9WbE7uu9piNoTCftYwhp3/5q4rBvwo4RpCI0NDp2lGm1p4FcRJ1gT9aLIY+OHSVrpXe5KcSXoXaKJ1SO+dejeJ0FfXN7WbNiQcUP5OjYUboywW0Pem/QtmYzRtYygcso1sipcyF2X2vZnbCMThb19ndCWKWSTvpspCI0dOdjzzKem/5PP4KTAeBEDHnLPz7FkYnpMfP+hbNrpngWFGL7zWhPkEqnrj/e+wqHx0vbfnRyquJyfrz3lWnb6qW3lu7k8rY0eurcjHosYZdRrY1pr2sU5fZ3Slilkk76bCT+jGB43yHu3jnatOVtHOw/Xl+8K+PVKirct2BWd+Z4zfGxIxPcvWt6DfGPX3QGY0cmSjoBgDu+/1zgmcHCOT1kg+6dWeb81X30dgXP9yfr3xT4zb/SrfSG9x3iz+9/ctr8ZsHLP3B4nL/49k+mTf/U+jdx04Y1Feuxz+RWfs24X0KYZdRqY9pvmRnV9nfSbR4r6aTPRuLPCMKGVjLAb7z9VN7xhqX8zp27pn0bBu+Wlu9726lsvvgN0+LrxTn+fXN7uWfnSOB6Fs3pqZniWWx07CgndWdLbp1ZrDeb4a82rmH9mmVsf/pFPvT3j/Na0TfzOT1Zzlo2P/C9lVIjh0YO+iGd0n3QnQ1fCmNOb5azls9nzYoFFW8nOdPUzGbcqrLWMsK0Me23zIxi+5OSttspn41EdwTV0inLZTJwyS+cwmlL5jJVIWY0lXPHD2bxAQ06uJVCOtVCPYXXigcma6W9WgbW/lwfAGcum0++bNxjqih9s/yDWOnUdWDFAnIuoAR33oUuhVFpXvDO0oZGDrKq76TQp86VBmubUY+l2jKCtm085+3TsMtIg2ZvfyNhlbgWGeyEz0ZiO4Li+GIYU3nY8IUfsGntyuPpgkDJiH8un+eR4f2h4pSrl85j09qV3PH90nGAwjf+Sq+Vx0U3DvaTK/qDyJg3CBuU7lk4Ff343UMUMmWn8nDxZ7Yzu7trWpy1Umrk6qXzuGnDGn73a0Ml23TOqoU1S2GUp1iWb8/gqQtL0m+L014rldtoZ6y4eNvA+zyYc6y/ZUfHxaw7Sb2lmDt9PKHdEpk+OtNy0w987AIWzunhyedf4YO3/6jkGoFG0j+HRg4ysGLBtLBP+Wth2t3bZXz7I+eXhKGKBd2Gs1itVNKwpSqClC8n7HG457fOpbsrW3F72p2CB4Vbl/5LSekRpUlGL8y3/Lh8RuIo1SUmZlpuemjkIBsGVzB/dg+9XV1M5E7E5+uNU65eOq/iP87y18K0uyfrjUWsqRBiCroNZ7Gg9gedutYzjlFpOWGPwzMHXmPD4IqK2xOHWPGRiRy9XdkZfRakfmHCKnH5jHSyRHYE9ZRvDrLwpG6+8sh/0tuVmVYyYXwqx5yebOC336Dn5YPIQYrnrdXuiVzueCnoJ59/FXCcuWz+8WXXGlMIm77WyBgHTB/fCHMcqi2zWSl41b5ZhvnWGaYd1c7+6m2ThBf3NM1OOM6J7AjK44uvTeYq3noyyNV37Kr4msOLuReXZdj4ln7u3jVaEte/e+doSeloIDBuGTQmcPfO0mV99bHnjsf88w4+88DTfPWHJ6Z1Z42/unINlw4sZ8fw/pIxhWLFpS9qqTXGESRMSYpTF83m3/cdKXnfT154tepZxkxv21ctftys8uFhrwsJ0yapTzM+I1HplOOcyDGCAu9b8yv85h07p5WXaIegEs9Bsc37rz2vpNx0tZh/QW9Xhm9/ZHo57PJ5Hr2+vrhp2G+5YUpSzOnJ8u6/2TFtW8LcZrPRb1W12lVvbLkZ4ymKaUcjbt+843CcUz1GUNA3t5f5s3voyU6vDtoO5XHLSrHN4jGA3SMHq8b8C7IZCyyHXaynwnUA1VQb4yhWLU5bKEexe+Rg4MVxYW6z2WgKXrV2FR7PtHx4veMpimlHI25pmp10nBN9ZfGBw+O8cnRixrenbJbyuOWcnizjU9NLUExO5Y7XIAobZ8/lHQMrFtRVYjqMsLVewtzack5PNjBEl3PRlapklHgKAAAKOElEQVSu1q56YsvV6kLVO54S95i2NEcnHefEnhEUx+aaWVOoO2tMltUtMiCbgdndXRydnCq5jWSxt556Ig+/0L5MxiDn6M0aljEGT13Ihi/84Ph7Nq1dycbB/pJY/fmr+/j+zw6UjBHctMHL/y+OlR6byuGcK7mGoFmx9XKV4rTlt7Z8z1v7ueuHI8f3YVcG3jO4IrJS1bXix2Fiy7Xi//WOp8Q5pi3N00nHOZFjBDO5jqDLoNpwQgYCAy892Qyf3riGT/zD7qphnMI1CuXt6+nK8Plff3PgQHVPlmnXMtx/7Xk8/8oxyrOGoHaJ6TAajW+GubWl1/ajgFW8/Wez46iNZg3VE/9X1pAEaedxjvUYgZldDHwGyAK3OeduaObyZ3IdQTbj3Ry+krK7WJZMPzaZqxnPHxo5yOlL501rX282w+7RV4LXWdb9FMYRLjhjSeD8YUpg1NJofLN43ZVubem1/XVV52l2HLVa/Ljaa/XE/8OOp4RZryRHJxznlo8RmFkW+BzwLuBNwK+Z2ZuauY6ZXEfganQe1eqA1orRF+apFDu84PTFodrUijhjq8o8xz2O2uj1FCKdpB2DxecAw865nznnJoCvAZc1cwWF2Fxx+dfzV/dNm68r48Vyi+f7yysHvFtBBjDg0++Z/no2UxqjLyyvPEOmEDcOat+NV5zN4Gl905a9ae1K/vLKgZaXsq3UxnrWG2YZzVhPlArx/2K1rqcQ6TQtHyMwsw3Axc65D/rPPwC8zTl3baX3zOQ6guLY3PC+Q+wY3s+srgzLFp7EmctODrwqGE7Eexee1M3I2GssnjuLtT/XV/L6juH9LJ7bWzK9fL1jRyYqxo0rxQ6DYs3tijM2Y71h68XEOV5eb/xfJA7CjhG0oyO4EviVso7gHOfcR8rmuwa4BmDlypVvefbZZ1vaThGRThfnexaPAsUVxvqB58tncs7d6pwbdM4NLlkSPCgqIiIz146O4EfA6WZ2mpn1AO8F7mtDO0REhDakjzrnpszsWuCf8dJHv+ycm36DXBERaYm2XEfgnPsO8J12rFtEREolutaQiIjU1hElJszsJaDetKHFwP4ImhMn2sbkSMN2ahtb71TnXM1sm47oCBphZjvDpE11Mm1jcqRhO7WN8aXQkIhIyqkjEBFJuSR3BLe2uwEtoG1MjjRsp7YxphI7RiAiIuEk+YxARERCSFxHYGYXm9m/m9mwmV3f7vY0g5mtMLOHzOwpM3vSzD7qT19kZt8zs5/6vxe2u63NYGZZM/tXM7vff36amT3mb+fX/dIkHcvMFpjZPWb2b/4xXZu0Y2lmH/M/qz82s7vMbFYSjqOZfdnMXjSzHxdNCzx25vms/79oj5n9YvtaXl2iOoJW3PSmTaaATzjn3gicC3zY367rgQedc6cDD/rPk+CjwFNFz7cAN/vbOQZc3ZZWNc9ngH9yzr0BWIO3rYk5lma2HPhdYNA5dxZeKZn3kozj+BXg4rJplY7du4DT/Z9rgM+3qI11S1RHQAtuetMOzrkXnHOP+48P4f3jWI63bbf7s90OXN6eFjaPmfUD7wZu858bcCFwjz9LR2+nmZ0MXAB8CcA5N+GcO0jyjmUXMNvMuoCTgBdIwHF0zm0HXi6bXOnYXQbc4Tw/ABaY2SmtaWl9ktYRLAdGip6P+tMSw8xWAW8GHgOWOudeAK+zAF7XvpY1zV8D13HiJs19wEHn3JT/vNOP6euBl4D/44e/bjOzOSToWDrn9gJ/CTyH1wG8AuwiWcexWKVj1zH/j5LWEQTdUjgxaVFmNhe4F/g959yr7W5Ps5nZeuBF59yu4skBs3byMe0CfhH4vHPuzcAROjgMFMSPkV8GnAYsA+bghUnKdfJxDKNjPrtJ6whC3fSmE5lZN14ncKdz7hv+5H2FU03/94vtal+TrAMuNbNn8MJ6F+KdISzwQwzQ+cd0FBh1zj3mP78Hr2NI0rG8CPhP59xLzrlJ4BvA20nWcSxW6dh1zP+jpHUEibzpjR8n/xLwlHPu00Uv3Qdc5T++Ctja6rY1k3PuD5xz/c65VXjHbptz7n3AQ8AGf7aO3k7n3P8DRszs5/1J7wR+QrKO5XPAuWZ2kv/ZLWxjYo5jmUrH7j5gk589dC7wSiGEFDvOuUT9AJcATwP/AfxRu9vTpG06D++Ucg8w5P9cghc/fxD4qf97Ubvb2sRtfgdwv//49cAPgWHgH4Dedrdvhts2AOz0j+e3gIVJO5bAnwH/BvwY+DugNwnHEbgLb9xjEu8b/9WVjh1eaOhz/v+iJ/CyqNq+DUE/urJYRCTlkhYaEhGROqkjEBFJOXUEIiIpp45ARCTl1BGIiKRcV+1ZROLLzP4UOAycDGx3zj1QYb7Lgaedcz/xn78B76I1B2xwzv1Ha1ocnpkNAMucc99pd1sk2XRGIIngnPtUpU7AdzleRdri51udc28u7gT8i3/i8ncxgHe9yDRFV+iKzJiuI5COY2Z/BGzCK+j1El5Bs7PwLkC7x8xuAC7FK9/9XbwSB/fjFT97Bfgi8CdADu/iw/8O/CPela9r8TqJtwN/iHdR0Ledc5v9dR/Gu0joIrxSyn8I3AisxKsBdZ9fDn0L8Ct4ZxxfdM79jZm9Bfg0MBfYD/yGc+4FM3sYr4jgLwEL8C5SegzvwqvZwF7gfwFvxKvds8p//weAG/AuvusFPuec+8LM97CkTruvaNOPfur5Ad6Cd5XmSXjhoGHgk3h14jcAi4B/58SXnAX+76/ghYAKy/lT4JP+41V41U7P9Z8vwyuTsAQvfLoNuNx/zQHv8h9/E6+j6ca7r8CQP/238epCdfnPF/nzPAos8ae9B/iy//hh4K/8x5cAD/iPfwO4pazNu4DZ/vNrgD/2H/fiXa18WruPkX4670enl9Jpzge+6Zx7DcDMymtJvQocA24zs2/jnQmE8azzasYDvBV42Dn3kr+OO/HuIfAtYAL4J3++J4Bx59ykmT2B16GAd7bwt84vueyce9nMzsI7a/meV36HLF6pgoJCIcFdRcsJcp9z7qj/+JeBs82sUL9nPt5NUP4z5DaLABosls5UMZ7pnJsys3PwCp29F7gWr4ppLUeKHgeVDy6YdM4V1p8Hxv315ovi9hbQRgOedM6trbDccf93jup/l+Xt/Ihz7p+rzC9SU1wGxUTC2g78VzObbWbzgF8tftG/Z8N852Xa/B7egCvAIWBeyHU8BvwXM1vsx/t/Dfi/dbTxu8CHCh2DmRXCVUvMbK0/rdvMzqyxnFpt/mfgt/0S5ZjZGf5NbkTqoo5AOorzbtn5dbwKrPcC/1I2yzzgfjPbg/fP+2P+9K8Bv+/fFeznaqzjBeAP8AaPdwOPO+fqKZl8G94Ywx4z2w38uvNunboB2OJPG8IbkK7mIeBNZjZkZu+psJ6fAI/7N1P/AjrLlwYoa0hEJOV0RiAiknLqCEREUk4dgYhIyqkjEBFJOXUEIiIpp45ARCTl1BGIiKScOgIRkZT7/33Lxlvfg6MwAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fb9aea39e10>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df[(~df['gp_coords'].isnull()) & (df.risk > 0) & (df.pt == 'case')].plot(x='distfromcentre', y='risk', kind='scatter')"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x7fb9a532af60>"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAHQBJREFUeJzt3X+UXXV57/H3M5PJJCWRkMmYAkma2MldV9AQlyNNOujF4K2oCLQGrIrELmxuW+3S3noTauuPsuqqwXvF29YqGF3Gu7TAJUqyqFWQkKJYohMZAshVovwIISsJA8EEyWR+PPePvU9yZubMnH1+7L3P3vvzWmvWnLNnn72/e8+Z88x+nv39fs3dERGR4mpLuwEiIpIuBQIRkYJTIBARKTgFAhGRglMgEBEpOAUCEZGCUyAQESk4BQIRkYJTIBARKbgZaTcgigULFvjSpUvTboaISKbs3r37WXfvrrZeJgLB0qVL6e/vT7sZIiKZYmZPRllPqSERkYKL9YrAzJ4AjgKjwIi795rZfOAWYCnwBHCluz8fZztERGRqSVwRvNHdV7p7b/j8WuBud18O3B0+FxGRlKSRGroM2BI+3gJcnkIbREQkFHcgcOBOM9ttZuvDZQvd/QBA+P3llV5oZuvNrN/M+g8fPhxzM0VEiivuu4b63P0ZM3s5cJeZ/b+oL3T3m4CbAHp7exOfPWfw2BBPP/8Si86YTdeczqR3LyKSmFgDgbs/E34/ZGbfAs4HDprZme5+wMzOBA7F2YZ6bBvYz8ate+hoa2N4bIzr37GCS1eenXazRERiEVtqyMxOM7O5pcfA7wEPA9uBdeFq64BtcbWhHoPHhti4dQ/Hh8c4OjTC8eExNmzdw+CxobSbJiISizivCBYC3zKz0n6+4e7fMbMfA7ea2TXAU8AVMbahZk8//xIdbW0cZ+zkso62Np5+/iWliEQkl2ILBO7+S+C8CssHgYvi2m+jFp0xm+GxsXHLhsfGWHTG7JRaJCISL/UsnqBrTifXv2MFszramNs5g1kdbVz/jhW6GhCR3MrEWENJu3Tl2fT1LNBdQyJSCAoEU+ia06kAICKFoNSQiEjBKRCIiBScAoGISMEpEIiIFJwCgYhIwSkQiIgUnAKBiEjBKRCIiBScAoGISMEpEEhTDB4b4sF9RzRct0gGaYgJaZgm8hHJNl0RSEM0kY9I9ikQSENKE/mUK03kIyLZoEAgDdFEPiLZp0AgDdFEPiLZp2KxNEwT+YhkmwKBNIUm8hHJLqWGREQKToFARKTgFAhERApOgUBEpOAUCERECk6BQESk4BQIREQKToFARKTgFAgaoDH4RSQP1LO4ThqDX0TyQlcEddAY/CKSJwoEddAY/CKSJ7EHAjNrN7MHzOyO8PkyM9tlZo+Z2S1mNjPuNjRqYi0g72Pwq/YhUixJ1Ag+BDwKvCx8vgm4wd1vNrMvAtcAX0igHXWZqhZw/TtWsGHC8jyMvqnah0jxmLvHt3GzRcAW4FPAfwfeDhwGftPdR8xsNfBJd3/zdNvp7e31/v7+2No5lcFjQ/Rt2sHx4VP//c/qaOO+jWvomtPJ4LGhXI3BX+14RSRbzGy3u/dWWy/u1NDngA1A6ZOlCzji7iPh86eBiv9umtl6M+s3s/7Dhw/H3MzKqtUCuuZ0ct7iebn5kFTtQ6SYYgsEZnYJcMjdd5cvrrBqxUsSd7/J3Xvdvbe7uzuWNlaTpVpAM/L6WTpeEWmeOK8I+oBLzewJ4GZgDcEVwjwzK9UmFgHPxNiGhmRlPt5tA/vp27SDqzbvom/TDrYP7K9rO1k5XhFprlhrBCd3YnYh8BF3v8TM/i+wtaxYvMfd/3m616dVIyhp5VpAHHn9Vj5eEYkuao0gjZ7FG4GbzezvgAeAL6fQhpq08ny8pbz+cU4FglJev942t/LxikjzJRII3H0nsDN8/Evg/CT2WwTK64tIo9SzOEZJdMxSXl9EGqVB52KSZMesS1eeTV/PAuX1RaQuCgQxKB+UrpS737B1D309C2L7kFZeX0TqpdRQDNQxS0SyJNeBIK3B01TAFZEsyW1qKM3B00oF3DwOSici+ZPLQJBGjn4iFXBFJCtyGQji6GRVDxVwRSQLclkjUI5eRCS6XAYCdbISEYkul6khUI5eRCSq3AYCUI5eRCSKXAeCLIs6FHT5ekDqV0AawlokexQIWlDUPhDl6700PIKZMWtGe2qTzmvie5FsymWxOMvK+0AcHRrh+PAYG7bumdQ7euJ6I2MwPOrTvqYV2i0irUeBoMVEHaeo0nrVXhMnja8kkl0KBC0mah+ISutVe02c1HdDJLsUCFpM1D4QE9eb0QYd7ZZavwn13RDJrkQmr29U2pPXp0F3DYlIo1p58nqJIGofiInrpf3hq74bItmj1JCISMEpEIiIFJwCgYhIwSkQiIgUnAJBwtKaR1lEZCq6ayhBGotHRFqRrggSorF4RKRVKRAkRGPxiEirUiBISFHG4lENRCR7VCNISGksng0TagR56oWrGohINikQJCjP8yiX10COE1z5bNi6h76eBbk6TpE8UiBIWF7H4inVQEpBAE7VQPJ4vCJ5EluNwMxmmdmPzOxBM3vEzP42XL7MzHaZ2WNmdouZzYyrDZKcotRARPIozmLxELDG3c8DVgIXm9kqYBNwg7svB54HromxDZIQzUcgkl2xpYY8mOjgWPi0I/xyYA3w7nD5FuCTwBfiaockJ881EJE8i7VGYGbtwG6gB/g88AvgiLuPhKs8DVS8rcTM1gPrAZYsWRJnM6WJ8loDEcmzWPsRuPuou68EFgHnA6+stNoUr73J3Xvdvbe7uzvOZoqIFFoiHcrc/QiwE1gFzDOz0pXIIuCZJNogIiKVxXnXULeZzQsfzwbeBDwK3AOsDVdbB2yLqw0iIlJdnDWCM4EtYZ2gDbjV3e8ws58CN5vZ3wEPAF+OsQ0iIlJFnHcN7QFeU2H5LwnqBSIi0gI06JyISMEpEIiIFJwCgYhIwSkQiIgUnAJBAWiyGBGZjoahzjlNFiMi1eiKIMfKJ4s5OjTC8eExNmzdoysDERlHgSDHSpPFlCtNFiMiUhIpEJjZpDkDzOzTzW+ONJMmixGRKKJeEaw1s/eUnpjZPwMaErTFabIYEYkiarH4D4DtZjYGvAV4zt3/LL5mSbNoshgRqWbaQGBm88uevh+4HbgPuM7M5rv7c3E2TppDk8WIyHSqXRHsJpg4xsq+vy38cuAVsbZOpILBY0O6whFpomkDgbsvS6ohIlGoX4RI80W9a+gKM5sbPv4bM/ummU0aYlokTuoXIRKPqHcNfczdj5rZBcCbgS3AF+Nrlshk6hchEo+ogWA0/P424Avuvg2YGU+TJAvSGL9I/SJE4hE1EOw3sxuBK4Fvm1lnDa+VnNk2sJ++TTu4avMu+jbtYPvA/kT2q34RIvEwd6++ktlvABcDD7n7Y2Z2JvBqd78z7gYC9Pb2en9/fxK7kioGjw3Rt2kHx4dP/Wc+q6ON+zauSewDWXcNiURjZrvdvbfaetX6EbzM3X8FzAJ2hsvmA0OAPpkLqJSnP86pQFDK0yf1oax+ESLNVa0fwTfM7O3As8ATBP0IStSPoICUpxfJn2nz/O5+iQe5owF3f4W7Lyv7UhAoIOXpRfIn6lhDPzSz17n7j2NtjWSCxi8SyZeogWAN8Cdm9iTwIuGQE+6+IraWSUtTnl4kP6IGgrfE2goREUlNpEDg7k/G3RAREUmHOoWJiBScAoGISMEpEIiIFJwCgUgVaQywJ5KkqHcNiRSSJsKRItAVgcgUNBGOFEVsgcDMFpvZPWb2qJk9YmYfCpfPN7O7zOyx8PsZcbVBpBGaCEeKIs4rghHgL939lcAq4ANmdg5wLXC3uy8H7g6f55Jyy9nWjAH29B6QLIitRuDuB4AD4eOjZvYocDZwGXBhuNoWguGtN8bVjrQot5x9pQH2Nkz4PUYdWkPvAcmKSBPTNLwTs6XAvcCrgKfcfV7Zz55392nTQ1mbmKYVJm+R5qlnIhy9B6QVRJ2YJvZisZnNAbYCHw4nuYn6uvVm1m9m/YcPH46vgTFQbjlfuuZ0ct7ieTV9gOs9IFkSayAwsw6CIPB1d/9muPhgONUl4fdDlV7r7je5e6+793Z3d8fZzKbT5C2NyUNeXe8ByZI47xoy4MvAo+7+2bIfbQfWhY/XAdviakNaNHlL/bYN7Kdv0w6u2ryLvk072D6wP+0m1UXvAcmS2GoEZnYB8H3gITg5we1HgV3ArcAS4CngCnd/brptZa1GUKJJ1muTx7y63gOSpqZMXt8Id/8B4+c4LndRXPttJZq8pTalvPpxTgWCUl49q+dR7wHJAvUslpahvLpIOhQIpGUory6SDg06Jy3l0pVn09ezQHl1kQQpEEjLUV5dJFlKDTUgD/e7i4joiqBOGkdGRPJCVwR10Dj1IpInCgR10DgyIpInCgR10P3uIpInCgR10P3uIpInKhbXSfe7i0heKBA0QPe7i0geKDWUcerLICKN0hVBhqkvg4g0g64IMkp9GUSkWRQIMkp9GUSkWRQIMkp9GdKjuozkjWoEGVXqy7BhQo1AdzHFS3UZySMFggxTX4ZklddlStNpbti6h76eBTr3kmkKBBmnvgzJyeOcyiKgGkGs9h48ym39+9h78Oi06ynnnA2qy0he6YogJh+//SG+dv9TJ59fvXoJ11326knrKeecHarLSF4pEMRg78Gj44IAwNf+4ymuXrWUnoVzTy5Tzjl7VJeRPFJqKAYD+45EWq6+ANnUNaeT8xbPUxCQ3FAgiMHKxfMiLVfOWURagQJBDHoWzuXq1UvGLbuy92xePDE6riCseQ1EpBWYu6fdhqp6e3u9v78/7WbUbO/BowzsO8Lgiye44Xs/n7IgPHhsSDlnEWk6M9vt7r3V1lOxOEY9C+dyxmkz6du0Y9qCsPoCiEialBqKmQrCItLqFAhipoKwiLQ6BYKYqSAsIq1ONYIEqBOSiLSy2K4IzOwrZnbIzB4uWzbfzO4ys8fC72fEtf9Wo05IItKq4kwNfRW4eMKya4G73X05cHf4vGUUefC3Ih+7SNHFlhpy93vNbOmExZcBF4aPtwA7gY1xtaEWRR78rcjHLiLJF4sXuvsBgPD7yxPef0VFngi+yMcuIoGWvWvIzNabWb+Z9R8+fDjWfRX5Xv8iH7uIBJIOBAfN7EyA8PuhqVZ095vcvdfde7u7u2NtVJbv9W80t1/vsaumIJIfSQeC7cC68PE6YFvC+68oq/f6bxvYT9+mHVy1eRd9m3awfWB/zduo59ibsV8RaR2xDTpnZv9CUBheABwEPgHcDtwKLAGeAq5w9+eqbSupQeeyNPjb4LGhk2MYlczqaOO+jWvqanvUY2/2fkUkPqkPOufu75riRxfFtc9GZWnwt2ZPpB712DWBu0j+tGyxOIuSzJs3u64Rte1ZrqeISGUaYqJJkr4Xv5kTqdfSdk3gLpI/mpimCdLMmzda16i37Vmqp4gUVeo1giJJM2/eaF2j3rZnqZ4iItNTjaBMvTn+NPPmafUjaFQj7U67D0Pa+xdpNl0RhBrJ8aeVN29GXSKNtjfS7rTHRUp7/yJxUI2A5uX4k8ybp9WPoFGNtDvtPgxp71+kVlFrBEoN0bzxdpKcc6DZYwQl1fZG2p32uEhp718kLgoEJD/eTjNyzFm9n7+Rdqd9zGnvXyQuCgQkO95Os8bpyer4SI20O+1jTnv/InFRjaBM3OPtxJFjzur9/I20O+1jTnv/IlGpH0Ed4h5vJ47+Blm9n7+Rdqd9zGnvX6TZcp0aiut+70q54hOjY7zw0vC4fZXvf/DYEC+8dIKhkZFxr2s0x5zHe9rzeEzNoPMiccntFUGc93tPvPf++Mgoo2NjfODrPzm5L4eT+39peAQzow04MTp+W1f2Lqr7v8s83tOex2NqBp0XiVMuawRJ3e89eGyIR555gT/+Wj9DI6fOY+eMNsDHLZtK54w2fnht7e3K4z3teTymZtB5kXoVuh9BUvd7d83p5PTZM5nZ3j5ueXub0W7RTm17m9XVrjze057HY2oGnReJWy5TQ0ne711pX6NjDkS70hod87ralcd72vN4TM2g8yJxy+UVQZL3e3fN6eRjl5zDzBltnNbZzqyONj6zdgWfWXveyf3PaIOOdqOj3ca9tqPd+Mzayu2qVhjM4z3tUY+paEXTPP6upbXkskZQksT93qUiXrsZw6NjfOLt5/KeVb81bv8P73+B6+54hHZrY2RslGsueAWrf7uLc886vWK7aikM5vGe9umOqchF0zz+riVeUWsEuQ4EcYtSxKu10KfC4NR0bkRqU+hicVKiFPFqLfSpMDg1nRuReCgQNKBaES/oRDbMidHohb7ptlm03PhERS+aFv33L/HJ5V1DSema08mVr13E1+5/6uSyUgex8lz2iZHxH17TdSKbaqKYH+x9trC58ZK0JgBqBUWujUj8VCNowFQ56zs+eAGX/NMPxi0vV+vk8IBy42WKVjRVbUTqpRpBAqbKWQ/sOzJp+cR1quW1yyeKUW58vCQnAGoF+v1L3BQIIqqUnz1tZjtDI+MHD/r1iRFGRsc4MTo6cRMnTcxrV8v9tmLdYO/Bo9zWv4+9B482fdv9jw/y2Tt/Rv/jg03ZXtZz60nWRrJ+rqQ+Sg1FUCk/WxpUDqiYAmpvMwxndseMk4POzZrRPim/GzX3u31g/6TcePnAdknmjT9++0Pj6iJXr17CdZe9uinbvmrz/fxg76kA8PqeLv7P+1fVvb285NYr/f6bfRx5OVdyivoRNEml/GznDAOMoZHKNYBT67Xxpat7OfeslwFMymvX08cg7brB3oNHedMN905a/r2/eAM9C+c2tO3+xwdZe+P9k5bf9t9W0busq+bt5S23HmdtJG/nSgKqETRJpfxsu7XR3mZTvKJsvTbj9NkdJycymZjXrjX32wp1g4F9R2paXot7H3u2puXV5C23HmdtJG/nSmqT60DQjFxzpfzsyNhoOLDc9KoNKFdr7rc8f9vsvHHU3PDKxfNqWl6LNyxfEGl51LYWvd9BLXSuWk+S9ZrcBoKrNt/P2hvv5x927GXtjffz3s2TUw5RdM3pZMn88X8Mw6PwztctYlZHG7M6Kp/C9rapB5Qr33bUwcQmTnp/395nmzYQ2cRtbx/YP+W6PQvncvXqJeOWXb16ScNpIYDeZV28vmd8Cuj1PV3j0kK1tFWDtUWnc9VaanmfN0MuawTNzDVPta2Z7W184/3n8+7NuzgxOvkcds4wfnjtRZH+kKrlfqfL38Lk2kMt6s0N7z14lIF9R1i5eF5TgkC5/scHufexZ3nD8gXjfl/1trVo/Q4aoXOVvmbWa1p68nozuxj430A7sNndP93M7U+Xa641EEydn3aeGPw1nTPaOTE6MumnM9vbI09KX20y9OkmvW80Zzzdtqfbbs/CuU0PACW9y7oq/p7qbasmm49O5yp99b7PG5F4asjM2oHPA28BzgHeZWbnNHMfUXPNjWwLjJWL503Kq5Y0M78aZ/42S7nhLLVVpF5pvM/TqBGcD+x191+6+wngZuCyZu4gSq65kW21GfzPK1bQs3DuybxqZzjpTKlu0Mz8apz52yzlhrPUVpF6pfE+T7xGYGZrgYvd/f3h8/cCv+PuH5zqNfX2I5gq11yP/scH+c4jB+npPo3/eu5vjvullPKqp81s58UTo7HlV+O+jzwrueEstVWkXs14n7dshzIzuwJ484RAcL67//mE9dYD6wGWLFny2ieffDLRdoqIZF0rdyh7Glhc9nwR8MzEldz9Jnfvdffe7u7uxBonIlI0aQSCHwPLzWyZmc0E/hDYnkI7RESEFG4fdfcRM/sg8F2C20e/4u6PJN0OEREJpNKPwN2/DXw7jX2LiMh4uR1iQkREosnEEBNmdhiodtvQAqC+YSqLRecpGp2naHSeokvjXP2Wu1e92yYTgSAKM+uPcptU0ek8RaPzFI3OU3StfK6UGhIRKTgFAhGRgstTILgp7QZkhM5TNDpP0eg8Rdey5yo3NQIREalPnq4IRESkDrkIBGZ2sZn9zMz2mtm1abenVZjZV8zskJk9XLZsvpndZWaPhd/PSLONrcDMFpvZPWb2qJk9YmYfCpfrXJUxs1lm9iMzezA8T38bLl9mZrvC83RLOHRM4ZlZu5k9YGZ3hM9b9jxlPhAkMdFNhn0VuHjCsmuBu919OXB3+LzoRoC/dPdXAquAD4TvIZ2r8YaANe5+HrASuNjMVgGbgBvC8/Q8cE2KbWwlHwIeLXvesucp84GABCa6ySp3vxd4bsLiy4At4eMtwOWJNqoFufsBd/9J+PgowR/v2ehcjeOBY+HTjvDLgTXAbeHywp8nADNbBLwN2Bw+N1r4POUhEJwN7Ct7/nS4TCpb6O4HIPgABF6ecntaipktBV4D7ELnapIw3TEAHALuAn4BHHH30sTd+vsLfA7YACcnHu6ihc9THgKBVVimW6GkZmY2B9gKfNjdf5V2e1qRu4+6+0qCeUTOB15ZabVkW9VazOwS4JC77y5fXGHVljlPqYw+2mSRJrqRkw6a2ZnufsDMziT4z67wzKyDIAh83d2/GS7WuZqCux8xs50ENZV5ZjYj/G9Xf3/QB1xqZm8FZgEvI7hCaNnzlIcrAk10U5vtwLrw8TpgW4ptaQlh/vbLwKPu/tmyH+lclTGzbjObFz6eDbyJoJ5yD7A2XK3w58nd/8rdF7n7UoLPox3u/h5a+DzlokNZGHk/x6mJbj6VcpNagpn9C3AhwaiHB4FPALcDtwJLgKeAK9x9YkG5UMzsAuD7wEOcyul+lKBOoHMVMrMVBEXOdoJ/Im919+vM7BUEN2nMBx4ArnL3ofRa2jrM7ELgI+5+SSufp1wEAhERqV8eUkMiItIABQIRkYJTIBARKTgFAhGRglMgEBEpuDx0KJMCM7NPAscIOu3c6+7fm2K9y4Gfu/tPw+f/meBWPgfWuvsvkmlxdGa2EjjL3b+ddlsk33RFILng7h+fKgiELicYnbb8+TZ3f015ELBAq/xdrATeWukHZqZ/4qRp1I9AMsfM/hq4mmCwwcPAbuBVwB3ufpuZfRq4lGB46TuBbwJ3AC+EX18CPgaMAj8H/gj4N4Ken6sJgsTvEnQqM+Bf3X1juO9jBMOev4lgKOGPAtcTdDr7sLtvD4dG3wS8meCK40vu/o9m9lrgs8Ac4FngfeHwFTsJOq+9EZhHMDzxLmAvMBvYD/w9wbg+ZwFLw9e/F/g0QafBTuDz7n5j42dYCsfd9aWvzHwBryXoAfwbBOmgvcBHCOZeWEvQa/NnnPonZ174/asEKaDSdj5J0OMTgg/WMWBV+Pwsgp7E3QTp0x3A5eHPHHhL+PhbBIGmAzgPGAiX/ynBuEUzwufzw3V+CHSHy95J0AseYCfwv8LHbwW+Fz5+H/BPE9q8G5gdPl8P/E34uBPoB5al/TvSV/a+dHkpWfN64Fvu/msAM5s4rtSvgOPAZjP7V4IrgSiedPf7w8evA3a6++FwH18H3kAwPMcJ4Dvheg8BQ+4+bGYPEQQUCK4WvujhkMPu/pyZvYrgquWuYGgj2oEDZfsvDXS3u2w7lWx395fCx78HrDCz0vg1pwPLgccjHrMIoGKxZNOU+Ux3HzGz84GLCAb8+iDBhCDVvFj2uNKQwSXD7l7a/xjBrF24+1hZ3t4qtNGAR9x99RTbLY05M8r0f5cT2/nn7v7dadYXqapVimIiUd0L/L6ZzTazucDby38Yzilwugd32nyYoOAKcBSYG3Efu4D/YmYLwnz/u4B/r6GNdwJ/UgoMZlZKV3Wb2epwWYeZnVtlO9Xa/F3gT8MhtDGz/2Rmp9XQThFAgUAyxoMpJW8BBgjy8N+fsMpc4A4z20Pw4f0X4fKbgf8RTib+21X2cQD4K4Li8YPAT9y9liGDNxPUGPaY2YPAuz2YRnUtsClcNkBQkJ7OPcA5ZjZgZu+cYj8/BX5iZg8DN6KrfKmD7hoSESk4XRGIiBScAoGISMEpEIiIFJwCgYhIwSkQiIgUnAKBiEjBKRCIiBScAoGISMH9f1GswbTgFFERAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fb9a3b8d2e8>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df[(~df['gp_coords'].isnull()) & (df.risk > 0) & (df.pt == 'control')].plot(x='distfromcentre', y='risk', kind='scatter')"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"case 310\n",
"control 129\n",
"Name: pt, dtype: int64\n"
]
},
{
"data": {
"text/plain": [
"False 378\n",
"True 62\n",
"Name: gp_coords, dtype: int64"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"print(df.drop_duplicates(subset='patient_id')['pt'].value_counts())\n",
"df['pt'] = df['pt'].fillna('case') # this should be removed \n",
"df.drop_duplicates(subset='patient_id').gp_coords.isnull().value_counts()"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[<matplotlib.axes._subplots.AxesSubplot object at 0x7fb9a4b4e390>,\n",
" <matplotlib.axes._subplots.AxesSubplot object at 0x7fb9aae08470>],\n",
" [<matplotlib.axes._subplots.AxesSubplot object at 0x7fb9ab1fb630>,\n",
" <matplotlib.axes._subplots.AxesSubplot object at 0x7fb9afa79128>]],\n",
" dtype=object)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAJUCAYAAAC/nuosAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XuwpXdZJ/rvY9oMd4LQ3BLajdgKKIMXJoV6RhmCcomYWAUzeBQiFQerDiC3CTTWTME55UDwOIOcM+OMGQIERQHjJZGOKIPADKMGCGZEiJzEEJoQLkEJBNSBwHP+WG/DpunuvXpfev3WXp9P1a5e7/3ZK3s9+a7f+653VXcHAIDF+4ZFFwAAwIxgBgAwCMEMAGAQghkAwCAEMwCAQQhmAACDEMx2map6aVU9e9F1bKeq+vWqevEW9/HcqvqFbSoJWEK7sT9uhn44tnIfs92jqvYmuTrJt3b33y+6niNV1Q1J7pXkS+tmf1t337TBdr+e5LrufvEWjn2HJNcm+cfd/Teb3Q+wnEbuj1X1uXWTd0jyv/LVPvmz3f26bT6efjgwI2a7y08nuWK0pnOEx3f3ndb9HDeUbZfu/rskf5TkySfjeMBwfjqD9sf1PTHJoXxtn/y6UFZVe7Z4PP1wYILZ7vLYJO9YP6Oqzqmqq6vqs1X111X1mGn+U6vqmqq6taqur6qfXbfNParqTVV1S1X9bVX996r6hmnZfavqt6vq5qr6UFX93FaLrqpvqKpLq+rj0zHfXlUPOsa696yqK9bV9t/WLTujqn53XW1PP2Lztyc5e6v1AktpKfvjtN9fqKo3VNVvVtWtSX7qyEs8qupR01mJw9P64ZISzHaXhyT54OGJqjozyWuTXJDktCQ/mOSGafEnk/xokrskeWqSl1fV90zLnpfkxiR7Mzv1+PNJemo+v5/kfyY5PclZSZ5dVY/ehtrflGR/knsn+cskv3aM9S5Icv1U272T/Jvpdz1l2se7p9p+OMkFVXXWum2vSfLQbagVWD7L3B+T5MeT/EaSuyZ5w/FW1A+Xm2C2u5yW5NZ10+cneVV3v6W7v9zdH+3uv0qS7j7Y3X/dM+/IbFj7n07bfTHJfZJ8c3d/sbv/e88uRvwnSfZ29//V3V/o7uuT/JckTzqBGn9veqd5S1X93lTLl7v7Nd19a3f/Q5IXJ/neqrrjUbb/YpL7Jtk31XD4HfDDk9ylu18yzb8uycVH1Hbr9BwBq2cZ+uPxvLO7f3+qdaPTsfrhEhPMdpdPJ7nzuun7Jfnro61YVY+tqj+bhuJvSfK4JPeYFv/fSa5L8kfTMP6Baf43J7nvumB1S2bvFu91AjWe292nTT/nTrWcUlW/OB3rs9Oxs66e9S5M8uEkb51OPVywrrZ9R9T2/MxG1Q67c5JbTqBWYPdYhv54PB85gXX1wyW2pQsIGc5fJPm2zIavk9kL+QFHrlRV/yjJbyd5SpLLuvuL0+hVJUl335rZcP3zquo7krytqt497e9D3b1/m+t+SmaN75GZha67J7n5cD3rdfdnkzwnyXOq6iFTbe+aaru2u496bdrkQZmdZgBWz7L2x8OOvIXC5zP7BOdh60OXfrjEjJjtLlck+aF10xcneWpVnTVdYH96VT0wyalJ/lFm4ee2qnpskh85vFFV/WhVfWtVVZLPZvax7S8leVeSz1bVC6rq9tNI13dW1T+ZtntEVW3m/it3zuzj4X+TWaP5t8dasaoeX1UPmGr7zLra/jTJF6rqeVV1u6m2h1TV967b/IeS/MEm6gOW37L2x2O5OsnZVXW3qrpPkvUfNNAPl5hgtru8Nsnjqur2SdLd78p04WpmIeYdmV0XcWtmL+I3Zja8/78nuXzdfvYn+a9JPpfZC/xXuvvt3f2lJI9P8l1JPpTkU0lemdnFqMns1MCfbqLuVye5afp5f5I/Oc66357kj6fa/keSV3T3O7v7tsxG3c7M7ALeTyX51cwu3s30nDwms+cIWD3L2h+P5TWZXcD/4SRvTvL6wwv0w+XmBrO7TFW9JMknu/uXF3DsVyb5re7+w5N97I1U1XMyuzD35xddC7AY+uNXatEPByaYAQAMwqlMAIBBCGYAAIMQzAAABiGYAQAM4qTeYPYe97hHr62tncxDAifZVVdd9anu3rvoOpaBngir4UT64kkNZmtra3nPe95zMg8JnGRV9eFF17As9ERYDSfSF53KBAAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMYumC2dqBg1k7cHDRZQAAbLulC2YAALuVYAYAMAjBDABgEIIZAMAgBDMAgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwiD2LLmBevoYJANjtjJgBAAxCMAMAGIRgBgAwCMEMAGAQghkAwCAEMwCAQQhmAACDEMwAAAYhmAEADGKuYFZVz6mq91fVX1bVb1bV7arq/lV1ZVVdW1VvqKpTd7pYAIDdbMNgVlWnJ/m5JA/r7u9MckqSJyV5WZKXd/f+JJ9Ocv5OFgoAsNvNeypzT5LbV9WeJHdI8rEkj0xy6bT8kiTnbn95AACrY8Ng1t0fTfJLSQ5lFsg+k+SqJLd0923TajcmOX2nigQAWAXznMq8W5Jzktw/yX2T3DHJY4+yah9j+6dN16HdfOjQoa3UCrD09ETgeOY5lfmoJB/q7pu7+4tJfifJ9yc5bTq1mSRnJLnpaBt390Xdvb+79+7bt29bigZYVnoicDzzBLNDSR5eVXeoqkpyVpIPJHlbkidM65yX5LKdKREAYDXMc43ZlZld5P/eJO+btrkoyQuSPLeqrkty9yQX72CdAAC73p6NV0m6+0VJXnTE7OuTnLntFQEArCh3/gcAGIRgBgAwCMEMAGAQghkAwCAEMwCAQQhmAACDEMwAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBDABgEIIZAMAghg1mawcOZu3AwUWXAQBw0gwbzAAAVs2eRRewEaNmAMCqMGIGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBDABgEIIZAMAgBDMAgEEIZgAAgxDMAAAGMVcwq6rTqurSqvqrqrqmqr6vqr6pqt5SVddO/95tp4sFANjN5h0xe0WSN3f3A5M8NMk1SQ4keWt370/y1mkaAIBN2jCYVdVdkvxgkouTpLu/0N23JDknySXTapckOXenigQAWAXzjJh9S5Kbk7y6qv68ql5ZVXdMcq/u/liSTP/e82gbV9XTptOdNx86dGjbCgdYRnoicDzzBLM9Sb4nyX/q7u9O8vmcwGnL7r6ou/d39959+/ZtskyA3UFPBI5nnmB2Y5Ibu/vKafrSzILaJ6rqPkky/fvJnSkRAGA1bBjMuvvjST5SVd8+zToryQeSXJ7kvGneeUku25EKAQBWxJ4513tmktdV1alJrk/y1MxC3Rur6vwkh5I8cWdKBABYDXMFs+6+OsnDjrLorO0tBwBgdbnzPwDAIHZtMFs7cDBrBw4uugwAgLnt2mAGALBsBDMAgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwCMEMAGAQghkAwCAEMwCAQQhmAACDEMwAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABjE0geztQMHs3bg4LatBwCwKEsfzAAAdgvBDABgEIIZAMAgBDMAgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwCMEMAGAQghkAwCDmDmZVdUpV/XlVvWmavn9VXVlV11bVG6rq1J0rEwBg9zuREbNnJblm3fTLkry8u/cn+XSS87ezMACAVTNXMKuqM5KcneSV03QleWSSS6dVLkly7k4UCACwKuYdMfvlJM9P8uVp+u5Jbunu26bpG5OcfrQNq+pp0+nOmw8dOrSlYgGWnZ4IHM+GwayqfjTJJ7v7qvWzj7JqH2377r6ou/d39959+/ZtskyA3UFPBI5nzxzr/ECSH6uqxyW5XZK7ZDaCdlpV7ZlGzc5IctPOlQkAsPttGMy6+4VJXpgkVfWIJP+qu3+yqn4ryROSvD7JeUku28E657Z24OCiSwAA2JSt3MfsBUmeW1XXZXbN2cXbUxIAwGqa51TmV3T325O8fXp8fZIzt78kAIDV5M7/AACDEMwAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABnFCX8m0jHypOQCwLIyYAQAMQjADABiEYAYAMAjBDABgEIIZAMAgBDMAgEEIZlu0duCgW3IAANtCMAMAGMTKBTMjXADAqFYumAEAjEowAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBB7Fl0AZvllhcAwG5jxAwAYBCCGQDAIAQzAIBBCGYAAIMQzE6Q79oEFulw/9GHYHcSzAAABrG0t8vYKYffhd5w4dlzrbfV/QAAHGbEDABgEIIZAMAgNgxmVXW/qnpbVV1TVe+vqmdN87+pqt5SVddO/95t58sFANi95hkxuy3J87r7QUkenuTpVfXgJAeSvLW79yd56zQNAMAmbRjMuvtj3f3e6fGtSa5JcnqSc5JcMq12SZJzd6pIAIBVcELXmFXVWpLvTnJlknt198eSWXhLcs/tLg4AYJXMHcyq6k5JfjvJs7v7syew3dOm69BuPnTo0GZqBNg19ETgeOYKZlX1jZmFstd19+9Msz9RVfeZlt8nySePtm13X9Td+7t77759+7ajZoClpScCxzPPpzIrycVJrunuf79u0eVJzpsen5fksu0vb36b/aqknf6KJV/hBOwUvQV2n3nu/P8DSZ6c5H1VdfU07+eTXJjkjVV1fpJDSZ64MyUCAKyGDYNZd78zSR1j8VnbWw4AwOpy538AgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwCMEMAGAQghkAwCAEs4mvTgJ2Oz0OxieYAQAMYp7vytyVvHP86nNww4Vnb2o5ALC9jJgBAAxCMAMAGIRgtov5QAMALBfBDABgEILZBg6POhl9Akaxvhcdfny0ecDyEcwAAAYhmG0TI2rAstCrYFyCGQDAIASzY1jFEbBV/J2Br6cPwOIIZgAAgxDMAAAGIZidZE4Xbp7nDrbXIl5PXsNwfIIZAMAg9iy6gN1m3neDh9e74cKzd7IcAGCJGDEDABiEYLYkXF8FHMuxesOJjuBv97rAiRPMAAAGIZgBrJDNjHgZJYOTRzADABiEYAYAMAjBbEkd/jDAPKcYfHAAVsNOv871Edh5ghkAwCAEM1gyRkBX17wj5Jvd9mQbsabj8drjZBDMAAAGIZjxdY71rtC7RVgeh1+rG71u1y87cr3jLdvomMfazyKd6O8Di7ClYFZVj6mqD1bVdVV1YLuKAgBYRZsOZlV1SpL/mOSxSR6c5Ceq6sHbVRgAwKrZyojZmUmu6+7ru/sLSV6f5JztKYvNcKoRmNdGpxuPdVry8LzjLd9o/8e7XGIrdmKfmzkmbMVWgtnpST6ybvrGaR4AAJtQ3b25DauemOTR3f0z0/STk5zZ3c88Yr2nJbkgyWlJ7pTk/XPsfl+SQ5sqbAzqXyz1L9Z3dPftF13EqDbZE5Pl/7tQ/2Kpf7Hm7otbCWbfl+TF3f3oafqFSdLdL93UDr923zd3996t7mdR1L9Y6l+sZa9/VMv+vKp/sdS/WCdS/1ZOZb47yf6qun9VnZrkSUku38L+1rtlm/azKOpfLPUv1rLXP6plf17Vv1jqX6y569+z2SN0921V9Ywkf5jklCSv6u55h+Q38plt2s+iqH+x1L9Yy17/qJb9eVX/Yql/seauf9PBLEm6+4okV2xlH8dw0Q7s82RS/2Kpf7GWvf5RLfvzqv7FUv9izV3/pq8xAwBge/lKJgCAQQhmAACDEMwAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBDABgEIIZAMAgBDMAgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwCMEMAGAQghkAwCAEMwCAQQhmAACDEMwAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBDABgEIIZAMAgBDMAgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwCMEMAGAQghkAwCAEMwCAQQhmAACDEMwAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBDABgEIIZAMAgBDMAgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwCMEMAGAQghkAwCAEMwCAQQhmAACDEMwAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBbAlV1Uur6tmLrmMjVXVDVT1qh4/x61X14i3u47lV9QvbVBKwQPrj1uiHiyeYLZmq2pvkKUl+dYePc1KbxnS8v6+qz637ue9JOvx/TvLUqrr7SToesAN2Y388oid++Yg++ZM7cEj9cMEEs+Xz00mu6O6/X2QRVbVnB3b7+O6+07qfm3bgGF+nu/8uyR8lefLJOB6wY346u6w/ru+JSQ7la/vk67b72Prh4glmy+exSd6xfkZVnVNVV1fVZ6vqr6vqMdP8+1bV5VX1t1V1XVX9y3XbvLiq3lhVr62qW6vq/VX1sGnZryXZl+T3p3dlz6+qtarqqjq/qg4l+eNp3R+btr2lqt5eVQ/azl+2qr6hqi6tqo9vdIyqumdVXTGt97dV9d/WLTujqn63qm6uqg9V1dOP2PztSc7eztqBk26l+uN0jF+oqjdU1W9W1a1JfurISzyq6lFVdcO6af1wYILZ8nlIkg8enqiqM5O8NskFSU5L8oNJbpgW/2aSG5PcN8kTkrykqs5at68fS/L6abvLk/yHJOnuJ+dr35n94rptfijJg5I8uqq+bTrGs5PsTXJFZs3q1G38fZPkTUn2J7l3kr9M8mvHWO+CJNdPtdw7yb9Jkqo6ZdrHu5OcnuSHk1xwxHNxTZKHbnPdwMm1iv0xSX48yW8kuWuSNxxvRf1wfILZ8jktya3rps9P8qrufkt3f7m7P9rdf1VV90vyvyV5QXf/Q3dfneSV+drh6Xd29xXd/aXMws48L8QXd/fnp1MF/yLJwenYX0zyS0lun+T7N/m7/d70zvKWqvq9JJl+p9d0963d/Q9JXpzke6vqjkfZ/ouZNdl93f2F7j78zvnhSe7S3S+Z5l+X5OIkT1q37a2ZPbfA8trN/fF43tndvz/9jhudxtUPByeYLZ9PJ7nzuun7Jfnro6x33yR/293rm9SHM3uHdNjH1z3+uyS3m+P6hI8ccYwPH57o7i9Py08/cqM5ndvdp00/5yazd3dV9YtVdX1VfTbJddO69zjK9hdO9bx1OmVxwTT/m5PsWxf6bkny/MxG1Q67c5JbNlk3MIbd3B/nPe5G9MPB7cQF3Oysv0jybZkNQyezF+QDjrLeTUm+qaruvK757Evy0TmP03PMvymzUwdJkqqqzBrhvMeYx1OSPC7JIzNrcndPcnOS+rrCuj+b5DlJnlNVD0nytqp6V2bP0bXdfbzrOx6U5H9uY93Aybdq/fFY9Xw+yR3WTa8PXfrh4IyYLZ8rMruO4bCLM/to81nThfKnV9UDu/sjSf4kyUur6nZV9Y8zG9b/uk/xHMMnknzLBuu8McnZ07G/Mcnzkvyv6bhfo6oeUVXHambHc+dpn3+TWaP5t8dasaoeX1UPmBrgZ5J8afr50yRfqKrnTc/FKVX1kKr63nWb/1CSP9hEfcA4Vq0/HsvV07HvVlX3SfJz65bph4MTzJbPa5M8rqpunyTd/a4kT03y8szCyDsyG6pOkp9IspbZO7ffTfKi7n7LnMd5aZJ/PQ11/6ujrdDdH0zyU0n+3ySfSvL4zC6I/cJRVr9fZg3hRL06s/pvSvL+HKWprfPtmX0a6nNJ/keSV3T3O7v7tsxG3c7M7MLfT2V2n6O7JMn0XD4ms+cWWF6r1h+P5TWZXcD/4SRvzuxDDIfr0g8HV93bGdI5GarqJUk+2d2/vOha5lVVr0zyW939h4uu5UhV9Zwke7v75xddC7A1+uOWa9EPF0wwAwAYhFOZAACDEMwAAAYhmAEADEIwAwAYxEm9wew97nGPXltbO5mHBE6yq6666lPdvXfRdSwDPRFWw4n0xZMazNbW1vKe97znZB4SOMmq6sMbr0WiJ8KqOJG+6FQmAMAgBDMAgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwCMEMAGAQghkAwCAEMwCAQQhmAACDEMwAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBDABgEIIZAMAgBDMAgEEsVTBbO3AwawcOLroMAIAdsVTBDABgNxPMAAAGMVcwq6rnVNX7q+ovq+o3q+p2VXX/qrqyqq6tqjdU1ak7XSwAwG62YTCrqtOT/FySh3X3dyY5JcmTkrwsycu7e3+STyc5fycLBQDY7eY9lbknye2rak+SOyT5WJJHJrl0Wn5JknO3vzwAgNWxYTDr7o8m+aUkhzILZJ9JclWSW7r7tmm1G5OcvlNFAgCsgnlOZd4tyTlJ7p/kvknumOSxR1m1j7H906br0G4+dOjQVmoFWHp6InA885zKfFSSD3X3zd39xSS/k+T7k5w2ndpMkjOS3HS0jbv7ou7e39179+3bty1FAywrPRE4nj0br5JDSR5eVXdI8vdJzkryniRvS/KEJK9Pcl6Sy3aiQDeUBQBWxTzXmF2Z2UX+703yvmmbi5K8IMlzq+q6JHdPcvEO1gkAsOvNM2KW7n5RkhcdMfv6JGdue0UAACvKnf8BAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBDABgEIIZAMAgBDMAgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwCMEMAGAQghkAwCAEMwCAQQhmAACDEMwAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBDABgEHMFs6o6raouraq/qqprqur7quqbquotVXXt9O/ddrpYAIDdbN4Rs1ckeXN3PzDJQ5Nck+RAkrd29/4kb52mAQDYpA2DWVXdJckPJrk4Sbr7C919S5JzklwyrXZJknN3qkgAgFWwZ451viXJzUleXVUPTXJVkmcluVd3fyxJuvtjVXXPo21cVU9LckGS0/bu3Tt3YWsHDs69LsCy2GxPBFbDPKcy9yT5niT/qbu/O8nncwKnLbv7ou7e39179+3bt8kyAXYHPRE4nnmC2Y1JbuzuK6fpSzMLap+oqvskyfTvJ3emRACA1bBhMOvujyf5SFV9+zTrrCQfSHJ5kvOmeecluWxHKgQAWBHzXGOWJM9M8rqqOjXJ9Umemlmoe2NVnZ/kUJIn7kyJAACrYa5g1t1XJ3nYURadtb3lAACsLnf+BwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBDABgEIIZAMAgBDMAgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwCMEMAGAQSxnM1g4czNqBg4suAwBgWy1lMAMA2I0EMwCAQQhmAACDEMwAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBDABgEHMHs6o6par+vKreNE3fv6qurKprq+oNVXXqzpUJALD7nciI2bOSXLNu+mVJXt7d+5N8Osn521kYAMCqmSuYVdUZSc5O8sppupI8Msml0yqXJDl3JwoEAFgV846Y/XKS5yf58jR99yS3dPdt0/SNSU4/2oZV9bTpdOfNhw4d2lKxAMtOTwSOZ8NgVlU/muST3X3V+tlHWbWPtn13X9Td+7t77759+zZZJsDuoCcCx7NnjnV+IMmPVdXjktwuyV0yG0E7rar2TKNmZyS5aefKBADY/TYcMevuF3b3Gd29luRJSf64u38yyduSPGFa7bwkl+1YlQAAK2Ar9zF7QZLnVtV1mV1zdvH2lAQAsJrmOZX5Fd399iRvnx5fn+TM7S8JAGA1ufM/AMAgBDMAgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwCMEMAGAQghkAwCAEMwCAQQhmAACDEMwAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABiGYAQAMQjADABiEYAYAMAjBDABgEIIZAMAgBDMAgEEIZgAAgxDMAAAGIZgBAAxCMAMAGIRgBgAwCMEMAGAQghkAwCAEMwCAQQhmAACDEMwAAAaxYTCrqvtV1duq6pqqen9VPWua/01V9Zaqunb69247Xy4AwO41z4jZbUme190PSvLwJE+vqgcnOZDkrd29P8lbp2kAADZpw2DW3R/r7vdOj29Nck2S05Ock+SSabVLkpy7U0UCAKyCE7rGrKrWknx3kiuT3Ku7P5bMwluSe253cQAAq2TuYFZVd0ry20me3d2fPYHtnjZdh3bzoUOHNlMjwK6hJwLHM1cwq6pvzCyUva67f2ea/Ymqus+0/D5JPnm0bbv7ou7e39179+3btx01AywtPRE4nnk+lVlJLk5yTXf/+3WLLk9y3vT4vCSXbX95AACrY88c6/xAkicneV9VXT3N+/kkFyZ5Y1Wdn+RQkifuTIkAAKthw2DW3e9MUsdYfNb2lgMAsLrc+R8AYBC7LpitHTiYtQMHF10GAMAJ23XBDABgWQlmAACDEMwAAAYhmAEADEIwAwAYhGBIgmjzAAAIUElEQVQGADCIee78P6z1t8W44cKzF1gJAMDWGTEDABjEUo+YreemsgDAsjNiBgAwCMEMAGAQKxXMfI8mADCylQpmAAAjE8wAAAYhmAEADEIwAwAYhGAGADAIwQwAYBCCGQDAIAQzAIBBCGYAAIMQzAAABrFn0QXsNF/BBAAsCyNmAACDEMwAAAYhmAEADGLXXmM2z7Vl69e54cKzT3jfJ7INwHZZO3BQ/4FdyogZAMAgBDMAgEHs2lOZx3Myb6Gx2dOlLJbT1YzO3yjsTkbMAAAGIZgB7CJuqg3LTTADABiEYLYFawcOHvPd6fGWwbLwd7w8/HeC3WFLwayqHlNVH6yq66rqwHYVBQCwijYdzKrqlCT/Mcljkzw4yU9U1YO3qzAAgFWzlRGzM5Nc193Xd/cXkrw+yTnbUxYAR3O0U5abnXe8SzGAxdhKMDs9yUfWTd84zQMAYBOquze3YdUTkzy6u39mmn5ykjO7+5lHrPe0JBckOS3JnZK8f47d70tyaFOFjUH9i6X+xfqO7r79oosY1SZ7YrL8fxfqXyz1L9bcfXErwez7kry4ux89Tb8wSbr7pZva4dfu++bu3rvV/SyK+hdL/Yu17PWPatmfV/UvlvoX60Tq38qpzHcn2V9V96+qU5M8KcnlW9jferds034WRf2Lpf7FWvb6R7Xsz6v6F0v9izV3/Zv+rszuvq2qnpHkD5OckuRV3T3vkPxGPrNN+1kU9S+W+hdr2esf1bI/r+pfLPUv1tz1b+lLzLv7iiRXbGUfx3DRDuzzZFL/Yql/sZa9/lEt+/Oq/sVS/2LNXf+mrzEDAGB7+UomAIBBCGYAAIPY0jVm26WqHpjZtwacnqST3JTk8u6+ZqGFASyIvgiraeHXmFXVC5L8RGZf6XTjNPuMzG6/8fruvnBRtcHJUlX3yrr/AXf3JxZc0tyqqjL7irb1AeJdvejmssT0RVbdKvfEEYLZ/5fZHXG/eMT8U5O8v7v3L6ay+VTVXZO8MMm5SQ7fPO6TSS5LcmF3L8W9V1b5RbBIVfVdSf5zkrsm+eg0+4zM7nnzf3T3exdV2zyq6keS/EqSa/O19X9rZvX/0aJqW2b64hj0xZNPTxzjVOaXk9w3yYePmH+fadno3pjkj5M8ors/niRVde8k5yX5rSQ/vMDaNnSsF0FVLf2LoKqWIRi8JsnPdveV62dW1cOTvDrJQxdR1Al4RZJHdfcN62dW1f0zu5XOgxZR1C6gLy6QvrhQr8mK98QRRswek+Q/ZPYHdPhL0fdlli6f0d1vXlRt86iqD3b3t5/oslFU1dU59ovgV7t76BdBVV2T5LHHehF099DBoKquPdboR1Vd193ferJrOhFVdW2SB3X3bUfMPzXJB0avf1T64mLpi4ujJw4wYtbdb66qb8tXh1wrs2sq3t3dX1pocfP5cFU9P8klh4e5p+Hvn85XG+rI7nhk80mS7v6zqrrjIgo6QXvy1Wtw1vtokm88ybVsxh9U1cEkr81X/17ul+QpSYb+n+/kVUneXVWvz9fW/6QkFy+sqiWnLy6cvrg4K98TFz5ituyq6m5JDmT26al7ZXYu/xOZfW/oy7r7bxdY3oaq6v9J8oAc/UXwoe5+xqJqm0dVvTDJP8/sIukjXwRv7O6XLqq2eVXVY/PVT98d/h/w5dM3awyvqh6Uo9f/gYUWxsLoi4u17H1x1XuiYLbNquqfZvYu932Dn8f/ilV/EQA7S188+fTF5SWYbVFVvau7z5we/0ySpyf5vSQ/kuT3fayd41n36bVzktxzmr00n16rqsccvt5p+l3+XWb/A/7LJM9Zpk+xsX30RTZLT3Tn/+2w/nz9zyb5ke7+PzNrQD+5mJLmV1V3raoLq+qaqvqb6eeaad5pi65vI9NF0ocf37WqXllVf1FVvzFd0zK6Nyb5dJJ/1t137+67J/lnmX00/LcWWtl8XrLu8b9L8vEkj0/y7iS/upCKGIG+uEBL3hdXvicKZlv3DVV1t6q6e2YjkDcnSXd/Psltx990CCv/Iliwte5+2eFbCiRJd398GlHYt8C6NuNh3f2vu/vD3f3yJGuLLoiF0RcXa5n74sr3xIV/KnMXuGuSqzI7h99Vde/u/nhV3WmaN7q17n7Z+hnTC+LCqnrqgmrarId193dNj19eVecttJr5LPun1+5ZVc/N7G/9LlVV625g6Y3f6tIXx7FsfXHle6JgtkXdvXaMRV9O8uMnsZTNWvkXwYL9i8w+vfaO6Xlf/+m1f77Iwub0X5LceXp8SZJ7JLm5ZjcTvXphVbFQ+uLCLXNfXPme6OL/FXfEx9oPX2h5+EVwYXd/elG1zaOqXnTErF/p7sMvgl/s7qcsoq4TUbMvqz4jyZ919+fWzf/KRaQjm+o/PcmVy1g/HElfXKxV74mCGcdUVU/t7lcvuo7NWob6q+rnMvvE2jVJvivJs7r7smnZe7v7exZZ30aq6plJnpElrR9O1DL0leMZvX49UTDjOKrqUHcv28WWX7EM9VfV+5J8X3d/rqrWklya5Ne6+xVV9efd/d0LLXADy14/nKhl6CvHM3r9y95TtqN+15ituKr6i2MtyuyO3UNb9vqTnHJ4qLu7b6iqRyS5tKq+OctxkfSy1w9fZ9n7ypLXv+w9Zcv1C2bcK8mjM/to+HqV5E9OfjknbNnr/3hVfVd3X50k07usH83s+9YestjS5rLs9cPRLHtfWeb6l72nbLl+wYw3JbnT4T+i9arq7Se/nBO27PU/JUfc16m7b0vylKoa/X5DyfLXD0ez7H1lmetf9p6y5fpdYwYAMIjR72cCALAyBDMAgEEIZgAAgxDMAAAGIZgBAAzi/wd4ctUmtEAPnAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fb9a3bb5e80>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df[~df['gp_coords'].isnull()].hist(column='distfromcentre', by=['pt', 'exposed'], sharex=True, sharey=True, bins=100, figsize=(10,10))"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pt exposed\n",
"case False 18.40\n",
" True 13.15\n",
"control False 6.50\n",
" True 6.50\n",
"Name: distfromcentre, dtype: float64"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df[~df['gp_coords'].isnull()].groupby(['pt', 'exposed']).distfromcentre.median()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pt exposed\n",
"case False 27.00\n",
" True 20.09\n",
"control False 12.49\n",
" True 9.31\n",
"Name: distfromcentre, dtype: float64"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df[~df['gp_coords'].isnull()].groupby(['pt', 'exposed']).distfromcentre.mean()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"results = []\n",
"def sig(df):\n",
" for i in range(0,10000):\n",
" sample = df[df['pt'] == 'case'].sample(n=80, random_state=i)\n",
" result = sample[sample.exposed == True].patient_id.count() / sample.patient_id.count() * 100\n",
" results.append(result)\n",
" return results\n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"results = sig(df)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([1.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 1.800e+01, 0.000e+00,\n",
" 0.000e+00, 0.000e+00, 4.100e+01, 0.000e+00, 0.000e+00, 0.000e+00,\n",
" 1.040e+02, 0.000e+00, 0.000e+00, 0.000e+00, 1.590e+02, 0.000e+00,\n",
" 0.000e+00, 0.000e+00, 2.730e+02, 0.000e+00, 0.000e+00, 0.000e+00,\n",
" 4.760e+02, 0.000e+00, 0.000e+00, 0.000e+00, 7.140e+02, 0.000e+00,\n",
" 0.000e+00, 0.000e+00, 8.910e+02, 0.000e+00, 0.000e+00, 0.000e+00,\n",
" 9.960e+02, 0.000e+00, 0.000e+00, 0.000e+00, 1.090e+03, 0.000e+00,\n",
" 0.000e+00, 0.000e+00, 1.123e+03, 0.000e+00, 0.000e+00, 0.000e+00,\n",
" 1.064e+03, 0.000e+00, 0.000e+00, 0.000e+00, 8.750e+02, 0.000e+00,\n",
" 0.000e+00, 0.000e+00, 6.990e+02, 0.000e+00, 0.000e+00, 0.000e+00,\n",
" 5.230e+02, 0.000e+00, 0.000e+00, 0.000e+00, 3.670e+02, 0.000e+00,\n",
" 0.000e+00, 2.510e+02, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n",
" 1.720e+02, 0.000e+00, 0.000e+00, 0.000e+00, 8.200e+01, 0.000e+00,\n",
" 0.000e+00, 0.000e+00, 5.400e+01, 0.000e+00, 0.000e+00, 0.000e+00,\n",
" 1.500e+01, 0.000e+00, 0.000e+00, 0.000e+00, 6.000e+00, 0.000e+00,\n",
" 0.000e+00, 0.000e+00, 3.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n",
" 1.000e+00, 0.000e+00, 0.000e+00, 2.000e+00]),\n",
" array([ 7.5 , 7.8125, 8.125 , 8.4375, 8.75 , 9.0625, 9.375 ,\n",
" 9.6875, 10. , 10.3125, 10.625 , 10.9375, 11.25 , 11.5625,\n",
" 11.875 , 12.1875, 12.5 , 12.8125, 13.125 , 13.4375, 13.75 ,\n",
" 14.0625, 14.375 , 14.6875, 15. , 15.3125, 15.625 , 15.9375,\n",
" 16.25 , 16.5625, 16.875 , 17.1875, 17.5 , 17.8125, 18.125 ,\n",
" 18.4375, 18.75 , 19.0625, 19.375 , 19.6875, 20. , 20.3125,\n",
" 20.625 , 20.9375, 21.25 , 21.5625, 21.875 , 22.1875, 22.5 ,\n",
" 22.8125, 23.125 , 23.4375, 23.75 , 24.0625, 24.375 , 24.6875,\n",
" 25. , 25.3125, 25.625 , 25.9375, 26.25 , 26.5625, 26.875 ,\n",
" 27.1875, 27.5 , 27.8125, 28.125 , 28.4375, 28.75 , 29.0625,\n",
" 29.375 , 29.6875, 30. , 30.3125, 30.625 , 30.9375, 31.25 ,\n",
" 31.5625, 31.875 , 32.1875, 32.5 , 32.8125, 33.125 , 33.4375,\n",
" 33.75 , 34.0625, 34.375 , 34.6875, 35. , 35.3125, 35.625 ,\n",
" 35.9375, 36.25 , 36.5625, 36.875 , 37.1875, 37.5 , 37.8125,\n",
" 38.125 , 38.4375, 38.75 ]),\n",
" <a list of 100 Patch objects>)"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAD89JREFUeJzt3X3M3WV9x/H3ZyA+ZpaHG8LabsXZ+BAzkXTI5rIQcI4HY1kiicZpY7p0S9DhcNPqP2wuSzDZRE0Wkg6QmhiVoBmNkhlSMG5/yCzCEKwLHTK4paO34UE3o4753R/nqh7a+6E9pz3nPlzvV3Ln/H7X7zr3+f5ytffnXNfvPKSqkCT155emXYAkaToMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnTpx2Acs57bTTasOGDdMuQ5Jmyt133/39qppbqd+qDoANGzawZ8+eaZchSTMlyX8eST+XgCSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVOr+p3A0rG2YfuXn7X/8DWXTqkSafqcAUhSpwwASeqUASBJnTIAJKlTBoAkdcpXAWmmDb+qx1f0SEfHGYAkdcoAkKROGQCS1CkDQJI65UVgaQVeaNZzlTMASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6tWIAJLkxyYEk9w+1nZLk9iQPttuTW3uSfDLJviT3JTln6D5bWv8Hk2w5PqcjSTpSRzIDuAm46JC27cDuqtoI7G77ABcDG9vPNuA6GAQGcDXweuBc4OqDoSFJmo4VA6CqvgY8cUjzZmBn294JXDbU/uka+DqwJsmZwO8Dt1fVE1X1JHA7h4eKJGmCRr0GcEZV7Qdot6e39rXAo0P95lvbUu2HSbItyZ4kexYWFkYsT5K0kmP9WUBZpK2WaT+8sWoHsANg06ZNi/bRc5efuyNNzqgzgMfb0g7t9kBrnwfWD/VbBzy2TLskaUpGDYBdwMFX8mwBbh1qf1d7NdB5wNNtiegrwJuSnNwu/r6ptUmSpmTFJaAknwXOB05LMs/g1TzXADcn2Qo8Alzeut8GXALsA34EvBugqp5I8tfAN1q/j1TVoReWJUkTtGIAVNXblzh04SJ9C7hiid9zI3DjUVUnSTpufCewJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASerUsf5KSOlZ/IpHafVyBiBJnTIAJKlTLgFJx5jLXpoVzgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdWqsAEjyZ0keSHJ/ks8meUGSs5LcleTBJJ9PclLr+/y2v68d33AsTkCSNJqRAyDJWuBPgU1V9RrgBOBtwEeBa6tqI/AksLXdZSvwZFW9HLi29ZMkTcm4S0AnAi9MciLwImA/cAFwSzu+E7isbW9u+7TjFybJmI8vSRrRyAFQVd8D/hZ4hMEf/qeBu4GnquqZ1m0eWNu21wKPtvs+0/qfeujvTbItyZ4kexYWFkYtT5K0gnGWgE5m8Kz+LOBXgBcDFy/StQ7eZZljv2io2lFVm6pq09zc3KjlSZJWMM4S0BuB71bVQlX9L/BF4LeBNW1JCGAd8FjbngfWA7TjLwWeGOPxJUljGCcAHgHOS/KitpZ/IfBt4E7gra3PFuDWtr2r7dOO31FVh80AJEmTMc41gLsYXMz9JvCt9rt2AB8Erkqyj8Ea/w3tLjcAp7b2q4DtY9QtSRrTWF8JWVVXA1cf0vwQcO4ifX8MXD7O40mSjh3fCSxJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpU2N9FIT6s2H7l3++/fA1l06xEknjcgYgSZ0yACSpUy4BSVPmspqmxRmAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSp8YKgCRrktyS5DtJ9ib5rSSnJLk9yYPt9uTWN0k+mWRfkvuSnHNsTkGSNIpxZwCfAP6pql4JvBbYC2wHdlfVRmB32we4GNjYfrYB14352JKkMYwcAEl+Gfhd4AaAqvppVT0FbAZ2tm47gcva9mbg0zXwdWBNkjNHrlySNJZxZgAvAxaATyW5J8n1SV4MnFFV+wHa7emt/1rg0aH7z7c2SdIUjBMAJwLnANdV1euA/+EXyz2LySJtdVinZFuSPUn2LCwsjFGeJGk54wTAPDBfVXe1/VsYBMLjB5d22u2Bof7rh+6/Dnjs0F9aVTuqalNVbZqbmxujPEnSckYOgKr6L+DRJK9oTRcC3wZ2AVta2xbg1ra9C3hXezXQecDTB5eKJEmTN+6Xwr8X+EySk4CHgHczCJWbk2wFHgEub31vAy4B9gE/an0lSVMyVgBU1b3ApkUOXbhI3wKuGOfxJEnHju8ElqROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkTo37URCSJmzD9i//fPvhay6dYiWadc4AJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ3yjWCd801FUr+cAUhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqfGDoAkJyS5J8mX2v5ZSe5K8mCSzyc5qbU/v+3va8c3jPvYkqTRHYsZwJXA3qH9jwLXVtVG4Elga2vfCjxZVS8Hrm39JElTMlYAJFkHXApc3/YDXADc0rrsBC5r25vbPu34ha2/JGkKxp0BfBz4APCztn8q8FRVPdP254G1bXst8ChAO/506y9JmoKRAyDJm4EDVXX3cPMiXesIjg3/3m1J9iTZs7CwMGp5kqQVjDMDeAPwliQPA59jsPTzcWBNkoMfM70OeKxtzwPrAdrxlwJPHPpLq2pHVW2qqk1zc3NjlCdJWs7IAVBVH6qqdVW1AXgbcEdVvQO4E3hr67YFuLVt72r7tON3VNVhMwBJ0mQcj/cBfBC4Ksk+Bmv8N7T2G4BTW/tVwPbj8NiSpCN0TL4RrKq+Cny1bT8EnLtInx8Dlx+Lx5Mkjc93AktSp/xOYOk5zu991lKcAUhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOuUXwjyH+MUfko6GMwBJ6pQBIEmdMgAkqVMGgCR1yovAkp71AgLwRQS9cAYgSZ0yACSpUwaAJHXKAJCkTo0cAEnWJ7kzyd4kDyS5srWfkuT2JA+225Nbe5J8Msm+JPclOedYnYQk6eiNMwN4Bnh/Vb0KOA+4Ismrge3A7qraCOxu+wAXAxvbzzbgujEeW5I0ppEDoKr2V9U32/YPgb3AWmAzsLN12wlc1rY3A5+uga8Da5KcOXLlkqSxHJNrAEk2AK8D7gLOqKr9MAgJ4PTWbS3w6NDd5lubJGkKxg6AJC8BvgC8r6p+sFzXRdpqkd+3LcmeJHsWFhbGLU+StISxAiDJ8xj88f9MVX2xNT9+cGmn3R5o7fPA+qG7rwMeO/R3VtWOqtpUVZvm5ubGKU+StIxxXgUU4AZgb1V9bOjQLmBL294C3DrU/q72aqDzgKcPLhVJkiZvnM8CegPwTuBbSe5tbR8GrgFuTrIVeAS4vB27DbgE2Af8CHj3GI8tSRrTyAFQVf/C4uv6ABcu0r+AK0Z9PEnSseWngUoam19HOpv8KAhJ6pQzgFXMZ1WSjidnAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI65UdBTJAf7SBpNXEGIEmdMgAkqVMGgCR1ymsAkibO62GrgzMASeqUASBJnTIAJKlTBoAkdcqLwJJWPS8aHx/OACSpU84AxuCzEkmzzBmAJHXKAJCkThkAktSpiV8DSHIR8AngBOD6qrpm0jUsxTV9ST2ZaAAkOQH4e+D3gHngG0l2VdW3J1mHpOc2n8wdmUkvAZ0L7Kuqh6rqp8DngM0TrkGSxOSXgNYCjw7tzwOvn3ANkvQso8wYjscsY9Izl1TVcX+Qnz9Ycjnw+1X1R23/ncC5VfXeoT7bgG1t9xXAv0+swKN3GvD9aRcxplk/h1mvHzyH1WDW64dnn8OvVdXcSneY9AxgHlg/tL8OeGy4Q1XtAHZMsqhRJdlTVZumXcc4Zv0cZr1+8BxWg1mvH0Y7h0lfA/gGsDHJWUlOAt4G7JpwDZIkJjwDqKpnkrwH+AqDl4HeWFUPTLIGSdLAxN8HUFW3AbdN+nGPk5lYqlrBrJ/DrNcPnsNqMOv1wwjnMNGLwJKk1cOPgpCkThkARyjJjUkOJLl/qO2UJLcnebDdnjzNGpezRP1/meR7Se5tP5dMs8aVJFmf5M4ke5M8kOTK1j4T47BM/TMzDklekORfk/xbO4e/au1nJbmrjcHn24s8VqVlzuGmJN8dGoezp13rcpKckOSeJF9q+0c9BgbAkbsJuOiQtu3A7qraCOxu+6vVTRxeP8C1VXV2+1nt12aeAd5fVa8CzgOuSPJqZmcclqofZmccfgJcUFWvBc4GLkpyHvBRBuewEXgS2DrFGley1DkA/MXQONw7vRKPyJXA3qH9ox4DA+AIVdXXgCcOad4M7GzbO4HLJlrUUVii/plSVfur6ptt+4cM/vGvZUbGYZn6Z0YN/HfbfV77KeAC4JbWvmrHAJY9h5mRZB1wKXB92w8jjIEBMJ4zqmo/DP5zA6dPuZ5RvCfJfW2JaFUunSwmyQbgdcBdzOA4HFI/zNA4tKWHe4EDwO3AfwBPVdUzrcs8qzzYDj2Hqjo4Dn/TxuHaJM+fYokr+TjwAeBnbf9URhgDA6Bv1wG/zmAavB/4u+mWc2SSvAT4AvC+qvrBtOs5WovUP1PjUFX/V1VnM3gn/7nAqxbrNtmqjs6h55DkNcCHgFcCvwmcAnxwiiUuKcmbgQNVdfdw8yJdVxwDA2A8jyc5E6DdHphyPUelqh5v/xF+BvwDg//Mq1qS5zH44/mZqvpia56ZcVis/lkcB4Cqegr4KoPrGWuSHHxf0WEf8bJaDZ3DRW2JrqrqJ8CnWL3j8AbgLUkeZvCJyhcwmBEc9RgYAOPZBWxp21uAW6dYy1E7+Eez+QPg/qX6rgZtnfMGYG9VfWzo0EyMw1L1z9I4JJlLsqZtvxB4I4NrGXcCb23dVu0YwJLn8J2hJxFhsH6+Ksehqj5UVeuqagODj9O5o6rewQhj4BvBjlCSzwLnM/jEvceBq4F/BG4GfhV4BLi8qlblhdYl6j+fwbJDAQ8Df3xwLX01SvI7wD8D3+IXa58fZrCOvurHYZn6386MjEOS32BwgfEEBk8gb66qjyR5GYNno6cA9wB/2J5JrzrLnMMdwByD5ZR7gT8Zuli8KiU5H/jzqnrzKGNgAEhSp1wCkqROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXq/wEXlrMUf/ClbwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fb9b21304a8>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.hist(results, bins=100)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x7fb9bdfd3eb8>"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl8VfWd//HXJwlJ2BJICFtC2JFdlgBqFZdWiytWoWJtxXnYop3admo3nc5Y68/OtNP5lU6n/jpatWotIm6VtliqotQFaVjCvoU1IUAChACBQJbP7497aO/EhNxAknuTvJ+Px33knO/5nnM/5zxy7+d+v+ec7zF3R0REJC7aAYiISGxQQhAREUAJQUREAkoIIiICKCGIiEhACUFERAAlBBERCSghiIgIoIQgIiKBhGgH0Bg9evTwAQMGRDsMEZFWZeXKlQfdPaOheq0qIQwYMIAVK1ZEOwwRkVbFzHZHUk9dRiIiAighiIhIQAlBREQAJQQREQkoIYiICKCEICIiASUEEREBlBBERCSghCAiIkAru1NZpLnMW76nzvLPTclu4UhEoieiFoKZTTOzLWaWb2YP1LF8qpmtMrMqM5sRVn6lmeWFvSrM7OZg2TNmtjNs2bim2y0REWmsBlsIZhYPPAZcDRQCuWa20N03hlXbA9wFfCt8XXd/BxgXbCcNyAf+HFbl2+7+8vnsgIiINI1IuowmA/nuvgPAzOYD04G/JQR33xUsqznLdmYAb7j7iXOOVkREmk0kXUaZQEHYfGFQ1lizgBdqlf3QzNaa2VwzSzqHbYqISBOJJCFYHWXemDcxsz7AGGBxWPGDwHBgEpAGfLeedeeY2QozW1FSUtKYtxURkUaIpMuoEOgXNp8FFDXyfT4LvObulWcK3H1fMHnKzH5NrfMPYfWeAJ4AyMnJaVQiEqnr6iFdOSRSt0haCLnAUDMbaGaJhLp+FjbyfW6nVndR0GrAzAy4GVjfyG2KiEgTajAhuHsVcB+h7p5NwAJ332Bmj5jZTQBmNsnMCoGZwONmtuHM+mY2gFALY2mtTf/WzNYB64AewKPnvzsiInKuIroxzd0XAYtqlT0UNp1LqCuprnV3UcdJaHe/qjGBiohI89LQFSIiAmjoCpFG04lqaavUQhAREUAJQUREAkoIIiICKCGIiEhACUFERAAlBBERCSghiIgIoIQgIiIBJQQREQGUEEREJKCEICIigMYyklZIYwmJNA+1EEREBFBCEBGRgBKCiIgASggiIhJQQhAREUAJQUREAhElBDObZmZbzCzfzB6oY/lUM1tlZlVmNqPWsmozywteC8PKB5rZcjPbZmYvmlni+e+OiIicqwYTgpnFA48B1wIjgdvNbGStanuAu4B5dWzipLuPC143hZX/GJjr7kOBUuDuc4hfRESaSCQthMlAvrvvcPfTwHxgengFd9/l7muBmkje1MwMuAp4OSh6Frg54qhFRKTJRZIQMoGCsPnCoCxSyWa2wsw+MrMzX/rpwBF3rzrHbYqISBOLZOgKq6PMG/Ee2e5eZGaDgCVmtg44Guk2zWwOMAcgO1vDE4iINJdIWgiFQL+w+SygKNI3cPei4O8O4F1gPHAQ6GZmZxJSvdt09yfcPcfdczIyMiJ9WxERaaRIEkIuMDS4KigRmAUsbGAdAMysu5klBdM9gE8AG93dgXeAM1ckzQZeb2zwIiLSdBpMCEE//33AYmATsMDdN5jZI2Z2E4CZTTKzQmAm8LiZbQhWHwGsMLM1hBLAj9x9Y7Dsu8D9ZpZP6JzCU025YyIi0jgRDX/t7ouARbXKHgqbziXU7VN7vQ+BMfVscwehK5hERCQG6E5lEREBlBBERCSghCAiIoAeoSkxoq7HYoIejSnSktRCEBERQAlBREQC6jISaUZ1dYWpG0xilVoIIiICKCGIiEhACUFERAAlBBERCSghiIgIoIQgIiIBJQQREQGUEEREJKCEICIigBKCiIgElBBERARQQhARkYASgoiIABEmBDObZmZbzCzfzB6oY/lUM1tlZlVmNiOsfJyZLTOzDWa21sxuC1v2jJntNLO84DWuaXZJRETORYPDX5tZPPAYcDVQCOSa2UJ33xhWbQ9wF/CtWqufAO50921m1hdYaWaL3f1IsPzb7v7y+e6EiIicv0iehzAZyHf3HQBmNh+YDvwtIbj7rmBZTfiK7r41bLrIzIqBDOAIIiISUyLpMsoECsLmC4OyRjGzyUAisD2s+IdBV9JcM0tq7DZFRKTpRJIQrI4yb8ybmFkf4DfAP7j7mVbEg8BwYBKQBny3nnXnmNkKM1tRUlLSmLcVEZFGiCQhFAL9wuazgKJI38DMUoA/Av/i7h+dKXf3fR5yCvg1oa6pj3H3J9w9x91zMjIyIn1bERFppEgSQi4w1MwGmlkiMAtYGMnGg/qvAc+5+0u1lvUJ/hpwM7C+MYGLiEjTajAhuHsVcB+wGNgELHD3DWb2iJndBGBmk8ysEJgJPG5mG4LVPwtMBe6q4/LS35rZOmAd0AN4tEn3TEREGiWSq4xw90XAolplD4VN5xLqSqq93vPA8/Vs86pGRSoiIs1KdyqLiAighCAiIgElBBERAZQQREQkoIQgIiKAEoKIiASUEEREBIjwPgSRczVv+Z6PlX1uSnYUIhGRhqiFICIigBKCiIgE1GUkEgPq6loDda9Jy1ILQUREACUEEREJKCGIiAighCAiIgElBBERAZQQREQkoIQgIiKAEoKIiASUEEREBIgwIZjZNDPbYmb5ZvZAHcunmtkqM6sysxm1ls02s23Ba3ZY+UQzWxds8+dmZue/OyIicq4aTAhmFg88BlwLjARuN7ORtartAe4C5tVaNw34PjAFmAx838y6B4t/CcwBhgavaee8FyIict4iaSFMBvLdfYe7nwbmA9PDK7j7LndfC9TUWvfTwJvuftjdS4E3gWlm1gdIcfdl7u7Ac8DN57szIiJy7iJJCJlAQdh8YVAWifrWzQymz2WbIiLSDCJJCHX17XuE269v3Yi3aWZzzGyFma0oKSmJ8G1FRKSxIkkIhUC/sPksoCjC7de3bmEw3eA23f0Jd89x95yMjIwI31ZERBorkoSQCww1s4FmlgjMAhZGuP3FwDVm1j04mXwNsNjd9wHHzOyi4OqiO4HXzyF+ERFpIg0mBHevAu4j9OW+CVjg7hvM7BEzuwnAzCaZWSEwE3jczDYE6x4G/g+hpJILPBKUAXwZeBLIB7YDbzTpnomISKNE9MQ0d18ELKpV9lDYdC7/uwsovN7TwNN1lK8ARjcmWBERaT66U1lERAAlBBERCSghiIgIoIQgIiIBJQQREQGUEEREJBDRZacibVHRkZPk7jrM4fLTfJB/iH5pHRnUowvxcRqJXdonJQRpVyoqq1mQW8CrqwtZvvMwXmsErY4d4hnXrxtXDu9JlyR9PKR90X+8tHnVNc624mOs3nOETfuOUlXjDEjvxD99chjXjOpFr5RkXs/by/bictYXlbF85yFWF5Ry1fBefDYni4R49axK+6CEIG3WgaMV/Gn9PlbuLqX8dDWdEuPJGdCd70wbzvh+3Qh/SF9SQjwj+6Ywsm8Klw/LYNG6fSxat4+9pSf4z5kXMrRX1yjuiUjLUEKQNqf4aAVz39rKKyv3Ulldw4g+KUzI7s6w3l1IiItjQnb3s67fKyWZuy4ZwLq9Zfx54wGu//n7fO/6Edx5cX/0pFdpy5QQpM2oqq7h2WW7mfvmVk5X1fDZSVn06ppMepekRm/LzBib1Y1vffoCvvvyWr6/cAPbio/x/RtHNUPkIrFBCUHahA+3H+ThhRvYeuA4lw/L4OGbRjGwR2fmLd9zXtvt0SWJX92Zw48Xb+bxpTvYc/gknxreU+cVpE1SQpBW7ciJ0yxav5/1e8vI6t6Rx78wkWtG9mrSrp24OOPBa0cwIL0zD766jmMnK5kxMStq3Ud1JbnPTcmOQiTS1ighSKtUVVPDX7YeZOnWYtzhG58axj2XDyK5Q3yzveftk7MpOXaKn765lfQuSVw1vGezvZdINCghSKtTfKyCBbkFFJVVMKpvCteN6cNXrhzSIu/91auGsGRzMW9tOkDvlGRG9k1pkfcVaQlKCNKqvLa6kMfeyadDfByfn9K/xb+QzYxbxmdSfLSC11YXkp0+TDewSZuhM2PSaizILeD+BWvo170TX/vk0Kj9Ok+Ij2NGTj8qqmp4PW8vXvt2Z5FWSglBWoUFuQV899W1TB2awexLBpCS3CGq8fROSebqEb3YUHSUNYVHohqLSFNRQpCY99bGA3z31bVcNjSDx78wkQ4xcsnnpUN7kJ3WiYVrithfVhHtcETOW0SfLDObZmZbzCzfzB6oY3mSmb0YLF9uZgOC8jvMLC/sVWNm44Jl7wbbPLNMl2zIx2wsOsrX5q9mTGYqj39+YrNeRdRYcWbMnJhFdY3znVfWqutIWr0GE4KZxQOPAdcCI4HbzWxkrWp3A6XuPgSYC/wYwN1/6+7j3H0c8AVgl7vnha13x5nl7l7cBPsjbcixikq++GwuKckdePLOHDomxk4yOCO9SxLXju7DX7aWMO+v53cTnEi0RXJ5xGQg3913AJjZfGA6sDGsznTg4WD6ZeAXZmb+v38y3Q68cN4RS9S1xI1RldU1PP/RbkpPVPLSvRfTMyW5SbfflKYMTKP0xGl++MdNXDqkB/3TO0c7JJFzEkmXUSZQEDZfGJTVWcfdq4AyIL1Wndv4eEL4ddBd9K9Wz22fZjbHzFaY2YqSkpIIwpXWzt15ZVUhBaUnmXvbhYzOTI12SGdlZvz41rHExxnfemkN1TXqOpLWKZKEUNcXde3/+LPWMbMpwAl3Xx+2/A53HwNcFry+UNebu/sT7p7j7jkZGRkRhCut3ZItxawtLOOakb2YNrpPtMOJSN9uHXn4xlHk7irl6fd3RjsckXMSSUIoBPqFzWcBRfXVMbMEIBU4HLZ8FrVaB+6+N/h7DJhHqGtK2rm8giO8vamY8f26cfmw1vUD4JYJmVwzshc/+fMWth44Fu1wRBotkoSQCww1s4Fmlkjoy31hrToLgdnB9AxgyZnzB2YWB8wE5p+pbGYJZtYjmO4A3ACsR9q1nQfLeWVVIQPSO/OZ8Zmt7tkDZsa/3TKGLkkJfHPBGiqra6IdkkijNJgQgnMC9wGLgU3AAnffYGaPmNlNQbWngHQzywfuB8IvTZ0KFJ45KR1IAhab2VogD9gL/Oq890ZarYPHTvH8R7vp3imRz1+U3WqHl+7RJYl/+8xo1u0t47F38qMdjkijRDQIi7svAhbVKnsobLqCUCugrnXfBS6qVVYOTGxkrNJGHTp+imeW7cIMZl/cn06JrXtsoGmj+/CZ8Zn8Ykk+90wdTGb3jtEOSSQirfNnmLQZFZXVzPnNSo6erOTOi/qf09PNYtHDN46iR5ckXlpZQJW6jqSVUEKQqHF3Hnx1HSt3lzIzpx/Zbej6/dROHfjRrWMoPnaKd7bonktpHZQQJGp+u3wPr63ey/1XD2NMjN9rcC6uuKAnE7K7sXRrCUVHTkY7HJEGKSFIVKzfW8Yjv9/I5cMyuK+FHm4TDdeN6UOnxAReXVWoG9Yk5ikhSIs7WlHJP/52FeldEpl72zji4lrX5aWN0SkxgZsu7EtRWQXvbdOd9hLblBCkRbk7//zqOvYeOckvPjeetM6J0Q6p2Y3OTGV0Zipvby6m+KiGyZbYpYQgLWrVnlL+sHYf9189jIn906IdTou5cWwfEuPjeEVdRxLDlBCkxRQfq2DhmiIuGZzOvZcPjnY4LaprcgduGNuHgtKTPPPhrmiHI1InJQRpEZXVNbyYW0CH+Djm3jaO+DZ83qA+4/p144JeXfnJ4s3sPlQe7XBEPkYJQVrE4g372VdWwYwJWfSK4WcbNCcz4+bxmXSIi+OBV9ZRo64jiTFKCNLsNu87yofbD3Hx4HSG90mJdjhRldqxAw9eN4JlOw7xQq6esCaxRQlBmtWJ01W8snovfVKTmTaqd7TDiQm3T+7HJYPT+fdFm9mrG9YkhighSLP60/r9nDxdxYyJWXRopSOYNrUzT1ircefbL61R15HEjNY9rKTEtGXbD7FidylTh2bQJ1Ujfobrl9aJf7l+JP/82jp+89HuZkuWLfH8a2k79JNNmkVFZTXfe20daZ0TuWp4z2iHE5Nun9yPKy7I4N/f2MTBY6eiHY6IEoI0j+c/2s2Og+VMv7AviQn6N6vLma6jpIR4Xl29l+AhgyJRo0+qNLljFZU89k4+lw3twdBeXaMdTkzrlZLMg9cOZ9ehclbtKY12ONLOKSFIk3vq/Z2Unqjk25++INqhtAqfzelH//ROvLF+P+WnqqIdjrRjSgjSpA6Xn+bJ93YybVRvxmZ1i3Y4rUJcnHHzuEwqKqt5Y/3+aIcj7ZgSgjSpx5du58TpKr55zbBoh9Kq9EpJ5rKhGazaU8qOkuPRDkfaqYgSgplNM7MtZpZvZg/UsTzJzF4Mli83swFB+QAzO2lmecHrf8LWmWhm64J1fm5m7W9wmzbmyInTPP/Rbm68sK/OHZyDKy/oSVrnRH6XV6TnMEtUNJgQzCweeAy4FhgJ3G5mI2tVuxsodfchwFzgx2HLtrv7uOB1b1j5L4E5wNDgNe3cd0NiwXPLdlN+upovX9G+RjJtKokJcdx0YV8OHj/FUj1MR6IgkhbCZCDf3Xe4+2lgPjC9Vp3pwLPB9MvAJ8/2i9/M+gAp7r7MQ9faPQfc3OjoJWacOF3Frz/YyVXDezK8d/ser+h8DOvVlTGZqSzdUqJ7E6TFRZIQMoGCsPnCoKzOOu5eBZQB6cGygWa22syWmtllYfULG9gmAGY2x8xWmNmKkhL9aopVL+YWUHqikn9U6+C8XT+2Dwnxxu/ydG+CtKxIEkJdv/Rr/5fWV2cfkO3u44H7gXlmlhLhNkOF7k+4e46752RkZEQQrrS06hrnV3/ZwaQB3ckZ0H6egtZcUpI78OlRvdlxsJxXVu2NdjjSjkSSEAqBfmHzWUBRfXXMLAFIBQ67+yl3PwTg7iuB7cCwoH5WA9uUVmJNwRGKyir4xyuGRDuUNmPSgDT6p3Xi0T9u5NBxdR1Jy4gkIeQCQ81soJklArOAhbXqLARmB9MzgCXu7maWEZyUxswGETp5vMPd9wHHzOyi4FzDncDrTbA/0sJq3Fm6rYThvbtyxQVqwTWVuOBhOuWnqnjkDxujHY60Ew0mhOCcwH3AYmATsMDdN5jZI2Z2U1DtKSDdzPIJdQ2duTR1KrDWzNYQOtl8r7sfDpZ9GXgSyCfUcnijifZJWtDmfUcpOXaKL18xGF053LR6pSTzlSuH8HpeEQvXqAEtzS+i4a/dfRGwqFbZQ2HTFcDMOtZ7BXilnm2uAEY3JliJLe7Ou1tLSOucyPVj+kQ7nDbpviuHsHRrCf/y2jpy+nenbzcNIy7NR3cqyznbcbCcwtKTXDa0Bwl6+E2zSIiP42e3jaOqxvnmAj1MR5qXHpAjQN0PUoH6H6bi7ry9qZiuyQlMyO7enKG1e/3TO/PwjaP4zitrefL9HXRJ6hDtkKSN0s86OSc7Dpaz61A5lw/L0KMxW8DMnCymjerNTxZvoUjPYZZmok+yNJq789bGA6QkJzBJ9x20CDPj328ZQ/dOiSxYUUClxjqSZqCEII2WX3yc3YdPcMUFPdU6aEHdOyfynzMvpPjYKRat2xftcKQN0qdZGsXdeWvTAVI7diCnv84dtLSpwzK4dEgPlu88zLq9ZdEOR9oYnVSWRnlj/X4KSk9yy/hMXVkUJZ8e1Zvdh8p5dVUhfVOTSe+S1Kj167qAoL6LB6R90SdaIna6qob/+NNmenZNYoJaB1ETH2fMmpxNnBkv5O7RsxOkySghSMTmLd/NrkMnuHZ0b+J0V3JUde+UyK0Tsig6UqHHbkqTUUKQiBytqOTnS/K5ZHA6w/Q0tJgwsm8KnxiczrIdh3hDJ5mlCSghSEQeX7qdw+WnefDaERqzKIZ8enRvsrp35DuvrGXXwfJohyOtnBKCNGhf2UmefG8n08f1ZUxWarTDkTAJcXHMmpRNfJxx97O5lJ2sjHZI0oopIUiDfvrnrbjDt665INqhSB3SOifyP5+fyJ7DJ7hv3iqdZJZzpoQgZ7V5/1FeXlXI7Ev60y+tU7TDkXpcNCidH948hve2HeQHv9fzE+Tc6D4EqZe78+gfNtE1KYGvXKmnocW6z07qx/aS4zz+lx0M6dmF2ZcMiHZI0sooIUi9NhQd5f38gzx840i6dUqMdjgSge9MG86Og+X84Pcb6J+uFp00jrqMpE6nq2pYtG4fw3t35fMX9Y92OBKh+DjjZ7eNY3jvFL46bzX7j1ZEOyRpRZQQpE5Lt5Zw5GQlP7hplIaoaGU6JyXw1F05dEyM57kPd3GsQlceSWT0SZePKT5WwXvbShiblcqUQenRDkfOQZ/Ujjx91yTKT1fx3LLdnK7SlUfSsIgSgplNM7MtZpZvZg/UsTzJzF4Mli83swFB+dVmttLM1gV/rwpb591gm3nBq2dT7ZScuxp3Xl21lw7xcXpOcis3OjOVWZOyKTpykgUrCqhxPX5Tzq7Bk8pmFg88BlwNFAK5ZrbQ3cOvbbsbKHX3IWY2C/gxcBtwELjR3YvMbDSwGMgMW+8Od1/RRPsiTeDD7YfYc/gEMydm0TVZj2ps7Ub0SeH6sX34w9p9/Gn9fq5rgiSv0VLbrkhaCJOBfHff4e6ngfnA9Fp1pgPPBtMvA580M3P31e5eFJRvAJLNrHFj9UqLOXT8FG9u3M/w3l0Z169btMORJnLJ4B5cPCid9/MP8tGOQ9EOR2JYJJedZgIFYfOFwJT66rh7lZmVAemEWghn3AqsdvdTYWW/NrNq4BXgUXe1aZtapL/mqmpqeHFFAfFxxvRxmRqvqI25fmwfSk+c5vdriuiclMCYTA1BIh8XSQuhrm+G2l/cZ61jZqMIdSPdE7b8DncfA1wWvL5Q55ubzTGzFWa2oqSkJIJw5Vy8ueEAhaUnuXVCFqkd1VXU1sSZMWtSNtlpnViQW8C2A8eiHZLEoEgSQiHQL2w+Cyiqr46ZJQCpwOFgPgt4DbjT3befWcHd9wZ/jwHzCHVNfYy7P+HuOe6ek5GREck+SSO9s6WY9/IPMmVgGqP66pdjW5WYEMedFw+gZ0oSzy/fzY6Dx6MdksSYSBJCLjDUzAaaWSIwC1hYq85CYHYwPQNY4u5uZt2APwIPuvsHZyqbWYKZ9QimOwA3AOvPb1fkXBQcPsH9L+bROyW5SU44SmzrmBjPXZcMoFunRJ75YBeb9x+NdkgSQxpMCO5eBdxH6AqhTcACd99gZo+Y2U1BtaeAdDPLB+4Hzlyaeh8wBPjXWpeXJgGLzWwtkAfsBX7VlDsmDTt5upp7frOS6hrnjinZdNANaO1C1+QOzLlsEL1Sknn+o93kFRyJdkgSIyIay8jdFwGLapU9FDZdAcysY71HgUfr2ezEyMOUpubuPPjqWjbtP8rTsyexr0xDHLQnnZMSuPvSgfzmo90sWFFA75Qk/ulTw4iL08UE7Zl+ErZT/+/d7fwur4hvfGoYVw7XPYHtUXKHeP7hkgFMzO7Oz5fk85V5qzTMRTunhNAOvZ63l58s3sL0cX356lUa1ro9S4iP45YJmXzvuhEs3rCfG//7fdYVlkU7LIkSJYR2ZkfJcb790lqmDEzjP2aM1f0GgpnxpamDmD/nYioqa7jllx/w9Ps70W1B7Y8SQjuy62A5zy3bTXZ6J574Qg5JCfHRDkliyOSBabzx9cuYOjSDR/6wkS89t5LS8tPRDktakBJCO7H7UDnPLNtFSscOzPviFFI76eYz+bjunRN5cnYO/3rDSJZuLeba/3qPpVt1Q2h7oSemtUKNHVxs076jvJhbQErHBL542UB6piQ3Z3jSypkZd186kMkD0vjGgjxmP/1Xbp+czfeuH0GXJH1ltGVqIbRh7s7720p4/qPdZHRN4ouXDSJFI5hKhMZkpfKHr17KPVMHMT93D9N+9heWbdfgeG2ZEkIbVXaykq/Nz2PR+v2M7JvCl5QM5Bwkd4jnwetG8NI9F5MQZ9z+q49YtG4fVTV64E5bpPZfG/TXnYf5xot57D9awadG9OKKCzKI09VEch5yBqSx6OuX8aM3NvPcst3sOlTOrEnZpHVOjHZo0oSUENqQshOV/OhPm3nhr3vITuvES/dezOZ9GtVSmkanxAQemT6aqmrn1dWF/OKdbdw6IeusAyLqYTqtixJCG1BVXcOzH+7iv5ds43D5ab546UC+cfUwOiclKCFIkxudmUrfbh2Zn7uH3y7fw8WD0rl1YqYuY24DlBBaseoaZ9XuUpZsKabsZCWTB6bxzD+MZLQefiLNLK1zInOmDmLx+v18sP0QM365jF98bjz90ztHOzQ5D0oIrVB1jbOm4AhLthRzuPw0/bp35LHPTeATQ9J157G0mIS4OK4f25eBPTqzcE0RN/z8ff7tljHcMLaP/g9bKV1l1IrU1DgL1xTxX29v4+VVhSQnxDH74v7ce/lgLh3aQx9CiYqRfVP549cuY3DPLnz1hdV88dkVFBw+Ee2w5ByohdAKuDuLNxxg7ptb2XLgGL1SkrhjSjYj+6QoCUhM6JfWiZfvvZhff7CLuW9t5eq5S7nz4gGkd06kawSXO9d18hl0ArqlKSHEMHfnrU3F/OytrWwoOsqgjM78/PbxHD1ZqctIJeYkxMfxpamDuG5sH37yp808+d4O4uOMCdndmTQgjb7dOkY7RGmAEkKMCP+F5O5s3n+MtzcdoKisgv7pnfjPmRdy87i+JMTH1ftrSiQWZHbryM9mjedrnxzKt15ay8rdpSzfeZg+qcmM7JPCmMxURmeqdRuLlBBiSHWNs6GojL9sLaGorIK0zon8ZMZYPjM+kwQ93lJamUEZXZgxMYvrx/Qhr6CUvIIjLNlczNubi+mVksRVw3ty1fBeXDI4PdqhSkAJIQbkFx/nzY0HWLH7MMcqqkjrnMitE7IY168bM3P6RTs8kfPSMTGeiwf34OLBPTh+KvT//famAyzMK+KFvxaQGB9HdnonhvXqygW9utKjS6JaD1GihBAF1TVsSO0KAAAIB0lEQVTOmsIjvLulhDc3HmDTvqMYMLRXFz4zPp1hvbrqHIG0SV2SEpgxMYsZE7M4VVXNil2lvLulmN/lFbFo3T4WrdtHWudERvTuyog+KVRV16h13IIiSghmNg34LyAeeNLdf1RreRLwHDAROATc5u67gmUPAncD1cDX3H1xJNtsK+Yt38PxU1UUHj7BntITFBw+wYGjpzh+qoo4gwnZ3fn+jSM5VVWjweekXUlKiOcTQ3rwiSE9GNijC6Xlp9ly4Bib9x/lo52H+WD7IV5aWcgVF2Rw2dAMpgxMo19ap2iH3aY1mBDMLB54DLgaKARyzWyhu28Mq3Y3UOruQ8xsFvBj4DYzGwnMAkYBfYG3zGxYsE5D22xV3J2S46fILz7O9uLj5BcfJ7/kOOsKyzhaUQVAnEHv1GQ+Mz6TSQPTmDq0B906hQYH04liae+6d07kokHpXDQonVNV1Ww7cJzT1TW8s7mY1/OKAOidkky3Th3ok5pMWuckunfqQLdOidxz+SA6qCVx3iJpIUwG8t19B4CZzQemA+Ff3tOBh4Ppl4FfWKgTcDow391PATvNLD/YHhFss0W4OzUONe7UuON/mw79ra52jp+qovx0FeWnqjhWUUXZyUr2lVWw78hJ9pVVsP9oBbsOlv/tix9CTePBGZ0ZnNGFXinJ9EvrRGa3jiQmxOnaapEGJCXEMzozlc9NyaamxtlWfJxl2w+SV3CED7cfYsv+Y4Q/8fmnb26hZ9dQskjp2IHU4NUlKYHkDvEkJcSd9W98nJEQb8THGfFmf58/Mx0XR1xc6O7s+DjDDOLMMMAs9FAhMzCCcgMj+BtMx52pF7ZOrIkkIWQCBWHzhcCU+uq4e5WZlQHpQflHtdbNDKYb2maTufc3K1m6taTWF37oS/98dE1OoG9qR3qnJjN2XCpDMrowpGdXhvTsQq+UJMxMv/xFzlNcnHFB765c0LsrEGpNV9XUUHaiktITlRw5cZqs7h3ZV1ZB2clKjpyspODwCdafrORw+Wmqapzq8/2wN5M6E0swDaFEcqbe7796KYMzujRvPO5nP1BmNhP4tLt/MZj/AjDZ3b8aVmdDUKcwmN9OqCXwCLDM3Z8Pyp8CFhEaMuOs2wzb9hxgTjB7AbClkfvYAzjYyHWiobXECYq1ObSWOEGxNofmjrO/u2c0VCmSFkIhEH7tYxZQVE+dQjNLAFKBww2s29A2AXD3J4AnIoizTma2wt1zznX9ltJa4gTF2hxaS5ygWJtDrMQZyVmYXGComQ00s0RCJ4kX1qqzEJgdTM8Alnio6bEQmGVmSWY2EBgK/DXCbYqISAtqsIUQnBO4D1hM6BLRp919g5k9Aqxw94XAU8BvgpPGhwl9wRPUW0DoZHEV8BV3rwaoa5tNv3siIhKpiO5DcPdFhPr+w8seCpuuAGbWs+4PgR9Gss1mcs7dTS2stcQJirU5tJY4QbE2h5iIs8GTyiIi0j7oTg4REQHacEIws11mts7M8sxsRbTjCWdmT5tZsZmtDytLM7M3zWxb8Ld7NGM8o55YHzazvcGxzTOz66IZYxBTPzN7x8w2mdkGM/t6UB5zx/UsscbUcTWzZDP7q5mtCeL8QVA+0MyWB8f0xeDCkKg6S6zPmNnOsGM6LtqxQmgECDNbbWZ/COZj4pi22YQQuNLdx8XC5Vy1PANMq1X2APC2uw8F3g7mY8EzfDxWgLnBsR0XnA+Ktirgm+4+ArgI+EowdEosHtf6YoXYOq6ngKvc/UJgHDDNzC4iNDTN3OCYlhIauiba6osV4NthxzQveiH+L18HNoXNx8QxbesJISa5+18IXY0VbjrwbDD9LHBziwZVj3pijTnuvs/dVwXTxwh92DKJweN6llhjioccD2Y7BC8HriI0RA3EzjGtL9aYY2ZZwPXAk8G8ESPHtC0nBAf+bGYrg7udY10vd98HoS8MoGeU42nIfWa2NuhSino3TDgzGwCMB5YT48e1VqwQY8c16NrIA4qBN4HtwBF3PzNwV/hwNFFVO1Z3P3NMfxgc07kWGpk52n4GfAeoCebTiZFj2pYTwifcfQJwLaEm+dRoB9SG/BIYTKhpvg/4v9EN5+/MrAvwCvBP7n402vGcTR2xxtxxdfdqdx9HaDSBycCIuqq1bFR1qx2rmY0GHgSGA5OANOC7UQwRM7sBKHb3leHFdVSNyjFtswnB3YuCv8XAa/x9lNVYdcDM+gAEf4ujHE+93P1A8OGrAX5FjBxbM+tA6Av2t+7+alAck8e1rlhj9bgCuPsR4F1C5zy6BUPUwFmGnYmWsFinBd1zHoy4/Guif0w/AdxkZruA+YS6in5GjBzTNpkQzKyzmXU9Mw1cA6w/+1pRFz78x2zg9SjGclZnvmADnyEGjm3QD/sUsMndfxq2KOaOa32xxtpxNbMMM+sWTHcEPkXofMc7hIaogdg5pnXFujnsx4AR6peP6jF19wfdPcvdBxAa0WGJu99BjBzTNnljmpkNItQqgNDd2POCO6Zjgpm9AFxBaITDA8D3gd8BC4BsYA8w092jfjK3nlivINSt4cAu4J4z/fTRYmaXAu8B6/h73+w/E+qbj6njepZYbyeGjquZjSV0gjOe0I/HBe7+SPD5mk+oC2Y18PngF3jUnCXWJUAGoW6ZPODesJPPUWVmVwDfcvcbYuWYtsmEICIijdcmu4xERKTxlBBERARQQhARkYASgoiIAEoIIiISUEIQERFACUFERAJKCCIiAsD/B1Ej2tpnT+3dAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fb9be40d240>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"sbn.distplot(results)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"# Do diff cut-offs analysis, whilst being blind to case-control status, as per coggon suggestion"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" cases controls odds ratio p-value\n",
"ever 268 110 1.08 0.88\n",
"never 43 19 1.00 0.88\n"
]
},
{
"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>cases</th>\n",
" <th>controls</th>\n",
" <th>odds ratio</th>\n",
" <th>p-value</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>ever</th>\n",
" <td>268</td>\n",
" <td>110</td>\n",
" <td>1.08</td>\n",
" <td>0.88</td>\n",
" </tr>\n",
" <tr>\n",
" <th>never</th>\n",
" <td>43</td>\n",
" <td>19</td>\n",
" <td>1.00</td>\n",
" <td>0.88</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" cases controls odds ratio p-value\n",
"ever 268 110 1.08 0.88\n",
"never 43 19 1.00 0.88"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ever_vs_never(df)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" cases controls odds ratio p-value\n",
"ever 74 64 2.51 0.09\n",
"never 6 13 1.00 0.09\n"
]
},
{
"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>cases</th>\n",
" <th>controls</th>\n",
" <th>odds ratio</th>\n",
" <th>p-value</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>ever</th>\n",
" <td>74</td>\n",
" <td>64</td>\n",
" <td>2.51</td>\n",
" <td>0.09</td>\n",
" </tr>\n",
" <tr>\n",
" <th>never</th>\n",
" <td>6</td>\n",
" <td>13</td>\n",
" <td>1.00</td>\n",
" <td>0.09</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" cases controls odds ratio p-value\n",
"ever 74 64 2.51 0.09\n",
"never 6 13 1.00 0.09"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ever_vs_never(df[df['distfromcentre'] < 10].copy())"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" cases controls odds ratio p-value\n",
"Non-construction high risk occupations 93 48 0.86 0.75\n",
"Construction 61 21 1.28 0.57\n",
"Medium risk industrial 59 16 1.63 0.24\n",
"Low risk industrial 55 25 0.97 1.00\n",
"Office and other low risk 43 19 1.00 1.00\n",
"\n",
"\n",
" cases controls odds ratio \\\n",
"Non-construction high risk occupations 93 48 0.86 \n",
"Carpenters 5 4 0.55 \n",
"Plumber, electrician, painter or decorator 41 13 1.39 \n",
"Other construction 15 4 1.66 \n",
"Medium risk industrial 59 16 1.63 \n",
"Low risk industrial 55 25 0.97 \n",
"Office and other low risk 43 19 1.00 \n",
"\n",
" p-value \n",
"Non-construction high risk occupations 0.75 \n",
"Carpenters 0.46 \n",
"Plumber, electrician, painter or decorator 0.53 \n",
"Other construction 0.56 \n",
"Medium risk industrial 0.24 \n",
"Low risk industrial 1.00 \n",
"Office and other low risk 1.00 \n"
]
},
{
"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>cases</th>\n",
" <th>controls</th>\n",
" <th>odds ratio</th>\n",
" <th>p-value</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>Non-construction high risk occupations</th>\n",
" <td>93</td>\n",
" <td>48</td>\n",
" <td>0.86</td>\n",
" <td>0.75</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Carpenters</th>\n",
" <td>5</td>\n",
" <td>4</td>\n",
" <td>0.55</td>\n",
" <td>0.46</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Plumber, electrician, painter or decorator</th>\n",
" <td>41</td>\n",
" <td>13</td>\n",
" <td>1.39</td>\n",
" <td>0.53</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Other construction</th>\n",
" <td>15</td>\n",
" <td>4</td>\n",
" <td>1.66</td>\n",
" <td>0.56</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Medium risk industrial</th>\n",
" <td>59</td>\n",
" <td>16</td>\n",
" <td>1.63</td>\n",
" <td>0.24</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Low risk industrial</th>\n",
" <td>55</td>\n",
" <td>25</td>\n",
" <td>0.97</td>\n",
" <td>1.00</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Office and other low risk</th>\n",
" <td>43</td>\n",
" <td>19</td>\n",
" <td>1.00</td>\n",
" <td>1.00</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Construction</th>\n",
" <td>61</td>\n",
" <td>21</td>\n",
" <td>1.28</td>\n",
" <td>0.57</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" cases controls odds ratio \\\n",
"Non-construction high risk occupations 93 48 0.86 \n",
"Carpenters 5 4 0.55 \n",
"Plumber, electrician, painter or decorator 41 13 1.39 \n",
"Other construction 15 4 1.66 \n",
"Medium risk industrial 59 16 1.63 \n",
"Low risk industrial 55 25 0.97 \n",
"Office and other low risk 43 19 1.00 \n",
"Construction 61 21 1.28 \n",
"\n",
" p-value \n",
"Non-construction high risk occupations 0.75 \n",
"Carpenters 0.46 \n",
"Plumber, electrician, painter or decorator 0.53 \n",
"Other construction 0.56 \n",
"Medium risk industrial 0.24 \n",
"Low risk industrial 1.00 \n",
"Office and other low risk 1.00 \n",
"Construction 0.57 "
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bycat(df)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" cases controls odds ratio p-value\n",
"Non-construction high risk occupations 27 28 2.09 0.28\n",
"Construction 18 16 2.44 0.16\n",
"Medium risk industrial 18 8 4.88 0.02\n",
"Low risk industrial 11 12 1.99 0.35\n",
"Office and other low risk 6 13 1.00 1.00\n",
"\n",
"\n",
" cases controls odds ratio \\\n",
"Non-construction high risk occupations 27 28 2.09 \n",
"Carpenters 2 4 1.08 \n",
"Plumber, electrician, painter or decorator 12 8 3.25 \n",
"Other construction 4 4 2.17 \n",
"Medium risk industrial 18 8 4.88 \n",
"Low risk industrial 11 12 1.99 \n",
"Office and other low risk 6 13 1.00 \n",
"\n",
" p-value \n",
"Non-construction high risk occupations 0.28 \n",
"Carpenters 1.00 \n",
"Plumber, electrician, painter or decorator 0.11 \n",
"Other construction 0.41 \n",
"Medium risk industrial 0.02 \n",
"Low risk industrial 0.35 \n",
"Office and other low risk 1.00 \n"
]
},
{
"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>cases</th>\n",
" <th>controls</th>\n",
" <th>odds ratio</th>\n",
" <th>p-value</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>Non-construction high risk occupations</th>\n",
" <td>27</td>\n",
" <td>28</td>\n",
" <td>2.09</td>\n",
" <td>0.28</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Carpenters</th>\n",
" <td>2</td>\n",
" <td>4</td>\n",
" <td>1.08</td>\n",
" <td>1.00</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Plumber, electrician, painter or decorator</th>\n",
" <td>12</td>\n",
" <td>8</td>\n",
" <td>3.25</td>\n",
" <td>0.11</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Other construction</th>\n",
" <td>4</td>\n",
" <td>4</td>\n",
" <td>2.17</td>\n",
" <td>0.41</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Medium risk industrial</th>\n",
" <td>18</td>\n",
" <td>8</td>\n",
" <td>4.88</td>\n",
" <td>0.02</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Low risk industrial</th>\n",
" <td>11</td>\n",
" <td>12</td>\n",
" <td>1.99</td>\n",
" <td>0.35</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Office and other low risk</th>\n",
" <td>6</td>\n",
" <td>13</td>\n",
" <td>1.00</td>\n",
" <td>1.00</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Construction</th>\n",
" <td>18</td>\n",
" <td>16</td>\n",
" <td>2.44</td>\n",
" <td>0.16</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" cases controls odds ratio \\\n",
"Non-construction high risk occupations 27 28 2.09 \n",
"Carpenters 2 4 1.08 \n",
"Plumber, electrician, painter or decorator 12 8 3.25 \n",
"Other construction 4 4 2.17 \n",
"Medium risk industrial 18 8 4.88 \n",
"Low risk industrial 11 12 1.99 \n",
"Office and other low risk 6 13 1.00 \n",
"Construction 18 16 2.44 \n",
"\n",
" p-value \n",
"Non-construction high risk occupations 0.28 \n",
"Carpenters 1.00 \n",
"Plumber, electrician, painter or decorator 0.11 \n",
"Other construction 0.41 \n",
"Medium risk industrial 0.02 \n",
"Low risk industrial 0.35 \n",
"Office and other low risk 1.00 \n",
"Construction 0.16 "
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bycat(df[df['distfromcentre'] < 10].copy())"
]
}
],
"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.5.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment