ioannidis_2005.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:cff6f2ccecc7e0e8ced0822c29c9dc08287d3ce60b6a8397379f16a5dc41ed44"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The argument in \"Why most published research findings are false\""
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I spent some time trying to understand the argument in this paper:\n",
      "\n",
      "* Ioannidis, John PA. 2005. \u201cWhy most published research findings are false.\u201d *PLoS medicine* 2 (8): e124.\n",
      "\n",
      "These papers were useful for understanding the argument:\n",
      "\n",
      "* Goodman, Steven, and Sander Greenland. 2007. \u201cAssessing the unreliability of the medical literature: a response to \u2018why most published research findings are false.\u2019\u201d *Johns Hopkins University, Dept. of Biostatistics Working Papers*.\n",
      "* Kass, Robert E., and Adrian E. Raftery. 1995. \u201cBayes factors.\u201d *Journal of the american statistical association* 90 (430): 773\u2013795.\n",
      "* Wacholder, Sholom, Stephen Chanock, Montserrat Garcia-Closas, Nathaniel Rothman, and others. 2004. \u201cAssessing the probability that a positive report is false: an approach for molecular epidemiology studies.\u201d *Journal of the National Cancer Institute* 96 (6): 434\u2013442."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "\"The practice of science is profoundly broken\".  Discuss? - no - model and test!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The point that Ioannnides makes is:\n",
      "\n",
      "We know that the scientific process is flawed in a variety of ways.  We assume that these flaws do not have a large effect on the outcome.  But, if we model some of the flaws, we see that their effect can be catastrophic, in the sense that a large proportion of scientific findings are likely to be wrong.\n",
      "\n",
      "We scientists commit ourselves to rational thinking.  In this case, rational thinking is asking, \"how likely is it that we are getting the answers wrong\"?.  We have to ask this question in a rational way.  This is what Ioannidis sets out to do in this paper."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Different ways of phrasing the argument"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The basis of Ioannidis' argument comes from [Wacholder et al 2004](http://jnci.oxfordjournals.org/content/96/6/434.long) (see appendix table 1).  [Goodman and Greenland 2007](http://www.plosmedicine.org/article/info%3Adoi%2F10.1371%2Fjournal.pmed.0040168) explain Ioannidis in terms of Bayes theorem.\n",
      "\n",
      "Both Ioannidis and Goodman & Greenland use odds ratios rather than probability values in their exposition.  I found it easier to think in terms of probabilities."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Some terms"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We've done an experiment, and we have conducted a statistical test:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "* $H_A$ - alternative hypothesis\n",
      "* $H_0$ - null hypothesis\n",
      "* $\\alpha$ : false positive rate - probability for test to reject $H_0$ when $H_0$ is true ($H_A$ is false)\n",
      "* $\\beta$ : false negative rate - probability for test to accept $H_0$ when $H_A$ is true ($H_0$ is false)\n",
      "* $1 - \\beta$ : power - probability we will reject $H_0$ if $H_A$ is true ($H_0$ is false)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's say that the test can either be \"significant\" (test gives $p \\le \\alpha$) or \"not significant\" ($p > \\alpha$).  Denote \"test is significant\" by $T_S$, \"test not significant\" by $T_N$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Before we did the experiment there were two possibilities - $H_A$ is true, or $H_0$ is true.  After we have four possibilities: \n",
      "\n",
      "* $H_A \\land T_S$ : $H_A$ is true, test is significant; \n",
      "* $H_A \\land T_N$ : $H_A$ is true, test is not significant;\n",
      "* $H_0 \\land T_S$ : $H_0$ is true ($H_A$ is false) - test is significant;\n",
      "* $H_0 \\land T_N$ : $H_0$ is true ($H_A$ is false) - test is not significant."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "What does a \"significant\" statistical test result tell us?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this section we work up slowly to Ioannidis table 1."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First we need to load functions for symbolic mathematics from the Sympy library:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sympy import symbols, Eq, solve, simplify, lambdify, init_printing, latex\n",
      "init_printing(use_latex=True, order='old')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Define a function to display the probability tables:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.display import HTML\n",
      "# Code to make HTML for a probability table\n",
      "def association_table(assocs, title):\n",
      "    latexed = {'title': title}\n",
      "    for key, value in assocs.items():\n",
      "        latexed[key] = latex(value)\n",
      "    latexed['s_total'] = latex(assocs['t_s'] + assocs['f_s'])\n",
      "    latexed['ns_total'] = latex(assocs['t_ns'] + assocs['f_ns'])\n",
      "    return \"\"\"<h3>{title}</h3>\n",
      "              <TABLE><TR><TH><TH>$T_S$<TH>$T_N$\n",
      "              <TR><TH>$H_A$<TD>${t_s}$<TD>${t_ns}$\n",
      "              <TR><TH>$H_0$<TD>${f_s}$<TD>${f_ns}$\n",
      "              <TR><TH>Total<TD>${s_total}$<TD>${ns_total}$\n",
      "              </TABLE>\"\"\".format(**latexed)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If we do not take any prior probabilities into account, then we have the following probabilities:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sympy.abc import alpha, beta # get alpha, beta symbolic variables\n",
      "assoc = dict(t_s = 1 - beta, # H_A true, test significant = true positives\n",
      "             t_ns = beta, # true, not significant = false negatives\n",
      "             f_s = alpha, # false, significant = false positives\n",
      "             f_ns = 1 - alpha) # false, not sigificant = true negatives\n",
      "HTML(association_table(assoc, 'Not considering prior'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$\\newcommand{Frac}[2]{\\frac{\\displaystyle #1}{\\displaystyle #2}}$\n",
      "\n",
      "Some new terms:\n",
      "\n",
      "* $Pr(H_A)$ - prior probability of $H_A$ - probability of $H_A$ before the experiment was conducted.\n",
      "* $Pr(H_0)$ - prior probability of $H_0$ = $1 - Pr(H_A)$ - probability of null hypothesis before experiment conducted\n",
      "\n",
      "We are interested in updating the probability of $H_A$ and $H_0$ as a result of a test on some collected data.  This updated probability is $Pr(H_A | T)$ - the probability of $H_A$ given the test  result $T$. $Pr(H_A | T)$ is called the *posterior* probability because it is the probability after the test result.\n",
      "\n",
      "The test result $T$ is assumed to have arisen under either $H_A$ or $H_0$.\n",
      "\n",
      "$Pr(T) = Pr(T | H_A) Pr(H_A) + Pr(T | H_0) Pr(H_0)$\n",
      "\n",
      "Also the probability of a *signficant* result of the test $T_S$ is from the same formula:\n",
      "\n",
      "$Pr(T_S) = Pr(T_S | H_A) Pr(H_A) + Pr(T_S | H_0) Pr(H_0)$\n",
      "\n",
      "(From Kass & Rafferty 1995)\n",
      "\n",
      "Remembering [Bayes theorem](http://en.wikipedia.org/wiki/Bayes'_theorem#Derivation):\n",
      "\n",
      "$P(A | B) = \\Frac{P(B | A) P(A)}{P(B)}$\n",
      "\n",
      "Bayes theorem gives:\n",
      "\n",
      "$P(H_A | T) = \\Frac{Pr(T | H_A) Pr(H_A)}{Pr(T)} = \\Frac{Pr(T | H_A) Pr(H_A)}{Pr(T | H_A) Pr(H_A) + Pr(T | H_0) Pr(H_0)}$\n",
      "\n",
      "Consider only the test result $T_S$ (the test is signficant).  What is the posterior probability of $H_A$ given that the test is signficant?\n",
      "\n",
      "$P(H_A | T_S) = \\Frac{Pr(T_S | H_A) Pr(H_A)}{Pr(T_S | H_A) Pr(H_A) + Pr(T_S | H_0) Pr(H_0)}$\n",
      "\n",
      "We have $Pr(T_S | H_A)$, $Pr(T_S | H_0)$ from the first column of the table above. Substituting into the equation:\n",
      "\n",
      "$P(H_A | T_S) = \\Frac{(1 - \\beta) Pr(H_A)}{(1 - \\beta) Pr(H_A) + \\alpha Pr(H_0)}$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To make this a little less cluttered, define:\n",
      "\n",
      "$\\pi := Pr(H_A)$\n",
      "\n",
      "So\n",
      "\n",
      "$1 - \\pi = Pr(H_0)$\n",
      "\n",
      "and:\n",
      "\n",
      "$P(H_A | T_S) = \\Frac{(1 - \\beta) \\pi}{(1 - \\beta) \\pi + \\alpha (1 - \\pi)}$\n",
      "\n",
      "Let's put that formula into Sympy for later:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sympy.abc import pi # get symbolic variable pi\n",
      "post_prob = (1 - beta) * pi / ((1 - beta) * pi + alpha * (1 - pi))\n",
      "post_prob"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A table shows the new probabilities that take the prior into account:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "assoc = dict(t_s = pi * (1 - beta),\n",
      "             t_ns = pi * beta,\n",
      "             f_s = (1 - pi) * alpha,\n",
      "             f_ns = (1 - pi) * (1 - alpha))\n",
      "HTML(association_table(assoc, r'Considering prior $\\pi := Pr(H_A)$'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This table is equivalent to Ioannidis table 1. The first column of the table gives the probabilities in the case we're interested in, of $T_S$.  The posterior probability is the first row, first column - $Pr(T_S | H_A)$, divided by the total row, first column - $Pr(T_S)$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ioannidis uses \"positive predictive value\" (PPV) for the posterior probability $P(H_A | T_S)$.  Goodman and Greenland complain, reasonably enough, that \"positive predictive value\" is a confusing new term for a standard concept."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ioannidis also prefers his equations in terms of $R$ - the *prior odds ratio*.  $R := \\Frac{Pr(H_A)}{Pr(H_0)}$. Also (from $\\pi := Pr(H_A)$ and $Pr(H_0) = 1 - Pr(H_A)$): $R = \\Frac{\\pi}{1 - \\pi}$. \n",
      "\n",
      "Ioannidis' formula for PPV is $\\Frac{(1 - \\beta) R}{R - \\beta R + \\alpha}$.  This is the same as our formula above, only rephrased in terms of $R$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "R = pi / (1 - pi)\n",
      "ppv = (1 - beta) * R / (R - beta * R + alpha)\n",
      "# Is this the same as our formula above?\n",
      "simplify(ppv - post_prob) == 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The posterior probability is our estimate of whether $H_A$ is true, given our prior knowledge $Pr(H_A) = \\pi$ combined with the new information from the test result."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "What is a finding that is likely to be true?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A finding that is likely to be true is one for which the posterior probability $Pr(H_A | T_S) > 0.5$.  That is, the likelihood of the tested hypothesis, given the reported significant test result, is greater than $0.5$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Whether a finding is likely to be true depends on the power of the experiment"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Assume that $\\alpha = 0.05$ (standard significance threshold for null hypothesis test)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First we load the plotting and numerical libraries."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's have a look at the posterior probability as a function of power and prior probability:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make symbolic expression into a function we can evaluate\n",
      "post_prob_func = lambdify((alpha, beta, pi), post_prob)\n",
      "# Make set of pi values for x axis of plot\n",
      "pi_vals = np.linspace(0, 1, 100)\n",
      "for power in (0.8, 0.5, 0.2):\n",
      "    beta_val = 1 - power\n",
      "    plt.plot(pi_vals, post_prob_func(0.05, beta_val, pi_vals), label='power={0}'.format(power))\n",
      "plt.xlabel('Prior probability $Pr(H_A)$')\n",
      "plt.ylabel('Posterior probability $Pr(H_A | T_S)$')\n",
      "plt.legend()\n",
      "plt.title(\"Posterior probability for different priors and power levels\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The posterior probability depends on the power.  If the power is low and $H_A$ is true, the likelihood of getting a significant test result is small.  Assuming $\\pi = Pr(H_A) = 0.5$, our posterior probability is given by $\\Frac{(1 - \\beta)}{(1 - \\beta) + \\alpha}$.  As the chance of finding a true positive $= 1-\\beta$ drops towards the chance of finding a false negative $= \\alpha$, our confidence in the truth of the significant result must drop too.  \n",
      "\n",
      "The posterior likelihood also depends on the prior. When the prior $Pr(H_A)$ drops then we become more wary of the (apriori more unlikely) true positive compared to the (apriori more likely) false positive.\n",
      "\n",
      "As you can see from the figure.  When power is 0.2, and the prior probability is less than around 0.2, then even if there is a significant test result, the null is still more likely than the $H_A$ (posterior < 0.5)."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Quantifying the effect of bias"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Working scientists know that working scientists have many sources of bias in data collection and analysis.\n",
      "\n",
      "We tend to assume that the effect of this bias is relatively minor.  Is this true?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ioannidis quantifies bias with a parameter $u$.  $u$ is the proportion of not-significant findings that become significant as a result of bias.  Put another way, the effect of bias is the result of taking the second column in the probablity table above (the not-significant findings) and multiplying by $u$. We add this effect to the first column (significant findings) and subtract from the second column (not-significant findings). Before applying the priors, this looks like:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "u = symbols('u')\n",
      "bias_assoc_noprior = dict(t_s = (1 - beta) + u * beta,\n",
      "                          t_ns = beta - u * beta,\n",
      "                          f_s = alpha + u * (1 - alpha),\n",
      "                          f_ns = (1 - alpha) - u * (1 - alpha))\n",
      "HTML(association_table(bias_assoc_noprior, r'Bias $u$ not considering prior'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "After applying the prior:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bias_assoc = bias_assoc_noprior.copy()\n",
      "bias_assoc['t_s'] *= pi\n",
      "bias_assoc['t_ns'] *= pi\n",
      "bias_assoc['f_s'] *= 1 - pi\n",
      "bias_assoc['f_ns'] *= 1 - pi\n",
      "HTML(association_table(bias_assoc, r'Bias $u$ with prior $\\pi := Pr(H_A)$'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The first cell in the table is $Pr(T_S | H_A) Pr(H_A)$.  The total for the first column gives $Pr(T_S)$.  Therefore the posterior probability $Pr(H_A | T_S)$ is:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "post_prob_bias = bias_assoc['t_s'] / (bias_assoc['t_s'] + bias_assoc['f_s'])\n",
      "post_prob_bias"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Same as Ioannidis formulation?\n",
      "# This from Ioannidies 2005:\n",
      "ppv_bias = ((1 - beta) * R + u * beta * R) / (R + alpha - beta * R + u - u * alpha + u * beta * R)\n",
      "# Is this the same as our formula above?\n",
      "simplify(ppv_bias - post_prob_bias) == 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "What effect does bias have on the posterior probabilities?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Formula as a function we can evaluate\n",
      "post_prob_bias_func = lambdify((alpha, beta, pi, u), post_prob_bias)\n",
      "pi_vals = np.linspace(0, 1, 100)\n",
      "fig, axes = plt.subplots(3, 1, figsize=(8, 16))\n",
      "for i, power in enumerate((0.8, 0.5, 0.2)):\n",
      "    beta_val = 1 - power\n",
      "    for bias in (0.05, 0.2, 0.5, 0.8):\n",
      "        pp_vals = post_prob_bias_func(0.05, beta_val, pi_vals, bias)\n",
      "        axes[i].plot(pi_vals, pp_vals, label='u={0}'.format(bias))\n",
      "        axes[i].set_ylabel('Posterior probability $Pr(H_A | T_S)$')\n",
      "    axes[i].set_title('Power = {0}'.format(power))\n",
      "    axes[-1].set_xlabel('Prior probability $Pr(H_A)$')\n",
      "axes[-1].legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As we'd expect, as bias increases to 1, the result of the experiment has less and less effect on our posterior estimate.  If the analysis was entirely biased, then our posterior estimate is unchanged from the prior (diagonal line on the graph)."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The effect of multiple studies"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ioannidis makes the point that when a field is particularly fashionable, there may be many research groups working on the same question.\n",
      "\n",
      "Given publication bias for positive findings, it is possible that only positive research findings will be published.  If $n$ research groups have done the same experiment, then the probability that *all* the $n$ studies will be not significant, given $H_A$ is true, is $\\beta^n$. Conversly the probability that there is at least one positive finding in the $n$ tests is $1 - \\beta^n$.  Similarly the probability that all $n$ studies will be not significant, given $H_0$, is $(1 - \\alpha)^n$.  The probability of at least one false positive is $1 - (1 - \\alpha)^n$.\n",
      "\n",
      "The probability table (without the priors) is:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = symbols('n')\n",
      "multi_assoc_noprior = dict(t_s = (1 - beta ** n),\n",
      "                          t_ns = beta ** n,\n",
      "                          f_s = 1 - (1 - alpha) ** n,\n",
      "                          f_ns = (1 - alpha) ** n)\n",
      "HTML(association_table(multi_assoc_noprior,\n",
      "                       r'$n$ replications with positive publication bias, not considering prior'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Considering the prior:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "multi_assoc = multi_assoc_noprior.copy()\n",
      "multi_assoc['t_s'] *= pi\n",
      "multi_assoc['t_ns'] *= pi\n",
      "multi_assoc['f_s'] *= 1 - pi\n",
      "multi_assoc['f_ns'] *= 1 - pi\n",
      "HTML(association_table(multi_assoc,\n",
      "                       r'$n$ replications with positive publication bias, with prior $\\pi$'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Giving posterior probability of:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "post_prob_multi = multi_assoc['t_s'] / (multi_assoc['t_s'] + multi_assoc['f_s'])\n",
      "post_prob_multi"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Formula from Ioannidies 2005:\n",
      "ppv_multi = R * (1 - beta ** n) / (R + 1 - (1 - alpha) ** n - R * beta ** n)\n",
      "# Is this the same as our formula above?\n",
      "simplify(ppv_multi - post_prob_multi) == 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Formula as a function we can evaluate\n",
      "post_prob_multi_func = lambdify((alpha, beta, pi, n), post_prob_multi)\n",
      "pi_vals = np.linspace(0, 1, 100)\n",
      "fig, axes = plt.subplots(3, 1, figsize=(8, 16))\n",
      "for i, power in enumerate((0.8, 0.5, 0.2)):\n",
      "    beta_val = 1 - power\n",
      "    for n_studies in (1, 5, 10, 50):\n",
      "        pp_vals = post_prob_multi_func(0.05, beta_val, pi_vals, n_studies)\n",
      "        axes[i].plot(pi_vals, pp_vals, label='n={0}'.format(n_studies))\n",
      "        axes[i].set_ylabel('Posterior probability $Pr(H_A | T_S)$')\n",
      "    axes[i].set_title('Power = {0}'.format(power))\n",
      "axes[-1].set_xlabel('Prior probability $Pr(H_A)$')\n",
      "axes[-1].legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Putting it together"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Considering analysis bias and positive publication bias together:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "multi_bias_assoc_noprior = dict(t_s = (1 - beta ** n) + u * beta ** n,\n",
      "                                t_ns = beta ** n - u * beta ** n,\n",
      "                                f_s = 1 - (1 - alpha) ** n + u * (1 - alpha) ** n,\n",
      "                                f_ns = (1 - alpha) ** n - u * (1 - alpha)**n)\n",
      "HTML(association_table(multi_bias_assoc_noprior,\n",
      "                       r'$n$ replications, $u$ analysis bias, not considering prior'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "multi_bias_assoc = multi_bias_assoc_noprior.copy()\n",
      "multi_bias_assoc['t_s'] *= pi\n",
      "multi_bias_assoc['t_ns'] *= pi\n",
      "multi_bias_assoc['f_s'] *= 1 - pi\n",
      "multi_bias_assoc['f_ns'] *= 1 - pi\n",
      "HTML(association_table(multi_bias_assoc,\n",
      "                       r'$n$ replications, $u$ analysis bias, with prior'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "post_prob_multi_bias = multi_bias_assoc['t_s'] / (multi_bias_assoc['t_s'] + multi_bias_assoc['f_s'])\n",
      "post_prob_multi_bias"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's assume that two groups are doing more or less the same study, and only the positive study publishes ($n = 2$).  There is an analysis bias of 10% ($u= 0.1$).  We take the power from the Button et al estimate for neuroimaging studies = 0.08.  Therefore $\\beta = 1 - 0.08 = 0.92$:\n",
      "\n",
      "* Button, Katherine S., John PA Ioannidis, Claire Mokrysz, Brian A. Nosek, Jonathan Flint, Emma SJ Robinson, and Marcus R. Munaf\u00f2. 2013. \u201cPower failure: why small sample size undermines the reliability of neuroscience.\u201d *Nature Reviews Neuroscience*."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pp_mb_func = lambdify((alpha, beta, pi, u, n), post_prob_multi_bias)\n",
      "pp_vals_nobias = pp_mb_func(0.05, 0.92, pi_vals, 0, 1)\n",
      "pp_vals_bias = pp_mb_func(0.05, 0.92, pi_vals, 0.1, 2)\n",
      "plt.plot(pi_vals, pp_vals_nobias, label='no analysis or publication bias')\n",
      "plt.plot(pi_vals, pp_vals_bias, label='with analysis and publication bias')\n",
      "plt.plot(pi_vals, pi_vals, 'r:', label='$T_S$ not informative')\n",
      "plt.ylabel('Posterior probability $Pr(H_A | T_S)$')\n",
      "plt.xlabel('Prior probability $Pr(H_A)$')\n",
      "plt.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This graph tells us that, for a study with average power in neuroimaging, with some mild analysis bias and positive publication bias, the significant finding $T_S$ does not change our posterior very much from our prior.\n",
      "\n",
      "If we do some study with an hypothesis that is suitably unlikely apriori - say $Pr(H_A) = 0.25$ - then our posterior probability is:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pp_mb_func(0.05, 0.92, 0.25, 0.1, 2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "What if the result was significant at $p < 0.01$?:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pp_mb_func(0.01, 0.92, 0.25, 0.1, 2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So, even if our result is significant at $p < 0.01$, the probability that $H_A$ is correct is still less than $0.5$."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}