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 ≥ δ*

and

*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 ≥ δ*

and

*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 ≥ δ*

and

*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.

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.