{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from lifelines.datasets import load_rossi\n",
"rossi = load_rossi()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" week | \n",
" arrest | \n",
" fin | \n",
" age | \n",
" race | \n",
" wexp | \n",
" mar | \n",
" paro | \n",
" prio | \n",
"
\n",
" \n",
" \n",
" \n",
" 0 | \n",
" 20 | \n",
" 1 | \n",
" 0 | \n",
" 27 | \n",
" 1 | \n",
" 0 | \n",
" 0 | \n",
" 1 | \n",
" 3 | \n",
"
\n",
" \n",
" 1 | \n",
" 17 | \n",
" 1 | \n",
" 0 | \n",
" 18 | \n",
" 1 | \n",
" 0 | \n",
" 0 | \n",
" 1 | \n",
" 8 | \n",
"
\n",
" \n",
" 2 | \n",
" 25 | \n",
" 1 | \n",
" 0 | \n",
" 19 | \n",
" 0 | \n",
" 1 | \n",
" 0 | \n",
" 1 | \n",
" 13 | \n",
"
\n",
" \n",
" 3 | \n",
" 52 | \n",
" 0 | \n",
" 1 | \n",
" 23 | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
"
\n",
" \n",
" 4 | \n",
" 52 | \n",
" 0 | \n",
" 0 | \n",
" 19 | \n",
" 0 | \n",
" 1 | \n",
" 0 | \n",
" 1 | \n",
" 3 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" week arrest fin age race wexp mar paro prio\n",
"0 20 1 0 27 1 0 0 1 3\n",
"1 17 1 0 18 1 0 0 1 8\n",
"2 25 1 0 19 0 1 0 1 13\n",
"3 52 0 1 23 1 1 1 1 1\n",
"4 52 0 0 19 0 1 0 1 3"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rossi.head()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# let's b-spline age\n",
"cph = CoxPHFitter().fit(rossi, \"week\", \"arrest\", formula=\"fin + bs(age, df=4) + wexp + mar + paro + prio\")"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" week | \n",
" arrest | \n",
" fin | \n",
" age | \n",
" race | \n",
" wexp | \n",
" mar | \n",
" paro | \n",
" prio | \n",
"
\n",
" \n",
" \n",
" \n",
" 0 | \n",
" 52.0 | \n",
" 0.0 | \n",
" 0.5 | \n",
" 17.000000 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 0.0 | \n",
" 1.0 | \n",
" 2.0 | \n",
"
\n",
" \n",
" 1 | \n",
" 52.0 | \n",
" 0.0 | \n",
" 0.5 | \n",
" 17.551020 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 0.0 | \n",
" 1.0 | \n",
" 2.0 | \n",
"
\n",
" \n",
" 2 | \n",
" 52.0 | \n",
" 0.0 | \n",
" 0.5 | \n",
" 18.102041 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 0.0 | \n",
" 1.0 | \n",
" 2.0 | \n",
"
\n",
" \n",
" 3 | \n",
" 52.0 | \n",
" 0.0 | \n",
" 0.5 | \n",
" 18.653061 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 0.0 | \n",
" 1.0 | \n",
" 2.0 | \n",
"
\n",
" \n",
" 4 | \n",
" 52.0 | \n",
" 0.0 | \n",
" 0.5 | \n",
" 19.204082 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 0.0 | \n",
" 1.0 | \n",
" 2.0 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" week arrest fin age race wexp mar paro prio\n",
"0 52.0 0.0 0.5 17.000000 1.0 1.0 0.0 1.0 2.0\n",
"1 52.0 0.0 0.5 17.551020 1.0 1.0 0.0 1.0 2.0\n",
"2 52.0 0.0 0.5 18.102041 1.0 1.0 0.0 1.0 2.0\n",
"3 52.0 0.0 0.5 18.653061 1.0 1.0 0.0 1.0 2.0\n",
"4 52.0 0.0 0.5 19.204082 1.0 1.0 0.0 1.0 2.0"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# now we need to \"extend\" our data to plot it\n",
"# we'll plot age over it's observed range\n",
"age_range = np.linspace(rossi['age'].min(), rossi['age'].max(), 50)\n",
"\n",
"# need to create a matrix of variables at their means, _except_ for age. \n",
"x_bar = cph._central_values\n",
"df_varying_age = pd.concat([x_bar] * 50).reset_index(drop=True)\n",
"df_varying_age['age'] = age_range\n",
"\n",
"df_varying_age.head()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAloUlEQVR4nO3deXxU9b3/8dcnCUlIgCQkIWEJJAiyyCoRXHABl+JS0YoWqtVWre1t6Xq7eOv93fZ621uX22p7a1tx12oRtVZaF4qKCyhKkF32sCUCCYQlEEhI+Pz+mIEbQ5Blkpxk5v18POYxc77nnJnPeTDkPed8zzlfc3dERCR2xQVdgIiIBEtBICIS4xQEIiIxTkEgIhLjFAQiIjEuIegCTkZWVpbn5+cHXYaISJsyf/78be6e3bC9TQZBfn4+RUVFQZchItKmmNmGxtp1aEhEJMYpCEREYpyCQEQkxikIRERiXJMEgZmNM7OVZrbGzG5vZP59ZrYw/FhlZjvrzaurN296U9QjIiLHL+KzhswsHngAuBgoAeaZ2XR3//jQMu7+/XrLfxsYXu8t9rn7sEjrEBGRk9MUewQjgTXuXuzuNcBUYPxnLD8J+EsTfK6IiDSBpgiC7sCmetMl4bYjmFkvoAB4s15zspkVmdlcM7uqCeo5qpcXb+bPcxs9jVZEJGa19AVlE4Hn3b2uXlsvdy81s97Am2a2xN3XNlzRzG4DbgPo2bPnSX34K0s2M2ftNiaM6EFyu/iTeg8RkWjTFHsEpUBeveke4bbGTKTBYSF3Lw0/FwNv8en+g/rLTXH3QncvzM4+4grp4zJpZE92Vh1gxrItJ7W+iEg0aoogmAf0NbMCM0sk9Mf+iLN/zKw/kAG8X68tw8ySwq+zgHOAjxuu21TOPiWTnp1T+MuHG5vrI0RE2pyIg8Dda4HJwAxgOTDN3ZeZ2Z1mdmW9RScCU/3TY2MOAIrMbBEwC7ir/tlGTS0uzvjiGXnMLa6guHxPc32MiEibYm1xzOLCwkI/2ZvOlVXu5+xfvcnNowv46WUDmrgyEZHWy8zmu3thw/aYu7K4S8dkLhqQw/PzS6iurTv2CiIiUS7mggBg0qieVOytYebHW4MuRUQkcDEZBOf2yaJ7ent1GouIEKNBEBdnTDwjjzlrtrNh+96gyxERCVRMBgHAtYV5xMcZU+dtOvbCIiJRLGaDIDctmbH9u/Bc0SZqag8GXY6ISGBiNggAJo3MY9ueGt5Yrk5jEYldMR0E55/aha5pyfxFh4dEJIbFdBDEh680fnd1OZsqqoIuR0QkEDEdBADXFeZhwLPaKxCRGBXzQdAtvT0X9OvCtKJNHKhTp7GIxJ6YDwKAL5/Zi7LKal5ZsjnoUkREWpyCADj/1Gx6Z6Xy6Ox1tMWb8ImIREJBQOhK46+ek8+ikl18tHFH0OWIiLQoBUHYF07vQafkBB6dvT7oUkREWpSCICw1KYFJo3ry6tLNlOzQqaQiEjsUBPXceFY+ZsaT728IuhQRkRajIKine3p7xg3K5S8fbmRvdW3Q5YiItAgFQQO3jC6gcn8tL3xUEnQpIiItQkHQwOk9MxiWl85jc9Zz8KBOJRWR6KcgaMTNowtYt20vs1aWBV2KiEiza5IgMLNxZrbSzNaY2e2NzP+KmZWb2cLw49Z6824ys9Xhx01NUU+kLh2US9e0ZB6dsy7oUkREml3EQWBm8cADwKXAQGCSmQ1sZNFn3X1Y+PFweN3OwM+AUcBI4GdmlhFpTZFqFx/HjWflM2fNdpZv3h10OSIizaop9ghGAmvcvdjda4CpwPjjXPdzwEx3r3D3HcBMYFwT1BSxSSPzaN8unse0VyAiUa4pgqA7UP8eziXhtoauMbPFZva8meWd4LqY2W1mVmRmReXl5U1Q9mdLT0nkmhHd+dvCTyivrG72zxMRCUpLdRb/Hch39yGEfvU/caJv4O5T3L3Q3Quzs7ObvMDG3HxOAQfqDqqvQESiWlMEQSmQV2+6R7jtMHff7u6HflY/DIw43nWD1Du7A5cP7sqT761nZ1VN0OWIiDSLpgiCeUBfMysws0RgIjC9/gJm1rXe5JXA8vDrGcAlZpYR7iS+JNzWakwe24e9NXU8Nmd90KWIiDSLiIPA3WuByYT+gC8Hprn7MjO708yuDC/2HTNbZmaLgO8AXwmvWwH8F6EwmQfcGW5rNfrnduKSgTk8NmcdlfsPBF2OiEiTs7Y4EEthYaEXFRW12OctKdnF538/mx99rh/fGtOnxT5XRKQpmdl8dy9s2K4ri4/D4B5pXNAvm0dmr6OqRjejE5HooiA4Tt8e24eKvTU888HGoEsREWlSCoLjNKJXZ87qncmUd4rZf6Au6HJERJqMguAEfPvCPpRVVvNc0aZjLywi0kYoCE7AWb0zGdErgz+9XUxN7cGgyxERaRIKghNgZkwe24fSnfv424JWc92biEhEFAQn6IJTsxncPY0H3lpDbZ32CkSk7VMQnKBDewUbtlfx0sJPgi5HRCRiCoKTcPGAHAZ178R9r6+iulZnEIlI26YgOAlxccbt4wZQsmMff56r6wpEpG1TEJyk0X2zOLdvFr9/czW7dQ8iEWnDFAQR+Mm4/uyoOsCUt4uDLkVE5KQpCCIwqHsaVw7txsOziynbvT/ockREToqCIEI/vKQfdQed+99YHXQpIiInRUEQoZ6ZKVw/qhfPztvE2vI9QZcjInLCFARNYPLYPiQnxPE/M1YGXYqIyAlTEDSBrA5J3HbeKby6dAsfbdwRdDkiIidEQdBEbj23gKwOSdz16gra4qhvIhK7FARNJDUpge9e1JcP11Uwa2VZ0OWIiBw3BUETmnhGHgVZqfzi5eW6TbWItBlNEgRmNs7MVprZGjO7vZH5PzCzj81ssZm9YWa96s2rM7OF4cf0pqgnKO3i4/jZ5wdSXL6XR2avC7ocEZHjEnEQmFk88ABwKTAQmGRmAxsstgAodPchwPPAPfXm7XP3YeHHlZHWE7QL+nXhkoE5/O6N1Xyyc1/Q5YiIHFNT7BGMBNa4e7G71wBTgfH1F3D3We5eFZ6cC/Rogs9ttf7fFQNxnF+8/HHQpYiIHFNTBEF3oP4gviXhtqO5BXi13nSymRWZ2Vwzu+poK5nZbeHlisrLyyMquLnldU5h8pg+vLJkC++sat21ioi0aGexmd0AFAL31mvu5e6FwJeA+83slMbWdfcp7l7o7oXZ2dktUG1kvnZeb/IzU/j59GUas0BEWrWmCIJSIK/edI9w26eY2UXAHcCV7l59qN3dS8PPxcBbwPAmqClwSQnx/PzK0yjetpeH31XHsYi0Xk0RBPOAvmZWYGaJwETgU2f/mNlw4EFCIVBWrz3DzJLCr7OAc4CoObB+Qb8ujDstl9+/uYZSdRyLSCsVcRC4ey0wGZgBLAemufsyM7vTzA6dBXQv0AF4rsFpogOAIjNbBMwC7nL3qAkCgP/3+XDH8T+iarNEJIokNMWbuPsrwCsN2v6j3uuLjrLee8Dgpqihteqe3p5vj+3LvTNW8vaqcs4/tfX3b4hIbNGVxS3g1nMLKMhK5T9eWsq+GnUci0jroiBoAUkJ8fz31YPZsL2Ke2asCLocEZFPURC0kLNOyeQrZ+fz2Jz1zC3eHnQ5IiKHKQha0I/H9aNXZgo/fn4xe6trgy5HRARQELSolMQE7p0wlE07qrj7NR0iEpHWQUHQwkYWdOarZxfw5PsbeG/NtqDLERFREAThR5/rR0FWKj9+YTF7dIhIRAKmIAhA+8R4/ufaIZTu3MevXlkedDkiEuMUBAEZ0aszXzu3N09/sJF3V+sOpSISHAVBgH5w8amckp3KT55fzK6qA0GXIyIxSkEQoOR28fzmumGU76nmh88vwt2DLklEYpCCIGBD89K5/dIBzPx4K4/OWR90OSISgxQErcDN5+RzycAcfvXKchZs3BF0OSISYxQErYCZce+EoeSmJTP5mQXsrKoJuiQRiSEKglYiLaUdD3zpdMoq9/PD5xarv0BEWoyCoBUZmpfOTy8bwOvLt/LIbA1vKSItQ0HQynzl7Hw+d1oOd726go/UXyAiLUBB0MqYGfdMGErX9GQmP/0RFXvVXyAizUtB0AqltQ/1F2zbW8M3/jyfmtqDQZckIlFMQdBKDemRzr0ThvDhugrueHGJOo9FpNk0SRCY2TgzW2lma8zs9kbmJ5nZs+H5H5hZfr15/xZuX2lmn2uKeqLF+GHd+c6FfXlufglT3ikOuhwRiVIRB4GZxQMPAJcCA4FJZjawwWK3ADvcvQ9wH3B3eN2BwETgNGAc8Ifw+0nY9y7sy+VDunLXayuYsWxL0OWISBRqij2CkcAady929xpgKjC+wTLjgSfCr58HLjQzC7dPdfdqd18HrAm/n4TFxRm/vnYoQ3qk872pC1lauivokkQkyjRFEHQHNtWbLgm3NbqMu9cCu4DM41wXADO7zcyKzKyovDy2btuc3C6eh748gvSUdtz6RBFbd+8PuiQRiSJtprPY3ae4e6G7F2ZnZwddTovr0imZh28qZPf+A3ztySL21dQFXZKIRImmCIJSIK/edI9wW6PLmFkCkAZsP851Jey0bmn8duJwlpTuYvIzH3GgTqeVikjkmiII5gF9zazAzBIJdf5Ob7DMdOCm8OsJwJseOh9yOjAxfFZRAdAX+LAJaopaFw/M4b/GD+KNFWX8YNoi6g7qtFIRiUxCpG/g7rVmNhmYAcQDj7r7MjO7Eyhy9+nAI8BTZrYGqCAUFoSXmwZ8DNQC33J3HfM4hhvO7EXl/lrufm0FHZLi+e+rBxPqexcROXHWFi9UKiws9KKioqDLCNy9M1bwwKy1fO3cAn562QCFgYh8JjOb7+6FDdsj3iOQ4Pzwkn7s2V/LQ++uo2NyO75zYd+gSxKRNkhB0IaZGT/7/GlUVtfym5mr6JCUwM2jC4IuS0TaGAVBGxcXZ9xzzRCqquu48x8fk5IYz8SRPYMuS0TakDZzHYEcXUJ8HL+dNIzzT83m9r8u4fE5GtRGRI6fgiBKJCXEM+XGEVwyMIef//1jHpi1JuiSRKSNUBBEkaSEeP5w/elcPbw7985YyV2vrtDtq0XkmNRHEGUS4uP49bVDSUmM509vr2VvdS3/eeVpxMXp1FIRaZyCIArFxRm/uGoQHZISePCdYvbW1HLPNUNIiNcOoIgcSUEQpcyM2y/tT4ekBH49cxV7q2u5/4vDaZ+o4R5E5NP0EzGKmRnfvrAvP/v8QP758VYmPjSXskrdwlpEPk1BEAO+ek4BD94wglVbKrn6gfdYsWV30CWJSCuiIIgRl5yWy3PfOIvagweZ8Mf3eWtlWdAliUgroSCIIYO6p/HSt0bTKzOFmx+fx5Pvrw+6JBFpBRQEMSY3LZlpXz+Lsf1z+I+XlvHz6cuo1QA3IjFNQRCDUpMSePDLI7h1dAGPv7eeLz30AWUaB1kkZikIYlR8nPHvVwzk/i8OY0npLi773WzeW7st6LJEJAAKghh31fDuTJ98DmntE7jh4Q94YNYaDmr4S5GYoiAQ+uZ0ZPrk0Vw+pBv3zljJrU8WsbOqJuiyRKSFKAgECPUb/G7iMO4cfxrvri7n8t/NZv6GiqDLEpEWoCCQw8yMG8/K57lvnI0ZXPun97nntRXU1OqsIpFoFlEQmFlnM5tpZqvDzxmNLDPMzN43s2VmttjMvlhv3uNmts7MFoYfwyKpR5rGsLx0Xv3uuVw7Io8/vLWW8Q/MYeWWyqDLEpFmEukewe3AG+7eF3gjPN1QFXCju58GjAPuN7P0evN/5O7Dwo+FEdYjTaRjcjvunjCEh24spLxyP5//39lMeWctdepIFok6kQbBeOCJ8OsngKsaLuDuq9x9dfj1J0AZkB3h50oLuXhgDjO+dx5j+mfz36+sYNJDc9mwfW/QZYlIE4o0CHLcfXP49RYg57MWNrORQCKwtl7zL8OHjO4zs6TPWPc2Mysys6Ly8vIIy5YTkdkhiT/dMIJfXzuU5Z/s5pL73uH3b66murYu6NJEpAnYsYYyNLPXgdxGZt0BPOHu6fWW3eHuR/QThOd1Bd4CbnL3ufXathAKhynAWne/81hFFxYWelFR0bEWk2awZdd+/usfH/Pyks2ckp3KL64azFmnZAZdlogcBzOb7+6FDduPuUfg7he5+6BGHi8BW8N/zA/9UW/0lpZm1gl4GbjjUAiE33uzh1QDjwEjT27zpKXkpiXzwPWn89hXz6Cm7iCTHprLD6YtZNue6qBLE5GTFOmhoenATeHXNwEvNVzAzBKBF4En3f35BvMOhYgR6l9YGmE90kLG9OvCzO+fz+Qxffj7ok+48Ndv89T76zmgG9iJtDnHPDT0mSubZQLTgJ7ABuA6d68ws0LgG+5+q5ndQOjX/rJ6q37F3Rea2ZuEOo4NWBheZ8+xPleHhlqXNWWV/PvfljK3uILe2an8ZFx/LhmYQyjfRaS1ONqhoYiCICgKgtbH3Xl9eRl3vbqcteV7OSM/g3+7bACn92y0y0hEAnDSfQQix8PMDp9q+surB7FuWxVf+MN7fPPp+azbptNNRVoz7RFIs9hbXctD7xYz5Z1iqmsPMn5oN745pg99unQIujSRmKVDQxKIssr9PPh2MU9/sIHq2oNcNrgrk8f0YUDXTkGXJhJzFAQSqO17qnlk9jqefH8De6pruWhADpPH9mFYXnrQpYnEDAWBtAq7qg7w+HvreXTOOnbtO0BhrwxuOjufcYNyaRevLiuR5qQgkFZlT3UtUz/cyFNzN7BhexU5nZK4flQvJo3sSXbHo95pREQioCCQVungQeetVWU8/t4G3llVTmJ8HJcP6crEM/I4I78zcXG6FkGkqRwtCBKCKEbkkLg4Y2z/HMb2z2Ft+R6een8Dz88v4cUFpeR1bs8XhvfgmtN70DMzJehSRaKW9gik1amqqWXGsi28ML+UOWu34Q4j8ztzzYjuXDq4K52S2wVdokibpEND0iZ9snMfLy4o5YX5JRRv20u7eOOcPllcOiiXiwfm0jk1MegSRdoMBYG0ae7Owk07eWXJZl5duoWSHfuIMxhVkMmlg3O5eGAOXdPaB12mSKumIJCo4e4s+2Q3ry3dwqtLN7O2PHQLi1NzOnBe32zOOzWbkQWdSW4XH3ClIq2LgkCi1uqtlcxaWcbbq8qZt24HNXUHSW4Xx6iCTM7tm8WogkwGdO1Igq5TkBinIJCYUFVTywfFFby9qpx3VpdTHN5bSE2M5/ReGYzM78wZBZ0ZlpeuPQaJOTp9VGJCSmICY/p3YUz/LkBoaM0P11cwb10FH66r4NczVwHQLt7on9uJwT3SGNI9jcE90jg1p6OubpaYpD0CiSk7q2ooWr+Dog07WFK6k8Ulu6jcXwtAUkIcA7t1YkDXTvTL6Ui/3I70y+lIhs5MkiihQ0MijTh40NlYUcXi0l0sKQkFw4otlezad+DwMtkdk+if25FTsjvQOzuVgqzQo1tae135LG2KDg2JNCIuzsjPSiU/K5Urh3YDQmcllVVWs2JLJau2VLJyayUrt1QyrWgTVTV1h9dNTIgjPzOFXpmp5GWk0COjPXmdU8jr3J4eGSl0SNJ/L2kb9E0VacDMyOmUTE6nZM4/Nftwu7tTXllN8ba9rAs/isv3sn7bXmav3sa+A3Wfep/0lHZ0TWtPt7RkctOS6ZbentxOyXRNTz78/goLaQ30LRQ5TmZGl07JdOmUzJm9Mz81z93ZvreGkh372FRRxaYdVZTu2MfmXfv5ZNd+Ptq4gx1VB454z5TEeHI6JZPdMYmcTsl06ZhEdscksjuEn8OPjJRE4nUYSppJREFgZp2BZ4F8YD1wnbvvaGS5OmBJeHKju18Zbi8ApgKZwHzgy+5eE0lNIkEwM7I6JJHVIemog+3sq6ljy+79bN65j7LKarbu3v9/z7urWVyyk22V1eytqTti3TiDzPD7Z3dMIqtD4hGB0aVjKFA6JSdgptCQ4xdRZ7GZ3QNUuPtdZnY7kOHuP2lkuT3ufsRgtWY2Dfiru081sz8Bi9z9j8f6XHUWSzTbW13Ltj3VlFeGHmWV1Yen67dv21NDTd3BI9ZPSoijS6ckcjqGDj916ZQUPhQVeu6W1p7ctGRdRxGDmuWsITNbCVzg7pvNrCvwlrv3a2S5I4LAQj9ZyoFcd681s7OAn7v75471uQoCkdDhqF37DhwOhvI91ZTtPvS8n627q9laGdrb2FNde8T6WR0SQ30Y6cl0TQt1dPcMd3bnZaSQqv6LqNNcZw3luPvm8OstQM5Rlks2syKgFrjL3f9G6HDQTnc/9A0tAbof7YPM7DbgNoCePXtGWLZI22dmpKckkp6SSN+cjp+57J7qWrbu3s/WcJ/F5p37+GTXPj7ZuZ/i8lBnd8NDUlkdEsnrnEJ+ZiqnZKfSO7sDp2R3oFdmivYmoswxg8DMXgdyG5l1R/0Jd3czO9ruRS93LzWz3sCbZrYE2HUihbr7FGAKhPYITmRdkVjXISmBDuE/5I1xd3ZWHWBjRdXhx6aKKjZsr2Ju8XZeXFB6eNk4gx4ZKfTp0oH+uR0Z0LUTA7p2JD8zVfdzaqOOGQTuftHR5pnZVjPrWu/QUNlR3qM0/FxsZm8Bw4EXgHQzSwjvFfQAShtbX0Sal5mRkZpIRmoiQxvp7N5bXcu6bXtZW76HteWh5zVb9/DOqnJqD4Z+lyUlxNEvtyMDcjsxJC+NoT3S6Z+rm/21BZEeGpoO3ATcFX5+qeECZpYBVLl7tZllAecA94T3IGYBEwidOdTo+iISvNSkBAZ1T2NQ97RPtVfX1rGmbA8rNleyYstulm+uZObyrTxbtAmA5HZxDO4eCoWheemMLOhMTqfkIDZBPkOkncWZwDSgJ7CB0OmjFWZWCHzD3W81s7OBB4GDQBxwv7s/El6/N6EQ6AwsAG5w9+pjfa46i0VaL3dnU8U+FpbsZOHGnSwq2cnS0l1U14bOcCrISuXM3p05s3cmowoyyU1TMLQU3WtIRAJzoO4gyzfv5sN1Fcwt3s4H6yoO3+yvICuVc/pkMqZfF84+JYv2ieqIbi4KAhFpNeoOOss372Zu8XbeX7ud94u3U1VTR2JCHGf1zmRMv2zG9s+hZ2ZK0KVGFQWBiLRa1bV1zFu3gzdXlPHWyjKKt4UGFOrbpQOXDe7KFUO6HvMUWTk2BYGItBnrt+1l1soyZizbwgfrKnCHfjkduXxIVy4b3JU+XRo/DVY+m4JARNqkst37eW3ZFv6xeDPz1odCYUDXTkwY0YOrh3enswYOOm4KAhFp87bu3s+rSzbz4oJSFpXsIjE+josH5nDdGXmM7pOlO7Qeg4JARKLKii27eXbeJl5cUMrOqgN0T2/PhBE9+NKonrpW4SgUBCISlapr65j58VaenbeJ2Wu2kRBnXDGkG7eMLjjiArhYpyAQkai3YfteHpuz/vCwoqMKOnPL6AIuHJCjw0YoCEQkhuzad4Bn523kifc2ULpzH70yU/j6eacwYUQPEhNi995HCgIRiTm1dQd5bdkWHnqnmEUlu+ie3p7JY/swYUQP2sXgzfAUBCISs9ydt1aVc//MVSwq2UWPjPZMHtOHa2IsEBQEIhLz3J23VpZz3+urWBwOhO9fdCpXD+9OXAz0IRwtCGInCkUk5pkZY/p34aVvncOjXykkIyWRf31uEeMfmMOH6yqCLi8wCgIRiTlmxtj+Obz0rXO474tD2banmusefJ9vPj2fjdurgi6vxSkIRCRmxcUZVw/vwZv/egHfv+hUZq0o56LfvM2vXl1O5f4DQZfXYhQEIhLz2ifG892L+jLrhxfw+aHdePDtYsb8z9v8Y/EntMV+1BOlIBARCctNS+bX1w1l+uRz6JqWzORnFnDrE0V8snNf0KU1KwWBiEgDQ3qk8+I3z+aOywbw3trtXPybt3l8zjrqDkbn3oGCQESkEQnxcXztvN788/vncXqvDH7+94+55o/vsWLL7qBLa3IKAhGRz5DXOYUnbx7J/V8cxsaKKq743Wz++NbaqNo7iCgIzKyzmc00s9Xh54xGlhljZgvrPfab2VXheY+b2bp684ZFUo+ISHMwM64a3p3Xf3A+Fw/M4e7XVnD9w3Ojpu8g0j2C24E33L0v8EZ4+lPcfZa7D3P3YcBYoAr4Z71FfnRovrsvjLAeEZFm0zk1kT9cfzr3XDOExSW7uPS37/Ly4s1BlxWxSINgPPBE+PUTwFXHWH4C8Kq7x94VGyISFcyM687I4+XvnEt+VirfeuYjfvTcIvZU1wZd2kmLNAhy3P1QHG4Bco6x/ETgLw3afmlmi83sPjNLOtqKZnabmRWZWVF5eXkEJYuIRK4gK5Xnv3EW3x7bhxc+KuHy373L0tJdQZd1Uo550zkzex3IbWTWHcAT7p5eb9kd7n5EP0F4XldgMdDN3Q/Ua9sCJAJTgLXufuexitZN50SkNflwXQXfnbqAir01/PLqwUwY0SPokhp10jedc/eL3H1QI4+XgK3hP+aH/qiXfcZbXQe8eCgEwu+92UOqgceAkSe6YSIiQRtZ0Jm/f3s0p/fM4IfPLeI/XlpKTe3BoMs6bpEeGpoO3BR+fRPw0mcsO4kGh4XqhYgR6l9YGmE9IiKByOqQxFO3jOS283rz5PsbmPTQXLbu3h90Wccl0iC4C7jYzFYDF4WnMbNCM3v40EJmlg/kAW83WP9pM1sCLAGygF9EWI+ISGAS4uP46WUD+N9Jw1m+eTdX/O9s5q1v/be31sA0IiLNYOWWSr7+VBElO/Zx5/hBfGlUz6BL0sA0IiItqV9uR16aPJrRfbP46YtLuPu1FRxspVcjKwhERJpJWvt2PHxjIV8a1ZM/vrWW7z67kOrauqDLOkJC0AWIiESzhPg4fnnVIPIyUrj7tRVs3b2fKV8eQXpKYtClHaY9AhGRZmZm/MsFp/DbicNYuHEn1/zxPTZVtJ4bLCgIRERayPhh3XnqlpFs21PD1X+Yw+KSnUGXBCgIRERa1KjembzwL2eT3C6eSVPm8kHx9qBLUhCIiLS0Pl068Pw3ziY3LZmbHvuQ2au3BVqPgkBEJAC5ack8+/WzyM9M5eYn5vHG8q2B1aIgEBEJSFaHJKbedib9czvy9afm88qSYMY2UBCIiAQoPSWRP986iqF56Ux+5iNeXFDS4jUoCEREAtYpuR1P3jySM3tn8oNpi5j64cYW/XwFgYhIK5CalMCjXzmD80/N5va/LuGF+S23Z6AgEBFpJZLbxfOnG0ZwTp9MfvT8Il5toT4DBYGISCuS3C6eh24sZHjPDL4zdQFvrfys8b6ahoJARKSVSUkMHSY6NSd0NlFzX3SmIBARaYXS2oc6kHtktOeWJ4pYtGlns32WgkBEpJXK7JDE07eeSUZqO2567ENWbqlsls9REIiItGK5ack8c+uZJCXEcf3DH7B+294m/wwFgYhIK5fXOYWnbx3FwG6d6Jjc9MPIaGAaEZE2oE+Xjjx588hmee+I9gjM7FozW2ZmB83siAGR6y03zsxWmtkaM7u9XnuBmX0Qbn/WzFrPkD0iIjEi0kNDS4EvAO8cbQEziwceAC4FBgKTzGxgePbdwH3u3gfYAdwSYT0iInKCIgoCd1/u7iuPsdhIYI27F7t7DTAVGG9mBowFng8v9wRwVST1iIjIiWuJzuLuwKZ60yXhtkxgp7vXNmhvlJndZmZFZlZUXl7ebMWKiMSaY3YWm9nrQG4js+5w95eavqTGufsUYApAYWGht9TniohEu2MGgbtfFOFnlAJ59aZ7hNu2A+lmlhDeKzjULiIiLaglDg3NA/qGzxBKBCYC093dgVnAhPByNwEttochIiIhkZ4+erWZlQBnAS+b2YxwezczewUg/Gt/MjADWA5Mc/dl4bf4CfADM1tDqM/gkUjqERGRE2ehH+Zti5mVAxtOcvUsYFsTltNWaLtjS6xuN8Tuth/Pdvdy9+yGjW0yCCJhZkXuftSL36KVtju2xOp2Q+xueyTbrXsNiYjEOAWBiEiMi8UgmBJ0AQHRdseWWN1uiN1tP+ntjrk+AhER+bRY3CMQEZF6FAQiIjEupoLgaOMiRBsze9TMysxsab22zmY208xWh58zgqyxOZhZnpnNMrOPw+NkfDfcHtXbbmbJZvahmS0Kb/d/httjYrwPM4s3swVm9o/wdNRvt5mtN7MlZrbQzIrCbSf9PY+ZIDjGuAjR5nFgXIO224E33L0v8EZ4OtrUAv/q7gOBM4Fvhf+No33bq4Gx7j4UGAaMM7MziZ3xPr5L6K4Fh8TKdo9x92H1rh046e95zAQBRxkXIeCamoW7vwNUNGgeT2jMB4jSsR/cfbO7fxR+XUnoj0N3onzbPWRPeLJd+OHEwHgfZtYDuBx4ODwdy+OcnPT3PJaC4GjjIsSKHHffHH69BcgJspjmZmb5wHDgA2Jg28OHRxYCZcBMYC0nMN5HG3Y/8GPgYHj6hMY5acMc+KeZzTez28JtJ/091+D1Mcjd3cyi9rxhM+sAvAB8z913h34khkTrtrt7HTDMzNKBF4H+wVbU/MzsCqDM3eeb2QUBl9PSRrt7qZl1AWaa2Yr6M0/0ex5LewRHGxchVmw1s64A4eeygOtpFmbWjlAIPO3ufw03x8S2A7j7TkK3dz+L8Hgf4VnR+H0/B7jSzNYTOtQ7Fvgt0b/duHtp+LmMUPCPJILveSwFQaPjIgRcU0uaTmjMB4jSsR/Cx4cfAZa7+2/qzYrqbTez7PCeAGbWHriYUP9IVI/34e7/5u493D2f0P/nN939eqJ8u80s1cw6HnoNXAIsJYLveUxdWWxmlxE6phgPPOruvwy2ouZhZn8BLiB0W9qtwM+AvwHTgJ6EbuF9nbs37FBu08xsNPAusIT/O2b8U0L9BFG77WY2hFDnYDyhH3fT3P1OM+tN6JdyZ2ABcIO7VwdXafMJHxr6obtfEe3bHd6+F8OTCcAz7v5LM8vkJL/nMRUEIiJypFg6NCQiIo1QEIiIxDgFgYhIjFMQiIjEOAWBiEiMUxCIiMQ4BYGISIz7/8dGLAQBW+rjAAAAAElFTkSuQmCC\n",
"text/plain": [
"