Example 3
Beta-Binomial Bayesian Estimation of a Small Proportion

As stated in the opening to Example 1, most research questions are best addressed by estimating a "focal parameter" that relates as directly as possible to each question and by quantifying the statistical uncertainty about that estimate.

Suppose that Streptococcus heftans (fictitious) is an uncommon oral bacteria that in had been recently hypothesized to cause heart infections affecting chambers and valves, especially in the elderly and those with other heart disease. The prevalence for S. heftans, here denoted as P, had not yet been rigorously estimated, so Maya Cardiya, PhD, sought to estimate it and understand the statistical uncertainty of that estimate. She originally planned to test for the presence of S. heftans in 500 dental patients recruited from 10 clinics that served different economic/racial/ethnic populations.

Statistical planning. The study's progressive and astute statistician was Anne Nova, MS. Her Bayesian design called for using the standard beta-binomial method in which P ~ beta(a, b). If N patients are tested for the presence of S. heftans, the number of positives, m1, is taken to be m1 ~ binomial(N, Ptrue). Denoting the prior as P0 ~ beta(a0, b0), the posterior is beta(a1 = a0 + m1, b1 = b0 + N - m1).

When planning the study, Dr. Cardiya believed the true P was around 0.005 (0.5% of the population), but could be as high as 0.015, and was "virtually certain" to be less than 0.03. Given those values, Ms. Nova tailored the prior distribution to have median and 99% quantile of approximately 0.015 and 0.03. Using BetaBinomial(), she easily found that that P0 ~ beta(8.90, 562.75) fits this to four decimals.

> prior <- BetaBinomial(P0.50=0.015, Q=0.99, P0.Q=0.03,
+ CP.points=c(0.005, 0.01, 0.02, 0.03))

Solution for prior: P0 ~ beta(a0 = 8.90, b0 = 562.75)

Fit of beta(a0, b0) to requested P0.50 and P0.99.
*************************************************************
   P0.50    P0.99
As requested: 0.01500    0.03000
As fit: 0.01500    0.03000
*************************************************************

Median, 95% PI, and its spread, P0 ~ beta(8.90, 562.75).
********************************************************
Median    Prob. Interval PI Spread
P 0.0150 [0.00712, 0.0272]    0.0201
********************************************************

100q% quantiles of P0 ~ beta(8.90, 562.75).
*************************************************************************
   q: 0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
P(q)    0.0054 0.0071 0.0081 0.0118 0.0150 0.0187 0.0249 0.0272 0.0320
*************************************************************************

Prob[P0 < p], where P0 ~ beta(8.90, 562.75).
*********************************************
p: 0.0050 0.0100 0.0200 0.0300
Prob[P < p]    0.003 0.131 0.814 0.990
*********************************************

Ms. Nova then used BB.WhatIf() to assess whether using P0 ~ beta(8.90, 552.75) and data gathered from N = 500 is likely to yield a satisfactory posterior distribution. (Example 2 summarizes this methodology.)

> prior <- BetaBinomial(Median=0.01, P.Q=0.05, Q=0.99,
+ CP.points=c(0.005, 0.01, 0.025, 0.05))

Solution for prior: P ~ beta(a0 = 1.408, b0 = 106.751)

Fit of beta(a0, b0) to requested Median (P_0.50) and P_0.99.
*************************************************************
   Median    P_0.99
As requested: 0.0100    0.0500
As fit: 0.0102    0.0500
*************************************************************

Prior: P ~ beta(a0 = 1.41, b0 = 106.75)

Median, 95% PI and its spread, P ~ beta(1.41, 106.75).
******************************************************
Median Prob. Interval PI Spread
P 0.010 [0.001, 0.041] 0.040
******************************************************

100q% quantiles of P ~ beta(1.41, 106.75).
****************************************************************
   q: 0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
P(q) 0.000 0.001 0.001 0.005 0.010 0.018 0.034 0.041 0.057
****************************************************************

Prob[P < p], where P ~ beta(1.41, 106.75).
*******************************************
p: 0.005 0.010 0.025 0.050
Prob[P < p] 0.246 0.493 0.873 0.990
*******************************************

