Phase II Vaccine Power Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Phase II Power Analysis\n",
    "\n",
    "Using the sample information from the phase I trial, we can calculate power for a phase II trial.\n",
    "\n",
    "![Table 2](http://d.pr/i/O9Y6+)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proper Bayesian power analysis\n",
    "\n",
    "We wish to identify sample sizes for high dose and standard dose (in a 2:1 ratio) that will be able to detect a positive difference in the probability $p$ of a $\\ge 4$ fold titer rise to H3N2 influenza. That is, we wish to be able to detect:\n",
    "\n",
    "$$\\delta = p^{(HD)} - p^{(SD)} > 0$$\n",
    "\n",
    "with 80% power. This analysis will be conducted using the *proper Bayesian* approach (Spiegelhalter et al. 2004). The advantage of this approach for power analysis over a frequentist power analysis is that it does not assume that the underlying probabilities are known. Failure to do so will tend to underestimate variance and result in underpowered experiments.\n",
    "\n",
    "A sample leading to a successful trial (one in which a difference is detected) under a particular sample size $n$ is:\n",
    "\n",
    "$$E_n = \\left\\{\\mathbf{x}; Pr( \\delta > 0 | \\mathbf{x} ) \\ge 0.95\\right\\}$$\n",
    "\n",
    "so the probability of a successful trial is:\n",
    "\n",
    "$$\\psi(n) = P_{X^{(n)}}(X^{(n)} \\in E_n)$$\n",
    "\n",
    "which is the **expected Bayesian power**.\n",
    "\n",
    "In our case, the power is:\n",
    "\n",
    "$$\\psi(n) = \\Phi\\left(\\frac{\\pi_{hd} - \\pi_{sd}}{\\sqrt{\\pi_{hd}(1-\\pi_{hd})/n_{hd} - \\pi_{sd}(1-\\pi_{sd})/n_{sd}}}\\right)$$\n",
    "\n",
    "where the $\\pi_i$ are the posterior probabilities, parameterized by a beta distribution with a binomial likelihood for the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the sample data from the phase I trial (GiaQuinta et al. 2014) to construct prior distributions for $p^{(HD)}$ and $p^{(SD)}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import pymc3 as pm\n",
    "import theano.tensor as tt\n",
    "import matplotlib.pyplot as plt\n",
    "from theano.tensor.shared_randomstreams import RandomStreams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x78ee8af53240>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VNeZ5/HvqSrtu9CCQIDEIhAIY0BsNraxWQ1e4i2xE9uxQ0wn7U7S6ZnJTGfppHs83Ul3J3EnnXTidhLHjpM4tuMlNovB4AUbMMIsEgYkEFqR0C60S6U688eVsMCASlJV3bq33s/z1FNCKm69F8SPo/eee47SWiOEEMI6HGYXIIQQYmQkuIUQwmIkuIUQwmIkuIUQwmIkuIUQwmIkuIUQwmIkuIUQwmIkuIUQwmIkuIUQwmJc/jhoSkqKzsrK8sehhRDClg4cONCgtU715rV+Ce6srCwKCgr8cWghhLAlpVS5t6+VVokQQliMBLcQQliMBLcQQliMBLcQQliMBLcQQliMBLcQQliMBLcQQliMX+Zx24rW0NUMnU3Q1QSdjcbHnY0QlwFz7walzK5SCBFCJLiH89rX4cBvLv913Q/z7g1cPUKIkCfBfSVnPzJCe86dkLMOosdBdJLxHJkIf7gXNv8vmHINJE42u1ohRIiQ4L6Sd/4NwmNhww8hOvmTX7/jF/Bfy+GlL8PnXwWHM/A1CiFCjlycvJy643D0JVi86dKhDZCUBTf/AMp3w56fBbQ8IUTokuC+nHf+DcKiYdnfXPl1V38WZt0CO/8v1BYFpjYhREiT4L6U+mIoehEWfxFixl35tUrBrf9h9Lz/vAncPYGpUQgRsiS4L+Xdf4ewKLjmq969PiYFbv9PqDsKOx/zb21CiJAnwX2xhpNQ+Dws2mgEsrdy1sLCh+H9n0LZbv/VJ4QIeRLcF3v3h+CM8H60PdSaxyA5G176EvS0+b42IYRAgvtCTaVw5DnI/wLEpo3890fEwu0/h9ZKOPR739cnhBBIcF/onR+CMwyuHcVoe9CUZTBhPhT8xrhdXgghfEyCe1BzGRz+Ayx8COLGj+1YCx+G+mNQ+YEvKhNCiAtIcA9694fgcMG1fzv2Y+XdBeFxV17jRAghRkmCG6C/DwpfgHmfgfiMsR8vIhau+rRx52VX89iPJ4QQQ0hwA9Qcgb5OmHaT746Z/zC4u+HwH313TCGEQILbULHHeJ68zHfHHD8XJubLRUohhM9JcIMR3EnZY78oebH8h6HhxMf/MQghhA9IcGttBKsvR9uD5twJEQnGqFsIIXxEgrvxpLEN2eSlvj92eLRxwfOjV4ztzoQQwgckuAfbGFOu8c/xFz4E/T1yJ6UQwmckuMv3GFuRjZvun+Onz4HMxXDgKblIKYTwCQnuwf62P3dqz38YGktk1UAhhE+EdnC31ULzaf/0t4eacwdEJhijbiGEGCOvg1sp5VRKHVRKvebPggLq/PxtP/W3B4VFwbz74Nir0NHo3/cSQtjeSEbcXwOO+asQU1TsBVcUZFzl//da+DD09xrLxgohxBh4FdxKqUxgA/Ckf8sJsIo9kJlvLOXqb2mzIH2uMeoWQogx8HbE/TjwDcDjx1oCq/sc1Bb6bxrgpcy+zRjlt50N3HsKIWxn2OBWSt0C1GmtDwzzuk1KqQKlVEF9fb3PCvSbqv2gPf6/MDlU7q2AhuP2uUwghAg8b0bc1wK3KaXKgD8CNymlfnfxi7TWT2it87XW+ampqT4u0w8q9oJyQOaiwL1n6ixjvvixvwTuPYUQtjNscGut/15rnam1zgLuBXZqre/3e2X+VrEHxl8FEXGBe0+ljFF32btyC7wQYtRCcx63uxeqCvyzsNRwcm8FjxuKtwX+vYUQtuAayYu11m8Bb/mlkkCqPQLursD2twdNWADxE412ydX3Bf79xYg0dfRysKKZDyuaqWru4pHrppI3McHsskSIG1Fw24Y/Nk7w1mC75MBT0NNubHMmgkZ7j5tXDlVzoLyZgxUtnG7oAMDpUESHOdlSVMt3b53NZxdPRvlzmQQhriA0g7t8DyRPhbh0c94/9zbY9ws4ud24HV4EhR53P194aj8fnG4iJTac+ZOT+HT+JBZMTuSqzEQ6e918/U+H+dZLRXxwuol/vmMuMRGh+U9ImCv0vusGN06YebN5NUxeCtEpRrtEgjsoaK35+xcL+eB0Ez+8Zx53Lpj4iRF1VLiTpx5axM92neTHO4opqm7l559byMzxAbzALQSheHGyoQS6mszpbw9yOGHWBuMCZV+3eXWI83668yR/PljN363O4a6FmZdtgzgciq+snMHvNi6htcvN7T/bzQsHqgJcrQh1oRfcFe8bz/5eWGo4ubdBbzucftvcOgSvHKrmR9uLuXPBRL5yk3frsl8zPYXNX13OvMxE/ufzh3nqvdN+rlKIj4VgcO812hTjpplbR/b1EBEva5eYbH9ZE//r+SMszk7mX+6cO6ILjmnxkTz7xSXckJPKv79RTN05+elJBEYIBvceo01i9owAVzjkrIPjm6HfbW4tIaqsoYNNTxeQmRTFEw8sJMLlHPExXE4H/3jbHHrdHv55s70WzxTBK7SCu70OmsvMmQZ4Kbm3Gv328vfMriTkNHf08vBT+wH49UOLSIwOH/WxslJi2HT9VF4+dIZ9pbLeuvC/0Aru+uPG8/g8c+sYNH2VsR64rF0ScN98qZDq5i6eeDCfrJSYMR/v0RunMzExiu++ehR3v30W0RTBKbSCu6HYeE7JMbeOQeHRMGOVsVqgR/6xB0pRdStbimr58oppLMpK9skxo8KdfOeWXI7XtvHM3nKfHFOIywmx4C6B8FiIyzC7ko/l3gZtNVB9xVVzhQ89vqOE+EgXX1ie7dPjrp0znutmpPCjN4qpb+vx6bGFGCrEgrsYUmaYf2FyqBlrwBEGx14xu5KQUFjVyo5jZ/nidVNJiPLtzkdKKb532xy63f38YOtxnx5biKFCLLhLYNwMs6u4UFQiZF8HJ7aaXUlIeHxHMQlRYTx8bZZfjj8tNZaNy6fywoEqDpQ3++U9hAid4O7tgNbK4OlvDzVjLTSWQOMpsyuxtcOVLbx5vI5HrssmLtJ/+4x+5abpjI+P5LuvFtHv0X57HxG6Qie4G08azylBNuIGyFljPJe8YW4dNvf4jmISo8P4/DVZfn2fmAgX39qQS1H1Of7wQYVf30uEptAJ7oYS4zkYR9zJUyFlpmyu4EcHK5rZdaKeR66b6tfR9qBbrsogf0oSv3j7lIy6hc+FUHAXG3tMJk81u5JLy1kDZbuhp83sSmzp8R0lJAVgtD1IKcVD12ZR1dzFO8UW2DxbWEpoBXfiFAiLNLuSS8tZB54+OLXL7Eps50B5M28X17Pp+mnEBnD97DWzx5MSG8HvZF638LEQCu6S4GyTDJq0BCISoETaJb72+I5ikmPCeXDZlIC+b7jLwb2LJrHzRB1VzZ0BfW9hb6ER3J5+4+JkMF6YHOQMg+krofgNuYvShw6UN/FuSQN/df1UU3aruXfxJAD++EFlwN9b2FdoBHdrJbi7g3vEDZCzFjrqoOaQ2ZXYxuM7SkiJDeeBAI+2B2UmRXPTzDT+uL+SXrf8hyx8IzSCO5hnlAw1fTWgZHaJj5Q1dPBuSQMPXZNFdLh5u/Tdv3QKDe09bP/orGk1CHsJkeAOssWlLidmHGQukj63j7xwoAqHgrsXTjK1jutzUpmYGCUXKYXPhE5wRyUbwRjsctbCmYPQVmt2JZbW79G8cKCKG3JSGZ9g7kwip0Px2SWT2VPayMm6dlNrEfYQIsEd5DNKhspZazyXbDe3Dot7p6Se2nPdfGaRuaPtQZ9ZNIkwp+LZfTLqFmMXIsFdHNwzSoZKz4P4iVAsi06NxfMFlSTHhHPTrHSzSwEgJTaCdXkZvHigiq7efrPLERZn/+DubIKOeuuMuJUylno9tQvcsqbzaDR19LL9o7PcMX8i4a7g+Rb/3JLJnOt285cjZ8wuRVhc8HxX+8v5xaUsEtxg3EXZ1yF7UY7Sywer6evXfDo/ONokg5ZkJzMjLZZn5SKlGCP7B/f5GSUWaZUAZF8PrkiZFjgKWmv+VFDJvMwEZo6PM7ucCyil+NySyRyuaqWwqtXscoSFhUBwlxg7zCSacwPGqIRHG+FdvBW0rCw3EkXV5zhe28Y9QTbaHnTnwkyiwpxykVKMSWgE97hp4DTvBoxRmbEGmss+vnlIeOW5ggoiXA5unTfB7FIuKT4yjFvnZfCXw2fo7pOLlGJ0QiC4LTSjZKjz0wKlXeKt7r5+Xjl0hvVzM3y+n6Qv3TZvIh29/bx1QpZ7FaNj7+Du74Pm09a6MDkocTKkzZY+9whsO1pLW7ebe/IzzS7lipZOTSY5JpzXC2vMLkVYlL2Du+k0eNzWDG4wRt0Ve6CrxexKLOFPBZVMSo5iaXZw3yHrcjpYlzeeN4+dlXaJGJVhg1spFamU+kApdVgpdVQp9Y+BKMwnrDijZKgZa43/eE7tNLuSoFfZ1Ml7Jxu5Z+EkHA5ldjnD2jA3g87eft46UWd2KcKCvBlx9wA3aa3nAVcD65RSS/1blo8MBvc4iwZ35iKISpJNhL3wwoEqlIK7FgZ3m2TQkuxkxsWE89oRaZeIkRs2uLVhcGWcsIGHNeaoNZRAXAZExptdyeg4XTB9lRHcHvmR+nI8AwtKLZ+ewsTEKLPL8crH7ZI6uQVejJhXPW6llFMpdQioA7Zrrff5tywfseqMkqFy1kFnI1R/aHYlQWt/WRPVLV3cbZHR9qANczPo6pN2iRg5r4Jba92vtb4ayAQWK6XyLn6NUmqTUqpAKVVQXx8E05y0ttaqgJcz7SZjd3pZdOqyNhfWEOFysCo3OBaU8tbi7GRSYsN5TWaXiBEa0awSrXUL8Baw7hJfe0Jrna+1zk9NTfVReWPQXgc9rdYP7uhkmLRU5nNfRr9Hs7molptmpZmyp+RYDLZLdkq7RIyQN7NKUpVSiQMfRwGrgOP+LmzMrD6jZKicNVBbCK3VZlcSdArKmqhv62H93AyzSxmVDXMn0NXXzy5pl4gR8GbEnQHsUkodAfZj9Lhf829ZPmCV7cq8kTPwA47MLvmE1wtriAxzcNOsNLNLGRWjXRLB6zK7RIzAsD9baq2PAPMDUItvNZRAWAzEBeeaFSOSOgsSJht3UeY/bHY1QaPfo9lSVMuNM63XJhnkdChuzhvP8wcq6ex1m7qpsbAO+9452VAMKdPBYYNTVMq4i/L029DXbXY1QWP/QJtkw1XWbJMM2nBVBt19HnYdD4KL+sISbJBql2GHGSVD5ayFvk4o2212JUFjs8XbJIMWZSWTGhfB64WyM47wjj2Du7cTWivsFdxZy8EVJdMCB/R7NJsLjdkkVm8vDLZLdh6vo6PHbXY5wgLsGdzntyuzwYySQWFRMHWFMS1QNldgf1kTDe3WnU1ysQ1zjXbJzuMyu0QMz97BPW66uXX4Ws4aaKmA+uCfjelvrx+xR5tkUH5WMmlxEWyWm3GEF+wZ3C0D20IlZZlahs/NGNhcIcTX6B6cTWKHNskgaZeIkbBncDeXQ1QyRATXZrFjljAR0ueGfHB/cNpok2yYa4OpnkOsy8ugx+3h3ZIGs0sRQc6ewd1SbuwgY0c5a6FyH3Q2mV2JaQZnk9w4KwiWVvChRVlJJESFsf2js2aXIoKcTYO7ApIstKv7SOSsBd0fspsrDLZJVs5Kt02bZJDLafTsdx4/S79HLkCLy7NfcHs8RnDbdcQ9cSFEjwvZdslgm8Qus0kutio3nebOPj6saDa7FBHE7Bfc7WehvxcSbTridjhh+mo4uR36Q+8i1uuFZ4gKc9quTTLo+pwUwpxK2iXiiuwX3HadUTLUzHXQ1Wz0ukNIv0ez1WazSS4WFxnGsmkp7JDgFldgw+CuMJ7t2ioBYzszZzic2Gx2JQG173QjDe29tm2TDFqdm0ZpQwen6tuHf7EISfYL7uaBEbedgzsizriL8vhrIXUX5baiWiJc9ptNcrGVAzv5yKhbXI79grulHGLSjFvE7Wzmemgug7qPzK4kIDwezbajZ7khJ9W2bZJBExKjyJsYL31ucVn2DG67TgUcauZ6QMHx182uJCCOVLdSe66bdXnjzS4lIFblpnOgopnG9h6zSxFByH7B3Wzjm2+GikuHzEVGuyQEbC2qxeVQrJxlrQ2BR2tVbjpaI4tOiUuyV3D3u+FctX2nAl5s1gaoOQytVWZX4ldaa7YdrWXZtHEkRIeZXU5AzJkQz4SESGmXiEuyV3C3nQGPOzRG3ACzbjGej9t7dklJXTunGzpYMyc02iQASilWzU7n3ZIGuvtkB3hxIXsF9+BUwFDocYOxNVvKTNu3S7YW1aIUrJ0dGm2SQaty0+nq6+f9U7LolLiQvYL7/FTAEAluMNolZbuNG3JsatvRWhZMTiItPtLsUgJqydRkYiNc0i4Rn2Cv4G6pABQkZJpdSeDMusVYdKr4DbMr8YvKpk6OnjnH2jmhNdoGiHA5uSEnlR3H6vDIolNiCJsFdznETwBXhNmVBM6E+RCXYdt2ybajtQCsDaH+9lCrZ6dT39bDkepWs0sRQcRmwW3jVQEvx+Ew5nSffBP6us2uxue2Ha0lNyOeKeNizC7FFCtmpuJ0KLmLUlzAXsHdXB5a/e1Bs9ZDXwecftvsSnyqrq2bgvLmkGyTDEqMDmdRVpL0ucUF7BPc7t6BOdwhNuIGyLoeIuJt1y7Z/tFZtCZk7pa8nNWzx3PibBsVjZ1mlyKChH2C+1wVoENnKuBQrnCYsRpObAGPfeb8bjt6lqxx0cxMt9neoSO0enDRqWMy6hYG+wR3KE4FHGrWBuioh6r9ZlfiE61dfbx/soG1c8ajlDK7HFNNHhdNTnqsBLc4zz7BHQrrcF/J9NXgCLPNolM7j5/F7dGsDfE2yaBVuensO91Ea2ef2aWIIGCj4C4H5YT4iWZXYo7IeJh6g23W6N5aVEt6fARXZyaaXUpQWDU7nX6P5q1iWXRK2Cq4KyBhIjjtvVbzFc3aAE2lUH/c7ErGpKu3n7eL61kzezwOR2i3SQZdnZlISmw4O45JcAs7BXeoTgUcauYGUA44+rLZlYzJ28X1dPd5Qn42yVCOgSVt3zpRR6/bY3Y5wmT2Ce6WCgnuuHTIWg5FL1i6XbK1qIbE6DAWZyebXUpQWZmbRlu3m/1lTWaXIkxmj+Du64L22tCcCnixvLug8STUHjG7klHpcffz5rE61sxOJ8xpj29PX1k+I4UIl0NuxhHDB7dSapJSapdS6phS6qhS6muBKGxEWiqN51AfcQPk3gYOFxS9aHYlo/LeyQbaetzcnGfvndxHIzrcxfLpKew4dhZt4Z+oxNh5M6RxA/9Da50LLAUeVUrN9m9ZIxTqUwGHik6GaSuh6M/gsV4vdHNhLXGRLq6ZPs7sUoLSqtnpVDV3ceJsm9mlCBMNG9xa6xqt9YcDH7cBx4DgmnPXUmY8S6vEkHcXtFZa7macXreHN47Wsjo3nQiX0+xygtLKWWkAvCmzS0LaiJqISqksYD6w7xJf26SUKlBKFdTX1/umOm+1VIAzHGJlFgJgLDrlirRcu2RPaSPnut3cPFfaJJeTFh/JvEmJ0ucOcV4Ht1IqFngR+Fut9bmLv661fkJrna+1zk9NTfVljcNrLoeEScYSpwIi4iBnLRx9ydhA2SK2FNYQE+7kuhkpZpcS1FbnpnGosoW6Nvst4yu841XSKaXCMEL7Wa31n/1b0iiE4jrcw8m7CzrqoHy32ZV4xd3vYdvRWlbmphMZJm2SK1k1sPfmTmmXhCxvZpUo4FfAMa31j/xf0ii0lEt/+2Iz1kB4nGXaJftON9Hc2cf6udLuGs7M9Dgyk6Jk0akQ5s2I+1rgAeAmpdShgcd6P9flvZ526GyUEffFwqKMW+A/etVYqzzIbS6sISrMyQ05aWaXEvSUUqzKTefdkga6eu2zjK/wnjezSnZrrZXW+iqt9dUDj82BKM4r56cCyoj7E/Lugu4WOLXT7EquqN+j2Xb0LDfOSiUqXNok3lg9O50et4fdJxvMLkWYwPpX81oG1uFOyjK1jKA0dQVEJQV9u6SgrImG9h656WYEFmcnExfpkr0oQ5QNgltuvrksVzjMvt1Yo7s3eLe92lJUS4TLwY2zpE3irTCngxUz03jz+Fk8HrmLMtRYP7iby8EVBTEBnoJoFXl3GRsJl2wzu5JL8ng0W4pquCEnldiIEF6SdxRW5abR0N7LoaoWs0sRAWb94G4pN0bbIb691WVNuRZi06HwBbMruaSDlc2cPdfDernpZsRWzEzD5VC8cVTaJaHGHsEtUwEvz+GEOXdAyXbobjW7mk/YXFhLuNPBTbnSJhmphKgwlk0bx9aiGll0KsTYILjl5pth5d0N/T1Btx+l1pothTVcNyOF+Mgws8uxpJvzMihr7OR4rSw6FUqsHdxdLcYoUqYCXllmPiRlw8Fnza7kAoerWjnT2i1rk4zBmjnpOJRxgVeEDmsH9+CMEmmVXJlSsOBB4/b3hhKzqzlvS2ENLodidW662aVYVkpsBIuzk9lSWGN2KSKALB7cA3O4pVUyvKs/Z2yw8OFvza4EMGaTvHakhmunp5AQLW2Ssbg5L4OSunZO1km7JFRYO7ibB4NbRtzDikuHmTfDod+Du8fsaigob6a6pYs75gfX0u5WNLip8pZCaZeECmsHd0uFsZBSVJLZlVjDgoeMdV2C4CLlSweriQpzsnq2tEnGKj0+koVTkqTPHUKsHdzNZZCcJXO4vTXtRmPdcpPbJT3ufl4/coa1c9KJkZtufOLmvPF8VHOO8sYOs0sRAWD94JY1SrzncML8B6D0LWg6bVoZu47Xc67bzaekTeIz59slMuoOCdYNbo9n4OabLLMrsZb594NywIdPm1bCK4eqSYkNZ/l02enGVzKTopmXmSDBHSKsG9ztZ8HdLcE9UgkTjU0WDj0L/X0Bf/vWrj7ePFbHrfMm4HJa99svGK3Ly+BwZQvVLV1mlyL8zLr/cprLjGcJ7pFb8HnjP77iwC88taWwht5+j8wm8YObB9olW2XUbXs2CO5sU8uwpBlrIC4DDjwV8Ld+6WA1U1NimDsxIeDvbXdZKTHkZsTLzTghwOLBrYxZEmJknC6j131yB7RUBuxtq1u62He6iU/Nn4iSmUB+cXPeeA5UNHP2nOwAb2fWDu6ETGOzADFy8x8wng/+LmBv+eqhMwB86mppk/jL+rnj0Rq2HZV2iZ1ZO7ilvz16SVOMed0HnwGP/zec1Vrz0sEqFk5JYvK4aL+/X6ianhbH9LRYuYvS5iwe3HKr+5gsfAjOVRstEz87VtNG8dl2mbsdAOvzxrPvdCON7eYvbSD8w5rB3dsJ7bUy4h6rnJuNLd/2/8rvb/XKoWpcDsUGWcLV79blZeDR8IZsJGxb1gzu88u5yoySMXGFQ/5GYz/KumN+e5t+j+aVQ2dYMTOV5Bi5JuFvuRlxZI2L5vUjMrvErqwZ3DKH23cWb4KwaHjvP/z2FvtKG6k91y1tkgBRSnH71RN571QDNa1yM44dSXCHuphxxg05hc/7bWrgSweriY1wsUo2TAiYuxZkojX8+cNqs0sRfmDd4A6PhehxZldiD8seNZ73/Mznh+7ocbOlqJZ1eeOJDHP6/Pji0iaPi2ZJdjIvHqiSjYRtyLrBnZQly7n6SuIkmHuPsdxrR6NPD/3SwWrae9x8donsUhRody/MpLShgw8rWswuRfiYtYNb+M61X4O+TvjgCZ8dUmvN7/aWM2dCPPMnJfrsuMI7N8/NICrMyQsHqswuRfiY9YJbawluf0jLhZnr4YNfQq9vFuPfX9bM8do2Hlw2RW5xN0FshIub547ntcNn6O7z/01WInCsF9ztdeDukuD2h+Vfh65mn63V/czecuIjXdw2T2aTmOXuhZm09bjlFnibsV5wy4wS/5m0GCZfA+//55jX6q5r62ZrUQ335E8iKlwuSpplafY4JiZG8aLMLrEVCW5xoeVfh3NVUPjCmA7zxw8q6evX3L9UliUwk8OhuGthJrtL6qltlRUD7cKiwS3LufrNjNWQNgfee9zYHm4U3P0efr+vgutmpJCdEuPjAsVI3bVgIh4Nfz4oFyntYtjgVkr9WilVp5QqCkRBw2oug/gJEBZpdiX2pJQx6q4/DsVbR3WIHcfOUnuumwdktB0UpoyLYXFWMi/InG7b8GbE/RSwzs91eE9mlPjfnDsgcTLs/pExi2eEntlbzsTEKFbKnZJB4+6FmZTWd3CwUuZ028Gwwa21fgdoCkAt3pHg9j+nC5b/HVTth2Ovjui3nqxr572TjXx2yWScDpkCGCzWXyVzuu3EWj3uvm5oOyPBHQjzHzB63W98x/hz99Lv9pYT5lR8ZpFcgwgmsREubs4bz19kTrct+Cy4lVKblFIFSqmC+vp6Xx32QueXc83yz/HFx5wuWPv/oKUc9v3Cq9/S0ePmxQNVrJ+bQUpshJ8LFCN198JM2rrdbJd1ui3PZ8GttX5Ca52vtc5PTU311WEvJFMBA2vajcZmC+/8u3Hj0zBePlRNW4+bB5fJRclgtHSqMaf7eWmXWJ61WiUS3IG35jHjTtVd/++KL9Na88yecnIz4lkwOSlAxYmRcDgUdy/M5N2Sekrr280uR4yBN9MB/wDsAWYqpaqUUhv9X9ZlNJcZi/7H+GlELz4pZbqx2cKHT0Pt5WeEvneyUdYlsYAHlk0hzOngv989bXYpYgy8mVVyn9Y6Q2sdprXO1Fr7f4PCy2k+Lcu5muGGb0BkAmz75iWnB2qt+eH2E2QkRHLnAlmXJJilxEZw98JMXvywiro2uZPSqqzXKpE2SeBFJcGKb8Lpt+HElk98+a0T9RysaOFvbppOhEvWJQl2j1w3lb5+D799v8zsUsQoWSe4ZTlXc+U/DCk58Ma3wd17/tNaa360vZjMpCjuWShTAK0gOyWGdXPG88yectp73GaXI0bBOsHdUW8s9C/BbQ5nGKz9Z2g6Bfv/+/ynt390lsLqVr66cgbhLut8O4W6TddP5Vy3m+f2+2efUeFf1vmXJjNKzDdjNUxbCW//ANrO4vEepaL5AAAQXklEQVQYo+3slBjulB3cLWX+5CSWZCfzq3dL6esf3WJiwjwS3GJk1n0f3D3wyl+ztaiG47VtfG3lDFxO63wrCcOXbpjGmdZuXjtyxuxSxAhZ51/bYHAnyqazpkrNMeZ2n9xB6eYfMz0tllvnTTC7KjEKK2amMjM9jl++XSqrBlqMtYI7LgPCosyuRCz6IrXp1/NI12/49mJkMSmLUkqx6fqpHK9t4+1iPy1TIfzCWsEtbZKg4PZoHm3bSJcjmhsK/95onQhLunXeBMbHR/LLt0vNLkWMgAS3GLGXD53hQFMYpdf8AHX2KLz5T2aXJEYp3OVg4/Js9pQ2cljW6rYMawR3Xzeck+Vcg0Ffv4efvFnC3IkJzF91L+RvhD3/Cad2mV2aGKV7F08iLtLFE+/IqNsqrBHcrZWAluAOAk/vKaeiqZO/W51jrEmy5jHjxpyXvwydwbPfhvBeXGQY9y+dwpaiGorPtpldjvCCNYJbpgIGhbKGDv5t23FumpXGipkDC32FR8NdT0JHA/zla6Pa6kyY75HrphIXGcZ3Xi6SGSYWIMEtvOLxaL7x4hHCnA7++Y65F64AmDEPbvq2sc3Z3p+bV6QYteSYcL6xbib7TjfxyiGZ1x3srBPcrkiIlc1nzfLM3nI+ON3EP9wym/EJkZ98wTVfhdxbjRUEC18IfIFizO5dNJl5kxJ57PVjtHb1mV2OuALrBLcs52qaisZOvr/lOCtmpnL3wsxLv8jhgDufhCnXwktfglM7A1ukGDOnQ/HY7Xk0dvTw4+3FZpcjriD4g1trqCqA9DlmVxKSjBbJYVwOxb/cOffKmySERcK9vzcuVj73AJw5GLhChU/MzUzggaVTeHpPGUXVrWaXIy4j+IO77hi018K0m8yuJCQ9u6+cvaVNfOeW2WQkeHHXalQi3P8iRCXD7+6GxlP+L1L41P9YM5PkmHC+/XIRHo9cqAxGwR/cgz9yT73R3DpCUGVTJ/+y5TjX56RyT/5lWiSXEp8BD/wZtAd+dye0ya7iVpIQFcY31+dyqLKF5wpk2ddgFPzBXboLUmZCgiwbGkgej+YbLxzBoRTfH65FcikpM+Bzzxu7wz97F3Sf80+hwi/umD+RxdnJ/GDrcZo6eof/DSKggju4+7qh7D2YJqPtQPvN+2XsKW3k2xtymZA4yoW9MvPh008b7a7f3grnanxbpPAbpRSPfSqP9m43P9hy3OxyxEWCO7gr94G7S/rbAba5sIbHXv+I1bPT+cyiMW5HNmM1fOZZaCiBJ1dCzRHfFCn8Lic9jo3Ls3muoJL3TzaYXY4YIriD+9ROcIQZU8xEQLx3soG//eMhFk5O4if3zh95i+RSZq6DL2w1Pv71uktuOCyC01dXzmBaagxffvZDTtW3m12OGBDcwV26CyYtgYhYsysJCUeqWtj0dAHZKTH86vOLiAr34Y7tGVfBIzuNjRj+cB+8/59ye7wFxES4eOrhxYQ5FQ//Zj8N7bKEbzAI3uDuaICawzBthdmVhITS+nYe+s1+EqPDeXrjYhKiw3z/JnHj4aHNkHsLvPEteO1voV/u0At2k5KjefLzi6hr6+aLvy2gu6/f7JJCXvAGd+lbxrP0t/2utrWbB371AQp4ZuNi0uMvcUu7r4RHwz1Pw/Kvw4GnjIuWDSX+ez/hE1dPSuTxz8zncFULX3/ukMzvNlnwBvepXRCZCBlXm12JrbV09vLgr/fR2tXHUw8vZmpqANpSDges+h7c+d9Q9xH817Xw9r+BW6adBbN1eeP51vpcthTV8v2tMtPETMEZ3Fob/e2pN4DDh31WcYHyxg4+9+Q+yho6eeKBhczNTAhsAVd9Gh7dD7PWw67H4JfXQ+UHga1BjMjG5dk8uGwKT7xTyjN7y80uJ2QFZ3A3FMO5ammT+NHrR2q45Se7qWzq5JcPLOSa6SnmFBKXDvc8Bfc9Bz1t8Ks18Pr/lBt2gpRSin+4ZTYrZ6Xx3VeKeP2IzM03Q3AG9+A2WHKbu8919/XznZeLePT3HzItLZbNX7uOG2elmV2WMWXw0b2w5K9g/5Pw04Ww+8fQLQsdBRuX08FP7pvPvEmJPPr7D/nHvxyl1+0xu6yQEqTBvROSp0HSFLMrsZXTDR3c+fP3eWZvOY9cl82f/moZmUnRZpf1sYg4uPkH8MU3jdUgd3wPfpwH2/8B2mrNrk4MERPh4o+blvLQNVn85r0y7vnF+1Q2dZpdVsgIvuB290LZbrnN3Yc8Hs1LB6u45Sfvcqa1i199Pp9vbZhNuCv4/voByFwID74Mm96G6avg/Z/C43Ph1a/IDJQgEuFy8r3b5vCL+xdQ2tDB+p+8y5ZCaZ0EgvLH/nL5+fm6oKBgdL+5bDc8tcFY13nWBt8WFmLc/R5eL6zh57tOceJsGwunJPHT++aPfu0RszSVGuF98Fno74GJ+TDnDph9OySO8ZZ84ROVTZ38ze8/5HBVKw8um8I31+cSGSYTC0ZCKXVAa53v1WuDLrjf/CfY/Tj879MQGeBZDjbR4+7nzx9W84u3T1He2MmMtFj++sZp3HrVBFzOIB1le6O9Dg4+A0dfhtqBNU8m5sOcTw2E+GRz6wtxvW4PP9h6nF/tPk3WuGi+sDybuxZkEhPhMrs0S7B2cD9xIzjDYeM23xYVAqpbuth8pIYnd5dy9lwPV2Um8OiN01mdm47DYbNt3xpPwUcvXxji42bA5CUweZnxSJ4q292ZYNeJOh7fUcLhyhbiIl3ct3gyDy6bElzXU4KQz4NbKbUO+A/ACTyptf7+lV4/6uDubIJ/nQor/o/xEFfU4+5n/+lm3i6u460T9ZTUGYsALZs6jkdvnM6108f5ZpGoYNd4Co6/BuV7oHIvdDUbn49JhclLjV3oU2YaW6olTwVXuLn1hogPK5r59e7TbCmqRWvNurzx3L9kCgumJEkb5RJ8GtxKKSdQDKwGqoD9wH1a648u93tGHdxHX4LnH4KN22HS4pH/fhvr7HVTWt9BaUMHp+raKapu5f1TjXT19RPudLA4O5kVM1NZMTOV6WlxZpdrHo/HuA+gYo+xLHDFHmOz6UEOFyRlQ+pMSM6GuAnGjj1xEyB+grGeitMP67SEsDMtXTy9p5w/fFBBa1cfYU5F3sQE8qcksXBKMvlZSaTERphdpul8HdzLgO9prdcO/PrvAbTW/3K53zPq4H71K3D0FfhGKTjt2RfTWtPXr+l299Pd109Pn4fuvn7Odbtp7uilqbOX5o5emjv7aO7o5UxrF6X1HVS3dJ0/hlKQNS6G5dNTWDEzlWXTxhEdbs8/L5/oaYfGEqgvhoYTUH/CCPfmcuNi5wUURI+D6GRjyYWopCGPRAiPgbDoIc/REBYDrgjj4YwwRvTOIQ+Ha+DhDOnWTWevm90lDRwob6agvJnCqlZ6+43535OSo5iUFE16fCTp8ZGMj49gfEIkafGRxEe6iAp3ERXmJDrcSYTLYcufJEcS3N78a58IDN14rgpYMprCrkhrzh7aSokjl3/8j/d8fvhh3/5KXxvyn9sFr9PGr7XWeDRoNFobd+x7tMbt0Xg8xnP/wKPH3Y836/OEORVJ0eGkx0eyKCuJe1MnMTU1lmlpMWSNi5EfNUciIhYmzDceQ2lttFXOVRu787SdMZ7ba6Grxfhaey3UHzN+3eODuzkHQ1w5QTmMdVvUJR6ogZBXA59jyOe46OtDQ2zIx1cMNx8F3wgCNBpYM/AA8GRoetweunr76XH3467VuM9o3B7PBSv+eoCOgccgh1IXvLVS6hNnZEa2dzoTmP0t/+eXN8F9qdP/RPQopTYBmwAmTx7F1X13N6fj8imOnM+MeHPW3/7kX/0FX7zUhyilcCjjc+e/eRS4HArn4EMpnA4HTocx9zUyzEFkmDFyiAhzEhnmJC7SRXJ0OMkx4STFhBMT7rTlqCKoKGWMrKOTYfzc4V/v6YfeDujrvPC5twP6e8Hdc+Fzf6+xbK3Hbfxez+DHA7/Wngsfg59DD3yOgY/1x5+HgXXM9cfPgy746fkKowOfTUgY23EcQNTA44KjaujtN34S7e7z4PZ4Ph4AaWMw1OsZHCTpgcHTwAN9/o/N+wp9N0HDHRbvs2NdSXC1SoQQIkSNpFXizaTe/cAMpVS2UiocuBd4dSwFCiGEGL1hWyVaa7dS6m+AbRjTAX+ttT7q98qEEEJckldTEbTWm4HNfq5FCCGEFyx8/7MQQoQmCW4hhLAYCW4hhLAYCW4hhLAYCW4hhLAYvyzrqpSqB0a7BXQK0ODDcqxAztn+Qu18Qc55pKZorVO9eaFfgnsslFIF3t49ZBdyzvYXaucLcs7+JK0SIYSwGAluIYSwmGAM7ifMLsAEcs72F2rnC3LOfhN0PW4hhBBXFowjbiGEEFdgWnArpdYppU4opU4qpT6xM7BSKkIp9dzA1/cppbICX6XveHG+f6eU+kgpdUQp9aZSaooZdfrScOc85HV3K6W0UsryMxC8OWel1KcH/q6PKqV+H+gafc2L7+3JSqldSqmDA9/f682o01eUUr9WStUppYou83WllPrJwJ/HEaXUAp8XobUO+ANjedhTwFQgHDgMzL7oNX8N/GLg43uB58yoNYDneyMQPfDxl618vt6e88Dr4oB3gL1Avtl1B+DveQZwEEga+HWa2XUH4JyfAL488PFsoMzsusd4ztcDC4Ciy3x9PbAFY2OspcA+X9dg1oh7MXBSa12qte4F/gjcftFrbgd+O/DxC8BKZd29vIY9X631Lq1158Av9wKZAa7R17z5Owb4v8C/At2BLM5PvDnnR4Cfaa2bAbTWdQGu0de8OWcNDO7plQCcCWB9Pqe1fgdousJLbgee1oa9QKJSKsOXNZgV3JfagHji5V6jtXYDrcC4gFTne96c71AbMf7HtrJhz1kpNR+YpLV+LZCF+ZE3f885QI5S6j2l1F6l1LqAVecf3pzz94D7lVJVGOv6fyUwpZlmpP/eR8yrjRT8wJsNiL3apNgivD4XpdT9QD5wg18r8r8rnrNSygH8GHgoUAUFgDd/zy6MdskKjJ+q3lVK5WmtW/xcm794c873AU9prX84sIftMwPn7PF/eabwe3aZNeKuAiYN+XUmn/zx6fxrlFIujB+xrvTjSTDz5nxRSq0CvgXcprXuCVBt/jLcOccBecBbSqkyjF7gqxa/QOnt9/UrWus+rfVp4ARGkFuVN+e8EfgTgNZ6DxCJsaaHXXn1730szApubzYgfhX4/MDHdwM79UDn34KGPd+BtsEvMULb6n1PGOactdatWusUrXWW1joLo69/m9a6wJxyfcKb7+uXMS5Eo5RKwWidlAa0St/y5pwrgJUASqlcjOCuD2iVgfUq8ODA7JKlQKvWusan72Dildn1QDHGFelvDXzunzD+8YLxl/s8cBL4AJhq9tVkP5/vDuAscGjg8arZNfv7nC967VtYfFaJl3/PCvgR8BFQCNxrds0BOOfZwHsYM04OAWvMrnmM5/sHoAbowxhdbwS+BHxpyN/xzwb+PAr98X0td04KIYTFyJ2TQghhMRLcQghhMRLcQghhMRLcQghhMRLcQghhMRLcQghhMRLcQghhMRLcQghhMf8fFOw4BBBObMEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xvals = np.linspace(0, 1)\n",
    "\n",
    "hd_dist = pm.Beta.dist(alpha=13, beta=11).logp(xvals)\n",
    "sd_dist = pm.Beta.dist(alpha=3, beta=14).logp(xvals)\n",
    "\n",
    "plt.plot(xvals, np.exp(hd_dist.eval()))\n",
    "plt.plot(xvals, np.exp(sd_dist.eval()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Standard normal from complimentary error function\n",
    "Φ = lambda x: 0.5*tt.erfc(-x/tt.sqrt(2))\n",
    "\n",
    "def power_analysis(n):\n",
    "    \n",
    "    # Allocate samples 2:1\n",
    "    n_hd = int(2*n/3)\n",
    "    n_sd = n - n_hd\n",
    "    \n",
    "    with pm.Model() as model:\n",
    "\n",
    "        # probabilities of >= 4-fold-rise taken from Beta priors constructed from phase I data\n",
    "        π_hd = pm.Beta('π_hd', 12, 10)\n",
    "        π_sd = pm.Beta('π_sd', 2, 13)\n",
    "        \n",
    "        # Draw binomial samples from probabilities\n",
    "        rng = RandomStreams(seed=20090425)\n",
    "        #x_hd = pm.Binomial('x_hd', n_hd, π_hd)\n",
    "        x_hd = rng.binomial(n=n_hd, p=π_hd)\n",
    "        #x_sd = pm.Binomial('x_sd', n_sd, π_sd)\n",
    "        x_sd = rng.binomial(n=n_sd, p=π_sd)\n",
    "\n",
    "        \n",
    "        # Empirical probabilities\n",
    "        π_hd_post = x_hd/n_hd\n",
    "        π_sd_post = x_sd/n_sd\n",
    "        \n",
    "        # Convert difference in probabilities to standard normal\n",
    "        d = pm.Deterministic('d', (π_hd_post - π_sd_post)\n",
    "                             / tt.sqrt(π_hd_post*(1-π_hd_post)/n_hd + π_sd_post*(1-π_sd_post)/n_sd))\n",
    "        \n",
    "        # P-value\n",
    "        p = 1 - Φ(d)\n",
    "        \n",
    "        # Successful detection\n",
    "        power = pm.Deterministic('power', (p<0.05).astype('int32'))\n",
    "        \n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use model to estimate power for a range of sample sizes.\n",
    "\n",
    "n=60:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [π_sd, π_hd]\n",
      "Sampling 2 chains: 100%|██████████| 6000/6000 [00:04<00:00, 1284.20draws/s]\n"
     ]
    }
   ],
   "source": [
    "iterations = 1000\n",
    "tune = 2000\n",
    "\n",
    "mod60 = power_analysis(60)\n",
    "with mod60:\n",
    "    tr = pm.sample(iterations, tune=tune, cores=2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>mc_error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>power</th>\n",
       "      <td>0.891</td>\n",
       "      <td>0.006796</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        mean  mc_error\n",
       "power  0.891  0.006796"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pm.summary(tr, varnames=['power'])[['mean', 'mc_error']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "n=48:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [π_sd, π_hd]\n",
      "Sampling 2 chains: 100%|██████████| 6000/6000 [00:05<00:00, 1148.71draws/s]\n",
      "The acceptance probability does not match the target. It is 0.8859333589363031, but should be close to 0.8. Try to increase the number of tuning steps.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>mc_error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>power</th>\n",
       "      <td>0.8605</td>\n",
       "      <td>0.007625</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         mean  mc_error\n",
       "power  0.8605  0.007625"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod48 = power_analysis(48)\n",
    "with mod48:\n",
    "    tr = pm.sample(iterations, tune=tune, cores=2)\n",
    "    \n",
    "pm.summary(tr, varnames=['power'])[['mean', 'mc_error']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "n=42:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [π_sd, π_hd]\n",
      "Sampling 2 chains: 100%|██████████| 6000/6000 [00:05<00:00, 1085.89draws/s]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>mc_error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>power</th>\n",
       "      <td>0.8485</td>\n",
       "      <td>0.008528</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         mean  mc_error\n",
       "power  0.8485  0.008528"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod42 = power_analysis(42)\n",
    "with mod42:\n",
    "    tr = pm.sample(iterations, tune=tune, cores=2)\n",
    "    \n",
    "pm.summary(tr, varnames=['power'])[['mean', 'mc_error']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "n=36:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [π_sd, π_hd]\n",
      "Sampling 2 chains: 100%|██████████| 6000/6000 [00:05<00:00, 1162.65draws/s]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>mc_error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>power</th>\n",
       "      <td>0.834</td>\n",
       "      <td>0.007513</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        mean  mc_error\n",
       "power  0.834  0.007513"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod36 = power_analysis(36)\n",
    "with mod36:\n",
    "    tr = pm.sample(iterations, tune=tune, cores=2)\n",
    "    \n",
    "pm.summary(tr, varnames=['power'])[['mean', 'mc_error']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "n=27:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [π_sd, π_hd]\n",
      "Sampling 2 chains: 100%|██████████| 6000/6000 [00:04<00:00, 1244.64draws/s]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>mc_error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>power</th>\n",
       "      <td>0.7505</td>\n",
       "      <td>0.009151</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         mean  mc_error\n",
       "power  0.7505  0.009151"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod27 = power_analysis(27)\n",
    "with mod27:\n",
    "    tr = pm.sample(iterations, tune=tune, cores=2)\n",
    "    \n",
    "pm.summary(tr, varnames=['power'])[['mean', 'mc_error']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power analysis suggests that a sample size in the low- to mid-30s would be sufficient to detect a difference between high and standard dose with 80% power in a similar population. A sample size over 60 would be required for 90% power."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "1. GiaQuinta, S., Michaels, M. G., McCullers, J. A., Wang, L., Fonnesbeck, C., O'Shea, A., et al. (2014). Randomized, double-blind comparison of standard-dose vs. high-dose trivalent inactivated influenza vaccine in pediatric solid organ transplant patients. Pediatric Transplantation, 19(2), 219–228. http://doi.org/10.1111/petr.12419\n",
    "1. Spiegelhalter, D. J., Abrams, K. R., and Myles, J. P. (2004). Bayesian approaches to Clinical Trials and Health Care Evaluation. Wiley, New York."
   ]
  }
 ],
 "metadata": {
  "_draft": {
   "nbviewer_url": "https://gist.github.com/c6c1dc8c112d785331684e2ccd2575c6"
  },
  "gist": {
   "data": {
    "description": "Vaccine Power Analysis-Copy1.ipynb",
    "public": true
   },
   "id": "c6c1dc8c112d785331684e2ccd2575c6"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}