{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.seed:Setting seed to 1245502385\n" ] } ], "source": [ "%load_ext autoreload\n", "%autoreload 2\n", "%matplotlib inline\n", "import random\n", "random.seed(1100038344)\n", "import survivalstan\n", "import numpy as np\n", "import pandas as pd\n", "from stancache import stancache\n", "from matplotlib import pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulate survival data " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to demonstrate the use of this model, we will first simulate some survival data using `survivalstan.sim.sim_data_exp_correlated`. As the name implies, this function simulates data assuming a constant hazard throughout the follow-up time period, which is consistent with the Exponential survival function.\n", "\n", "This function includes two simulated covariates by default (`age` and `sex`). We also simulate a situation where hazard is a function of the simulated value for `sex`. \n", "\n", "We also center the `age` variable since this will make it easier to interpret estimates of the baseline hazard.\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:sim_data_exp_correlated: cache_filename set to sim_data_exp_correlated.cached.N_100.censor_time_20.rate_coefs_54462717316.rate_form_1 + sex.pkl\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:sim_data_exp_correlated: Loading result from cache\n" ] } ], "source": [ "d = stancache.cached(\n", " survivalstan.sim.sim_data_exp_correlated,\n", " N=100,\n", " censor_time=20,\n", " rate_form='1 + sex',\n", " rate_coefs=[-3, 0.5],\n", ")\n", "d['age_centered'] = d['age'] - d['age'].mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Aside: In order to make this a more reproducible example, this code is using a file-caching function `stancache.cached` to wrap a function call to `survivalstan.sim.sim_data_exp_correlated`. *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Explore simulated data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is what these data look like - this is `per-subject` or `time-to-event` form:" ] }, { "cell_type": "code", "execution_count": 3, "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", "
sexageratetrue_tteventindexage_centered
0male540.0820851.0138551.013855True0-1.12
1male390.0820854.8905974.890597True1-16.12
2female450.0497874.0934044.093404True2-10.12
3female430.0497877.0362267.036226True3-12.12
4female570.0497875.7122995.712299True41.88
\n", "
" ], "text/plain": [ " sex age rate true_t t event index age_centered\n", "0 male 54 0.082085 1.013855 1.013855 True 0 -1.12\n", "1 male 39 0.082085 4.890597 4.890597 True 1 -16.12\n", "2 female 45 0.049787 4.093404 4.093404 True 2 -10.12\n", "3 female 43 0.049787 7.036226 7.036226 True 3 -12.12\n", "4 female 57 0.049787 5.712299 5.712299 True 4 1.88" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*It's not that obvious from the field names, but in this example \"subjects\" are indexed by the field `index`.*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot these data using `lifelines`, or the rudimentary plotting functions provided by `survivalstan`." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4VFX+x/H3mfTeE0JCSEIPHUK30AREir2uws+C2Nay7i7rqquurrqWtawNFAVFBcXCIqKAIL33DgkBAqRDSALp5/fHHSCEhCRkZm4y+b6eZ55pd+795s7kk5Mz956jtNYIIYRwLhazCxBCCGF7Eu5CCOGEJNyFEMIJSbgLIYQTknAXQggnJOEuhBBOSMJdCCGckIS7EEI4IQl3IYRwQq5mbTg0NFTHxsaatXkhhGiUNmzYkKW1DqtpOdPCPTY2lvXr15u1eSGEaJSUUgdrs5x0ywghhBOScBdCCCck4S6EEE7ItD53IYQAKCkpITU1lcLCQrNLaVA8PT2Jjo7Gzc3tkl4v4S6EMFVqaip+fn7ExsailDK7nAZBa012djapqanExcVd0jpq7JZRSk1VSmUopbZX87xSSr2jlNqvlNqqlOpxSZUIIZqkwsJCQkJCJNgrUEoREhJSr/9matPn/hkw4iLPXw20sV4mAB9ccjVCiCZJgv1C9d0nNXbLaK2XKqViL7LIWGC6NubrW62UClRKRWqtj9WrsmrsWvMLGTuXs7XZ9ZS4eNtjEwD4ebrxfwNicXWR75yFEI2PLfrco4DDFe6nWh+7INyVUhMwWvfExMRc0sZy967gyoPvEHHge+4veYJDRFzSei7mzLSyPVoG0rNlsM3XL4RoeN555x0++OADevTowYwZM2y+/ueeew5fX1+efPJJm6+7Kg79QlVrPRmYDJCYmHhJM3P3vfMFSBpI+2/v5nf9PNz4CbQeatM6dx07ydVvLyPjZJFN1yuEaLjef/99Fi5cSHR0tNml2IQt+hyOAC0q3I+2PmY/rQbDfYvBPwpm3ATL3zrX3LaBMD8PADLyJNyFaAomTpxIcnIyV199NS+99BJ33303vXv3pnv37vz4448AfPbZZ1x77bVcddVVxMbG8t///pc333yT7t2707dvX3JycgCYMmUKvXr1omvXrtxwww2cOnXqgu0lJSUxYsQIevbsyeWXX87u3btt/jPZouU+B3hYKfU10AfItVd/+3mC4+DeBfDjQ7DwH3BsM4x9D9x96r3qIG93XCyKTAl3IRzq+f/tYOfRkzZdZ0Jzf/4xuuNFl/nwww+ZP38+ixcv5s0332Tw4MFMnTqVEydO0Lt3b4YONXoHtm/fzqZNmygsLKR169a8+uqrbNq0iccff5zp06fz2GOPcf3113PfffcB8PTTT/PJJ5/wyCOPnLe9CRMm8OGHH9KmTRvWrFnDgw8+yG+//WbTn7vGcFdKfQUMBEKVUqnAPwA3AK31h8A8YCSwHzgF/J9NK7wYdx+48VOI7AoLn4eM3XDDx9CsU71W62JRhPi4S7gL0QT9+uuvzJkzh9dffx0wDtU8dOgQAIMGDcLPzw8/Pz8CAgIYPXo0AJ07d2br1q2A8Qfg6aef5sSJE+Tn5zN8+PDz1p+fn8/KlSu56aabzj5WVGT7rKnN0TK31fC8Bh6yWUV1pRRc9jg06wLfT4Qpg+Gq56HPROO5SxTq60FmvoS7EI5UUwvbEbTWzJ49m3bt2p33+Jo1a/Dw8Dh732KxnL1vsVgoLS0FYPz48fzwww907dqVzz77jCVLlpy3nvLycgIDA9m8ebNdfw7nOc6v9RB4YCXED4T5k2DGjZCfccmrC/F1J7ug2GblCSEah+HDh/Puu++ird/jbdq0qU6vz8vLIzIykpKSkiqPuvH39ycuLo5vvvkGMP6YbNmypf6FV+I84Q7gGwa3z4SRr0PKcni/H+z95ZJWFerrQba03IVocp555hlKSkro0qULHTt25JlnnqnT6//5z3/Sp08fBgwYQPv27atcZsaMGXzyySd07dqVjh07nv3S1paUtuFRJnWRmJio7TpZR8YumH0vpG+Hfg/D8Jfq9PJ/zt3Jl2sOseufFzs5VwhRX7t27aJDhw5ml9EgVbVvlFIbtNaJNb3WuVruFYV3gHsXQeebYNV/IS+9Ti8P8XXndEmZtN6FEI2S84Y7gJsntLW2vAtP1Oml8aHGIZX9XvmNB2dsYP72NApLymxdoRBC2IXzD/nrGWBcF+bW6WUjOkXy3YP9mbP5KHO3HmXetjT8PF0Z2SmSsd2a0yc+BBeLDHYkhGiYmlC41/3EiB4xQfSICeLpazqwIimbHzcfYe7Wo8xcf5gIfw9Gd2nO2G5RdIryl1HthBANShMI90DjOnsftLm0MWhcXSxc2TaMK9uGcfraMhbtTufHzUeZtiqFj5cfID7Mh7FdoxjTrTlxofU/Q1YIIerL+cM9pDXE9IPF/4L2oyCwRc2vuQgvdxdGdWnOqC7NOXGqmJ+3p/Hj5iO8tWgv/1m4l67RAYztFsUtvVrg4+H8u1cI0TA59xeqABYLXPch6HL44QEoL7fZqgO93bmtdwxfT+jHykmDeWpke0rLNS/M3cm0VSk2244QouFasmQJo0aNMruMCzh/uAMExcKIVyBlGax+3y6biAzwYsIVrfjpj5cT4e9BcmaBXbYjhBC10TTCHaD7H6DdNbDoeeMEJztqGezDoewLh/kUQjRMKSkptG/fnvHjx9O2bVvuuOMOFi5cyIABA2jTpg1r165l7dq19OvXj+7du9O/f3/27NlzwXoKCgqqHC7YDE2nU1gpGP02vNEOtn0LQ+p2SnFdtAj2Zvn+TLutXwin9fMkSNtm23U26wxXv1LjYvv37+ebb75h6tSp9OrViy+//JLly5czZ84c/vWvfzF9+nSWLVuGq6srCxcu5KmnnmL27NnnreOll16qcrhgHx/HH2jRdMIdjLFnAqLhxEG7bqZliDezNxZxqrgUb/emtYuFaKzi4uLo3LkzAB07dmTIkCEopejcuTMpKSnk5uYybtw49u3bh1KKkpKSC9ZR3XDBZgyv0PSSJ6glHE+x6ybahPsC0P+V3xjcPpxhCRFc0TZMgl6ImtSihW0vNQ3n+8wzzzBo0CC+//57UlJSGDhw4AXrqG64YDM0nT73M4Ji4bh9W+7DOzbjozt7Mrh9OIt2ZTDxi410f2EB905bx6x1h8mS8WqEaHRyc3OJiooCjCn3qlLf4YJtqek1JYPjoSADNn8JXW+r14Qe1bFYFMM7NmN4x2aUlpWzNiWHBTvT+XVHOgt3ZaAUJLYMYlhCM65KiCBWTnwSosH7y1/+wrhx43jxxRe55pprqlzmmWee4bHHHqNLly6Ul5cTFxfH3LlzHVypwXmH/K1OQTbMuhMOroCO18Go/4BXkEM2rbVm57GT/LojnQU709l5zBgSoW2EL1clRDAsoRmdowKwyJg1ogmRIX+rV58hf5teuAOUl8GKt2HxS+AbAdd9BHGXO7yMwzmnWLDTCPq1KTmUlWua+XsyNCGcYQnNuKx1qAS9cHoS7tWTcL9URzYaE3rkJMNlj8HAp8DV3ZRSjhcU89vuDBbsTOf3vZmcLiljdNfmvHVLNxl9Ujg1Cffq1Sfcm16fe0VRPWDiMpj/N1j+H0haDHd8axwy6WBBPu7c0DOaG3pGU1hSxpSlybyxYC/uLhZeu7GLtOCFU9Nay8iqldS34d30jpapzN0HxrwDN0+HY5th2yyzK8LTzYVHhrTh8aFtmb0xlWd+3F7vN1qIhsrT05Ps7Gz5jFegtSY7OxtPT89LXkfTbrlXlDAWvEPtPjRBXfxxSGsKS8v4YEkSHq4uPDOqg7RuhNOJjo4mNTWVzEw5q7siT09PoqOjL/n1Eu4VhXeAzN1mV3GWUoq/DG9HYUkZU1ccwNPNwp+Ht5OAF07Fzc2NuLg4s8twOtItU1FYe8jcAw3o30OlFM+OSuD2PjG8vySJtQdyzC5JCNEISLhXFJEARSdhziOQc8Dsas5SSjG+fywA2QXF5hYjhGgUJNwr6nob9LoXts6Cd3vCd/dD5l6zqxJCiDqTcK/IzQuueQMe3QJ9H4Bdc+C93vDNeEjbbnZ1QIPqMRJCNGAS7lXxj4ThL8Fj2+DyJ2DfQvhwAHx1G6RuMLs6IYSokYT7xfiEwpBn4fFtxtmrB1fCx4Ph8+sgfadDS3G1nsS0JfWEHA8shKiRhHtteAXBwL/C49th6PNwdDNMvhKWv2WMU+MAcaE+XNc9islLk3nl590S8EKIi5JwrwsPP2MMmofWQpthsPAf8OlIY2waO1NK8cZNXbmrX0s+WprMpNnbKCuXgBdCVE3C/VL4hsEtX8B1k40zWj8YAOs+tvu3nRaL4vkxHfnj4NbMXH+YR77aSFGpY/5zEEI0LhLul0op6HoLPLgKWvSBn/4EX1wPuUfsvFnFE8Pa8fQ1HZi3LY17p62noKjUrtsUQjQ+Eu71FRAFd35vHEJ5aDVMGQzFBXbf7L2Xx/PvG7uwYn8Wj35t3lReQoiGScLdFpQyTn664xvIT4NNMxyy2ZsTW/DokLYs3JXB/ow8h2xTCNE4SLjbUuxlEN0bVv3XYUfR3NE3BncXC5+vsu+k30KIxkXC3dYG/BFOHDTObnWAUF8PrukSyeyNR8iXvnchhFWtwl0pNUIptUcptV8pNamK52OUUouVUpuUUluVUiNtX2oj0W4kBMfDinccNlbAnf1akl9UyvcbUx2yPSFEw1djuCulXID3gKuBBOA2pVRCpcWeBmZprbsDtwLv27rQRsPiAv0ehqMb4eAKh2yye4tAOkcF8NbCffxr3i7WHsihtKzcIdsWQjRMtWm59wb2a62TtdbFwNfA2ErLaMDfejsAOGq7EhuhbreDT5hxBqsDKKX413WdSWjuz6crDnDzR6vo9dJCnpi1mXnbjkl3jRBNUG1mYooCDle4nwr0qbTMc8CvSqlHAB9gqE2qa6zcvKDP/fDbi8Zoks062X2TnaMD+PyePuQVlrB0bxaLdqXz2+4Mvtt4BHcXC31bhTC0QzhDOkQQFehl93qEEOZSNY1RopS6ERihtb7Xev9OoI/W+uEKyzxhXdcbSql+wCdAJ611eaV1TQAmAMTExPQ8eNCJj/A4fRze7AgdRsH1k00pobSsnA0Hj7NodwYLdqZzIMs4/r5DpD9XdQhnaEIEnZoHYLHItH1CNBZKqQ1a68Qal6tFuPcDntNaD7fe/xuA1vrlCsvswPgDcNh6Pxnoq7XOqG69iYmJev369bX5WRqv+U/Bmg+NAcf8m5tdDUmZ+Szalc7CnRmsP5hDuYZwPw+GdIjgqoRw+rcKxdPNxewyhRAXYctwdwX2AkOAI8A64Hat9Y4Ky/wMzNRaf6aU6gAsAqL0RVbeJMI9Owne7WEMG3z5n8yu5jw5BcUs2ZPBwl3p/L4nk4LiMrzcXHh+bEduTmxhdnlCiGrUNtxr/EJVa10KPAz8AuzCOCpmh1LqBaXUGOtifwLuU0ptAb4Cxl8s2JuMkFbQ8jLY9EWDm0Ip2Med63tE8/4dPdn47FVMv7s3Ph4uLN2baXZpQggbqM0Xqmit5wHzKj32bIXbO4EBti3NSXT/A/ww0ZjoI7Zh7iIPVxeuaBtGdJA3uadLzC5HCGEDcoaqvSWMAXc/o/XewAV5u3HilIS7EM5Awt3e3H2MgN/zE5Q17OPNA73dOX6q2OwyhBA2IOHuCG2GQWEuHGnYXyAHeruRcbKIdxftY11KDsWlcparEI1VrfrcRT3FDwTlAvsWQExfs6up1oiOzViVlM0bC/bCAvB0s5DYMpi+8cH0axVC56hA3F2lPSBEY1DjoZD20iQOhaxo6ggoOQ33/252JTU6XlDMmgM5rE7OZnVyNrvTjLHivdxcSIwNom98CH3jQ+gSHYCbi4S9EI5U20MhpeXuKK0Gw+KX4PQJ8Ao0u5qLCvJxZ0SnZozo1AwwjolfeyCbVUnZrE7O4bVf9gDg7e5CYqy1ZR8fQueoAFwl7IVoECTcHcUv0rguymvw4V5ZsI87IzpFMqKT8TNk5Rex1tqyX5WUzb/nG2HvYw37fq2Mln2n5v4S9kKYRMLdUVw9jevSInPrsIFQXw9Gdo5kZGcj7DPzjLBflZzF6uQcXvl5NwC+Hq5c2S6Mf9/QBR8P+agJ4UjyG+corh7GdWmhuXXYQZifMRvUNV3Ohf3q5GxWJmXx1drDNA/w5O/XVJ4CQAhhTxLujnIm3Msaf8u9JmF+Hozu2pzRXZsDiqkrUhjbLYpOUQFmlyZEkyEdoo7iHWpc5x4xtw4HmzSiPUHe7jz1/TbKyhvW+DpCODMJd0eJ6AgWVzi6yexKHCrA241nRyewNTWX6atSzC5HiCZDumUcxc0TwhOMuVWbmNFdIpm9IZXXf9lDWbmmT1wIHSL95EgaIexIwt2RmneHnT8Yw/+qpjP7kVKKF6/txD3T1vHiT7sA47DJnrHB9IkLpndcMF2iA/BwlYlChLAVCXdHat4dNk6DE4cgqKXZ1ThUi2Bvfn38StJyC1mbksPaA9msO3D87AlR7q4WurUIPBv2PWKC5PBJIepBfnscyTfCuD6d0+TC/YxmAZ6M6dqcMV2NaQePFxSzLiWHtQdyWJuSw/tLknj3t/24WBSdmvvTOy6Y3nEh9IoNItDb3eTqhWg8JNwdycPPuC7KM7eOBiTIx51hHZsxrKMx1EF+USkbDx43wv5ADtNWHWTKsgMAtG/mR6/YYGvgBxPh72lm6UI0aBLujiThXiNfD1euaBvGFW3DACgsKWNrai5rD2Sz5kAO321M5fPVBwGIDfGmf+tQJl3dHn9PNzPLFqLBkXB3pDPhnrXP3DoaEU83l7Mt9YeB0rJydh47ydoDOaw5kMNXaw8R5uvB41e1NbtUIRoUORbNkQKiIaw9LPwHzLgZMveYXVGj4+pioUt0IPdeHs+UuxK5sm0YX609REmZTCwiREUS7o7k6gETfoerXoBDq+H9fjD3ccjPMLuyRuvOvi3JyCtiwc50s0sRokGRcHc0N08Y8Cj8cRP0ugc2Tod3esDS143JPESdDGwXTlSgF5+vOmh2KUI0KBLuZvEJgZGvwYOrIe4K+O2f8G5P2PI1lEsXQ225WBR39I1hVXI2v+1Ox6yZxYRoaCTczRbaBm77Esb/BD5h8P398NEVsPcX40xWUaNbElsQ5ufB3Z+tZ/Abv/PBkiQy8pxvaGUh6kLmUG1Iysth+2xY/CIcT4GYfjDkH9Cyn9mVNXinikuZty2NWesOszYlBxeLYnD7cG5JbMHAdmEyjo1wGrWdQ1XCvSEqLYZN0+H3f0N+OrQZBkOehWadza6sUUjKzGfW+sPM3nCErPwiwv08uLFnNDcntiA21Mfs8oSoFwl3Z1B8CtZ+BMv/A4W50OlGGPQUhLQyu7JGoaSsnN92ZzBr3WEW78mgXEPf+GBu7RXDiE7N8HSTgcpE4yPh7kxOH4cV78CaD6GsGHqMg8FPg3ew2ZU1Gmm5hXy74TCz1qdyKOcU/p6ujO0WxSODWxMuwxiIRkTC3RnlpcPSf8P6T8HT3+iq6TEOLNICra3ycs3qA9nMXHeYn7enEe7nwRf39JHuGtFoSLg7s/SdMO/PcHA5RHaDa96A6Brfa1HJ1tQTjJu6FheLhel39yahub/ZJQlRo9qGuxxC0BhFJMD4uXDDJ5CXBh8PgR8fhoIssytrVLpEB/LNxP64uShumbyKdSk5ZpckhM1IuDdWSkHnG+GR9dD/EdjylXES1NopUF5mdnWNRutwX759oD9hvh7c+ckaFu+WoSCEc5Bwb+w8/GDYizBxhXGo5Lwn4Ze/m11VoxIV6MU3E/vROtyX+6av5+u1hygrlxPIROMm4e4swtvDuP9B4j2w5gNIWWF2RY1KiK8HX93Xl8TYICZ9t40r/r2YD5YkkVNQbHZpQlwS+ULV2RQXwAf9AQUPrAR3b7MralRKy8pZuCudaSsPsio5G3dXC2O6Nueufi3pEh1odnlCyNEyTdqBZTBtFPR9EEa8bHY1jdbe9Dymr0rhu41HOFVcRrcWgYzr35KRnSPxcJXDT4U5JNybup+ehHUfw//9LGPT1NPJwhK+25DK9FUHSc4qIMTHndt6x3B7nxiaB3qZXZ5oYiTcm7qifGMyEDcvmLgcXN3NrqjRKy/XrEjKYvqqgyzalY5Siqs6RDBxYCu6tZAuG+EYcpx7U+fha5zclLUHVr5tdjVOwWJRXN4mjCl3JfL7nwdx3+XxrDmQzS0frSIzr8js8oQ4T63CXSk1Qim1Rym1Xyk1qZplblZK7VRK7VBKfWnbMsUlaTsMEsYaszxlJ5ldjVNpEezNpKvb892DAyguK+fTFQfMLkmI89QY7kopF+A94GogAbhNKZVQaZk2wN+AAVrrjsBjdqhVXIoRr4LFDX76k0z+YQdxoT6M7BTJ56sOcrKwxOxyhDirNi333sB+rXWy1roY+BoYW2mZ+4D3tNbHAbTWcppfQ+EfaQwwlrzYmMpPAt7mHhjYiryiUmasPmR2KUKcVZtwjwIOV7ifan2sorZAW6XUCqXUaqXUiKpWpJSaoJRar5Ran5mZeWkVi7rrda8xeuSyN2D+3yTgbaxTVACXtwnlk+UHyD0lrXfRMLjacD1tgIFANLBUKdVZa32i4kJa68nAZDCOlrHRtkVNLBYY/Ta4+8Dq96GkAEa9JUMF29BDg1pz6+TVdP/nr3SOCqBvqxD6xYfQKzYYHw9b/ZoJUXu1+dQdAVpUuB9tfayiVGCN1roEOKCU2osR9utsUqWoP6Vg+L+MgF/6mjHL03Ufgoub2ZU5hb7xIXz3YH+W7M5gVXI2U5cf4KPfk3G1KLpEB9CvVQj94kPp2TIIL3f5oyrsrzbhvg5oo5SKwwj1W4HbKy3zA3Ab8KlSKhSjmybZloUKG1DKmMHJzRsWPQ+lhXDTNHCRlqUt9IgJokdMEGBM2L3h4HFWJWWzKjmbD39P5r3FSbi5KLq1CKRffAh9W4XQIyZIpvsTdlGrk5iUUiOBtwAXYKrW+iWl1AvAeq31HKWUAt4ARgBlwEta668vtk45iclkqz+A+ZNkiAIHyS8qZV1KDqutYb/9SC7lGtxdLfSICaRffCj9WoXQrUUg7q5y+omonpyhKmr28yRjBMnRb0PP8WZX06Tkni5h3YEcViVnsyopm11pJ9EaPN0sJLYMpl+rEPrGB5MQGSDdOOI8Eu6iZmWl8NUtkLwE7vwB4i43u6Im68SpYlYn57DaGvZ70vMAoyetZbA3bSP8aN/Mj7bNjOvYEB9cXaSF3xRJuIvaKcyFj6+Cggy4dxGEtDK7IgFk5xexLuU4u9NOsjc9j91peaRkFXBmDhF3Fwutwn2NwLcGf7tmfkQGeGL0kgpnJeEuai8nGaYMNlryve+Dfg+BT6jZVYlKCkvK2J+Rz560vLOBvyctj7SThWeX8fN0pV3EuRZ+uwgj9AO9ZeA4ZyHhLuomcy8seRl2fA+unpB4tzE3q3+k2ZWJGuSeKmFPep5xSTvJnjQj+PMKS88uE+HvUaGF70+7CD/aRPjKkTqNkIS7uDRZ+2DZm7B1pnGSU/c7YcCjENTS7MpEHWitSTtZyB5r635PmhH++zLyKS4tB8CioGWIz9nWfbtmfiRE+hMb6mNy9eJiJNxF/RxPgeVvwaYvAA1dboXLn5A++UautKyclOxTFbp1TrI3PZ+U7IKzo1L8dUR7Hhgo73NDJeEubCP3CKx8BzZ8BmXF0HsCXP2q2VUJGztdbPTnf/D7fuZtS2PynT0Z1rGZ2WWJKshkHcI2AqKMMH9sG7QbCWs+hKI8s6sSNubl7kLn6ADevLkbXaMDeGzmZnYePWl2WaIeJNxF7fiGQ4cxxu28dHNrEXbj6ebClLsS8fd0477p62WGqUZMwl3Unp/13/S8o+bWIewq3N+TKXclkl1QxMQvNlBUWmZ2SeISSLiL2vNvblznpZlbh7C7ztEBvHFTNzYcPM7Vby9jytJksvOlFd+YSLiL2vOLBBSkbTO7EuEA13SJ5IM7ehDk7c5L83bR9+VFPDRjI0v3ZlJeLtMxNHRytIyom5l/gANL4fEd4OFndjXCQfam5zFz3WG+25jK8VMlRAV6cUuvFtyUGE1kgJfZ5TUpciiksI/UDfDxYBj2onEGq2hSikrL+HVHOjPXHWb5/iwsCga2C+eWXi0Y3D4cNxnMzO4k3IX9TBttnMn66BZw9TC7GmGSQ9mnmLX+MN9sOEz6ySLC/Dy4sWc0tyS2kLNc7UiOcxf2c9njkHcM1n9qdiXCRDEh3jw5vB0r/jqYT8Yl0q1FIJOXJjPw9SX8uLnyTJzC0STcRd3FDzIuv/wNtsw0uxphMlcXC0M6RDDlrkRWThpMbIg3325INbusJk/CXdSdUnDrDIi9DL6/3zr+jBAQ4e/J0A4RrEnO4VRxac0vEHYj4S4ujbsP3D4LWg2CHx+SLhpx1sB24RSXlbMqKdvsUpo0CXdx6dy84NavoM0wmPsYrJ1idkWiAegVF4SXmwu/7800u5QmTcJd1I+bJ9zyBbS7BuY9CV/fYQwXLJosD1cX+rcK4dcd6eSeKjG7nCZLwl3Un6sH3DwNhjwLSb/Be31g8b+g+JTZlQmT3H9lK7ILirhv+noKS2RsGjNIuAvbcHGDy/8ED6+H9tfA768aIb9zDph0LoUwT++4YN68uRtrU3J47OvNlMlwBQ4n4S5sKyAKbpwK438yhieYdSd8fi1k7jG7MuFgo7s259lRCczfkcZzc3Zg1gmTTZWr2QUIJxV7Gdy/FNZPhcUvwgf9IWEsNO8OzTpDsy7gHWx2lcLO7r4sjvSThXy0NJkIfw8eHtzG7JKaDAl3YT8urtBnAnS63uiD3/MzbJ997nn/KGvQV7gExoJF/qF0Jn8d0Z6MvCJe/3Uvob4e3No7xuySmgQJd2F/PqEw6k3jUpBlDBlc8bJvAWjrl27uftCs0/mBH9bBOCpHNEoWi+LVG7qQU1DM377fho+HK6O7Nje7LKcnA4cJ85Wchoxd5wd++nYozjeeVy4Q1u78wI/oDD4h5tYt6uR0cRnjpq5l46HwHrjiAAASzklEQVTjTL6rJ4PbR5hdUqMko0KKxq28HI4fuLCVX3GKv8rdOq0GyxjzDdzJwhJun7Kafen5TLu7N33j5Q90XUm4C+dUVbdO1h7Q5dBqCNz5ndkVihrkFBRz80erOHbiNF/e15euLQLNLqlRkXAXTUfJaVj8Eqx8Fx7bDoEtzK5I1CAtt5CbPlpJ7qkSBrYLJy7Uh/gwH+JDfYkN9cbP083sEhus2oa7fKEqGj83L0i8xwj3bbOMk6lEg9YswJMZ9/Tlhbk72XT4OP/bevS8c93C/DyItwZ+XKgPcaG+xIf50CLIG3dXOZqqNqTlLpzHJ8Pg9Al4aI0xLLFoNApLyjiUc4rkzAKSs/I5kFnAgawCkrMKyCkoPruci0XRIsiL+DBfa+j7WP8I+BLh74FqAu+7tNxF09PlFvjpCTi6CaJ6mF2NqANPNxfaRvjRNuLCL8RPnCo2gt4a+AeyCkjKzGdlUhaFJeVnl/N2dyE25Ez3jg9x1m6euDAf/JtgN4+03IXzOH0c/tMZWg+Gm6ebXY2ws/JyTdrJQmvw55NsDf7kzAJSj5+i4nA2ob7u1la+EfZxoT60CvOhRbA3Hq4u5v0Ql0Ba7qLp8QqCPvfDstchbbtxMpRwWhaLonmgF80DvRjQOvS854pKyzicc4qkM619a3fPot3pZK0/181jURAd5E2bcF/6xAfTv1UoCZH+WCyNv3tHWu7CuZzKgbe7QvyVxjjzQlSSe7rE2r1j9O0nZRWw69hJkjMLAAj0dqNvXAgDWofQv3Uo8aE+DaovX1ruomnyDoa+DxhDDqcsh5j+MlaNOE+AlxvdWgTSrdLx9Wm5haxKzmLF/mxW7s9i/o40ACL8PRjQKpR+rYywjwr0MqPsOqtVy10pNQJ4G3ABPtZav1LNcjcA3wK9tNYXbZZLy13YzekTRuu98AS4+UB4B4joCBGdICIBwhNkREpxUVprDmafYmVSNiuSsliVlH32qJ3YEG/6tQplQOsQ+sWHEOLr4dDabHYSk1LKBdgLXAWkAuuA27TWOyst5wf8BLgDD0u4C1OdOATJSyB9h/Wy3fjC9Qz/KCPkIzqeu4S0AVd300oWDVd5uWZPeh4rk4xW/ZoDOeQXlQLQvpkfA1qH0r9VCL3jgu1+ApYtw70f8JzWerj1/t8AtNYvV1ruLWAB8GfgSQl30aBoDXlpRtBn7DgX+pl7oNw6z6fFDULbnh/4ER3BL1KOmxfnKS0rZ+uRXFYlZbMyKYt1KccpLi3HxaLoEh3A40PbckXbMLts25Z97lHA4Qr3U4E+lTbWA2ihtf5JKfXnOlUqhCMoBf6RxqXN0HOPl5VA1r7zQ//gSuNM1zO8giD8TNgnGN074R3A3cfxP4doEFxdLPSICaJHTBAPDWpNYUkZGw8dZ+X+bOZsOcojX21i4RNXEubn2C6b82qs7wqUUhbgTWB8LZadAEwAiImRAftFA+DiZg3sBOCmc4+fPg7pOyFjp9Glk74DNs84NwwxCoJioe0IGPw0ePiaULxoKDzdXOjfKpT+rUIZ260517yznBfm7uTd27qbVlNtwv0IUHEkpmjrY2f4AZ2AJdbDhZoBc5RSYyp3zWitJwOTweiWqUfdQtiXVxDEDjAuZ5SXw4mD1lb+Tji2BdZ8CHt/hus+gpi+5tUrGow2EX48NKg1/1m4l2u7NWdIB3PGra9Nn7srxheqQzBCfR1wu9Z6RzXLL0H63EVTcXAlfD/R+AJ3wKMw6ClwNe9fcdEwFJeWM+rdZeQVlrLgiSvx9bDdUee17XOv8QBgrXUp8DDwC7ALmKW13qGUekEpNab+pQrRiLXsDw+sgB53wYq3YPIgY4x50aS5u1p4+foupJ0s5OV5u0ypQc5QFcJW9v4CPz5s9NcPegr6/9GYJFw0WS/9tJMpyw7w3OgExg+Is8k6bdZyF0LUUtvh8OBqaD8SFj0PHw+Go5vNrkqYaNLVHbgqIYLn5+7kp63HHLptCXchbMknBG6aBjd+ahxXP2UQzH8KivJrfq1wOi4Wxbu3dadHTBCPz9zMqqRsh21bwl0IW1MKOl0PD62FHuNg9Xvwfl/YM9/syoQJPN1c+GRcIjEh3kyYvp5dx046ZLsS7kLYi1cgjH4L7v7FOOHpq1tg1l1Gi140KYHe7ky7uzc+Hq6M/3QtqcdP2X2bEu5C2FtMX7h/mXGy05758N9esPlLs6sSDhYV6MW0u3tzuriMZfuy7L49OVpGCEfKToI5j8DBFdDlVrjmdfC4cGo54byy84vqNZKkHC0jREMU0grG/Q8GPmWMX/PRlcaZrqLJcNQQwRLuQjiaxQUG/tUI+ZJT8PFQWPORMXKlEDYi4S6EWWIvg4krIH4Q/PwX+PoOyM8wuyrhJCTchTCTTwjcPhOGvwz7foU3Oxghv+dnKCs1uzrRiMm50UKYTSno9yC0GQYbP4MtM2H3XPAJh663Qvc/QFg7s6sUjYwcLSNEQ1NWAvsXwqYvYO98KC+F6F7Q7Q7j5CjPALMrFCay2TR79iLhLkQt5GfC1plG0GfuAlcvSBhjtOZbXgYW6VltaiTchXAmWsPRjbBpBmz7FopyISAGontCcDwEt7Jex4NvuMz56sRsOYeqEMJsSkFUT+My/CXY/RNsn20cI79zDuiyc8u6+RghHxJ/LvDPBn8zae03ERLuQjQ2bl7Q+UbjAkYffe5hyE6GnAqX9J2wex6Ul5x7rauXNejjLgx+/ygJfici4S5EY+fidi6gKysvM4L/bOgfMK6z98O+BVBWVGE9HpVCP+5cl09AtHHylWg0JNyFcGYWFwiKNS6tBp//XHk5nDxyfmv/zCVpMZSerrAeN2MdFVv64e3lS90GTMJdiKbKYoHAFsYl/srznysvh/y0C0M/OxlSlkNJgbFcZFcY+tyFfziE6STchRAXsljAv7lxib3s/Oe0NoZJSFoEi1+Gz6+DuCuNkI/qYUa1ogry/5QQom6UAr8I6HY7PLIeRrwC6duNKQVnjYOs/WZXKJBwF0LUh6sH9H0AHt0CV04yvqR9rzf87zGZccpkchKTEMJ28jNg6Wuw/lOwuELsAONIG/9o4zogGgKijMMuXR0zrrmzkZOYhBCO5xsOI1+Dvg/C8jfh2FbjRKuCzCqWjbAGfxQEtDgX/AHRxn3vUDkSpx4k3IUQthccB2PePXe/5DScPAq5qRUuh41DMTN3GwOllVSaNNrFw/hC90zYnw3+Cv8JePg69udqRCTchRD25+ZlTDEY0qrq57WG08fPD/+TFW4fWAp5R0GXn/86z8Dqgz8gGvwiwaVpxlzT/KmFEA2LUuAdbFwiu1S9TFkp5B2zBv8Ro+Wfmwq5R4zrQ6uh8ESl9VqMgD8b/NYuoOA4YwYsJw5+5/3JhBDOxcX13ElX1SnKrzr4cw/D0U2w639QVmwsG9beGISt9VDH1O9gEu5CCOfh4WvMWlXdzFXl5XAqCw6uhIXPwRc3QOurjJB3stmu5KtoIUTTYbEYR/R0vBYeWgPDXoTDa+H9fjDvL3Aqx+wKbUbCXQjRNLl6QP9H4I8boed4WDcF3ukOqz8whlFu5CTchRBNm08ojHoTJq6A5t1h/iTjLNvFLxvH6Zt0omd9yRmqQghxhtaw71dY/hYcWgVoCIyB9qOMS0xf08e1lzlUhRCiPvIzYe/PsGsuJC8xJjbxDoF2VxtBHz/QOH7fwSTchRDCVoryjLNod/8Ee3+BopPGXLWthxhB33YYeAU5pBQZW0YIIWzFww86XmdcSoshZZkR9Lt/gl1zrIOkXWYEfddbjeVNJi13IYS4VOXlcHSjcXLU7rnG3LThHeH2mRc/2aoeattyl6NlhBDiUlksEJ0IVz0Pj2yAP8w2zob9eAgc2WhuaaZuXQghnEnroXDPr8aIlp+ONFr0JqlVuCulRiil9iil9iulJlXx/BNKqZ1Kqa1KqUVKqZa2L1UIIRqB8A5w3yKI6Agz74SV75pyrHyN4a6UcgHeA64GEoDblFIJlRbbBCRqrbsA3wL/tnWhQgjRaPiGw/i5kDAGfn0alr3u8BJq03LvDezXWidrrYuBr4GxFRfQWi/WWp8ZaX81EG3bMoUQopFx84IbP4O4K2DL1w7ffG3CPQo4XOF+qvWx6twD/FyfooQQwilYLMbJTtn7jclIHLlpW65MKfUHIBF4rZrnJyil1iul1mdmVjGnohBCOJso61GLRzY4dLO1CfcjQMUDNqOtj51HKTUU+DswRmtdVNWKtNaTtdaJWuvEsLCwS6lXCCEal+bdAQWpjj2vpzbhvg5oo5SKU0q5A7cCcyouoJTqDnyEEewZti9TCCEaKU9/aNYZtn/n0KGEawx3rXUp8DDwC7ALmKW13qGUekEpNca62GuAL/CNUmqzUmpONasTQoimZ9BTkLUHVr/vsE3K8ANCCOEIX91mjC758Dpjsu5LJMMPCCFEQzLiFeNkpvl/c8jmJNyFEMIRglrCFU8ao0juW2j3zcmQv0II4Sj9HzEm5Hb1sPumJNyFEMJRXD3gjlkO2ZR0ywghhBOScBdCCCck4S6EEE5Iwl0IIZyQhLsQQjghCXchhHBCEu5CCOGEJNyFEMIJmTZwmFIqEzh4iS8PBbJsWI6tSF11I3XVXUOtTeqqm/rU1VJrXeOEGKaFe30opdbXZlQ0R5O66kbqqruGWpvUVTeOqEu6ZYQQwglJuAshhBNqrOE+2ewCqiF11Y3UVXcNtTapq27sXlej7HMXQghxcY215S6EEOIiGnS4K6VGKKX2KKX2K6UmVfG8h1JqpvX5NUqpWAfU1EIptVgptVMptUMp9WgVywxUSuVaJwvfrJR61t51WbebopTaZt3mBRPUKsM71v21VSnVwwE1tauwHzYrpU4qpR6rtIzD9pdSaqpSKkMptb3CY8FKqQVKqX3W66BqXjvOusw+pdQ4O9f0mlJqt/V9+l4pFVjNay/6ntuptueUUkcqvF8jq3ntRX9/7VDXzAo1pSilNlfzWrvss+qywbTPl9a6QV4AFyAJiAfcgS1AQqVlHgQ+tN6+FZjpgLoigR7W237A3irqGgjMNWGfpQChF3l+JPAzoIC+wBoT3tM0jON0TdlfwBVAD2B7hcf+DUyy3p4EvFrF64KBZOt1kPV2kB1rGga4Wm+/WlVNtXnP7VTbc8CTtXivL/r7a+u6Kj3/BvCsI/dZddlg1uerIbfcewP7tdbJWuti4GtgbKVlxgLTrLe/BYYopZQ9i9JaH9Nab7TezgN2AVH23KYNjQWma8NqIFApFenA7Q8BkrTWl3ryWr1prZcCOZUervg5mgZcW8VLhwMLtNY5WuvjwAJghL1q0lr/qrUutd5dDUTbYlt1Vc3+qo3a/P7apS5rBtwMfGWr7dWypuqywZTPV0MO9yjgcIX7qVwYomeXsf4i5AIhDqkOsHYDdQfWVPF0P6XUFqXUz0qpjg4qSQO/KqU2KKUmVPF8bfapPd1K9b9wZuyvMyK01sest9OAiCqWMXPf3Y3xH1dVanrP7eVha5fR1Gq6GczcX5cD6VrrfdU8b/d9VikbTPl8NeRwb9CUUr7AbOAxrfXJSk9vxOh66Aq8C/zgoLIu01r3AK4GHlJKXeGg7dZIKeUOjAG+qeJps/bXBbTxP3KDOYRMKfV3oBSYUc0iZrznHwCtgG7AMYwukIbkNi7earfrPrtYNjjy89WQw/0I0KLC/WjrY1Uuo5RyBQKAbHsXppRyw3jzZmitv6v8vNb6pNY633p7HuCmlAq1d11a6yPW6wzge4x/jSuqzT61l6uBjVrr9MpPmLW/Kkg/0z1lvc6oYhmH7zul1HhgFHCHNRQuUIv33Oa01ula6zKtdTkwpZptmvJZs+bA9cDM6pax5z6rJhtM+Xw15HBfB7RRSsVZW323AnMqLTMHOPOt8o3Ab9X9EtiKtT/vE2CX1vrNapZpdqbvXynVG2M/2/WPjlLKRynld+Y2xhdy2ystNge4Sxn6ArkV/l20t2pbU2bsr0oqfo7GAT9WscwvwDClVJC1G2KY9TG7UEqNAP4CjNFan6pmmdq85/aoreL3NNdVs83a/P7aw1Bgt9Y6taon7bnPLpIN5ny+bP2NsS0vGEd37MX41v3v1sdewPjAA3hi/Ju/H1gLxDugpssw/q3aCmy2XkYCE4GJ1mUeBnZgHCGwGujvgLrirdvbYt32mf1VsS4FvGfdn9uARAe9jz4YYR1Q4TFT9hfGH5hjQAlGv+Y9GN/TLAL2AQuBYOuyicDHFV57t/Wzth/4PzvXtB+jD/bMZ+zMUWHNgXkXe88dsL8+t35+tmIEV2Tl2qz3L/j9tWdd1sc/O/O5qrCsQ/bZRbLBlM+XnKEqhBBOqCF3ywghhLhEEu5CCOGEJNyFEMIJSbgLIYQTknAXQggnJOEuhBBOSMJdCCGckIS7EEI4of8HOxEhzr3d4vsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "survivalstan.utils.plot_observed_survival(df=d[d['sex']=='female'], event_col='event', time_col='t', label='female')\n", "survivalstan.utils.plot_observed_survival(df=d[d['sex']=='male'], event_col='event', time_col='t', label='male')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## model1: original spec" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "model_code = '''\n", "functions {\n", " // Defines the log survival\n", " vector log_S (vector t, real shape, vector rate) {\n", " vector[num_elements(t)] log_S;\n", " for (i in 1:num_elements(t)) {\n", " log_S[i] = gamma_lccdf(t[i]|shape,rate[i]); \n", " }\n", " return log_S;\n", " }\n", " \n", " // Defines the log hazard\n", " vector log_h (vector t, real shape, vector rate) {\n", " vector[num_elements(t)] log_h;\n", " vector[num_elements(t)] ls;\n", " ls = log_S(t,shape,rate);\n", " for (i in 1:num_elements(t)) {\n", " log_h[i] = gamma_lpdf(t[i]|shape,rate[i]) - ls[i];\n", " }\n", " return log_h;\n", " }\n", " \n", " // Defines the sampling distribution\n", " real surv_gamma_lpdf (vector t, vector d, real shape, vector rate) {\n", " vector[num_elements(t)] log_lik;\n", " real prob;\n", " log_lik = d .* log_h(t,shape,rate) + log_S(t,shape,rate);\n", " prob = sum(log_lik);\n", " return prob;\n", " }\n", "}\n", "\n", "data {\n", " int N; // number of observations\n", " vector[N] y; // observed times\n", " vector[N] event; // censoring indicator (1=observed, 0=censored)\n", " int M; // number of covariates\n", " matrix[N, M] x; // matrix of covariates (with n rows and H columns)\n", "}\n", "\n", "parameters {\n", " vector[M] beta; // Coefficients in the linear predictor (including intercept)\n", " real alpha; // shape parameter\n", "}\n", "\n", "transformed parameters {\n", " vector[N] linpred;\n", " vector[N] mu;\n", " linpred = x*beta;\n", " for (i in 1:N) {\n", " mu[i] = exp(linpred[i]);\n", " }\n", "}\n", "\n", "model {\n", " alpha ~ gamma(0.01,0.01);\n", " beta ~ normal(0,5);\n", " y ~ surv_gamma(event, alpha, mu);\n", "}\n", "'''" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we are ready to fit our model using `survivalstan.fit_stan_survival_model`. \n", "\n", "We pass a few parameters to the fit function, many of which are required. See ?survivalstan.fit_stan_survival_model for details. \n", "\n", "Similar to what we did above, we are asking `survivalstan` to cache this model fit object. See [stancache](http://github.com/jburos/stancache) for more details on how this works. Also, if you didn't want to use the cache, you could omit the parameter `FIT_FUN` and `survivalstan` would use the standard pystan functionality.\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:Step 1: Get compiled model code, possibly from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:StanModel: cache_filename set to anon_model.cython_0_29_2.model_code_14429915565770599621.pystan_2_18_1_0.stanmodel.pkl\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:StanModel: Loading result from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:Step 2: Get posterior draws from model, possibly from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:sampling: cache_filename set to anon_model.cython_0_29_2.model_code_14429915565770599621.pystan_2_18_1_0.stanfit.chains_4.data_36753546383.iter_5000.seed_9001.pkl\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:sampling: Loading result from cache\n" ] } ], "source": [ "testfit = survivalstan.fit_stan_survival_model(\n", " model_cohort = 'model 1',\n", " model_code = model_code,\n", " df = d,\n", " time_col = 't',\n", " event_col = 'event',\n", " formula = '~ age_centered + sex',\n", " iter = 5000,\n", " chains = 4,\n", " seed = 9001,\n", " FIT_FUN = stancache.cached_stan_fit,\n", " drop_intercept = False,\n", " )\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# 0:01:33.270480 elapsed" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " mean se_mean sd 2.5% 50% 97.5% Rhat\n", "lp__ -272.845825 0.022499 1.435829 -276.502302 -272.511241 -271.047870 1.000226\n", "alpha 1.023739 0.002147 0.152227 0.755217 1.014877 1.349258 1.001638\n", "beta[1] -3.029740 0.004090 0.274891 -3.614153 -3.011143 -2.532143 1.001077\n", "beta[2] 0.618876 0.003042 0.239763 0.156590 0.614498 1.108694 1.000019\n", "beta[3] -0.003055 0.000162 0.013851 -0.030628 -0.002910 0.023814 1.000772\n" ] } ], "source": [ "survivalstan.utils.print_stan_summary([testfit], pars=['lp__', 'alpha', 'beta'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## model2: alternate version of surv_gamma_lpdf" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "model_code2 = '''\n", "functions {\n", " // Defines the log survival\n", " real surv_gamma_lpdf (vector t, vector d, real shape, vector rate) {\n", " vector[num_elements(t)] log_lik;\n", " real prob;\n", " for (i in 1:num_elements(t)) {\n", " log_lik[i] = d[i] * (gamma_lpdf(t[i]|shape,rate[i]) - gamma_lccdf(t[i]|shape,rate[i]))\n", " + gamma_lccdf(t[i]|shape,rate[i]);\n", " }\n", " prob = sum(log_lik);\n", " return prob;\n", " }\n", "}\n", "data {\n", " int N; // number of observations\n", " vector[N] y; // observed times\n", " vector[N] event; // censoring indicator (1=observed, 0=censored)\n", " int M; // number of covariates\n", " matrix[N, M] x; // matrix of covariates (with n rows and H columns)\n", "}\n", "parameters {\n", " vector[M] beta; // Coefficients in the linear predictor (including intercept)\n", " real alpha; // shape parameter\n", "}\n", "transformed parameters {\n", " vector[N] mu;\n", " {\n", " vector[N] linpred;\n", " linpred = x*beta;\n", " mu = exp(linpred);\n", " }\n", "}\n", "\n", "model {\n", " alpha ~ gamma(0.01,0.01);\n", " beta ~ normal(0,5);\n", " y ~ surv_gamma(event, alpha, mu);\n", "}\n", "'''" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:Step 1: Get compiled model code, possibly from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:StanModel: cache_filename set to anon_model.cython_0_29_2.model_code_9177012762674257483.pystan_2_18_1_0.stanmodel.pkl\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:StanModel: Loading result from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:Step 2: Get posterior draws from model, possibly from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:sampling: cache_filename set to anon_model.cython_0_29_2.model_code_9177012762674257483.pystan_2_18_1_0.stanfit.chains_4.data_36753546383.iter_5000.seed_9001.pkl\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:sampling: Loading result from cache\n" ] } ], "source": [ "testfit2 = survivalstan.fit_stan_survival_model(\n", " model_cohort = 'model 2',\n", " model_code = model_code2,\n", " df = d,\n", " time_col = 't',\n", " event_col = 'event',\n", " formula = '~ age_centered + sex',\n", " iter = 5000,\n", " chains = 4,\n", " seed = 9001,\n", " FIT_FUN = stancache.cached_stan_fit,\n", " drop_intercept = False,\n", " )\n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [], "source": [ "# 0:01:20.742172 elapsed" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " mean se_mean sd 2.5% 50% 97.5% Rhat\n", "lp__ -272.859001 0.023962 1.470993 -276.580657 -272.522315 -271.028212 1.000202\n", "alpha 1.022005 0.002037 0.150233 0.749135 1.014273 1.341307 1.000132\n", "beta[1] -3.032481 0.003923 0.273614 -3.612430 -3.020702 -2.543998 1.000088\n", "beta[2] 0.614401 0.003031 0.241061 0.152624 0.610390 1.103668 0.999922\n", "beta[3] -0.003198 0.000161 0.013956 -0.031153 -0.003048 0.023837 0.999972\n" ] } ], "source": [ "survivalstan.utils.print_stan_summary([testfit2], pars=['lp__', 'alpha', 'beta'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## model3: use `log_mix` inside surv_gamma_lpdf" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "model_code3 = '''\n", "functions {\n", " // Defines the log survival\n", " real surv_gamma_lpdf (vector t, vector d, real shape, vector rate) {\n", " vector[num_elements(t)] log_lik;\n", " real prob;\n", " for (i in 1:num_elements(t)) {\n", " log_lik[i] = log_mix(d[i], gamma_lpdf(t[i]|shape,rate[i]), gamma_lccdf(t[i]|shape,rate[i]));\n", " }\n", " prob = sum(log_lik);\n", " return prob;\n", " }\n", "}\n", "data {\n", " int N; // number of observations\n", " vector[N] y; // observed times\n", " vector[N] event; // censoring indicator (1=observed, 0=censored)\n", " int M; // number of covariates\n", " matrix[N, M] x; // matrix of covariates (with n rows and H columns)\n", "}\n", "parameters {\n", " vector[M] beta; // Coefficients in the linear predictor (including intercept)\n", " real alpha; // shape parameter\n", "}\n", "transformed parameters {\n", " vector[N] linpred;\n", " vector[N] mu;\n", " linpred = x*beta;\n", " mu = exp(linpred);\n", "}\n", "\n", "model {\n", " alpha ~ gamma(0.01,0.01);\n", " beta ~ normal(0,5);\n", " y ~ surv_gamma(event, alpha, mu);\n", "}\n", "'''" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:Step 1: Get compiled model code, possibly from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:StanModel: cache_filename set to anon_model.cython_0_29_2.model_code_1293841621968646714.pystan_2_18_1_0.stanmodel.pkl\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:StanModel: Loading result from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:Step 2: Get posterior draws from model, possibly from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:sampling: cache_filename set to anon_model.cython_0_29_2.model_code_1293841621968646714.pystan_2_18_1_0.stanfit.chains_4.data_36753546383.iter_5000.seed_9001.pkl\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:sampling: Loading result from cache\n" ] } ], "source": [ "testfit3 = survivalstan.fit_stan_survival_model(\n", " model_cohort = 'model 3',\n", " model_code = model_code3,\n", " df = d,\n", " time_col = 't',\n", " event_col = 'event',\n", " formula = '~ age_centered + sex',\n", " iter = 5000,\n", " chains = 4,\n", " seed = 9001,\n", " FIT_FUN = stancache.cached_stan_fit,\n", " drop_intercept = False,\n", " )\n" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# 0:00:42.036498 elapsed" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " mean se_mean sd 2.5% 50% 97.5% Rhat\n", "lp__ -272.836212 0.023769 1.461660 -276.474541 -272.493666 -271.040873 1.001796\n", "alpha 1.023676 0.001992 0.149010 0.746153 1.017081 1.331976 0.999766\n", "beta[1] -3.030396 0.004034 0.267378 -3.606247 -3.013798 -2.553470 1.000366\n", "beta[2] 0.619250 0.003194 0.242402 0.156037 0.617779 1.106182 1.000445\n", "beta[3] -0.003286 0.000155 0.013740 -0.030794 -0.003000 0.023025 1.000123\n" ] } ], "source": [ "survivalstan.utils.print_stan_summary([testfit3], pars=['lp__', 'alpha', 'beta'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## model4: vectorize surv_gamma_lpdf" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "model_code4 = '''\n", "functions {\n", " int count_value(vector a, real val) {\n", " int s;\n", " s = 0;\n", " for (i in 1:num_elements(a)) \n", " if (a[i] == val) \n", " s = s + 1;\n", " return s;\n", " }\n", "\n", " // Defines the log survival\n", " real surv_gamma_lpdf (vector t, vector d, real shape, vector rate, int num_cens, int num_obs) {\n", " vector[2] log_lik;\n", " int idx_obs[num_obs];\n", " int idx_cens[num_cens];\n", " real prob;\n", " int i_cens;\n", " int i_obs;\n", " i_cens = 1;\n", " i_obs = 1;\n", " for (i in 1:num_elements(t)) {\n", " if (d[i] == 1) {\n", " idx_obs[i_obs] = i;\n", " i_obs = i_obs+1;\n", " }\n", " else {\n", " idx_cens[i_cens] = i;\n", " i_cens = i_cens+1;\n", " }\n", " }\n", " print(idx_obs);\n", " log_lik[1] = gamma_lpdf(t[idx_obs] | shape, rate[idx_obs]);\n", " log_lik[2] = gamma_lccdf(t[idx_cens] | shape, rate[idx_cens]);\n", " prob = sum(log_lik);\n", " return prob;\n", " }\n", "}\n", "data {\n", " int N; // number of observations\n", " vector[N] y; // observed times\n", " vector[N] event; // censoring indicator (1=observed, 0=censored)\n", " int M; // number of covariates\n", " matrix[N, M] x; // matrix of covariates (with n rows and H columns)\n", "}\n", "transformed data {\n", " int num_cens;\n", " int num_obs;\n", " num_obs = count_value(event, 1);\n", " num_cens = N - num_obs;\n", "}\n", "parameters {\n", " vector[M] beta; // Coefficients in the linear predictor (including intercept)\n", " real alpha; // shape parameter\n", "}\n", "transformed parameters {\n", " vector[N] linpred;\n", " vector[N] mu;\n", " linpred = x*beta;\n", " mu = exp(linpred);\n", "}\n", "model {\n", " alpha ~ gamma(0.01,0.01);\n", " beta ~ normal(0,5);\n", " y ~ surv_gamma(event, alpha, mu, num_cens, num_obs);\n", "}\n", "'''" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:Step 1: Get compiled model code, possibly from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:StanModel: cache_filename set to anon_model.cython_0_29_2.model_code_16881928540873162731.pystan_2_18_1_0.stanmodel.pkl\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:StanModel: Loading result from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:Step 2: Get posterior draws from model, possibly from cache\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:sampling: cache_filename set to anon_model.cython_0_29_2.model_code_16881928540873162731.pystan_2_18_1_0.stanfit.chains_4.data_36753546383.iter_5000.seed_9001.pkl\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:stancache.stancache:sampling: Loading result from cache\n" ] } ], "source": [ "testfit4 = survivalstan.fit_stan_survival_model(\n", " model_cohort = 'model 4',\n", " model_code = model_code4,\n", " df = d,\n", " time_col = 't',\n", " event_col = 'event',\n", " formula = '~ age_centered + sex',\n", " iter = 5000,\n", " chains = 4,\n", " seed = 9001,\n", " FIT_FUN = stancache.cached_stan_fit,\n", " drop_intercept = False,\n", " )\n" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "# 0:00:16.703755 elapsed" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " mean se_mean sd 2.5% 50% 97.5% Rhat\n", "lp__ -272.823550 0.022981 1.461548 -276.488239 -272.492426 -271.021688 1.000197\n", "alpha 1.023161 0.002054 0.150942 0.749636 1.011858 1.346196 1.000033\n", "beta[1] -3.030795 0.003943 0.275100 -3.610841 -3.013843 -2.531514 0.999949\n", "beta[2] 0.615333 0.002998 0.240180 0.152485 0.613556 1.096432 0.999838\n", "beta[3] -0.003085 0.000159 0.013548 -0.030261 -0.002951 0.023048 1.000245\n" ] } ], "source": [ "survivalstan.utils.print_stan_summary([testfit4], pars=['lp__', 'alpha', 'beta'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## compare coefficient estimates for each model spec" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhMAAAEKCAYAAACoktfqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xt8k+X5P/DPlTb0QJFKW0VbMZQ2YpXiABHn1g0PDIq4KZtz4qwTh3NK/Q5Utq8tAwpfRcR9VzwyUYpf3UScQwVU/M3DdDpWEAER0oItFC2QKofSlh5y/f5I0rUlaZImaZLm83698mry3Pfz5MrTQq/ez/3cl6gqiIiIiHrKEOoAiIiIKLIxmSAiIiK/MJkgIiIivzCZICIiIr8wmSAiIiK/MJkgIiIivzCZICIiIr8wmSAiIiK/MJkgIiIiv8SGOoBol5qaqiaTKdRhEBFFlM2bN1tVNc3PY5wRGxv7NIALwT+uu2MDsKO1tfW20aNHH3LVgclEiJlMJpSXl4c6DCKiiCIi1f4eIzY29unBgwefn5aW9o3BYGBtCTdsNpscPnw4p7a29mkA17jqw0yMiIii1YVpaWnHmEh0z2AwaFpa2lHYR3Bc9+nFeIiIiMKJgYmEdxznyW3OwGSCiIiI/MI5E0RERACm3zHzvLpvjvYL1PFSTh/YvOKJZbsDdTxP0tPTR5SXl39+1llntfraZ+bMmekvvfRSyrFjx2IaGho+8fW9mUwQEREBqPvmaL/951wZsGQC+98O2KGC7Uc/+tGRe+6559D555/vdl5Ed5hMUEQqLS1FZWVlt31qamoAABkZGW77ZGVlobCwMKCxERF5Y/fu3f0mTpyYPWrUqBObN29Oys3NPXHrrbdaFyxYkF5XVxe7cuXKvePHj284ePBgzLRp00z79u2LS0hIsC1fvrz6kksuaaytrY2ZOnVq5sGDB/uNHj26XvU/0z8ef/zxQU888cSZLS0tMmrUqBOrVq2qjo11/yv/iiuuOOHPZ2EyQRGpsrISn2zfCVviILd9DA1HAQAHT7r+MTc0fB2U2IiIvLV///74F198ce/o0aOrcnNzz3/++edTysvLd73wwgvJixYtOmv8+PF77rvvvrNHjhzZ8Pbbb+959dVXBxQUFAzdtWvXzt/+9rdnX3rppfUPP/zwV3/5y18Grl69OhUAtmzZEr9mzZpB5eXlu+Li4vSmm24a8uSTT6bcdddddcH6HEwmKGLZEgehKedqt+3xO18HALd9nO1ERKGSnp5+cuzYsY0AYDabGy+//PJjBoMBo0aNali4cOHZALBp06YBL7/8ciUAXHPNNcdnzJgR+/XXXxs+/vjjAX/9618rAeCGG244evvtt7cBwBtvvDFgx44diSNHjjwfAJqamgxnnHGG23kUgcBkggKutLQUAHj5ADwXRNS9fv36tV+bMBgMiI+PVwCIiYlBW1ub9OSYqio/+clP6h577LEDgYrTEyYTUcxqtWL+/PmYN28eVLX9eUpKCqxWK4qKigAAs2fPxtKlS9HU1IQvv/wSycnJqK2tRVZWltt5C/Hx8fwFCmDLli3Yu3cv1qxZA6PRiNbWVgwdOhRLly5FSkpKqMMjoghwySWXHH/22WdTlixZ8tXrr78+4PTTT28dNGiQbdy4ccdXrlyZ8tBDD321evXq044dOxYDABMnTjx23XXXZf33f//3wfT09NaDBw/GHD16NMZsNjcHK0YmE1GsrKwM27ZtQ1lZGVS1/fmsWbNQVlaGnTt3AgBKSkpQVVXVvl9tbS0AdDsBsqmpKaixR4qO562lpQUAsHfv3vbzTEThI+X0gc2BvAMj5fSBAfnlvXjx4i+nTZtmMpvNOQkJCbaVK1d+AQAPPvjgl1OnTs3Mysq6YMyYMfVnnXVWMwCMHj26qaio6MAVV1xhttlsMBqNWlpauq+7ZOJXv/pVxiuvvDKoqanJcOaZZ+ZOmzbN+sgjj3zpbYzScfYn9b4xY8ZoKGpzWK1W3HDDDWhubka/fvY7oZqbmxEXF4fHH38cv/rVr9p/+fXUiBEj8NhjjwUi3FMUFhZi855av+dMjB42uP1SRKBt2rQJ99xzj8u22NhYvPTSSxydIOohEdmsqmP8Ocann35aNXLkSGugYurrPv3009SRI0eaXLVxZCJKOUcjAHRKGmw2G0pKStDa6v9cne3btwftUkdFRQWk2b9EWJqOoaLieNBi3LZtm9u21tZWjk4QUZ/B5bSj1MaNG9uTCFXtlFhUVVWBI1b+s9ls3ba/9dZbvRQJEVFwcWQiSl111VVYv349WlpaIGKfMKyqMBqNSE9PR3V1dUASimBdQnBe5vCHxp+G7CBe5sjPz0d9fb3b9gkTJgTlfYmIehtHJqJUQUFBexJhNBphNBoB2G9NKi4uRncrpXlrxIgRfh8jks2bN89tW2xsLAoKCnovGCKiIGIyEaVSU1MxadIkiAjy8/Pbn0+aNAnZ2dnIz89v72symXr0HsGafBkpxo4d256kdXX11Vdz8iUR9RlMJqJYQUEBcnNzUVBQ0Om5sy0nJwc5OTkoLi5GTk4OMjMzER8fj8GDBwOw17VwJ9pHJZweeOCBTiNAIoLMzEyOShBRn8I5E1EsNTUVy5Yta3/d8XlqaiqefPLJ9tcdn5P3xo4di/feey/UYRCRF2b9evp59UesAasampSc2vzI4yvCvgT58ePHDVOmTMmsrq6Oi4mJwYQJE448/vjjPq2eyWSCIlJNTQ0Mx+uQWF7muoOtDVAbAHHfp60VNTVBXa6eiCJI/RFrv6LhewOWTCzcFagjBd/s2bMPTpky5XhTU5Ncdtll5tWrV592/fXXH/N2f17moIjVP7YN5w9ocPlINLTCKIpB/Vrd9ukf2xbqj0BEUWz37t39hg4desHUqVNNJpPpwmuuuWbo3/72twGjRo0afu655174zjvvJALAwYMHY6688sphZrM5Z+TIkcP/9a9/JQBAbW1tzGWXXZadlZV1wU9/+tNzu5YgHzFixPnDhw/PufHGG8/tbu2gAQMG2KZMmXIcAOLj4zU3N7dh//79PiVVHJmgiJSRkYGm1q9QNMb1rZcLy5MAwG27s098RkZQ4iMi8ka4lSC3Wq0xGzduTL733nsP+vI5ODJBREQUIs4S5DExMaeUIK+pqYkD7CXIp0+fXgfYS5AfOXKkvQT5rbfeWgfYS5Cfdtppp5QgHz58eM4HH3xw2t69e+M8xdLS0oLrrrsuc8aMGQdzcnJ8qivCkYk+pmMl0K63HlosFhQWFiIjIwOLFy9GSkoKLBYL7r77bixbtgzJycmYPXs2qqurMXfuXLzwwguorq5u399gMGDRokVYunQpDhzwPDdn/vz5GD9+fMA/YyQpLS3Fvn37sGnTJgDAsGHD8PDDD7dXZnX3vSKi6BBOJchvvPFGU2ZmZtPcuXMP+fqeHJnoYzpWAu1q4cKFaGhogMViaW9fuHAhTpw4gQULFqCsrAxffPEFbDYbFi5cCIvFgpMnT7Y/GhsbMXfuXK8SCcBebTTaVVZWticSALBnz572c9/d94qIyMlZghwAXJUgB4CuJchff/310w8cOBAL2OdcWCyWbudAFBYWnn3s2LGYFStW7O9JjByZ6EOsVis2bNgAVcWGDRtQUFDQ/hevxWLpVA573bp1+O53v9u+raqqqtMohLvJOt0tD91Va2sr3nnnnagenfj6669P2bZu3TpMmTLF7feKiEIjKTm1OZB3YCQlp0ZECfI9e/YYly1bdtbQoUObLrjgghwAmDFjxqFZs2Z5XVGVJchDLJAlyJcuXdpeb8NoNGLy5MntVSlvvvnmTskEACQlJfmUHPTURRddFPBjVlRUoF/LMSzLc33nkjcTMGe+fxqajachOzs74PE5bd261eV2k8mEAwcOuPxeEZFnLEHe+7orQc7LHH1Ix0qgLS0tnapSdk0kAN9GGSiwqqqq3H6viIgiDS9z9CEdK4EajcZOVSlNJlNIRiZiY2ODUpWzsLAQTVX/9usYZybaEG/KDlrVUADIy8tzub3ryAQriBJRJOPIRB/SsRKowWDoVP+hqKioU1+j0Yj58+d32ubcN5CKi4sDfsxIMmTIkFO2GY1GFBcXu/1eERFFGiYTfUjHSqCTJk3qNKHPbDZ3qv45efJkXHzxxe3bTCYTrrnmmvZ2dyXIk5KSvI4nNjY2qidfAsCgQYNO2TZ58mRkZ2e7/V4REUWaiEwmRMTj2LyI/JeIJPZGPOH0/l2rf3ZUVFSExMREmM3m9vaioiL0798fc+fORUFBAYYOHQqDwYCioiKYzWbExcW1PxISErBgwQKkp6d7FUu0j0oA9sqqY8eObX89bNiwTpVZ3X2viIgiSUTezSEi9ara7Z/IIlIFYIyqen9ri0iMqgakYIO37x/IuzmiiXPOhN/LaZsuDuqcCSIKDt7N0fu6u5sjoidgisj3AcwDYAVwIYDNAG4CMBPA2QDeERGrqo4XkQkA5gOIA7AHwC9Utd7xS/9FAFcBeEhEygE8CSANQBuAn6jqHhG5F8D1jv1fUdXfi4gJwBuO9x0F4DMANwO4rev7B/lURKV99THtSUNX1cdjAMBtu3N/c1AiI6JIdNudt51Xd7QuYFVDUwamND/92NNhX4IcAL773e9mHzp0yNjW1iZjx449vmrVqn3uLne7EtHJhMO3AFwA4EsAHwK4TFVLRWQWgPGqahWRVABFAK5U1RMiMgfALAALHMeoU9VRACAi/wLwoKq+IiLxAAyORCQbwFgAAuBVEckDsA/AeQCmq+qHIvIMgF+r6sMd37+XzkNUqa+vh8QNwP4W1+0tchJtbW2orO+HuDjXS9InJicgKysriFESUSSpO1rX79AlhwKWTOBfATtS0K1du3bPoEGDbDabDZMmTRr2zDPPnD5jxoxvvN0/IudMdLFJVWtU1QZgKwCTiz7jAOQA+FBEtgIoAHBuh/YXAUBEBgBIV9VXAEBVm1S1AcAEx+MTAFsADIc9uQCA/ar6oeP5/wH4jqeARWSGiJSLSPnhw4d9+rBkl5SUhPrmetQnun60nN4CW6wNLdLiuk9zPTIyMlBYWBjqj0JEUSpcSpADwKBBg2wA0NLSIi0tLeLr3X19IZk42eF5G1yPtgiAjap6keORo6rTO7Sf8PAeAuCBDvtnqeoKR1vXSSceJ6Go6nJVHaOqY9LS0jx1J3eSAdv3bW4fSHbfB8mhDp6IyF6CfM6cOQf37NmzY8+ePfHOEuSLFi2qWbRo0VkA4CxBbrFYdpaUlBwoKCgYCgDOEuSVlZWfXXvttUe++uqrfkDnEuS7du3aaTAY9Mknn/R4y9h3vvOd7LS0tJH9+/dv+8UvfuH1qATQN5IJd44DGOB4/jGAy0QkCwBEpL+InHK5XFWPA6gRkR85+sU57sh4E8CtIpLk2J4uImc4dhsiIpc6nt8I4AMX7x91SktLo35iI88BEXkSTiXIP/jgg4ra2tpPm5ubDa+99tppvnyOvjBnwp3lAN4QkS8dEzBvAfBnEXGe0CIAFhf7/RzAUyKyAEAL7BMw3xKR8wF85Bj6qYd9omcbgN0A7nTMl9gJ4AlX7x+cjxh4VqsVc+bMQU1NDc466ywA6FQltK2tzW0RsK769+8f1ZcRPv/8c+zZswcvv/wyHnnkEYwePbpTyXfO1yCicCpBDgCJiYk6ZcqUI6+88krytdde67r4kQsROTLhvC1UVd9V1as7bL9LVVc6ni9T1fOcv8hV9e+qerGq5joerzq2mzpOklTVClW93NFntKrudWz/o6qOcDwuVdU9jl1aVfUmVT1fVac65lic8v6RoqysDBUVFWhsbMTevXuxd+/eTmXIvU0kAODECU9Xj/q22tpaNDU1QVXb19zoWPKdiMgbwS5BfvToUUN1dbURsNcK2rBhw8Dhw4c3+hJjXx6ZIB9ZrVasW7cuoMe877778NBDDwX0mJHAarV2Kj9eX1+PtWvXdir5XllZydEJojCSMjClOZB3YKQMTImIEuTHjh0zTJ48Oau5uVlUVb797W8fu/fee326OyAiF63qS8Jp0aqlS5di7dq1AT9usEqQ19vqYbva5raP4V37wJvt+6f2MbxuQJIhKWjlx/fv34+6urpO20QEHf+9mUwmrFq1KijvT9TXcdGq3scS5OSVjRs3hjqEPuObb06dCN01cXdVFp6IKBLxMge1u+qqq4IyMhGsEuSfHPik5wdIArLTg1d+3NUoj6uRCSKivoAjE9SuoKDAbbXQnho3blxAjxcpOpaDd5o1a1an13Pnzu3NkIiIgobJBLVLTU3F5MmTA3rMaJx8CdjP5dChQ9tfJyUl4Yc//GGnku+cfElEfQWTCeqkoKAA2dnZSEhIQGZmJjIzMzuVIfdl5CJaRyWcHn74YWRmZkJEUFJSAqBzyXcior6Ccyaok9TUVKxYscJzR/IoNTUVK1eu7LTNbDZjw4YNoQmIiChImExQRKqpqQHqAMPfuhlcc1QUddmnFahoqEBpaWlUr9JJRP9x922/PO94nTVgVUMHpKQ2//HpP0VECXKnyy+/PGv//v1xFRUVn/ny3kwmKGIlqGCwmxLkALDf8fUcN31qT5xAZWVlwOMiosh0vM7ar6Du64AlE2WBOlAvKSsrS+7fv39bT/ZlMkERKSMjA0esVkyH+6XrVzgKuLrrs8JzgVcioqDZvXt3v4kTJ2aPGjXqxObNm5Nyc3NP3HrrrdYFCxak19XVxa5cuXLv+PHjGw4ePBgzbdo00759++ISEhJsy5cvr77kkksaa2trY6ZOnZp58ODBfqNHj67vWoL8iSeeOLOlpUVGjRp1YtWqVdXdzXk7evSoobS09Mzly5dX33DDDcN8/SycgElERBQi4VKCfNasWel33333waSkJPfLCneDyUQEs1qtmDlz5inLNnu73V0/ANi0aRO+973v4aabbsItt9yCiRMn4sEHH0ReXp7Xj1dffTV4Hz4ClJaWYvHixZg+fTomTJiAH/zgBygvL8fMmTNRUVHh9twTUfQIhxLk//znPxO++OKLuJtvvvlITz8Hk4kIVlZWhm3btqGsrKxH2931A4B58+ZBVbFv3z7s3bsXDQ0NWL9+vU/xLV261MdP1LdUVlbi448/RkVFBZqamtDY2Ii5c+di27ZtKCkpcXvuiSh6BLME+a5du3bu2rVrZ1VV1Y5HHnnkS3f9//GPfyTt2LEjMT09fUReXt7wqqqquLFjx57ny3symYhQVqsVGzZsgKpiw4YNnUYbvNleUVHhsh9gH5Wor6/3O0ZVjerRiZaWllNGHurr66GqqKqqcnnuiYi6CnYJ8jlz5hw+dOjQtgMHDmx///33d5lMppObNm3y6S4UTsCMUGVlZe11Hmw2G8rKyjBr1iyvt5eUlLjsB9hHJQLl4Ycfxttvvx2w4zlVVFT4nQnXAThcURG0W0N37/b8b7HruSei0BmQktocyLHCASmpEVGCPBBYgjzEelqCfOLEiWhoaGh/nZiYiDfeeMPr7V05+wFAXl6ez/F0J1glyA0nTuA+P+7meAgKW//+QStDvnXrVq/6dTz3ROQdliDvfd2VIOfIRIS66qqrsH79erS0tMBoNGLChAk+bU9PT8eBAwdO6QfY60gE4jIHYK+UGayqoUe8/GXtTgqA5OzgVQ699tprPV7C6HruiYgiEedMRKiOVSkNBgMKCgp82l5cXOyyHxDYyxyzZ88O2LEizeDBgz326XruiYgiEZOJCJWamopJkyZBRDBp0iSkpKT4tD07O9tlPwAYO3YskpKS/I5RRHDNNdf4fZxIZTQaO51XwD7qIyIwmUwuzz0R9SqbzWbr0R0T0cZxntyuQcFkIoIVFBQgNzf3lL9svd3urh9gH50QEQwZMgSZmZlITExEfn6+T/FF86gEAGRlZWHcuHHIzs5GfHw8EhISsGDBAuTm5qK4uNjtuSeiXrPj8OHDA5lQdM9ms8nhw4cHAtjhrg8nYIZYTydgRjvnnAl/l9NOvuiioM2ZIKLgCcQEzM2bN58RGxv7NIALwT+uu2MDsKO1tfW20aNHH3LVgRMwKWLVovv6Gl85vrrrUwsgOeBREVGkcPxijN5rsQHEZIIiUlZWlsc+9TU1AIDkjAyX7cleHoeIiLrHZIIiUrAWmiIiIt/xGhERERH5hckEERER+YXJBBEREfmFyQQRERH5hckEERER+YXJBBEREfmFyQQRERH5hckEERER+YXJBBEREfmFyQQRERH5hckEERER+YXJBBEREfmFhb6IiAKotLQUlZWVLttqHJVsM9xUss3KymIRO4pITCaIiAKosrISn2zfCVvioFPaDA1HAQAHT576X6+h4eugx0YULEwmiIgCzJY4CE05V5+yPX7n6wDQbRtRJOKcCSIiIvILRyaIKGxZrVbMnz8f8+bNQ0pKSo/2nzNnDmpqatC/f39YrdZT+px77rl47rnnAhFu2LnttttgsVg89ktLS8Px48chIrjvvvuwZMkSZGRkYPHixT067xR9ODJBRGGrrKwM27ZtQ1lZWY/3r6ioQGNjo8tEAgCqq6v9CTGseZNIAMDhw4fR1NSExsZGLFq0CA0NDbBYLD0+7xR9mEwQUViyWq3YsGEDVBUbNmxAXV2dz/uvX7/eq75LlizpSYhh7amnnurRfq2tre3P161b5/N5p+jEyxxEFJbKysqgqgAAm82GsrIyzJo1y6f9W1pavOr72muvYf/+/T2Ks6uKigpIs/q8nzQdQ0XF8YDdGrp161a/j9HS0uLzeafoxJEJIgpLGzdubE8GWlpa8NZbb/m8P/nP1/NO0YkjE0QUlq666iqsX78eLS0tMBqNmDBhgs/7r1271uv+paWlvoboUmFhITbvqfV5P40/DdnDBgcsjry8vIAcx9fzTtGJIxNEFJYKCgogIgAAg8GAgoICn/c3Go1e9Z0yZYrP8YW7adOm+X0Mo9Ho83mn6MRkgojCUmpqKiZNmgQRwaRJk3y+RTE1NRX5+fle9b333nt7EmJYu/3223u0X2zsfwasJ0+ezFtDyStMJogobBUUFCA3N7fHfx0XFBQgOzsbCQkJSE1Nddnn3HPP9SfEsGY2m73ql5aWhvj4eCQkJOD+++9HYmIizGYzRyXIax7nTIjImQD+B8DZqjpJRHIAXKqqK4IeHRFFtdTUVCxbtsyv/VesiN7/qp5++uke7XfFFVcEOBLq67yZgLkSwLMA7ne8tgB4EUD0/gslIuqGoeFrl7U2DA32NRtct30NYHCwQyMKCm+SiVRVXS0ivwMAVW0VkbYgx0VEFHFKS0tRU1ODpHgjYDt+SnuDzb4gVJKLNsQbUV9fH+wQiYLCm2TihIikAFAAEJFxAI4GNSoioghUWVmJhiOHMCTJ9d9b1YYYAMA5xiOntO2rj0FSUnZQ4yMKFm+SiVkAXgUwTEQ+BJAG4MdBjYqIKEINSWpD0RjXIwwLy5MAwGW7s40oEnlMJlR1i4h8D8B5AATAblX1bo1aIqIw5VwcKlDLV4ejaPiMFB7cJhMicp2bJrOIQFX/GqSYiIj85ixfXlhYiNLSUtx888343e9+h+bm5vY+IoL8/HxkZWWFMNLg2bZtGywWC9asWXNKW0JCAhobG9tfx8XFIT09vX2didjYWCxatKjH60x4Kh/vb3l5Ci/drTMxpZvH1cEPjYio55zly0tKSrBt2zb8/ve/75RIAICqYsGCBSGKMPi6K6/eMZEAgJMnT2Lv3r2wWCywWCzYuXOnXyXIPZWP97e8PIUXt8mEqv6im8etvRkkEZEvOpYvr6qqgqq6vVOiqqoKlZWVvRxh8FksFpw8edKvY6xfv75HJcg9lY/3t7w8hR9vFq1KAfB7AN+B/Y6ODwAsUFV+94koLHUsX+6NO+64A+eff77f71tRUYF+LT1bWPhggwHNFRUBm9/w+eef+32MnpYg91Q+3t/y8hR+vPmp/wuAwwCmwn4Xx2HYF60iIgpLHcuXe8Pfv+DDUSA+k6r2qAS5p/Lx/paXp/Djza2hZ6lqSYfXC0Xkp8EKiIjIXx3Ll3vDZDIFpPR3YWEhmqr+3aN9z0y0Id6UHbAS5DfffDOqqqr8OoaI9KgEuafy8f6Wl6fw483IxFsicoOIGByP6wG8GezAiIh6qmP5cm/MnTs3iNGERlFRkd/H6GkJck/l4/0tL0/hx20yISLHReQYgF8CeAFAs+PxFwAzeic8IiLfdSxfbjKZICJISnK9KJTJZOqTt4aazWYMGDDAr2Pk5+f36LZNT+Xj/S0vT+Gnu7s5BqjqaY6vBlWNdTwMqnpabwZJROQrZ/ny4uJi5ObmYv78+ejXr1+nPnFxcX1yVMLpD3/4g9sRmoSEhE6v4+LikJmZCbPZDLPZjJycHL9GDDyVj/e3vDyFF2/mTEBETgeQDSDeuU1V3w9WUOFGRPoDWA0gA0AMgBIAlQAeAZAEwArgFtgnp34E4F5VfVdEHgBgU9X7XR2XiIKnY/ly59e33347lCH1OrPZjPfeey8k7+2pfLy/5eUpvHhza+htAO6G/RfpVgDjYP+FeXlwQwsrEwF8qaqTAUBEBgLYAOCHqnrYMSF1kareKiK3AFgjIjMd+10SqqCJqHfV1NTg6yOx+OU7A122N7bZRwlctZ9sEwyqqQlqfETB4s3IxN0ALgbwsaqOF5HhAP4nuGGFne0AlorIYgCvA/gGwIUANjqGEGMAfAUAqvqZiDzn6HepqjZ3PZiIzIBj3smQIUN65QMQUe9oMwINyW4mfzpW52kY5KL91EKiRBHDm2SiSVWbRAQiEqequ0TkvKBHFkZU1SIiowDkA1gI4O8APlPVS93sMgL2/xrOcHO85QCWA8CYMWO8X1mHiMJaRkYGDsth2L5vc9lueNc+Tc1Vu+FdAzLSM4IaH1GweHNraI2IJAP4G+x/ia8F4H7B9z5IRM4G0KCq/wdgCeyXLtJE5FJHu1FELnA8vw7AIAB5AJY5zh0REVGf5U0J8msdT+eJyDsABgJ4I6hRhZ8RAJaIiA1AC4A7ALQCKHXMn4gF8L8ichDAgwCuUNX9IvIogD8C4HRloiBzVqH80Y9+hJKSEvzmN7/BY489BpvNhpbbV6e4AAATAElEQVSWFpfLaw8bNgzPPvtsCKLtPYsXL8bf//73Uwp7AfY7OH77299iyZIlyMjIwOLFi5GSktKpomddXR1mzpyJc845B3PmzEFpaSnmzZsHVWXVT2rXXQny01T1mIgM6rB5u+NrEoCvgxpZGFHVN+F6oa48F9vMHfYLzFJ2ROSRswrlZ599BpvNhkceecRjfY49e/b0UnSh8/HHH7tMJAD7ktuLFi1Ca2srLBZLe42MjhU9t27disbGRlgsFpSUlKC6urq9toazD+tqUHeXOV5wfN0MoNzFVyKisNCxCmVraysAeF3o69FHHw1maCFltVo9VuR0ni8AWLduHSoqKtrP5bp16zotye2swLp+/XqsX7+eVT+pnduRCVW9Wuy3KnxPVff1YkxERD7xtUpoR6tXr4bFYglIHBUVFYDruZee1dv3D1TVUADYv3+/T/1bWlpQUlLSfi7d1TbpuJ1VPwnwMAFT7T9R63opFiKiHvG1Smi0+Oabb3zep6qqyuO5VNVOCQerfpI3t4ZuEZGLVbVnpfCIiILM1yqhXQWqUmdhYSE+OfBJz3ZOArLTA1c1FACWLl2KtWvX+rSPyWTCgQMHuj2XziW6VZVVPwmAd7eGXgLgIxHZIyLbRGS7iGwLdmBERN7ytUpoR9dff32Aowkfvta9MBqNKC4ubj+XRqPRbb/YWPvfoqz6SYB3ycQPAAyDffnsKQCudnwlIgoLHatQOn/JeZtc3HXXXcEMLaRSU1M93rbpPF8AMHnyZGRnZ7efy8mTJ8NkMrW3Oyuw5ufnIz8/n1U/qZ3HZEJVq1W1GkAjAO3wICIKG84qlPfffz8MBgNmzZqF+Ph49OvXz21iMWzYsF6OsveNGzfulAqhTnFxcbj//vuRmJgIs9ncPsLQsaJnUVEREhISYDab2yuwFhQUsOondSKeZkCLyDUAlgI4G8AhAOcC+FxVLwh+eH3fmDFjtLycd9oS9QXOORM9XU77W+nfCuicib5MRDar6phQx0F23kzALIG9UujbqvotERkP4KbghkVEFKGO/CdpcNUGuGk/AiA9aFERBZU3yUSLqtaJiEFEDKr6joj8b9AjIyKKMFlZWaipqUFjg+sVJxtaGwAAiQ2Jpzb2A+rr64MZHlHQeJNMHBGRJAD/APC8iBwCcCK4YRERRZ7CwkJUVlZi19atGOyi/SvH17QTp/4XWgsgKSkpmOERBY03yYSzuNfdsF/eGAhgQTCDIiKKZIMBTMepkz5XOOaud9dGFIm8uTU0FsBbAN4FMADAi6rKhdiJKKKVlpb2+cmO0fAZKTx4U4J8PoD5IpIL4KcA3hORGlW9MujRERH5yFk+u7CwEEuXLkVjYyMOHDiA5uZml/3T0tLws5/9rJej7B1btmzB3r17sWbNGpftsbGxaG1txezZs/HSSy9h3759GDJkCP74xz9CVVFUVISmpiZ8+eWXGDJkCObMmdN+Tg8ePIhHH30UWVlZvfypKBx5c5nD6RDsl/XqAJwRnHCIiPzjLJ9dUlLSqeKlO0888USfTSY8fX5nxdCO5dr37dvXXjht586d7X2dJcg7HnPBggVYtWpVwOOmyOPxMoeI/FpE3gXw/wCkAPilquYGOzAiIl91LEXuTSLh9Oc//zl4QYXIpk2bYLN5V8K063pDr732GtatO7XGY9dzWlVVhcrKyh7HSH2HNyMT5wD4L1XdGuxgiIj80dNS5E888QQ++uijgMRQUVHh1WS0ruoAHA5gCfJt23peQqmtrc3rvhydIMC7ORO/641AiIj8xVLk/+HtqIS/fBkBor7LlzkTRERhzZ9S5IEsQ35kq+8DuSkAkrMDV4I8Pz+/VxbB6lgIjKJXT0bjiIjCUk9Lkd9xxx1BiCa05s2b1+N9Y2JiOlUT7c7cuXN7/D7UdzCZIKI+o2Mpcl/+Yu6Ld3OMHTsWRqPRq75dE7ApU6Zg8uTJp/Trek5NJhNvDSUATCaIqI9xlsYuLi5GTk4Ohg4din79+rnt3xdHJZweeOCBbtudow+zZs3CkCFDAABDhgxpLzGek5ODzMxMxMfHt5cgd57TxMREjkpQO48lyCm4WIKcqG9xzpnoyXLayRddxBUrvcQS5OGFEzCJiAKsFq5rbTgLfblqqwWQHNSoiIKHyQQRUQB1N4egvqYGAJCckXFKW7KHfYnCGZMJIqIACtSiU0SRhBMwiYiIyC9MJoiIiMgvTCaIiIjIL0wmiIiIyC9MJoiIiMgvTCaIiIjIL0wmiIiIyC9MJoiIiMgvTCaIiIjIL0wmiIiIyC9MJoiIiMgvTCaIiIjILyz0RURRqbS0FO+++y6Sk5PxzDPPhDocoojGZIKIolJlZSWsVisaGxtDHQpRxONlDiIiIvILkwkiijqlpaWoqakBAJw8eRKlpaUhjogosjGZIKKoYrFYsGbNGlitVgBAa2sr1qxZg7y8PLzzzjshjo4oMjGZIKKosnDhQrdtJSUlvRgJUd/BZIKIoobFYkFVVZXb9tbWVo5OEPUAkwkiihrdjUo4cXSCyHdMJogoanQ3KuHU2toa/ECI+hgmE0QUNUwmk8c+sbFcfofIV0wmiChqFBUVeexTXFzcC5EQ9S1MJogoapjN5m5HJ2JjYzF+/PjeC4ioj2AyQURRpbvRCY5KEPUMLw4SUVQxm814//33UVhYiK1bt6J///7YsGFDqMMiimgcmSCiqORcTruhoSHEkRBFPiYTRERE5BcmE0QUlTIyMgAAiYmJIY6EKPIxmSAiIiK/MJkgoqjy1FNPIS8vD1u3bgUAnDhxAnl5ecjLy8OKFStCHB1RZGIyQURR5fnnn3fbVlZW1ouREPUdTCaIKGo89dRTHvtwdILId0wmiChqdDcq4cTRCSLfMZkgIiIivzCZICIiIr8wmSCiqDFt2jSPfQoKCnohEqK+hckEEUWN22+/3WOf6dOn90IkRH0LkwkiiirdjU5wVIKoZ0RVQx1DyIjI9wE0q+o/g3T8lQBeV9U17vqMGTNGy8vLg/H2RNQNVg2NbCKyWVXHhDoOsov2kYnvA/i2LzuICMu2ExERdRDUZEJE/iYim0XkMxGZ4dg2XUQsIrJJRP4kIo86tqeJyMsi8m/H47JujpskIs+KyHYR2SYiUx3bJ4jIRyKyRUReEpEkx/YqEZnv2L5dRIaLiAnArwD8RkS2ish33cUgIvNE5DkR+RDAcyISIyJLHH22icjtjn4iIo+KyG4ReRvAGcE7u0Tkj/r6eogIWlpaQh0KUcQL9l/Zt6rq1yKSAODfIrIOQDGAUQCOA/g7gE8dff8I4A+q+oGIDAHwJoDz3Ry3GMBRVR0BACJyuoikAigCcKWqnhCROQBmAVjg2MeqqqNE5NcA7lHV20TkSQD1qvqw4zgvdBNDDoDvqGqjIzE6qqoXi0gcgA9F5C0A3wJwnqPvmQB2Anima/CO/WcAwJAhQ3w6oUQUGElJSVBVGI3GUIdCFPGCnUwUisi1jufnAPg5gPdU9WsAEJGXAJgd7VcCyBER576niUiSqta7OO6VAG5wvlDVb0Tkath/iX/oOEY/AB912Oevjq+bAVznJl6XMTiev6qqjY7nEwDkisiPHa8HAsgGkAfgz6raBuBLEfm7qzdR1eUAlgP2ORNuYiEiIooIQUsmHJMbrwRwqao2iMi7AHbB/WiDAcA4VW3q6VsC2KiqP3PTftLxtQ3uP7fLGBzJxYku7zVTVd/s0i/f16CJqPeVlpaipqYGAHDy5EmUlpaisLAwxFERRa5gzpkYCOAbRyIxHMA4AP0BfM9xWSIWwNQO/d8CMNP5QkQu6ubYGwHc2aHv6QA+BnCZiGQ5tvUXEbOb/Z2OAxjQgxjeBHCHiBgd/cwi0h/A+wB+6phTcRaA8R7en4h62XPPPYc1a9bAarUCAFpbW7FmzRrk5eXhzjvv9LA3EbkSzGTiDQCxIvI5gAdh/2V/AMD/ANgE4EMAVQCOOvoXAhjjmNC4E/bJke4sBHC6iOwQkU8BjFfVwwBuAfBnEdkG+yWO4R5ifA3Atc4JmD7E8DTs8yG2iMgOAE/BPtrxCoAKR9sqdL7MQkRh4E9/+pPbtu3bt/diJER9R6+vM+GcB+EYmXgFwDOq+kqvBhFGuM4EUe957rnnuk0mAGDEiBF47LHHeiki6imuMxFeQrHOxDwR2QpgB4AvAPwtBDEQURTylEgAHJ0g6oleX4BJVe/xtq+I/ALA3V02f6iqvLBJREQUJsJ6NUdVfRbAs6GOg4iIiNyL9uW0iSiK/PKXv/TYZ8SIEb0QCVHfwmSCiKLGz3/+c499OPmSyHdMJogoqnQdneiw4i1HJYh6KKpLkIcD3hpKFBosQR7ZeGtoeOHIBBFFJedy2g0NDSGOhCjyMZkgIiIivzCZIKKolJGRAQBITEwMcSREkY/JBBEREfklrBetIiIKNIvFgttuu6399YkTJ5CXlweAdTmIeoojE0QUVRYuXOi2jXU5iHqGyQQRRQ2LxYKqqqpu+9x5J0v/EPmKyQQRRY3uRiWcODpB5DsmE0QUNTyNShBRzzCZIKKoYTKZQh0CUZ/EZIKIokZRUZHHPqzPQeQ7JhNEFDXMZrPH0QneGkrkOyYTRBRVuhud4KgEUc9w0Soiiipmsxnvv/8+q4YSBRBHJoiIiMgvHJkgoqiUlZWFmpoaJCcnhzoUoognqhrqGKLamDFjtLy8PNRhEBFFFBHZrKpjQh0H2fEyBxEREfmFyQQRERH5hckEERER+YXJBBEREfmFyQQRERH5hckEERER+YW3hoaYiBwGUB3qOLyUCsAa6iB8FIkxA5EZdyTGDERm3JEYMxDYuM9V1bQAHYv8xGSCvCYi5ZF2X3ckxgxEZtyRGDMQmXFHYsxA5MZNnvEyBxEREfmFyQQRERH5hckE+WJ5qAPogUiMGYjMuCMxZiAy447EmIHIjZs84JwJIiIi8gtHJoiIiMgvTCbIZyIyW0RURFJDHYs3RKRERLaJyFYReUtEzg51TJ6IyBIR2eWI+xURiYg62SLyExH5TERsIhLWs/ZFZKKI7BaRShH5bajj8YaIPCMih0RkR6hj8ZaInCMi74jITsfPxt2hjokCj8kE+UREzgEwAcC+UMfigyWqmquqFwF4HcDcUAfkhY0ALlTVXAAWAL8LcTze2gHgOgDvhzqQ7ohIDIDHAEwCkAPgZyKSE9qovLISwMRQB+GjVgCzVTUHwDgAd0bIuSYfMJkgX/0BwH0AImayjaoe6/CyPyIgdlV9S1VbHS8/BpARyni8paqfq+ruUMfhhbEAKlV1r6o2A/gLgB+GOCaPVPV9AF+HOg5fqOpXqrrF8fw4gM8BpIc2Kgq02FAHQJFDRH4I4ICqfioioQ7HJyKyCMDNAI4CGB/icHx1K4AXQx1EH5MOYH+H1zUALglRLFFDREwAvgXgX6GNhAKNyQR1IiJvAxjsoul+AP8N+yWOsNNd3Kq6VlXvB3C/iPwOwF0Aft+rAbrgKWZHn/thHyZ+vjdj6443cRN1JSJJAF4G8F9dRgupD2AyQZ2o6pWutovICABDAThHJTIAbBGRsapa24shuuQubheeB7AeYZBMeIpZRG4BcDWAKzSM7uH24VyHswMAzunwOsOxjYJARIywJxLPq+pfQx0PBR6TCfKKqm4HcIbztYhUARijqmFfbEhEslW1wvHyhwB2hTIeb4jIRNjnpnxPVRtCHU8f9G8A2SIyFPYk4gYAN4Y2pL5J7H99rADwuao+Eup4KDg4AZOiwYMiskNEtsF+mSYSbk17FMAAABsdt7Q+GeqAvCEi14pIDYBLAawTkTdDHZMrjsmtdwF4E/YJgatV9bPQRuWZiPwZwEcAzhORGhGZHuqYvHAZgJ8DuNzxs7xVRPJDHRQFFlfAJCIiIr9wZIKIiIj8wmSCiIiI/MJkgoiIiPzCZIKIiIj8wmSCiIiI/MJkgohOISL1oY6BiCIHkwkiIiLyC5MJoiggIg+KyJ0dXs8TkSIR+X8iskVEtjsKuXXd7/si8nqH1486lvmGiIwWkfdEZLOIvCkiZ/XKhyGisMNkgig6vAjg+g6vrwdQBuBaVR0FeyXVpeJlOVhHrYVlAH6sqqMBPANgUWBDJqJIwdocRFFAVT8RkTNE5GwAaQC+AVAL4A8ikgfABntZ7jMd2z05D8CFsC/3DQAxAL4KRuxEFP6YTBBFj5cA/Bj28uEvApgGe2IxWlVbHMXb4rvs04rOI5jOdgHwmapeGtSIiSgi8DIHUfR4EfbqmD+GPbEYCOCQI5EYD+BcF/tUA8gRkTgRSQZwhWP7bgBpInIpYL/sISIXBP0TEFFY4sgEUZRQ1c9EZACAA6r6lYg8D+A1EdkOoBwuSrOr6n4RWQ1gB4AvAHzi2N4sIj8GUCoiA2H/v+R/AYR95U0iCjxWDSUiIiK/8DIHERER+YXJBBEREfmFyQQRERH5hckEERER+YXJBBEREfmFyQQRERH5hckEERER+YXJBBEREfnl/wP5P4c2ObsVxAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "survivalstan.utils.plot_coefs([testfit, testfit2, testfit3, testfit4])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "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.6.7" } }, "nbformat": 4, "nbformat_minor": 2 }