{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "%config InlineBackend.figure_format = 'retina'\n", "\n", "from matplotlib import pyplot as plt\n", "from lifelines import CoxPHFitter\n", "import numpy as np\n", "import pandas as pd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Testing the proportional hazard assumptions\n", "\n", "This Jupyter notebook is a small tutorial on how to test and fix proportional hazard problems. \n", "\n", "The proportional hazard assumption is that _all_ individuals have the same hazard function, but a unique scaling factor infront. So the _shape_ of the hazard function is the same for all individuals, and only a scalar infront changes. \n", "\n", "$$h_i(t) = a_i h(t)$$\n", "\n", "At the core of the assumption is that $a_i$ is not time varying, that is, $a_i(t) = a_i$.. In this tutorial we will test this non-time varying assumption, and look at ways to handle violations. \n" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from lifelines.datasets import load_rossi\n", "rossi = load_rossi()\n", "cph = CoxPHFitter()\n", "\n", "cph.fit(rossi, 'week', 'arrest')" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " duration col = 'week'\n", " event col = 'arrest'\n", "number of subjects = 432\n", " number of events = 114\n", " log-likelihood = -658.748\n", " time fit was run = 2019-02-19 17:35:13 UTC\n", " model = untransformed variables\n", "\n", "\n", "---\n", " coef exp(coef) se(coef) z p -log2(p) lower 0.95 upper 0.95\n", "fin -0.379 0.684 0.191 -1.983 0.047 4.398 -0.755 -0.004\n", "age -0.057 0.944 0.022 -2.611 0.009 6.791 -0.101 -0.014\n", "race 0.314 1.369 0.308 1.019 0.308 1.698 -0.290 0.918\n", "wexp -0.150 0.861 0.212 -0.706 0.480 1.058 -0.566 0.266\n", "mar -0.434 0.648 0.382 -1.136 0.256 1.965 -1.182 0.315\n", "paro -0.085 0.919 0.196 -0.434 0.665 0.589 -0.469 0.299\n", "prio 0.091 1.096 0.029 3.194 0.001 9.476 0.035 0.148\n", "---\n", "Concordance = 0.640\n", "Log-likelihood ratio test = 33.266 on 7 df, -log2(p)=15.370\n" ] } ], "source": [ "cph.print_summary(model=\"untransformed variables\", decimals=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Checking assumptions with `check_assumptions`\n", "\n", "New to lifelines 0.16.0 is the `CoxPHFitter.check_assumptions` method. This method will compute statistics that check the proportional hazard assumption, produce plots to check assumptions, and more. Also included is an option to display advice to the console. Here's a breakdown of each information displayed:\n", "\n", " - Presented first are the results of a statistical test to test for any time-varying coefficients. A time-varying coefficient imply a covariate's influence _relative to the baseline_ changes over time. This implies a violation of the proportional hazard assumption. For each variable, we transform _time_ four times (these are common transformations of time to perform). If _lifelines_ rejects the null (that is, _lifelines_ rejects that the coefficient is not time-varying), we report this to the user.\n", " - Some advice is presented on how to correct the proportional hazard violation based on some summary statistics of the variable. \n", " - As a compliment to the above statistical test, for each variable that violates the PH assumption, visual plots of the the _scaled Schoenfeld residuals_ is presented against the four time transformations. A fitted lowess is also presented, along with 10 bootstrapped lowess lines (as an approximation to the confidence interval of the original lowess line). Ideally, this lowess line is constant (flat). Deviations away from the constant line are violations of the PH assumption. \n", " \n", "#### Why the _scaled Schoenfeld residuals_?\n", " \n", "This section can be skipped on first read. Let $s_{t,j}$ denote the scaled Schoenfeld residuals of variable $j$ at time $t$, $\\hat{\\beta_j}$ denote the maximum-likelihood estimate of the $j$th variable, and $\\beta_j(t)$ a time-varying coefficient in (fictional) alternative model that allows for time-varying coefficients. Therneau and Grambsch showed that. \n", "\n", "$$E[s_{t,j}] + \\hat{\\beta_j} = \\beta_j(t)$$\n", "\n", "The proportional hazard assumption implies that $\\hat{\\beta_j} = \\beta_j(t)$, hence $E[s_{t,j}] = 0$. This is what the above proportional hazard test is testing. Visually, plotting $s_{t,j}$ over time (or some transform of time), is a good way to see violations of $E[s_{t,j}] = 0$, along with the statisical test. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "The ``p_value_threshold`` is set at 0.05. Even under the null hypothesis of no violations, some covariates\n", "will be below the threshold (i.e. by chance). This is compounded when there are many covariates.\n", "\n", "Similarly, when there are lots of observations, even minor deviances from the proportional hazard\n", "assumption will be flagged.\n", "\n", "With that in mind, it's best to use a combination of statistical tests and eyeball tests to\n", "determine the most serious violations.\n", "\n", "\n", "\n", " test_name = proportional_hazard_test\n", " null_distribution = chi squared\n", "degrees_of_freedom = 1\n", "\n", "---\n", " test_statistic p -log2(p)\n", "age km 11.03 <0.005 10.12\n", " rank 11.09 <0.005 10.17\n", "fin km 0.02 0.89 0.17\n", " rank 0.02 0.90 0.16\n", "mar km 0.60 0.44 1.19\n", " rank 0.67 0.41 1.27\n", "paro km 0.12 0.73 0.45\n", " rank 0.14 0.71 0.49\n", "prio km 0.02 0.88 0.18\n", " rank 0.02 0.88 0.18\n", "race km 1.44 0.23 2.12\n", " rank 1.46 0.23 2.14\n", "wexp km 7.48 0.01 7.32\n", " rank 7.18 0.01 7.08\n", "\n", "\n", "1. Variable 'age' failed the non-proportional test: p-value is 0.0009.\n", " Advice: try binning the variable 'age' using pd.cut, and then specify it in `strata=['age']` in the call in `.fit`. See documentation in link [A] and [B] below.\n", " Advice: try adding an interaction term with your time variable. See documentation in link [A] and specifically link [C] below.\n", "\n", "2. Variable 'wexp' failed the non-proportional test: p-value is 0.0063.\n", " Advice: with so few unique values (only 2), you can try `strata=['wexp']` in the call in `.fit`. See documentation in link [A] and [B] below.\n", "\n", "---\n", "[A] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html\n", "[B] https://lifelines.readthedocs.io/en/latest/Survival%20Regression.html#checking-the-proportional-hazards-assumption\n", "[C] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Option-2:-introduce-time-varying-covariates\n", "\n" ] } ], "source": [ "cph.check_assumptions(rossi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively, you can use the proportional hazard test outside of `check_assumptions`:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " test_name = proportional_hazard_test\n", " time_transform = rank\n", " null_distribution = chi squared\n", "degrees_of_freedom = 1\n", " model = untransformed variables\n", "\n", "---\n", " test_statistic p -log2(p)\n", "age 11.094 0.001 10.173\n", "fin 0.017 0.896 0.158\n", "mar 0.666 0.414 1.271\n", "paro 0.138 0.711 0.493\n", "prio 0.023 0.881 0.183\n", "race 1.462 0.227 2.141\n", "wexp 7.180 0.007 7.084\n" ] } ], "source": [ "from lifelines.statistics import proportional_hazard_test\n", "\n", "results = proportional_hazard_test(cph, rossi, time_transform='rank')\n", "results.print_summary(decimals=3, model=\"untransformed variables\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the advice above, we can see that `wexp` has small cardinality, so we can easily fix that by specifying it in the `strata`. What does the `strata` do? Let's go back to the proportional hazard assumption.\n", "\n", "In the introduction, we said that the proportional hazard assumption was that \n", "\n", "$$ h_i(t) = a_i h(t)$$\n", "\n", "In a simple case, it may be that there are two subgroups that have _very_ different baseline hazards. That is, we can split the dataset into subsamples based on some variable (we call this the stratifying variable), run the Cox model on all subsamples, and compare their baseline hazards. If these baseline hazards are _very_ different, then clearly the formula above is wrong - the $h(t)$ is some weighted average of the subgroups' baseline hazards. This ill fitting average baseline can cause $a_i$ to have time-dependent influence. A better model might be:\n", "\n", "$$ h_{i |i\\in G}(t) = a_i h_G(t)$$\n", "\n", "where now we have a unique baseline hazard _per_ subgroup $G$. Because of the way the Cox model is designed, inference of the coefficients is identical (expect now there are more baseline hazards, and no variation of the stratifying variable within a subgroup $G$). \n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " duration col = 'week'\n", " event col = 'arrest'\n", " strata = ['wexp']\n", "number of subjects = 432\n", " number of events = 114\n", " log-likelihood = -580.89\n", " time fit was run = 2019-02-19 17:30:49 UTC\n", "\n", "---\n", " coef exp(coef) se(coef) z p -log2(p) lower 0.99 upper 0.99\n", "fin -0.38 0.68 0.19 -1.99 0.05 4.42 -0.87 0.11\n", "age -0.06 0.94 0.02 -2.64 0.01 6.91 -0.12 -0.00\n", "race 0.31 1.36 0.31 1.00 0.32 1.65 -0.49 1.10\n", "mar -0.45 0.64 0.38 -1.19 0.23 2.09 -1.44 0.53\n", "paro -0.08 0.92 0.20 -0.42 0.67 0.57 -0.59 0.42\n", "prio 0.09 1.09 0.03 3.16 <0.005 9.33 0.02 0.16\n", "---\n", "Concordance = 0.61\n", "Log-likelihood ratio test = 188.99 on 6 df, -log2(p)=124.17\n" ] } ], "source": [ "cph.fit(rossi, 'week', 'arrest', strata=['wexp'])\n", "cph.print_summary()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "The ``p_value_threshold`` is set at 0.05. Even under the null hypothesis of no violations, some covariates\n", "will be below the threshold (i.e. by chance). This is compounded when there are many covariates.\n", "\n", "Similarly, when there are lots of observations, even minor deviances from the proportional hazard\n", "assumption will be flagged.\n", "\n", "With that in mind, it's best to use a combination of statistical tests and eyeball tests to\n", "determine the most serious violations.\n", "\n", "\n", "\n", " test_name = proportional_hazard_test\n", " null_distribution = chi squared\n", "degrees_of_freedom = 1\n", "\n", "---\n", " test_statistic p -log2(p)\n", "age km 11.29 <0.005 10.32\n", " rank 4.62 0.03 4.99\n", "fin km 0.02 0.90 0.16\n", " rank 0.05 0.83 0.28\n", "mar km 0.53 0.47 1.10\n", " rank 1.31 0.25 1.99\n", "paro km 0.09 0.76 0.40\n", " rank 0.00 0.97 0.05\n", "prio km 0.02 0.89 0.16\n", " rank 0.02 0.90 0.16\n", "race km 1.47 0.23 2.15\n", " rank 0.64 0.42 1.23\n", "\n", "\n", "1. Variable 'age' failed the non-proportional test: p-value is 0.0008.\n", " Advice: try binning the variable 'age' using pd.cut, and then specify it in `strata=['age']` in the call in `.fit`. See documentation in link [A] and [B] below.\n", " Advice: try adding an interaction term with your time variable. See documentation in link [A] and specifically link [C] below.\n", "\n", "---\n", "[A] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html\n", "[B] https://lifelines.readthedocs.io/en/latest/Survival%20Regression.html#checking-the-proportional-hazards-assumption\n", "[C] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Option-2:-introduce-time-varying-covariates\n", "\n" ] } ], "source": [ "cph.check_assumptions(rossi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since `age` is still violating the proportional hazard assumption, we need to model it better. From the residual plots above, we can see a the effect of age start to become negative over time. This will be relevant later. Below, we present two options to handle `age`. \n", "\n", "#### Option 1: bin variable and stratify on it\n", "\n", "\n", "The first option proposed is to bin the variable into equal-sized bins, and stratify like we did with `wexp`. There is a trade off here between estimation and information-loss. If we have large bins, we will lose information (since different values are now binned together), but we need to estimate less new baseline hazards. On the other hand, with tiny bins, we allow the `age` data to have the most \"wiggle room\", but must compute many baseline hazards each of which has a smaller sample size. Like most things, the optimial value is somewhere inbetween." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ageage_strata
027(24, 27]
118(15, 18]
219(18, 21]
323(21, 24]
419(18, 21]
\n", "
" ], "text/plain": [ " age age_strata\n", "0 27 (24, 27]\n", "1 18 (15, 18]\n", "2 19 (18, 21]\n", "3 23 (21, 24]\n", "4 19 (18, 21]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rossi_strata_age = rossi.copy()\n", "rossi_strata_age['age_strata'] = pd.cut(rossi_strata_age['age'], np.arange(0, 80, 3))\n", "\n", "rossi_strata_age[['age', 'age_strata']].head()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# drop the orignal, redundant, age column\n", "rossi_strata_age = rossi_strata_age.drop('age', axis=1)\n", "cph.fit(rossi_strata_age, 'week', 'arrest', strata=['age_strata', 'wexp'])" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " duration col = 'week'\n", " event col = 'arrest'\n", " strata = ['age_strata', 'wexp']\n", "number of subjects = 432\n", " number of events = 114\n", " log-likelihood = -392.443\n", " time fit was run = 2019-02-19 17:30:49 UTC\n", " model = stratified age and wexp\n", "\n", "\n", "---\n", " coef exp(coef) se(coef) z p -log2(p) lower 0.99 upper 0.99\n", "fin -0.395 0.674 0.197 -2.004 0.045 4.472 -0.903 0.113\n", "race 0.280 1.324 0.313 0.895 0.371 1.431 -0.527 1.088\n", "mar -0.194 0.824 0.392 -0.494 0.621 0.687 -1.202 0.815\n", "paro -0.163 0.849 0.200 -0.818 0.413 1.275 -0.678 0.351\n", "prio 0.080 1.084 0.028 2.854 0.004 7.857 0.008 0.153\n", "---\n", "Concordance = 0.582\n", "Log-likelihood ratio test = 565.874 on 5 df, -log2(p)=396.379\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvkAAAIPCAYAAADpd8F8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmYZWldH/Dvj31QaBhUdu0OEWdMDM0MiiwyRasIKAiKio2GQstA3Bi3BC2iA6aCuw2KAS1Jo1KAARQ1YTOXblncZoYiJhkRtVqUGRhwoNmGdd78cU5BTVHVXV3bvffM5/M89zlTZ3nP78w9fet7T73nPdVaCwAAMBw3G3cBAADA7hLyAQBgYIR8AAAYGCEfAAAGRsgHAICBEfIBAGBghHwAABgYIR8AAAZGyAcAgIER8gEAYGCEfAAAGBghHwAABkbIBwCAgRHyAQBgYIR8AAAYGCEfAAAG5hbjLmAaVNVKktsnOTXmUgAAGLaDST7QWju0k0aE/K25/XnnnXf+hRdeeP64CwEAYLiuuuqqXH/99TtuR8jfmlMXXnjh+VdcccW46wC4SVhcXEySzM3NjbkSgP118cUX58orrzy103b0yQdg4oxGo4xGo3GXATC1hHwAABgYIR8AAAZGyAcAgIER8gEAYGCEfAAAGBhDaAIwcQ4d2tEzYABu8oR8ACbOwsLCuEsAmGq66wAAwMAI+QAAMDC66wAwcY4ePZokWVpaGnMlk+PSSy/N8vLyOW93+PDhHDt2bA8qAiaZkA/si+PHj+fUqVM5ePBgZmdnx10OTJ3l5eWcPHly3GXATcIQfmcJ+cC+OH78eE6ePJlLLrlkaj8wYRIcOHAghw8fPut6y8vLOX369D5UBMMzhN9ZQj4ATJHDhw/nxIkTZ11vZmbGlX+4Cdu1G2+r6mBVtao6XlX3rqqXVtW1VXVDVc1U1cVV9eyqemtVXVdVH62qt1fVL1bVHc/Q7rdV1f9as82pqnpxVd1vg3W/vapeX1Xv79e9qqqeXlW33q3jBACASbcXV/LvleTPk/xNkhclOS/JB5L8uySPTXIyyR+n+4JxcZIfTvKIqrp/a+2Dq41UVSX5b0memOS9SV6R5D1J7pHkoUneluTyNeu/IMmTkvxTkpcneX+Sr0zy00m+uqq+trX2yT04XgAAmCh7EfIfnORZrbWfWDuzqp6V5Ptaa59aN/+7kywm+d4kP7tm0fekC/h/meRrW2un12xz8yRfsObn2XQB//eSPKG1dv2aZZcl+akk35fk2WcqvKqu2GTRBWfaDti65eXlzMzMjLsMJtxVV12VJM6VNbYzss7qdv4/wrnZ7r+3SbIXIf/dSZ6xfmZr7R82Wf8FSX4pydflxiH/B/rpk9cG/L6tTyW5Zs2spyb5ZJLvWhvwez+d5PuTPCFnCfnA3jt9+rR+wmzZtddeO+4Spp5/c3DTtBch/62ttY+tn1lVt0zy5CSPT/KlSQ7kxvcE3H3Nup+T5F8neXdr7S1n2llV3TbJfdJ16bm06+XzWT6W5MKzFd5au3iTfVyR5KKzbQ+c3VZHBgFubLuj5fg3B+duCKNT7UXIf9cm81+ark/+3yd5Zb/e6peBS5OsvTn2Dv30nVvY3x2TVJLPT9ctB5hgWx0ZBLix7Y6W498cnLshjE61FyG/rZ/Rj4Tz2HQ33D5i7Q2wVXWzJP9h3Sbv76d3z9mtfs16S2vN1XaAARiNRkmSI0eOjLkSgOm0X+Pk/8t++gcbjHDzFelG4Pm01tqHq+r/JPnXVXXfM3XZaa19qKr+b5J/VVXnt9au29XKAdh3i4uLSYR8gO3ar5B/qp/OJPmV1ZlV9QVJnrvJNs9J8utJnt8Pf7l2dJ2bJblza2315ttfSvKbSV5QVbOttfevbagfh/9Qa+3KXTgWABibrY6WM4TRQYDt26+Q/5dJ3pTkm6rqzUnemOTOSR6Rbrz7qzfYZjHJVyX5ziRvr6pXphsn/25JjqQbleeyJGmtvaCqLk43DOffVdVrkrwjyflJDiV5SLox95+yR8cHnMXs7GxmZmZy8ODBcZcCU81oObD3hvA7a19CfmvtU1X16CT/Ockjk/xguptqF/t5/2+DbVqSf9sH9n+X5FvT3Zx7TZI3JPmDdet/X1W9Kl2Q/5p0N+9ely7s/3yS39mTgwO2ZHZ2dtwlwFTb7gg5RtaBczeE31m7FvJba6fSjXKz2fLr0l1p38jBM2z3onRPzt1KDX+U5I+2si4ATJNjx46NuwRgitzs7KsAAADTRMgHAICB2a8bbwFgy5aWlsZdAsBUcyUfAAAGRsgHAICBEfIBmDjz8/OZn58fdxkAU0uffAAmzsrKyrhLAJhqruQDAMDACPkAADAwQj4AAAyMkA8AAAMj5AMAwMAYXQeAiXPkyJFxlwAw1YR8ACbO3NzcuEsAmGq66wAAwMAI+QBMnJWVFQ/EAtgB3XUAmDjz8/NJkqWlpTFXAjCdXMkHAICBEfIBAGBghHwAABgYIR8AAAZGyAcAgIER8gEAYGAMoQnAxFlYWBh3CQBTTcgHYOIcOnRo3CUATDXddQAAYGCEfAAmzuLiYhYXF8ddBsDUEvIBmDij0Sij0WjcZQBMLSEfAAAGRsgHAICBEfIBAGBghHwAABgYIR8AAAbGw7AAmDgehgWwM0I+ABNnYWFh3CUATDXddQAAYGCEfAAAGBghH4CJc/To0Rw9enTcZQBMLSEfAAAGRsgHAICBEfIBAGBghHwAABgYIR8AAAZGyAcAgIHxxFsAJs7c3Ny4SwCYakI+ABPnyJEj4y4BYKrprgMAAAMj5AMwcUajUUaj0bjLAJhauusAMHEWFxeT6LYDsF2u5AMAwMAI+QAAMDBCPgAADIyQDwAAAyPkAwDAwAj5AAAwMIbQBGDiLC0tjbsEgKnmSj4AAAyMkA8AAAMj5AMwcebn5zM/Pz/uMgCmlj75AEyclZWVcZcAMNVcyQcAgIER8gEAYGCEfAAAGBghHwAABkbIBwCAgTG6DgAT58iRI+MuAWCqCfkATJy5ublxlwAw1XTXAQCAgRHyAZg4KysrHogFsAO66wAwcebn55MkS0tLY64EYDq5kg8AAAMj5AMAwMAI+QAAMDBCPgAADIyQDwAAAyPkAwDAwBhCE4CJs7CwMO4SAKaakA/AxDl06NC4SwCYarrrAADAwAj5AEycxcXFLC4ujrsMgKkl5AMwcUajUUaj0bjLAJhaQj4AAAyMkA8AAAMzdaPrVFVLcrK1NjPuWgCYTJdeemmWl5fPebvDhw/n2LFje1ARwP6aupAPwE3H8ePHc+rUqRw8eDCzs7Nb3m55eTknT57cu8Ky/doA9sM0hvwLk3xk3EUAsPeOHz+ekydP5pJLLtlWkD5w4EAOHz581vWWl5dz+vTpfa0NYC9NXchvrf31uGsAYG+tPgzr6quv3lE7hw8fzokTJ8663szMzJ5f+QfYT3t6421VHayqVlXHq+qCqvr9qrquqj5cVW+sqoetW3+2X3+2qh5eVSeq6nTfD391nVZVJzbY14GqelZVva2qPlpV76uq11TV1+zlMQKw+xYWFrKwsDDuMgCm1n6NrnMoyZ8mOT/J85P89yQXJ3lVVX3bBus/LskfJflgkucleemZGq+qOyR5c5KnJTmd5FiSlyd5QJLXVtWTd+cwAABg8u1Xd52HJPmF1tqPrc6oql9NF/yfV1Wvaq19YM36j0zyyNbaq7fY/s8m+dIkv57kKa211u/jZ5NcnuQ5VfWa1tqpMzVSVVdssuiCLdYBwB5YXl7OzMzMOa2/1/vZ7j4A9sN+hfzTSZ65dkZr7fKqelGSJyZ5bJIXrln8yq0G/Kq6VZLvSPKhJD++GvD7fby9qp6T5OlJ/u36GgCYTEePHr3Rz6dPn96XPvP7tR+AvbZfIf/K1toHN5h/Il3Iv29uHPL/4hza/pIkt03yptbadRssH6UL+fc9W0OttYs3mt9f4b/oHGoCYBdtdZScVdsZLedc97PdfQDsh/0K+e/eZP67+umBTeZvxeq212yyfHX+Hc6hTQAmyFZHyVm13dFyzmU/RuQBJtl+3Xh7503m36Wfrr8U0taveAar295lk+V33WQfAAAwSPsV8i+qqtttMH+mn75lB22/Ld3Dse7Tj7Kz3kP76ZU72AcAAEyN/equcyDJTyZZO7rO/ZI8Id0V9t/bbsOttY/3N/B+T5KfTvIDa/ZxryQ/mOQTSX57u/sAYDptdbQcI+UAQ7NfIf9PksxV1f2TvCldF5pvS/eXhCevGz5zO56W5KuSfH9VfXmS1yf5vCTfmuR2Sb6/tbayw30AsM9mZ2czMzOTgwcPbmv7vRwtZ6e1Aeyl/Qr5K0mekuRn+umt03WfeWZr7TU7bby1dl1VPSDJjyf5piQ/nOT6dKP0/Hxr7bU73QcA+2dubi5JcuTIkW1tfy4j8Wx3u9nZ2W3tA2A/1Jph5Xe/8aqD6QL+C1trs3u2oz1WVVdcdNFFF11xxWbPygIAgJ27+OKLc+WVV1652dDuW7VfN94CAAD7RMgHYOKMRqOMRqNxlwEwtfarTz4AbNni4mKS7ffJB7ip29OQ31o7laT2ch8AAMCN6a4DAAADI+QDAMDACPkAADAwQj4AAAyMkA8AAANjCE0AJs7S0tK4SwCYaq7kAwDAwAj5AAAwMEI+ABNnfn4+8/Pz4y4DYGrpkw/AxFlZWRl3CQBTzZV8AAAYGCEfAAAGRsgHAICBEfIBAGBghHwAABgYo+sAMHGOHDky7hIAppqQD8DEmZubG3cJAFNNdx0AABgYIR+AibOysuKBWAA7oLsOABNnfn4+SbK0tDTmSgCmkyv5AAAwMEI+AAAMjJAPAAADI+QDAMDACPkAADAwQj4AAAyMITQBmDgLCwvjLgFgqgn5AEycQ4cOjbsEgKmmuw4AAAyMkA/AxFlcXMzi4uK4ywCYWkI+ABNnNBplNBqNuwyAqSXkAwDAwAj5AAAwMEI+AAAMjJAPAAADI+QDAMDAeBgWABPHw7AAdkbIB2DiLCwsjLsEgKmmuw4AAAyMkA8AAAMj5AMwcY4ePZqjR4+OuwyAqSXkAwDAwAj5AAAwMEI+AAAMjJAPAAADI+QDAMDACPkAADAwnngLwMSZm5sbdwkAU03IB2DiHDlyZNwlAEw13XUAAGBghHwAJs5oNMpoNBp3GQBTS3cdACbO4uJiEt12ALbLlXwAABgYIR8AAAZGyAcAgIER8gEAYGCEfAAAGBghHwAABsYQmgBMnKWlpXGXADDVXMkHAICBEfIBAGBghHwAJs78/Hzm5+fHXQbA1NInH4CJs7KyMu4SAKaaK/kAADAwQj4AAAyMkA8AAAMj5AMAwMAI+QAAMDBG1wFg4hw5cmTcJQBMNSEfgIkzNzc37hIAppruOgAAMDBCPgATZ2VlxQOxAHZAdx0AJs78/HySZGlpacyVAEwnV/IBAGBghHwAABgYIR8AAAZGyAcAgIER8gEAYGCEfAAAGBhDaAIMzKWXXprl5eVz3u7w4cM5duzYHlR07hYWFsZdAsBUE/JhII4fP55Tp07l4MGDmZ2dHXc5jNHy8nJOnjw57jJ25NChQ+Mu4SbHZwgMi5APA3H8+PGcPHkyl1xyiV/QJEkOHDiQw4cPn3W95eXlnD59eh8qYpL5DIFhEfIBBurw4cM5ceLEWdebmZmZuCv/i4uLSZK5ubkxVwIwnXZ8421VHayqVlXHq+qCqvr9qrquqj5cVW+sqoetW/9AVf1YVY2q6p+q6uNV9Z6q+oOqesAm+2hVdaKq7lJVi1X1zqr6VFXNrlnnrlX13Ko6tabNV1TVxTs9RgD212g0ymg0GncZAFNrN6/kH0ryp0n+Ksnzk9w1ybcleVVVHW2tvbRf78IkC0n+JMn/SPK+JF+Y5NFJHlFVj2qtvXqD9s9P8mdJPpTkFUluSPLuJKmqQ0nemORuSUZJXpzknkm+JcnXV9U3t9b+aBePFQAAJtZuhvyHJPmF1tqPrc6oql9NF/yfV1Wvaq19IMlVSe7WWnvv2o2r6h5J/iLJLyfZKOR/WZLfTvJdrbVPrlv2vHQB/+mttU8PyVBVv5buy8QLq+qLWmsfOtMBVNUVmyy64EzbwSRZXl7OzMzMuMtgjLYzss7qdpNy7lx11VVJMjH13BRs97wBJtNuhvzTSZ65dkZr7fKqelGSJyZ5bJIXttY2vLurtfZPVfWyJD9QVV/YWnvHulU+nuRH1wf8/svBw5K8I8nPrWvzzVX14iTfkeSbkvzWto8OpsTp06cnrn8102ESz51rr7123CUATKXdDPlXttY+uMH8E+lC/n2TvDBJqupBSZ6a5AFJviDJrdZtc/d0oX2tU621jT7t79tP39Ba+8QGy0fpQv59c5aQ31rbsP9+f4X/ojNtC5NiqyOqMFzbHS1nks6d1Sv5F1544ZgruekwyhIMy26G/HdvMv9d/fRAklTVY5O8LMlHk7wuyd8l+XC6PvYzSS5JcusztLPegX56zSbLV+ffYZPlMChbHVGF4druaDmTdO4cPXo0SbK0tDTmSm46JnGUJWD7djPk33mT+Xfpp6uXB346Xdeb+7XWrlq7YlU9P13I30jbZP5qu3fZZPld160HwITzMCyAndnNkH9RVd1ugy47M/30Lf30Xyb5vxsE/JslefA29rva7oOr6hYb3JT70H565TbaBmAMFhYWzr4SAJvazZB/IMlPJlk7us79kjwh3VX03+tnn0ryxVV1t9ba1f16leSyJF96rjvtb9h9XZKvTXJpkl9Ys//7JzmabpjO39u4BYBh2upoOUZVARie3Qz5f5Jkrg/Wb8pnxsm/WZIn98NnJt0Qmc9L8paqenmSTyR5ULqA/4dJHrWNfT+l3+fP9w/fujyfGSf/hiRP2uSmYBiM2dnZzMzM5ODBg+MuhQkxiaPlMLl8hsCw7GbIX0kXtn+mn946XReZZ7bWXrO6Umvt+VX1sXRX3Z+Y5Pokb0jypCTfnG2E/Nba3/d/NXh6kkem6yL0gXTj7S+01v5y+4cF02F2dnbcJTAhtjtCzqSMrJO48XYcfIbAsOxmyE/fz/4bt7De8STHN1j0V+m67axfv7bQ5juT/PuzrQcwdMeOHRt3CQCM2c3GXQAAALC7hHwAABgYIR8AAAZmx33yW2unkpy1zzwAALA/XMkHAICB2dXRdQBgN8zNzY27BICpJuQDMHGOHDky7hIAppruOgAAMDBCPgATZzQaZTQajbsMgKmluw4AE2dxcTGJbjsA2+VKPgAADIyQDwAAAyPkAwDAwAj5AAAwMEI+AAAMjJAPAAADYwhNACbO0tLSuEsAmGqu5AMAwMAI+QAAMDBCPgATZ35+PvPz8+MuA2Bq6ZMPwMRZWVkZdwkAU82VfAAAGBghHwAABkbIBwCAgRHyAQBgYIR8AAAYGKPrADBxjhw5Mu4SAKaakA/AxJmbmxt3CQBTTXcdAAAYGCEfgImzsrLigVgAO6C7DgATZ35+PkmytLQ05koAppMr+QAAMDBCPgAADIyQDwAAAyPkAwDAwAj5AAAwMEI+AAAMjCE0AZg4CwsL4y4BYKoJ+QBMnEOHDo27BICpprsOAAAMjJAPwMRZXFzM4uLiuMsAmFpCPgATZzQaZTQajbsMgKkl5AMAwMAI+QAAMDBCPgAADIyQDwAAAyPkAwDAwHgYFgATx8OwAHZGyAdg4iwsLIy7BICpprsOAAAMjJAPAAADI+QDMHGOHj2ao0ePjrsMgKkl5AMAwMAI+QAAMDBCPgAADIyQDwAAAyPkAwDAwAj5AAAwMJ54C8DEmZubG3cJAFNNyAdg4hw5cmTcJQBMNd11AABgYIR8ACbOaDTKaDQadxkAU0t3HQAmzuLiYhLddgC2y5V8AAAYGCEfAAAGRsgHAICBEfIBAGBghHwAABgYIR8AAAbGEJoATJylpaVxlwAw1VzJBwCAgRHyAQBgYIR8ACbO/Px85ufnx10GwNTSJx+AibOysjLuEgCmmiv5AAAwMEI+AAAMjJAPAAADI+QDAMDACPkAADAwRtcBYOIcOXJk3CUATDUhH4CJMzc3N+4SAKaa7joAADAwQj4AE2dlZcUDsQB2QHcdgClz6aWXZnl5+Zy3O3z4cI4dO7YHFe2++fn5JMnS0tKYKwGYTkL+hDt+/HhOnTqVgwcPZnZ2dtzlABNgeXk5J0+eHHcZwBSTL4ZPyJ9wx48fz8mTJ3PJJZf4RwjcyIEDB3L48OGzrre8vJzTp0/vQ0XAtJAvhk/IB5hShw8fzokTJ8663szMjCv/ADcxe3bjbVUdrKpWVcer6l5V9bKq+ueq+mBVvbaq/nW/3udX1a9X1TVV9dGq+suqeui6tu5WVT9ZVW+qqndV1cer6uqqWqqqLz3Lvu9dVS+tqmur6oaqmtmrYwYAgEmwH1fyDyb58yRXJTne//zYJCeq6gFJXp3kA0lemuT8JI9P8qqqundr7R19Gw9J8rQkr0/y8iQfSvLFSR6X5NFV9aDW2ls32Pe9+n3/TZIXJTmv3xcAAAzWfoT8S5I8vbW2sDqjqv5TkmemC+C/m+R7W2s39Mtel+S3kvxQ/0qSUZI7t9Y+uLbhqrpPkjcl+Zkkj9hg3w9O8qzW2k9spdCqumKTRRdsZfu9tLy8nJmZmXGXAUyA7Yyss7rdtHyOXHXVVUkyNfXCtNnu5wjTYz9C/ql0IXytF6YL+bdO8mOrAb+3lOQFST59N1lr7dqNGm6tvbWqRkkeVlW3bK19Yt0q707yjJ2VPxlOnz6tTy2wI9P4OXLttRt+/ANwFvsR8pdba59aN+/qfvo366/Ot9Y+VVXvTnKPtfOr6uuTPCXJ/ZJ8Xj679s9Lcs26eW9trX1sq4W21i7eaH5/hf+irbazF7Y6igYwfNsdLcfnCLDKqFvDtx8h/7POoNbaJ6tqw2W9Tya55eoPVfXUJMeSvC/J65K8I8lHkrQkj0lyn3R/FVjvXTspfJJsdRQNYPi2O1qOzxFglVG3hm/ih9CsqlskuSxdYL+otXbNuuUPOMPmbQ9LA2CPLC4uJknm5ubGXAnAdNqzITR30ecluUOSN28Q8D83Y+5GA8DuG41GGY1G4y4DYGpN/JX8JNem65pzcVV9bmvtQ0lSVbdM8ux0XwIAbnK2OlqOUTQAbnomPuS31m6oquekGyf/r6rqlUluleSh6cbVf33/34M0OzubmZmZHDx4cNylABNmGkfLASaDfDF8Ex/ye/8pyXuSzCV5crobdl+X5OkZyBCZm5mdnR13CcCE2e4IOUbWAVbJF8O3ZyG/tXYqSZ1h+ZmWHVz38yeT/FL/Wm+2f2153wDT7NixY+MuAYAJNw033gIAAOdgWrrrAHATcujQoXGXADDVhHwAJs7CwsK4SwCYarrrAADAwAj5AAAwMEI+ABPn6NGjOXr06LjLAJhaQj4AAAyMkA8AAAMj5AMAwMAI+QAAMDBCPgAADIyQDwAAA+OJtwBMnLm5uXGXADDVhHwAJs6RI0fGXQLAVNNdBwAABkbIB2DijEajjEajcZcBMLV01wFg4iwuLibRbQdgu1zJBwCAgRHyAQBgYIR8AAAYGCEfAAAGRsgHAICBEfIBAGBgDKEJwMRZWloadwkAU82VfAAAGBghHwAABkbIB2DizM/PZ35+ftxlAEwtffIBmDgrKyvjLgFgqrmSDwAAAyPkAwDAwAj5AAAwMEI+AAAMjJAPAAADY3QdACbOkSNHxl0CwFQT8gGYOHNzc+MuAWCq6a4DAAADI+QDMHFWVlY8EAtgB3TXAWDizM/PJ0mWlpbGXAnAdHIlHwAABkbIBwCAgRHyAQBgYIR8AAAYGCEfAAAGRsgHAICBMYQmABNnYWFh3CUATDUhH4CJc+jQoXGXADDVdNcBAICBEfIBmDiLi4tZXFwcdxkAU0vIB2DijEajjEajcZcBMLWEfAAAGBghHwAABkbIBwCAgRHyAQBgYIR8AAAYGA/DAmDieBgWwM4I+QBMnIWFhXGXADDVdNcBAICBEfIBAGBghHwAJs7Ro0dz9OjRcZcBMLWEfAAAGBghHwAABkbIBwCAgRHyAQBgYIR8AAAYGCEfAAAGxhNvAZg4c3Nz4y4BYKoJ+QBMnCNHjoy7BICpprsOAAAMjJAPwMQZjUYZjUbjLgNgaumuA8DEWVxcTKLbDsB2uZIPAAADI+QDAMDACPkAADAwQj4AAAyMkA8AAAMj5AMAwMAYQhOAibO0tDTuEgCmmiv5AAAwMEI+AAAMjJAPwMSZn5/P/Pz8uMsAmFr65AMwcVZWVsZdAsBUcyUfAAAGRsgHAICBmdjuOlX1g0mekuRQktsk+aEkv5zkZGttZoylAVPi0ksvzfLy8jlvd/jw4Rw7dmwPKgKA/TGRIb+qHp/k2UnekuRYko8l+bOxFsXUOX78eE6dOpWDBw9mdnZ23OUwBsvLyzl58uS4y2CK+NwAhmIiQ36Sb1idttauXp1ZVRcm+ch4SmLaHD9+PCdPnswll1zil/VN3IEDB3L48OGzrre8vJzTp0/vQ0VMKp8bwFBMasi/W5KsDfj9z389nnKAaXb48OGcOHHirOvNzMy48j8hjhw5Mu4SAKbaRN14W1WXVVVL8tD+57b6WvPziY22qaqZqnpcVf1FVX2kqq6rqpdU1d33/0gA2Im5ubnMzc2NuwyAqTVpV/JP9NPZJF+U5BnnsO33Jnl0kj9IcjLJ/ZN8W5L7VNXh1trHdq9MAACYXBMV8ltrJ5KcqKqZJF/UWrvsHDZ/eJIvb6391eqMqlpK8u1JvjHJ756tgaq6YpNFF5xDHUyY5eXlzMzMjLsMxmA7I+usbuecGa/rr78+SXLeeeft6363e84ATJqJCvk79Jy1Ab/3G+lC/ldkCyGfYTp9+rR+1pwT5wwA025IIf/yDeb9Yz+941YaaK1dvNH8/gr/RdusizHb6sgqDM92R8txzozfVVddlSS58MIL93W/RlgChmJIIf/9G8z7ZD+9+X4WwmTZ6sgqDM92R8txzozf0aNHkyRLS0v7ul8jLAFDMVGj6wAAADsn5AMAwMAMqbsOwIa2OlqOkVUAGAohn8GanZ3NzMxMDh48OO5SGDOj5bBv4uxKAAATYklEQVRVPjeAoRDyGazZ2dlxl8CYbXeEHCPrjN/CwsJY9utzAxiKiQz5rbWZTebXBvMuS3LZJuufSvJZ2wA3DceOHRt3CWzToUOHxl0CwFRz4y0AAAyMkA/AxFlcXMzi4uK4ywCYWkI+ABNnNBplNBqNuwyAqSXkAwDAwAj5AAAwMEI+AAAMjJAPAAADI+QDAMDATOTDsAC4afMwLICdEfIBmDgLCwvjLgFgqlVrbdw1TLyq+ufzzjvv/AsvvHDcpQAAMGBXXXVVrr/++utaa3faSTtC/hZU1UqS2yc5NeZShuKCfvrXY62CSeKcYD3nBOs5J1hvqOfEwSQfaK3tqN+ikM++q6orkqS1dvG4a2EyOCdYzznBes4J1nNOnJnRdQAAYGCEfAAAGBghHwAABkbIBwCAgRHyAQBgYIyuAwAAA+NKPgAADIyQDwAAAyPkAwDAwAj5AAAwMEI+AAAMjJAPAAADI+QDAMDACPnsmaq6ZVU9tar+W1UtV9XHq6pV1dwO2nxgVf3Pqrquqq6vqv9dVZdW1c13s3b21m69j/35tNnrz/aqfranqu5RVS+oqqur6mNVdaqqjlXVHc+xnfP77U717Vzdt3uPvaqdvbEb50RVnTjLZ8Ft9vIY2B1V9biq+pWqekNVfaB/735nm23tymfNtLvFuAtg0D4nybH+v9+d5F1J7rndxqrqG5O8PMlHk7w0yXVJHpXkl5M8KMm37KRY9scevI//kOT4BvP/aftVstuq6l5J3pzkC5K8MslfJ/mKJE9N8vCqelBr7Z+30M6d+nbunWSU5CVJLkjypCRfX1UPaK39/d4cBbtpt86JNZ6xyfxP7qhQ9svTk9wnyYfSfX5fsJ1G9uC8ml6tNS+vPXkluVWSRyS5a//zZUlakrlttHX7JNcm+ViS+62Zf5t0/5hbkseP+5i99vd97Nc/Me7j8trSe/Wa/v36gXXzf6mf/7wttvP8fv1fXDf/B/v5rx73sXrt+zlxoosz4z8mrx2dDw9N8sVJKslMfw78zjba2ZXzaggv3XXYM621j7fWXtVau2YXmntcks9P8pLW2uVr9vHRdN/+k+Tf78J+2Fvex5ug/sraw5KcSvLcdYt/KsmHk3xnVX3OWdr53CTf2a9/2brFv5rurzpfV1X/YudVs5d265xgOFprr2+tvb31iXw7nFc3JuQzLY7001dvsOxPknwkyQOr6tb7VxLbsBfv4x2q6ruq6ieq6vuq6it3XCW77aH99LWttRvWLmitfTDJm5LcNsnZ3ruvTHJekjf1261t54Z0V/DW7o/JtVvnxKdV1bdV1dOq6oer6hF+H9wk7fp5Nc2EfKbFl/TTv1m/oLX2ySQr6e4xcQVvsu3F+3ifJL+ZZCHd1dw/7W/0/rId1sru2fR97729n957n9ph/PbivXxJkmcl+cUk/zPJO6rqcdsrjynlM2INIZ9pcaCfnt5k+er8O+xDLWzfbr+Pv5TuZt3PT3K7JF+e5GXpgv+oqu6+zTrZXbv1vvscGI7dfC9fme7m/Xuk+0vPBenC/h2SvLSqHr6DOpkuPiPWEPI5o37YqTMNTbb+ta3hrpgek3ROtNZ+pLX25tbae1trH2qtXd5a+5Z0o/d8XpIf3at9A5OhtfbLrbU/aq29s7X20dba21prP5HkR9LlnGeNuUQYC0NocjZ/l26ow626eo/qWP32fWCT5avz379H++czdnJO7Nf7+Lwk35zkITtsh92xW++7z4Hh2I/3cjHd0LyHq+p26+/jYJB8Rqwh5HNGrbWvHncNvbcluV+6fnRXrF1QVbdIcijdWMjGx95jOzwn9ut9fE8/vUmMoDAF3tZPN+sH+8X9dLN+tLvdDuO35+9la+2jVfXBJHdM91kg5A+fz4g1dNdhWoz66UZ9Kx+S7m75N7fWPrZ/JbEN+/U+ro6c4EvfZHh9P31YVd3o905V3S7dfRUfSXK2pxT/WZLrkzyo325tOzdLN3Te2v0xuXbrnNhUVX1JuoD/wSTv3W47TJU9P6+miZDPRKmqA1V1QVXddd2il6X7kH58Vd1vzfq3SfKf+x//6z6Vyfad8/tYVbftz4kvXDf/31TVLdfvoKr+TbqRdpLEPSIToLX2d0lem+Rgku9bt/gZ6a6y/nZr7cOrM/v3/EZPvGytfSjJb/frX7aune/v239N88Tbibdb50RVHaqq89e3X1Wfn+S/9T++pB+9i4Goqlv258O91s7fznk1ZLWDZw7AWVXV0/KZR1MfTjfqyZvzmWGs3thaW1yz/my6D+YXttZm17X1mHQh8aPphkq7Lsmj0w2Z9bIk37qTh2iwP871fayqmXRXZ0621mbWzD+ebkSNNyT5x3RP0b0g3V8Jbp7kN5I82TkxGTZ41PxVSe6fblzrv0nywLbmUfNV1T3SuLVa186d+nbune4vQ3+R5MIk35juacoP7H/RM+F245zof2c8L8kb0/3l7rokX5jkken6X1+e5GtbazeJPtjTrP/d8Jj+x7sk+bp07+kb+nnvba39aL/uwXRDLv9Da+3gunbO6bwatHE/ctdr2K/0jxs/w+v4uvVnN5q/ZvmD0o1//L50f7b/qyQ/lOTm4z5Wr3M6L7b8PuYzjzc/sW7+Y5K8IsnfJvlAko8nuSbJHyZ59LiP0WvD9/2e6b7EX9O/X/+Q5FiSO26wbut+RW3YzvlJnt1vv/q+vyDJPcZ9jF77e04k+bIkx/vPkH9O8ol0Qf8NSX4gya3GfYxeWz4XLjtLXji1Zt2D6+dt97wa8suVfAAAGBh98gEAYGCEfAAAGBghHwAABkbIBwCAgRHyAQBgYIR8AAAYGCEfAAAGRsgHAICBEfIBAGBghHwAABgYIR8AAAZGyAdYp6oOVlWrquNjrOHeVfXxqvoP46qhr6Oq6q1V9YYdtDERxzItquruVXV9Vf3ncdcCTC8hH2Ay/VKSf07yq2tnVtXx/gvI7GYbVtVl/TqXbTJ/7eujVfW3VfXrVXVwfVuttZbkJ5M8uKoet5vH0td0UVX996p6d/9F4B1V9WtVdedNjq2q6nuq6s+r6kNV9eGquryqnlJVn/U7rapu37f3T1X1z1X1h1V1r03anquqT1TVfbd5nKmqe1bVz1TVFVX1vr69a6vqj6vqqVV1YN36p/r34eDqvNbaO5M8L8kPV9U9t1sLcNMm5ANMmKp6YJKvT/IrrbWP7MEuTiZ5Rv/6jSQfS/I9Sa6sqi9ev3Jr7ZVJrkqyUFV1Ljs607FU1Tck+bMk39xPn53krUmekuTyqvrCDZr8nSS/nuRgkhcnWUxy2yT/NcnxDdY/nuTfJfnjJC9N8tVJ/ldV3XZdLXdP8gtJfra19pZzOcY1bcwleXuS/5jk5n19P5fkFUnukuRYkr/bYnM/n+RWSf7TdmoBuMW4CwDgs3xfkhuS/NYetX+itXbZ6g/9FfA/TPLIJD+R5EkbbPPCJD+TLiT/8Tnsa8NjqarbpAvot0zyza21V6xZ9u1JltJd+X/0mvmPTXI0yUqSr2itvbeff6skL0/ynVX1+6tt9X8NeGySn2qtPbOf9+fpgv83JPndNSU9L8k7kzzzHI5t7fE8Id0Xpvf1x/M/NljnQUmeu5X2WmtXV9Xrkhytqh9rrZ3eTl3ATZcr+QDnoKruWlXP7btZfLyq3lNVr6iqizdZ/0BVHeu7i3y0qv66qn64qv7FRv3+q+r2SR6X5M2ttX/ah0NKa+2GfOYq+JdvstpL+ul3b7XdsxzLA5PcOcnlawN+X8+L013R/4aq+qI1ix7bT39xNeD36388n7ni/f1r1l/d9i/WzPuLdctSVd+R7gvOd/VtnZOqul2S5/Q/Pn6jgN/X+aYk9z+Hpl+S5HOSPP5cawIQ8gG2qKoOJbk8yfem63bxi0lek647ypv77idr179NklGSpya5Nl13lBNJ5vttN/KQdN003rj7R7Aln9hoZmvtH9Jd6f6ac+iyc6ZjuUs//ftNtv37JJXkyBa3WZ33Vf2V/SR5Rz9d+wXsfv30H5JPX+0/luSXW2t/vkktZ/O4JOcn+bPW2mvPtGJr7WPn0O6b+unXbrMu4CZMdx2ArXtekrsleXprbWF1ZlX9WpI/SfLCqvqi1tqH+kU/luSidFdkj/Y3saaqFpJcuck+HtxPLz9LLY/Z6EbZ3sxZtr2Rqrp5PnOF/kxfLv4yyWOSXJjk/22h6TMdy+qV+EObbPsv+umXbHGb1fVv0f/3X7fW3lVVr0zyU/3Nth9N8sR04X/1avtzk1yXnfV9Xz3O/7WDNj5La+1vq+r96b4sAZwTIR9gC6rqHkkeli4g/tzaZa21N1fVi5N8R5Jvymf6nz8xXX/0H18N+P36/1hVx5JsNETi6s2m15ylpG/sX9sxs2bknfPTXSm+IF1w/+kzbPeuNTVuJeSf6VjelOT9Sb68qr6xv7k3SVJV35rkPv2Pd1yzzf9I8u3pRp15SWvtun79W6a7iTgbbPPEdO/Xo5LcJt1fUi5trX24Hy3om5JckuSGqvqVJE9I8rlJ3pzke1trWznOu/bTvehe9a4kF1TVbVprH92D9oGB0l0HYGtWh1V8Q2ttoy4to7Xr9f3R75Xkna21Uxusv9kV8zv10/edpZ4ntdZqo1duHHg3ckmSn+pfP5Au4C8neWBr7dozbHddP/28s7S/atNjaa19OF03ppbkFVX1e1X1c1X1B+n+8rHcr3rDms1ekq571L2S/L+qen5VPbtf96vyme45n96mtXa6tfbk1trdWmvnt9a+vrX29qo6P92Nvb/WWntDkp9NNwrPZelu9r1Tklf3Xa7G6Vz/nwMkEfIBtmp1fPPNrrCvzr9DP719P333JutvNv/6frqX4fIZ/ZeBm6e72v6cJIeT/O5GY82vcV4/vf4M66x1xmNprf1WutF6Xpeui9FT+3pmk7yoX+3aNet/Kt0V+acleU+6q/RPTDds5QOTfHD9NmfwnL6+p1XV5yT590l+u7X2nNbaq9Pdd3HPdKP5nM3qe3/3Lax7rs71/zlAEiEfYKtWhzC8yybL77puvQ/00w0f6nSG+asB9U6bLN81rbUbWmv/2Fp7apKXpeuO9P1n2GS1pq2E6LXrbXosrbXXt9Ye3lq7Y2vt1q21w334/zf9Kn+5bv1PtNZ+trX2Za2127TW7tBae0ySU0m+OMl7W2srZyqqqr4+Xbec7+nvn7hXuhuE194ncUU//VdbOM7Vv8p89RbWPVd3SvLJfOaKPsCWCPkAW7P6gKQHV9VG9zM9tJ9emSSttQ+kG/Hl7pvcIPvgDeYlyf/upxdsr8xt+5F0D8X6yb6r0UYuSNcV5q+22Oa2jqWq7pDuiv170l3l34rHpwvqLz5L2weSPD/Jb7bW1o/3f+s1/30uf0l5WboQ/oCq+pqz7P/WZ1q+bt3PTffXgf+99p4OgK0Q8gG2oB/n/XXpnrR66dplVXX/dN063pfk99Ys+q10n7PPWjvsZFXdc30ba5zop1+5G3VvVWvtHeke5nSndIH/RvpwejjJW1pr799isyf66YbH0o8vv37ebdM9eOsOSX5y/ZCTG30BqarD6Z4Q+750D+w6k9WhS9ce498l+Xi6B2StelQ//b9naS+ttQ8m+cH+x5dW1ddttF5VfWWSPz1be2t8ebouVa8/h20AkhhdB+BcPCXdqDA/X1UPSzc05D2TfEu6K9xP6gPfqp9LN+Tk45N8SVW9Nl3f/m9NN+TmY3LjG0vTWvs/VfW2JF9dVTfv+6Hvl/+SbijNH6qqX1n7wKl0feZXnyy7JVs4lidW1Y+k+zJwTbovGI9K1/Xp2a21523Q7Ouq6vok/yddH/wL0z2n4Pokj2qtXb1ZPf1V9u/u1/v0E2T7kXae2x/3q5P8bbqn/v5juifvbuVYX1RV56W7mffVVbWcboSe9/XH9YB0Iwa9d/NWPsvD+umW/58DrHIlH2CLWmt/n+5hSs9LN377jyZ5RJJXJ3nQ2mEg+/WvT9eN51fS9eX/of7n/5LkWf1qH8hn+6/9+g/bYNmeaa1d0+/7dkl+fN3iJ6a72v2b59jsmY7l8iRXJXl4uv+X35xupJxHttY2+0vHy/r6viPJD6fru//rSb60tXZysyL6ri+/keRFrbU/2mCVH0/3sLKLk8wl+fMkDz+XYStba4vp7gv4uXSjBj0hyX9M97Cs96V7/++1lbb6G6C/I8lbW2vncvUfIElSuvkB7L+q+p504fQprbXnr1t2+3RdSN7cWtvuWPi7pqq+IN2NrUuttblz3HaijmVaVNWjkvxBku9srf3OuOsBpo+QD7CHqupu67uQVNUXphuR5a5JvmijLiZV9b3pnsZ6v9baFeuX76f+wV3fneTe/dX+c91+Yo5lGvT3b1yR5FNJvsJNt8B26JMPsLde3j+R9Yp0T3g9mO4Gz9umexLuZn3In5/u5tPNhuzcF33gvCbdFeVzDvi9iTiWKXKXdFfxf1/AB7bLlXyAPdRfxf7OdH21DyT5ULrhOH+1tfaKcdYGwHAJ+QAAMDBG1wEAgIER8gEAYGCEfAAAGBghHwAABkbIBwCAgRHyAQBgYIR8AAAYGCEfAAAGRsgHAICBEfIBAGBghHwAABgYIR8AAAZGyAcAgIH5/8PrsTqZFJ6IAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 263, "width": 380 }, "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "cph.print_summary(3, model=\"stratified age and wexp\")\n", "cph.plot()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Proportional hazard assumption looks okay.\n" ] } ], "source": [ "cph.check_assumptions(rossi_strata_age)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Option 2: introduce time-varying covariates\n", "\n", "Our second option to correct variables that violate the proportional hazard assumption is to model the time-varying component directly. This is done in two steps. The first is to transform your dataset into _episodic format_. This means that we split a subject from a single row into $n$ new rows, and each new row represents some time period for the subject. It's okay that the variables are static over this new time periods - we'll introduce some time-varying covariates later.\n", "\n", "See below for how to do this in _lifelines_:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
stopstartarrestagefinidmarparoprioracewexp
01.00.00270001310
12.01.00270001310
23.02.00270001310
34.03.00270001310
45.04.00270001310
56.05.00270001310
67.06.00270001310
78.07.00270001310
89.08.00270001310
910.09.00270001310
1011.010.00270001310
1112.011.00270001310
1213.012.00270001310
1314.013.00270001310
1415.014.00270001310
1516.015.00270001310
1617.016.00270001310
1718.017.00270001310
1819.018.00270001310
1920.019.01270001310
201.00.00180101810
212.01.00180101810
223.02.00180101810
234.03.00180101810
245.04.00180101810
\n", "
" ], "text/plain": [ " stop start arrest age fin id mar paro prio race wexp\n", "0 1.0 0.0 0 27 0 0 0 1 3 1 0\n", "1 2.0 1.0 0 27 0 0 0 1 3 1 0\n", "2 3.0 2.0 0 27 0 0 0 1 3 1 0\n", "3 4.0 3.0 0 27 0 0 0 1 3 1 0\n", "4 5.0 4.0 0 27 0 0 0 1 3 1 0\n", "5 6.0 5.0 0 27 0 0 0 1 3 1 0\n", "6 7.0 6.0 0 27 0 0 0 1 3 1 0\n", "7 8.0 7.0 0 27 0 0 0 1 3 1 0\n", "8 9.0 8.0 0 27 0 0 0 1 3 1 0\n", "9 10.0 9.0 0 27 0 0 0 1 3 1 0\n", "10 11.0 10.0 0 27 0 0 0 1 3 1 0\n", "11 12.0 11.0 0 27 0 0 0 1 3 1 0\n", "12 13.0 12.0 0 27 0 0 0 1 3 1 0\n", "13 14.0 13.0 0 27 0 0 0 1 3 1 0\n", "14 15.0 14.0 0 27 0 0 0 1 3 1 0\n", "15 16.0 15.0 0 27 0 0 0 1 3 1 0\n", "16 17.0 16.0 0 27 0 0 0 1 3 1 0\n", "17 18.0 17.0 0 27 0 0 0 1 3 1 0\n", "18 19.0 18.0 0 27 0 0 0 1 3 1 0\n", "19 20.0 19.0 1 27 0 0 0 1 3 1 0\n", "20 1.0 0.0 0 18 0 1 0 1 8 1 0\n", "21 2.0 1.0 0 18 0 1 0 1 8 1 0\n", "22 3.0 2.0 0 18 0 1 0 1 8 1 0\n", "23 4.0 3.0 0 18 0 1 0 1 8 1 0\n", "24 5.0 4.0 0 18 0 1 0 1 8 1 0" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from lifelines.utils import to_episodic_format\n", "\n", "# the time_gaps parameter specifies how large or small you want the periods to be. \n", "rossi_long = to_episodic_format(rossi, duration_col='week', event_col='arrest', time_gaps=1.)\n", "rossi_long.head(25)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each subject is given a new id (but can be specified as well if already provided in the dataframe). This id is used to track subjects over time. Notice the `arrest` col is 0 for all periods prior to their (possible) event as well. \n", "\n", "Above I mentioned there were two steps to correct `age`. The first was to convert to a episodic format. The second is to create an interaction term between `age` and `stop`. This is a time-varying variable.\n", "\n", "Instead of `CoxPHFitter`, we must use `CoxTimeVaryingFitter` instead since we are working with a episodic dataset. " ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "rossi_long['time*age'] = rossi_long['age'] * rossi_long['stop']" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from lifelines import CoxTimeVaryingFitter\n", "ctv = CoxTimeVaryingFitter()\n", "\n", "ctv.fit(rossi_long, \n", " id_col='id', \n", " event_col='arrest', \n", " start_col='start', \n", " stop_col='stop', \n", " strata=['wexp'])" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " event col = 'arrest'\n", " strata = ['wexp']\n", "number of subjects = 432\n", " number of periods = 19809\n", " number of events = 114\n", " log-likelihood = -575.080\n", " time fit was run = 2019-02-19 17:30:53 UTC\n", "\n", "---\n", " coef exp(coef) se(coef) z p -log2(p) lower 0.95 upper 0.95\n", "age 0.073 1.075 0.040 1.830 0.067 3.893 -0.005 0.151\n", "fin -0.386 0.680 0.191 -2.018 0.044 4.520 -0.760 -0.011\n", "mar -0.397 0.672 0.382 -1.039 0.299 1.743 -1.147 0.352\n", "paro -0.098 0.907 0.196 -0.501 0.616 0.698 -0.481 0.285\n", "prio 0.090 1.094 0.029 3.152 0.002 9.267 0.034 0.146\n", "race 0.295 1.343 0.308 0.955 0.340 1.558 -0.310 0.899\n", "time*age -0.005 0.995 0.002 -3.337 0.001 10.203 -0.008 -0.002\n", "---\n", "Log-likelihood ratio test = -1150.160 on 7 df, -log2(p)=-0.000\n" ] } ], "source": [ "ctv.print_summary(3)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzQAAAIPCAYAAACluJutAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xt8LVV99/HPT0HAqoeLiiJiUity+liJ52gteDkh9W6901qj1mhjpV5PbX0ea6g9XvLSPm01WmtBowaVeClotT5Fsd0kIhSRA9tai/dERQEvSEC5w3r+mIlsNjvZyU52Zq/k83699mtOZtas+e1kkpxvZs2aSCkhSZIkSTm6Q9UFSJIkSVKnDDSSJEmSsmWgkSRJkpQtA40kSZKkbBloJEmSJGXLQCNJkiQpWwYaSZIkSdky0EiSJEnKloFGkiRJUrYMNJIkSZKyZaCRJEmSlC0DjSRJkqRsGWgkSZIkZctAI0mSJClbBhpJkiRJ2TLQSJIkScrWPlUXoO6KiDngbsB8xaVIkiRpc+sDrkop9W/kQQ00m9/dDjjggIO3b99+cNWFSJI6Nzc3B0B//4b+P0GSVuziiy/m2muv3fDjGmg2v/nt27cfvHfv3qrrkCStwfDwMADT09MVVyJJre3cuZMLL7xwfqOP6z00kiRJkrJloJEkSZKULQONJEmSpGwZaCRJkiRly0kBJEnKwOjoaNUlSFJPMtBIkpSBoaGhqkuQpJ7kkDNJkiRJ2TLQSJKUgVqtRq1Wq7oMSeo5DjmTJCkDk5OTgEPPJKmZV2gkSZIkZcsrNJIkacvYvXs39Xp91fsNDAwwMTHRhYokrZWBRpLUkampKebn5+nr62NkZKTqcqQVqdfrzM7OVl2G1FVb7eezgUaS1JGpqSlmZ2fZtWvXlviFqc1l27ZtDAwMtG1Xr9dZWFjYgIqk9bPVfj4baCRJ0pYzMDDAzMxM23aDg4Ne0ZF6nJMCABHRFxEpIqYi4siI+FhE/DgibomIwYjYGRHviIivRMQVEXFdRHwrIv4+Ig5apt9nR8R/NOwzHxEfiYiHtmj7nIg4KyKuLNteHBEnRsR+3X33kiRJUr68QnNb9we+BHwTOBU4ALgK+BPgGcAs8O8UQXAn8GrgiRHx8JTS1YudREQAHwBeAPwU+ATwE+Bw4DjgG8AFDe3fD7wQuAQ4HbgS+B3gTcDvRsRjU0o3de1dS5J63vT0dNUlSFJPMtDc1iOBt6SUXte4MiLeArwspXRz0/o/BiaBlwJ/07DpxRRh5svAY1NKCw373BG4Z8PHIxRh5pPAc1NK1zZs2wP8NfAy4B3LFR4Re5fYdNRy+0nSWtXrdQYHB6suQ1qRTmY4W9zP81y56PQ8z5WB5rYuB97QvDKl9L0l2r8feBvweG4baF5RLl/SGGbKvm4GLm1Y9SrgJuBFjWGm9Cbg5cBzaRNoJKkqCwsL3mOgTc/zXOpdBprb+kpK6frmlRGxL/AS4A+B3wS2cdv7j+7T0PbXgAcBl6eULlruYBFxZ+BoimFpu4uRardzPbC9XeEppZ1LHGMvsKPd/pLUqZXOFqW1mZubA6C/v7/iSvLW6axlnufKyVabnc9Ac1uXLbH+YxT30HwX+FTZbjH47AYab9w/sFz+cAXHOwgI4B4UQ8skKTsrnS1KazM8PAx4L81adTprmee5crLVZucz0NxWal5Rzkj2DIrJAJ7YeHN+RNwB+N9Nu1xZLu9De4vR+aKUkldRJEmSpFVy2ub2fqNcfrrFTGO/TTET2q+klH4J/DdwaEQ8ZLmOU0q/AL4G/K+IOHid6pUkSZK2DK/QtDdfLgeBf1hcGRH3BP5xiX3eCbwHOLmccrlxlrM7AIemlBYnBngb8D7g/RExklK6srGj8jk3/SmlC9fhvUiSJFY+a9lWmy1KypGBpr0vA+cAz4yIc4EvAocCT6R4nsyPWuwzCTwKeD7wrYj4FMVzaA4DhihmR9sDkFJ6f0TspJj6+TsR8Tng+8DBQD/waIpn2pzQpfcnSR0ZGRlhcHCQvr6+qkuRVs1Zy7SZbbWfzwaaNlJKN0fEU4E3A08CXklxw/9kue5/WuyTgD8qw8mfAH9AMXHApcDZwKeb2r8sIs6gCC2PoZhY4AqKYPO3wIe78uYkaQ1GRkaqLkFatU5nKnOGM+Vkq/18NtAAKaV5itnGltp+BcUVlFb6ltnvVODUFdbwGeAzK2krSdp6hoaGqi5hU5iYmKi6BEnrzEAjSVIGRkdHqy5BknqSs5xJkiRJypaBRpKkDMzNzTE3N1d1GZLUcxxyJklSBsbGxgCYnp6uuBJJ6i1eoZEkSZKULQONJEmSpGwZaCRJkiRly0AjSZIkKVsGGkmSJEnZMtBIkiRJypbTNkuSlIHx8fGqS5CknmSgkSQpA/39/VWXIEk9ySFnkiRJkrJloJEkKQOTk5NMTk5WXYYk9RwDjSRJGajVatRqtarLkKSeY6CRJEmSlC0DjSRJkqRsGWgkSZIkZctAI0mSJClbBhpJkiRJ2fLBmpIkZcAHa0pSawYaSZIyMD4+XnUJktSTHHImSZIkKVsGGkmSJEnZMtBIkpSB4eFhhoeHqy5DknqOgUaSJElStgw0kiRJkrJloJEkSZKULQONJEmSpGwZaCRJkiRly0AjSZIkKVv7VF2AJElqb3R0tOoSJKknGWgkScrA0NBQ1SVIUk9yyJkkSZKkbBloJEnKQK1Wo1arVV2GJPUch5xJkpSByclJwKFnktTMKzSSJEmSsmWgkSRJkpQtA40kSZKkbBloJEmSJGXLQCNJkiQpWwYaSZIkSdly2mZJkjIwPT1ddQmS1JO8QiNJkiQpWwYaSZIkSdky0EiSlIGxsTHGxsaqLkOSeo730EiSlIG5ubmqS5CknuQVGkmSJEnZMtBIkiRJypaBRpIkSVK2vIdmA0VEAmZTSoNV1yJJ0u7du6nX66veb2BggImJiS5UJEmrZ6CRJGkTmZqaYn5+nr6+PkZGRpZtW6/XmZ2drbwOSVoLA83G2g5cU3URkqT8DA0Nrajd1NQUs7Oz7Nq1a8VBYtu2bQwMDLRtV6/XWVhY6FodktQJA80GSil9veoaJEl5Gh0d7VrfAwMDzMzMtG03ODjYtSs6ktQpJwVYQkT0RUSKiKmIOCoi/iUiroiIX0bEFyPicU3tR8r2IxHxhIiYiYiF8r6ZxTYpImZaHGtbRLwlIr4REddFxM8j4nMR8ZgNeKuSJElStgw07fUD/wkcDJwM/DOwEzgjIp7dov3xwGeAq4GTgI8t13lEHAicC7wWWAAmgNOBY4AzI+Il6/M2JEk5m5ub8+GaktSCQ87aezTwdyml1yyuiIh3UYSckyLijJTSVQ3tnwQ8KaX02RX2/zfAbwLvAU5IKaXyGH8DXAC8MyI+l1KaX66TiNi7xKajVliHJKmHjY2NATA9Pb2i9vV6ncHBwbZtOtHNviVptQw07S0Ab2xckVK6ICJOBV4APAM4pWHzp1YaZiLiTsDzgF8Af7kYZspjfCsi3gmcCPxRcw2SJC1nYWGha/e7dLNvSVotA017F6aUrm6xfoYi0DyE2waa81fR9wOBOwPnpJSuaLG9RhFoHtKuo5TSzlbryys3O1ZRkyRpE1jJzGWrmbVso/qWpNUy0LR3+RLrLyuX25ZYvxKL+166xPbF9Qeuok9JklY0c1mns5Z1s29JWi0nBWjv0CXW36tcNv/5KTU3XMbivvdaYvu9lziGJEmSJAw0K7EjIu7aYv1gubxoDX1/g+JBm0eXs501O65cXriGY0iSJEmblkPO2tsGvB5onOXsocBzKa6cfLLTjlNKN5STC7wYeBPwioZj3B94JXAj8KFOjyFJUjsrmbVssZ0k9RoDTXtfAEYj4uHAORTDwJ5NcXXrJU1TNnfitcCjgJdHxMOAs4C7A38A3BV4eUrJBw9I0hY3Pj6+onYjIyMMDg7S19e34r67MWtZJ3VIUicMNO3NAScAby2X+1EMAXtjSulza+08pXRFRBwD/CXwTODVwLUUs6X9bUrpzLUeQ5KUv/7+/hW1GxkZWXGf7WYqW8t+q6lDktbCQLMCKaWLgae1aTMFTLVpE0usvxL4P+VLkqQNMTExUXUJkrRmTgogSVIGJicnmZycrLoMSeo5BhpJkjJQq9Wo1WpVlyFJPcdAI0mSJClb3kOzhJTSPNDynhdJkiRJvcErNJIkSZKyZaCRJEmSlC0DjSRJkqRseQ+NJEkZWOmDNSVpqzHQSJKUgfHx8apLkKSe5JAzSZIkSdky0EiSJEnKloFGkqQMDA8PMzw8XHUZktRzDDSSJEmSsmWgkSRJkpQtA40kSZKkbBloJEmSJGXLQCNJkiQpWwYaSZIkSdnap+oCJElSe6Ojo1WXIEk9yUAjSVIGhoaGqi5BknqSQ84kSZIkZctAI0lSBmq1GrVareoyJKnnOORMkqQMTE5OAg49k6RmXqGRJEmSlC0DjSRJkqRsGWgkSZIkZctAI0mSJClbBhpJkiRJ2TLQSJIkScqW0zZLkpSB6enpqkuQpJ7kFRpJkiRJ2TLQSJIkScqWgUaSpAyMjY0xNjZWdRmS1HO8h0aSpAzMzc1VXYIk9SSv0EiSJEnKloFGkiRJUrYMNJIkSZKyZaCRJEmSlC0DjSRJkqRsOcuZJEkZGBoaqroESepJBhpJkjIwOjpadQmS1JMcciZJkiQpWwYaSZIyMDc358M1JakFh5xJkpSBsbExAKanpyuuRJJ6i1doJEmSJGXLQCNJkiQpWwYaSZIkSdky0EiSJEnKloFGkiRJUrYMNJIkSZKy5bTNkiRlYHx8fE377969m3q9vur9BgYGmJiYWNOxJambDDSSJPW4qakp5ufn6evro7+/v6M+6vU6s7Oz61zZ8hrrHhkZ2dBjS9o6DDSSJPW4qakpZmdn2bVr15qDwbZt2xgYGGjbrl6vs7CwsKZjrWfdkrQUA40kSZm49NJL19zHwMAAMzMzbdsNDg5u+BUdSeqEkwK0EREjEXF6RHw3Iq6NiKsi4pyIeN4S7R8WEWdGxNVl23+PiGMiYk9EpIgYbLHPURExFRE/iIgbIuLyiJiOiAd2/Q1KkrJx5ZVXVl2CJPUcr9C090/A14AvAJcChwBPAj4UEQ9MKf3VYsOIeDRwJnBH4BPAd4DfAs4Caq06j4gnlG33Bf4V+DZwOPBM4MkRcVxK6cLuvDVJkiQpbwaa9h6UUvpO44qIuBNwBvDaiDgppfTDiLgD8D5gP+BJKaUzGtqfQBGMaOrnIOAjwDXAo1NK/9Ow7UHAecAksKNdkRGxd4lNR7XbV5KUh5///OcMDg52tG8nM5wt7rfRx5Sk1TDQtNEcZsp1N0TEPwJDwO8CHwSOBX4DOKsxzJTeA/wZcGTT+j8CDgRe3hhmymP8d0S8F9gdEb/ZvF2StPXceOONG35fy8LCgvfSSOppBpo2IuII4P9QBJcjgAOamtynXD6kXH6xuY+U0i0RcS63DzTHlMujI2JPi8Mvtt8OLBtoUko7l6h/Lyu4wiNJ6n377rsvxx57bEf7djpr2UpnRVvPY0rSahholhERvw6cDxwEnE1xf8wCcDPQB7yAYogZwLZyefkS3bVaf0i5fHGbUu6ysoolSZvZQQcdtKIZylrpdNaylc6Ktp7HlKTVMNAs79UUoeOFKaWpxg0R8RyKQLPoqnJ56BJ9tVq/+Gero1NK/7WGOiVJW8D+++9fdQmS1HOctnl5v1EuT2+xbVfTxxeVy0c2NywnDGg1RuC8cvmojqqTJG0p/f39VZcgST3HKzTLmy+XgxRTKgMQEY8HRpvankMxTfNxEfHEpokB/oTb3z8D8AFgDPjriPhySun8xo1lEHp0SmlmDe9BkqRfWemsZc5QJikXBprlvRt4IfDPEXEa8CPgQcATgI8Dz15sWN74Pwp8Fvh0RJxOEXAeDDyWYprnJwK3NOzzs4g4HvgkcF5E/AfFM28ScF+KSQMOARxjIElb2MjICIODg/T19a25r42ctWw965akpRholpFS+q+IOA54M/Bkis/XVygeenklDYGmbD8TEbsa2gN8CTgOeG758VVN+/xHRDwY+Avg8RTDz26gCE81Wg93kyRtISMjIwwPD/PNb36TkZGRjvrodKayTvcDOq5VklbDQNNGSulciufNtBIt2n+J4orMbRtG/B3F7GjfarHPPPDyNRUqSdIyJiYmqi5BkrrCSQHWUUTcOSIObLF+hGJSgDNTSr/c8MIkSZKkTcorNOvrCOCiiPg88G2Kz+9DKGY+uxL48wprkyRJkjYdA836uhw4lWJK5+MoHrp5GcVsZuMppe9UWJskSZK06Rho1lFK6efcfjpnSZIkSV3iPTSSJEmSsuUVGkmSMjA66gAASWrFQCNJUgaGhpZ6goAkbW0OOZMkSZKULQONJEkZqNVq1Gq1qsuQpJ7jkDNJkjIwOTkJOPRMkpp5hUaSJElStgw0kiRJkrJloJEkSZKULQONJEmSpGwZaCRJkiRly0AjSZIkKVtO2yxJUgamp6erLkGSepJXaCRJkiRly0AjSZIkKVsGGkmSMjA2NsbY2FjVZUhSz/EeGkmSMjA3N1d1CZLUk7xCI0mSJClbBhpJkiRJ2TLQSJIkScqWgUaSJElStgw0kiRJkrLlLGeSJGVgaGio6hIkqScZaCRJysDo6GjVJUhST3LImSRJkqRsGWgkScrA3NycD9eUpBYcciZJUgbGxsYAmJ6errgSSeotXqGRJEmSlC0DjSRJkqRsGWgkSZIkZctAI0mSJClbBhpJkiRJ2TLQSJIkScqW0zZLkpSB8fHxqkuQpJ5koJEkKQP9/f1VlyBJPckhZ5IkSZKyZaCRJCkDk5OTTE5OVl2GJPUcA40kSRmo1WrUarWqy5CknmOgkSRJkpQtA40kSZKkbBloJEmSJGXLQCNJkiQpWwYaSZIkSdnywZqSJGXAB2tKUmsGGkmSMjA+Pl51CZLUkxxyJkmSJClbXb1CExGDwFnAG1JKe7p5LEmScrd7927q9fqq9xsYGGBiYqILFUlS71tzoImIPmAOOCWlNLLW/jbSYu0ppai4FEmSqNfrzM7OVl2GJGWl2/fQnA9sB37a5eOsWETcF7gkpZRabNsXuHtK6dKNr0ySpMK2bdsYGBho265er7OwsLABFUlS7+pqoEkpXQN8vZvH6MAkcLeIeDnws8WVEfFY4J3ALHBCRbVJksTAwAAzMzNt2w0ODnpFR9KWt6ZJASJiD8VwM4AXRERqeI1ExGD57z1N+82U6/eNiNdHxHci4rqI+EZEvLih3QkR8dWIuDYiLomIN0REy5oj4uERcVpEXBYRN0TEDyLi5Ig4rKnpk4EPAZ8G3lbuexrwbuANwJ829HlkRLw1Ii6IiJ9ExPUR8b2IeE9EHL5EHftFxJ6I+G7Zfi4i3lyuTxEx02KffSLipRFxXkRcFRHXRMRFEfHypd6vJEmSpLVfoZkBDgReBXwF+JeGbfVy23I+Cjwc+DfgRuB44D0RcSPwYOAFwGeA/wCeCrweuAb4m8ZOIuJFwHuA6ymCyg+ABwCjwFMi4ndSSt8HSCndBLy7DDFnl13cCzi6vKLU6JkUV2vOAs4FbgD+V0O/D00p/bChjgBOpwhN3wLeBewLjJT73U45zO1fgccD3wCmgeuA44B/KD8/z1/2syhJkiRtUWsKNCmlmYiYpwg09eaZzMpZzpZzBPCglNKVZfu/pxii9nbgSuDBi4GhvMrzbeAvIuLvy2BCRBwJnATMA7uaAsbvAmcC7wCeUa7bB3gxcCLwJeBI4DLgKxFxIvDxhvtrPgS8PaV0fdP7ehxwRtnHnzZseh5FmDkbeExK6Yay/euB85b4HIxRhJl3AbtTSjeX+9yRIqS9KCJOSyl9arlPZETsXWLTUcvtJ0mSJOWs6uFMr10MMwAppe8CX6S4svOmxnBStvtX4O7AfRr6+FOKqyCvamxf7vMfFFdsnhIRdy1X/z+KKz9PBV5dtjseeCnFkLN/atj/h81hplx/JvA1iiDS6AXl8sTFMNNQ+5ua+ymHk72CIlD92WKYKfe5GfhzIAHPbd5XkiRJUvdnOWvnghbrflQuW11xWAwshwPfK/99TLncFREPa7HPPYE7UlyJ2UsxXOySlFIqp20GIKX0+Yj4LYrABPxqCNlzKYaMHQ0cVPa16AZu6yHALRTD05p9scW6I4GDKYannVgc7naupZgpblkppZ2t1pdXbna021+SJEnKUaWBJqXUaq7Jm8rlctv2bVh3SLl8TZvD3aU85g+WqedGoHHK5rcBu8t1n6MIVNeW20aA+zV1sQ24YnE4XJPLW6xbrP0BwF+3q12SJEnSbVV9hWY9LAafbSmlq1azY0ppHmh5WSQi7gm8Evhv4NiU0tVN25/TYrergIMjYp8WoebQZWr/ZErpmaupXZIkSdL63EOzeN/HHZdt1T2LN9s/ap37/XWKz8+ZLcLM4eX2ZheV+xzbYtsjW6z7OsXkB79TznYmSZIkaRXW4wrNzyluXD9iHfrqxLuAPwHeHhHfSil9s3FjRNwJeHhK6eyWey9tvlw+MiLu2DD72F2A99L6c/dBYAh4c0Q0znK2Dfir5sYppZsi4h/Kbe+MiFenlK5tbBMR9wYOSin9zyrrlyRlql6vMzg4uKJ2krTVrTnQpJR+ERFfAh4VEacC36S4avPptfa9wuN/vXwOzfuBr0XEZ8sa9qUIWY8CfsIqpy9OKV0WER8F/hCoR8SZFPfIPJbiOTF1YKBptw+W7Z8A/HdEfLqs41nAl4EHUkwa0OhNFBMOnEAxG1uN4l6de1LcW/MIiqmdDTSStEUsLCwwOztbdRmSlIX1uofm+RTPjnkC8ByK+1Iu4darHF2VUvpwRHyFYprj44DHAb+kmDHtNOBjHXb9x8B3gWcDL6MIRp+meMDn6S3qSBHxDOB1FJ+TV1BMKHAK8G7g6RT32TTuc2NEPJ3iGTYjwO9RTALwE2CO4urNqR3WL0nKyMBA89/JbnXxxRcDsH377Se+XG4/Sdrs4tZnSKqbIuKxFA/5fGtK6S838Lh7d+zYsWPv3qWeuylJysHw8DAA09PTFVciSa3t3LmTCy+88MKlHifSLVU/WHPTiYjDWqw7BHhr+eEnN7YiSZIkafPaDNM295q3RcTRFA/X/AnFQ0CfSPEAzZNTSudXWZwkSZK0mRho1t8nKJ458xTgQIoJBL4GvK98SZIkSVonBpp1llL6OPDxquuQJG0u3jsjSa15D40kSZKkbBloJEmSJGXLQCNJUgbGxsYYGxurugxJ6jneQyNJUgbm5uaqLkGSepJXaCRJkiRly0AjSZIkKVsGGkmSJEnZMtBIkiRJypaBRpIkSVK2nOVMkqQMDA0NVV2CJPUkA40kSRkYHR2tugRJ6kkOOZMkSZKULQONJEkZmJub8+GaktSCQ84kScrA2NgYANPT0xVXIkm9xSs0kiRJkrJloJEkSZKULQONJEmSpGwZaCRJkiRly0AjSZIkKVsGGkmSJEnZctpmSZIyMD4+XnUJktSTDDSSJGWgv7+/6hIkqSc55EySJElStgw0kiRlYHJyksnJyarLkKSeY6CRJCkDtVqNWq1WdRmS1HMMNJIkSZKyZaCRJEmSlC0DjSRJkqRsGWgkSZIkZctAI0mSJClbPlhTkqQM+GBNSWrNQCNJUgbGx8erLkGSepJDziRJkiRly0AjSZIkKVsGGkmSMjA8PMzw8HDVZUhSzzHQSJIkScqWgUaSJElStgw0kiRJkrJloJEkSZKULQONJEmSpGwZaCRJkiRla5+qC5AkSe2Njo5WXYIk9SQDjSRJGRgaGqq6BEnqSQ45kyRJkpQtA40kSRmo1WrUarWqy5CknuOQM0mSMjA5OQk49EySmhloJElawu7du6nX66veb2BggImJiS5UJElqZqCRpA00NTXF/Pw8fX19jIyMVF2O2qjX68zOzlZdhrrE70dpczDQSNIGmpqaYnZ2ll27dvkfqIxs27aNgYGBtu3q9ToLCwsbUJHWg9+P0uZgoJEkqY2BgQFmZmbathscHPSKjiRtsC0/y1lE9EVEioipiDgqIv4lIq6IiF9GxBcj4nFN7bdFxGsiohYRl0TEDRHxk4j4dEQcs8QxUkTMRMS9ImIyIn4YETdHxEhDm3tHxD9GxHxDn5+IiJ1d/hRIkiRJ2fIKza36gf8EvgqcDNwbeDZwRkQMp5Q+VrbbDowDXwD+H/Bz4AjgqcATI+IpKaXPtuj/YOA84BfAJ4BbgMsBIqIf+CJwGFADPgLcF/h94MkR8ayU0mfW/R1LkiRJmTPQ3OrRwN+llF6zuCIi3kURck6KiDNSSlcBFwOHpZR+2rhzRBwOnA+8HWgVaH4L+BDwopTSTU3bTqIIMyemlMYb+nw3RXA6JSLul1L6xVLFR8TeJTYdtdQ+kqpTr9cZHBysugy10ckMZ4v7devr63mzfjr9+krqLQaaWy0Ab2xckVK6ICJOBV4APAM4JaXU8m7PlNIlEXEa8IqIOCKl9P2mJjcAf9EcZsog9Djg+8D/berz3Ij4CPA84JnABzt+d5J6ysLCgvdabGJ+fSVp4xhobnVhSunqFutnKALNQ4BTACLiEcCrgGOAewJ3atrnPhQBpdF8SunHLfp/SLk8O6V0Y4vtNYpA8xCWCTQppZb32pRXbnYstZ+kaqx01ixVq9NZy/z65sFZ6aTNwUBzq8uXWH9ZudwGEBHPAE4DrgM+D3wH+CXFPTGDwC5gv2X6abatXF66xPbF9QcusV1ShlY6a5aq1emsZd34+o6NjQEwPj7epqVWylnppM3BQHOrQ5dYf69yufgnnDdRDB97aErp4saGEXEyRaBpJS2xfrHfey2x/d5N7SRJW9Dc3FzVJUhST9ry0zY32BERd22xfrBcXlQufwP4nxZh5g7AIzs47mK/j4yIVgHzuHJ5YQd9S5IkSZuaV2hutQ14PdA4y9lDgedSXB35ZLl6HnhARByWUvpR2S6APcBvrvag5WQCnwceC+wG/q7h+A8Hhimmhv5k6x4kSd220lnLnDVLkjaegeZWXwBGyxBxDrc+h+YOwEvKKZuhmJb5JOCiiDgduBEZWcDfAAAdrElEQVR4BEWY+VfgKR0c+4TymH9bPsjzAm59Ds0twAuXmLBAUmZGRkYYHBykr6+v6lK0Cs5atjn5/ShtDgaaW81RBIu3lsv9KIZ5vTGl9LnFRimlkyPieoqrKS8ArgXOBl4IPIsOAk1K6bvl1aATgSdRDHO7iuJ5NuMppS93/rYk9ZKRkZGqS9AqdDpTmTOc5cHvR2lzMNA0KO+LedoK2k0BUy02fZVi6Flz+1hBnz8E/rRdO0nSxpmYmKi6BElSGwYaSZIyMDQ0VHUJktSTDDSSJGVgdHS06hIkqSc5bbMkSZKkbG35KzQppXmg7T0ukiRVafHBmv39/RVXIkm9ZcsHGkmScjA2NgbA9PR0xZVIUm9xyJkkSZKkbBloJEmSJGXLQCNJkiQpWwYaSZIkSdky0EiSJEnKloFGkiRJUractlmSpAyMj49XXYIk9SQDjSRJGfCBmpLUmkPOJEmSJGXLQCNJUgYmJyeZnJysugxJ6jkGGkmSMlCr1ajValWXIUk9x0AjSZIkKVsGGkmSJEnZMtBIkiRJypaBRpIkSVK2DDSSJEmSsuWDNSVJyoAP1pSk1gw0kiRlYHx8vOoSJKknOeRMkiRJUrYMNJIkSZKyZaCRJCkDw8PDDA8PV12GJPUcA40kSZKkbBloJEmSJGXLQCNJkiQpWwYaSZIkSdky0EiSJEnKloFGkiRJUrb2qboASZLU3ujoaNUlSFJPMtBIkpSBoaGhqkuQpJ7kkDNJkiRJ2TLQSJKUgVqtRq1Wq7oMSeo5DjmTJCkDk5OTgEPPJKmZV2gkSZIkZctAI0mSJClbBhpJkiRJ2TLQSJIkScqWgUaSJElStgw0kiRJkrLltM2SJGVgenq66hIkqSd5hUaSJElStgw0kiRJkrJloJEkKQNjY2OMjY1VXYYk9RzvoZEkKQNzc3NVlyBJPckrNJIkSZKyZaCRJEmSlC2HnHVJRLwSOAHoB/YH/gx4OzCbUhqssDRJm9Du3bup1+ur3m9gYICJiYkuVCRJ0sYw0HRBRPwh8A7gImACuB44r9KipDWYmppifn6evr4+RkZGqi5HLdTrdWZnZ6suQ+us8XtPktSagaY7fm9xmVL60eLKiNgOXFNNSVLnpqammJ2dZdeuXQaaHrdt2zYGBgbatqvX6ywsLGxARVqLxu+9ww47rOpyJKknGWi64zCAxjBTfvz1asqRtFUMDAwwMzPTtt3g4KBXdDIzNDRUdQmS1JOcFGAdRcSeiEjAceXHafHV8PFMq30iYjAijo+I8yPimoi4IiI+GhH32fh3IknqNaOjo4yOjlZdhiT1HK/QrK+ZcjkC3A94wyr2fSnwVODTwCzwcODZwNERMZBSun79ypQkSZI2BwPNOkopzQAzETEI3C+ltGcVuz8BeFhK6auLKyJiGngO8DTg48vtHBF7l9h01CpqkJZVr9cZHBysugy10MkMZ4v7+TXtXY1f18UHa/b391dVjiT1JANN73hnY5gpvZci0Pw2bQKNtBEWFha872KT8Wuaj7GxMQCmp6crrkSSeouBpndc0GLdD8rlQe12TintbLW+vHKzYw11Sb+y0hm0tPE6nbXMr2lvczY6SWrPQNM7rmyx7qZyeceNLERaykpn0NLG63TWMr+mvc3Z6CSpPWc5kyRJkpQtA40kSZKkbDnkTJI2kZXOWtbprGiSJPUaA42ktkZGRhgcHKSvr6/qUtSGs5ZtLo3fe2eeeWbV5UhSTzLQSGprZGSk6hLURqczlTnDWW9r/N7btWtXdYVIUg8z0HRBSmlwifXRYt0eYM8S7eeB2+0jSc0mJiaqLkFd5gM1Jak1JwWQJEmSlC0DjSRJGZicnGRycrLqMiSp5xhoJEnKQK1Wo1arVV2GJPUcA40kSZKkbBloJEmSJGXLQCNJkiQpWwYaSZIkSdky0EiSJEnKlg/WlCQpAz5YU5JaM9BIkpSB8fHxqkuQpJ7kkDNJkiRJ2TLQSJIkScqWgUaSpAwMDw8zPDxcdRmS1HMMNJIkSZKyZaCRJEmSlC0DjSRJkqRsGWgkSZIkZctAI0mSJClbBhpJkiRJ2dqn6gIkSVJ7o6OjVZcgST3JQCNJUgaGhoaqLkGSepJDziRJkiRly0AjSVIGarUatVqt6jIkqec45EySpAxMTk4CDj2TpGZeoZEkSZKULQONJEmSpGwZaCRJkiRly0AjSZIkKVsGGkmSJEnZMtBIkiRJypbTNkuSlIHp6emqS5CknuQVGkmSJEnZMtBIkiRJypaBRpKkDIyNjTE2NlZ1GZLUc7yHRpKkDMzNzVVdgiT1JK/QSJIkScqWgUaSJElStgw0kiRJkrJloJEkSZKULQONJEmSpGw5y5kkSRkYGhqqugRJ6kkGGkmSMjA6Olp1CZLUkxxyJkmSJClbBhpJkjIwNzfnwzUlqQWHnEmSlIGxsTEApqenK65EknqLV2gkSZIkZctAI0mSJClbBhpJkiRJ2TLQSJIkScqWgUaSJElStpzlTJIytXv3bur1+qr3GxgYYGJiogsVSZK08Qw06pqpqSnm5+fp6+tjZGSk6nKkTaderzM7O1t1Gdog4+PjVZegNvy9J1XDQKOumZqaYnZ2ll27dvmDXeqibdu2MTAw0LZdvV5nYWFhAypSN/T391ddgtrw955UDQONJGVuYGCAmZmZtu0GBwe9oiNJ2nScFKCFiOiLiBQRUxFx/4g4LSJ+FhFXR8SZEfGgst09IuI9EXFpRFwXEV+OiOOa+josIl4fEedExGURcUNE/CgipiPiN9sc+8iI+FhE/DgibomIwQ36FEiSeszk5CSTk5NVlyFJPcdAs7w+4EvAocAUcCbwGGAmIh4AnAc8DPgY8HHgaOCMiDiioY9HA68FrgROB95e7nc8cH5EHL3Ese9fHrsPOBV4D3DVur0zSVJWarUatVqt6jIkqec45Gx5u4ATU0q/uhMzIv4KeCNF2Pg48NKU0i3lts8DHwT+rHwB1IBDU0pXN3ZcBplzgLcCT2xx7EcCb0kpvW4lhUbE3iU2HbWS/bupXq8zODhYdRnSptPJDGeL+/k9mZ+LL74YwK9dD+v0e1LS2hholjdPETganUIRaPYDXrMYZkrTwPuBX92dm1L6cauOU0pfiYga8LiI2DeldGNTk8uBN6yt/N6wsLDguH2ph/g9mbcf/7jlrxVJ2rIMNMurp5Rublr3o3L5zearLimlmyPicuDwxvUR8WTgBOChwN25/ef97sClTeu+klK6fqWFppR2tlpfXrnZsdJ+umGlMzBJWp1OZy3zezJPi1dotm/fXnElWoozCUrVMNAs73Y/lVJKN0VEy22lm4B9Fz+IiFcBE8DPgc8D3weuARLwdIr7bvZr0c9laym8l6x0BiZJq9PprGV+T+ZpeHgYgOnp6Yor0VKcSVCqhoGmiyJiH2APRTjZkVK6tGn7McvsnrpYmiRJkrQpGGi66+7AgcAnWoSZu1DxUDBJUj58sKYktWag6a4fUwwv2xkRd0kp/QIgIvYF3kEReCRpTVY6a5kzMOVtfHy8fSNJ2oIMNF2UUrolIt5J8Ryar0bEp4A7AccBBwNnlf/elEZGRhgcHKSvr6/qUqRNzVnLpN7g7z2pGgaa7vsr4CfAKPASiskEPg+cyCaZlnkpIyMjVZcgbWqdzlTmDGdSd/h7T6pGpOS955tZROzdsWPHjr17l3rupiQpB85yJqnX7dy5kwsvvPDCpR4n0i132MiDSZIkSdJ6MtBIkiRJypaBRpIkSVK2DDSSJEmSsmWgkSRJkpQtA40kSZKkbPkcGkmSMjA6Olp1CZLUkww0kiRlYGhoqOoSJKknOeRMkiRJUrYMNJIkZaBWq1Gr1aouQ5J6TqSUqq5BXRQRPzvggAMO3r59e9WlSJLWYG5uDoD+/v6KK5Gk1i6++GKuvfbaK1JKh2zkcQ00m1xEzAF3A+YrLiUHR5XLr1dahTYLzyetN88prTfPKa23o4GbU0r7beRBnRRgk0sp+ae8FYqIvQAppZ1V16L8eT5pvXlOab15Tmm9LZ5TG817aCRJkiRly0AjSZIkKVsGGkmSJEnZMtBIkiRJypaBRpIkSVK2nLZZkiRJUra8QiNJkiQpWwYaSZIkSdky0EiSJEnKloFGkiRJUrYMNJIkSZKyZaCRJEmSlC0DjSRJkqRsGWi0pUTEvhHxqoj4QETUI+KGiEgRMbqGPo+NiH+LiCsi4tqI+K+I2B0Rd1zP2tXb1us8KM/HpV7ndat+bbyIODwi3h8RP4qI6yNiPiImIuKgVfZzcLnffNnPj8p+D+9W7epN63FORcRMm59D+3fzPah3RMTxEfEPEXF2RFxVfv0/3GFf6/Lzbin7rEcnUkZ+DZgo/305cBlw3047i4inAacD1wEfA64AngK8HXgE8PtrKVZ56MJ58D1gqsX6SzqvUr0kIu4PnAvcE/gU8HXgt4FXAU+IiEeklH62gn4OKfs5EqgBHwWOAl4IPDkijkkpfbc770K9ZL3OqQZvWGL9TWsqVDk5ETga+AXF75+jOumkC+fm7aWUfPnaMi/gTsATgXuXH+8BEjDaQV93A34MXA88tGH9/uU3bgL+sOr37Ku7r/U+D8r2M1W/L1/dfQGfK7/Wr2ha/7Zy/Ukr7Ofksv3fN61/Zbn+s1W/V18b81rHc2qm+O9h9e/JV7Uv4DjgAUAAg+V59OEO+lmXc3O5l0POtKWklG5IKZ2RUrp0Hbo7HrgH8NGU0gUNx7iO4q8aAH+6DsdRb/M80KqUf618HDAP/GPT5r8Gfgk8PyJ+rU0/dwGeX7bf07T5XRRX+h4fEb++9qrVy9brnJIapZTOSil9K5XpoxMbdW4aaKTODZXLz7bY9gXgGuDYiNhv40pSBbpxHhwYES+KiNdFxMsi4nfWXKV6yXHl8syU0i2NG1JKVwPnAHcG2n3dfwc4ADin3K+xn1so/iraeDxtXut1Tv1KRDw7Il4bEa+OiCf6u0wdWvdzsxUDjdS5B5bLbzZvSCndBMxR3KfmX0c3t26cB0cD7wPGKf7S/p/lJBa/tcZa1RuWPGdK3yqXR25QP8pfN86FjwJvAf4e+Dfg+xFxfGflaQvbkJ9TBhqpc9vK5cIS2xfXH7gBtag6630evI1iIoF7AHcFHgacRhFyahFxnw7rVO9Yr3PGn0FatJ7nwqcoJjU5nOIK4FEUweZA4GMR8YQ11KmtZ0N+ThlolJ1yqr/lppRsfnU0xaC2jl46p1JKf55SOjel9NOU0i9SSheklH6fYha1uwN/0a1jS1JK6e0ppc+klH6YUroupfSNlNLrgD+n+H/jWyouUbodp21Wjr5DMT3uSv2oS3Us/lVh2xLbF9df2aXja/2s5ZzaqPPgJOBZwKPX2I+qt17njD+DtGgjzoVJiqnoByLirs33bUlL2JCfUwYaZSel9LtV11D6BvBQinGfexs3RMQ+QD/FfP0+A6LHrfGc2qjz4Cfl0lmK8veNcrnUmPEHlMulxpyvdz/KX9fPhZTSdRFxNXAQxc8hA41WYkN+TjnkTOpcrVy2Gk/8aIpZO85NKV2/cSWpAht1HizOAGNAzt9Z5fJxEXGb38MRcVeKe6iuAc5r0895wLXAI8r9Gvu5A8VUqY3H0+a1XufUkiLigRRh5mrgp532oy2n6+cmGGiktiJiW0QcFRH3btp0GsUP9T+MiIc2tN8feHP54T9tUJmqzqrPg4i4c3lOHdG0/sERsW/zASLiwRQzngF4T1jmUkrfAc4E+oCXNW1+A8Vfvz+UUvrl4sryfLnNU7pTSr8APlS239PUz8vL/j+XUjIEb3LrdU5FRH9EHNzcf0TcA/hA+eFHyxkcpV+JiH3Lc+r+jes7OTc7Ov4anpUjZSkiXksxawvAAMXsUedy69SBX0wpTTa0H6H4QX5KSmmkqa+nU/yH9jqKKS6vAJ5KMU3hacAfrOWBVMrDas+DiBik+KvVbEppsGH9FMXsQmcDPwCupzhXnwDcEXgv8BLPqfyVv/TPBe5JMavUxcDDKZ7Z8E3g2JTSzxraJ4CUUjT1c0jZz5EUVwvPB7YDTwN+XPbznW6/H1VvPc6p8vfdScAXKa4GXwEcATyJ4l6HC4DHppS8L2sLKH+3Pb388F7A4ynOi7PLdT9NKf1F2baP4jEF30sp9TX1s6pzs6Na/b2orSYiZoBdyzS5TXBZLtCU2x8BjAHHAPsD3wbeD7wzpXTzetWt3raa82CZQPN04I+AB1P84N8f+BnFfyLem1L6dNffiDZMRNwXeCNFYD0EuBT4JPCGlNLPm9q2DDTltoMpnrj9dODeFOfMGcDrU0qXdPM9qLes9Zwqn3X158BO4DDgbhRDzL4GfBw4OaV0Q/ffiXpBROyh+NmylF+Fl+UCTbl9xedmR7UaaCRJkiTlyntoJEmSJGXLQCNJkiQpWwYaSZIkSdky0EiSJEnKloFGkiRJUrYMNJIkSZKyZaCRJEmSlC0DjSRJkqRsGWgkSZIkZctAI0mSJClbBhpJkiRJ2TLQSNImEhF9EZEiYqrCGo6MiBsi4n9XVUNZR0TEVyLi7DX00RPvJRcRceeIuCwiPlx1LZK2DgONJGm9vQ34GfCuxpURMVWGrZGldoyIPWWbPUusb3xdFxHfjoj3RERfc18ppQS8HnhkRBy/nu+lrGlHRPxzRFxehp7vR8S7I+LQJd7bfIv3sPi6rEX7O0XEmyNiLiIWIuKsiNixRN+PKfv5vQ7fJxFxSET8VUScGxE/jYgbI+JnEXF2RLyu+X1FxEx5zMHFdSmla4C3AMMR8bBOa5Gk1din6gIkSZtHRBwLPBkYK/9zu95mgZny34cAQ8CLgeMj4uEppW81Nk4pfSoiLgbGI+L0MuSsyHLvpQwOn6D4PfqvwDeBo4ATgKdExCNSSt9v0e0CMNFi/S9arHsr8GfA6cAlwPOBsyLiqJTSpQ213AV4L3BqSukzK31/Ld7Ph4FtwLeBTwI/Lj9+OPBm4HUR8RsppduFryYnA38NjAOP66QeSVoNA40kaT29DLgF+GCX+p9JKe1Z/CAi7kARKJ4EvA54YYt9TqEIB78L/PsqjtXyvUTE/sAksC/wrJTSJxq2PQeYprii89QWfV7ZWP9SIiKAlwAfSCm9qFz3SYow93zg/zY0fytwAPCqFb6v5mPtoggwN1F8/k5pDn4R8VvAO4D92/WXUrouIj4GvCQiHtAcMiVpvTnkTJK2iIi4d0T8Yzn06YaI+ElEfCIidi7RfltETETEJeXwrq9HxKsj4tdb3acTEXcDjgfOTSldsgFviZTSLcBiHUsNcfpoufzjlfbb5r0cCxwKXNAYZsp6PgJ8Bfi9iLjfSo/Xwj2AOwPnN6xb/Pev+o2IRwEvBV6eUvrZag9SBsKTKf7A+aqU0lSrq1gppa8CjwF+uMKuPwoE8KLV1iRJq+UVGknaAiKiH/gicBhQAz4C3Bf4feDJEfGsxuFK5VWIGrADuAg4lWL40RjwqCUO82jgTuVxqnBjq5Uppe9FxA+Bx0RErHDY2XLv5V7l8rtL7Ptd4GiK4XAfaNq2X0Q8DzgC+CXwX8AXUko3N7X7KXAt0Bg2H1ouvwcQEQcA7wM+kVI6rd0bWsIu4IEUQeV9yzUsw+MtK+z3fIqvx2OBv+ywNklaEQONJG0NJ1GEmRNTSuOLKyPi3cAXgFMi4n4ppcV7OV5DEWY+CgwvhoCIGAcuXOIYjyyXF7Sp5emtbuIvDbbZ9zYi4o7ceuVluSD1ZeDpwHbgf1bQ9XLv5aflsn+JfX+9XD6wxbZ7AR9qWjcXES9MKc0urkgp3RIR7wFeGRHbKALH84GrKMIlwJso7iN62XJvpI3F9znTIlR1LKV0bUR8DXhIRNw1pXT1evUtSc0MNJK0yUXE4RQ3Z3+f2957QUrp3Ij4CPA84Jncer/ICyj+Gv+XjVc0Uko/iIgJipvEmx1RLi9tsa3R08pXJwYbZkA7mOIKwFEUIeVNy+y3eCP7Eaws0Cz3Xs4BrgQeFhFPSyl9anFDRPwBxdUZgIOa9vsAcDbwNeBqiuDzcuBPgDMi4piU/n979xdqWVUHcPz7y0pUZhwwKpsMwaJ/EEaQUpJGMEylERX90WwMkhmjNBujDObq9KAwFaRjNFPRg0kPgQoaMTjkjFmBWI2SUD5UYkyTaWWTNYyUvx5+69jh3L3v2fc2d5wz8/3AZXH2XnvvtfZ9uOd311q/lQ+O1f88NUrzEaqve4CrMnNvRLwZ+AxwCfB4eyfrqalqDwKXZ+ZPB/Tz1FYuxxTBPwFnAquB3yzD/SUJcA2NJB0L3tjKezOza1rW3eP12vqRM4C9mflIR/2+kZBTWvm3Ke35eGZG1w+wecq151IZtK4BPk0FMw8Ab8nMPy9w3V9b+aIp9x/p7Utm/pNagJ/AbRFxe0RsiYg7qBGtB1rVZyau25yZd2fmY5n5r8x8KDM3UKmhTwCunah/MDOvzszTM3NlZp6bmfdHxAup4GhHZt4CXN7ex3bgnVSwtKMvffRhtNh3LklLYkAjSUe/k1vZN3IyOr6qlStb+VhP/b7jB1o5NRPW/2FzC3yOo0ZRbqRGAb7fFrj3OaGVBxaoM27BvmTmzVTWtJ3UNLkrWnsu4X9TwhYKsMZta+XbBtafo0Y91rfPnwN+lJnXZuZO4GPASQybijb63a8e+OzFWOw7l6QlccqZJB39/t7Kl/acP3Wi3v5W9v2Hv+/46Av8KT3nD5m2QP0PwBUR8TIqI9mnqACny6hNQ4OMqX3JzF3ArsnjETGatnf/wGc93sqTplWMiDOpqWiXtalnK6m1UaMgisx8NCKeAF4/4Nmj0bbzIuK4Q7mOhsW/c0laEkdoJOnot6eV50RE1z+y3t7KXwJk5n4qU9fqnsX753Qcg8rYBTUN7HDaCBwE5toX/C6voaaA/WrgPZfUl4hYBVxABSk7B152div7sqaN7v18aqrZ7sz89sTp4yc+Dx0luwd4GHg53Xv4jD//eRHxgoH3hUqK8BeWZ32OJD3LgEaSjnJtH5WdwOnUQvJnRcRZwIXUWpHbx07dTP2NuL5t8jiqf9rkPcbsbuXZPeeXRWY+CnyLGhHYOHk+Io6npqXtycwnB952dys7+xIRKzqOnUht4rkKmMvMg2PnXhsR80ZgWsB4U/t4y5Q2fQF4JXDp6EALPvcCa0fBatsocwWVfGBBbaRrPbWp5o0R8dHx3/dYO18H3MXAqWktTfhLqOBrSJpsSVoyp5xJ0rFhA5Wd68sRsYZKRzzah+YZaqH+eGrdLVSa4w8Dr46Iu6i1OB+k0jy/l/mL3h+KiIeBdyzD9KVprqPSN18ZEVsz84mxc+dRe8rcOvRmA/qyLiI2UoHPPiqYuoCavndDZm6bqP8hYGNE/JjaR+YfVOKFd1OjKT8EvtLXnhZQbAI2diRq2ALcANwbEfcBFwFPAV8f2Nd7IuJ9VDrp7wKbImI3Ncp0MrX/zVnUvjlD18OsaeXgdy5JS+UIjSQdAzLzd9QX023UVKCrqIxYO4C3jqcebvUPUFPRtlJrb65sn68Drm/V9jPfN1r9NR3nlk1m7mvPXsH8jRzXAU8zZePIDgv15efAr4G11Lt8P5Xd7F2Z2TWCtQv4ARXEXAh8lsrY9pPWvvMz8+muRrS9dr4D3Ed3kLKVCnZOAy4Dfg+szcy+5A3zZOadrW1z1DSxD1BrdS6isrnNAWcs4p7rqIDIgEbSsgtHgiVJixERlwLfBDZk5vaJcyuB3wI/y8yl7jVzyETEi4FHgO9l5icWee0R1ZdZERFvoPbC2ZSZXfsVSdIh5QiNJKlTyx42eewV1GjAv4E7J8+3NR3XAO+JiDcteyOn+yLwH6rNi3IE9mVWfInKQPfV57ohko4NrqGRJPW5tWW1+gXwJJVU4HzgRODqzPxjz3XbqYXxfWmiD4u2uH0fcHGbkrYUR0RfZkVLjLAH+FqbtihJy84pZ5KkThHxSeBi4FXU4vCnqC+rN2Xmbc9l2yRJGjGgkSRJkjSzXEMjSZIkaWYZ0EiSJEmaWQY0kiRJkmaWAY0kSZKkmWVAI0mSJGlmGdBIkiRJmlkGNJIkSZJmlgGNJEmSpJllQCNJkiRpZhnQSJIkSZpZBjSSJEmSZpYBjSRJkqSZZUAjSZIkaWb9FwMMwSbEoIAyAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 263, "width": 410 }, "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ctv.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the above scaled Schoenfeld residual plots for `age`, we can see there is a slight negative effect for higher time values. This is confirmed in the output of the `CoxTimeVaryingFitter`: we see that the coefficient for `time*age` is -0.005.\n", "\n", "#### Conclusion\n", "\n", "The point estimates and the standard errors are very close to each other using either option, we can feel confident that either approach is okay to proceed. " ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.2" } }, "nbformat": 4, "nbformat_minor": 2 }