Methodology for Examples 3 and 4:

Bayesian Analysis of the Incidence Rate Ratio (IRR)

Examples 3 and 4 use BetaBinomial() to reanalyze data that were obtained in two Phase 3 clinical trials and widely reported in the U.S. news media in late 2020 and early 2021. Both were two-group randomized designs that compared vaccine to placebo. The Pfizer-BioNTech trial used the beta-binomial approach delineated here. The Janssen (Johnson and Johnson) used a frequentist method not discussed here.

Let IR.vaccine and IR.placebo be incidence rates (here, per subject-year) in the vaccine and placebo groups. The incidence event (IE) might be as a confirmed a first-time COVID infection, a hospitalization, or a death. The incidence rate ratio is IRR = IR.vaccine/IR.placebo. For example, IRR = 0.70 indicates that vaccine subjects have a 30% lower incidence rate than placebo subjects. This 30% number is the vaccine's efficacy, VE = 1 - IRR.

The Pfizer-BioNTech reports state that they used the transform P = IRR/(1 + IRR) and modeling it as P ~ beta(a, b). (P was denoted by θ (theta) in the reports cited above.) But as shown below, the transform was actually P = OSTR*IRR/(1 + OSTR*IRR), where OSTR (overall surveillance time ratio) is the ratio of the groups' surveillance times, vaccine vs. placebo. The Pfizer-BioNTech reported 2214 (vaccine) and 2222 (placebo) surveillance years , so OSTR = 2214/2222 = 0.996. For the Janssen trial, OSTR = 3113.9/3089.1 = 1.008. These are so close to 1.0, because the randomization was 1:1 and the numbers of confirmed COVID cases were so small relative to the numbers randomized. Example 4 deals with the AstraZeneca trial, but those data have not yet been released. Because it used a randomization of 2:1, its OSTR will likely be close to 2.0.

The rationale behind using P = OSTR*IRR/(1 + OSTR*IRR) is straightforward, but it unnecessarily oversimplifies. Let NAR.vaccine[j] and NAR.placebo[j] be the numbers of subjects in each group who were at risk for the IE at the time of the j-th IE (either group, ordered by calendar time). Then the probability that the j-th IE would have come from the vaccine group is

P[j] = NAR.vaccine[j]*IR.vaccine/(NAR.vaccine[j]*IR.vaccine + NAR.placebo[j]*IR.placebo)
= NARR[j]*IRR/(NARR[j]*IRR + 1)

where NARR[j] = NAR.vaccine[j]/NAR.placebo[j] is the number at-risk ratio when the j-th IE occurred. The reverse transform is

IRR = P[j]/{(1 - P[j])*NARR[j]}.

However, if we don't have the individual subject data needed to compute the NARR[j] values, we must compromise and use the reported group surveillance times, and substitute their ratio, OSTR for every NARR[j] value. This simplifies P[j] to

P = OSTR*IRR/(OSTR*IRR + 1),

with the reverse transform being

IRR = P/[(1 - P)*OSTR)].

Such a model pretends that all events have the same NARR, an unnecessary simplification if you have all event and censoring times (as would be in the dataset for any such clinical trial).

This leads to the transform functions P.IRR(IRR, OSTR), which transforms IRR to P, and IRR.P(P, OSTR)), which transforms P to IRR.

P.IRR <- function(IRR, OSTR) { return(OSTR*IRR/(1 + OSTR*IRR)) }

IRR.P <- function(P, OSTR) { return(P/((1 - P)*OSTR)) }

These transforms are well behaved (monotonic and differentiable), so the probability density function for IRR ~ betaIRR(a, b | OSTR) is easily derived from that of P ~ beta(a, b). In base R, qbeta(q, a, b) computes 100q% quantiles, pbeta(p, a, b) computes cumulative probabilities, Prob[P < p], and dbeta(p, a, b) computes densities at P = p. Corresponding functions for IRR ~ betaIRR(a, b | OSTR) have been written and are given in the Rfuncs code:

qbetaIRR(q, a, b, OSTR) computes the 100q% quantile of IRR ~ betaIRR(a, b | OSTR).
pbetaIRR(irr, a, b, OSTR) computes Prob[IRR < irr | OSTR].
dbetaIRR(irr, a, b, OSTR) computes the density of betaIRR(a, b | OSTR) at IRR = irr.

