Big Data, Plainly Spoken (aka Numbers Rule Your World)

Pfizer modeling in Python

Here is code for those who want to play around with the Bayesian analysis used in the Pfizer vaccine trial study (Related blog post here). I got a request for Python code so here it is. The charts I published the other day were made in R.

# Python code for Pfizer vaccine trial analysis

import pandas as pd
import numpy as np
from scipy.stats import beta
import matplotlib.pyplot as plt


# functions

def ve2vsc (ve):
	if isinstance(ve, pd.Series) is False:
		ve = pd.Series(ve)
	return ((1-ve)/(2-ve)).round(3)

def vsc2ve (vsc):
        if isinstance(vsc, pd.Series) is False:
                vsc = pd.Series(vsc)
	return ((1-2*vsc)/(1-vsc)).round(3)

ve = pd.Series(np.linspace(0,1,11))
vsc = ve2vsc(ve)


# line chart showing VSC and VE relationship

plt.plot(ve, vsc, color="darkgray", lw=3)
plt.title("Relationship between Vaccine Efficacy and Vaccine's Share of Cases")
plt.xlabel("ve")
plt.ylabel("vsc")

ax = plt.gca()
ax.spines['right'].set_color("none")
ax.spines['top'].set_color("none")
ax.spines['bottom'].set_position(('data',0))
ax.spines['left'].set_position(('data',0))

xlabels = (ve*100).astype('int').astype('str')+"%"
xlabels[list(range(1,11,2))] = ""
plt.xticks(ve, xlabels) 
plt.yticks(pd.Series(np.linspace(0,0.5,6)), (pd.Series(np.linspace(0,50,6))).astype('int').astype("str")+"%")

ax.hlines(ve2vsc(0.5), 0, 0.5, color="brown", ls="dotted")
ax.vlines(0.5, 0, ve2vsc(0.5), color="brown", ls="dotted")
ax.hlines(ve2vsc(0.9), 0, 0.9, color="brown", ls="dotted")
ax.vlines(0.9, 0, ve2vsc(0.9), color="brown", ls="dotted")

plt.text(0.02,ve2vsc(0.5)[0]+0.02,(ve2vsc(0.5)[0]*100).astype("int").astype("str")+"%", color="brown", fontsize=8)
plt.text(0.02,ve2vsc(0.9)[0]+0.02,(ve2vsc(0.9)[0]*100).astype("int").astype("str")+"%", color="brown", fontsize=8)

plt.show()

# plotting the posterior

a1 = 0.700102	        # From protocol
b1 = 1			# From protocol
data_x = 8		# Number of cases in vaccine arm
data_n = 94		# Total number of cases

def a2 (a1, x):
	return a1 + x

def b2 (b1, n, x):
	return b1 + n - x

def vsc_posterior(samp, data_x):
	return pd.Series(beta.pdf(samp, a2(a1,data_x), b2(b1,data_n, data_x)))

vsc_sample = pd.Series(np.linspace(0,0.5,72))

# single posterior curve vs vsc

plt.plot(vsc_sample, vsc_posterior(vsc_sample, 8), lw=3)
plt.title("Predicted Prob. of Vaccine's Share of Cases")
plt.xlabel("VSC")
plt.ylabel("Probability")

ax = plt.gca()
ax.spines['right'].set_color("none")
ax.spines['top'].set_color("none")
ax.spines['bottom'].set_position(('data',0))
ax.spines['left'].set_position(('data',0))

plt.xticks(pd.Series(np.linspace(0,0.5,6)), (pd.Series(np.linspace(0,50,6))).astype('int').astype('str')+"%") 
plt.yticks(pd.Series(range(0,21,4)), [""]*6)

plt.show()

# posterior vs ve

plt.plot(vsc2ve(vsc_sample), vsc_posterior(vsc_sample, 8), lw=3)

xlabels = (pd.Series(range(0,110,10))).astype("str")+"%"
xlabels[list(range(1,11,2))] = ""
xlabels

plt.xticks(pd.Series(np.linspace(0,1.0,11)), xlabels) 
plt.yticks(pd.Series(range(0,21,4)), [""]*6)
plt.title("Predicted Prob. of Vaccine Efficacy")
plt.xlabel("ve")
plt.ylabel("probability")

plt.show()


# three posterior curves vs ve

plt.plot(vsc2ve(vsc_sample), vsc_posterior(vsc_sample, 3), color="blue", lw=1.5, linestyle="dashed", alpha=0.4)
plt.plot(vsc2ve(vsc_sample), vsc_posterior(vsc_sample, 47), color="purple", lw=1.5, linestyle="dashed", alpha=0.4)
plt.plot(vsc2ve(vsc_sample), vsc_posterior(vsc_sample, 8), 
	color="darkgray", lw=3)

