top of page

Example 1
Does inhaled anosine improve lung function in severe ZOVID?

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.

When the research question is motivated by a hypothesis, then we can go on to assess how much that focal estimate supports or undermines the hypothesis.

​

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.

 

In frequentism, we obtain estimates and confidence intervals for that parameter, and, when the research question is a bona fide hypothesis, we can compute and try to interpret a p-value.

 

In Bayesianism, we set a sound prior distribution for the focal parameter and then let the data revise that prior to give us a posterior distribution. From that posterior, we can compute compute a posterior mean or median, a set of posterior quantiles (giving us, say, 80%, 90%, 95% and/or 99% probability/credible intervals), and posterior probabilities of our choosing. Critically, plotting the prior and posterior density functions gives us rich images, which tap into our native visual processing abilities.

 

Note. This example is fiction, but mirrors a real trial I collaborated in designing in 2021 to treat COVID-19.

​

Investigators. Brandon Pike, MD, the study's lead investigator, an expert in pulmonary O2 and CO2 exchange. Cybill Varing, MS, an innovative and modern Bayesian statistical scientist.

​

Subjects. Patients so severely ill with ZOVID that they had acute respiratory distress syndrome (ARDS), a condition with a mortality rate that exceeds 45% under high-quality standard care. They were put on ventilators immediately and treated within four hours under study protocol.

​

Treatment question. Does aerosolized inhaled anosine (fictitious) improve lung function in these patients? Whether the patient survives is not the whole story. A survivor may still have permanently weakened lung function, a sequela in "long ZOVID."

​

Primary outcome measure. The PaO2/FiO2 ratio (real) is a "hard" measure of oxygen uptake; let's just call it the O2ratio. Higher the better. Dr. Pike was adamant that by measuring the O2ratio at baseline and 5 days later, a straightforward algorithm reliably and validly distinguishes patients as having markedly improved (+) or not markedly improved (-) over this super-critical time period.

​

Statistical endpoint. So, this boils down to a binary outcome. Ms. Varing is typically wary about dichotomizing a continuous measure, but comes to endorse it in this case.

​

Research goal. Demonstrate "sufficient" clinical efficacy: P = Prob[+ | treated] > 0.60. If this study used a frequentist structure, the null hypothesis would be H0: P < 0.60.

​

Planning scenarios. To shape her statistical planning, Ms. Varing engaged Dr. Pike in an "elicitation session." She asked: If this study had a control group in which patients were given high-quality standard care, what proportion would have "markedly improved" O2ratios over the first five days? Let's call this P0 = Prob[+ | not treated]. This generated considerable discussion and clarification, but a scenario emerged. Given the new Omega variant of ZOVID, there was little good data so far, but P0 = 0.30 "sounds about right," and P0 = 0.45 is "almost surely" close to the upper limit. Dr. Pike also said that he believed that the inhaled anosine treatment was likely 80% efficacious, Prob[+ |  treated] = 0.80, a conjecture necessary for the sample-size analysis discussed below.

​

Bayesian approach. The study conforms neatly with the ordinary beta-binomial model. For a given prior distribution, P0 ~ beta(a0, b0), and m1 patients who markedly improve and m2 who do not, the posterior distribution is P1 ~ beta(a0+m1, b0+m2).

 

Two priors, both legitimate but markedly different, will be used:

Non-informative

Prior

​

​

​​

​

Tailored

Prior

P0 ~ beta(1, 1) is the uniform(0, 1) distribution, which gives every value of P0 the same density, In effect, it says "We know nothing whatsoever about Prob[+ | treated]", so let the posterior distribution, P1 ~ beta(1+m1, 1+m2) be dominated by the data. The resulting posterior median and 95% probability (credible) interval (PI) most closely resemble the ordinary frequentist estimate and 95% confidence interval.

​