A core component of this methodology is that any quantile of P ~ beta(a, b) transforms directly to the corresponding quantile of IRR ~ betaIRR(a, b | OSTR), and vice-versa. Formally, qbetaIRR(q, a, b,OSTR) = IRR.P(P=qbeta(q, a, b), OSTR), and qbeta(q, a, b) = P.IRR(IRR=qbetaIRR(q, a, b, OSTR), OSTR). Thus, we can easily and simply summarize a betaIRR(a, b | OSTR) distribution by reporting its median, qbetaIRR(0.50, a, b, OSTR) and, say, its 95% probability (credible) interval, [qbetaIRR(0.025, a, b, OSTR), qbetaIRR(0.975, a, b, OSTR). The 95% spread, PIspread(0.95) = qbetaIRR(0.975, a, b, OSTR) - qbetaIRR(0.025, a, b, OSTR) gives a single measure of variability. (For a Normal variate, PIspread(0.95)/4 is the standard deviation.).

Two ancillary functions within BetaBinomial() describe the characteristics of a given IRR ~ betaIRR(a, b | OSTR) distribution:

DescribeIRR() tables the central tendencies (means, medians, modes) and the 95% PI and spread of IRR ~ betaIRR(a, b | OSTR).
DP.IRR() tables densities and cumulative probabilities of IRR ~ betaIRR(a, b | OSTR).

In typical beta-binomial analyses, P ~ beta(1.0, 1.0) is considered to be a non-informative prior. All P have density 1.0. The mean and median for P are both 0.50. That's all fine, but P is not the variable of interest here.

Using DP.IRR() and DescribeIRR() reveals that IRR ~ betaIRR(1.0, 1.0 | OSTR=1) has a maximum density of 1.00 at IRR = 0.00, then it falls to an asymptotic limit of 0.00. Its extreme right skewness gives a mean over 13 versus a median of 1.0. Critically, unless OSTR = 1, IRR ~ betaIRR(1.0, 1.0 | OSTR) has unacceptable medians. Thus, IRR ~ betaIRR(1.0, 1.0 | OSTR) cannot be said to be a non-informative prior for analyzing IRR.

> DP.IRR(a=1, b=1, OSTR=1,"M.1",
+    irr=c(0,0.001,0.05,0.5,0.7,1,2,10^10),
+    names.irr=c("0.000", "0.001", "0.050", "0.500",
+    "0.700", "1.000" ,"2.000", "10^10"))

Table M.1. Densities & Cum. Prob's of IRR ~ betaIRR(1.00, 1.00, OSTR=1).
*************************************************************************
   irr: 0.000 0.001 0.050 0.500 0.700 1.000 2.000 10^10
   Density(irr)    1.000 0.998 0.907 0.444 0.346 0.250 0.111 0.000
Prob[IRR < irr]    0.000 0.001 0.048 0.333 0.412 0.500 0.667 1.000
*************************************************************************

> DescribeIRR(parms=matrix(c(1, 1, 1/2,
+    1, 1, 2/3,
+    1, 1, 1,
+    1, 1, 1.5,
+    1, 1, 2),
+    ncol=3, byrow=T), table="M.2")

Table M.2. Descriptive measures, IRR ~ betaIRR(a, b | OSTR).
******************************************************************************
   a b OSTR Mean Median Mode    PI95 Spread95
Case 1 1.00 1.00 0.500 26.785 2.000 0.000 [0.051, 78.000] 77.949
Case 2 1.00 1.00 0.667 20.089 1.500 0.000 [0.038, 58.500] 58.462
Case 3 1.00 1.00 1.000 13.393 1.000 0.000 [0.026, 39.000] 38.974
Case 4 1.00 1.00 1.500 8.928 0.667 0.000 [0.017, 26.000] 25.983
Case 5 1.00 1.00 2.000 6.696 0.500 0.000 [0.013, 19.500] 19.487
******************************************************************************

Setting priors for IRR. The one-to-one correspondence between the quantiles of P and IRR undergirds the following methodology.

Formulate the problem directly in terms of a prior median for IRR, IRR_0.50, and a q100% quantile, IRR_Q.
Transform those IRR quantiles to corresponding P quantiles, specifically, P_0.50 = P.IRR(IRR_0.50, OSTR) and P_Q = P.IRR(IRR.Q, OSTR).
Use BetaBinomial()'s Median, P.Q, and Q arguments to specify the prior median and q100% quantile for P ~ beta(a0, b0).
BetaBinomial() fits a P ~ beta(a0, b0) distribution that has quantiles of P_0.50 and P_Q.
As specified by the arguments Q.points and CP.points, BetaBinomial() computes a set of quantiles and cumulative probabilities for the fitted P ~ beta(a0, b0) distribution. However, using the Trans and TransNames arguments, BetaBinomial() will transform those quantiles and cumulative probabilities for P into the corresponding values for IRR.

For example, let us set a diffuse, null prior for IRR as one having a median of 1.00 and a 5% quantile of 0.10. Such prior quantiles would be set in the study protocol. Because a0 and b0 also depend on the overall surveillance time ratio (OSTR), a0 and b0 cannot be set until the dataset is frozen for analysis. Thus, the median for P ~ beta(a0, b0) is P.IRR(1.00, OSTR) and its 5% quanitle is P.IRR(0.10, OSTR). The inverse transform is IRR(P) = P/((1-P)*OSTR). Here we use two overall surveillance time ratios, 1.03 and 2.05.

OSTR = 1.03. Suppose the randomization ratio was 1:1, which led to OSTR = 1.03, meaning that the total surveillance times for the vaccine group was 3% greater than for the placebo group. The fitted prior is P ~ beta(a0=1.43, b0=1.40). The fitted IRR median and 5% quantiles are 0.999 (versus 1.000 requested) and 0.1000 (versus 0.1000 requested).

> OSTR=1.03
> BB_1.03 <- BetaBinomial(Median=P.IRR(1.00, OSTR),
+ P.Q=P.IRR(0.10, OSTR), Q=0.05,
+ Trans="P/((1-P)*OSTR)", TransName="IRR")

Transformation: t(P) = IRR(P) = P/((1-P)*OSTR)

Solution for prior: P ~ beta(a0 = 1.434, b0 = 1.402)

Fit of beta(a0, b0) to requested Median (IRR_0.50) and IRR_0.05.
*****************************************************************
   Median    IRR_0.05
As requested: 1.000    0.1000
As fit: 0.999    0.1000
*****************************************************************

Prior: P ~ beta(a0 = 1.43, b0 = 1.40)

Median, 95% PI and its spread, IRR(P) = P/((1-P)*OSTR).
*******************************************************
Median Prob. Interval PI Spread
IRR    0.999 [0.059, 17.437]    17.4
*******************************************************

100q% quantiles of t(P) = IRR(P) = P/((1-P)*OSTR).
*******************************************************************
   q: 0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
t(q) 0.018 0.059 0.100 0.413 0.999 2.427 10.185 17.437 57.424
*******************************************************************

What are the densities, etc. for IRR ~ betaIRR(1.434, 1.402) | OSTR=1.03)? This distribution is unimodal, but the values for its mode (0.176), median (0.999), and mean (3.4) indicate right skewness. Of course, beta(1, 1) | OSTR=1.03) is not unimodal, and although its median of 0.971 is tolerable, it is extremely right skewed.

