bootstrap_logrank_stratification.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Required packages and example dataset\n",
    "from lifelines.datasets import load_rossi\n",
    "from lifelines.statistics import logrank_test\n",
    "from sklearn.utils import resample\n",
    "import statistics\n",
    "from lifelines import KaplanMeierFitter\n",
    "kmf = KaplanMeierFitter()\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import argparse\n",
    "from pylab import rcParams      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stratify_plot(df,variable,time,event,iterations):\n",
    "        '''\n",
    "        Given a dataframe and a desired variable, and time and event,stratify based on a log-rank test\n",
    "        '''\n",
    "        rcParams['figure.figsize'] = 10, 8\n",
    "        minp=np.int(len(df)/5)\n",
    "        for f in range(0,iterations):\n",
    "            boot = resample(df, replace=True, n_samples=len(df))\n",
    "            T = boot[time]\n",
    "            E = boot[event]\n",
    "            start = boot[variable].min()\n",
    "            stop = boot[variable].max()\n",
    "            step = np.abs((start-stop))/100\n",
    "            cutoff_list = np.arange(start, stop, step)\n",
    "            p_list = []\n",
    "            num_list = []\n",
    "            teststast_list = []\n",
    "            resdic={}\n",
    "            minvalue=[]\n",
    "            i = 0\n",
    "            for Cutoff in cutoff_list:\n",
    "                low = (boot[variable] < Cutoff )\n",
    "                high=(boot[variable] >= Cutoff )\n",
    "                if len(low) >=minp and len(high)>=minp:\n",
    "                    results = logrank_test(T[low], T[~low], E[low], E[~low], alpha=.99)\n",
    "                    p_list.append(results.p_value)\n",
    "                    num_list.append(T[low].size/low.size)\n",
    "                    teststast_list.append(results.test_statistic)\n",
    "                    if results.p_value==results.p_value:\n",
    "                        resdic[Cutoff]=results.p_value\n",
    "            minvalue.append(min(resdic, key=resdic.get))\n",
    "        T=df[time]\n",
    "        E=df[event]\n",
    "        bestvalue=statistics.mean(minvalue)\n",
    "        low = df[variable] < bestvalue\n",
    "        results = logrank_test(T[low], T[~low], E[low], E[~low], alpha=.99)\n",
    "        ax1 = plt.subplot(111,alpha=0.4)\n",
    "        kmf.fit(T[low], event_observed=E[low], label=\"low {} \".format(variable))\n",
    "        kmf.plot(ax=ax1,show_censors=True, censor_styles={'ms': 6, 'marker': 's'})\n",
    "        kmf.fit(T[~low], event_observed=E[~low], label=\"high {}\".format(variable))\n",
    "        kmf.plot(ax=ax1,show_censors=True, censor_styles={'ms': 6, 'marker': 's'})\n",
    "        lvalue=-np.log2(results.p_value)\n",
    "        ax1.title.set_text(\"p<{}  | -log2(p)={}  | Cuttoff={}\".format(round(results.p_value,5),lvalue,round(bestvalue,3)))\n",
    "        plt.savefig('myimage.svg', format='svg', dpi=1200)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAHwCAYAAABkCCFxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucnGV99/HvL8kmG5KQ0+ZQcmCzAUIIhwRjUCJhIahBqZRK1aiVFBWKgoqHPtjKKYiipU+1VB4LSrFaEQoqgQLSViggWBMEgQSskAOEQ0g2yZLDLslufs8f970wmczuXLu55rD3fN6v174yM/c111xzyOx37+tk7i4AAADsvwGVbgAAAEBWEKwAAAAiIVgBAABEQrACAACIhGAFAAAQCcEKAAAgEoIVUICZNZrZ2l6Ubzaz9SVsUu5jjTOz35tZfUDZCWb2tJkNKUfbUNvM7EYzW1LpdvTEzL5qZpvM7JX0+hlm9oKZbTezOZVuH/o/ghX6xMwWmtkzZrbTzO4zs4N7KNuYltmZ3ueUvOMXmtkrZtZqZjfkhoCe7mtmZ5nZo2b2mpmtN7NvmtmgnOPb8346zeya2K9FTGY23sxuMrOX0tfjV2Z2XF6xiyT9s7u3F6vP3TdIuk/SOX1oy5Fm9ov0l9A+C96Z2flmtsLMXjezG4vUZekvtBfT53W/mc0qUG6MmW00s4dybvtI3vu408zczN6SHr/MzHbnlWkqUPdZ6f0+kXNbt/c1s8PM7Pa0PZvT12JGXn09ff5mmtkv0+f7rJmdkdeeA8zs2vT1bTWzB/KOH2tmD6Rt2mBmny3wnE5Mn9NXc25bkn7Wc59Tc87xK8zsSTPrMLPL8up7r5k9ZGZb0/+T15vZiPzH7SszG5y+5n8wsx1mtjb9P98YcN/G9LnmvsZLcj8rAXVMkfQFSUe4+8T05qslne/uw939scB6hpjZ981snZltM7PHzOzUvOd5a/r8PPf176a+MWb2s/Q1WWdmH8451mxme/Lez7NCnzPKj2CFIGZWZ2Zj0ssNkn4q6WJJYyStkHRzD3e/SdJjksZK+htJt5rZuLSudysJCgslNUpqknR5yH0lHSDpc5IaJB2X1vHFrjumX5TD3X24pAmS2iT9W99egbIZLmm5pLcoeW1/IOnfzWy4lHyhSzpL0o96Uee/Sjq3D23ZLekWSR/v5vhLkr4q6YaAuv5M0tmSTlDyvB6R9MMC5b4h6encG9z9X/Pey09JWi3ptznFbs4t4+6rc+sws9GSvixpZYHH7O6+oyQtkzRDyefnN5Juz7lft5+/9Jf/7ZLuTJ/vOZJ+ZGaH5dz/uvTYzPTfC3Pa2yDpHkn/pOSzf4ike/OeU52kb0v6nwLP6ZG853R/zrFnJf2VpH8vcL+RSt7Tg9J2TZb0twXK9dWtkt4n6cPpYx0j6VElr105HCypxd1fzbut0OeiJ4MkvSDpRCXP42JJt+QFxIckfVTSKwH1fUfSLiWfs49I+n95f3i8lPd+/qCX7UU5uTs/NfQjaa2SXzCrJG2R9M+S6nsof6Skv5O0QdIH09vOkfRwTplhSkLL4QXuf5ik1yWNyLntQUl/mV7+saSv5RxbKOmVkPsWeKzPS7qjm2NnKfllbIGvU6Oktb14XZslrc+5PlPS/ZK2KvnSfl/OsbGS7pD0mpIQ9VVJD/VQ92uS3pJeXiDp2bzj90v6upJf/K1KfqGPyTk+SNJOSQf38TNzSPJV0e3xr0q6sUgd/0fSLTnXZ0lqzyvzdiWB6y+KvB73Sbo05/plkn5U5PG/qySQ3S/pE725b07ZMZJc0thin7/0/8323M+bkmB0RXp5Rvq+HthNXV+T9MMi7blI0jcl3Sjpqzm3L+np9csp9yNJlxUp86eSnuzFZ+VGSUu6OXaKku+JKT3cf62kUwq9P5KeT1//7enP2yW1S+pMr29Ny42U9C+SNkpaJ+krSk4idD3+nrT8Tem/LmmHpOf68v8jp61PSHp/gdvXS2ru4X7DlISqw3Ju+6Gkq9LLzcr5buGn+n84Y1WbPiLp3ZKmKwkvX8k9aGajzexTZrZcyS+DPZJOdveus1KzJP2uq7y775D0XHp7vlmSVrv7tpzbfpdTdq+60ssTzGxswH3zLVD3f3meJelfPP2mKqX0TMIdSl678ZIukPSvOd1I31HyRT4xbVe3p/XNbLakwUrOMkjSUZJ+X6Dox5ScETpIUoekf+g64O4d6f2PSev8cNrV093P1D4+9Z78RNIhafdanZLnfE/XQTMbqOR1OV/JL7qCLOlyXqDkF2euP06761aa2Xl595knaa6ScFVIt/fNs0BJ6G/p4XjX588KNV9J4JKSM1zrJF2edgU+aWbvzyn7NkmbzexhM3vVzO7IfV/S1+FsSUu7acuctN7/NbOLc7vPeqmn/1O9dYqk37j7C/vRFkka5clZm0ck/aXePDs3Kj1+jZJw1aTkjNLHJP2Fu/+npFP15tmfxZ6cAZWkY9x9uiSZ2Z09/N+4s1DDzGyCku/SvrxWh0nqdPf/zbkt/3tufNodvMbM/t7MhvXhcVAmBKva9I/u/oK7b5Z0paTFkmRmB5rZTyStUfJX0qVK/rr8krvnfmEMV3JmJFerpEJjMYqVzT/edXlEbx7HzP5CyS/Pqwscm6rkC7Zcp8/fpqTtV7n7Lnf/pZIuocVpgHi/kjMuO919VXftMrMDlfzlerm7d70OoyRtK1D8h+7+VBpyL5b0gfSxumxL7yt3/7G7j+rh5/n9fwn28bKSs42/V3LW4M+U0/Ul6TOS/sfdHy1Sz8ckPejua3Juu0XJGcJxkj4p6RIz6/pMD5R0raQL3H1Pgfq6vW8uM5usJPh9vlCjCnz+npH0qqQvpd3o71LyGTwgPT5ZSchqVRKGz5f0AzObmXP8LEmflTRVyf/Jm3Ie8h8kXezu2ws054G07vFKPmuLJX2pULt7YmbvTNtwSW/v242xSj4HJZO+3x+U9GV33+bua5Wccf/z0Drc/bQe/m+cVuAx65R0t//A3Z/pQ7OLfc89I2m2pD+SdLKSYQL/tw+PgzIhWNWm3L8Y1yn5YpekOiVfyJuV/MX0lLt3Frj/dkkH5t12oAr/wi9WNv941+VtoY9jZn8i6SpJp7r7pgJt+JiSrpE1BY71mplNzR1IWqDIQZJeyPtFvk7SJCW/wLvGZ3TZ5y94Mxuq5KzXr9396zmHtqhwgM1/T+uUjP3pMkJJt2SlXCrprZKmSKpXMo7ul5YM4D5ISbD6m4B6Pqa8IOruq9z9JXfvdPeHlYw7OjM9/ClJT6RnN/ZR5L6SklmYSs4+XuvuN+XXUejz5+67Jf2JpPcqGWPzBSUhrmvmaJuSMWxfTcP3fyvp4nxXzvGfuftyTyYpXC7peDMbaWZ/rKR7vOC4Rndf7e5r3H2Puz+p5KzWmYXKdsfM3qakm/7MvDMp+6NFSTgopQYlZ3jX5dzW9X8vOjMboOSPn11KwnFf9Pg95+6vpJ/TPel32F+pl+8nyotgVZum5FyeqmQQsty9xd2PVPIX3yRJv7VkVtOSrsHTqZVKu5UkKT0tPV2FT4OvlNRke88sOian7F51pZc3pN0txe4rM1sk6XpJf5z+Eilkn1/G+8Pdn/e9B1Pne0nSlPRLt8tUSS8qGffRoeSMRJfc96NrgPrP0/L5g86fUNJ1kC//Pd0taVNa3yAl46R+l17Pn2WX/1OKrsBjlAwSX+/uHe5+o6TRko6QNE/JL9xVlkyB/7akeZbMSnvjrJuZzVcSWm8t8liuN7viFko6I63rFUnHS/o7M/vHgPt2DXq/V9Iyd78yv3BPnz93f8LdT3T3se7+biVdU79JDz9R5Dk8ob27RLsuW/qc5uY8pw9K+pyZ3a7C9npOxViy5MAySWe7+3+F3i/Afyp5Xyf3UGaH3jyrJyXd5V0KdRHn37ZJyWc/d5Zy1/+9IGZ2dw//N+7OKWeSvq9kwPn70zDdF/8raZCZHZpz217fc3l69X6iAkoxcIuf6v1RMjj0SSW/2Mco6Z75Wjdlhyj50r5byUDbRent45Scqn6/krMP31ByZqW7x/y1ki6SeklnKDlzMi49tkjJX/RHKPlF+0ulgzYD7nuykr+CF/Tw2Mcr+bIeEfoapfdrVB8Hryv5i/k5JYOL69Jj25QO7lcyg/LHSn6BHK5kUO5D6bGu8Vk/lzSowOMMVhLOJuXcdr+SMyFHpHX+m6Qf570Gq/rwWbH0dT9CyZd5vaQhOccHpbd9Xclf7fWF2pyWvVTJLKkJSv6g+/P0fRmVfs4m5vx8VslMt4l5dVynZJxcft2np58dUxLSXpR0VnpsVF7dDyvpzhsZcN8DlQShf+zmOfX4+ZN0dPqaHKBktuCartcvfZ+fVdJtO0jS/LzPyMlKzk7OTsv+vZIuUCk5+5j7nG5Oj49Jj58qaUJ6+XBJT2nvwf51abt+rGTiQb2kgemxI5UzUaUPn5kb1c3g9fT4Mr0563VQ+lz+UkmIk5IutR+nbZyrJCh1DV4/QMlA9dxB3ouUfKcNzrntR5J+ltZ9sJKutE/k/z/NKe+SDunDc/2uku+n4d0cH5K+tuuVnImsVzeTZ5SMQbxJyUD2+Uq+X2fltHlq+hmdouTM5j/35f3hpzw/FW8AP2V+w/eeFbhVyZmcAwLud5CkaTnXT0m/sNqU/GJvzDn2XUnfzbnemJZpUzLG5pS8uj+ffpm/pmSW4pCQ+6ZfMB16c5bQdkl359X9Tyoyu6qb59uo/ZsVOEvSf6dfkKsknZFzbJySqe5dswK/Iem/0mMnpl/0O/Oe1wk59/9bSf8n5/r9enNW4GtKgllDzvHvSPpMH18Dz/tZm3P8sgLHL0uPTU3bPTW9Xp+24+W0jb9VGtQLPO4S5c1qS++/VdLCAuVvUhJwtqefyW6fq/adFdjtfZWML+qaMZb7XnQ9px4/f+n7tKXrduX98k4/I4+k9e/1GUmPn6ck6G1J39OCs+m076zAq5X8f9qhZCbsUkl1eeXz37cl6bF/1puz5rp+VvbiM3Ojeg5Wg5V0az6btm+dpO/lvKZNSkL1diX/R/5BObM20+eyMf0svC2t79+VDF/YlJYZrSRcbVTSRX6JpAGF/p+mt/U6WCkJbK5kVmLua/WRnDJrC7zOjemxv877rIxR8sfUDiV/aH0459jn08/BzvT5XKNe/qHIT3l/LH3jUCMsWU38E57MkEE3LFmP5n53byzDY31DydmZswLLj1NypnGOu7eZ2f1Kfvl8r0DZ8UoC3hwPWFAU2B+WLBR7vyddvUBN6usUXAB9ZGaHK/lL+0klA7o/LukTPd4ph7tvVNLFE1L2VSWz3gAAZUCwAgrbKulbJap7hJIuqIOUTMn/O+29ojfQX/1cSRcYULPoCgQAAIiE5RYAAAAiIVgBAABEUrExVg0NDd7Y2FiphwcAAAj26KOPbnL3ccXKVSxYNTY2asWKFZV6eAAAgGBmtq54KboCAQAAoiFYAQAAREKwAgAAiIQFQgEAyJjdu3dr/fr1am9nJ6veqq+v1+TJk1VXV9en+xOsAADImPXr12vEiBFqbGyUmVW6Of2Gu6ulpUXr16/XtGnT+lQHXYEAAGRMe3u7xo4dS6jqJTPT2LFj9+tMH8EKAIAMIlT1zf6+bgQrAAAQ3fDhwyvdhIogWAEAAERCsAIAACXj7vrSl76kI488UkcddZRuvvlmSdKnPvUpLVu2TJJ0xhln6Oyzz5Ykff/739dXvvKVfeo577zzNHfuXM2aNUuXXnrpG7ffddddOvzww/WOd7xDn/nMZ3TaaadJknbs2KGzzz5bb33rWzVnzhzdfvvtpX6qkpgVCABApl1+x0qteum1qHUecdCBuvSPZwWV/elPf6rHH39cv/vd77Rp0ya99a1v1YIFC7RgwQI9+OCDet/73qcXX3xRL7/8siTpoYce0oc+9KF96rnyyis1ZswYdXZ2auHChXriiSd02GGH6dxzz9UDDzygadOmafHixXuVP/nkk3XDDTdo69atmjdvnk455RQNGzYszovQDc5YAQCAknnooYe0ePFiDRw4UBMmTNCJJ56o5cuX64QTTtCDDz6oVatW6YgjjtCECRP08ssv65FHHtHxxx+/Tz233HKLjj32WM2ZM0crV67UqlWr9Mwzz6ipqemNpRFyg9W9996rq666SrNnz1Zzc7Pa29v1/PPPl/z5csYKAIAMCz2zVCruXvD2SZMmacuWLbrnnnu0YMECbd68WbfccouGDx+uESNG7FV2zZo1uvrqq7V8+XKNHj1aS5YsUXt7e7d1dz3ubbfdphkzZkR9PsVwxgoAAJTMggULdPPNN6uzs1MbN27UAw88oHnz5kmS3v72t+tb3/qWFixYoBNOOEFXX321TjjhhH3qeO211zRs2DCNHDlSGzZs0N133y1JOvzww7V69WqtXbtWkt4YvyVJ7373u3XNNde8Eb4ee+yxEj/TRNFgZWY3mNmrZvZUN8fNzP7BzJ41syfM7Nj4zQQAAP3RGWecoaOPPlrHHHOMTj75ZH3zm9/UxIkTJUknnHCCOjo6dMghh+jYY4/V5s2bCwarY445RnPmzNGsWbN09tlna/78+ZKkoUOH6tprr9WiRYv0jne8QxMmTNDIkSMlSRdffLF2796to48+WkceeaQuvvjisjxf6+k0miSZ2QJJ2yX9i7sfWeD4eyRdIOk9ko6T9G13P67YA8+dO9dXrFjRp0YDAIDuPf3005o5c2alm1EW27dv1/Dhw+Xu+vSnP61DDz1UF1544X7VWej1M7NH3X1usfsWPWPl7g9I2txDkdOVhC53919LGmVmf1SsXgAAgP11/fXXa/bs2Zo1a5ZaW1t17rnnVrQ9MQavT5L0Qs719eltL0eou89+fe0nNWLr00FlG4YP0YQR9cULHnWmNPcv9rNlAAAglgsvvHC/z1DFFCNYFdpUp2D/opmdI+kcSZo6dWqEh95/O3d16pXWdg0ssjfQyNZntHPzBv3h5Y6ide6uG66NExf0WGbE0DqdNGN8r9oKAACqW4xgtV7SlJzrkyW9VKigu18n6TopGWMV4bG79bZPXR9U7oP/9Ihatr+uK95zVI/ljrh3seraWjR0zEFF6xy1/QWN3fQfPZZ5rW23tGt0UBvLrn6kdOg7K90KAAD6nRjBapmk883sJ0oGr7e6e0W7AUtlyI4X1bjiiqLlWicery2TF/ZYZltHu371SvHVLg4YMkhzpowKbmMU2zeU9/EAAMiIosHKzG6S1CypwczWS7pUUp0kuft3Jd2lZEbgs5J2SsrkIKRNje9TXVtL0XL129ZJUtFgNT5kTJekrTt3BZUDAACVVzRYufviIsdd0qejtagCXtzapqV3rixS6mjV7fqKfODgHkt9bc8VUqv01w/3XNuJk6RTD+5dOwEA6A/Wrl2r0047TU89te8SmJdccokWLFigU045pdv7X3bZZRo+fLi++MUvlrKZJVHzW9qcPnuSWra/HrXOaXvW6Wtt3XcZtnVKQ5+VGjcWr2tXxx61Plm83IABphFDIr2dBx2bzIAEACCypUuXVroJJVXzwerDx03V8CEDNS6ga27smju0e+i4HsuMXn+87BWpqYcyq1vD2zd4UNiuQx2de8Ir7cnmNVLn7jh1AQBqVmdnpz75yU/q4Ycf1qRJk3T77bdr6NChWrJkiU477TSdeeaZuuuuu/T5z39eDQ0NOvbYY7V69WrdeeedkqRVq1apublZzz//vD73uc/pM5/5zD6Pcd5552n58uVqa2vTmWeeqcsvv1ySuq13x44duuCCC/Tkk0+qo6NDl112mU4//fSoz7vmg1VsWyYvLDq+qqub8Kqi67eG27pzl+Yf0rD/Fd1zkdTxuvTkrcXLMnsQAKrf3RdJrwR0ffTGxKOkU6/qscgf/vAH3XTTTbr++uv1gQ98QLfddps++tGPvnG8vb1d5557rh544AFNmzZNixfvPfLomWee0X333adt27ZpxowZOu+881RXV7dXmSuvvFJjxoxRZ2enFi5cqCeeeEKHHXZYt/VeeeWVOvnkk3XDDTdo69atmjdvnk455RQNGzYswouSIFhlxK7OPfrVs5uKlguaZThoiDR8QvEHZfYgAKAb06ZN0+zZsyVJb3nLW97YKLnLM888o6amJk2bNk2StHjxYl133XVvHH/ve9+rIUOGaMiQIRo/frw2bNigyZMn71XHLbfcouuuu04dHR16+eWXtWrVKu3Zs6fbeu+9914tW7ZMV199taQk3D3//PNRt/8hWGUEswwBAAUVObNUKkOGDHnj8sCBA9XW1rbX8WJ7Feffv6Nj7wW616xZo6uvvlrLly/X6NGjtWTJErW3t/dYr7vrtttu04wZM3rzVHqFYFUhq1uli4rMHJSYPQgAyKbDDz9cq1ev1tq1a9XY2Kibb765V/d/7bXXNGzYMI0cOVIbNmzQ3Xffrebm5h7rffe7361rrrlG11xzjcxMjz32mObMmRP1eRGsKuDESWHluga5xwxWxboMj2zbncwwDKlsd3vYWKzYGNsFAP3e0KFDde2112rRokVqaGjQvHnzenX/Y445RnPmzNGsWbPU1NSk+fPnF6334osv1uc+9zkdffTRcnc1Nja+MVg+Fit2Kq5U5s6d6ytWrKjIY+db9viL0WYFxtR1Ruuq48v2kGpccYU6Ovdo5Bl/V74H7a3tG1gOAgB68PTTT0cdN1Qq27dv1/Dhw+Xu+vSnP61DDz00yobK+1tvodfPzB5196LTzsLm8gMAAER2/fXXa/bs2Zo1a5ZaW1t17rnnVnW9IegKBAAAFXHhhRdGOUNVrnpDEKyqXCUGuQ/b8XyynlUxTc3SYYviPCgAABlAsKpilRjk3jrxeA3r3KORxQpuXpP8W4lgFTJongHuAGqcu8vMKt2Mfmd/x54TrKrYqQeHhaWQM1qhtkxeqDVjTii+invIGa1SGR3worB4KYAaVl9fr5aWFo0dO5Zw1QvurpaWFtXXh60NWQjBCvsIWcX9yLZkP8GnYq32HlvoUhCc2QKQQZMnT9b69eu1cePGSjel36mvr99nhffeIFhhHyGruA8amEwoHXXA4KJlK7Lae8hZLYkzWwAyqa6u7o0tXVBeBKuMYCX3PuLMFgAgIoJVBlRyJfd+jzNbAICICFYZUIlB7pJUv22dGldcUbTcro49an2yeH0DBphGDAn4SDY1s8wDAKAqEax6oWPwgaprKz4QcEBnu14fPqUMLaqc1onh++wMHhS2wH9H557ihSq1zANdhgCAAASrXmiddGJQubFr7ihxSypvy+SF2jJ5YdQ6t+7cVb3LPNBlCAAIwF6BAAAAkXDGqgRCuwzLzTpH6rntg/TlBzt6LHfyhO1652FjytQqAACyg2BVAqFdhuU2r32Ddj+3Sbt7KLOuZae0QXrnYWVrFgAAmUGwqiELZ07QwpkTeiyz9M6VUltbmVqUQexjCAA1jWCFqhG6lc6wHc9rx8++0GO5jePerg0Tm4s+ZvTtdtjHEABqGsEKVSNkK50dk96hQa883OMHt37bOg3aPECvN72raH0V2W6HpRsAILMIVuhXQpZ5CFm0tKJYugEAMotghUwKXRX+5THHSYecWYYWAQBqAcEKmRO6Knz9tnUa17lHUpUGK7oMAaDfIVghc0JXhW9ccYUUso1OpdBlCAD9DiuvAwAARMIZK9S0YTueD9t/sKm5/Bs/h6LLEACqBsFK0oihddq4rb2sj9ne0akpo4eV9TFDPbd9kC56OF59J06STg3s1Sqn1onHa2jHHg1o62kt+iR8tb3eoRHVGqzoMgSAqkGwknTSjPFlf8xlj79Y9scMMX96g+raWuQaHKW+1a3Jv9UYrHozFmtPNY/FAgBUDYIV9rJw5gT96YHPaNCu4ptID+hs1+vDp/RYJuaZr0rKRJchAKDkCFbYR+gm0mPX3FHillSH1onHa1jnHo0sVnDzmuRfghUA1CyCFVDElskL9fuR8zV4YM+TaI988utS2249VWS/w+j7EwIAqgbBqkJCB8xX8yD3UKtbw7oEq3WQuxS2j+GgNHiNOqDn8WkV2Z8QAFAWBKsKCR0wX62D3EOdOCmsXDUPcgcAIBTBCiV16sFhYSkrg9wBALWNldcBAAAi4YwVEFH9tnXJHoQ9eHnMcdIhFdj4OWSFdlZnB4D9QrACImmdeHzRMvXb1mlc5x5JFQhWISu0szo7AOwXghX6rGPwgaprK76QaAjrTFaJqmtrLVo2ZGHSSghZyb1xxRUSq7gDQGYRrKpcNS/LELqQaIjdK1dKklqmLShatr8vTMoq7gCQXQSrKlcryzLUisxs/AwAKIhgBZQRGz8DQLYRrIAqVZEuw5CZg73BLEMANYZghaqxrmWnlt65smi5uraR8oFxHrNat9Gp2MbPITMHe4NZhgBqDMEKVWH+9AZJPW9eHFs1b6OzZfJCrRlzguYf0tBzwZAzWpUUegaMM1sAMoJghaqwcOYELZw5Iajs2DV3aPfQcfv9mGyjUwahZ8A4swUgIwhWGVHNyzKgxDavYfkGAKgSBKuMqKVlGUIWJq3WRUSja2oOKxd7LBYAoCCCFfqdkIVJ+/siopK0q3OPfvVskXFnA+bqgFlv05wpo3ouV+1jsQAgIwhWQJUaP6I+qNzWnbtK3BIAQKgBlW4AAABAVnDGCpkUMg4r2fh5j6SwM0MAABRDsEImhYzD2r1yperaWsrQmtIKGYt1ZNtuDRhgGlGmNgFArSJY1ZiQZRlYkqF/CRmLNWjgAHVU896DIQuJsogogH6AYFVjQpZlyMKSDKGe2z4oaKHQat36JjNCFhJlEVEA/QDBCjVr/vQGDXq9VdbZ86y657YPknXu0qkHMxYLANAzghVqVug2OkvvXCllYCxWv8e+gwD6AYIVUCP2uIovOBrogCGDii9KGhv7DgLoBwhWQI0YPGiARh0wOEpdLEoKAIURrIAAz22v05cfDAkTe+QDex6LVamB8PXb1qlxxRU9lmmdeLy2TF5YphYBQPYQrIAi5k9vkCTtDihbbF2s1a3Jv+UOVq0Tjy9apn7bOkkiWAHAfiBYAUWEDnKXks2fdw8d1+3xkKXgSEE9AAAdDElEQVQdSmHL5IVFA1Oxs1m5gjaI7oWKjNkCgBIgWAHotdANokMxZgtAVrAJMwAAQCScsQKQLax3BaCCCFYAsoX1rgBUEMEKwBtClmTojdDlG0IGwzPAHUB/QLDCPkYMrdPGbe3R6mvv6NSU0cOi1VfNOgYfqLq2jd0et86RkvZIqr59B0OWZOiN3izfEDIYPvoAd7oMAZQAwQr7OGnG+Kj1LXv8xaj1VbPWSSf2eHz3ypVF17qqlJAlGXoj5pmvkqDLEEAJEKwAlExo1yIrvgPICoIVUGbPbR8UtFBopba+iSW0a7HqV3wP7TIMQbcikHkEK6CM5k9v0KDXW2WdPY8Xem77IFnnLp16cPWNxQoV2rWYmS7DEHQrAplHsALKKHR7nKV3rpSqdCwWAKB7BCsAKJeY3Yqx0U0JREGwQsmFLt9QS8syoEbF7FaMjW5KIIqgYGVmiyR9W9JASd9z96vyjh8s6QZJ4yRtlvRRd18fua3op0KXb6ilZRnQeyGLiPYGC44CKIWiwcrMBkr6jqR3SlovabmZLXP3VTnFrpb0L+7+AzM7WdLXJf15KRoM1AofMKjHxUa7DOhs1+vDp5ShRZUVsohob0RfcLS/i91NSdcialTIGat5kp5199WSZGY/kXS6pNxgdYSkC9PL90n6ecxGArWoY8hItUwrvmTB2DV3lKE1yLzY3ZR0LaJGhQSrSZJeyLm+XtJxeWV+J+n9SroLz5A0wszGujvTmgCgFoWcAeOsFjIoJFhZgds87/oXJf2jmS2R9ICkFyV17FOR2TmSzpGkqVOn9qqhALIrZIX22Kuzh47ZYixWH4WcAeOsFjIoJFitl5Q7gGOypJdyC7j7S5L+VJLMbLik97t7a35F7n6dpOskae7cufnhDEAfFNv4uUu1jsUKWaG9FKuzh47ZYiwWgN4ICVbLJR1qZtOUnIn6kKQP5xYwswZJm919j6QvK5khCGA/rGvZmSwUWsT86YcHLTparWOxQlZor/rV2QEgVTRYuXuHmZ0v6RdKllu4wd1XmtlSSSvcfZmkZklfNzNX0hX46RK2Gci8+dMbJBXvplrXslPSpqBg1d+xoTOA/iBoHSt3v0vSXXm3XZJz+VZJVbqcMND/9GrrmxqQmQ2dAWQeK68DqHqV3NA55sKkDIQHso9gBQA9iLkwKQPhgewjWAH9XOgg93eNrNe7Di1DgwCghhGsgH6sN4Pc72sbQrBCdWEbHWQQwQpVY8TQOm3c1l7Wx2zv6NSU0cPK+pgx9WqQe1tbGVqEnrAoaR620UEGEaxQNU6aMb7sj7ns8RfL/pioXSxKCmTfgEo3AAAAICsIVgAAAJHQFQjUiOe2D9JFD8ep68RJ0qmRh8cAQBYQrIAaMH96g+raWuQavN91rU63V6/WYJWFrW8Y5N5HsWcZhmAmIvIQrIAasHDmBH2g/jfaPXTcftcV66xXKWRl6xsGufdR7FmGIZiJiDwEK6BGdAw+UHVtG3ssM6CzXa8Pn1KmFsVXya1vAEAiWAE1o3XSiUXLjF1zRxlaUh2y0GWIKhDa/UiXYc0gWAHotdWtcbsEyz0YPitdhqgCod2PdBnWDIIVgF45cVLc+ioxGJ4uQwClQrAC0CunHhw3BFXzYHgpvMswBN2KQPYRrACgG6FdhiHoVgRqA8EKNS104+f+vlkz+ia0yzAE3YpAbSBYoaaFbvzMZs2oRiELibKIaJVg8dKaQbAC8IaQta56I3RdrJBZhlnYRifmeC1Jagwo09G5R1pZV7xgU7N02KL9axC6x+KlNYNgBeANIWtd9UbIulghswyrfRudEDHHa0W3eU3yL8EK2G8EKwAVFTLLsNpnDoaIOV6rN7bu3KX5hzT0XOiei8rTGKAGEKwAAMlZq5gBq6mZM2CoSQQrACUTOmarv+9R2O81Ncetj65F1DCCFRAgdFmGELW0dEPomK1a2qOwKh22KG4IomsRNYxgBQQIXZYhBEs3oCaEdC02NXNWC5lDsAIAxNXUXLwM3YXIKIIVAGRYyCKiUuSFREO6FukuLL3QRUlZSDQqghUAZNj4EfVB5bbu3FXilqDsQhclZSHRqAhWAPqFkNXZeyMLK7kDqD4EKwBVL2R19t7IwkruQDR0GUZFsAJQ9UJWZ++NLKzkDkRDl2FUBCsAQGUGuQMZRLACADDIHcWFdBnSXUiwAlB5IVvfsO0NUGEhXYZ0FxKsAFReyNY3sbe9CZ1lyOzBEmLjZ2QQwQpAzQmdZcjswRJqao5bHyu5VwdmGBKsgHIL3dC5ljZrLrfQWYbMHiwhNn7OJmYYEqyAcgvd0JnNmvcWMg5LYiwWgMoiWAHoF0LGYUnxx2IBKIHQLsNQVdS1SLACgB7E3EqHgfBAKrTLMFQVdS0SrACgGzG30mEgPFAbCFYA0I2YW+kwEL4MYi/fEKKpmZmI2AvBCgDQ/zU1l/8xWeIBBRCsAAD9X+zlG0KwxAMKIFgBANBXod2PTc2c2aoRBCsAKBO20cmYpuawcnQZ1hSCFQCUAdvoZFBo9yNdhjWFYAUAZcA2OkBtIFgBVYo9BQGg/yFYAVWKPQX7hj0FAVQSwQpApmRhT8FqHuS+q3OPfvXspqLlDhgySHOmjCpDi4DqQrACgCpS7YPcx4+oDyq3deeuErekn2FV+JpBsAJQk6q1y5BB7hnU1Fz+x2SJh4ohWAGoSVnoMkQ/warwNWVApRsAAACQFZyxAvo5lmUAgOpBsAL6OZZlKK3QsVghWOIByD6CFQD0IHQsVgjGawHZxxgrAACASAhWAAAAkdAVCABlEnO8lnWOlCTVtbVGqS/2+K+QFdpZnR1ZRLACakTo7MEQzDDsm5jjtXavXKl1LTv12ZVjeyw3f3qDFs6cULS+2OO/QlZoZ3X2Egtd7b2pmYVEIyJYATUidPZgCGYYVt786Q2Sej4jtK5lp6RNQcEKGdPUHFaOFdqjI1gBQD+0cOaEooFp6Z0ry9QaVJ3Q1d5ZoT06ghUAIOr4L4k1u/oVugyjIlgBAKKO/5LCxmyFDHCXGOReUk3NYeXoMgxGsAIAVETIAHeJQe4lRZdhdAQrAL0Wc4ahxCzDLArpWqS7EFlEsALQazFnGErMMsyikK5FtvjpZ0LGYjU113x3IcEKAFDVQsdihWLMVh80NRcvwzgsSQQrAECVCx2LFYoxW30QMhaLcViSCFYAACAWlm4gWAFAlq1r2Rm0UGjo1jdAt5qaw8plvMuQYAWg4kJmGTJzsPdCtr2R2PoGkbB0gySCFYAqEDLLkJmDvRey7Y3E1jdATAQrAABQXqFjsUIddKx01Jnx6tsPBCsAAFA+Tc1x69u8RurcHbfO/UCwAgAwyB3lEzoWK9Q9FxGsAADVg0HuQDwEKwCocZUa5B6yn2BvsPcgqkFQsDKzRZK+LWmgpO+5+1V5x6dK+oGkUWmZi9z9rshtBQBkSMh+gr3B3oOoBgOKFTCzgZK+I+lUSUdIWmxmR+QV+4qkW9x9jqQPSbo2dkMBAACqXdFgJWmepGfdfbW775L0E0mn55VxSQeml0dKeileEwEAAPqHkK7ASZJeyLm+XtJxeWUuk3SvmV0gaZikU6K0DgAAoJCbPyq1b33z+mUjk3+HjZe+9IfKtElhZ6yswG2ed32xpBvdfbKk90j6oZntU7eZnWNmK8xsxcaN8QYsAgCAGpMbqnLteLW87cgTEqzWS8qdZjFZ+3b1fVzSLZLk7o9IqpfUkF+Ru1/n7nPdfe64ceP61mIAAIAqFdIVuFzSoWY2TdKLSganfzivzPOSFkq60cxmKglWnJICEE3IRs29wabO2RO6fMOQ9t3S9s44D7q7XRp9cJy6kAlFg5W7d5jZ+ZJ+oWQphRvcfaWZLZW0wt2XSfqCpOvN7EIl3YRL3D2/uxAA+ixko+beYFPn7AldvmHjtnbpqElxHvTJW+PUg8wIWscqXZPqrrzbLsm5vErS/LhNAwAA6F9YeR0AUFPaOzqjnbGcvuE1zTqMLX4qon5U4QHsw+Ke3e4tghUAoKbEHFv32itDpe0bihdkLFZ8H/zRm5e3b5COOrNybclBsAIAoI82TlwQNl6LsVg1g2AFoCaFzjJk9iCA3iBYAahJobMMmT0IoDcIVgAA9FHoQHgGudcOghUAAH0U2k3c9kKkBUlR9UK2tAEAAEAAzlgBAFBiu+uGhy3LEPVBWeKhEghWANCDmHsUMsOwdgUvyxATSzxUBMEKAHoQc4/CLMwwXNeyU0vvXFm03PzpDVo4k8HaqD0EKwBAkPnTGyRtKlpuXctOSZsIVqhJBCsAQJCFMycEhaWQM1pAVjErEAAAIBKCFQAAQCR0BQJAmbA/Ye0KXaF9xNC6qBMmUH4EKwAoE/YnrF2hQTnW0h6oHLoCAQAAIiFYAQAAREJXIAAAWVQ/MmwbHba+iYpgBQBAFh36zrBybH0TFcEKAKpMFmYPhmx9w7Y3yCKCFQBUmf4+ezBk6xu2vUFWEawAAFGFbH3DtjdVhLFYURGsAACoZYzFiopgBQBAlQhZoZ3V2asbwQoAUBEhA9yl2hrkHjIZoWKrs4d0GdJdSLACAJRfyAB3iUHuVSWky5DuQoIVAKD8Qga4SwxyR//DljYAAACRcMYKAADEwdINBCsAABAJSzcQrAAA1S109mCoWppliPIjWAFAPxWyp2A17ycYInT2YChmGVaJ0C7DUIOr5zNOsAKAfipkkchq3U8wVOjswVDMMqwSoV2GodY8GLe+/UCwAgCgHwlZnV1ihfZKIVgBANCPhHbtVmyF9hrHOlYAAACREKwAAAAiIVgBAABEwhgrAAAyKHSQeygGw4chWAEAkEGx1y9jMHwYugIBAAAiIVgBAABEQrACAACIhDFWAJBhIfsJSv1/T0GUXshg+IoNcJ92QvkfsxsEKwDIsNBfcv19T0GUXkjwZoA7XYEAAADREKwAAAAiIVgBAABEQrACAACIhMHrAABmDwKREKwAAMweBCIhWAEAasq6lp1aeufKKHXNn96ghTMnRKkL2UCwAgDUjPnTGyRtilLXupadkjYRrLAXghUAoGYsnDkhWhCKddYL2cKsQAAAgEgIVgAAAJEQrAAAACJhjBUAAH0UOsOQ2YO1g2AFAEAfhM4wZPZgbSFYAQDQB6EzDJk9WFsYYwUAABAJwQoAACASghUAAEAkjLECAABRtHd0Bm3UPWJoXfDG3/0NwQoAAEQxZfSwoHIbt7WXuCWVQ1cgAABAJAQrAACASOgKBACgxFihvXYQrAAAKCFWaK8tBCsAAEqIFdprC2OsAAAAIiFYAQAAREJXIAAg2IihdUFrELV3dAavaQRkCcEKABAsdLXskNW3UbtCV2gPVU0ruROsAABAWcU+m1lNK7kzxgoAACASghUAAEAkBCsAAIBIgsZYmdkiSd+WNFDS99z9qrzjfy/ppPTqAZLGu/uomA0FACDrQra+Ydub6lY0WJnZQEnfkfROSeslLTezZe6+qquMu1+YU/4CSXNK0FYAADIrZOsbtr2pfiFnrOZJetbdV0uSmf1E0umSVnVTfrGkS+M0DwCA2hCy9Q3b3lS/kDFWkyS9kHN9fXrbPszsYEnTJP1y/5sGAADQv4QEKytwm3dT9kOSbnX3zoIVmZ1jZivMbMXGjRtD2wgAANAvhASr9ZKm5FyfLOmlbsp+SNJN3VXk7te5+1x3nztu3LjwVgIAAPQDIWOslks61MymSXpRSXj6cH4hM5shabSkR6K2EADQ74TsKch+gsiiosHK3TvM7HxJv1Cy3MIN7r7SzJZKWuHuy9KiiyX9xN276yYEANSIkH3b2E8QWRS0jpW73yXprrzbLsm7flm8ZgEAAPQ/rLwOAAAQCcEKAAAgkqCuQAAAUB1Ctr2R2PqmUghWAICKCJk5KDF7MFfItjcSW99UEsEKAFARITMHJWYP5grZ9kZi65tKYowVAABAJAQrAACASAhWAAAAkTDGCgCADAqdPRiKWYZhCFYAgKoWOnswVC3MMgydPRiKWYbhCFYAgKoWOnswVC3MMgydPRiKWYbhGGMFAAAQCcEKAAAgEroCAQBAUSGD4RngTrACAABFhAyGZ4B7gmAFAAB6FDIYngHuCcZYAQAAREKwAgAAiIRgBQAAEAnBCgAAIBKCFQAAQCQEKwAAgEgIVgAAAJEQrAAAACIhWAEAAERCsAIAAIiEYAUAABAJwQoAACASghUAAEAkBCsAAIBICFYAAACREKwAAAAiGVTpBgAAUE4jhtZp47b2KHW1d3RqyuhhUepCNhCsAAA15aQZ46PVtezxF6PVhWygKxAAACASzlgBAIAo1rXs1NI7VxYtN396gxbOnFCGFpUfwQoAAOy3+dMbJG0qWm5dy05JmwhWAAAA3Vk4c0JQWAo5o9WfEawAAOij0BmGzB6sHQQrAAD6KHSGIbMHawezAgEAACIhWAEAAERCsAIAAIiEYAUAABAJwQoAACASghUAAEAkLLcAAADKKnTrm1BHTxqp982eFK2+/UGwAgAAZRO69U2odS07tbtjT7T69hfBCgAAlE3o1jehlt65kmAFAEAtYeub+P7yR4+qtW33G9cbL/p3SVLD8MFa8ZV3VqpZBCsAAEqNrW/iyw1VuTZt31XmluyNWYEAAACREKwAAAAiIVgBAABEQrACAACIhGAFAAD6nZFD6wre3jB8cJlbsjdmBQIAgH7nux99yxuXN25rr5qV1zljBQAAEAnBCgAAIBKCFQAAQCSMsQIAoEqEbH3DtjfVjWAFAECVCNn6hm1vqhtdgQAAAJEQrAAAACIhWAEAAERCsAIAAIiEYAUAABAJwQoAACASghUAAEAkBCsAAIBICFYAAACREKwAAAAiYUsbAAD6kZD9BCX2FKwUghUAAP1IyH6CEnsKVgpdgQAAAJEQrAAAACIhWAEAAETCGCsAADIodJB7KAbDhyFYAQCQQaGD3EMxGD4MXYEAAACREKwAAAAiCQpWZrbIzH5vZs+a2UXdlPmAma0ys5Vm9uO4zQQAAKh+RcdYmdlASd+R9E5J6yUtN7Nl7r4qp8yhkr4sab67bzGzuB27AAAA/UDI4PV5kp5199WSZGY/kXS6pFU5ZT4p6TvuvkWS3P3V2A0FAACVEzLLkJmDYcFqkqQXcq6vl3RcXpnDJMnMfiVpoKTL3P2e/IrM7BxJ50jS1KlT+9JeAABQASGzDJk5GDbGygrc5nnXB0k6VFKzpMWSvmdmo/a5k/t17j7X3eeOGzeut20FAACoaiHBar2kKTnXJ0t6qUCZ2919t7uvkfR7JUELAACgZoQEq+WSDjWzaWY2WNKHJC3LK/NzSSdJkpk1KOkaXB2zoQAAANWuaLBy9w5J50v6haSnJd3i7ivNbKmZvS8t9gtJLWa2StJ9kr7k7i2lajQAAEA1CtrSxt3vknRX3m2X5Fx2SZ9PfwAAAGoSK68DAABEQrACAACIhGAFAAAQCcEKAAAgEoIVAABAJAQrAACASAhWAAAAkRCsAAAAIiFYAQAAREKwAgAAiIRgBQAAEAnBCgAAIBKCFQAAQCQEKwAAgEgIVgAAAJEMqnQDAABANowYWqeN29qLlmvv6NSU0cPK0KLyI1gBAIAoTpoxPqjcssdfLHFLKoeuQAAAgEgIVgAAAJEQrAAAACIhWAEAAERCsAIAAIiEYAUAABAJwQoAACASghUAAEAkBCsAAIBIWHkdAACUVejWN6EOGFI9caZ6WgIAAGpC6NY3oR55riVqffuDrkAAAIBICFYAAACREKwAAAAiIVgBAABEQrACAACIhGAFAAD6tbdPH1vpJryBYAUAABAJwQoAACASghUAAEAkBCsAAIBICFYAAACREKwAAAAiIVgBAABEQrACAACIhGAFAAAQCcEKAAAgEoIVAABAJAQrAACASAhWAAAAkRCsAAAAIiFYAQAAREKwAgAAiIRgBQAAEAnBCgAAIBJz98o8sNlGSetK/DANkjaV+DFQHO9DdeB9qDzeg+rA+1B5/fE9ONjdxxUrVLFgVQ5mtsLd51a6HbWO96E68D5UHu9BdeB9qLwsvwd0BQIAAERCsAIAAIgk68Hquko3AJJ4H6oF70Pl8R5UB96Hysvse5DpMVYAAADllPUzVgAAAGWT2WBlZovM7Pdm9qyZXVTp9tQKM7vBzF41s6dybhtjZv9hZn9I/x1dyTZmnZlNMbP7zOxpM1tpZp9Nb+d9KCMzqzez35jZ79L34fL09mlm9j/p+3CzmQ2udFuzzswGmtljZnZnep33oMzMbK2ZPWlmj5vZivS2TH4nZTJYmdlASd+RdKqkIyQtNrMjKtuqmnGjpEV5t10k6b/c/VBJ/5VeR+l0SPqCu8+U9DZJn04//7wP5fW6pJPd/RhJsyUtMrO3SfqGpL9P34ctkj5ewTbWis9KejrnOu9BZZzk7rNzllnI5HdSJoOVpHmSnnX31e6+S9JPJJ1e4TbVBHd/QNLmvJtPl/SD9PIPJP1JWRtVY9z9ZXf/bXp5m5JfKJPE+1BWntieXq1Lf1zSyZJuTW/nfSgxM5ss6b2SvpdeN/EeVItMfidlNVhNkvRCzvX16W2ojAnu/rKU/NKXNL7C7akZZtYoaY6k/xHvQ9mlXVCPS3pV0n9Iek7SVnfvSIvw3VR635L0V5L2pNfHivegElzSvWb2qJmdk96Wye+kQZVuQIlYgduY/oiaYmbDJd0m6XPu/lryhzrKyd07Jc02s1GSfiZpZqFi5W1V7TCz0yS96u6Pmllz180FivIelN58d3/JzMZL+g8ze6bSDSqVrJ6xWi9pSs71yZJeqlBbIG0wsz+SpPTfVyvcnswzszoloepf3f2n6c28DxXi7lsl3a9kzNsoM+v6o5bvptKaL+l9ZrZWyZCQk5WcweI9KDN3fyn991Ulf2TMU0a/k7IarJZLOjSd+TFY0ockLatwm2rZMklnpZfPknR7BduSeekYku9Letrd/2/OId6HMjKzcemZKpnZUEmnKBnvdp+kM9NivA8l5O5fdvfJ7t6o5PfAL939I+I9KCszG2ZmI7ouS3qXpKeU0e+kzC4QambvUfKXyUBJN7j7lRVuUk0ws5skNSvZuXyDpEsl/VzSLZKmSnpe0p+5e/4Ad0RiZu+Q9KCkJ/XmuJK/VjLOivehTMzsaCUDcgcq+SP2FndfamZNSs6ejJH0mKSPuvvrlWtpbUi7Ar/o7qfxHpRX+nr/LL06SNKP3f1KMxurDH4nZTZYAQAAlFtWuwIBAADKjmAFAAAQCcEKAAAgEoIVAABAJAQrAACASAhWACrKzEaZ2afSyweZ2a3F7hNY72Vm9sX08lIzOyVGvQDQE5ZbAFBR6X6Gd7r7kZHrvUzSdne/Oma9ANATzlgBqLSrJE03s8fN7N/M7ClJMrMlZvZzM7vDzNaY2flm9nkze8zMfm1mY9Jy083snnRz1wfN7PD8BzCzG83szPTyWjO73Mx+a2ZPdpVPV4e+wcyWp49xehlfAwAZQbACUGkXSXrO3WdL+lLesSMlfVjJvmJXStrp7nMkPSLpY2mZ6yRd4O5vkfRFSdcGPOYmdz9W0v9L7yNJf6Nky5O3SjpJ0t+m228AQLBBxYsAQMXc5+7bJG0zs1ZJd6S3PynpaDMbLul4Sf+WbJEoSRoSUG/XxtSPSvrT9PK7lGzY2xW06pVstfH0/j0FALWEYAWgmuXu37Yn5/oeJd9fAyRtTc929aXeTr35PWiS3u/uv+9jWwGArkAAFbdN0oi+3NHdX5O0xsz+TJIscUwf2/ELSRdYeurLzOb0sR4ANYxgBaCi3L1F0q/SQet/24cqPiLp42b2O0krJfV10PkVkuokPZG25Yo+1gOghrHcAgAAQCScsQIAAIiEYAUAABAJwQoAACASghUAAEAkBCsAAIBICFYAAACREKwAAAAiIVgBAABE8v8BboRkcupv9vUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "## Do stratification for the example data set.\n",
    "#variable of stratification : age, time : week, event: arrest, 10 bootstrap iterations\n",
    "rossi = load_rossi()\n",
    "stratify_plot(rossi,\"age\",\"week\",\"arrest\",10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}