P0 ~ beta(a0, b0) is fit so that it has median 0.30 and Prob[P0 > 0.60] = 0.05. This follows the principle of the Scientific Method, which has us postulate a believable state of nature contrary to our research hypothesis. We then challenge the study to produce data that convincingly refutes that null state in a way that supports the research hypothesis. Setting the median of P0 at 0.30 and Prob[P0 > 0.60] = 0.05 is a more formidable challenge than using the non-informative prior, which sets the median at 0.50 and Prob[P0 > 0.60] = 0.40. This fitted tailored prior (run Ex1.Tprior) is P0 ~ beta(2.49, 5.38), which has median 0.3008 and 0.95 quantile 0.6000.

>       Ex1.Tprior <- BetaBinomial(Median=0.30, P.Q=0.60, Q=0.95)


Solution for prior:  P0 ~ beta(a0 = 2.49, b0 = 5.38)

Fit of beta(a0, b0) to requested Median (P_0.50) and P_0.95.
*************************************************************
               Median      P_0.95
As requested:  0.3000      0.6000
      As fit:  0.3008      0.6000
*************************************************************


Median, 95% PI, and its spread, P0 ~ beta(2.49, 5.38).
******************************************************
    Median   Prob. Interval   PI Spread
P    0.301  [0.0657, 0.655]       0.589
******************************************************


100q% quantiles of P0 ~ beta(2.49, 5.38).
****************************************************************
     q:  0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
P(q)     0.033 0.066 0.089 0.197 0.301 0.421 0.600 0.655 0.752
****************************************************************

Visually comparing priors. The Rfunc PlotBetaDensities() can be used to easily graph the density functions of these two priors. (Also employed here is PlotDirector(), an Rfunc utility documented and downloadable elsewhere.)

     PlotDirector(PlotSize=list(w=8, h=5),CloseOld=T)  # opens graphics window
     
     PlotBetaDensities(
       Title3="Comparing Non-Informative and Tailored Priors",
       XLabel="P: Probability of Improvement",
       Parms1=c(1,1),
       Parms2=c(Ex1.Tprior$a0, Ex1.Tprior$b0),
       DensLabel=c("Non-informative", "Tailored"),
       LineType=c(1,1),
       LineColor=c("purple","red"),
       ShiftRt=c(.75,-0.05),
       ShiftUp=c(0.005,0.02))
     lines(c(0.60,0.60),c(0,1.6),col="blue")
     text(0.60, 1.85, "P = 0.60", pos=1, col="blue",cex=0.85,srt=0)
    

Non-informative, beta(1, 1), vs. tailored, beta(2.49, 5.38).

Efficacy analysis 1a: Non-informative, (0+, 0-) ⇒ (4+, 0-)

Prior: Non-informative.

​

Data: Under a typical Phase I protocol, four patients were treated and monitored closely with frequent blood, urine, sputum, and nasal swab samples. No serious adverse events or toxicology signals occurred. All four patients had markedly improved O2ratios after five days (4+ vs. 0-).

>     Ex1a <- BetaBinomial(A0=1, B0=1,
+                           CP.points=c(0.45, 0.60, 0.70, 0.80),
+                           M=c(4,0))

 

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

Posterior:  P1 ~ beta(a1 = 1.00 + 4 = 5.00, b1 = 1.00 + 0 = 1.00)

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.871   [0.478, 0.995]       0.517
**********************************************


100q% quantiles of P0 ~ beta(1.00, 1.00) and P1 ~ beta(5.00, 1.00).
*********************************************************************
          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.347 0.478 0.549 0.758 0.871 0.944 0.990 0.995 0.999
*********************************************************************


Prob[P0 < p] and Prob[P1 < p].
*****************************************
          p:  0.450 0.600 0.700 0.800
    Prior     0.450 0.600 0.700 0.800
Posterior     0.018 0.078 0.168 0.328
****************************************

