BayesPPD-package: Bayesian sample size determination using the power and...

Description Details References


The BayesPPD (Bayesian Power Prior Design) package provides two catogries of functions: functions for Bayesian power/type I error calculation and functions for model fitting. Supported distributions include normal, binary (Bernoulli/binomial), Poisson and exponential. The power parameter a_0 can be fixed or modeled as random using a normalized power prior.


Following Chen et al.(2011), for two group models (i.e., treatment and control group with no covariates), denote the parameter for the treatment group by μ_t and the parameter for the control group by μ_c. Suppose there are K historical datasets D_0 = (D_{01},\cdots, D_{0K})'. We consider the following normalized power prior for μ_c given multiple historical datasets D_0

π(μ_c|D_0,a_0) = \frac{1}{C(a_0)}∏_{k=1}^K ≤ft[L(μ_c|D_{0k})^{a_{0k}}\right]π_0(μ_c)

where a_0 = (a_{01},\cdots,a_{0K})', 0≤ a_{0k} ≤ 1 for k=1,\cdots,K, L(μ_c|D_{0k}) is the historical data likelihood, π_0(μ_c) is an initial prior, and C(a_0)=\int ∏_{k=1}^K [L(μ_c|D_{0k})^{a_{0k}}]π_0(μ_c)dμ_c. When a_0 is fixed, the normalized power prior is equivalent to the power prior

π(μ_c|D_0,a_0) = ∏_{k=1}^K ≤ft[L(μ_c|D_{0k})^{a_{0k}}\right]π_0(μ_c).

The power/type I error calculation algorithm assumes the null and alternative hypotheses are given by

H_0: μ_t - μ_c ≥ δ


H_1: μ_t - μ_c < δ,

where δ is a prespecified constant. To test hypotheses of the opposite direction, i.e., H_0: μ_t - μ_c ≤ δ and H_1: μ_t - μ_c > δ , one can recode the responses for the treatment and control groups. To determine Bayesian sample size, we estimate the quantity

β_{sj}^{(n)}=E_s[I\{P(μ_t-μ_c<δ|y^{(n)}, π^{(f)})≥ γ\}]

where γ > 0 is a prespecified posterior probability threshold for rejecting the null hypothesis (e.g., 0.975), the probability is computed with respect to the posterior distribution given the data y^{(n)} and the fitting prior π^{(f)}, and the expectation is taken with respect to the marginal distribution of y^{(n)} defined based on the sampling prior π^{(s)}(θ), where θ=(μ_t, μ_c, η) and η denotes any nuisance parameter in the model. Let Θ_0 and Θ_1 denote the parameter spaces corresponding to H_0 and H_1. Let π_0^{(s)}(θ) denote a sampling prior that puts mass in the null region, i.e., θ \subset Θ_0. Let π_1^{(s)}(θ) denote a sampling prior that puts mass in the alternative region, i.e., θ \subset Θ_1. Then β_{s0}^{(n)} corresponding to π^{(s)}(θ)=π_0^{(s)}(θ) is a Bayesian type I error, while β_{s1}^{(n)} corresponding to π^{(s)}(θ)=π_1^{(s)}(θ) is a Bayesian power. We compute n_{α_0} = \min\{n: β_{s0}^{(n)} ≤ α_0\} and n_{α_1} = \min\{n: β_{s1}^{(n)} ≥ 1-α_1\}. Then Bayesian sample size is max\{n_{α_0}, n_{α_1}\}. Choosing α_0=0.05 and α_1=0.2 guarantees that the Bayesian type I error rate is at most 0.05 and the Bayesian power is at least 0.8. To compute β_{sj}^{(n)}, the following algorithm is used:

Step 1:

Generate θ \sim π_j^{(s)}(θ)

Step 2:

Generate y^{(n)} \sim f(y^{(n)}|θ)

Step 3:

Compute P(μ_t < μ_c + δ|y^{(n)}, π^{(f)})

Step 4:

Check whether P(μ_t < μ_c + δ|y^{(n)}, π^{(f)}) ≥ γ

Step 5:

Repeat Steps 1-4 N times

Step 6:

Compute the proportion of times that \{μ_t < μ_c + δ|y^{(n)}, π^{(f)} ≥ γ\} is true out of the N simulated datasets, which gives an estimate of β_{sj}^{(n)}.

For positive continuous data assumed to follow exponential distribution, the hypotheses are given by

H_0: μ_t/μ_c ≥ δ


H_1: μ_t/μ_c < δ,

where μ_t and μ_c are the hazards for the treatment and the control group, respectively. The definition of β_{sj}^{(n)} and the algorithm change accordingly.

If there are covariates to adjust for, we assume the first column of the covariate matrix is the treatment indicator, and the corresponding parameter is β_1, which, for example, corresponds to a difference in means for the linear regression model and a log hazard ratio for the exponential regression model. The hypotheses are given by

H_0: β_1 ≥ δ


H_1: β_1 < δ.

The definition of β_{sj}^{(n)} and the algorithm change accordingly.

This implementation of the method does not assume any particular distribution for the sampling priors. The user is allowed to specify a vector or matrix of samples for θ (matrix if θ is of dimension >1) from any distribution, and the algorithm samples with replacement from the vector or matrix at each iteration of data simulation. In order to accurately approximate a joint distribution for multiple parameters, the number of iterations should be large (e.g., 10,000).

Gibbs sampling is used for normally distributed data. Slice sampling is used for all other data distributions. For two group models with fixed a_0, numerical integration using the RcppNumerical package is used.


Chen, Ming-Hui, et al. "Bayesian design of noninferiority trials for medical devices using historical data." Biometrics 67.3 (2011): 1163-1170.

BayesPPD documentation built on Sept. 8, 2021, 5:06 p.m.