P ~ beta(1.41, 106.751 fits the requested median and 99% quantile superbly. The 95% probability interval of [0.001, 0.041] is wide for this problem, thus one can argue that this prior is essentially "non-informative."

Interim analysis. Dr. Cardiya's protocol calls for one interim analysis near 50% of testing. Of the first 241 subjects, 7 test positive for S. heftans.

Using the common beta-binomial model, the interim posterior distribution is P ~ beta(prior$a + 7, prior$b + 234) = beta(1.408 + 7, 106.751 + 234). BetaBinomial() handles this easily.

> interim <- BetaBinomial(Median=0.01, P.Q=0.05, Q=0.99,
+ CP.points=c(0.005, 0.01, 0.025, 0.05),
+ M=c(7, 234))

Solution for prior: P ~ beta(a0 = 1.408, b0 = 106.751)

Fit of beta(a0, b0) to requested Median (P_0.50) and P_0.99.
*************************************************************
   Median    P_0.99
As requested: 0.0100    0.0500
As fit: 0.0102    0.0500
*************************************************************

Prior: P ~ beta(a0 = 1.41, b0 = 106.75)
Posterior: P ~ beta(a1 = 1.41 + 7 = 8.41, b1 = 106.75 + 234 = 340.75)

Median and 95% PI and its spread, P ~ beta().
*********************************************

Median Prob. Interval PI Spread

Prior 0.010 [0.001, 0.041] 0.040
Posterior 0.023 [0.011, 0.043] 0.032
*********************************************

100q% quantiles of P ~ beta(1.41, 106.75) and beta(8.41, 340.75).
*********************************************************************
q: 0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
Prior 0.000 0.001 0.001 0.005 0.010 0.018 0.034 0.041 0.057
Posterior 0.008 0.011 0.012 0.018 0.023 0.029 0.039 0.043 0.050
*********************************************************************

Prob[P < p], where P ~ beta(1.41, 106.75) or beta(8.41, 340.75).
*****************************************************************
p: 0.005 0.010 0.025 0.050
Prior 0.246 0.493 0.873 0.990
Posterior 0.000 0.017 0.588 0.995
*****************************************************************

These data, 7 positives vs. 234 negatives, have increased the median from 0.010 (prior) to 0.023 (posterior) and the 95% spread has tightened from 0.040 to 0.032. The prior Prob[P< 0.01] of 0.493 has plunged to 0.017. These results, although only interim, strongly suggest that S. heftans is more prevalent than was initially believed.

Final analysis. Suppose testing 247 more patients finds 5 positives. Thus, the final posterior distribution is P ~ beta(a0 + 7 + 5, b0 + 234 + 242). This prior.to.final run measures the change from the prior to the final analysis.

> prior.to.final <- BetaBinomial(Median=0.01, P.Q=0.05, Q=0.99,
+    CP.points=c(0.005, 0.01, 0.025, 0.05),
+    M=c(7+5, 234+242))

Solution for prior: P ~ beta(a0 = 1.408, b0 = 106.751)

Fit of beta(a0, b0) to requested Median (P_0.50) and P_0.99.
*************************************************************
   Median    P_0.99
As requested: 0.0100    0.0500
As fit: 0.0102    0.0500
*************************************************************

Prior: P ~ beta(a0 = 1.41, b0 = 106.75)
Posterior: P ~ beta(a1 = 1.41 + 12 = 13.41, b1 = 106.75 + 476 = 582.75)

Median and 95% PI and its spread, P ~ beta().
*********************************************
Median Prob. Interval PI Spread
Prior 0.010 [0.001, 0.041] 0.040
Posterior 0.022 [0.012, 0.036] 0.024
*********************************************

100q% quantiles of P ~ beta(1.41, 106.75) and beta(13.41, 582.75).
*********************************************************************
q: 0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
Prior 0.000 0.001 0.001 0.005 0.010 0.018 0.034 0.041 0.057
Posterior 0.010 0.012 0.013 0.018 0.022 0.026 0.033 0.036 0.041
*********************************************************************

Prob[P < p], where P ~ beta(1.41, 106.75) or beta(13.41, 582.75).
******************************************************************
p: 0.005 0.010 0.025 0.050
Prior 0.246 0.493 0.873 0.990
Posterior 0.000 0.006 0.687 1.000
******************************************************************

This interim.to.final run measures the change from the interim to the final analysis.

> interim.to.final <- BetaBinomial(A0 = round(interim$a1,3),
+    B0 = round(interim$b1,3),
+    CP.points=c(0.005, 0.01, 0.025, 0.05),
+    M=c(5, 242))

Prior: P ~ beta(a0 = 8.408, b0 = 340.751)
Posterior: P ~ beta(a1 = 8.408 + 5 = 13.4, b1 = 340.751 + 242 = 582.8)

Median and 95% PI and its spread, P ~ beta().
*********************************************
Median Prob. Interval PI Spread
Prior 0.023 [0.011, 0.043] 0.032
Posterior 0.022 [0.012, 0.036] 0.024
*********************************************

100q% quantiles of P ~ beta(8.408, 340.751) and beta(13.4, 582.8).
*********************************************************************
q: 0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
Prior 0.008 0.011 0.012 0.018 0.023 0.029 0.039 0.043 0.050
Posterior 0.010 0.012 0.013 0.018 0.022 0.026 0.033 0.036 0.041
*********************************************************************

Prob[P < p], where P ~ beta(8.408, 340.751) or beta(13.4, 582.8).
******************************************************************
p: 0.005 0.010 0.025 0.050
Prior 0.000 0.017 0.588 0.995
Posterior 0.000 0.006 0.687 1.000
******************************************************************

The final median has barely changed from the interim median, going from 0.023 to 0.022. But the 95% PI spread tightens by 25% (1 - 0.024/0.032). With Prob[P < 0.010 | all data] = 0.006, Dr. Cardiya is confident that the prevalence of S. heftans is considerably greater than has been believed for years, an important finding in terms of public health. She reports the posterior median of 0.022 and 95% PI of [0.012, 0.036].

Informative vs. non-informative priors. In Bayesian lingo, the above prior is informative, because it is based on "believing" that P, the true prevalence rate, is very likely (99%) to be 0.05. Let us compare the above results to those derived from three common approaches that make no "informative guesses" about P is likely to be:

Wilson's (frequentist) method.
Beta-binomial (Bayesian) method with uniform prior, P ~ beta(1, 1).
Beta-binomial method with Jeffrey's prior, P ~ beta(1/2, 1/2).

Wilson's method was touted in the influential article by Agresti and Coull (1998, Approximate is better than ”exact” for interval estimation of binomial proportions. The American Statistician, 52(2):119–126). It is computed by the built-in R function prop.test(). Here, the sample proportion of positive tests is 12/488 = 0.025. The 95% confidence interval is [0.014, 0.042].

> prop.test(x=12, n=488, correct=F)

1-sample proportions test without continuity correction

data: 12 out of 488, null probability 0.5
X-squared = 441.18, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
0.01412152 0.04248505
sample estimates:
p
0.02459016

P ~ beta(1, 1) is identical to the standard uniform distribution in which all values for P from 0 to 1 are equally likely. This is the classical non-informative prior for this problem and represents the Bayesian counterpart to the above frequentist analysis. But in this case, P has virtually no chance of being greater that 0.10, so Bayesians will argue there is no justification for giving this a 90% prior probability.

> beta_uniform <- BetaBinomial(A0=1, B0=1,
+ CP.points=c(0.005, 0.01, 0.025, 0.05, 0.10),
+ M=c(12, 488-12))

Prior: P ~ beta(a0 = 1, b0 = 1)
This prior distribution for P not unimodal.

Posterior: P ~ beta(a1 = 1 + 12 = 13, b1 = 1 + 476 = 477)

Median and 95% PI and its spread, P ~ beta().
*********************************************
Median Prob. Interval PI Spread
Prior 0.500 [0.025, 0.975] 0.950
Posterior 0.026 [0.014, 0.042] 0.028
*********************************************

100q% quantiles of P ~ beta(1, 1) and beta(13, 477).
*********************************************************************
q: 0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
Prior 0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
Posterior 0.011 0.014 0.016 0.021 0.026 0.031 0.039 0.042 0.049
*********************************************************************

Prob[P < p], where P ~ beta(1, 1) or beta(13, 477).
****************************************************
p: 0.005 0.010 0.025 0.050 0.100
Prior 0.005 0.010 0.025 0.050 0.100
Posterior 0.000 0.002 0.450 0.996 1.000
****************************************************

The Jeffrey's prior, P ~ beta(1/2, 1/2) is often promoted for analyzing a single proportion. It is symmetric but U-shaped (having infinite densities at both 0 and 1), so it seems inappropriate for most such applications, especially this one.

> beta_Jeffreys <- BetaBinomial(A0=1/2, B0=1/2,
+ CP.points=c(0.005, 0.01, 0.025, 0.05, 0.10),
+ M=c(12, 488-12))

Prior: P ~ beta(a0 = 0.5, b0 = 0.5)
This prior distribution for P not unimodal.

Posterior: P ~ beta(a1 = 0.5 + 12 = 12.5, b1 = 0.5 + 476 = 476.5)

Median and 95% PI and its spread, P ~ beta().
*********************************************
Median Prob. Interval PI Spread
Prior 0.500 [0.002, 0.998] 0.997
Posterior 0.025 [0.014, 0.041] 0.028
*********************************************

100q% quantiles of P ~ beta(0.5, 0.5) and beta(12.5, 476.5).
*********************************************************************
q: 0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
Prior 0.000 0.002 0.006 0.146 0.500 0.854 0.994 0.998 1.000
Posterior 0.011 0.014 0.015 0.020 0.025 0.030 0.038 0.041 0.047
*********************************************************************

Prob[P < p], where P ~ beta(0.5, 0.5) or beta(12.5, 476.5).
************************************************************
p: 0.005 0.010 0.025 0.050 0.100
Prior 0.045 0.064 0.101 0.144 0.205
Posterior 0.000 0.003 0.505 0.998 1.000
************************************************************

How do these methods compare?

Table 1.1. Four analyses of 12 positives vs. 476 negatives for S. heftans.
**************************************************************************
   Prior    Posterior
**************************    ***************************
Beta Prior    Median 95% PI (CI) Median 95% PI (CI)
(1.4, 106.8) 0.010    [0.001, 0.041] 0.022 [0.012, 0.036]
(1.00, 1.00) 0.500    [0.025, 0.975] 0.026 [0.014, 0.042]
(0.50, 0.50) 0.500    [0.002, 0.998] 0.025 [0.014, 0.041]
Wilson's method*    n/a n/a    0.025 [0.014, 0.042]
**************************************************************************
*Instead of the median and 95% probability interval, the estimate
associated with Wilson's frequentist method is the ordinary sample
proportion, 12/488, and the 95% confidence interval is computed by
inverting the score statistic without the continuity correction.

The beta(1, 1) and beta(0.50, 0.50) distributions are regarded as being non-informative priors for analyzing a single proportion. As a result, compared to P ~ beta(1.4, 106.8), their estimates lean a little towards 0.50 and they have wider statistical intervals. In a sense, the informative prior of P ~ beta(1.4, 106.8) is akin to having already run a no-cost, no-time, no-effort study of N = 1.4 + 106.8 = 108.2 pseudo subjects and observing 1.4 positives. The payoff? The spread of the 95% PI for P ~ beta(1.4, 106.8), [0.012, 0.036], is about 16% narrower than that of the other three methods, which are all around [0.014, 0.041].

The lack of a more striking difference between the informative prior versus the two non-informative priors and the frequentist method is due to having a sample size, N = 12 + 476, that largely dominates the prior distribution. As the sample size increases, Bayesian methods progressively "disregard" the prior. (The next example, which is real, shows how the prior can be critical when the sample size is relatively small.)

The same table constructed for the interim analysis (N = 7 + 234 = 241; output from BetaBinomial() runs not shown here), reveals a greater advantage for using the informative prior. Here, using P ~ beta((1.4, 106.8) yields a posterior 95% PI that is about 28% narrower than those of the three other methods.

> beta_uniform <- BetaBinomial(A0=1, B0=1,
+ Q.points=c(0.025, 0.50, 0.975),
+ M=c(7, 234))

> beta_Jeffreys <- BetaBinomial(A0=1/2, B0=1/2,
+ Q.points=c(0.025, 0.50, 0.975),
+ M=c(7, 234))

> prop.test(x=7, n=241, correct=F)

(Output from these runs not shown.)

Table 1.2. Four analyses of 7 positives vs. 234 negatives for S. heftans.
*************************************************************************
   Prior    Posterior
**************************    **************************
Beta Prior    Median 95% PI (CI) Median 95% PI (CI)
(1.4, 106.8) 0.010    [0.001, 0.041] 0.023 [0.011, 0.043]
(1.00, 1.00) 0.500    [0.025, 0.975] 0.032 [0.014, 0.059]
(0.50, 0.50) 0.500    [0.002, 0.998] 0.030 [0.013, 0.056]
Wilson's method*    n/a    n/a 0.029 [0.014, 0.059]
*************************************************************************
*Instead of the median and 95% probability interval, the estimate
associated with Wilson's frequentist method is the ordinary sample
proportion, 7/241, and the 95% confidence interval is computed by
inverting the score statistic without the continuity correction.

Example 3 Beta-Binomial Bayesian Estimation of a Small Proportion

Example 3
Beta-Binomial Bayesian Estimation of a Small Proportion