Comment. Under the beta(1, 1) prior, posterior median and 95% PI are 0.871 and [0.478, 0.995]. Prob[P > 0.60] = 0.40. Having observed 4+ vs. 0-, this progresses to Prob[P > 0.60] = 1 - 0.078 = 0.922. Based on only four subjects, this seems too favorable, but compare that to the frequentist estimate of 1.00 (which makes no sense!), 95% CI (score/Wilson's method) of [0.510, 1.00], and p-value of 0.130 for H0: P ≤ 0.60.

​

Prior and posterior density plots.

    path <- "/Users/ralphobrienCWRU_MBPro07/storeroom"
   PlotDirector(PlotSize=list(w=8, h=5), PlotType=".pdf", Directory=path)

   PlotBetaDensities(
     Title1="Analysis 1a",
     Title2="Non-informative prior: P0 ~ beta(1, 1)",
     Title3="First data: 4 improvements in 4 patients (4+, 0-)",
     XLabel="P: Probability of Improvement",
     Parms1=c(1,1),
     Parms2=c(Ex1a$a1, Ex1a$b1),
     LineType=c(1,1),
     LineColor=c("darkgreen","red"),
     ShiftRt=c(0,0),
     ShiftUp=c(0.01,0))
   lines(c(0.60,0.60),c(0,1.6),col="blue")
   text(0.60, 2, "P = 0.60", pos=1, col="blue",cex=0.85,srt=0)

   dev.off()
   setwd("/Users/ralphobrienCWRU_MBPro07")

Efficacy analysis 1a: Non-informative, (0+, 0-) ⇒ (4+, 0-).

Efficacy analysis 1b: Tailored, (0+, 0-) ⇒ (4+, 0-)

Prior: Tailored

Data: 4+, 0-.

> setwd("/Users/ralphobrienCWRU_MBPro07")
>     Ex1b <- BetaBinomial(Median=0.30, P.Q=0.60, Q=0.95,
+                           CP.points=c(0.45, 0.60, 0.70, 0.80),
+                          Q.points=c(.025,0.05,0.25,0.50,0.75,0.95,0.975,0.995),
+                           M=c(4,0))


Solution for prior:  P0 ~ beta(a0 = 2.49, b0 = 5.38)

Fit of beta(a0, b0) to requested Median (P_0.50) and P_0.95.
*************************************************************
               Median      P_0.95
As requested:  0.3000      0.6000
      As fit:  0.3008      0.6000
*************************************************************

Posterior:  P1 ~ beta(a1 = 2.49 + 4 = 6.49, b1 = 5.38 + 0 = 5.38)

Median and 95% PI, and its spread, P ~ beta().
**********************************************
          Median   Prob. Interval   PI Spread
    Prior  0.301   [0.066, 0.655]       0.589
Posterior  0.550   [0.273, 0.806]       0.532
**********************************************


100q% quantiles of P0 ~ beta(2.49, 5.38) and P1 ~ beta(6.49, 5.38).
*******************************************************************
          q:  0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
    Prior     0.066 0.089 0.197 0.301 0.421 0.600 0.655 0.752
Posterior     0.273 0.314 0.450 0.550 0.647 0.771 0.806 0.864
*******************************************************************


Prob[P0 < p] and Prob[P1 < p].
*****************************************
          p:  0.450 0.600 0.700 0.800
    Prior     0.797 0.950 0.987 0.998
Posterior     0.250 0.634 0.857 0.972
*****************************************

Comment. Under the P0 ~ beta(2.49, 5.38) prior, the posterior median and 95% PI are 0.550 and [0.273, 0.806]. While these are naturally markedly lower than those obtained using P0 ~ beta(1, 1), they are still encouraging given only four patients were treated. Moreover, while the tailored prior initially sets Prob[P0 > 0.60] = 0.05 to model healthy scientific skepticism, observing 4+ vs. 0- moves this to Prob[P1 > 0.60] = 1 - 0.634 = 0.366, a substantial trek forward in this endeavor. Finally, there is no frequentist counterpart to this analysis.

​

Prior and posterior density plots.

    path <- "/Users/ralphobrienCWRU_MBPro07/storeroom"
   PlotDirector(PlotSize=list(w=8, h=5), PlotType=".pdf", Directory=path)

   PlotBetaDensities(
     Title1="Analysis 1b",
     Title2=
       "Tailored prior: median(P0) = 0.30, Prob[P0 > 0.60] = 0.05",
     Title3="First data: 4 improvements in 4 patients (4+, 0-)",
     XLabel="P: Probability of Improvement",
     Parms1=c(Ex1b$a0, Ex1b$b0),
     Parms2=c(Ex1b$a1, Ex1b$b1),
     LineType=c(1,1),
     LineColor=c("darkgreen","red"))
   lines(c(0.60,0.60),c(0,2.7),col="blue")
   text(0.60, 2.95, "P = 0.60", pos=1, col="blue",cex=0.85,srt=0)

   dev.off()
   setwd("/Users/ralphobrienCWRU_MBPro07")

Efficacy analysis 1b: Tailored, (0+, 0-) ⇒ (4+, 0-).

Efficacy analyses 1c1 and 1c2: Non-informative, (0+, 0-) ⇒ (4+, 0-) and (4+, 0-) ⇒ (25+, 4-)

Initial prior: non-informative, P0 ~ beta(1, 1)
Additional data: 21 improvements in 25 patients.
Accumulated data: 4 + 21 = 25 improvements in 4 + 25 = 29 patients.

​

Analyses 1c1 and 1c2 produce the same posterior distribution, but arrive there using different priors and data. Analysis 1c1 uses the initial non-informative prior and the accumulated data. Analysis 1c2 uses only the additional data, so its prior is the posterior distribution produced by Analysis 1a, which used the initial non-informative prior and the (4+, 0-) data. This shows how the endeavor progressed from the 4th patient studied to the 29th.

​

For comparison, the function CI.P.score() computes a two-sided confidence interval for P using the score method.

​

Efficacy analysis 1c1: Non-informative, (0+, 0-) ⇒ (25+, 4-).

>   Ex1c1 <- BetaBinomial(A0=1, B0=1,
+                           CP.points=c(0.45, 0.60, 0.70, 0.80),
+                           M=c(25,4))


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

Posterior:  P1 ~ beta(a1 = 1.00 + 25 = 26.00, b1 = 1.00 + 4 = 5.00)

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.846   [0.693, 0.944]       0.251
**********************************************


100q% quantiles of P0 ~ beta(1.00, 1.00) and P1 ~ beta(26.00, 5.00).
*********************************************************************
          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.637 0.693 0.720 0.799 0.846 0.886 0.932 0.944 0.962
*********************************************************************


Prob[P0 < p] and Prob[P1 < p].
***************************************
          p:  0.450 0.600 0.700 0.800
    Prior     0.450 0.600 0.700 0.800
Posterior     0.000 0.002 0.030 0.255
***************************************

​

>     CI.P.score(Ns=25, Nf=4, CI.level=0.95)
Ns Nf  N   P.est CI.level                 CI
25  4 29 0.86207     0.95 [0.69441, 0.94503]

Comment. The P0 ~ beta(1, 1) prior lets the data dominate, so the posterior statistics are quite similar to their frequentist counterparts:

   Posterior Median, 95% PI        Frequentist estimate, 95% CI
  ************************        ****************************
     0.846, [0.693, 0.944]             0.862, [0.694, 0.945]

 

The advantage of the Bayesian approach is most evident by considering the simple and direct statements provided by the posterior probabilities Prob[P > 0.60] = 0.998 and Prob[P > 0.70] = 0.970.

​

Prior and posterior density plots.

    PlotDirector(PlotSize=list(w=8, h=5),CloseOld=T)  # opens graphics window

   PlotBetaDensities(
     Title1="Analysis 1c1",
     Title2="Non-informative initial prior: P0 ~ beta(1, 1)",
     Title3="Accumulated data: 25 improvements in 29 patients (25+, 4-)",
     XLabel="P: Probability of Improvement",
     Parms1=c(1,1),
     Parms2=c(Ex1c1$a1, Ex1c1$b1),
     DensLabel=c("Initial Prior","Posterior"),
     LineType=c(1,1),
     LineColor=c("darkgreen","red"),
     ShiftUp=c(0.01,0),
     ShiftRt=c(0,-0.1))
   lines(c(0.60,0.60),c(0,2),col="blue")
   text(0.60, 2.6, "P = 0.60", pos=1, col="blue",cex=0.85,srt=0)
   lines(c(0.60,0.60),c(2.5,6.5),col="blue")

Efficacy analysis 1c1: Non-informative, (0+, 0-) ⇒ (25+, 4-).

Efficacy analysis 1c2: Non-informative, (4+, 0-) ⇒ (25+, 4-).

>   Ex1c2 <- BetaBinomial(A0=Ex1a$a1, B0=Ex1a$b1,
+                           CP.points=c(0.45, 0.60, 0.70, 0.80),
+                           M=c(21,4))


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

Posterior:  P1 ~ beta(a1 = 5.00 + 21 = 26.00, b1 = 5.00 + 4 = 5.00)

Median and 95% PI, and its spread, P ~ beta().
**********************************************
          Median   Prob. Interval   PI Spread
    Prior  0.871   [0.478, 0.995]       0.517
Posterior  0.846   [0.693, 0.944]       0.251
**********************************************


100q% quantiles of P0 ~ beta(5.00, 5.00) and P1 ~ beta(26.00, 5.00).
*********************************************************************
          q:  0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
    Prior     0.347 0.478 0.549 0.758 0.871 0.944 0.990 0.995 0.999
Posterior     0.637 0.693 0.720 0.799 0.846 0.886 0.932 0.944 0.962
*********************************************************************


Prob[P0 < p] and Prob[P1 < p].
***************************************
          p:  0.450 0.600 0.700 0.800
    Prior     0.018 0.078 0.168 0.328
Posterior     0.000 0.002 0.030 0.255
***************************************

Comment. Even though Analyses 1c1 and 1c2 arrive at the same posterior distribution for P, Analysis 1c2 directly depicts how the 25 new observations changes our understanding of P.

​

Prior and posterior density plots.

  PlotDirector(PlotSize=list(w=8, h=5),CloseOld=T)  # opens graphics window

  PlotBetaDensities(
   Title1="Analysis 1c2",
   Title2="Current prior: posterior from non-informative prior & (4+, 0-)",
   Title3="Additional data: 21 improvements in 25 patients (21+, 4-)",
   XLabel="P: Probability of Improvement",
   Parms1=c(Ex1a$a1, Ex1a$b1),
   Parms2=c(Ex1c2$a1, Ex1c2$b1),
   DensLabel=c("Previous Posterior", "New Posterior"),
   LineType=c(1,1),
   LineColor=c("darkgreen","red"),
   ShiftUp=c(-0.06,0),
   ShiftRt=c(-.50,-0.19))
  lines(c(0.60,0.60),c(0,2),col="blue")
  text(0.60, 2.6, "P = 0.60", pos=1, col="blue",cex=0.85,srt=0)

Efficacy analysis 1c2: Non-informative, (4+, 0-) ⇒ (25+, 4-).

Efficacy analyses 1d1 and 1d2: Tailored, (0+, 0-) ⇒ (4+, 0-) and (4+, 0-) ⇒ (25+, 4-)

Initial prior: non-informative, P0 ~ beta(2.49, 5.38), which has median 0.30 and 0.95 quantile 0.60.
Additional data: 21 improvements in 25 patients.
Accumulated data: 4 + 21 = 25 improvements in 4 + 25 = 29 patients.

​

Analyses 1d1 and 1d2 produce the same posterior distribution, but arrive there using different priors and data. Analysis 1d1 uses the initial tailored prior and the accumulated data. Analysis 1c2 uses only the additional data, so its prior is the posterior distribution produced by Analysis 1b, which used the initial tailored prior and the (4+, 0-) data. This shows how the endeavor progressed from the 4th patient studied to the 29th.

​

There is no clear frequentist analog to this analysis.

​

Efficacy analysis 1d1: Tailored, (0+, 0-) ⇒ (25+, 4-).

>   Ex1d1 <- BetaBinomial(Median=0.30, P.Q=0.60, Q=0.95,
+                           CP.points=c(0.45, 0.60, 0.70, 0.80),
+                           M=c(25,4))


Solution for prior:  P0 ~ beta(a0 = 2.49, b0 = 5.38)

Fit of beta(a0, b0) to requested Median (P_0.50) and P_0.95.
*************************************************************
               Median      P_0.95
As requested:  0.3000      0.6000
      As fit:  0.3008      0.6000
*************************************************************

Posterior:  P1 ~ beta(a1 = 2.49 + 25 = 27.49, b1 = 5.38 + 4 = 9.38)

Median and 95% PI, and its spread, P ~ beta().
**********************************************
          Median   Prob. Interval   PI Spread
    Prior  0.301   [0.066, 0.655]       0.589
Posterior  0.750   [0.596, 0.871]       0.275
**********************************************


100q% quantiles of P0 ~ beta(2.49, 5.38) and P1 ~ beta(27.49, 9.38).
*********************************************************************
          q:  0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
    Prior     0.033 0.066 0.089 0.197 0.301 0.421 0.600 0.655 0.752
Posterior     0.544 0.596 0.622 0.700 0.750 0.796 0.854 0.871 0.899
*********************************************************************


Prob[P0 < p] and Prob[P1 < p].
***************************************
          p:  0.450 0.600 0.700 0.800
    Prior     0.797 0.950 0.987 0.998
Posterior     0.000 0.028 0.251 0.768
***************************************

Comment. The posterior median and 95% PI are now at 0.750 and [0.596, 0.871]. Before data, Prob[P0 > 0.60] = 0.05. At this point, Prob[P1 > 0.60] = 1 - 0.028 = 0.972. Given the skepticism built into this intitial prior, this is strong, if not compelling evidence, in favor of inferring that the treatment more than doubles the 0.30 improvement rate that Dr. Pike conjectured for those not treated with anosine. In addition, Prob[P1 > 0.70] = 1 - 0.251 = 0.749.

​

Prior and posterior density plots.

  PlotDirector(PlotSize=list(w=8, h=5),CloseOld=T)  # opens graphics window

  PlotBetaDensities(
   Title1="Analysis 1d1",
   Title2=
     'Initial Prior: "Tailored," median P of 0.30, Prob[P > 0.60] = 0.05',
   Title3="Accumulated data: 25 improvements in 29 patients (25+, 4-)",
   XLabel="P: Probability of Improvement",
   Parms1=c(Ex1d1$a0, Ex1d1$b0),
   Parms2=c(Ex1d1$a1, Ex1d1$b1),
   DensLabel=c("Initial Prior", "Posterior"),
   LineType=c(1,1),
   LineColor=c("darkgreen","red"),
   ShiftRt=c(-0.16,0))
  lines(c(0.60,0.60),c(0,2),col="blue")
  text(0.60, 2.6, "P = 0.60", pos=1, col="blue",cex=0.85,srt=0)

Efficacy analysis 1d1: Tailored, (0+, 0-) ⇒ (25+, 4-).

Efficacy analysis 1d2: Tailored, (4+, 0-) ⇒ (25+, 4-).

>   Ex1d2 <- BetaBinomial(A0=Ex1b$a1, B0=Ex1b$b1,
+                           CP.points=c(0.45, 0.60, 0.70, 0.80),
+                           M=c(21,4))


Prior:  P0 ~ beta(a0 = 6.49, b0 = 5.38)
Posterior:  P1 ~ beta(a1 = 6.49 + 21 = 27.49, b1 = 5.38 + 4 = 9.38)

Median and 95% PI, and its spread, P ~ beta().
**********************************************
          Median   Prob. Interval   PI Spread
    Prior  0.550   [0.273, 0.806]       0.532
Posterior  0.750   [0.596, 0.871]       0.275
**********************************************


100q% quantiles of P0 ~ beta(6.49, 5.38) and P1 ~ beta(27.49, 9.38).
*********************************************************************
          q:  0.005 0.025 0.050 0.250 0.500 0.750 0.950 0.975 0.995
    Prior     0.203 0.273 0.314 0.450 0.550 0.647 0.771 0.806 0.864
Posterior     0.544 0.596 0.622 0.700 0.750 0.796 0.854 0.871 0.899
*********************************************************************


Prob[P0 < p] and Prob[P1 < p].
***************************************
          p:  0.450 0.600 0.700 0.800
    Prior     0.250 0.634 0.857 0.972
Posterior     0.000 0.028 0.251 0.768
***************************************

Comment. Even though Analyses 1d1 and 1d2 arrive at the same posterior distribution for P, Analysis 1d2 directly depicts how the 25 new observations changes our understanding of P.

​

Prior and posterior density plots.

  PlotDirector(PlotSize=list(w=8, h=5),CloseOld=T)  # opens graphics window

  PlotBetaDensities(
   Title1="Analysis 1d2",
   Title2="Current prior: posterior from tailored prior & (4+, 0-)",
   Title3="Additional data: 21 improvements in 25 patients (21+, 4-)",
   XLabel="P: Probability of Improvement",
   Parms1=c(Ex1d2$a0, Ex1d2$b0),
   Parms2=c(Ex1d2$a1, Ex1d2$b1),
   DensLabel=c("Previous Posterior", "New Posterior"),
   LineType=c(1,1),
   LineColor=c("darkgreen","red"),
   ShiftRt=c(-0.24,0))
  lines(c(0.60,0.60),c(0,2.9),col="blue")
  text(0.60, 3.4, "P = 0.60", pos=1, col="blue",cex=0.85,srt=0)

Visual comparison of posteriors from non-informative and tailored priors

  PlotDirector(PlotSize=list(w=8, h=5),CloseOld=T)  # opens graphics window

  PlotBetaDensities(
   Title2=
     "Comparing posteriors from non-informative and tailored initial priors",
   Title3="Accumulated data: 25 improvements in 29 patients (25+, 4-)",
   XLabel="P: Probability of Improvement",
   Parms1=c(Ex1c1$a1, Ex1c1$b1),
   Parms2=c(Ex1d1$a1, Ex1d1$b1),
   DensLabel=c("Non-Informative", "Tailored"),
   LineType=c(1,1),
   LineColor=c("purple","blue"),
   ShiftRt=c(-0.055,-0.12),
   ShiftUp=c(0.017,0))
  lines(c(0.60,0.60),c(0,2.9),col="red")
  text(0.60, 3.45, "P = 0.60", pos=1, col="red",cex=0.85,srt=0)

Comparing the posteriors from (25+, 4-) data under non-informative and tailored priors of beta(1, 1) and beta(2.49, 5.38).

Comment. The plot illustrates concretely how much the scientific conservatism built-in to  the tailored prior affects the posterior distribution relative to letting the data dominate. The two priors are not in conflict, they merely view the (25+, 4-) from different perspectives, and—fortunately—the two posteriors support with the same conclusion.

A minimally sufficient statistical planning analysis

This will surely be improved after discussing it with Cyndi and Paul, and then turned be into webpage material.

# Sample-size analysis
# ********************
# For a given sample size, N, and true value, Ptrue = Prob[+ | treated], the
# the number of marked improvements is M1 ~ binomial(N, Ptrue). Consider five
# possible outcomes:
#       "very unlucky" ... M1 is 2.5% quantile of binomial(N, Ptrue)
#   "somewhat unlucky" ... M1 is 25% quantile of binomial(N, Ptrue)
#        "anticipated" ... M1 is median of binomial(N, Ptrue)
#     "somewhat lucky" ... M1 is 75% quantile of binomial(N, Ptrue)
#         "very lucky" ... M1 is 97.5% quantile of binomial(N, Ptrue)
#
# Using the tailored prior, P0 ~ beta(a0, b0), that has median 0.30 and
# Prob[P0 > 0.60] = 0.05, we can compute the following posterior results under
# each possible outcome.
#     Median ... posterior median of P
#     PI95 ... 95% posterior probability (credible) interval
#     Prob.P.gt.60 ... posterior Prob[P > 0.60]

# Scenario for statistical planning: 80% of patients will "markedly improve."
# Planned sample size: 30.
N <- 30
Ptrue <- 0.80
Q <- c(0.025, 0.25, 0.50, 0.75, 0.975)
Outcome <- c('"very unlucky"', '"somewhat unlucky"', '"anticipated"',
            '"somewhat lucky"','"very lucky"')
prior <- BetaBinomial(Median=0.30, P.Q=0.60, Q=0.95,
                     CP.points=0.60,Print=F)
M1 <- integer(length(Q))

Median <- PI95 <- Prob.P.gt.60 <- character(length(Q))
for (i in 1:length(Q)) {
  M1[i] <- qbinom(Q[i], N, Ptrue)
  results <- BetaBinomial(A0=prior$a0, B0=prior$b0,
                         CP.points=0.60,
                         M=c(M1[i],N-M1[i]),Print=T)
  Median[i] <- formatC(results$median[2],digits=3,format="f")
  PI95[i] <- paste("[",formatC(results$PI[2,1],digits=3,format="f"),", ",
                   formatC(results$PI[2,2],digits=3,format="f"),"]",sep="")
  Prob.P.gt.60[i] <- formatC(1-results$cprob[2,1],digits=3,format="f")
}
{
print(data.frame(N, Ptrue),row.names=F)
cat("\n")
print(data.frame(Outcome, Q, M1, Pcnt=round(M1/N,3),
                Median, PI95,Prob.P.gt.60),row.names=F)
}

#  N Ptrue
# 30   0.8
#
#            Outcome     Q M1  Pcnt Median           PI95 Prob.P.gt.60
#     "very unlucky" 0.025 19 0.633  0.569 [0.409, 0.719]        0.349
# "somewhat unlucky" 0.250 23 0.767  0.676 [0.518, 0.811]        0.833
#      "anticipated" 0.500 24 0.800  0.703 [0.547, 0.833]        0.905
#   "somewhat lucky" 0.750 26 0.867  0.757 [0.605, 0.874]        0.979
#       "very lucky" 0.975 28 0.933  0.811 [0.667, 0.913]        0.997

# Comment. Under the scenario of Ptrue = 0.80, with N = 30 and using the
# tailored prior (which builds in considerable scientific skepticism), even if
# the outcome is "very unlucky" (Q = 0.025) with only 19 markedly improved
# patients (63%), the posterior median is 0.569 with a 95% PI  of [0.409, 0.719]
# and a posterior Prob[P > 0.60] = 0.349. So the chance is 95% that these
# results would be better than that.  For "somewhat unlucky" (Q = 0.25), there
# is a 75% chance that the posterior results will be better than median 0.676,
# 95% PI of [0.518, 0.811], and Prob[P > 0.60] = 0.833. There is a 50% chance
# that the results will be better than the "anticipated" outcome, which has
# median 0.703, 95% PI of [0.547, 0.833], and Prob[P > 0.60] = 0.905.

# Dr. Varing advises the research team that they should plan on studying 30
# patients, but stresses that using a Bayesian approach does not lock them to a
# preset fixed sample size (which frequentist methods are technically required
# to do in order to properly assess Type I and Type II error rates and
# confidence interval coverage). This trial could stop sooner if the successive
# posterior distributions convincingly support concluding efficacy or futility.
# The trial can also be extended beyond N = 30 if the posterior results show
# substantial promise but are not yet convincing, as per the "very unlucky"
# outcome.

A safety analysis

I'm still creating this, again.

bottom of page