PKU Model.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Estimating the probability of IQ impairment from blood phenylalanine for phenylketonuria patients"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The association of blood phenylalanine levels with IQ was meta-analyzed using a hierarchical mixed-effects model, estimated using Markov chain Monte Carlo methods. In an effort to partially pool the information from the set of studies obtained in the literature search, we specified random effects for the intercept and slope parameters of a linear relationship between blood Phe level and IQ. Importantly, this allowed each study to have its own parameters, each sampled from a conceptual population of parameters. Those with smaller sample sizes were automatically shrunk towards the population means for each parameter, with larger studies influencing the estimate of the population mean more than being influenced by it. In turn, the magnitude of the effect (*i.e.* slope) was specified partly as a function of a fixed effect that accounted for whether measurements of Phe were carried out during the critical period. Hence, the overall model was a hierarchical mixed effects model. Bayesian hierarchical models are easily estimated using Markov chain Monte Carlo (MCMC) methods."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Data import and cleanup"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The following entries show how data were imported and cleaned prior to being input into the model."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import pandas as pd"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "observations = pd.read_csv(\"kq1-2.csv\", index_col=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "observations.head()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>Paper ID</th>\n",
        "      <th>Group ID</th>\n",
        "      <th>N</th>\n",
        "      <th>Age</th>\n",
        "      <th>Age SD</th>\n",
        "      <th>Age low</th>\n",
        "      <th>Age high</th>\n",
        "      <th>Concurrent</th>\n",
        "      <th>Critical period</th>\n",
        "      <th>Phe</th>\n",
        "      <th>Phe SD</th>\n",
        "      <th>Phe low</th>\n",
        "      <th>Phe high</th>\n",
        "      <th>IQ scale</th>\n",
        "      <th>IQ</th>\n",
        "      <th>IQ SD</th>\n",
        "      <th>IQ low</th>\n",
        "      <th>IQ high</th>\n",
        "      <th>Correlation</th>\n",
        "      <th>p</th>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>Observation ID</th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>1</th>\n",
        "      <td> 54</td>\n",
        "      <td>  1</td>\n",
        "      <td> 46</td>\n",
        "      <td>  7.500000</td>\n",
        "      <td> 3.3</td>\n",
        "      <td> 2.9</td>\n",
        "      <td> 15.5</td>\n",
        "      <td> 0</td>\n",
        "      <td> 1</td>\n",
        "      <td>  312</td>\n",
        "      <td> 132</td>\n",
        "      <td> 125</td>\n",
        "      <td> 852</td>\n",
        "      <td> Weschler</td>\n",
        "      <td> 104</td>\n",
        "      <td> 15</td>\n",
        "      <td> 68</td>\n",
        "      <td> 143</td>\n",
        "      <td>-0.17</td>\n",
        "      <td> 0.38</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2</th>\n",
        "      <td> 80</td>\n",
        "      <td> 12</td>\n",
        "      <td>  1</td>\n",
        "      <td>  7.666667</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>  NaN</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0</td>\n",
        "      <td>  704</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>    Raven</td>\n",
        "      <td>  91</td>\n",
        "      <td>NaN</td>\n",
        "      <td>NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>  NaN</td>\n",
        "      <td>  NaN</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>3</th>\n",
        "      <td> 80</td>\n",
        "      <td> 13</td>\n",
        "      <td>  1</td>\n",
        "      <td>  6.583333</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>  NaN</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0</td>\n",
        "      <td> 1418</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>    Raven</td>\n",
        "      <td> 128</td>\n",
        "      <td>NaN</td>\n",
        "      <td>NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>  NaN</td>\n",
        "      <td>  NaN</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>4</th>\n",
        "      <td> 80</td>\n",
        "      <td> 14</td>\n",
        "      <td>  1</td>\n",
        "      <td> 17.083333</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>  NaN</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0</td>\n",
        "      <td> 1402</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>    Raven</td>\n",
        "      <td> 112</td>\n",
        "      <td>NaN</td>\n",
        "      <td>NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>  NaN</td>\n",
        "      <td>  NaN</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>5</th>\n",
        "      <td> 80</td>\n",
        "      <td> 15</td>\n",
        "      <td>  1</td>\n",
        "      <td> 12.083333</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>  NaN</td>\n",
        "      <td> 1</td>\n",
        "      <td> 0</td>\n",
        "      <td> 1207</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>    Raven</td>\n",
        "      <td> 120</td>\n",
        "      <td>NaN</td>\n",
        "      <td>NaN</td>\n",
        "      <td> NaN</td>\n",
        "      <td>  NaN</td>\n",
        "      <td>  NaN</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "<p>5 rows \u00d7 20 columns</p>\n",
        "</div>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "                Paper ID  Group ID   N        Age  Age SD  Age low  Age high  \\\n",
        "Observation ID                                                                 \n",
        "1                     54         1  46   7.500000     3.3      2.9      15.5   \n",
        "2                     80        12   1   7.666667     NaN      NaN       NaN   \n",
        "3                     80        13   1   6.583333     NaN      NaN       NaN   \n",
        "4                     80        14   1  17.083333     NaN      NaN       NaN   \n",
        "5                     80        15   1  12.083333     NaN      NaN       NaN   \n",
        "\n",
        "                Concurrent  Critical period   Phe  Phe SD  Phe low  Phe high  \\\n",
        "Observation ID                                                                 \n",
        "1                        0                1   312     132      125       852   \n",
        "2                        1                0   704     NaN      NaN       NaN   \n",
        "3                        1                0  1418     NaN      NaN       NaN   \n",
        "4                        1                0  1402     NaN      NaN       NaN   \n",
        "5                        1                0  1207     NaN      NaN       NaN   \n",
        "\n",
        "                IQ scale   IQ  IQ SD  IQ low  IQ high  Correlation     p  \n",
        "Observation ID                                                            \n",
        "1               Weschler  104     15      68      143        -0.17  0.38  \n",
        "2                  Raven   91    NaN     NaN      NaN          NaN   NaN  \n",
        "3                  Raven  128    NaN     NaN      NaN          NaN   NaN  \n",
        "4                  Raven  112    NaN     NaN      NaN          NaN   NaN  \n",
        "5                  Raven  120    NaN     NaN      NaN          NaN   NaN  \n",
        "\n",
        "[5 rows x 20 columns]"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "HISTORICAL, CONCURRENT = 0,1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Extract historical observations only"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "observations = observations[observations.Concurrent == HISTORICAL]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Get list of set of unique papers"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "unique_papers = set(observations['Paper ID'])\n",
      "unique_papers"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "{54, 494, 794, 1054, 1156, 1222, 1420, 1429, 1431, 1432, 1448, 2939}"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Re-number papers"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "paper_lookup = dict(zip(unique_papers, range(len(unique_papers))))\n",
      "paper_id = observations['Paper ID'].replace(paper_lookup)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Unique group ID numbers"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "unique_groups = set(observations['Group ID'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "group_lookup = dict(zip(unique_groups, range(len(unique_groups))))\n",
      "group_id = observations['Group ID'].replace(group_lookup)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Individual-level data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "indiv = observations.N == 1\n",
      "obs_indiv = observations[indiv]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Group-summarized data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "obs_summ = observations[indiv - True]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Group and paper IDs for individual data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "group_id_indiv = group_id[indiv].values\n",
      "paper_id_indiv = paper_id[indiv].values"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Unique paper ID's for individual studies"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "unique_papers_indiv = set(paper_id_indiv)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Group and paper IDs for summarized data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "group_id_summ = group_id[indiv-True].values\n",
      "paper_id_summ = paper_id[indiv-True].values"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Unique paper IDs for summarized data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "unique_papers_summ = set(paper_id_summ)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "PKU Model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We developed two meta-analytic models. The first represents the relationship of blood Phe and IQ when Phe was measured historically with respect to IQ measurement (more than 12 months before IQ measurement), while in the second model, Phe and IQ were measured concurrently (Phe measured within 6 weeks prior to IQ measurement). Note that there were no studies in which Phe measurements were taken 6 weeks to 12 months prior to IQ measurements. We believed it to be unreasonable to combine historical and concurrent measurements in the same model, given the potential for a drastically different relationship between the two variables depending on how far apart they were taken. The example below illustrates the historical measurement model."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pymc import *"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For some studies, Phe was only reported as a range. In these cases, we can impute Phe simulaneously as we fit the rest of the model. Here, we model the missing values as uniform random variables across the reported range."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "phe_isnan = obs_indiv.Phe.isnull()\n",
      "phe_low_missing = obs_indiv['Phe low'][phe_isnan]\n",
      "phe_high_missing = obs_indiv['Phe high'][phe_isnan]\n",
      "phe_missing = Uniform('phe_missing', lower=phe_low_missing, upper=phe_high_missing, plot=False)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Linear model\n",
      "\n",
      "The core of each model is a linear relationship between the expected IQ ($\\mu$) and Phe ($x$):\n",
      "\n",
      "$$\\mu_i = \\beta_{0j[i]} + \\beta_{1i} x_i$$\n",
      "\n",
      "\n",
      "The subscript $j[i]$ denotes parameters for study $j$ corresponding to observation $i$. Hence, both the intercept $\\beta_0$ and slope $\\beta_1$ are allowed to vary by study. Note that by *observation* we refer here not to individuals, but to groups of individuals within a study that share a covariate. For example, within the same study, one group of individuals might have been measured for Phe in the critical period, and others not; these groups were considered separate observations in this analysis. One study reported a range of Phe measurements, rather than a single value, so we imputed values by randomly sampling at every iteration from a uniform distribution across the reported range.\n",
      "\n",
      "\n",
      "The intercept was modeled as a random effect, where each study is assumed to be an exchangeable (\\emph{i.e.} conditionally independent) sample from a population of PKU studies:\n",
      "\n",
      "$$\\beta_{0j[i]} \\sim N(\\mu_{\\beta}, \\tau_{\\beta})$$\n",
      "\n",
      "Thus, $\\mu_{\\beta}$ is the expected (population) mean and $\\tau_{\\beta}$ the associated variance, describing the expected differences within the population. This random effect was included to account for a variety of potential sources of variation in IQ, such as the IQ measurement scale in each study or geographic location."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Mean intercept\n",
      "mu_int = Normal('mu_int', mu=100, tau=0.01, value=100)\n",
      "\n",
      "\n",
      "# SD of intercepts\n",
      "sigma_int = Uniform('sigma_int', lower=0, upper=100, value=10.)\n",
      "tau_int = sigma_int**-2\n",
      "    \n",
      "# Intecepts by study\n",
      "beta0 = Normal('beta0', mu=mu_int, tau=tau_int, value=np.zeros(len(unique_papers)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The slope of the relationship included a study-level random effect and a fixed effect corresponding to whether the Phe measurement was taken during the critical period (via indicator function $I$):\n",
      "\n",
      "$$\n",
      "\\begin{aligned}\n",
      "    \\beta_{1i} &= \\alpha_{0i} + \\alpha_1 I(\\textbf{crit}_{j[i]}=1) \\\\\n",
      "    \\alpha_{0i} &= N(\\mu_{\\alpha},\\tau_{\\alpha})\n",
      "\\end{aligned}\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Mean slope\n",
      "mu_slope = Normal('mu_slope', mu=0, tau=0.1, value=-0.01)\n",
      "    \n",
      "# SD of slopes\n",
      "sigma_slope = Uniform('sigma_slope', lower=0, upper=10, value=1.)\n",
      "tau_slope = sigma_slope**-2\n",
      "\n",
      "# Random slopes by study\n",
      "alpha0 = Normal('alpha0', mu=mu_slope, tau=tau_slope, value=np.zeros(len(unique_papers)))\n",
      "alpha1 = Normal('alpha1', mu=0, tau=0.01, value=0.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Calculate Phe effect (slope) for each individual"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@deterministic\n",
      "def beta1_indiv(a0=alpha0, a1=alpha1, crit=obs_indiv['Critical period']):\n",
      "    \"\"\"Calculate Phe effect (slope) for each individual\"\"\"\n",
      "    return a0[paper_id_indiv] + a1*crit"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Calculate expected IQ"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = obs_indiv['Phe'].copy()\n",
      "\n",
      "@deterministic\n",
      "def mu_iq(b0=beta0, b1=beta1_indiv, phe_m=phe_missing):\n",
      "\n",
      "    # Insert values for missing phe\n",
      "    x[phe_isnan] = phe_m\n",
      "\n",
      "    return b0[paper_id_indiv] + b1*x.values"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 26
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Process noise (variance of observations about predicted mean)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sigma_iq = Uniform('sigma_iq', lower=0, upper=100, value=1)\n",
      "tau_iq = sigma_iq**-2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 27
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, the expected value of IQ was used to model the distribution of observed IQ values $IQ_i$, with error described by the variance $\\tau$\n",
      "\n",
      "$$IQ_i \\sim N(\\mu_i,\\tau)$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "iq_like = Normal('iq_like', mu=mu_iq, tau=tau_iq, value=obs_indiv['IQ'].values, observed=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Inference from studies with summarized data\n",
      "\n",
      "Twelve studies provided only summarized data, with no individual measurements of Phe or IQ. For studies that provided only data summaries, we were unable to directly estimate the quantities as specified above. Instead, we employed reported correlation coefficients to obtain inference regarding the relationship of these variables. Inference regarding the linear relationship (slope) between Phe and IQ can be obtained from the correlation coefficient ($\\rho$), using the Fisher transformation. Here, the hyperbolic function can be used to transform the correlation to a normally-distributed random variable:\n",
      "\n",
      "$$\\text{arctanh}(r_J) \\sim N \\left( \\text{arctanh}(\\rho_j), \\frac{1}{\\sqrt{n_j - 3}} \\right)$$\n",
      "    \n",
      "where $r_j$ is the reported Pearson correlation from study $j$, with a standard error that is solely a function of the corresponding sample size $n$ (for a Spearman correlation, the standard error is the inverse square root of $n-2$). This provides a measure of precision for the reported correlations, which in turn becomes a measure of precision for the slope of the relationship between Phe and IQ. The expected value of the slope was obtained in the model by converting $\\rho$ using the fundamental relationship:\n",
      "\n",
      "$$\\beta_{1j} = \\rho_j \\left( \\frac{s_{yj}}{s_{xj}} \\right)$$\n",
      "    \n",
      "where $s_{xj}$ and $s_{yj}$ are the reported standard deviations of the Phe levels and IQs, respectively, for study $j$, the availability of which was an inclusion criteria for the selected studies."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "stdev_phe = obs_summ['Phe SD'].values\n",
      "stdev_iq = obs_summ['IQ SD'].values\n",
      "mean_phe = obs_summ['Phe'].values\n",
      "mean_iq = obs_summ['IQ'].values"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 29
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Calculate mean of slope for summarized data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@deterministic(cache_depth=0)\n",
      "def beta1_summ(a0=alpha0, a1=alpha1, crit=obs_summ['Critical period']):\n",
      "    return a0[paper_id_summ] + a1*crit"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Likelihood for correlation coefficients of summarized data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@potential\n",
      "def r_like(b1=beta1_summ):\n",
      "    \n",
      "    # Convert slope to r\n",
      "    rho = b1*stdev_phe/stdev_iq\n",
      "\n",
      "    # Fisher transformation to allow for normality assumption\n",
      "    eps = np.arctan(rho) - np.arctan(obs_summ['Correlation'])\n",
      "    \n",
      "    # Difference should be mean-zero\n",
      "    return normal_like(eps, mu=np.zeros(len(obs_summ['N'])), tau=obs_summ['N']-3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 31
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Caclulating the probability of IQ dropping below threhsold values\n",
      "\n",
      "In order to evaluate the effect of particular levels of Phe on the likelihood of cognitive impairment, we chose a threshold value of IQ to bound the definition of impairment. We assume that for a standardized measure like IQ, a boundary of one standard deviation below the mean (IQ=85) was a reasonable choice. This threshold value was used to define indicator variables that were set to one if the value of the predicted IQ was below 85 during the current iteration of the MCMC sampler, and zero otherwise. Hence, for each combination of predictors, the total number of ones divided by the number of MCMC iterations represents a posterior predictive probability of observing IQ $< 85$. This corresponds to the integral of the posterior predictive distribution of IQ up to an 85 score. To illustrate the variation of this probability in response to Phe, this probability was calculated for a range of blood Phe levels from 200 to 3000 $\\mu$mol/L, in increments of 200.  This was done for critical period and non-critical period Phe measurement, under both the historical and concurrent measurement models. To assess the performance of the models for a lower threshold value, we also estimated the probability of IQ lower than 70 (two standard deviations below the population mean), characterizing those with intellectual disabilities."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Ranges of predictors\n",
      "crit_pred, phe_pred = np.transpose([[crit, phe] for crit in [0,1] for phe in np.arange(200,3200, 200)])\n",
      "\n",
      "@deterministic(cache_depth=0)\n",
      "def pred(a1=alpha1, mu_int=mu_int, tau_int=tau_int, mu_slope=mu_slope, tau_slope=tau_slope, tau_iq=tau_iq):\n",
      "    \"\"\"Estimate the probability of IQ<85 for different covariate values\"\"\"\n",
      "    \n",
      "    b0 = rnormal(mu_int, tau_int, size=len(phe_pred))\n",
      "    a0 = rnormal(mu_slope, tau_slope, size=len(phe_pred))\n",
      "\n",
      "    b1 = a0 + a1*crit_pred\n",
      "\n",
      "    iq = rnormal(b0 + b1*phe_pred, tau_iq)\n",
      "\n",
      "    return [iq<v for v in (70,75,80,85)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 32
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Run the model using MCMC (Metropolis-Hastings)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "iterations = 10000\n",
      "burn = 9000\n",
      "\n",
      "M = MCMC(locals())\n",
      "M.sample(iterations, burn)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------100%-----------------] 10000 of 10000 complete in 80.4 sec"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Visualization of the estimated probabilities of IQ less than either 85 (darker lines) or 75 (lighter lines) as a function of historical Phe measurement."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "historical = np.sum(M.trace(\"pred\")[:], 0)/float(iterations-burn)\n",
      "\n",
      "x = np.arange(200,3200, 200)\n",
      "\n",
      "colors = ('0.7', 'black')\n",
      "markers = (False, True)\n",
      "\n",
      "for i,p in enumerate((70,85)):\n",
      "    plt.plot(x, historical[i*-1][M.crit_pred==1], color=colors[i], marker='o'*markers[i])\n",
      "    plt.plot(x, historical[i*-1][M.crit_pred==0], color=colors[i], marker='o'*markers[i], linestyle=\"dashed\")\n",
      "plt.xlim(150, 3050)\n",
      "plt.ylim(0,1)\n",
      "\n",
      "plt.xlabel(\"Phe ($\\mu$mol/L)\")\n",
      "plt.ylabel(\"Pr(IQ < threshold)\")\n",
      "plt.text(1400, 0.9, \"critical period\", fontsize=12)\n",
      "plt.text(1700, 0.55, \"non-critical period\", fontsize=12)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 36,
       "text": [
        "<matplotlib.text.Text at 0x109b76d50>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEXCAYAAACzhgONAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8VfX5wPFPJiMJCVlkEbJIWEkgqAwrKir9IVJbKMcB\niiICdVXRgquOV6ugrWKlWsGFCpgeHEAFLSIoQ0ENQ0bCCBCycxMSbva49/z+CLlNyM3OzR0879cr\nL5Nzzv3eJyD3Oec7nq+TpmkIIYQQLXG2dgBCCCFsmyQKIYQQrZJEIYQQolWSKIQQQrRKEoUQQohW\nSaIQQgjRKteeeiNFUd4DpgAFqqrGt3DNkgvXVAB3qaqa1lPxCSGEMK8nnyjeB/6vpZOKotwIJKqq\nmgD8EVjVQ3EJIYRoRY8lClVVdwLFrVzyG+CDC9fuBXwURRnQE7EJIYRomS2NUYQCmY1+zgLCrBSL\nEEKIC2wpUQA4XfSz1BcRQggr67HB7HbIBgY2+jnswjGzvvnmG0kiQgjRCdddd93FN+WtsqVEsRF4\nAEhWFGUsUKKqan5rL0hKSuqRwIQQwlHs27evw6/pyemxHwNXA/6KomQCzwJuAKqqrlBVdbOiKBMU\nRTkElAN391RsQgghWuZkr2XGv/nmG02eKIQQomP27dvX4a4nWxvMFkIIYWMkUQghhGiVJAohhBCt\nkkQhhBCiVZIohBBCtEoShRBCiFZJohBCCNEqSRRCCCFaJYlCCCFEqyRRCNEO69atY/r06S2e/+GH\nHxgzZkyX3ycxMZHvvvuuy+20x6OPPsrf//73Tr126tSpfPTRR90ckbBVtlQUUAibNWPGDGbMmGH6\n2c/Pj5SUFCIiIgAYN24ce/fu7fL7ODk54eTUoeoKnfbKK690+rU9GafomC1btrBixQqqq6vp1asX\n8+fPZ9KkSV1qUxKFEG0wGAy4uLg0O26vddIAjEYjzs7SoeBotmzZwhNPPMHp06dNx86cOQPQpWQh\n/6eIS1ZWVhZ33nknsbGxxMTEsHjxYgDWrl3L5MmTefHFF0lISGDp0qWsXbuWG2+8EYApU6YAMGHC\nBMLDw1m/fj27du1ixIgRbbZ9+vRpbr75ZmJiYhg8eDDz589Hr9e3K97777+fJ554gjvvvJOYmBim\nTp1KVlaW6fzx48f53e9+R1RUFGPGjGH9+vVNXvvkk09y9913Ex0dzc6dO7n//vt54YUXTNd88MEH\nXHbZZURHRzNz5kzy8vJM57Zv386YMWMYPnw4f/3rXwH7TpTtsWXLFqZPn85NN93E9OnT2bJli823\nvWLFiiZJAur/n1u5cmWX2pUnCnFJMhgM3HbbbVx55ZUcPHgQZ2dnDhw4YDq/b98+xo8fz65du3B3\nd+ezzz4zndu0aRN+fn7s3LnT1PW0a9euVtvev3+/6fzChQsZP348x48f54EHHmDp0qW8+OKL7Yp7\n9erVrFu3jnfeeYc///nPzJs3j82bN1NeXs7vfvc75s6dy6pVq9i7dy9z585l6NChxMXFAfDRRx/x\n2muv8fbbb2MwGFBV1dR9tGPHDp577jk2btxIXFwcTz31FHPnzuWLL76gqKiI2bNn8/rrrzN16lTe\nfPNNfvzxRxRF6fSfv62z1J15W23fcMMNVFdXo9frzX6Vlpa2eLy0tNTUzsWqqqq6FLMkCnFJSklJ\nIS0tjQ0bNtCnTx+AJoPRrq6uLF68GHd3925pe+zYsQBERkYSGRkJwPDhw7n77rt577332t12QkKC\nqa17772XMWPGkJ2dzd69e/H09OSRRx4B6j/MJkyYwIYNG1i0aBEAo0ePNg3Iu7o2/ae/bt06rr/+\neuLj4wH4wx/+wOWXX05mZia7d+/Gx8eH3/72twDMnTu33YnN0rraH19aWkpOTg65ublNvtavX09R\nUVGTa0+fPs2dd96Jv78/Li4uODs7m/7b8NWe44cOHaK4uLhZ27NmzTKN/fTr16/Zl5eXF15eXvTr\n14/g4GDi4uKaXbNw4UK+//77Zr9n7969O/cHfIEkCnFJys7OJjo6Gl9fX7Pnhw8f3qkk0VbbBQUF\nPPHEE+zZs4eysjI0TcPHx6fdbTfu3oqJicHDw4O8vDyysrLIyMgwJSGof7IJDg42/Tx69OgW283P\nz2fChAmmn6OiovDy8iInJ4e8vDyGDx9uOtenTx8GDx7c7pgtpbU78+uvvx6dTkdubm6zRJCTk2M6\nZjQaCQkJITg42PQVGxuLv79/s0QB9Yn6vffew2g0YjAYmvy34aut40899VSzRAEwatQoNm7cSK9e\nvTr9Z/LQQw+Rm5vb5M8kIiKCefPmdbpNkEQhLlGhoaGkp6dz7tw5sx/oF99xm9NSH31rbf/lL39B\nr9fz/fff4+3tzYcffsjLL7/c7rgPHTpk+v7EiROUl5cTHBxMaGgoERER7Nmzp8XXtjZ4HRQU1KTr\nLT09ndLSUkJDQwkKCuLw4cOmc5WVlZw4caLdMVtKS/3xs2fPxmg04u3tTXBwcJNEMG7cONPPISEh\n9OvXz+zsra+++opjx441O96vXz/CwsK6FPeAAQM4evRos+NeXl5dShLwv26xlStXUlVVRe/evZk3\nb57MehKiMy677DKGDh3K008/zWOPPUZISAgHDx5s91qIwMBA9u/f3+QOvsHo0aNbbNvDwwNvb29c\nXV1JSUnhnXfe6VDchw4dYs+ePSQlJfHOO+8wZswYQkJCmDRpEs899xzLly9HURR8fX05dOgQnp6e\nxMbGttheQ7KbPn06d999N4cPH2bw4MGsWLGCsWPHEhYWxg033MCiRYvYsGEDU6ZM4d1337XaQLam\naRw8eJDNmze3OB15xIgRfPHFF1360J0/fz5nzpzp9jtzS7cN9cmiq4nhYjLrSVySnJ2dWbt2LaWl\npdxwww3Ex8ebZgmZWyNw8bHFixfz5JNPEhERwYYNG5qcd3FxabHt++67j6qqKhITE/nzn//MnDlz\nOrQe4Y477uDNN99k6NChHD16lBUrVgD1d6Off/45u3fv5le/+hVDhw7lL3/5C7W1tU1+h4s1HLv6\n6qt55plnmD17NsOGDSMrK8uUxPz8/Hj//fd54YUXSEhIoKSkpFsWF7ZXTU0N27dvZ9GiRcTHx3Pv\nvfdSVVXFkCFDzF7v7e3dLXfmS5YsYeLEiYwfP56JEyeydOnSbvkAtmTbliJ7ZgthJ+6//35CQkJ4\n6qmnrB2Kxen1er7++mu+/PJLvvnmGwYPHsyNN97I5MmTiY2NxcnJyewYRUREhM1/6FpbZ/bMlq4n\nIYRNyM7O5quvvmLTpk38/PPPjBs3jsmTJ/PXv/6VoKCgZtdbqj9eNCeJQgg7Yo9lM1qawqppGkeP\nHmXz5s18+eWXZGRkMGnSJO666y4+/PBDPD0922zbEv3xojlJFELYiTfeeMPaIXSYue6h1NRURo0a\nxdGjR9E0jcmTJ/P8888zbty4ds02E81pmkZlZSWlpaWUlZVRWlpKQEAAoaGh3dK+/K0IISzG3BTW\nvLw80tLSWL16NcOGDbPLpyRrq66uNiWEsrIyysrKcHNzw9PTEy8vLwIDA/Hw8Oi295NEIYSwiNTU\nVNLS0syeCw4ObrKIT7Ssrq7OlAwaEoPRaMTLywtPT09CQ0Px9PTEzc3NYjFIohBCdJvS0lI+++wz\nVq9eTU5OTosfXl0tKeGojEYjFRUVTZJCVVUVHh4eeHp64u/vT2RkJL169erRJzFJFEKILtE0jb17\n9/LRRx+xadMmrrrqKh577DGuu+46tm3bZnYKa3ctLrNnmqZRXV3dZFyhvLycXr16mZ4WgoOD6du3\nr9VLwkuiEEJ0Sn5+Pv/+979Zs2YNADNnzuTZZ58lMDDQdM2lMIVV0zTq6uo69eXi4mIq9hceHo6n\np6dNDujbXkRCCJtVV1fH1q1bWb16Nbt372bKlCn84x//YMyYMS12hdjbFFaDwUB1dbXpq60Pe4PB\ngKura4tf7u7u9O3b1+w5Z2dnuxjMl0QhhGhTeno6a9asITk5mbCwMGbNmsW//vUvvLy8rB1ah2ia\nRk1NTZNE0PDVcNxoNOLu7k6vXr3o1asXbm5uuLq60rt3b7Mf9i4uLnbxYd8VkiiEEGYXxf3qV79i\n48aNrF69muPHj6MoCp999lmLNZasraEL6OIP/sZftbW1uLm50atXL1My6NOnDz4+PqbE4Orq6vAf\n/B0liUKIS5y5RXE///wzBoOB8ePHM2/ePP7v//6v0/tzWIrBYODcuXMUFhZSWVlJdXU1Tk5Opg/8\nhq++ffuavnd3d7f6wLA9kkQhxCXMaDTy+uuvN1sUV1payvjx41FV1UqRmadpGufPn6egoIBz587h\n5eVFQEAAHh4epqcB0f3kT1UIO9CVLT81TSM3N5dTp06Rnp7OqVOnTN9nZGRQV1dn4ei7RtM0ysvL\n0el06HQ63N3dCQwMJCIiwuaechyVJAohbFxrW342JAtN09DpdE2SQcN/T58+TZ8+fYiKiiI6Opro\n6GimT59OdHQ0kZGRzJ49m+3btzd7X2sviquqqjIlB6PRSEBAACNGjKBv375WjetSJIlCCBvX0paf\njz/+OKqqmhKCi4sLUVFRxMTEEBUVxdSpU03JoV+/fi22b+kd1zqitraWwsJCdDodlZWV+Pv7ExMT\ng5eXlwwwW1GPJQpFUSYAr114z7dVVV1+0fk+wFtAAqAHXlVVdUNPxSeELaqoqCA/P9/subq6Om64\n4Qbmz59PdHS02b2/28Pai+IMBgPFxcUUFBSg1+vp378/YWFh+Pj4yMCzjeiRRKEoigvwHnA9kA38\npCjKVlVVUxtdNhsoV1V1lKIog4BtiqJsVFXVPrfgE6ITjEYjhw4dYvv27Xz77bfs27cPFxcXs9fG\nxsZyyy23dMv79vSiuIZBaZ1OR1FREZ6engQGBhIbGysD0jaop/5GrgBOqqp6BkBRlGTgZqBxojgP\neCmK4gb4AhWSJMSlICsry5QYduzYga+vL9dccw0LFizgyiuv5IcffnCIekmNB6ULCwtxc3MjICCA\n8PDwLu9xLSyrpxJFKJDZ6OcsoMnu7KqqfqwoylSg8EJc43ooNiF6lF6vZ/fu3Xz77bd8++23nDt3\njgkTJnDttdfy/PPPExYW1uR6a3cNdVVdXR15eXkUFBSYBqWHDx8ug9J2pKcSRZtPBoqiPADUAcFA\nPLBJUZRBqqoaLR2cEN2hpSmsdXV17Nu3z5QYDh8+TFJSEtdeey0rV64kPj6+zb54e6uXBPUD0zk5\nOeTl5dG/f38ZlLZjPZUosoGBjX4eSP1TRWMTgHdVVa0A9iqKkgPEAuZ3PhHChpibwvrLL78QGRnJ\niRMnCA0N5ZprruHRRx9l3LhxDn03XVNTQ3Z2NgUFBfj5+ZGYmGj1qbaia3oqUfwMDFYUJQLIAW4B\nbrvomm+AqYqifA1EAL6qqkqSEHbB3BTWoqIiBgwYwJ49exgwYICVIus5VVVVZGdnU1hYSEBAACNH\njpSxBwfRI3PPVFWtA+YAnwMpwHuqqqYqijJfUZT5Fy5LBgzUJ5V/AX/sidiE6KrCwkJOnDhh9pyP\nj4/DJ4nKykpOnDjBwYMHcXFxYdSoUURFRUmScCA9Ng9NVdXvgFEXHVvR6PvzSHIQdiQrK4t//vOf\nqKraYleSI3e5lJeXk5WVRUlJCcHBwSQlJVl032ZhPbKaRYgOOn78OPfffz8TJkzAzc2N3bt38+qr\nrxIZGdnkOnucwtoepaWlpKamcuTIETw8PBg9ejTh4eGSJByYrGwRop0OHDjAsmXL+OGHH5g7dy4p\nKSn0798fgODgYMB+p7C2h16vJzMzk4qKCkJDQ4mNjW1xMaBwLJIohGiFpmns3r2bZcuWkZaWxv33\n388bb7yBp6dns2vtcQprWxpWUGdmZlJdXU1YWBhDhw6V0hqXGEkUQphhNBr573//y7JlyyguLuah\nhx5CUZRLZoBW0zSKi4vJzMzEYDAQFhZGQECArIG4REmiEKKRuro6Pv/8c1577TVcXV15+OGH+c1v\nfnPJdLFomkZRURFZWfXLnMLCwvDz85ME0c26sr+INdqWRCEE9WsA1q5dy/LlywkNDeX555/nuuuu\nu2Q+II1GI4WFhWRlZeHi4kJ4eDj9+/e/ZH7/ntSe/UVsrW1JFOKScvHd1h133MHZs2d56623SExM\n5K233mLMmDFtN+QgjEYjBQUFZGVl0atXL6KiovD29pYEYUEt7S8yb948EhISuP3227n11lubve7b\nb79l9+7dpr2/G/YBT0pKIiEhodW2V65cKYlCiPYwd7f13XffMXbsWNatW8fw4cOtGF3PMhgM5Ofn\nk52djYeHB7Gxsa1ubiQ6pra2lrS0NPbv34+3tzc333yz6Vx1dbXZ1wwaNIjHHnusWVHIBg0Joqqq\nivPnz1NTU0NVVRUBAQGmRNFS21VVVV36fSRRiEuGubsto9FI7969L5kk0VDJNScnBy8vL4YOHWp2\nBpfoeF//6dOneeutt9i3bx+pqamEhYUxatQoJk+e3OS6liZEBAQEMGHChBbbHzduHOPGtV5Uu6W2\nu7rwUxKFuCRommYaoL1YV++27EFtbS25ubnk5ubi4+PD8OHD8fDwsHZYNqulvn6j0Uh8fDyhoaHN\nXuPm5kZ4eDg333wz8fHxeHl5mW3bklvPWqrtdiUKRVFcgeHAMCAcyKB+06GjqqrWdikCISzswIED\nPPHEE+Tm5po978hlNmpqasjJySE/Px9fX18SEhLo06ePtcOyeS319c+aNYuRI0eydevWZq8JCwvj\n/vvvb7NtS+4vYqm2W00UiqL0A/5AfQ0mZ+Ak9eXBpwIxgKYoyjLgX6qqlnYpEiG6WUFBAX/961/5\n+uuvefLJJ/H39+fpp5+2+53i2qO6uprs7Gx0Oh3+/v4OW+q7o91DdXV16HQ68vPzyc/PJy8vj8LC\nQhYuXNhkAL+lvv7LLruMr776qstxW3JxpiXabuuJYj/wKTC+YRvTxi6UDb8POABEd2tkQnRSTU0N\nK1eu5B//+Ae33nore/fuNQ3Uuri4OHSZjaqqKrKysigqKiIwMJBRo0bh7u5u7bAswlz3UFpaGrfc\ncgtPPvlks723NU0jKiqKvn37MmDAgCZfdXV1TWpVtdTXf6mO57SVKEaoqlrZ0skLyWORoijPdGtU\nQnTS119/zVNPPUVERASbN29m8ODBTc47YpkNgIqKCrKysiguLiYoKMimKrl2ZQFYw/oOX1/fZh/8\n5rqHcnNzWb16NQ8//HCzWVxOTk6cPn26XYsnLTmOYI9aTRStJYmLrnP80UBh006cOMFTTz3FmTNn\neOGFF7jhhhusHVKPKCsrIysrC71eT3BwMFFRUc0+UK2powvAXn75ZY4ePWoaeC8oKKBfv35s27at\n2bTRlrqHWpvq294V9va+T3l3a2uMYudFhzTAqdH3AKiq2vKcLiEsSK/X8/LLL5OcnMzDDz/M6tWr\nHbarpbHS0lIyMzMpLy8nJCSEwYMH22SZkZYGhZctW2b2QzciIoLo6GhCQ0MJDg4mKCioxW4gS00F\nbeCoT5+d0datx7uNvg+iflD7a+rHLpKAicByy4QmRMsMBgNr1qxhyZIlTJo0ie+//57AwEBrh2VR\nDZVcs7KyqKqqIjQ0lCFDhth0JdfDhw+bPd7Sym9FUdrdtnQP9Zy2up5WNXyvKMo3wPWqqh5pdGwY\n8Dqw1FIBCnGxPXv28MQTT9C7d28+/vhjRo4cae2QLKqyshKdTodOp8PJyYnQ0FACAgJsOkE0iI2N\nRafTNTve0o6AHSHdQz2nI52ZI4FTFx07zUXbmwphKVlZWTz//PP88MMPPP/880ybNs1haxLV1NRQ\nWFiITqejuroaf39/YmNj8fT0tJnfWdM09u/fz8aNG+nbty+LFi1qds2DDz5ITk6Oxe76pXuoZ3Qk\nUewGPlQU5V/UT4cdBSwAvrdEYOLSdfEsmbvuuou0tDRWrFjBnDlzeO211xxyVbHBYKCoqAidTkdp\naSm+vr6Eh4fj4+NjM8nBaDSSkpLCxo0b2bhxI+7u7tx8881MmTLF7PVy1+8YOpIo7gbeAL668Lo6\n4LMLx4XoFuZmyezYsYPRo0ezbds2wsPDrRhd99M0jZKSEnQ6HefOnaNfv34EBgYyZMgQmxycLi4u\n5pFHHmHKlCl8/PHHDB06tM0kJnf99q/diUJV1SLgVkVRXIAAQKeqqsFikYlLkrlZMgaDAU9PT4dJ\nEpqmUVZWhk6no7CwkF69ehEYGEhERITVZmxd/BQ3d+5cJk2a1CxZ+fn5sWvXLqvEKKynremxUa2c\nHtQwQ0FV1YvHLoTolJYK9DlC4b7Gg9JQXy00Pj7e6rWXWiq//uc//5k//vGPVoxM2Iq2nihOtqMN\nDbC9Z2Rhd1JSUvjll1/MnrPXOkW1tbUUFhZSUFBgs4PSy5cvN1t+fefOnZIoBND29Fjbn38n7F5F\nRQUvvvgin3zyCXPmzOGLL76w67nxBoOBc+fOodPp0Ov1pkFpb29vm5zSam76KjjGU5wtMxgMGI1G\nXF1dbeamoSUdXuuvKEo4EApkq6p6tvtDEpeSHTt28PDDD3P55Zeze/du/Pz8uPLKK+1ylozRaCQ7\nO5vs7Gy8vLwICAggLi7OJgelGwsNDeX48ePNjtvrU5y9yM/P5+zZsxiNRtzd3XFzc8Pd3Z3AwED8\n/PyaXa9pmtUSSrsThaIowUAycBWgB/opirILuFVV1RwLxScc1Pnz53nmmWfYtm0br7zySpNEYG+z\nZDRN49y5c5w+fRpPT09Gjhxpcx+ylZWVfPLJJ0yfPr3ZYjdZ4WwdISEhhISEYDAYqK2tpaamhtra\n2hZLk5w9e5b8/HxTQmlILn5+fi1uktRdOvJE8RZwGLhHVdWTiqIMpr6kx1vAbywRnHBMmzZtYtGi\nRUyePJndu3fb9V7NFRUVnD59murqamJiYvDx8bF2SE0UFBTw7rvvsmrVKpKSkrj22mubJQpZ62AZ\nmqah1+spKioiMjKyxacBFxcXXFxc2ry5CA8PJzg42JRQampqqKmpQdO0Vl/XHTqSKH4FBKuqWgOg\nquoJRVEeBeRpQrRLQUEBixcv5siRI7z99tuMHz/e2iF1Wl1dHVlZWeTn5xMWFkZwcLBNjT+kp6fz\n2muv8cUXXzBt2jS++OKLZiXXG7O3pzhbVltbS0FBAfn5+QAMGDCgW7qNnJycTE8SPa0jieIc9Vuh\nHmh0bAhQ3K0RCYejaRr//ve/efbZZ5k5cyZvvvmm1aeEdpamaeh0OjIyMvDx8bHZjYHy8/OJiIgg\nJSUFX19fa4dzyTh79iw5OTn4+voSExODl5eXzQ9Ut0dHEsXLwHZFUb4EUoDLgF8DT1giMOEYMjMz\neeSRR9DpdKiqSmJiorVD6rSysjJOnTqF0WgkLi7OprvMxo8fb9dPbPbK39+fkJAQm9oTpDu0+1lZ\nVdW3gelAJXAdUA7MUFV1hYViE3bMaDSycuVKrr32Wq688kq2bt1qt0mitraW9PR0jh49SmBgIImJ\niT2eJLZs2cL06dO56aabmD59Olu2bKGgoIClS5e2OL1VWIamaVRUVJg917dvX4dLEtDB6bGqqm4D\ntlkoFuEgjh07xh//+EecnZ358ssvW+0bt2WappGXl0dmZib+/v4kJSVZ5UPA3Mrpn376CaPRiKIo\nGAxSSacnVFdXm8YeevfuzfDhwx2iW6k9OjI9thdwF/XlxhtKdzoBmqqqd3Z/aMLe1NbW8vrrr/PW\nW2/x+OOPc/fdd9vUAG9H6PV6Tp06hYuLC8OHD7dqtVpz9a/Kysq46qqrePXVV60U1aWjuLiYvLw8\n9Ho9fn5+DBkyBE9PT2uH1aM6cnv0ATAG+C9N96Ww/NwsYXMuLiJ3ww03sGbNGoKDg9m+fXuz/Y3t\nRXV1NRkZGZw/f56IiAj8/f2tftfY0t7Q8iTRM86fP4+vry+xsbE2v3jSUjqSKH4NXKGq6onOvJGi\nKBOA1y6859uqqjbbQlVRlMuBZYAXUKyq6jWdeS9hWS0VkVuwYAF/+ctfrP7B2hlGo5GcnByys7MJ\nCgoiKSmpRz8UioqK2LlzJzt27ECn0/HRRx+Zzll6b2hR381YV1eHm5tbs3MRERE9H5CN6Ui/QBrQ\nqdVEF0qTvwdMA0YD9yiKMvSia3yAVdSv9E4Eft+Z9xKWZ64rxGg0kpaWZpdJori4mP3796PX60lI\nSGDQoEE9kiRqamp4+umnmTBhAklJSSQnJxMTE8Pjjz/e5Lr58+cTGRnZ5JisnO46o9FIcXExJ0+e\n5McffyQnR5aEtaStMuPX8b+upQ+AdxRF+YCmaykaBrlbcwVwUlXVMxfaTQZuBlIbXXM78KmqqlkX\n2ixs5+8gelh5ebnZ4/ZWRK6qqorTp09TUVFBZGRkl9cbXNwdN3/+/FYXsbm7uxMSEsLf/vY3kpKS\nzN7Ngqyc7m7V1dWcPXuWc+fO0adPH/z8/EhMTJQntFa01fX0Lk3HIJyAB81cF2nmWGOhQGajn7Oo\nH+9obDDgpijKTsAT+LuqqmvaaFf0IE3TUFWV/fv3mz1vL//QjEYjmZmZ5OXlERoaSlxcXJcH3c11\nxx07doxrrrmG3Nxcnn/+eUaMGNHsdffdd1+72peV093HxcUFDw8PwsPDW+zWE021VWY8opvepz0D\n3m7ANcD1QF/ga0VRPlNVtbKbYhBdkJ6ezmOPPUZxcTHPPfcc77zzjl0WkaupqSEtLQ03NzdGjhzZ\nbR8U5rrjcnJy2LFjB0uWLGHQoEHd8j6ifQwGAyUlJfj4+DTrRnR1dSUkJMRKkdmnTk8KVxRlGBCg\nqup37bg8GxjY6OeB1D9VNJYJfKmqat6F9n8GJlA/y0pYSU1NDcuXL+df//oXjzzyCPPnz8fV1ZXo\n6Gi76wopLy8nNTWVwMBABg4c2OHxlOrqatLT00lLSyMqKoqRI0c2OWdOeHg4U6ZM6VLcon0MBgPF\nxcUUFhZSUlKCp6cnHh4el+xMpe7UkXUUO4AnVFXdrSjKn4CnAaOiKM+rqvpaGy//GRisKEoE9UUE\nbwFuu+jzKaIzAAAgAElEQVSaDcCHiqL0BXoDo4Dd7Y1PdL89e/awcOFCwsPD2b59OwMH/i/X21tX\nyLlz5zhx4gQZGRksWbKk3eMI27Zt4/333+fYsWNkZmYSHh5OXFwcd9xxR5PrZGaSdWVnZ5OZmYmX\nlxd+fn5ER0e3OOYjOq4jTxTDgT0Xvl8AjKe+S+lt6qe9tkhV1TpFUeYAn/O/6bGpiqLMv3B+haqq\naYqivE99UulN/RhFWYd+G9Etzp8/z3PPPceWLVt44YUXuPnmm+1yNhPUj6s0THvNzc3lpZdeajaO\ncPPNN3PVVVfxf//3f81eHxgYyPTp04mLiyM6OrrFAoCOuqfDsmXLOHPmDP/4xz/Mnl+3bh3Jycl8\n+umnXXofPz8/UlJS2pyK2lIVVj8/PwIDA9uVHBRFYfr06dxyyy0djjMxMZHXX3+dq6++usOvtWdO\n7a1lrijKGSCB+oHp/6iqGqMoijNwXlVVy+6aYcY333yjJSUl9fTbOjRN0/j88895+umnufHGG3nm\nmWdsuvBdW4xGI6dOnaK0tJShQ4cyc+ZMtm/f3uy6wMBAXnjhBaZPn96l99uyZYvddcd1xNmzZxk1\nahQ6na7bV9y3lCg0TaO8vJySkhJKSkpwc3MjLi6uW9+7I0aOHMnrr7/OhAkTrBZDV+3bt4/rrruu\nQ3d+HXmi2Ar8C/CnfjYUQBQg01gdwNmzZ3nsscfIzs5m1apVXHHFFdYOqUtqa2s5duwYzs7OxMfH\n4+rq2uI4QkxMTJeTBNhfd1xbDAaD2f79ntgop6amhtOnT3P+/HlcXFzw8fEhODgYb2/vTrfZELe9\nPh1bU0duCx4Cvge+ABqeQ+MafS/sUF1dHcuXL2fixImMGzeOb7/91u6TRGVlJb/88gseHh4MHTrU\nVMjPlscREhMTef/995k0aRKRkZHcc889TRLbBx98wGWXXUZ0dDQzZ84kLy/PdM7Pz49169YxYcIE\noqKiWLRoUavvlZqayu9+9zuio6MZMmQIy5YtA2Dp0qXMnTuXRx99lKFDh7J27VqWLl3KggULAEyD\n8pGRkYSHh/PTTz+xdu1abrzxxjbbTklJYdKkSURERDB06FAWL15MbW1tizG6uLjg7e1NQkICzz33\nHOvXr+eee+4hLi6OWbNmUVJSYrr2p59+4te//jURERFMmDCB3bv/N7Q5depUXn31VaZNm0ZERARn\nzpxh6tSpppXvmqbx97//ncTEROLi4rjvvvvQ6/Wm1//73/8mISGBkSNHsnLlylb/XB1ZR8qMV6iq\n+oaqqstVVa24cGxTOwayhY3at28fEydOZPv27Xz99dc88sgjdj8AWFJSwqFDhwgNDW22/aQtr3B2\ncnLi/fffZ+nSpWzdupWUlBQ+/vhjAHbs2MFzzz3H+++/T2pqKiEhIcydO7fJ69esWcM777zD559/\nTnJyMt98843Z9yktLWXatGlce+21pKam8vPPPzfpRvnPf/5DXFwc+/fvZ8aMGU3+/DZv3gzAmTNn\nOHv2LJdffnm723Z1dWXJkiWkp6ezatUqvvzyS1555RUOHTpkNmG4uLgQFBRkSuLvvvsuS5cu5dCh\nQzg5OZlWr+fk5PD73/+e2bNnc+jQIe677z5mzpzJuXPnTG29++67LFiwgPT0dNNst4bfa82aNaxZ\ns4b//Oc/7Nu3D71ez+LFiwFIS0tj4cKFrFy5kh9//JFTp06Rm5vb1l+lQ2p3olAUpbeiKC8qinJK\nURT9hWOTFEV5wHLhCUvQ6/U8/vjjzJw5kwcffJBPP/202QeoPcrLy+P48ePExsYSFBTU7PykSZNY\nsmQJEydOZPz48UycOJGlS5faTHfRLbfcQlJSEtHR0UycOJFDhw4B9QPG119/PfHx8bi7u/OHP/yB\nH374gays/80wv+eee4iNjSUxMZErrriCw4cPm32PLVu2UFtbywMPPIC7uzuenp6MHj3adD40NJR5\n8+bRu3dvevfu3aSbqa0up9baTkxMZNCgQZw4UV8q7uqrr+ann34iNDS0zemrTk5OXH311QwZMoS+\nffty9913s379eoxGI+vWrWP06NHcfvvteHl5ceuttzJo0CC2bNliev3EiROZNGkSrq6uzcrEf/LJ\nJ0ybNo3w8HA8PDyYO3cun332GQaDgY0bNzJy5EjGjh2Lu7s7CxYsoK6urtVYHVVHxiiWUb/16ULq\nazIBHKF+xtM/uzcs0R3MlZSora1l8eLFXHvttezevdshtsnUNI3Tp09TXFxMfHy8aZvV06dPs3fv\nXm699VbTtbY8jhAfH2/6PjAwkIyMDKB+W9PGd/1RUVF4eXmRk5NjqtLb+LUDBgwwlVkZN24c2dnZ\nAKiqSnZ2NklJSS0ORjdOGm0pKytDr9dTU1NDRkYGBw8eJC4ujtLS0mZjCSdPnuRPf/oTR48epaqq\nCoPBwMiRI9v9/1/jVe0JCQnU1tZSVFREZmYme/bsaXKjYzAYKCgoAOqTTGu/U15eXpMNtUaOHEld\nXR0FBQXk5eU1ed+IiAi7ntzRFR1JFNOA8aqqpiuK0nBrkQvYZz1pB2eupMTevXvx9vZm5cqVXHnl\nlVaMrvvU1dVx/PhxjEYjiYmJpjvGlJQUZs2a1WZ/vT0ICgriwIH/lVdLT0+ntLS01dXFDXf/P/zw\nQ5Pj2dnZ7Nu3z+xAdXV1tanMel1dHbW1teh0OtNYSUN3TUPb58+fp6KiAqPRiLOzM+Hh4axevdps\n9+Vjjz1mOu/h4cELL7zQZCyhLQ1PVwAHDx7Ezc0Nf39/wsLCuPLKK1m3bl2Lr21ts6ng4GAOHDjA\nb37zGwD279+Pq6srAwYMICgoiG+//dZ07enTp5uMX1xKOjKYXUX9+obGRiKznmySuZISFRUVxMXF\nOUySqKqq4tChQ7i7uzNs2DDTB8KXX37JrbfeyrJly7j77rutHGXnNXwgT58+nW3btnH48GGqq6tZ\nsWIFY8eO7dSeH7/+9a9xd3fn+eefJzc3l9LSUlJSUoD6O3GDwYCTkxN9+/bFz88PLy8v0we/n58f\nzs7OplpfoaGhpnGEgQMHcsstt9C7d2+WLl3arG0PDw98fHxwdnZm586dJCcnd+jPYefOnRw7doyK\nigo++OAD09qeGTNm8NNPP5GcnExJSQlVVVXs2rWrSSXY1rrMpk2bxvr16zl79ixlZWW89957TJs2\nDWdnZ6ZOncqBAwfYs2cPNTU1vP322w65zWl7dCRRrAP+oSjKeABFUa4AXgQ+tkRgomtamgpaU1PT\nw5FYhl6v55dffmHAgAFER0ebulLeffddFi5cSHJystkFdPai8YDr1VdfzTPPPMPs2bMZNmwYWVlZ\nvPPOO02uNfd6czw9Pfn00085dOgQ48eP54orrjDd2Xt4eNCvXz/Cw8MJDg4mICCAPn36mJ48+vbt\ny8KFC7n11luJiori559/bhKnl5dXi20vXryYQ4cOMWLECJYvX869997bJMbWpqw6OTkxZ84cFi9e\nTHx8PAaDgaVLlwL1yerTTz9lzZo1JCUlkZCQwBtvvNEkObTW9qxZs7jtttuYMmUKSUlJeHp68tJL\nLwEwdOhQXn31VebNm8cVV1xBZGTkJVsjqiML7tyBl4B7qS/aV0n9quzFqqqa/1SyIFlw17rJkyez\nd+/eZscnTpzIJ598YoWIuk9BQQGnT59m8ODBTfq4i4uLmTlzJm+88YZDDM5bQsO+z2FhYXaznuA3\nv/kNiqIwa9Ysa4fiECy24O7CxkNPA09QP5gdABSqqmrscJTC4nbs2EFqair+/v4UFv6vZ9BWpoJ2\nlqZpnD17Fp1Ox4gRI5rtY92/f382bdpkNx+APamyspLs7GyKiooIDAzEaDTaVbG8nljkJ1rWrkSh\nqqpBUZTZwCuqqlYBBZYNS3TWBx98wIsvvshHH31EZWWlw5SUMBgMnDhxgpqaGhISElqsuSRJoqny\n8nIyMzM5f/48wcHBrW6QZMvk79W6OjIyswpYrCjKc6qqOkZHtwMxGAw8++yzbNmyhc2bNxMdHQ1g\nt4mhserqatLS0ujTpw8jRozo9jpDjqyqqgovLy9iYmLsdiB248aN1g7hkteR/3Ouo35L04cVRTkI\nNKw80VRVtd8KWQ6gtLSUefPmUVlZyZYtW/Dx6dTW5japrKyM1NRUgoKCmvSrHzp0iB9++MGuu9J6\ngp+fn7VDEA6gI7dm7wDzgPuAFdQXBmz4ElaSmZnJ5MmTCQoKYt26dQ6VJIqKijhy5AiRkZFNNhra\ntm0b06dPJzAw0MoR2gZN0ygsLMRgMFg7FOGg2v1EoarqKgvGITrhxx9/5K677uLBBx9kwYIFDtOP\nazAYyMzMRKfTMWzYMLy8/lfFfu3atTz//PN8+OGHjB071opRWp/RaKSgoIDs7Gzc3Nzw9PS0qwFq\nYT861GmpKMqvqV9k1zDdxIn6rqdnujsw0bpPPvmEJ598kn/+858OMQ7R4Pz585w8eRJPT08SExNN\ng9aapvHyyy+TnJzMf/7zH2JjY60cqfXU1dWRn59PTk4OHh4exMTEdKn8thBt6chWqP8E7qS+1HhD\nCUUn6ne5Ez3EaDSydOlSVFVl/fr1DBs2zNohdYu6ujrOnDlDcXExUVFRzfrWG6rCfvXVVwwYMMBK\nUdqGsrIy02ZMnp6e1g5HXAI68kRxO3CNqqr7LBWMaF1FRQUPPPAA2dnZfP311wQEBFg7pG5x7tw5\n0tPT6d+/P6NGjWLbtm3NihlOmjSJ1atXWztUm+Dj4+NQY1HC9nUkUZwFKiwViGhdXl4es2bNIjo6\nmg0bNtjEZjtdVVtby6lTpygrKyM2NhZvb2+zxQzPnDkDOMZU3/aqqKggPz+f0NDQFteMCNFTWk0U\niqJENfrxb8AbiqK8ARxofJ2qqqcsEJu44JdffmHmzJnMnj2bRx991O4HrTVNQ6fTcebMGQIDA4mJ\niTENwporZnj69GlWrlzp8InCYDBQWFhIfn4+VVVVMqtL2Iy2nihOmjl2rZljsgLKQjZt2sTDDz/M\n3/72N377299aO5wuq66uJj09nerqaoYOHWqa0VRcXMxXX33VpJx2Y1VVVT0ZZo8rKCjg1KlTeHt7\nExYWRv/+/e3+hkA4jlYThaqqkgCsRNM0li9fzooVK1BVlVGjRlk7pC7RNI28vDzOnj1LSEgIQ4YM\nwdnZmX379vGXv/yFffv2cfXVVxMcHExxcXGz1ztCV1trvL29SUpKkm4mYZM6shXq6y0clz2zu1lN\nTY1pi9ItW7bYfZKoqKjg8OHD6HQ64uPjGThwoKkMh6+vL/fccw+pqal8+OGHPPPMMza7r3VXaZrW\n4sY3vXr1kiQhbFZHBrPvBh5qfEBRFKcLxx/uzqAuZUVFRcyePZv+/fuzefPmZhVS7YnRaCQ7O5td\nu3Zx9uxZFi1a1Kw7JSIigoiICNPPDeMQjlLMEOortxYUFFBQUECvXr0YPny4LIwTdqXNRKEoyj0N\n1yqKMof/rZ1wApKAPMuF5/ga72ttMBg4ffo0t99+O08//bRNF78ztx93w4e50Whkx44dJCcns2vX\nLlxcXJg6dSq1tbXtumu25X2tO6KoqIicnBwqKysJCAhg+PDh9O3b19phCdFh7XmiuIP6xOB24fsG\nGpAPzLZAXJcEc1NB/f39GTt2rM0niZamsF533XVMmTKF/Px8pkyZwtq1a4mPj78kB2YNBgPBwcH4\n+vra9N+nEG3pyA53L6iq+pSF42k3R9jhbvr06Wzfvr3ZcVvfha6luCdMmMCf//xnNE1rdc8IIYT1\nWGyHOwBbShKOoqV9rW1xKmhdXR2pqamkpKRw5MgRs9cUFxcTERFxyZS2NhqN6HQ6SkpKiI2NvSSf\nmsSlwT53MnEABoPB1F1zMVubCrpp0yYWLFhAWFgYSUlJ+Pj4UFDQfJNDf3//SyJJ1NbWkpubS15e\nHh4eHoSGhlo7JCEsShKFFVRVVTFv3jz69++Pm5sbGRkZpnPdORW0tQHnxkpKSti3bx96vd7sor4J\nEyZw5MgR+vXrZ2r38ccfb5LoIiIimD9/frfEbcsyMzPJzs7G39+fESNGyOC0uCRIouhher2eWbNm\n4efnx9atW/nuu+8sMhW0tQHncePGkZycTEpKCvv27SMvL4/ExESuv/56s2013g8C4Oqrr2bevHls\n2LABZ2dn+vTpY/dTWNvLx8eHAQMGyPiLuKS0ezDb1tjjYHZ+fj4zZsxg7NixLFmyxKJz6VsbKF+1\nahVPP/00o0ePZvTo0cTFxbU7lvPnz3Ps2DFCQkIIDQ112H55TdMc9ncTlzaLDWYrijIAOA4MVFVV\n3+j4SeBOVVW/71Ckl6BTp07x+9//nttvv71HCvu1NlDu6enJa691bEG9pmnk5uaSlZXF4MGD6d+/\nf3eEaXNqa2vJy8vj3LlzJCQkSLIQgnaW8FBVNR/4Dril4ZiiKFcCRkkSbTt48CA33XQTDz30EI89\n9pjFP3xqa2vJysoye64zA+UGg4ETJ05QUFBAQkKCQyaJiooKTp48SUpKCtXV1QwePFiShBAXdGQV\n0PvAXY1+ng2s6s5gHNF3333HjBkzeOmll7jrrrss/n4ZGRnceOON+Pr6Eh4e3uRcZwbKq6qqOHTo\nEADx8fE2NyOrO5w5c4bDhw/j7u5OUlISMTExMkgtRCMdGcz+AnhLUZQYIAuYDiRaJCoHsX79ehYt\nWsT777/PlVdeafH327BhA3/605/44x//yB/+8Ae2bt3apYHykpISjh8/TlhYGMHBwQ57hx0UFER4\neLisnhaiBR1ZcFerKMoa6p8qfgF+VlXVfP+GGYqiTABeu/Ceb6uquryF6y4HfgAUVVU/a2/7tubd\nd9/l1Vdf5bPPPmPEiBE98p4ZGRkkJyfTMMjf2ZpJmqaRnZ1NTk4OcXFxeHt7d3eoNsURn5KE6E4d\nnR77PrAJOEQHup0URXEB3gOuB7KBnxRF2aqqaqqZ614CvqK+6KDd0TSNJUuW8Nlnn7F582YGDRrU\nY+/90EMPtX1RG+rq6jh58iTV1dUkJibSq1evbojMNpSVldGrVy/c3NysHYoQdqVDz9qqqh4CCoDx\nQEfu9q8ATqqqekZV1VogGbjZzHUPAp8Auo7EZSsMBgMLFy5k69atfPnllz2aJLpDZWUlv/zyC66u\nrsTHxztUktDpdBw5coTy8nJrhyKE3enMgrtFQICqqubnX5oXCmQ2+jkLGNP4AkVRQqlPHhOBy6mv\nTms3GlZbl5aWsmHDhmaL1LqTXq8nOzuboUOHdlub586d4+TJk4SHhxMUFNRt7VqbpmmcOXOGc+fO\nMWLECLve30MIa+nw6J2qqttUVf13B1/Wng/914DHVVVt2OvCbrqe9Ho9M2bMwM3NjeTkZIsmiX37\n9nHNNdewcePGbmlP0zTOnj1Leno6Q4cOdagkUVtby5EjR6ioqCAhIUGShBCd1FPTPLKBgY1+Hkj9\nU0Vjo4FkRVFOUz+j6k1FUX7TQ/F1Wl5eHjfddBPDhg3j7bfftlh3jdFoZPny5dx6660899xzLF68\nuMttNlSEPX/+PImJiRZNcNaQn5+Pp6cnw4YNk3EJIbqgp2o9/QwMVhQlAsihfuHebY0vUFU1quF7\nRVHeB/6jqmr33DZbSHp6Or///e+ZNWsWCxcutNj0UZ1Ox3333Yder2fr1q3N1kd0Rnl5OWlpafTv\n35+IiAiHnBrqyCVGhOhJPfLpoKpqHTAH+BxIAd5TVTVVUZT5iqLYZcnRAwcOMHXqVB5++GGLl+TY\nv38/8fHxfPHFF92SJAoLCzl8+DADBw4kKirKIZMEIElCiG4iRQE74dtvv2XevHksW7aMKVOmWCWG\nztA0jYyMDAoLCxkyZAienp7WDkkI0cMsusOdoii3AROon8G0H9gBfKuqqqFDUdqhxvs6nD9/nrNn\nz/Lxxx8zfvx4a4fWbrW1tRw7dgyAxMREh+qzLy8v59SpUwwdOhRXV6mcL0R3a/NflaIovYCtQByw\nGTgKRAALgG+A2xVFiQLGqKr6seVCtQ5z+zqEhIRQVlbWbe033lzot7/9LXfccUe3tN2grKyMtLQ0\n/P39GTRokEN1yRQWFpKenk5kZKQkCSEspD3/sl4A9ECkqqqm1UqKongCqqIoH1O/9uFRy4RoXStW\nrGiSJABycnJYuXJllzfqMZeEdu7cSf/+/bnpppu61DbUdzXl5+eTkZFBdHQ0/v7+XW7TVjRM69Xp\ndAwfPly60YSwoPaMYt4E3Nc4SQCoqloG3Ef9DKbFqqqutkB8VldRUWH2eFVVVZfbNpeE6urqWLVq\nVZfbrqys5MiRI+Tn5xMfH+9wSSI1NRW9Xk9iYqIkCSEsrD1PFIE0X/PQIBsoVVV1VbdFZENOnTpl\nKrF9se4oJKfX680e70oS0jSNnJwcsrKyCAsLIyQkxKG6mqB+NlNwcDDe3t4OO2NLCFvSnn9lx4Dr\nWjg3kfoxC4ezbds2Jk+ezG233UZkZGSTc53Z18GclhbndTYJlZeX88svv1BcXExCQoJDryPo37+/\nJAkhekh7niheAj5UFOUB4DNVVY2KojhTv3p6OfXdTw5D0zTeeOMN3nzzTVatWsW4cePYsmVLl/Z1\naMlDDz1Ebm5uk+6nziQho9FIZmYm+fn5DBo0iMDAQIdNEEKInteudRSKoswDlgJeQCHgD5QCT6qq\n+pZFI2yBJdZRVFZW8sgjj5CWlsbq1asJCwvrcpuaprF9+3Z8fX0ZOXJks/NdTUJ6vZ6TJ0/Sp08f\noqKiHKriK9SP2TTs8y2E6DqLraNQVXWloihrgQQgjPoxi0OqqpZ2PEzblJ2dzZ133klkZCSbN2/u\n8laYRqOR//73v7zyyiuUlZXx0ksvmb2us5sLGQwG0+K5qKgo/Pz8HO4poqKigtTUVPz9/SVRCGFF\nHdnhrgz43oKxWM3evXuZM2cO8+fP58EHH+zSB67BYGDDhg28+uqruLq6snDhQm666aZu7U8vLi4m\nPT0db29vRo0a5VCL54xGI6WlpRQXF1NQUMCgQYMYMGCAtcMS4pLWngV3O9u4RFNVdUI3xdPjPvzw\nQ/7617/yxhtvcMMNN3S5vYqKCpKTk3n22We5/vrru/Uuv7a2ljNnznD+/Hmio6Pp379/t7VtK44f\nP051dTU+Pj6MGDGiy092Qoiua88TxbttnLfLYlG1tbU89dRTfPfdd2zatInBgwd3S7teXl6oqtot\nbTXQNI2ioiJOnTqFv78/o0aNwsXFpVvfoycZDAYMBgPu7u7NzsXFxTlcF5oQ9q7NROGIayQKCwuZ\nM2cOffv25euvv6Zfv36tXn9xmY358+dz5ZVXkpubS0xMjEVjrampIT09ncrKSoYMGdJmrLZI0zQq\nKyspLi6mpKQEvV7PwIEDzU4WkCQhhO1pNVEoinIrsK61wn+KorgCv1dVNbm7g7OEw4cPM2vWLKZP\nn86TTz7Z5p25uTIbBw4cwGAwMGfOHJ555hmLxKlpGgUFBWRkZDBgwADi4uLsct1AWVkZqampODk5\n4ePjQ1BQEHFxcVKXSQg70ta/1huAJYqirAL2AOnUbzwUAsRQv+/1bOqLA9p8oli/fj1/+tOfeOml\nl5g2bVq7XmOuzEZxcTFjxoyxWJKoqqri5MmT1NXVMXz4cLvYwlPTNLNPA3369GH48OH06dNHnhaE\nsFOtJgpVVe9RFCUMWAi8TH1y6ANUUL9iexvwK1VVcywdaFcYjUaWLFmCqqp8+umnJCQktPu11dXV\nZo9bYozA3spv1NXVce7cOVN3UlJSUrOnHhcXFxmQFsLOtWeMIov6RMGFFdn+gE5VVbsYxNbr9SxY\nsAC9Xs8333zT4eJ43V1moyXl5eWcPHkSZ2dnEhIS6NOnT7e2353KysrIzc2lqKgIb29v+vfvT3h4\nuF12jQkh2taujuIL4xBpwHBVVQssG1L3SU9PZ+bMmVx11VV88MEHra430Ov1fP755/Tt25cZM2aY\njs+fP58zZ850ucxGSzRNIysri5ycHNOaAVt+ioD6yQB9+vQhKSnJ7MwlIYRjae/K7DpFUUqBIcBB\ny4bUPbZu3cr999/Pk08+yezZs81eYzQa2bFjBx9//DH//e9/ufrqq7n33nubXNOwatoStZ4qKys5\nceIEzs7OjBw50m7Kb0RERFg7BCFED+rI1JNlwHJFUd4G9gJ1DSdUVT3V3YF1VOMprPn5+RQVFbF2\n7VrGjh1r9vqioiKuueYa/Pz8uO2223jxxRfx8/Mze21ny2y0pGFG05kzZ2xyLELTNEpKSigrK2Pg\nwIHWDkcIYWUdSRSrLvz3Vxcd1wCrrv4yN4V14MCBLe73AODn58cnn3xCXFxcT4RoUltba1oXMWLE\nCJua0VRbW0tBQQF5eXm4uLgQHBxs7ZCEEDagI7WebHak0twU1szMTFasWEHfvn0ZNGiQ2Tvjnk4S\nJSUlnDhxAj8/P2JjY21q8PfUqVMUFBTg6+tLbGwsnp6eNvWUI4SwnvbUevIAngaGA/uBF1VVNT9n\n1Epyc3PNHt+5cyd5eXn8/e9/t2oXitFoNFV6HTx4MD4+PlaLpSW+vr4MHDjQoQoMCiG6R3ueKP5J\n/Q53XwJ3A37AA5YMqqPy8/PNHu/bty+7du2y6p1xeXk5x48fp0+fPowcOdLqH8RGo9Hsk4wtJi8h\nhG1oT9/HZGCqqqrzgZuBmywbUscFBQWZPW7NrUAbFs8dPnyYkJAQ4uLirJYkNE3j3LlzHDlyhLS0\nNKvEIISwX+1JFJ6qqh4EUFV1P+Bt2ZA6rqVEYa3B2Orqao4ePUphYSEJCQlWWxtRVVVFZmYmKSkp\nZGZmEhAQ0OPjMkII+9eeridnRVEmXvjeCXBt9DMAqqpu6/bIOsDSi+I6oqioiPT0dIKCghg4cKBV\nn2jS0tLw8vIiLi4OLy8vq8QhhLB/7UkUBTTdk6KI5ntURHZbRJ1gyUVx7WUwGDh16hR6vd4myoE7\nOSGNzsgAAA1CSURBVDmRmJgoM5eEEF3mpGl2UbKpmW+++UZLSkqydhgAlJaWcvz4cfr160dkZGSP\nldAuLy9Hp9PRu3fvFrvfhBCisX379nHdddd16A5SNgXoAk3TyMzMJC8vj+jo6BZXdnenmpoaCgsL\nKSgooLa2lsDAQLy9bW7YSAjhQCRRdFLjOk2JiYk9UqepsrKSgwcP4uvrS0REBN7e3tK1JISwOEkU\nHdS4TtPAgQMJDg7usQ/r3r17c/nll9v1ftlCCPsjiaIDeqJOU2VlJTqdjoCAgGZ7Ujg5OUmSEEL0\nOEkU7dRQp8nf37/b6zTV1dWZxh2qqqrw9/e3qTpQQohLmySKNli6TlNxcTHHjx/H29ubsLAwfHx8\nJEkIIWyKJIpWVFRUcPz4cXr16mWxOk0eHh42v/WpEOLS1qOJQlGUCcBrF973bVVVl190fiaw6MKP\nR6ivVHu4J2OE+gHr/Px8MjIyLL49qWwlKoSwdT3Wx6EoigvwHjANGA3coyjK0IsuOwVMUFU1Efgv\n8E5PxdegtraWtLQ08vLyiI+PJygoqMtJQtM0dDodpaWl3RSlEEL0nJ58orgCOKmq6hkARVGSqa9G\nm9pwgaqqPzS6fhPwQg/G12TAOi4urstjBZqmUVxcTEZGBs7OzkRFRXVTpEII0XN6MlGEApmNfs4C\nxrRy/Txgg0UjusBoNHL27Fl0Ol23DViXlJSQkZGB0WgkPDwcX19fWRwnhLBLPZko2l1USlGUa4FZ\nwHjLhVPPEgPWBoOBjIwMQkJC8Pf3lwQhhLBrPZkosoHG+5EOpP6poglFURKAlcBkVVVLLBVM4wHr\n8PDwbhmLaODi4kJiYmK3tCWEENbWk4niZ2CwoigRQA5wC3Bb4wsURQkHPgVmqap60lKB1NbWcvLk\nSaqrq4mPj6dv376dbqulrUWFEMJR9NgnnKqqdcAc4HMgBXhPVdVURVHmK4oy/8JlzwC+wFuKouxX\nFOXH7o6jpKSEAwcO0Lt3bxISEjqdJKqrqzl58iRHjx7t5giFEMK2XDL7UTQesI6JiaF///6det+a\nmhqysrLQ6XQEBQUREhJitb2whRCio2Q/ihZUVlZy7Ngx3N3duzRgnZOTY9p7etSoUbJYTghxSXDo\nRNG4JHh7BqwNBgNVVVU4OzubLanh4eFBYmIivXv3tmTYQghhUxw2UTQuCd7SgLVer6egoIDKykqq\nqqqoq6ujd+/eBAcHm00UspOcEOJS5JCJori4mBMnTuDp6UlAQABlZWVmE4WzszMeHh74+/vTu3dv\nevXqJWsehBDiIg6VKEpKSkhLS8NgMODm5oamadTU1LRYmdXT0xNPT88ejlIIIeyLQyWK7OxsPDw8\nGDx4sIwjCCFEN3GoRBEdHS3dR0II0c0cKlHIU4QQQnQ/qT0hhBCiVZIohBBCtEoShRBCiFZJohBC\nCNEqSRRCCCFaJYlCCCFEqyRRCCGEaJUkCiGEEK2SRCGEEKJVkiiEEEK0ShKFEEKIVkmiEEII0SpJ\nFEIIIVoliUIIIUSrJFEIIYRolSQKIYQQrZJEIYQQolWSKIQQQrRKEoUQQohWSaIQQgjRKkkUQggh\nWiWJQgghRKskUQghhGiVJAohhBCtkkQhhBD/3979x1xZ1nEcfzNMSbG5/BnxWKI2f0TIWigzUWGJ\noqY59k3bFI2aownkWhs/Mu0PTdcyJpaRQ0Hb0I81FVJqDaw0FWQkuqIUJwUiYOVvUlFOf1zXkcPh\nOTfPeXiew3NuPq+Ncf+47nNfX65xf8993fe5LivkRGFmZoWcKMzMrNA+rTpRRIwCZuVz3i5pdidl\nfgicC2wBLpf091bVz8zMOteSO4qI6A/cAVwEfB6YGBHH15UZBwyT9DlgKjCvFXUzM7Nirep6GgGs\nkbRW0lbgHuCCujJfBuYDSFoGHBQRh7eofmZm1kCrEsUngXU16+vztl2VGdzL9TIzs11oVaKodLFc\nv24eZ2ZmvaRVD7NfAjpq1jtIdwxFZQbnbQ2tXLmyRypnZmaNtSpRrACOjYhPAxuArwKX1JVZCFwF\n3BMRpwCvSdrU6APHjBlTf/dhZma9oF+l0prenYg4nR1fj70lIq4EkDQnl7mR9Hrs28AVkla3pHJm\nZtZQyxKFmZm1J/8y28zMCjlRmJlZoZYN4VFGEbEWeAP4ANgqaUREHAjcDQwBXgAulfRWLj8F+Abw\nPjBF0mN7pOIFIuIO0nOizZKG5m1Nx5R/eT8P+CiwSNLMFodSqEGc15FieSUXmyFpcd7XrnF2AHcB\nh5HimidpXpnatCDG6yhRe0bEAOCPwH7AO8C9kn7Sirb0HcXuqQBnSBouaUTedg3weB6K5EngewAR\ncQLwddIQJhcB8yKiL/773wmcXbetmZiqb6PNBybnY4ZHRP1n7mmdxVkBbs7tObzmotLOcW4FrpZ0\nIjAeuDFfJMrUpo1iLFV7SnoHOFPSScDppKGQjqUFbdkXL1Ttpv413Q+HIsl/X5iXLwAWSNoqaS2w\nhjS0SZ8i6VHg1brNzcR0ckR8AjhQ0vJc7q6aY/qEBnHCzu0J7R3nRklP5+V/A0+RRkEoTZsWxAjl\na88teXEg0B94lxa0pbuedk8FWBoR24CfSbodOLzm9x+bgOp4VYNI2b6qs2FM+qpmY9rKjj+ofIn2\niXVyREwEngC+I+k1ShJnRBwDnEiKpZRtWhPjE8CplKw9cy/EX0gxflvSvyKi19vSdxS751RJw4Cv\nATMi4rTanZIqFA9D0nbvJnchpnZ2G3AUMJL03OnHe7Y6PSciBpIG47y62n9dVZY2rYvxbUrYnpK2\n5WvOMcC3ImJ43f5eaUsnit0g6eX892rgflJX0qaIOAIg3+JtzsWbHqKkD2kmpvV5++C67X0+Vkmb\nJVUkvQ78lO1dg20dZ0R8BPg18EtJD+bNpWrTzmIsa3sC5K6kh0nPKnq9LZ0ouiki9s9vGxARhwLj\ngGdJQ5FMyMUmAA/k5YXAxRGxb0QcBRwLLKc9NBWTpI3AGxFxcn54dmnNMX1W/k9GROxDukt8Nu9q\n2zhzveYCf5U0q2ZXadq0UYxla8+IOCQiDsrLBwPn0I1rTnfi9C+zuyn/w9+fV/8DSNKcXbyqNpUd\nX1V7tPU1LxYRC0jfUg4mfTP5PvArmowpv3FxJ7A/8BtJ01scSqGaOA8h9eteC5wBnAS8B/wJuKna\n99vGcX6RFMszbO+SmA78mZK0aYMYZ5DGkytNe0bEUNLD6v7ARtI1Z253rjnNxulEYWZmhdz1ZGZm\nhZwozMyskBOFmZkVcqIwM7NCThRmZlbIicLMzAo5UZiZWSEnCjMzK+REYaUWEWsjYkwPft41ETGt\npz6vi+fsdgwRsSoihvR0nWzv4mHGre1FmmnwMGAb8BawGLgqjyDaY6NpRsQBwDdJw0K00k4xRMQg\nYJmkjhz/RElLOjl2DvBdYFKv19JKy3cUVgYV4DxJA4GxwBjyLF89bCKwRNJ/e+GzmzWOlBChOBne\nDYyPiENaUisrJd9RWKlIWhURi0kTu1R9JiKuB44DfgdMkPRu/lY+GxhFuhO5WdLsgo8fzfaLc6H8\nLf8m4DLSYG3zgR+QBmIbTZp8ZnweArs6h/FtwDDSkM/TJS0qOMU40sxkhSS9GREvkCbxeXBX5c06\n4zsKK4t+AHkil3NIF+Lq9knAFNJ8BCcDl+fhlReRxuf/LOmCPi0izio4x3Gk0Tm7okLqppoEfAm4\nEniE9A1/KHAgaQrL6lwKi4DfA4cCk4F783zIO8nlT8vlu2INcEIXy5rtxInCyqAf8EBEvEoa+n0h\ncEPeVwHmS3pK0vOkO4qTSEljCDBN0qY8/PJ9wMUF5+kANlRXIqJ/RDxWsz43Io6uKX+3pKclPQMs\nA9ZJWpQnvFpI6iIDOAU4Epgl6X1JjwArSXModGYUsCo/g+mK9cCnuljWbCfuerIyqAAXSFraYP/T\nNcsvk6aRPJI0Qf2GiKju60+at6CRf5LmFv5bXh+Zt1UnzxkpqfaOY1XN8ibg+Zr1zcCZeXkQ8Fzd\nhX9F3t6ZccBDBfWs18H2SXvMmuZEYXurdcCbwBGS3uviMauBo9ne5XM26Q4FYDi7vhj3a7B9A+k5\nygE1yeILNZ9d7xzgK12qcXI0oCbKm+3AicL2NtWL9XLgReCGiLiVlDiOBwZIWtHg2KWkZxw/z+tj\ngQV5+XxgSUSMk/Rwk3Vals8/NSJ+RHrwPBy4or5gnllxP0n/qNu1b0QMqFnfKumDiBhIShSPN1kn\nsw/5GYXtbSpARdI24DxS986TwCvAL4CPFRx7BzA6Ij6e50nvAMZGxEXA/4CjgNd3ce4d6gGQ72jO\nB87K9bgVuETSc518xrl03u30MLCl5s+1eftlwH2SXimol1khT4Vq1oSImAl8QHpAfLykmS0+/0PA\nbEm/7WL5VcCFkl7s3ZpZmbnryawJkq4HiIhbSL+NaLU/5D9dImlYr9XE9hq+ozAzs0J+RmFmZoWc\nKMzMrJAThZmZFXKiMDOzQk4UZmZWyInCzMwKOVGYmVkhJwozMyv0f+8ImB85dqMvAAAAAElFTkSu\nQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x109b55710>"
       ]
      }
     ],
     "prompt_number": 36
    }
   ],
   "metadata": {}
  }
 ]
}