{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import autograd.numpy as np\n", "import pandas as pd\n", "\n", "from autograd import grad, hessian\n", "from patsy import dmatrices\n", "from scipy.optimize import minimize\n", "from scipy.stats import norm\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Data pre-processing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read in the *NCCTG Lung Cancer* data from the R package [`survival`](https://cran.r-project.org/package=survival)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "ovarian = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/survival/ovarian.csv', index_col=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Replace dots with underscores in column names." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "ovarian.columns = ovarian.columns.str.replace('\\\\.', '_')" ] }, { "cell_type": "code", "execution_count": 4, "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", "
futimefustatageresid_dsrxecog_ps
159172.3315211
2115174.4932211
3156166.4658212
4421053.3644221
5431150.3397211
\n", "
" ], "text/plain": [ " futime fustat age resid_ds rx ecog_ps\n", "1 59 1 72.3315 2 1 1\n", "2 115 1 74.4932 2 1 1\n", "3 156 1 66.4658 2 1 2\n", "4 421 0 53.3644 2 2 1\n", "5 431 1 50.3397 2 1 1" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ovarian.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exploratory data analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bar plot of censored vs dead" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAACrpJREFUeJzt3V+MZnddx/HPl46NFjBg9gFht+tWbWqQmGAm/iPRhNpYlVAvvGgTTMUmc4WCsUIJicQ7Exv/JBrNBkpJbBaTirHhQmmKTWMsldlSoGWBEsSytLLT1Pj3otR+vdghlmF3Z57nnN3p/vb1SjYzzzlnnvO9mH3nt2fOma3uDgAXv5fs9wAAzEPQAQYh6ACDEHSAQQg6wCAEHWAQgg4wCEEHGISgAwxi7UKe7MCBA33kyJELeUqAi97x48ef7u7Fbsdd0KAfOXIkm5ubF/KUABe9qvqXvRznkgvAIAQdYBCCDjAIQQcYhKADDGLXoFfVHVV1qqoePcO+W6uqq+rA+RkPgL3aywr9ziTX79xYVVcmuS7JEzPPBMAKdg16dz+Q5Jkz7PrDJO9K4v+wA3gRWOnBoqp6S5Kvdfenq2q3YzeSbCTJ4cOHVzndBVd1+36PMJTuW/d7BLgkLP1D0aq6Isl7k/zOXo7v7qPdvd7d64vFrk+uArCiVe5y+YEkVyX5dFV9JcmhJA9X1ffOORgAy1n6kkt3fzbJq775ejvq69399IxzAbCkvdy2eCzJg0muqaqTVXXL+R8LgGXtukLv7pt22X9ktmkAWJknRQEGIegAgxB0gEEIOsAgBB1gEIIOMAhBBxiEoAMMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQQg6wCAEHWAQuwa9qu6oqlNV9egLtv1+VX2+qj5TVX9dVa84v2MCsJu9rNDvTHL9jm33Jnl9d/9Iki8mec/McwGwpF2D3t0PJHlmx7aPdfdz2y8/keTQeZgNgCWszfAev5bkL8+2s6o2kmwkyeHDh2c4HVy6qm7f7xGG0n3rfo8wq0k/FK2q9yZ5LsldZzumu49293p3ry8WiymnA+AcVl6hV9XNSd6c5Nru7vlGAmAVKwW9qq5P8u4kP9Pd/zPvSACsYi+3LR5L8mCSa6rqZFXdkuRPkrw8yb1V9UhV/fl5nhOAXey6Qu/um86w+QPnYRYAJvCkKMAgBB1gEIIOMAhBBxiEoAMMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQQg6wCAEHWAQgg4wCEEHGMSuQa+qO6rqVFU9+oJt31NV91bV49sfX3l+xwRgN3tZod+Z5Pod225Lcl93X53kvu3XAOyjXYPe3Q8keWbH5huSfGj78w8l+aWZ5wJgSateQ391dz+VJNsfXzXfSACs4rz/ULSqNqpqs6o2t7a2zvfpAC5Zqwb961X1miTZ/njqbAd299HuXu/u9cViseLpANjNqkG/J8nN25/fnORv5hkHgFXt5bbFY0keTHJNVZ2sqluS/F6S66rq8STXbb8GYB+t7XZAd990ll3XzjwLABN4UhRgEIIOMAhBBxiEoAMMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQQg6wCAEHWAQgg4wCEEHGISgAwxC0AEGMSnoVfWbVfVYVT1aVceq6jvnGgyA5awc9Ko6mOQ3kqx39+uTXJbkxrkGA2A5Uy+5rCX5rqpaS3JFkienjwTAKlYOend/LcntSZ5I8lSSf+/uj+08rqo2qmqzqja3trZWnxSAc5pyyeWVSW5IclWS1yZ5aVW9dedx3X20u9e7e32xWKw+KQDnNOWSy88m+efu3urubyT5SJKfmmcsAJY1JehPJPmJqrqiqirJtUlOzDMWAMuacg39oSR3J3k4yWe33+voTHMBsKS1KV/c3e9L8r6ZZgFgAk+KAgxC0AEGIegAgxB0gEEIOsAgBB1gEIIOMAhBBxiEoAMMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQUwKelW9oqrurqrPV9WJqvrJuQYDYDlrE7/+j5P8bXf/clVdnuSKGWYCYAUrB72qvjvJTyf51STp7meTPDvPWAAsa8oll+9PspXkg1X1qap6f1W9dKa5AFjSlKCvJfnRJH/W3W9I8t9Jbtt5UFVtVNVmVW1ubW1NOB0A5zIl6CeTnOzuh7Zf353Tgf8W3X20u9e7e32xWEw4HQDnsnLQu/tfk3y1qq7Z3nRtks/NMhUAS5t6l8uvJ7lr+w6XLyd52/SRAFjFpKB39yNJ1meaBYAJPCkKMAhBBxiEoAMMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQQg6wCAEHWAQgg4wCEEHGISgAwxC0AEGIegAg5gc9Kq6rKo+VVUfnWMgAFYzxwr9HUlOzPA+AEwwKehVdSjJLyZ5/zzjALCqqSv0P0ryriTPn+2Aqtqoqs2q2tza2pp4OgDOZuWgV9Wbk5zq7uPnOq67j3b3enevLxaLVU8HwC6mrNDfmOQtVfWVJB9O8qaq+otZpgJgaSsHvbvf092HuvtIkhuTfLy73zrbZAAsxX3oAINYm+NNuvv+JPfP8V4ArMYKHWAQgg4wCEEHGISgAwxC0AEGIegAgxB0gEEIOsAgBB1gEIIOMAhBBxiEoAMMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQaxctCr6sqq+vuqOlFVj1XVO+YcDIDlrE342ueS/FZ3P1xVL09yvKru7e7PzTQbAEtYeYXe3U9198Pbn/9nkhNJDs41GADLmeUaelUdSfKGJA+dYd9GVW1W1ebW1tYcpwPgDCYHvapeluSvkryzu/9j5/7uPtrd6929vlgspp4OgLOYFPSq+o6cjvld3f2ReUYCYBVT7nKpJB9IcqK7/2C+kQBYxZQV+huT/EqSN1XVI9t/fmGmuQBY0sq3LXb3PySpGWcBYAJPigIMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQQg6wCAEHWAQgg4wCEEHGISgAwxC0AEGIegAgxB0gEFMCnpVXV9VX6iqL1XVbXMNBcDyVg56VV2W5E+T/HyS1yW5qapeN9dgACxnygr9x5J8qbu/3N3PJvlwkhvmGQuAZa1N+NqDSb76gtcnk/z4zoOqaiPJxvbL/6qqL0w4J9/qQJKn93uI3VT99n6PwIXne3Ne37eXg6YEvc6wrb9tQ/fRJEcnnIezqKrN7l7f7zlgJ9+b+2PKJZeTSa58wetDSZ6cNg4Aq5oS9E8mubqqrqqqy5PcmOSeecYCYFkrX3Lp7ueq6u1J/i7JZUnu6O7HZpuMvXApixcr35v7oLq/7bI3ABchT4oCDELQAQYh6ACDmHIfOhdQVf1QTj+JezCn7/d/Msk93X1iXwcDXjSs0C8CVfXunP7VCpXkn3L6ltFKcswvRQO+yV0uF4Gq+mKSH+7ub+zYfnmSx7r76v2ZDM6tqt7W3R/c7zkuFVboF4fnk7z2DNtfs70PXqx+d78HuJS4hn5xeGeS+6rq8fz/L0Q7nOQHk7x936aCJFX1mbPtSvLqCznLpc4ll4tEVb0kp39l8cGc/otyMsknu/t/93UwLnlV9fUkP5fk33buSvKP3X2mf11yHlihXyS6+/kkn9jvOeAMPprkZd39yM4dVXX/hR/n0mWFDjAIPxQFGISgAwxC0AEGIegAg/g/vnkPINOqXx4AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ovarian['fustat'].value_counts().sort_index().plot.bar(color='darkblue')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Density plot of age" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAD8CAYAAABdCyJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XlYVeXax/HvLZOzmeKIphaOOQVankyzXk2ttEET07STpc2p1ck6p7I86dE6WlpvR83KIac00wanMssmBBU1B5LIEjVFxZFAgfv9Yy97OYSCwGbtDffnuvbF3ms/a+3fhg0361nPWo+oKsYYY0xBlXE7gDHGGP9mhcQYY0yhWCExxhhTKFZIjDHGFIoVEmOMMYVihcQYY0yhWCExxhhTKFZIjDHGFIoVEmOMMYUS6HaA4lC9enVt0KCB2zGMMcavbNiw4ZCqhubVrlQUkgYNGhAbG+t2DGOM8Ssi8kt+2lnXljHGmEKxQmKMMaZQrJAYY4wpFCskxhhjCsUKiTHGmEKxQmKMMaZQrJAYY4wplFJxHokxvmzXrhS++iqJAwdOUalSMFdeWZt27WohIm5HMyZfrJAY45JNmw7wxBNfsmbNr396rnnzaowdew29e1/mQjJjLox1bRlTzFSVCRPW067dHH744RDjxl3Dzp33kJY2nL177+ftt29ABG655UPuvns56ekZbkc25rxEVd3O4HWRkZFql0gxviArS7n//tVMn76Fvn0bM3VqN6pWLfundmfOZPLPf37Piy9+R6dOYSxffjvlywe5kNiUZiKyQVUj82pneyTGFBNV5bHH1jB9+haefvpK5s+/OdciAhAUFMALL1zNe+/dyNdf7+W225Zy+nRmMSc2Jn+skBhTTCZOjOX11zcxcmQEL73UkTJl8j6YfuedzZg2rRsrV+5m+PA1xZDSmAtnB9uNKQbr1iXx1FNfcfvt4bz88rUXNCJryJCW7Nx5mFdeiaVjxzDuvLOZF5Mac+Fsj8QYL0tJSSMq6mMaNbqIt9/unq89kZzGjr2Gq6+uy/33r2bPnuNeSGlMwVkhMcbLHn98LQcOnGLevBupXDmkQNsICgpg1qweZGZm8cADn1EaBskY/+HVQiIi3UUkXkQSRGRULs+HiMgC5/loEWngLG8vInHObbOI3Jptnd0istV5zoZiGZ+2atVu3nnnB558sh0REbUKta1GjS7in//syCefJPL++/FFlNCYwvPa8F8RCQB+BLoCSUAM0F9Vt2dr8yDQSlXvF5Eo4FZV7Sci5YHTqpohIrWBzUAd5/FuIFJVD+U3iw3/NW5IT8+gefN3CAoKIC5uEGXLFv6QZGZmFhERszl6NJ2dO+8pkm0acy6+MPy3PZCgqomqehqYD/TO0aY3MNO5vwi4XkREVVNV9exZWGUB2483fmfKlE0kJh5jypTriuwPfkBAGf7972v55ZfjvPbahiLZpjGF5c1CUhfYk+1xkrMs1zZO4TgGVAMQkStFZBuwFbg/W2FRYJWIbBCRoV7Mb0yBJSenMmbMd/Ts2ZCuXRsU6bavv/4Sbr75Ul56KZrk5NQi3bYxBeHNQpLb0JScexbnbKOq0araAmgHPC0iZ8/culpVrwB6AA+JSKdcX1xkqIjEikhscnJywd6BMQU0evS3nDp1hldeudYr2x8/vhMnT57mlVdivLJ9Yy6ENwtJElAv2+MwYN+52ohIIFAFOJK9garuAE4BlzuP9zlfDwJL8HSh/YmqTlPVSFWNDA0NLfSbMSa/EhOPMnXqZoYNa02zZtW88hrNmlUjKqopb7wRx6FDtldi3OXNQhIDhItIQxEJBqKAZTnaLAMGO/f7AGtUVZ11AgFE5BKgCbBbRCqISCVneQWgG/CDF9+DMRfspZe+JzCwDH//+1VefZ1//OMqUlPPMHGiHSsx7vJaIXGOaTwMrAR2AAtVdZuIvCgivZxmM4BqIpIAjATODhHuCGwWkTg8ex0POqO0agJfi8hmYD3wiaqu8NZ7MOZCJSYeZdas7Qwb1po6dSp69bWaN69O375NmDJlIykpaV59LWPOx67+a0wRuvfelcyZs53ExPu8XkgA4uIO0rbtLCZM6MSTT+bay2tMgfnC8F9jSpWffz7KzJnbimVv5Kw2bWrQpUs9Jk/exJkzdnVg4w4rJMYUkVdf3YgI/O1v7Yr1dUeMiCQp6QSLF+8q1tc15iwrJMYUgePH03nnnR/o168pdetWKtbXvvHGRoSHV2XixFi7BpdxhRUSY4rAzJnbOHHiNI8+2rbYX7tMGWH48CuIifmN77/fX+yvb4wVEmMKKStLmTJlEx061KFdu9quZLjrrhZUrBjE1KmbXXl9U7pZITGmkFas+Jldu1J49NErXMtQqVIwAwY0Z8GCeBsKbIqdFRJjCmny5I3UqVOR228PdzXHsGGtSEvLYPbs7Xk3NqYIWSExphB27jzMypW7eeCB1gQFBbiapW3bmrRvX4upUzfbQXdTrKyQGFMIr7++ieDgAIYObeV2FACGDWvN9u2H+eabvW5HMaWIFRJjCujYsXTefXcb/fs3pUaNCm7HAaBfvyZUrhzMtGlb3I5iShErJMYU0Dvv/MCpU2dcPcieU4UKwURFNWXx4h85ceK023FMKWGFxJgCyMzMYsqUjXTsWJcrrqjpdpz/MnhwC1JTM1i8+Ee3o5hSwgqJMQXw6ac/k5h4zKf2Rs7q0KEO4eFVmTlzm9tRTClhhcSYApg8eSNhYZW45ZbL3I7yJyLCoEHNWbt2D7t3H3M7jikFrJAYc4G2bz/EZ5/9woMPtnF9yO+53HVXcwA7p8QUCyskxlygKVM2ERISwH33tXQ7yjldckkVunSpx6xZ2+ycEuN1VkiMuQApKWnMmrWNAQOaUb16ebfjnNfgwS1ISDjKt9/uczuKKeG8WkhEpLuIxItIgoiMyuX5EBFZ4DwfLSINnOXtRSTOuW0WkVvzu01jvOntt7eSmprhkwfZc7r99sZUqBBkB92N13mtkIhIAPAG0ANoDvQXkeY5mg0BUlT1MmASMN5Z/gMQqaptgO7AVBEJzOc2jfGKzMwsXn99E507h9G6dQ234+SpYsVgbr01nEWLfuT0aZs90XiPN/dI2gMJqpqoqqeB+UDvHG16AzOd+4uA60VEVDVVVTOc5WWBs528+dmmMV7x0Uc/sXv3cb/YGzmrf/+mpKSksXLlbrejmBLMm4WkLrAn2+MkZ1mubZzCcQyoBiAiV4rINmArcL/zfH62aYxXTJ68kfr1K9Grl+8N+T2Xrl0voVq1csybt8PtKKYE82YhkVyW5Rw+cs42qhqtqi2AdsDTIlI2n9v0bFhkqIjEikhscnLyBcQ25s+2bk3miy/28NBDbQkM9J8xKkFBAfTp05ilSxM4dcoumWK8w5u/EUlAvWyPw4Ccw0f+aCMigUAV4Ej2Bqq6AzgFXJ7PbZ5db5qqRqpqZGhoaCHehjGeIb/lygVy772+O+T3XPr3b0pqagbLlv3kdhRTQnmzkMQA4SLSUESCgShgWY42y4DBzv0+wBpVVWedQAARuQRoAuzO5zaNKVKHD//OnDnbGTiwORdfXM7tOBfsmmvCqFu3IvPm7XQ7iimhvFZInGMaDwMrgR3AQlXdJiIvikgvp9kMoJqIJAAjgbPDeTsCm0UkDlgCPKiqh861TW+9B2MAZszYyu+/Z/DII23djlIgZcoI/fo1YcWKnzly5He345gSSErDWa+RkZEaGxvrdgzjhzIysrj00ulceulFrFnTz+04BRYb+xvt2s1h+vRu3Huvb0zCZXyfiGxQ1ci82vnPUUNjXLB0aQK//nqCxx6LcDtKoURE1CQ8vKp1bxmvsEJizHlMmrSBhg2rcNNNjdyOUigiQv/+Tfnii1/Zv/+k23FMCWOFxJhziInZzzff7OXRR68gIMD/f1X692+KKixYEO92FFPC+P9vhzFe8uqrG6lUKZh77rnc7ShFomnTarRpU4P58617yxQtKyTG5GLv3hMsXBjPkCEtqVw5xO04RSYqqgnR0fv5+eejbkcxJYgVEmNy8cYbcWRlKY8+6p9Dfs+lX7+mgHVvmaJlhcSYHFJTzzB16mZ6976Mhg0vcjtOkWrQoAodOtSx0VumSFkhMSaH2bO3c+RIGiNG+PeQ33OJimrKli3JbN9+yO0opoSwQmJMNpmZWfz737FERNSkY8eSeWHpvn0bU6aMWPeWKTJWSIzJ5oMPdrFrVwpPPdUekdwuNu3/ateuyLXX1mP+/J02n7spElZIjHGoKuPGRdO4cVVuuy3c7TheFRXVlB9/TGHTpoNuRzElgBUSYxyrVu1m06aDPPVU+xJxAuL53HZbOIGBZeycElMkSvZvizEXYNy4aMLCKjFwYHO3o3hdtWrl6NbtEhYs2ElWlnVvmcKxQmIM8O23e/nyyyQefzyS4OAAt+MUi/79m/Hrryf4/vtc54YzJt+skBgDjB79LdWrl+O++/xvBsSC6tXrUsqWDbRzSkyhWSExpd6XX+5h9epfePrpK6lQIdjtOMWmcuUQbryxIQsXxpORkeV2HOPHrJCYUk1VefbZb6hduwIPPNDa7TjFLiqqKQcPpvLll3vcjmL8mBUSU6qtXv0L69Yl8Y9/XEW5ckFuxyl2N97YiIoVg2z0likUrxYSEekuIvEikiAio3J5PkREFjjPR4tIA2d5VxHZICJbna/XZVtnrbPNOOdWw5vvwZRcWVnK3/++jksuqcyQIaXn2Eh25coFccst4SxevIvTpzPdjmP8lNcKiYgEAG8APYDmQH8RyTmucgiQoqqXAZOA8c7yQ8DNqtoSGAzMzrHeAFVt49zsjCpTIHPn7iA29gAvvPAXQkIC3Y7jmqioJqSkpLFq1W63oxg/5c09kvZAgqomquppYD7QO0eb3sBM5/4i4HoREVXdpKpnxyRuA8qKSMmZFMK4LjX1DE8/vY6IiJrcdVcLt+O4qmvXBlStWta6t0yBebOQ1AWyH8FLcpbl2kZVM4BjQLUcbW4HNqlqerZl7zjdWs9KSb0gkvGqV16JISnpBJMmdaFMmdL9EQoODqBPn8YsXZpAauoZt+MYP+TNQpLbb2fOU2jP20ZEWuDp7hqW7fkBTpfXNc7trlxfXGSoiMSKSGxycvIFBTclW1LSCcaPX0/fvo255powt+P4hKioppw8eYZPP010O4rxQ94sJElAvWyPw4Ccp9D+0UZEAoEqwBHncRiwBBikqj+dXUFV9zpfTwBz8XSh/YmqTlPVSFWNDA0NLZI3ZEqGRx9dQ1YWjB/fye0oPqNz5zBq1ixvJyeaAvFmIYkBwkWkoYgEA1HAshxtluE5mA7QB1ijqioiFwGfAE+r6jdnG4tIoIhUd+4HATcBP3jxPZgSZsmSXSxZsovRozuUuNkPCyMgoAx33NGETz5J5Pjx9LxXMCYbrxUS55jHw8BKYAewUFW3iciLItLLaTYDqCYiCcBI4OwQ4YeBy4BncwzzDQFWisgWIA7YC0z31nswJcuxY+k89NBntG4dysiRkW7H8Tn9+zcjPT2TpUsT3I5i/IyUholtIiMjNTY21u0YxmXDhq3irbe2Eh09gMjIWm7H8TmqSsOG02nRohqffHK723GMDxCRDaqa539ddma7KRWWLUtg2rQtjBwZYUXkHESEfv2asGrVLxw+/LvbcYwfsUJiSrx9+05yzz0radu2Bv/8Z0e34/i0qKimZGRksXjxj25HMX7ECokp0bKylLvvXk5q6hnmzr2xVJ/Bnh9t2tSgceOqdnKiuSBWSEyJ9sIL37J69S+8+moXmjbNea6ryUlE6N+/KWvX7mH//pNuxzF+wgqJKbEWL/6RF1/8jr/+9XLuu6+V23H8RlRUU1SxvRKTb1ZITIm0dWsygwcv56qravPmm/+DXUkn/5o2rUa7drWYOXOb21GMn7BCYkqcX389Ts+eH1C5cjCLF/e24yIFcPfdLdi8OZm4OLu4tsmbFRJTohw6lMoNNyzixInTrFjRhzp1KrodyS9FRTUlODiAd9+1C0eYvFkhMSXG0aNp3HjjB/z88zGWLbuFVq3sGmsFdfHF5ejV61Lee2+HTXhl8mSFxJQIBw+eokuXhWzadJCFC2+mU6d6ea9kzuvuu1tw6NDvLF/+s9tRjI+zQmL83q+/Hueaa+YTH3+Ejz66lV69LnM7Uolwww0NqVmzvHVvmTxZITF+LT7+CB07zuPAgVRWr+7LDTc0dDtSiREYWIaBA5vz8ceJJCenuh3H+LB8FRIRWSwiN4qIFR7jMzZuPEDHjvNIT89k7dp+XH11zgk4TWENHtyCjIwsm6fEnFd+C8ObwJ3ALhH5l4g09WImY/K0bl0SXbosoHz5QNati6JNmxpuRyqRWrYMJSKiJjNmbKU0XCncFEy+ComqfqaqA4ArgN3AahH5VkT+6kwwZUyx+fTTRLp1W0Tt2hX5+uv+NG58sduRSrT77mvFli3JREfvdzuK8VH57qoSkWrA3cC9wCbgNTyFZbVXkhmTiwULdtK794c0b16NdeuiqFevstuRSrw772xGxYpBTJ262e0oxkfl9xjJB8A6oDxws6r2UtUFqvoIYGd8mWIxdepm+vf/mA4d6rBmzR2EhpZ3O1KpUKlSMAMGNGfBgnhSUtLcjmN8UH73SN5S1eaqOk5V9wOISAhAfmbPMqaw/v3vGO6/fzU9ejRkxYrbqVIlxO1IpcqwYa34/fcMZs/e7nYU44PyW0j+mcuy7/JaSUS6i0i8iCSIyKhcng8RkQXO89Ei0sBZ3lVENojIVufrddnWiXCWJ4jIZLGr8ZV4//zndzzxxJf07duYJUtuoXx5OyxX3Nq2rUm7drWYOnWzHXQ3f3LeQiIitUQkAignIm1F5Arndi2ebq7zrRsAvAH0AJoD/UWkeY5mQ4AUVb0MmASMd5YfwtOF1hIYDMzOts6bwFAg3Ll1z/ttGn+kqvz97+t49tlvuOuu5sydexPBwQFuxyq1hg1rzfbth/nmm71uRzE+Jq89khuAV4AwYCLwb+c2Engmj3XbAwmqmqiqp4H5QO8cbXoDM537i4DrRURUdZOq7nOWbwPKOnsvtYHKqvqdev4tmgXckue7NH7pmWfWMXZsNPfe25J33+1BYKCdxuSmqKgmVK4czH/+YwfdzX8772+mqs5U1S7A3araJdutl6p+kMe26wJ7sj1Ocpbl2kZVM4BjQM5p7G4HNqlqutM+KY9tAiAiQ0UkVkRik5OT84hqfM348dH861/rGTasNVOndqNMGevBdFuFCsEMHtyChQvj+e23U27HMT4kr66tgc7dBiIyMuctj23n9pufs3P1vG1EpAWe7q5hF7BNz0LVaaoaqaqRoaF2FVh/Mn36FkaNWkdUVFPeeON6KyI+5JFHriAjI4s334xzO4rxIXn1FVRwvlYEKuVyO58kIPslWMOAfedqIyKBQBXgiPM4DFgCDFLVn7K1D8tjm8aPLV2awLBhq+jRoyEzZ/YgIMC6s3xJeHhVbrrpUt58M460tAy34xgfcd6p41R1qvP1hQJsOwYIF5GGwF4gCs9lVrJbhudg+ndAH2CNqqqIXAR8Ajytqt9ky7NfRE6IyFVANDAImFKAbMYHbd58kAEDPiEyshaLFvWyA+s+avjwCK6//ifmzdvBX//a0u04xgfk94TECSJSWUSCRORzETmUrdsrV84xj4eBlcAOYKGqbhORF0Wkl9NsBlBNRBLwHMA/O0T4YeAy4FkRiXNuZy+m9ADwFpAA/AQsz//bNb7qwIFT9Oq1hIsuCuHDD22Iry/r0qUeLVtW59VXN9pQYAOA5OeDICJxqtpGRG7FM0pqBPCFqrb2dsCiEBkZqbGxsW7HMOeQkZHFddctJDb2N9atiyIiopbbkUwe3n57K0OGrOTzz+/guuvqux3HeImIbMjPSef57YA+++9hT2Ceqh4pcDJjchg9+lvWrUti+vRuVkT8xJ13NiM0tBwvv7ze7SjGB+S3kHwkIjuBSOBzEQkF7KI7ptBWr97N2LHfM2RISwYMyHm+qvFVZcsGMmJEJCtW7GbjxgNuxzEuy+9l5EcBHYBIVT0DnOLPJxcac0EOHDjFwIGf0rx5NSZPvi7vFYxPefDBNlSpEsLYsd+7HcW47LyjtnJohud8kuzrzCriPKaUUFUeeOAzjh1LZ82aO+zguh+qUiWEhx9uy9ix37Njx2GaNct5LrEpLfI7ams2nkuldATaOTe76q8psPnzd7JkyS7GjLmaFi2qux3HFNDw4VdQrlwg48fbsZLSLL97JJFAc7WxfqYI/PbbKR5++HOuuqo2I0fa/yP+rHr18gwd2popUzby/PMdaNjwIrcjGRfk92D7D4ANpzFF4pFHPufUqTO88053O3O9BHjiiUgCA8vw4ot5zixhSqj8/hZXB7aLyEoRWXb25s1gpmRatWo3ixb9yD/+cRVNm1qfeklQt24lHn64LbNmbWf79kNuxzEuyO8JiZ1zW66qXxZ5Ii+wExJ9Q3p6Bi1bemYN2Lp1MCEhFzLWw/iyQ4dSadToLbp2vYTFi21AZ0lRpCckOgVjNxDk3I8BNhYqoSl1Xnklll27Upgy5TorIiVM9erlefzxSD74YBcxMfvdjmOKWX5Hbd2HZ+Kpqc6iusCH3gplSp7du4/x0kvfc/vt4dxwQ0O34xgvGDkykurVy/HMM1+7HcUUs/weI3kIuBo4DqCqu4Aa513DmGyeeWYdAJMmdXE5ifGWSpWC+fvfr+Kzz37h008T3Y5jilF+C0m6M10u8MfcITYU2ORLbOxvzJu3k5EjI6lXr7LbcYwXPfhgGxo3rsqIEV9w+nSm23FMMclvIflSRJ4ByolIV+B94CPvxTIlhary5JNfUr16Of72t3ZuxzFeFhwcwKRJXfjxxxRef32T23FMMclvIRkFJANb8Ux7+ynwD2+FMiXHp58msnbtHp5/vgOVK4e4HccUg549G9GjR0NeeOFbDh60ud1Lg/yO2srCc3D9QVXto6rT7Sx3k5fMzCyeeuorwsOrMmyYX0xdY4rIpEldSE3NsAPvpcR5C4l4jBaRQ8BOIF5EkkXkueKJZ/zZvHk72bbtMC+91JGgIJs2tzRp0uRiRoyIYMaMrXz55R634xgvy2uPZDie0VrtVLWaql4MXAlcLSIj8tq4iHQXkXgRSRCRUbk8HyIiC5zno0WkgbO8moh8ISInReT1HOusdbaZcwpe40MyMrJ48cXvaNUqlNtvb+x2HOOC55/vQIMGlRk2bDVpaRluxzFelFchGQT0V9Wfzy5Q1URgoPPcOYlIAPAG0ANoDvQXkZwzFw0BUlT1MmASMN5ZngY8Czxxjs0PUNU2zu1gHu/BuGDu3B3s2pXC6NF/oUwZcTuOcUGFCsFMndqN+PgjNmdJCZdXIQlS1T9dPEdVk/n/6XfPpT2QoKqJztDh+fx5MqzewEzn/iLgehERVT2lql9jszD6pYyMLMaM+Y42bWpwyy2XuR3HuKhbtwYMHNicf/1rPT/8kOx2HOMleRWS0wV8Djxnv2fvHE1yluXaRlUzgGNAfq7k947TrfWsiNi/uz5mzpztJCQcZfTov2A/HjNx4rVUqRLC3Xev4MwZO7ekJMqrkLQWkeO53E4ALfNYN7e/IDlHeuWnTU4DVLUlcI1zuyvXFxcZKiKxIhKbnGz/CRWXM2cyGTPmO664oia9el3qdhzjA0JDyzN1alc2bDjAmDF2qfmS6LyFRFUDVLVyLrdKqppX11YSUC/b4zBg37naOGfLVwGO5JFpr/P1BDAXTxdabu2mqWqkqkaGhobmEdUUlffe20Fi4jHbGzH/5bbbGjN4cAteeima77/P+WfA+DtvzioUA4SLSEMRCQaigJxzmCwDBjv3+wBrznd+iogEikh1534QcBOeSbeMD8jKUiZMWE/r1qHcdFMjt+MYH/Paa9cRFlaRu+76lFOn8uoZN/7Ea4XEOebxMLAS2AEsVNVtIvKiiPRyms0AqolIAjASzxn0AIjIbmAicLeIJDkjvkKAlSKyBYgD9gLTvfUezIVZvjyRHTuO8OST7WxvxPxJlSohzJrVk59+Osrw4V+4HccUoXxNbOXvbGKr4nHttfNJTDzGTz/daycgmnN65pl1jBsXzZw5PRkwIOcZAcaXFOnEVsbkJSZmP19+mcTw4RFWRMx5vfji1XTsWJdhw1YTH3/eQ6LGT1ghMUXilVdiqVIlhPvua+V2FOPjAgPLMG/eTZQrF0jfvsv4/fczbkcyhWSFxBRaYuJRFi36kWHDWlGpUrDbcYwfCAurxOzZPdm69RCPPWbHS/ydFRJTaK++uoGAAOGxxyLcjmL8SPfuDRk1qj3Tp29h1qxtbscxhWCFxBTK4cO/M2PGVgYMaEadOhXdjmP8zJgxHenSpR7Dhq1m06YDbscxBWSFxBTKm2/GkZqawRNP2OyH5sIFBpZh/vybqF69HLfdtpTDh393O5IpACskpsDS0jKYMmUTPXo0pEWL6m7HMX6qRo0KLF7ci337TnHnnZ+QmZnldiRzgayQmAKbPXs7Bw+m8uSTtjdiCqd9+9q8/vr1rFq1m+ee+8btOOYCWSExBZKVpfz737FERNTk2mvr5b2CMXm4775W3HtvS8aOjebDD3e5HcdcACskpkA++ugn4uOP8MQTdjkUU3SmTLmedu1qMWjQcjtZ0Y9YITEF8sorMTRoUJk+fWwaXVN0ypYNZPHiXpQtG8Ctt37IiRN2cUd/YIXEXLDvv9/H11/vZcSISAID7SNkila9epVZsOBm4uNT+Otfl1Margfo7+yvgLlgL78cQ9WqZbnnnsvdjmJKqC5d6jNhQicWL97FK6/EuB3H5MEKibkgCQkpLFmyiwceaE3FinY5FOM9I0dGcscdTRg1ah1ffPGr23HMeVghMRdk4sRYgoICeOSRK9yOYko4EWHGjBto3Lgq/ft/zG+/nXI7kjkHKyQm35KTU3nnnW0MGtScWrUquB3HlAIVKwazaFEvjh8/Tf/+H5ORYScr+iIrJCbf/vd/40hLy2DkyDznuTGmyLRoUZ3//Kcra9fu4fnn7WRFX2SFxORLauoZXn99EzfffCnNmlVzO44pZQYNavHqdzwZAAAUbUlEQVTHyYrLlye6Hcfk4NVCIiLdRSReRBJEZFQuz4eIyALn+WgRaeAsryYiX4jISRF5Pcc6ESKy1VlnstjZcMVi5sxtHDr0u10Oxbhm8uTraN06lIEDP+XXX4+7Hcdk47VCIiIBwBtAD6A50F9Eck7QPARIUdXLgEnAeGd5GvAs8EQum34TGAqEO7fuRZ/eZJeZmcXEibG0b1+Ljh3ruh3HlFLlygXx/vu9OHMmizvu+IjTpzPdjmQc3twjaQ8kqGqiqp4G5gO9c7TpDcx07i8CrhcRUdVTqvo1noLyBxGpDVRW1e/Uc5bSLOAWL74HAyxdmkBCwlGefNIuh2LcFR5elbffvoHo6P08++zXbscxDm8WkrrAnmyPk5xlubZR1QzgGHC+Dvi6znbOt01ThFSVl1+OoVGjKtx6a7jbcYyhT58mDBvWmgkTYvj881/cjmPwbiHJ7V/XnNc6yE+bArUXkaEiEisiscnJyefZpDmfb77Zy/ff72fkyEgCAmxshvENEydeS9OmFzNo0HIOHUp1O06p582/DElA9uuLhwH7ztVGRAKBKsD5LvmZ5GznfNsEQFWnqWqkqkaGhoZeYHRz1vjx66levRx//atdDsX4jvLlg5g790aSk1O5995Vdj0ul3mzkMQA4SLSUESCgShgWY42y4DBzv0+wBo9zydCVfcDJ0TkKme01iBgadFHNwDbth3i448TeeSRtpQvH+R2HGP+S9u2NfnXvzqxdGkC06ZtcTtOqea1QuIc83gYWAnsABaq6jYReVFEejnNZgDVRCQBGAn8MURYRHYDE4G7RSQp24ivB4C3gATgJ2C5t95DaffyyzGULx/IQw+1dTuKMbkaPjyCbt0aMGLEF+zYcdjtOKWWlIZdwsjISI2NjXU7hl9JSjpBw4bTefDBNrz22nVuxzHmnPbvP0mrVjOpW7ci0dEDCAkJdDtSiSEiG1Q1z0tZ2NFTk6tJk2JRVUaOjHA7ijHnVbt2Rd5+uzubNyczevS3bscplayQmD9JSUlj2rQtREU15ZJLqrgdx5g83XzzpQwZ0pIJE2L49tu9bscpdayQmD958804Tp48w9/+1t7tKMbk28SJ11KvXiUGD17OqVM2RW9xskJi/ktaWgavvbaR7t0b0KqVDZs2/qNy5RDefbc7CQlHeeqpr9yOU6pYITH/5Z13fuDgwVSeesr2Roz/ufba+gwfHsEbb8Tx2Wd21ntxsUJi/pCensG4cdFcfXVdOneul/cKxvigsWM70qTJxdxzzwqOHUt3O06pYIXE/GHmzG3s2XOC557rYBdnNH6rXLkgZs3qwb59J3nssTVuxykVrJAYAM6cyWTs2GiuvLI2Xbte4nYcYwqlffvaPP30lcycuY2lSxPcjlPiWSExAMyevZ1ffjlueyOmxHj22Q60bVuDoUNXkZxsF3b0JiskhoyMLF566XsiImrSo0dDt+MYUySCgwOYNasHR4+mc//9q+3Cjl5khcQwd+4OEhOP2d6IKXEuvzyUMWOu5oMPdjFnzna345RYVkhKudOnM3nhhW9p06YGN998qdtxjClyjz8eSceOdXn44c9trncvsUJSyk2fvoXExGOMG3eN7Y2YEikgoAwzZ/YgK0u5++7lZGVZF1dRs0JSip08eZoxY76jc+cwbrihgdtxjPGaRo0uYtKkLnzxxR6mTNnodpwSxwpJKfbqqxs4cCCVceM62d6IKfGGDGnJzTdfyqhR62zukiJmhaSUOnQolZdfjqF378vo0KGO23GM8ToRYfr0blSsGMRdd33KmTOZbkcqMayQlFJjx0Zz8uQZXnqpo9tRjCk2NWtWYNq0bmzYcIAxY75zO06JYYWkFIqPP8KUKZu4++4WtGhR3e04xhSrW28NZ/DgFowdG0109H6345QIXi0kItJdROJFJEFERuXyfIiILHCejxaRBtmee9pZHi8iN2RbvltEtopInIjY/LkFMGLEF5QvH8jYsde4HcUYV7z22nXUrVuRu+761OYuKQJeKyQiEgC8AfQAmgP9RaR5jmZDgBRVvQyYBIx31m0ORAEtgO7A/zrbO6uLqrbJz1zC5r998slPLF/+M88914GaNSu4HccYV1SpEsLMmT1ISEjhsce+cDuO3/PmHkl7IEFVE1X1NDAf6J2jTW9gpnN/EXC9eIYP9Qbmq2q6qv4MJDjbM4Vw+nQmI0aspUmTi3nkkSvcjmOMq669tj7PPHMVM2Zs5b337Kz3wvBmIakL7Mn2OMlZlmsbVc0AjgHV8lhXgVUiskFEhp7rxUVkqIjEikhscnJyod5ISTFxYiy7dqUwadK1BAcH5L2CMSXc6NF/oWPHutx//2p+/PGI23H8ljcLSW4nJuQ8pfRcbc637tWqegWeLrOHRKRTbi+uqtNUNVJVI0NDbcrYXbtSeOGF77jttnB69GjkdhxjfEJgYBnmzbuJkJBA+vX7mLS0DLcj+SVvFpIkIPs0e2HAvnO1EZFAoApw5HzrqurZrweBJViXV55UlaFDVxESEsCUKde7HccYnxIWVol33+1OXNxBnnhirdtx/JI3C0kMEC4iDUUkGM/B82U52iwDBjv3+wBr1HOt52VAlDOqqyEQDqwXkQoiUglARCoA3YAfvPgeSoS33/6BtWv3MGFCJ+rUqeh2HGN8zk03XcrIkZ653ufN2+F2HL8T6K0Nq2qGiDwMrAQCgLdVdZuIvAjEquoyYAYwW0QS8OyJRDnrbhORhcB2IAN4SFUzRaQmsMS5nEcgMFdVV3jrPZQEv/xyjMcfX0unTmHce28rt+MY47PGjevE+vW/MWTISpo2vZi2bWu6HclvSGmY7CUyMlJjY0vfKSeZmVl06bKQuLiDxMUNolGji9yOZIxPO3DgFJGRcyhTBmJiBlKjRukeIi8iG/JzmoWd2V6CTZgQw7p1SUyZcp0VEWPyoWbNCnz4YW8OHvydvn0/sutx5ZMVkhJq/fr9PPfcN9xxRxMGDWrhdhxj/EZERC1mzLiBr75K4sEHP7MpevPBa8dIjHsOHDjFbbctJSysIv/5T1e7RLwxF+jOO5uxffthXnrpe8LCKvH8839xO5JPs0JSwpw5k8kdd3zEkSNpfPvtnVStWtbtSMb4pTFjrmbv3hOMHv0ttWtXYOjQ1m5H8llWSEoQVWXkyLV89VUSc+b0pE2bGm5HMsZviQjTpnXj4MFUHnjgM0JDy3PrreFux/JJdoykBHn55Rhef30Tjz8eyYABOa+PaYy5UEFBASxceDPt29eiX7+PWLo0we1IPskKSQkxa9Y2nnrqK6KimjJhQme34xhTYlSoEMyKFX1o27YGffos48MPd7kdyedYISkB3n8/niFDVnL99fV5993ulCljB9eNKUpVqoSwalVfIiJq0rfvR8yfv9PtSD7FComfmzNnO1FRH9O+fS0++KA3ISF22MsYb/AUkz785S916N//YyZMWG9Dgx1WSPzYf/4Tx6BBn9K5cxgrV/ahcuUQtyMZU6JVrhzCypV96NevCU899RUPPfSZnbSIjdrySxkZWTz++FomT95Iz54NWbSoF+XKBbkdy5hSoWzZQObOvYn69Svz8ssxbNlyiPnzbyIsrJLb0VxjeyR+5uDBU9x44wdMnryR4cMjWLr0VisixhSzMmWECRM6M3fujcTFHaRt21msWPGz27FcY4XEjyxfnkjLljP58ss9TJvWjUmTuhAYaD9CY9zSv38zYmMHUqtWBXr0WMyQIStISUlzO1axs79CfiA5OZUhQ1bQs+cH1KhRnpiYgdx3n10S3hhf0LRpNdavH8CoUe2ZOXMbzZu/w9y5O8jKKj0H4q2Q+LC0tAwmT95I48YzmDVrO3/7WztiYgbSsqVNHWyMLylXLohx4zoREzOQunUrMmDAJ7RrN4fPPvvF7WjFwgqJDzp58jRTpmzk0kvf4rHH1hAZWYutWwczfnxnypa18RHG+Kq2bWuyfv1AZs/uyeHDv9O16/t06PAe778fT0ZGltvxvMYmtvIRqkpc3EGmT9/CnDk7OHHiNJ06hfH883+hS5d6dgVfY/xMWloGM2ZsZdKkDfz001Hq16/EXXe1YODAZjRtWs3tePmS34mtvFpIRKQ78BqeqXbfUtV/5Xg+BJgFRACHgX6qutt57mlgCJAJPKqqK/Ozzdz4aiFJS8tg/fr9fPTRTyxZksBPPx0lJCSAO+5owrBhrbn66rpuRzTGFFJmZhYff5zIm2/GsXr1L2RlKa1bh9KzZyO6d29Ahw51CAoKcDtmrlwvJCISAPwIdAWSgBigv6puz9bmQaCVqt4vIlHAraraT0SaA/OA9kAd4DOgsbPaebeZG18oJKdPZxIff4StWw+xZUsy3323j+jo/aSnZxIUVIbrrqvPrbeG07dvYy6+uJyrWY0x3vHbb6eYN28HS5cm8M03+8jIyKJcuUAiImpy5ZW1adeuFs2bV+Oyyy7yiWH9vlBIOgCjVfUG5/HTAKo6LlublU6b70QkEPgNCAVGZW97tp2z2nm3mZuiLiSqypkzWaSnZ/L772c4ejSdI0fSSElJIyUlnZSUNPbvP8WePcf59dcT7Nlzgl9+Of5HH2lgYBlatw6lc+cwOneuR+fO9ahSxc5KN6Y0OXYsnc8//4Wvvkpi/frf2LjxAOnpnrPkRaB+/cpceulF1K5dgTp1KlK7dgVq167IxReXpXLlYOcWQuXKwVSsGOyVa+zlt5B488htXWBPtsdJwJXnaqOqGSJyDKjmLP8+x7pn+3ny2maRufnmD9i58wjp6Zl/3NLSMv74YZ9PmTJC3boVqV+/Mu3a1eKOO5rQsmV1Lr+8Ok2aXExwsG/uyhpjikeVKiHcdltjbrvN09ly+nQm27cfZufOI/z44xHi41P4+edjfPPNXvbvP5Xn353AwDKEhAT8cQsO9nyNixvk9b0bbxaS3Mpjzt2fc7U51/LcRpnlukslIkOBoQD169c/d8rzCA+vSqVKwc4PJvCPH1DZsv//uGzZAC66qCxVq4ZQtWrZP27Vq5ezkwWNMfkWHBxAmzY1cp2QTlU5ejSd/ftPcvRoOsePn+bYMc/X48fTOXnyjPPPbsZ//eN7+nRmsRx/8WYhSQLqZXscBuw7R5skp2urCnAkj3Xz2iYAqjoNmAaerq2CvIGJE7sUZDVjjClSIvLHP6m+yJv/MscA4SLSUESCgShgWY42y4DBzv0+wBr1HLRZBkSJSIiINATCgfX53KYxxphi5LU9EueYx8PASjxDdd9W1W0i8iIQq6rLgBnAbBFJwLMnEuWsu01EFgLbgQzgIVXNBMhtm956D8YYY/JmJyQaY4zJVX5HbdnRYGOMMYVihcQYY0yhWCExxhhTKFZIjDHGFIoVEmOMMYVSKkZtiUgykJ8ZZqoDh7wcp6hY1qLnLznBsnqLv2QtrpyXqGqeM+mVikKSXyISm5+hbr7AshY9f8kJltVb/CWrr+W0ri1jjDGFYoXEGGNMoVgh+W/T3A5wASxr0fOXnGBZvcVfsvpUTjtGYowxplBsj8QYY0yhlNpCIiL1ROQLEdkhIttE5DFn+cUislpEdjlfq7qcs6yIrBeRzU7OF5zlDUUk2sm5wLmsvk8QkQAR2SQiHzuPfTKriOwWka0iEicisc4yn/r5nyUiF4nIIhHZ6XxmO/haVhFp4nwvz96Oi8hwX8t5loiMcH6nfhCRec7vmq9+Vh9zcm4TkeHOMp/5vpbaQoLn8vSPq2oz4CrgIRFpjme++M9VNRz43HnspnTgOlVtDbQBuovIVcB4YJKTMwUY4mLGnB4DdmR77MtZu6hqm2xDKX3t53/Wa8AKVW0KtMbz/fWprKoa73wv2wARQCqwBB/LCSAidYFHgUhVvRzPtBRR+OBnVUQuB+4D2uP52d8kIuH40vdVVe3mOU60FOgKxAO1nWW1gXi3s2XLWB7YiGee+kNAoLO8A7DS7XxOljA8H+rrgI/xTJvsq1l3A9VzLPO5nz9QGfgZ55imL2fNlq0b8I2v5gTqAnuAi/HMy/QxcIMvflaBvsBb2R4/C/zNl76vpXmP5A8i0gBoC0QDNVV1P4Dz9c8TKBczp6soDjgIrAZ+Ao6qaobTJAnPL4YveBXPhzzLeVwN382qwCoR2SAiQ51lPvfzBxoBycA7TpfhWyJSAd/MelYUMM+573M5VXUv8ArwK7AfOAZswDc/qz8AnUSkmoiUB3rimXLcZ76vpb6QiEhFYDEwXFWPu50nN6qaqZ7ugjA8u7fNcmtWvKn+TERuAg6q6obsi3Np6npWx9WqegXQA0/XZie3A51DIHAF8KaqtgVO4QPdQ+fiHFfoBbzvdpZzcY4n9AYaAnWACng+Bzm5/llV1R14utxWAyuAzXi65n1GqS4kIhKEp4i8p6ofOIsPiEht5/naePYCfIKqHgXW4jmmc5GInJ0qOQzY51aubK4GeonIbmA+nu6tV/HNrKjqPufrQTx9+e3xzZ9/EpCkqtHO40V4CosvZgXPH+SNqnrAeeyLOf8H+FlVk1X1DPAB8Bd897M6Q1WvUNVOeKYl34UPfV9LbSEREcEzZ/wOVZ2Y7allwGDn/mA8x05cIyKhInKRc78cnl+AHcAXQB+nmes5AVT1aVUNU9UGeLo21qjqAHwwq4hUEJFKZ+/j6dP/AR/7+QOo6m/AHhFp4iy6HtiOD2Z19Of/u7XAN3P+ClwlIuWdvwVnv6c+91kFEJEaztf6wG14vr++8311+0CSWzegI57d1i1AnHPriadP/3M8Ff9z4GKXc7YCNjk5fwCec5Y3AtYDCXi6EELc/p7myH0t8LGvZnUybXZu24C/O8t96uefLW8bINb5HHwIVPXFrHgGhBwGqmRb5nM5nVwvADud36vZQIgvfladrOvwFLrNwPW+9n21M9uNMcYUSqnt2jLGGFM0rJAYY4wpFCskxhhjCsUKiTHGmEKxQmKMMaZQrJAYY4wpFCskxhhjCsUKiTHGmEL5P1o+56TGYxY8AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ovarian['age'].plot.density(color='darkblue')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bar plot of [ECOG Performance Status](https://training.seer.cancer.gov/followup/procedures/dataset/ecog.html)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAACrFJREFUeJzt3VGMZndZx/Hf445ESzFg+lK122WrITUEjTUTFUnUUBtXbagXXrQJUrHJXAlorFhCInfGxEYl0agbWEu0KSYVA+FCadCmMdbG2bbQlgUhWMtCcaepUasXdePjBQNZht2ded9zdof97+eTNDPnvOd9/8/F9JuTM+fMVncHgEvfN+33AADMQ9ABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDGLtYi521VVX9eHDhy/mkgCXvOPHjz/X3YvdjruoQT98+HA2Nzcv5pIAl7yq+te9HOeSC8AgBB1gEIIOMAhBBxiEoAMMYtegV9WxqjpVVU+e5bU7q6qr6qoLMx4Ae7WXM/R7khzZubOqrk1yU5JnZp4JgBXsGvTufijJ82d56feSvCOJf8MO4BvASg8WVdUbk3yhuz9eVbsdu5FkI0kOHTq0ynIXXdXd+z3CULrv3O8R4LKw9C9Fq+qKJO9K8pt7Ob67j3b3enevLxa7PrkKwIpWucvle5Jcl+TjVfV0koNJHq2q75hzMACWs/Qll+5+Iskrv7K9HfX17n5uxrkAWNJeblu8L8nDSa6vqpNVdceFHwuAZe16ht7dt+3y+uHZpgFgZZ4UBRiEoAMMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQQg6wCAEHWAQgg4wCEEHGISgAwxC0AEGIegAgxB0gEHsGvSqOlZVp6rqyTP2/U5VfaqqPlFVf1VVL7+wYwKwm72cod+T5MiOfQ8keW13f3+Sf07yzpnnAmBJuwa9ux9K8vyOfR/t7tPbm/+Y5OAFmA2AJazN8Bm/lOQvzvViVW0k2UiSQ4cOzbAcXL6q7t7vEYbSfed+jzCrSb8Urap3JTmd5N5zHdPdR7t7vbvXF4vFlOUAOI+Vz9Cr6vYkNye5sbt7vpEAWMVKQa+qI0l+I8mPd/f/zDsSAKvYy22L9yV5OMn1VXWyqu5I8gdJXpbkgap6vKr++ALPCcAudj1D7+7bzrL7fRdgFgAm8KQowCAEHWAQgg4wCEEHGISgAwxC0AEGIegAgxB0gEEIOsAgBB1gEIIOMAhBBxiEoAMMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYxK5Br6pjVXWqqp48Y9+3V9UDVfWZ7a+vuLBjArCbvZyh35PkyI59dyX5WHe/OsnHtrcB2Ee7Br27H0ry/I7dtyR5//b370/yczPPBcCSVr2GfnV3P5sk219fOd9IAKzigv9StKo2qmqzqja3trYu9HIAl61Vg/5vVfWdSbL99dS5Duzuo9293t3ri8VixeUA2M2qQf9wktu3v789yYfmGQeAVe3ltsX7kjyc5PqqOllVdyT57SQ3VdVnkty0vQ3APlrb7YDuvu0cL9048ywATOBJUYBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQQg6wCAEHWAQgg4wCEEHGISgAwxC0AEGIegAgxB0gEEIOsAgBB1gEIIOMAhBBxjEpKBX1a9W1VNV9WRV3VdV3zLXYAAsZ+WgV9U1Sd6WZL27X5vkQJJb5xoMgOVMveSyluRbq2otyRVJvjh9JABWsXLQu/sLSe5O8kySZ5P8R3d/dOdxVbVRVZtVtbm1tbX6pACc15RLLq9IckuS65J8V5KXVtWbdh7X3Ue7e7271xeLxeqTAnBeUy65/GSSf+nure7+3yQfTPKj84wFwLKmBP2ZJD9SVVdUVSW5McmJecYCYFlTrqE/kuT+JI8meWL7s47ONBcAS1qb8ubufneSd880CwATeFIUYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQQg6wCAEHWAQgg4wCEEHGISgAwxC0AEGIegAgxB0gEEIOsAgBB1gEIIOMAhBBxiEoAMMYlLQq+rlVXV/VX2qqk5U1evmGgyA5axNfP97kvx1d/98Vb0kyRUzzATAClYOelV9W5IfS/KLSdLdLyZ5cZ6xAFjWlEsu351kK8mfVtVjVfXeqnrpTHMBsKQpQV9L8oNJ/qi7b0jy30nu2nlQVW1U1WZVbW5tbU1YDoDzmRL0k0lOdvcj29v358uB/xrdfbS717t7fbFYTFgOgPNZOejd/aUkn6+q67d33Zjkk7NMBcDSpt7l8tYk927f4fK5JG+ZPhIAq5gU9O5+PMn6TLMAMIEnRQEGIegAgxB0gEEIOsAgBB1gEIIOMAhBBxiEoAMMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQQg6wCAEHWAQk4NeVQeq6rGq+sgcAwGwmjnO0N+e5MQMnwPABJOCXlUHk/xskvfOMw4Aq5p6hv77Sd6R5P/OdUBVbVTVZlVtbm1tTVwOgHNZOehVdXOSU919/HzHdffR7l7v7vXFYrHqcgDsYsoZ+uuTvLGqnk7ygSRvqKo/n2UqAJa2ctC7+53dfbC7Dye5NcnfdvebZpsMgKW4Dx1gEGtzfEh3P5jkwTk+C4DVOEMHGISgAwxC0AEGIegAgxB0gEEIOsAgBB1gEIIOMAhBBxiEoAMMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQawc9Kq6tqr+rqpOVNVTVfX2OQcDYDlrE957OsmvdfejVfWyJMer6oHu/uRMswGwhJXP0Lv72e5+dPv7/0pyIsk1cw0GwHJmuYZeVYeT3JDkkbO8tlFVm1W1ubW1NcdyAJzF5KBX1ZVJ/jLJr3T3f+58vbuPdvd6d68vFoupywFwDpOCXlXfnC/H/N7u/uA8IwGwiil3uVSS9yU50d2/O99IAKxiyhn665P8QpI3VNXj2//9zExzAbCklW9b7O6/T1IzzgLABJ4UBRiEoAMMQtABBiHoAIMQdIBBCDrAIAQdYBCCDjAIQQcYhKADDELQAQYh6ACDEHSAQQg6wCAEHWAQgg4wCEEHGISgAwxC0AEGIegAg5gU9Ko6UlWfrqrPVtVdcw0FwPJWDnpVHUjyh0l+OslrktxWVa+ZazAAljPlDP2Hkny2uz/X3S8m+UCSW+YZC4BlrU147zVJPn/G9skkP7zzoKraSLKxvflCVX16wpp8rauSPLffQ+ym6tf3ewQuPj+b83rVXg6aEvQ6y77+uh3dR5McnbAO51BVm929vt9zwE5+NvfHlEsuJ5Nce8b2wSRfnDYOAKuaEvR/SvLqqrquql6S5NYkH55nLACWtfIll+4+XVW/nORvkhxIcqy7n5ptMvbCpSy+UfnZ3AfV/XWXvQG4BHlSFGAQgg4wCEEHGISgAwxC0IHJqup7q+rGqrpyx/4j+zXT5UjQB1BVb9nvGbh8VdXbknwoyVuTPFlVZ/5Np9/an6kuT25bHEBVPdPdh/Z7Di5PVfVEktd19wtVdTjJ/Un+rLvfU1WPdfcN+zrgZWTK33LhIqqqT5zrpSRXX8xZYIcD3f1CknT301X1E0nur6pX5ex/84kLRNAvHVcn+akk/75jfyX5h4s/DnzVl6rqB7r78STZPlO/OcmxJN+3v6NdXgT90vGRJFd+5X+aM1XVgxd/HPiqNyc5feaO7j6d5M1V9Sf7M9LlyTV0gEG4ywVgEIIOMAhBBxiEoAMM4v8B4g8O22F90f8AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ovarian['ecog_ps'].value_counts().sort_index().plot.bar(color='darkblue')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Modelling" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the negative log-likehood function of the Weibull model." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def nll(parameters, scale, X, Y):\n", " if scale is None:\n", " beta, sigma = parameters[:-1], np.exp(parameters[-1])\n", " else:\n", " beta, sigma = parameters, scale \n", " z = (np.log(Y[:, 0]) - np.dot(X, beta)) / sigma\n", " return np.sum(np.exp(z) - Y[:, 1] * (z - np.log(sigma)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define gradient and Hessian of the negative log-likehood w.r.t. `parameters` using [Autograd](https://github.com/HIPS/autograd)." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "nll_grad = grad(nll, 0)\n", "nll_hess = hessian(nll, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define a function to perform model fitting and hypothesis testing on the parameters." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def fit(X, Y, scale=None, x0=None):\n", " if x0 is None:\n", " p = X.shape[1]\n", " if scale is None:\n", " p += 1\n", " x0 = np.zeros(p, dtype=np.double)\n", " solution = minimize(nll, x0, args=(scale, X, Y), method='Newton-CG', jac=nll_grad, hess=nll_hess)\n", " if not solution.success:\n", " raise RuntimeError(solution.message)\n", " results = {}\n", " results['coef'] = solution.x\n", " results['vcov'] = np.linalg.inv(nll_hess(solution.x, scale, X, Y))\n", " results['se'] = np.sqrt(np.diag(results['vcov']))\n", " results['z'] = results['coef'] / results['se']\n", " results['p'] = 2 * norm.cdf(-np.abs(results['z']))\n", " return results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define a function to specify models using [`patsy`](https://patsy.readthedocs.io/) formulas." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def fit_formula(formula, data={}, scale=None):\n", " Y, X = dmatrices(formula, data=data)\n", " results = fit(np.asarray(X), np.asarray(Y), scale)\n", " index = X.design_info.column_names\n", " if scale is None:\n", " index.append('Log(scale)')\n", " return pd.DataFrame(results, index=index, columns=['coef', 'se', 'z', 'p'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Estimate a Weibull model." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefsezp
Intercept6.8966931.1775475.8568294.717874e-09
ecog_ps-0.3850430.527002-0.7306294.650061e-01
rx0.5286450.5292020.9989493.178196e-01
Log(scale)-0.1234420.252213-0.4894356.245341e-01
\n", "
" ], "text/plain": [ " coef se z p\n", "Intercept 6.896693 1.177547 5.856829 4.717874e-09\n", "ecog_ps -0.385043 0.527002 -0.730629 4.650061e-01\n", "rx 0.528645 0.529202 0.998949 3.178196e-01\n", "Log(scale) -0.123442 0.252213 -0.489435 6.245341e-01" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fit_formula('futime + fustat ~ ecog_ps + rx', ovarian)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare with the output of `survreg` in R:\n", "```\n", "Call:\n", "survreg(formula = Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian)\n", " Value Std. Error z p\n", "(Intercept) 6.897 1.178 5.857 4.72e-09\n", "ecog.ps -0.385 0.527 -0.731 4.65e-01\n", "rx 0.529 0.529 0.999 3.18e-01\n", "Log(scale) -0.123 0.252 -0.489 6.25e-01\n", "\n", "Scale= 0.884\n", "\n", "Weibull distribution\n", "Loglik(model)= -97.1 Loglik(intercept only)= -98\n", "\tChisq= 1.74 on 2 degrees of freedom, p= 0.42\n", "Number of Newton-Raphson Iterations: 5\n", "n= 26\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fix `scale` to 1, corresponding to an exponential model." ] }, { "cell_type": "code", "execution_count": 13, "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", "
coefsezp
Intercept6.9618381.3218775.2666291.389517e-07
ecog_ps-0.4331350.586994-0.7378864.605834e-01
rx0.5815030.5869940.9906463.218586e-01
\n", "
" ], "text/plain": [ " coef se z p\n", "Intercept 6.961838 1.321877 5.266629 1.389517e-07\n", "ecog_ps -0.433135 0.586994 -0.737886 4.605834e-01\n", "rx 0.581503 0.586994 0.990646 3.218586e-01" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fit_formula('futime + fustat ~ ecog_ps + rx', ovarian, scale=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare with the output of `survreg` in R:\n", "```\n", "Call:\n", "survreg(formula = Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian,\n", " dist = \"exponential\")\n", " Value Std. Error z p\n", "(Intercept) 6.962 1.322 5.267 1.39e-07\n", "ecog.ps -0.433 0.587 -0.738 4.61e-01\n", "rx 0.582 0.587 0.991 3.22e-01\n", "\n", "Scale fixed at 1\n", "\n", "Exponential distribution\n", "Loglik(model)= -97.2 Loglik(intercept only)= -98\n", "\tChisq= 1.67 on 2 degrees of freedom, p= 0.43\n", "Number of Newton-Raphson Iterations: 4\n", "n= 26\n", "```" ] } ], "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.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }