Phase II Power Analysis

Using the sample information from the phase I trial, we can calculate power for a phase II trial.

Table 2

Proper Bayesian power analysis

We wish to identify sample sizes for high dose and standard dose (in a 2:1 ratio) that will be able to detect a positive difference in the probability $p$ of a $\ge 4$ fold titer rise to H3N2 influenza. That is, we wish to be able to detect:

$$\delta = p^{(HD)} - p^{(SD)} > 0$$

with 80% power. This analysis will be conducted using the proper Bayesian approach (Spiegelhalter et al. 2004). The advantage of this approach for power analysis over a frequentist power analysis is that it does not assume that the underlying probabilities are known. Failure to do so will tend to underestimate variance and result in underpowered experiments.

A sample leading to a successful trial (one in which a difference is detected) under a particular sample size $n$ is:

$$E_n = \left\{\mathbf{x}; Pr( \delta > 0 | \mathbf{x} ) \ge 0.95\right\}$$

so the probability of a successful trial is:

$$\psi(n) = P_{X^{(n)}}(X^{(n)} \in E_n)$$

which is the expected Bayesian power.

In our case, the power is:

$$\psi(n) = \Phi\left(\frac{\pi_{hd} - \pi_{sd}}{\sqrt{\pi_{hd}(1-\pi_{hd})/n_{hd} - \pi_{sd}(1-\pi_{sd})/n_{sd}}}\right)$$

where the $\pi_i$ are the posterior probabilities, parameterized by a beta distribution with a binomial likelihood for the data.

We will use the sample data from the phase I trial (GiaQuinta et al. 2014) to construct prior distributions for $p^{(HD)}$ and $p^{(SD)}$.

In [1]:
%matplotlib inline
import numpy as np
import seaborn as sns
import pymc3 as pm
import theano.tensor as tt
plt = sns.plt
In [2]:
# Standard normal from complimentary error function
Φ = lambda x: 0.5*tt.erfc(-x/tt.sqrt(2))

def power_analysis(n):
    
    # Allocate samples 2:1
    n_hd = int(2*n/3)
    n_sd = n - n_hd
    
    with pm.Model() as model:

        # probabilities of >= 4-fold-rise taken from Beta priors constructed from phase I data
        π_hd = pm.Beta('π_hd', 12, 10)
        π_sd = pm.Beta('π_sd', 2, 13)
        
        # Draw binomial samples from probabilities
        x_hd = pm.Binomial('x_hd', n_hd, π_hd)
        x_sd = pm.Binomial('x_sd', n_sd, π_sd)
        
        # Empirical probabilities
        π_hd_post = x_hd/n_hd
        π_sd_post = x_sd/n_sd
        
        # Convert difference in probabilities to standard normal
        d = pm.Deterministic('d', (π_hd_post - π_sd_post)
                             / tt.sqrt(π_hd_post*(1-π_hd_post)/n_hd + π_sd_post*(1-π_sd_post)/n_sd))
        
        # P-value
        p = 1 - Φ(d)
        
        # Successful detection
        power = pm.Deterministic('power', p<0.05)
        
    return model

Use model to estimate power for a range of sample sizes.

n=60:

In [3]:
iterations = 30000
keep = 20000

mod60 = power_analysis(60)
with mod60:
    tr = pm.sample(iterations, step=pm.Metropolis(), init=None)
    
pm.df_summary(tr[-keep:], varnames=['power'])[['mean', 'mc_error']]
100%|██████████| 30000/30000 [00:07<00:00, 3802.60it/s]
Out[3]:
mean mc_error
power 0.8896 0.006621

n=48:

In [4]:
iterations = 30000
keep = 20000

mod48 = power_analysis(48)
with mod48:
    tr = pm.sample(iterations, step=pm.Metropolis(), init=None)
    
pm.df_summary(tr[-keep:], varnames=['power'])[['mean', 'mc_error']]
100%|██████████| 30000/30000 [00:08<00:00, 3689.98it/s]
Out[4]:
mean mc_error
power 0.8685 0.0075

n=42:

In [5]:
mod42 = power_analysis(42)
with mod42:
    tr = pm.sample(iterations, step=pm.Metropolis(), init=None)
    
pm.df_summary(tr[-keep:], varnames=['power'])[['mean', 'mc_error']]
100%|██████████| 30000/30000 [00:07<00:00, 3967.43it/s]
Out[5]:
mean mc_error
power 0.836 0.008768

n=36:

In [6]:
mod36 = power_analysis(36)
with mod36:
    tr = pm.sample(iterations, step=pm.Metropolis(), init=None)
    
pm.df_summary(tr[-keep:], varnames=['power'])[['mean', 'mc_error']]
100%|██████████| 30000/30000 [00:07<00:00, 3774.80it/s]
Out[6]:
mean mc_error
power 0.8212 0.008525

n=27:

In [7]:
mod27 = power_analysis(27)
with mod27:
    tr = pm.sample(iterations, step=pm.Metropolis(), init=None)
    
pm.df_summary(tr[-keep:], varnames=['power'])[['mean', 'mc_error']]
100%|██████████| 30000/30000 [00:08<00:00, 3413.36it/s]
Out[7]:
mean mc_error
power 0.74735 0.009091

The power analysis suggests that a sample size in the low- to mid-30s would be sufficient to detect a difference between high and standard dose with 80% power in a similar population. A sample size over 60 would be required for 90% power.

References

  1. GiaQuinta, S., Michaels, M. G., McCullers, J. A., Wang, L., Fonnesbeck, C., O'Shea, A., et al. (2014). Randomized, double-blind comparison of standard-dose vs. high-dose trivalent inactivated influenza vaccine in pediatric solid organ transplant patients. Pediatric Transplantation, 19(2), 219–228. http://doi.org/10.1111/petr.12419
  2. Spiegelhalter, D. J., Abrams, K. R., and Myles, J. P. (2004). Bayesian approaches to Clinical Trials and Health Care Evaluation. Wiley, New York.