ax = plt.gca()
ax.spines['right'].set_color("none")
ax.spines['top'].set_color("none")
ax.spines['bottom'].set_position(('data',0))
ax.spines['left'].set_position(('data',0))

xlabels = (pd.Series(np.linspace(0,100,11))).astype("int").astype("str")+"%"
xlabels[list(range(1,10,2))] = ""
xlabels

plt.xticks(pd.Series(np.linspace(0,1.0,11)), xlabels) 
plt.yticks(pd.Series(range(0,21,4)), [""]*6)
plt.title("Predicted Prob. of Vaccine Efficacy")
plt.xlabel("ve")
plt.ylabel("probability")

plt.text(0.02,vsc_posterior(vsc_sample, 47).max().round(2)+0.2,"47-47", color="purple", fontsize=8)
plt.text(0.02,vsc_posterior(vsc_sample, 47).max().round(2)+1,"No Effect", color="purple", fontsize=8)
plt.text(0.88,vsc_posterior(vsc_sample, 3).max().round(2)-0.2,"3-91", color="blue", fontsize=8)
plt.text(0.88,vsc_posterior(vsc_sample, 3).max().round(2)+0.5,"BEST", color="blue", fontsize=8)
plt.text(0.82,vsc_posterior(vsc_sample, 8).max().round(2)+0.2,"8-86", color="darkgray", fontsize=8)
plt.text(0.82,vsc_posterior(vsc_sample, 8).max().round(2)+1,"DATA", color="darkgray", fontsize=8, weight="bold")

plt.show()


# prior and posterior vs ve

vsc_prior = pd.Series(beta.pdf(vsc_sample, a1, b1))

plt.plot(vsc2ve(vsc_sample), vsc_prior, color="magenta", lw=1.5, ls="dashed", alpha=0.5)
plt.plot(vsc2ve(vsc_sample), vsc_posterior(vsc_sample, 8), 
	color="darkgray", lw=3)

ax = plt.gca()
ax.spines['right'].set_color("none")
ax.spines['top'].set_color("none")
ax.spines['bottom'].set_position(('data',0))
ax.spines['left'].set_position(('data',0))

xlabels = (pd.Series(np.linspace(0,100,11))).astype("int").astype("str")+"%"
xlabels[list(range(1,10,2))] = ""
xlabels

plt.xticks(pd.Series(np.linspace(0,1.0,11)), xlabels) 
plt.yticks(pd.Series(range(0,21,4)), [""]*6)
plt.title("Prior and Posterior Prob. of Vaccine Efficacy")
plt.xlabel("ve")
plt.ylabel("probability")

plt.text(0.15, 2,"prior", color="magenta", fontsize=8)
plt.text(0.76,12,"posterior", color="darkgray", fontsize=8)

plt.show()


# compute probability VE > 30%, 70% and 90% given 8-86 split of cases

beta.cdf(ve2vsc(0.3), a2(a1,8), b2(b1,data_n, 8))[0].round(3)
beta.cdf(ve2vsc(0.7), a2(a1,8), b2(b1,data_n, 8))[0].round(3)
beta.cdf(ve2vsc(0.9), a2(a1,8), b2(b1,data_n, 8))[0].round(3)

Good vaccine news without glorifying faux precision

The vaccine trials are slipping into parody with the way the scientific results are being dangled over a hungry press. It's clear that the news is great but let's not pump irrational exuberance over magical precision.

The most responsible thing to say at this point is that we have multiple highly promising vaccine candidates that have self-reported short-term efficiacy up to about 90 percent; if a vast majority of the world gets vaccinated, the world may successfully defeat the coronavirus.

***

When reading vaccine trial results, bear in mind the following:

