## Article in Quillette + Technical Details

Followers of my blog who haven’t already seen it might enjoy this opinion piece I wrote for Quillette magazine. I have wanted to write something like that for a while, but never thought I’d get around to it. I’m glad I proved myself wrong!

Here are the technical details of the calculations I did to get the results in the article. If you’re comfortable with an analytical derivation, see this one by Jared Tobin.

For both the Bayesian and frequentist calculations I used a binomial conditional prior for the data (aka likelihood) for the number of recoveries $x$:

$p(x | \theta) \sim \textnormal{Binomial}(100, \theta)$.

For the Bayesian analysis the prior for $\theta$ was a 50/50 mixture of a uniform prior from 0 to 1, and a delta function at $\theta=0.7$ (the old drug’s effectiveness):

$p(\theta) = \frac{1}{2}\delta(\theta - 0.7) + \frac{1}{2}$

where $\theta \in [0, 1]$. Is this prior debatable? Yes, just like the prior for the data. The conclusion of an argument can change if you change the premises.

Here’s R code to calculate the posterior:

# Theta values, prior, and likelihood function
# theta[7000] is basically 0.7
theta = seq(0, 1, length.out=10001)
likelihood = dbinom(83, prob=theta, size=100)

# Prior (discrete approx)
prior = rep(0.5/10000, 10001)
prior[7000] = 0.5

# Posterior
posterior = prior*likelihood/sum(prior*likelihood)

# Probability new drug is worse, same, better
print(sum(posterior[1:6999]))
print(posterior[7000])
print(sum(posterior[7001:10001]))

Here is R code for the (two-sided) p-value:

# Possible data sets
x = seq(0, 100)

# p(x | H0)
p = dbinom(x, prob=0.7, size=100)

# P(x <= 57 | H0)
prob1 = sum(p[x <= 57])

# P(x >= 83 | H0)
prob2 = sum(p[x >= 83])

# Print p-value
print(prob1 + prob2)