> DP.IRR(a=BB_1.03$a0, b=BB_1.03$b0, OSTR=1.03,"M.3",
+    irr=c(0,0.001,0.05,0.5,0.7,1,2,10^10),
+    names.irr=c("0.000", "0.001", "0.050", "0.500",
+    "0.700", "1.000" ,"2.000", "10^10"))

Table M.3. Densities & Cum. Prob's of IRR ~ betaIRR(1.43, 1.40, OSTR=1.03).
****************************************************************************
   irr: 0.000 0.001 0.050 0.500 0.700 1.000 2.000 10^10
   Density(irr)    0.000 0.114 0.542 0.523 0.422 0.308 0.130 0.000
Prob[IRR < irr]    0.000 0.000 0.020 0.298 0.392 0.500 0.702 1.000
****************************************************************************

> DescribeIRR(parms=matrix(c(BB_1.03$a0, BB_1.03$b0, 1.03,
+    1, 1,    1.03),
+    ncol=3, byrow=T), table="M.4")

Table M.4. Descriptive measures, IRR ~ betaIRR(a, b | OSTR).
******************************************************************************
   a b OSTR Mean Median Mode    PI95 Spread95
Case 1 1.43 1.40 1.030 3.388 0.999 0.176 [0.059, 17.437] 17.378
Case 2 1.00 1.00 1.030 13.003 0.971 0.000 [0.025, 37.864] 37.839
******************************************************************************