• No clinical trial can yield an efficacy estimate as precise as 90% or 95% or 94.5%. When a statistician say 90%, we mean 90% around a margin of error. To have highly precise estimates, we must have tiny margins of error. (See also end note.)
  
  
• The precision of a statistical estimate depends on the sample size. The larger the sample size, the more precise is the estimate. In accelerated vaccine trials, decisions are being made based on only 100-150 cases, so the margins of error are not tiny.
  
  
• The difference between 90% and 95% is meaningless. The difference between 94.5% and 95% is especially pointless; it may win PR points but it suggests lack of seriousness to statisticians.
  
  
• Participants in any clinical trial do not constitute a random selection from the entire population - there are self-selection bias and exclusion criteria. The trial results are extrapolated to the entire population, and that's another reason the efficacy estimates are imprecise.
  
  
• Subgroup analysis is shaky at this stage. When filtered to subgroups, the sample sizes get even smaller, and the error bars are wider around the reported number. For example, Pfizer  stated efficacy for people over 65 was 94%. This number is essentially the same as the average efficacy for all age groups. But the error bar for the over-65 people is much wider than the error bar for the average participant.
  
  
• Press releases may highlight particular effectiveness for some subgroup; remember that this also means it is less effective for some other subgroup. AstraZeneca, for example, disclosed 90% efficacy for the 1.5-dose treatment. SInce the average efficacy was 70%, we know that the efficacy for the 2-dose treatment must be lower than 70%.

