Binomial Distributions and Evidence

In this demonstration, we consider the binomial distribution. We proved in class that the probability of exactly j successes in a series of for a series of n Bernoulli trials with success probability p is given by

	b(n,p,j) = (({{n}, {j}}))(((p^j)(1-p))^(n-j))

For a given p and n this is called the Binomial Distribution (over j ∈ 1,...,n). Mathematica has the binomial distribution built in:

BinomialDistribution[n,p]

refers to this distribution. We can see that the PDF function is exactly what we expect:

PDF[BinomialDistribution[n,p]][j]

We can plot this distribution nicely with DiscretePlot

DiscretePlot[PDF[BinomialDistribution[100,0.1]][k],{k,0,100},PlotRange->All]

Suppose that we have some experimental data which shows that we have j successes in n trials, but we don’t know anything at all about p. What can we say about p?

Definition. Given some evidence X and a model described by a parameter θ, the probability mass function for observing X given θ is called the likelihood function p(X|θ).

In this case, we know that 

	p(X|θ) = b(n,θ,j)

We can contrast this with what we’d really like to know: the probability mass function for θ given evidence X. To define this posterior probability distribution, we need one more ingredient: a probability mass function p(θ) summarizing our beliefs, if any, about θ before we saw the evidence X.

Definition. The posterior probability mass function of parameter value θ given evidence X and a prior distribution p(θ) is given by

	p(θ|X) = K(x) p(X |θ) p(θ)

where K is a normalizing constant defined to make the integral of p(θ|X) over all possible θ equal to 1.

In our example, let’s suppose our prior belief about p is simply that it is uniformly distributed between 0 and 1. (We’ll see later that this is the least knowledge we can have about p.) Thus p(θ) = 1.

To work out the constant of proportionality, we just have to integrate

Integrate[PDF[BinomialDistribution[n,x]][j],{x,0,1}]

We can now write down the posterior mass function:

μ[n_,j_,x_] := (n+1)PDF[BinomialDistribution[n,x],j]

So suppose we observe 40 successes in 100 trials. It’s intuitive that the most likely value for p is 0.4, but what do we actually know? From this trial, we know that the mass function for p looks like:

Plot[μ[100,40,x],{x,0,1},PlotRange->Full,Filling->Axis]

which is certainly centered at 0.4, but has quite a wide distribution. We can find the probability that 0.3<p< 0.5 by integrating

N@Integrate[μ[100,40,x],{x,0.3,0.5}]

We would then describe [0.3,0.5] as a “96% credible interval” for p, meaning that, given the data, the probability that p ∈ [0.3,0.5] > 0.96. This approach is called a Bayesian approach.

We now consider the data from the largest trial of the Moderna COVID vaccine.  According to https://www.modernatx.com/covid19vaccine-eua/providers/clinical-trial-data in the Phase III trial, 14073 adults were given doses of a placebo and 185 cases of COVID were observed, while 14134 adults were given doses of the vaccine and 11 cases of COVID were observed.

Assume that x is your chance of contracting COVID (during the study period). Then our two posterior distributions for x (without vaccine and with vaccine) are given by

Plot[{μ[14073,185,x],μ[14134,11,x]},{x,0,1},PlotRange->{{0,0.1},{0,1500}},Filling->Axis,PlotLegends->{"No vax","Vax"}]

We can see that the distribution of probabilities for the “Vax” cases is centered on a much smaller number and there is little overlap. From this, we can conclude that it is likely that the vaccine reduces your chances of contracting COVID.

But this is a probability class, so... how likely? We can answer this as a question about the joint distribution of two variables

	x = probability of COVID with no vax
	y = probability of COVID with vax

We wish to know the probability of the event x > y. We can find this by integrating, but we’ll need one further trick.

Theorem. The binomial distribution b(n,p,j) is well approximated by

NormalApproximation[n_,p_,j_] := PDF[NormalDistribution[n p,Sqrt[n p (1 - p)]],j]

We can check this by plotting:

ComparisonPlot[n_,p_] := Show[{Plot[NormalApproximation[n,p,j],{j,0,n},Filling->None,PlotRange->All],DiscretePlot[PDF[BinomialDistribution[n,p],j],{j,0,n}]},PlotRange->All]

ComparisonPlot[1000,0.5]

Now we can define a posterior approximation as before by multiplying by (n+1):

μa[n_,j_,p_] := (n+1)NormalApproximation[n,p,j]

We can check that this looks like our previous data by plotting

Plot[{μa[14073,185,x],μa[14134,11,x]},{x,0,1},PlotRange->{{0,0.02},{0,1500}},Filling->Axis,PlotLegends->{"No vax","Vax"}]

Now we can approximate the joint distribution by

CovidJointDistribution[pNoVax_,pVax_] := μa[14073,185,pNoVax]*μa[14134,11,pVax]

We can plot this (after some fooling around with parameters), at least over the range [0,.03]x[0,.03]

Plot3D[CovidJointDistribution[pNoVax,pVax],{pNoVax,0,0.03},{pVax,0,0.03},PlotPoints->{100,100},PlotRange->Full]

Now we can find the probability that the vaccine helps by integrating this joint probability distribution over the region where pNoVax > pVax. (These are the possibilities where the vaccine really does reduce your chance of COVID.)

NIntegrate[CovidJointDistribution[pNoVax,pVax],{pVax,0,0.03},{pNoVax,pVax,0.03}]

Now we’re going to consider a case where the outcome is less clear. This link takes you to a Canadian study about mammography and breast cancer risk. In the study, Impact of screening mammography on mortality from breast cancer before age 60 in women 40 to 49 years of age, the authors divided 50436 women into two equal sized groups. One group received regular mammograms from age 40-49 and the other group received no mammograms during the same period.

Among the 25218 women who received mammograms, 132 died of breast cancer before 60. Among the 25218 women who did not, 122 died of breast cancer before age 60. (This is so surprising I just want to say it’s not a typo-- more women who got mammograms died.) 

If we assume a Bernoulli model, the distribution of pMammogram and pNoMammogram are

Plot[{μ[25218,132,x],μ[25218,122,x]},{x,0,0.01},PlotRange->All,PlotLegends->{"Mammograms","No Mammograms"},Filling->Axis]

We can see from the overlap that it’s not clear what we can conclude from this experiment.  Let’s look at the joint distribution.

MammogramJointDistribution[pNoMammogram_,pMammogram_] := μa[25218,132,pMammogram]*μa[25218,122,pNoMammogram]

Plot3D[MammogramJointDistribution[pNoMammogram,pMammogram],{pNoMammogram,0,0.01},{pMammogram,0,0.01},PlotPoints->100,PlotRange->Full]

We can now integrate to the compute the probability of the event “pNoMammogram > pMammogram” (that is, the probability that mammograms work at all):

NIntegrate[MammogramJointDistribution[pNoMammogram,pMammogram],{pMammogram,0,0.1},{pNoMammogram,pMammogram,0.1}]

We see that the odds are against the hypothesis that mammograms work. Based on the evidence from this study, with a Bernoulli model there’s only a ~26% chance that mammograms reduce risk, and a ~74% chance that they increase risk.

Of course, there are a lot of disclaimers. It’s not true that everyone’s risk for breast cancer is the same, so the Bernoulli assumption is debatable from the start. If the study happened to enroll more high-risk participants in the mammogram group, then it’s still possible that mammograms reduced their risk (but did not reduce it below the risk of the not-mammogram group).