OSTR = 2.05. Suppose the randomization ratio was 2:1, which led to OSTR = 2.05. Then the fitted prior is P ~ beta(a0 = 1.43, b0 = 1.40). The fitted IRR median and 5% quantiles are 0.999 (versus 1.00) and 0.100 (versus (0.100).

> OSTR=2.05
> BB_2.05 <- BetaBinomial(Median=P.IRR(1.00, OSTR),
+ P.Q=P.IRR(0.10, OSTR), Q=0.05,
+ Trans="P/((1-P)*OSTR)", TransName="IRR")

Transformation: t(P) = IRR(P) = P/((1-P)*OSTR)

Solution for prior: P ~ beta(a0 = 1.7002050, b0 = 1.0001000)

Fit of beta(a0, b0) to requested Median (IRR_0.50) and IRR_0.05.
*****************************************************************
   Median    IRR_0.05
As requested: 1.000 0.100
As fit: 0.969 0.101
*****************************************************************

Prior: P ~ beta(a0 = 1.70, b0 = 1.00)

Median, 95% PI and its spread, IRR(P) = P/((1-P)*OSTR).
*******************************************************
Median Prob. Interval PI Spread
IRR    0.969 [0.063, 32.502]    32.4
*******************************************************

100q% quantiles of t(P) = IRR(P) = P/((1-P)*OSTR).
********************************************************************
   q: 0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
t(q) 0.023 0.063 0.101 0.387 0.969 2.645 15.921 32.502 165.124
********************************************************************

The fitting algorithm is best gleaned by studying the code. Note that the solution here is not exact not exact (for the median, 0.969 vs. 1.00). In some cases, it will not succeed at all (but with an appropriate execution stop and message). The reasons:

Given the observed OSTR, there may not be a P ~ beta(a0, b0) distribution that has nearly the requested two quantiles, P_0.50 and P_Q. In such cases, BetaBinomial() stops and prints an appropriate message.
The fitting is iterative, terminating when the fitting "error" for P ~ beta(a0, b0) is tolerably close to 0.00, but rarely actually 0.00.
The fitting pertains to P, which is a probability ranging between 0 to 1, whereas IRR is a ratio ranging from 0 to infinity. The difference between the requested and fitted quantiles of IRR ~ betaIRR(a0, b0 ! OSTR) can seem magnified compared to those of P ~ beta(a0, b0).

IRR ~ betaIRR(1.7002, 1.0001) | OSTR=2.05) is unimodal, but the values for its mode (0.171), median (0.969), and mean (11.2) indicate substantial right skewness. betaIRR(1.0, 1.0 | OSTR=2.05), with a median of 0.488, is so far from being a null prior as to render it unacceptable.

> DP.IRR(a=BB_2.05$a0, b=BB_2.05$b0, OSTR=2.05,"M.5",
+    irr=c(0,0.001,0.05,0.5,0.7,1,2,10^10),
+    names.irr=c("0.000", "0.001", "0.050", "0.500",
+    "0.700", "1.000" ,"2.000", "10^10"))

Table M.5. Densities & Cum. Prob's of IRR ~ betaIRR(1.70, 1.00, OSTR=2.05).
****************************************************************************
   irr: 0.000 0.001 0.050 0.500 0.700 1.000 2.000 10^10
   Density(irr)    0.000 0.045 0.543 0.528 0.406 0.284 0.115 0.000
Prob[IRR < irr]    0.000 0.000 0.018 0.314 0.407 0.509 0.690 1.000
****************************************************************************

> DescribeIRR(parms=matrix(c(BB_2.05$a0, BB_2.05$b0, 2.05,
+ 1, 1,    2.05),
+    ncol=3, byrow=T), table="M.6")

Table M.6. Descriptive measures, IRR ~ betaIRR(a, b | OSTR).
******************************************************************************
   a b OSTR Mean Median Mode    PI95 Spread95
Case 1 1.70 1.00 2.050 11.229 0.969 0.171 [0.063, 32.502] 32.439
Case 2 1.00 1.00 2.050 6.533 0.488 0.000 [0.013, 19.024] 19.012
******************************************************************************

The above computations had no connection to any subject-matter problem, real or fictitious. Examples 2, 3, and 4 all stem directly from published data related to seeking the first FDA approvals to begin emergency vaccinations to stem the COVID-19 pandemic that began near the beginning of 2020.