%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# Table from Wacholder et al
from IPython.core.display import Image
Image('pics/h0_table.png')
def ppv(pi, beta, alpha=0.05):
""" Positive predictive value from Ioannidis 2005
This is the probability that the finding is true if the paper reports a
significant result
"""
true_sig = (1 - beta) * pi
false_sig = alpha * (1 - pi)
return true_sig / (true_sig + false_sig)
Let's say there's a 50% prior probability that the finding was true (fairly safe bet) but the power is low (20%):
ppv(0.5, 0.8, 0.05)
0.8
Maybe the bet is not so safe - and the prior probability is 20% (page 336 of Button et al).
ppv(0.2, 0.8, 0.05)
0.4999999999999999
So, if the prior odds were 20%, and the power is 20%, then a significant finding has less than 50% chance of reflecting an actual true finding.
Put another way - a positive finding in a low powered study is only likely to be correct if it was reasonably likely before you did the study.
Add bias. $u$ is the proportion of non-significant results that have been declared significant as a result of bias.
def biased_ppv(pi, beta, u, alpha=0.05):
""" Positive predictive value with bias, from Ioannidis 2005
This is the probability that the finding is true if the paper reports a
significant result, in the presence of bias
"""
true_sig = (1 - beta) * pi + u * beta * pi
false_sig = alpha * (1 - pi) + u * (1 - alpha) * (1 - pi)
return true_sig / (true_sig + false_sig)
Reproduce figure 1 from Ioannidis 2005
R = np.linspace(0, 1)
p = R / (R + 1)
for u in [0.05, 0.2, 0.5, 0.8]:
plt.plot(R, biased_ppv(p, 0.2, u), label='u={0}'.format(u))
plt.legend()
<matplotlib.legend.Legend at 0x101fee710>
The figure is not exactly the same as the Ioannidis figure:
from IPython.display import Image
Image('pics/ioannidis_figure1.png')
Check : the formula is:
Image('pics/bias_formula.png')
Check that looks like our figure, not the paper figure:
def biased_ppv_paper(alpha, beta, R, u):
return (
((1 - beta) * R + u * beta * R) /
(R + alpha - beta * R + u - u * alpha + u * beta * R))
for u in [0.05, 0.2, 0.5, 0.8]:
plt.plot(R, biased_ppv_paper(0.05, 0.2, R, u), label='u={0}'.format(u))
plt.legend()
<matplotlib.legend.Legend at 0x105ac0f90>