Survival Example.ipynb

{
 "metadata": {
  "name": "Survival Example"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# Cox regression in PyMC\n",
      "\n",
      "This example implements the Cox model using the counting process notation introduced by Andersen and Gill (1982) "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Leukemia data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Survival times in weeks\n",
      "obs_t = (1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12, 12, 15, 17, 22, 23, 6, \n",
      "    6, 6, 6, 7, 9, 10, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 32, 32, 34, 35)\n",
      "# Failure (1) or censored (0)\n",
      "fail = (1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, \n",
      "    0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0)\n",
      "# Indicators for placebo or treatment\n",
      "Z = (0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, \n",
      "    0.5, 0.5, 0.5, 0.5, 0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, \n",
      "    -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5)\n",
      "# Unique failure times\n",
      "t = (1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 15, 16, 17, 22, 23, 35)\n",
      "N = len(obs_t)\n",
      "T = len(t) - 1"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Calculate risk set"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Y = np.array([[int(obs >= time) for time in t] for obs in obs_t])"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Counting process. Jump = 1 if $\\text{obs}_t \\in [ t_j, t_{j+1} )$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dN = np.array([[Y[i,j]*(t[j+1] >= obs_t[i])*fail[i] for i in range(N)] for j in range(T)])"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Model using Poisson trick:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pymc import Normal, Lambda, Poisson, MCMC\n",
      "\n",
      "beta0 = Normal('beta0', 0, 0.001, value=np.zeros(T))\n",
      "beta1 = Normal('beta1', 0, 0.00001, value=0)\n",
      "\n",
      "# Poisson trick: independent log-normal hazard increments\n",
      "Idt = Lambda('Idt', lambda b0=beta0, b1=beta1: [[Y[i,j]*np.exp(b0[j] + b1*Z[i]) for i in range(N)] for j in range(T)])\n",
      "\n",
      "# Likelihood\n",
      "dN_like = Poisson('dN_like', Idt, value=dN, observed=True)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Initialize MCMC sampler."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M = MCMC([beta0, beta1, Idt, dN_like])"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M.sample(10000, burn=5000)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "('\\r', <pymc.progressbar.ProgressBar instance at 0x10c8ae680>)\n"
       ]
      }
     ],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "beta1.summary()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pymc import Matplot\n",
      "\n",
      "Matplot.plot(beta1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [],
     "prompt_number": "*"
    }
   ],
   "metadata": {}
  }
 ]
}