R/BayesPPD-package.R

In BayesPPD: Bayesian Power Prior Design

## usethis namespace: start
#' @importFrom Rcpp sourceCpp
#' @useDynLib BayesPPD,.registration = TRUE
## usethis namespace: end
NULL

#' Bayesian sample size determination using the power and normalized power prior for generalized linear models
#'
#' The \pkg{BayesPPD} (Bayesian Power Prior Design) package provides two categories 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 \eqn{a_0} can be fixed or modeled as random using a normalized power prior.
#' 
#' @details
#' 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 \eqn{\mu_t}
#' and the parameter for the control group by \eqn{\mu_c}. Suppose there are \eqn{K} historical datasets \eqn{D_0 = (D_{01},\cdots, D_{0K})'}. We consider the following normalized power prior
#' for \eqn{\mu_c} given multiple historical datasets \eqn{D_0}
#' \deqn{\pi(\mu_c|D_0,a_0) = \frac{1}{C(a_0)}\prod_{k=1}^K \left[L(\mu_c|D_{0k})^{a_{0k}}\right]\pi_0(\mu_c)}
#' where \eqn{a_0 = (a_{01},\cdots,a_{0K})'}, \eqn{0\le a_{0k} \le 1} for \eqn{k=1,\cdots,K}, \eqn{L(\mu_c|D_{0k})} is the historical data likelihood,
#' \eqn{\pi_0(\mu_c)} is an initial prior, and \eqn{C(a_0)=\int \prod_{k=1}^K [L(\mu_c|D_{0k})^{a_{0k}}]\pi_0(\mu_c)d\mu_c}. When \eqn{a_0} is fixed,
#' the normalized power prior is equivalent to the power prior
#' \deqn{\pi(\mu_c|D_0,a_0) = \prod_{k=1}^K \left[L(\mu_c|D_{0k})^{a_{0k}}\right]\pi_0(\mu_c).}
#' By default, the power/type I error calculation algorithm assumes the null and alternative hypotheses are given by
#' \deqn{H_0: \mu_t - \mu_c \ge \delta} and \deqn{H_1: \mu_t - \mu_c < \delta,} where \eqn{\delta} is a prespecified constant. To test hypotheses of
#' the opposite direction, i.e., \eqn{H_0: \mu_t - \mu_c \le \delta} and \eqn{H_1: \mu_t - \mu_c > \delta} , one can set the parameter \code{nullspace.ineq} to "<".
#' To determine Bayesian sample size, we estimate the quantity \deqn{\beta_{sj}^{(n)}=E_s[I\{P(\mu_t-\mu_c<\delta|y^{(n)}, \pi^{(f)})\ge \gamma\}]}
#' where \eqn{\gamma > 0} is a prespecified posterior probability threshold for rejecting the null hypothesis (e.g., \eqn{0.975}), the probability is computed with respect to the posterior distribution given the data
#' \eqn{y^{(n)}} and the fitting prior \eqn{\pi^{(f)}}, and the expectation is taken with respect to the marginal distribution of \eqn{y^{(n)}}
#'  defined based on the sampling prior \eqn{\pi^{(s)}(\theta)}, where \eqn{\theta=(\mu_t, \mu_c, \eta)} and \eqn{\eta} denotes any nuisance parameter in the model.
#'  Let \eqn{\Theta_0} and \eqn{\Theta_1} denote the parameter spaces corresponding to \eqn{H_0} and \eqn{H_1}.
#' Let \eqn{\pi_0^{(s)}(\theta)} denote a sampling prior that puts mass in the null region, i.e., \eqn{\theta \subset \Theta_0}.
#' Let \eqn{\pi_1^{(s)}(\theta)} denote a sampling prior that puts mass in the alternative region, i.e., \eqn{\theta \subset \Theta_1}.
#' Then \eqn{\beta_{s0}^{(n)}} corresponding to \eqn{\pi^{(s)}(\theta)=\pi_0^{(s)}(\theta)} is a Bayesian type I error,
#' while \eqn{\beta_{s1}^{(n)}} corresponding to \eqn{\pi^{(s)}(\theta)=\pi_1^{(s)}(\theta)} is a Bayesian power.
#' We compute \eqn{n_{\alpha_0} = \min\{n: \beta_{s0}^{(n)} \le \alpha_0\}} and \eqn{n_{\alpha_1} = \min\{n: \beta_{s1}^{(n)} \ge 1-\alpha_1\}}.
#' Then Bayesian sample size is max\eqn{\{n_{\alpha_0}, n_{\alpha_1}\}}. Choosing \eqn{\alpha_0=0.05} and \eqn{\alpha_1=0.2}
#' guarantees that the Bayesian type I error rate is at most \eqn{0.05} and the Bayesian power is at least \eqn{0.8}.
#' 
#' To compute \eqn{\beta_{sj}^{(n)}}, the following algorithm is used:
#' \describe{
#' \item{Step 1:}{Generate \eqn{\theta \sim \pi_j^{(s)}(\theta)}}
#' \item{Step 2:}{Generate \eqn{y^{(n)} \sim f(y^{(n)}|\theta)}}
#' \item{Step 3:}{Compute \eqn{P(\mu_t < \mu_c + \delta|y^{(n)}, \pi^{(f)})}}
#' \item{Step 4:}{Check whether \eqn{P(\mu_t < \mu_c + \delta|y^{(n)}, \pi^{(f)}) \ge \gamma}}
#' \item{Step 5:}{Repeat Steps 1-4 \eqn{N} times}
#' \item{Step 6:}{Compute the proportion of times that \eqn{\{\mu_t < \mu_c + \delta|y^{(n)}, \pi^{(f)} \ge \gamma\}} is true out of the \eqn{N} simulated datasets, which gives an estimate of \eqn{\beta_{sj}^{(n)}}.}
#' }
#'
#' For positive continuous data assumed to follow exponential distribution, the hypotheses are given by
#' \deqn{H_0: \mu_t/\mu_c \ge \delta} and \deqn{H_1: \mu_t/\mu_c < \delta,} where \eqn{\mu_t} and \eqn{\mu_c} are the hazards for the treatment and the control group, respectively.
#' The definition of \eqn{\beta_{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 \eqn{\beta_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
#' \deqn{H_0: \beta_1 \ge \delta} and \deqn{H_1: \beta_1 < \delta.}
#' The definition of \eqn{\beta_{sj}^{(n)}} and the algorithm change accordingly. 
#' 
#' By default, the package assumes the historical data is
#' composed of control group subjects only. If the user wants to use historical data to inform treatment effect, one can set \code{borrow.treat=TRUE} 
#' and include the treatment indicator in the historical covariate matrix. 
#'
#' 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 \eqn{\theta} (matrix if \eqn{\theta} 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 \eqn{a_0},
#' numerical integration using the \pkg{RcppNumerical} package is used.
#'
#' @references Chen, Ming-Hui, et al. "Bayesian design of noninferiority trials for medical devices using historical data." Biometrics 67.3 (2011): 1163-1170.
#' @docType package
#' @name BayesPPD-package
NULL
#> NULL