All three vaccines - Pfizer, Moderna, AstraZeneca - look promising from the short-term efficacy angle (even 60% isn't a bad result). It's hard to say more until the scientific reports are published. Efficacy or other health metrics are not the only consideration; cost, ease of transportation and storage, the schedules of availability, quantities, etc. - very practical matters - also must be evaluated.

 

***

The AstraZeneca result that just came out contains a mystery. In the design protocol, there isn't a 1.5-dose treatment so where did that come from? What the researchers appear to be saying is that an execution error - affecting a quarter of participants receiving the vaccine shots - resulted in a lower-than-expected first dose so the vaccine arm were split into two subgroups for analysis.

Most reporters have not noticed this abberation, or if they did, they didn't think it is material. This accident is immaterial only if the execution error affected a random subset of participants. In a clinical trial, we use randomization to ensure the two arms differ only by whether participants receive the vaccine or the placebo. If a third arm is added, it too must be randomly assigned. In the AstraZeneca/Oxford trial, we know the vaccine vs. placebo assignment is explicitly randomized but not so the split between 1.5 and 2 doses within the vaccine arm.

Execution errors sometimes have systematic causes, and may not be random. It's always risky to assume random, as I explained in this post.

What will happen is the researchers will give a story explaining the execution error in the scientific report, and argue why they can assume the errors were randomly scattered among the vaccine arm. This will be hard to prove. One common analysis is to show that the demographic and prior health characteristics of the group that received 1.5 doses are statistically identical to those of the group that received 2 doses. This is necessary but not sufficient; nevertheless, it will be considered good enough.

For those of us who run lots of testing, our worst nightmare is to have execution errors. What happens now is you have to trust the people, instead of trusting the (pre-registered) process. You have to trust our explanations about the errors, and the analyses justifying the assumption of random error.

I'm also curious how cutting the sample size by 75% has not affected the statistical significance of the result. That will come out in the scientific report.

That said, the average efficacy of about 70 percent is still a good result. I don't want you to think AstraZeneca's vaccine is bad. In fact, it has several advantages in terms of cost, storage, and so on.

 

P.S. The relationship between uncertainty and precision may not feel intuitive. The following probability curve is reproduced from my prior post. This curve shows the probability of vaccine efficacy ("posterior", in the sense of after looking at the data from the clinical trial).

The peak of this curve shows that the most likely value of vaccine efficacy is about 91% based on 8 cases in the vaccine group among a total of 94 cases.

The curve also shows we believe the vaccine efficacy to have a value between (roughly) 70% and 100%. Since the total probability of all possible values is 1, the total area under the curve between 70% and 100% should be exactly 1.

If you draw a vertical line under the peak of 91%, the area of that line is zero (height x 0 = 0). In other words, the probability of a precise estimate is zero. A statement such as efficacy is 91% has zero probability. But we can make other statements, e.g. efficacy is between 90% and 92%: the probability of this is the area under the curve between the values of 90% and 92%. This has positive probability.

The wider the margin of error, the larger the area under the curve, the higher the probability.

 

More notes on Pfizer's Bayesian model

I kept the previous post about Pfizer's Bayesian analysis of the vaccine trial data to a reasonable length. And then, I found some loose ends.

Bayesian analysis has the benefit of flexibility

Bayesians like to address the question head on. If our question is about vaccine efficacy (VE), we directly model that, resulting in a probability curve for VE. When we have such a curve (which I reprinted below), we have the ability to answer all kinds of questions.

In the last post, I told you that the Pfizer protocol disclosed that the FDA requires above 30% efficacy - not 50% as has been reported. That means we want at least 95% of the probability mass situated to the right of the VE = 30% point on the curve. With the results from that trial, this is almost certainly true.

What if the FDA has set 70% efficacy as the standard? We can use the same probability curve. Now we look for the probabilities to the right of 70%. This is still close to 100%.

 

Bayesian updating

I mentioned but never addressed one issue in the last post - why did the statisticians switch from the more intuitive vaccine efficacy (VE) to vaccine's share of cases (VSC)?

Here's one way to think about Bayesian models. We first set down our "prior" expectation, i.e. the pink line shown in the chart below (also reprinted from the last post):

This "prior" represents our subjective belief - prior to running the vaccine trial - of what the vaccine efficacy might be. The pink line is relatively flat across all values of VE with a slight lift towards 100%, indicating these researchers believe the most likely value of efficacy is close to 100% (!)

We can think of a Bayesian model as averaging a set of simulations. (This should remind you of the election forecasting models, which are also Bayesian.) For one such simulation, we select randomly a number from the pink curve (the "prior"). Given the flatness of the curve, pretty much any value of VE is equally likely, with a slight tip towards numbers close to 100%. Let's say we pick 75%.

In this particular simulation, the VE is now assumed to be 75%. The statistician now asks: if the vaccine's efficacy were 75%, and if I ran the trial to 94 cases, what would be the split of cases between vaccine and placebo?

If we flip a coin 94 times, how many heads will we get? That depends on the probability of the coin showing "heads". I didn't say the coin is fair so this probability does not have to be 50%. Now substitute "head" with "vaccine" and "tail" with "placebo", and you can see that the two questions are mathematically equivalent.

What is the probability of the coin showing "head" (i.e. "vaccine")? Is it 75%, the assumed VE for this scenario? Not really! The VE of 75% means the vaccine's case rate is a quarter of the placebo's. What we need is how many of the 94 cases turn out to come from the vaccine group. This number is the VSC - the proportion of cases coming from the vaccine group.

This is why the Bayesian model works with VSC instead of VE. As I showed in the previous post, for each value of VSC, we can derive the value of VE so the model can be interpreted in terms of VE, which is more intuitive.

Now back to the simulations. Each simulation starts with a different randomly-picked value of the VSC. This value of VSC is then used in our 94 coin tosses, after which we tally up how many "vaccine" sides turn up. So for that simulation, we have a split of cases between vaccine and placebo. We repeat these simulations, and average the results. This leads to the posterior probability curve (the black line) shown in the first chart of this post.

 

P.S. So far, I have ignored one small complication. In all these charts, I limited the VE to between 0% and 100% (from no effect to completely vanquishing the virus). However, VE can be negative! The vaccine could do more harm than good. In the VSC world, the share of cases can only go from 0% to 100%, and there are no other possible values. Note that when VSC exceeds 50%, i.e. the vaccine group accounts for more than half of the infections, VE is actually negative. The negative side of the VE scale is tricky - as in theory, VE can be negative infinity (all cases coming from the vaccine group).

[Added 11/18/2020] Here is the chart of the prior probability lifted from Sebastian Kranz's post. The right half of the probability density corresponds to theta > 0.5, meaning there are more cases in the vaccine group than the placebo group, i.e. VE is negative. Do they really believe that the vaccine has such a high chance (38%) of causing more cases than the placebo?

 

P.P.S. The analysis methodology is standard. As with most of statistics, the devil is in the definitions, the counting, and the exclusions.

 

 

 

 

How Pfizer plans to analyze the vaccine trial data

While we eagerly await the scientific report on the coronavirus vaccine trial from Pfizer/BioNTech, some of us in the blogosphere have checked out the Pfizer protocol in which they have pre-registered the method of analysis. This post is intended for those who are curious how the data from such trials are analyzed.

In a vaccine trial, the key metric is vaccine efficacy (VE). Pfizer previewed the "interim" analysis in a press release in which they claimed their vaccine has a VE of "above 90 percent". (Regardless of the optimism of this announcement, science by press release is an undesirable development. Such pre-announcements with few details prevents experts from offering any critique - the only reasonable response is to accept the headline tentatively, allowing it to take root.) 

VE of 90% means the vaccine cut infections by 90 percent. It does not mean 90 percent of those who are vaccinated will not get sick (link). That's because the coronavirus isn't spreading as quickly as you might presume. Inferring from Pfizer's press release, I think about 1 percent of the unvaccinated participants got infected within 30 days. So cutting infections by 90 percent means the case rate in the vaccinated group is 0.1 percent. Throughout this discussion, I refer to cases instead of infections because it's not clear from the protocol that Pfizer will detect every infection.

Since the participants are equally split between the vaccine and placebo groups, the average case rate is (1+0.1)%/2 = 0.55%. Out of every 10,000 participants, 55 have gotten sick. There were 10 times fewer sick people in the vaccine group than the placebo group. So 5 out of 55 are from the vaccine group, and 50 out of 55 are from the placebo group. The vaccine group accounts for 5/55 = 1/11 = 9 percent of the total cases. I'm going to call this the Vaccine Share of Cases (VSC) metric.

The Pfizer analysis does not work directly with the VE metric. It focuses on the VSC metric. (For those who've read the protocol, VSC is the entity known as theta.) Notice that the higher the VE, the more powerful is the vaccine, the smaller is the VSC. Those two metrics are inversely related. I did the calculation above to show you that any value of VE corresponds to some other value of VSC. Maximizing VE is the same as minimizing VSC. (The analytical reason for using VSC instead of VE will be explained below.)

The following chart shows the inverse relationship between VE and VSC.

The claimed VE of 90% implies that 9% of the cases were traced back to the vaccine group. The supposed minimum standard of 50% efficacy corresponds to cases split 1/3, 2/3 between the vaccine and placebo groups.

***

When the statisticians receive the data from the trial, they can estimate the VSC in the trial population. The Pfizer press release said at the time of the interim analysis, they recorded 94 cases in total. We just computed that 90% VE corresponds to 9% VSC, so 9% of 94 or roughly 8 cases came from the vaccinated group and 86 cases came from the placebo group. If the vaccine is completely useless, the cases would split evenly, 47 cases each.

That analysis is inadequate, though, because it does not take into account sampling error. Imagine re-running history with a different set of trial participants. It's quite possible that the VSC in this trial would be around 9 percent but it is highly unlikely that it would end up at exactly 9 percent.

Instead of concluding that VSC is 9%, as observed in the trial, Pfizer statisticians assign probabilities for the value of VSC in a Bayesian-style analysis. Like this:

In this chart (a probability density curve), the value of VSC runs along the horizontal axis, and the height represents the probability. The peak of this density curve gives the most likely value of VSC, which is 8.4%.

Statisticians focus on the middle 95% of the probability, i.e. the VSC range from 4% to 16%. The FDA scores the trial based on the chance that vaccine's share of cases is below 41%, which, if you glance back at the first chart of the post, corresponds to VE at least 30 percent. (This means that the FDA does not actually require at least 50 percent vaccine efficacy as frequently reported.)

Below is the same probability curve in terms of VE rather than VSC:

This chart is more intuitive. As you might suspect, the shape of the curve is influenced by the specific outcomes of the vaccine trial. The above chart assumes the cases split 8 and 86 between vaccine and placebo groups.

In the next chart, you can see how the probability curve changes if the trial had returned different results.

Some experts suggest that the vaccine cases could be as low as 3. For that trial outcome, the Pfizer analysis yields the blue dashed line shown below. The probability mass shifts to the right (higher VE) and also becomes tighter (lower uncertainty).

I include also the purple dashed line for the scenario in which the vaccine and placebo groups both produced 47 cases. As expected, the probability mass is now concentrated on the opposite end of VE. Even though the measured efficacy in the trial is 0%, the Bayesian model assigns probabilities of VE to be up to 50%.

If you're interested in even more technical details, Sebastian's post is your next stop.

***

As I explained the other day, the Bayesian analysis requires a "prior." The prior expresses the analyst's subjective expectation of the answer without the benefit of the data from the trial. Once the data become available, the prior is modified to yield the probability curves shown above.

The following chart compares the prior and "posterior" density curves.

You can see how far the prior has been pushed by the strong data observed in the Pfizer trial.

The shape of the prior (pink line) is rather curious though. The Pfizer statisticians seem to believe that the most likely VE is 100%, which is shockingly optimistic. They hedge by spreading the rest of the probability rather evenly across all possible values. The data from the vaccine trial pull the peak away from 100% (there were confirmed infections among the vaccinated), while concentrating all probabilities to between 70% and 100%.

Recall that the FDA wants to know the chance that VE is larger than 30%. The black line tells us that the strong data from the trial indicate an almost certainty. What was the researchers' prior expectation? That's asking how much of the probability mass of the pink line lies to the right of 30%. According to my software, it is 57%54%. [Oops, I forgot to switch from VE to VSC.]

In plain English: prior to running the vaccine trial, Pfizer believes their vaccine has a 65 in 10 chance of reaching at least 30% efficacy (the FDA minimum requirement), and based on the strong result in the trial, they now think it is almost certain that the vaccine has at least 30% efficacy.

 

P.S. I tied up a couple of loose ends in this next post.

 

[11/26/2020: In a new post, I put up some Python code that can be used to generate the charts shown on this post.]