A Unifying Bayesian Approach for Sample Size Determination Using Design and Analysis Priors

The classical power and sample size analysis problem in the frequentist setting is widely applied in diverse settings. For example, we can formulate a hypothesis concerning the population mean. A one-sided hypothesis test for the null $\displaystyle H_{0}:\theta=\theta_{0}$ against the alternative $H_{a}:\theta>\theta_{0}$, where $\theta$ is the population mean, is decided on the location of $\theta_{0}$ with respect to the distribition of the sample mean $\bar{y}$. Assuming a known value for the population variance $\sigma^{2}$ and that the sample mean’s distribition is (approximately) Gaussian, we reject $H_{0}$ if $\bar{y}>\theta_{0}+\frac{\sigma}{\sqrt{n}}Z_{1-\alpha}$, where $\alpha$ is the specified Type-I error and $Z_{1-\alpha}$ is the corresponding quantile of the Gaussian distribution. The statistical “power” of the test is $1-\beta$, where $\beta=P(\mbox{Rej }H_{0}\,|\,H_{a})$ is the Type-II error. Straightforward algebra yields the ubiquitous sample size formula, $\displaystyle n=(Z_{\alpha}+Z_{\beta})^{2}\left(\frac{\sigma}{\Delta}\right)^{2}$, derived from the power:

where $\Delta=\theta_{1}-\theta_{0}$ is the critical difference. Given any fixed value of $\Delta/\sigma$, the power curve is a function of sample size and can be plotted as in Figure 1. The sample size required to achieve a specified power can then be read from the power curve.

More generally, power curves, such as the one shown in Figure 1, may not be available in closed form but can be simulated for different sample sizes. They quantify our degree of assurance regarding our ability to meet our analysis objective (rejecting the null) for different sample sizes. Formulating sample size determination as a decision problem so that power is an increasing function of sample size, offers the investigator a visual aid in helping deduce the minimum sample size needed to achieve a desired power.

From a Bayesian standpoint, there is no need to condition on a fixed alternative. Instead, we determine the tenability of a hypothesis given the data we observe. A joint probability model is constructed for the parameters and the data using a prior distribition for the parameters and the likelihood function of the data conditional on the parameters. Inference proceeds from the posterior distribition of the parameters given the data.

In the design stage we have not observed the data. Therefore, we will need to formulate a data generating mechanism and, subsequently, consider the posterior distribution given the realized data to evaluate the tenability of a hypothesis regarding our parameter. We then use the probability law associated with the data generating mechanism to assign a degree of assurance of our analysis objective. As a sufficiently simple example, consider a situation where we seek to evaluate the tenability of $H:\theta>\theta_{0}$ given data from a Gaussian population with mean $\theta$ and a known variance $\sigma^{2}$. Assuming the prior $\theta\sim N\left(\theta_{1},\frac{\sigma^{2}}{n_{0}}\right)$, where $n_{0}$ reflects the precision of the prior relative to the data, and the likelihood of the sample mean $\displaystyle\bar{y}\sim N\left(\theta,\frac{\sigma^{2}}{n}\right)$, the posterior distribution of $\theta$ can then be obtained by multiplying the prior and the likelihood as

Our analysis objective is to ascertain whether $P(\theta>\theta_{0}\,|\,\bar{y})>1-\alpha$, where $\alpha$ is a fixed number. The posterior distribution in (2) gives us $P(\theta>\theta_{0}\,|\,\bar{y})$ and we define

as a Bayesian counterpart of statistical power, which we refer to as Bayesian assurance. The inner probability defines our analysis objective, while the outer probability defines our chances of meeting the analysis objective under the given data generating mechanism.

Our current manuscript intends to explore Bayesian assurance and subsequent sample size calculations as described through (2) and (3) in the general context of conjugate Bayesian linear regression. Of particular emphasis will be the data generating mechanism and providing motivation behind quantifying separate prior beliefs at the design and analysis stage of clinical trials (O’Hagan and Stevens,, 2001). The balance of this manuscript proceeds as follows. Section 2 presents the Bayesian sample size determination problem embedded within a conjugate Bayesian linear model. This offers us analytic tractability through which we motivate the need for different prior distributions for designing and analyzing a study. Section 2.3 considers the cases of known and unknown variances and offers corresponding algorithms. Section 3 casts the cost-effectiveness problem explored by O’Hagan and Stevens, (2001) in our framework, which we revisit using known and unknown variances, and we also compare with other Bayesian alternatives for sample size calculations including for inference using proportions. We close the paper in Section 4 with additional points of discussion and future aims.

Consider a proposed study where a certain number, say $n$, of observations, which we denote by $y_{n}=(y_{1},y_{2}\ldots,y_{n})^{\top}$, are to be collected in the presence of $p$ controlled explanatory variables, say $x_{1},x_{2},\ldots,x_{p}$, that will be known to the investigator for any unit $i$ at the design stage. Consider the usual normal linear regression setting such that $y_{n}=X_{n}\beta+\epsilon_{n}$, where $X_{n}$ is $n\times p$ with $i$-th row corresponding to $x_{i}^{\top}$ and $\epsilon_{n}\sim N(0,\sigma^{2}V_{n})$, where $V_{n}$ is a known $n\times n$ correlation matrix. Both $X_{n}$ and $V_{n}$ are assumed known for each sample size $n$ through design and modeling considerations. A conjugate Bayesian linear regression model specifies the joint distribution of the parameters $\{\beta,\sigma^{2}\}$ and the data as

Inference proceeds from the posterior distribution derived from (4),

where $a^{*}_{\sigma}=a_{\sigma}+n/2$, $b^{*}_{\sigma}=b_{\sigma}+(1/2)\left\{\mu_{\beta}^{\top}V_{\beta}^{-1}\mu_{\beta}+y_{n}^{\top}V_{n}^{-1}y_{n}-m_{n}^{\top}M_{n}m_{n}\right\}$, $M_{n}^{-1}=V_{\beta}^{-1}+X_{n}^{\top}V_{n}^{-1}X_{n}$ and $m_{n}=V_{\beta}^{-1}\mu_{\beta}+X_{n}^{\top}V_{n}^{-1}y_{n}$. Sampling from the joint posterior distribution of $\{\beta,\sigma^{2}\}$ is achieved by first sampling $\sigma^{2}\sim IG(a^{*}_{\sigma},b^{*}_{\sigma})$ and then sampling $\beta\sim N(M_{n}m_{n},\sigma^{2}M_{n})$ for each sampled $\sigma^{2}$. This yields marginal posterior samples from $p(\beta\,|\,y_{n})$, which is a non-central multivariate $t$ distribution, though we do not need to work with its complicated density function. See Gelman et al., (2013) for further details on the conjugate Bayesian linear regression model and sampling from its posterior.

We wish to decide whether our realized data will favor $H:u^{\top}\beta>0$, where $u$ is a $p\times 1$ vector of fixed contrasts. In practice, a decision on the tenability of $H$ is often based on the $100(1-\alpha)\%$ posterior credible interval,

obtained from the conditional posterior predictive distribution $p(\beta\,|\,\sigma,y_{n})$. The data would favor $H$ if $y_{n}$ belongs to the following set

This is equivalent to $0$ being below the two-sided $100(1-\alpha)\%$ credible interval for $u^{\top}\beta$. Practical Bayesian designs will seek to assure the investigator that the above criterion will be achieved with a sufficiently high probability through the Bayesian assurance,

which generalizes the definition in (3). Given the assumptions on the model, the fixed values of the parameters $\{\mu_{\beta},V_{\beta},\sigma,V_{n}\}$ and the fixed vector $u$ that determines the hypothesis being tested, the Bayesian assurance function evaluates the probability of rejecting the null hypothesis under the marginal probability distribution of the realized data corresponding to any given sample size $n$. Choice of sample size will be determined by the smallest value of $n$ that will ensure $\delta(n;\sigma,u,\mu_{\beta},V_{\beta},V_{n})>\gamma$, where $\gamma$ is a specified number.

Let us consider the special case when $X_{n}=1_{n}$ so that $\beta$ is a scalar with prior distribution $\beta\sim N(\beta_{1},\sigma^{2}/n_{0})$, where $n_{0}>0$ is a fixed precision parameter (sometimes referred to as “prior sample size”), $V_{n}=I_{n}$ and $H:\beta>\beta_{0}$. We decide in favor of $H$ if the data lies in

where the expression on the right reveals a convenient condition in terms of the sample mean. As $n_{0}\to 0$, i.e., the prior becomes vague, $S_{\alpha}(n;y,\sigma,\beta_{0},\beta_{1},n_{0})$ collapses to the critical region in classical inference for testing $H_{0}:\beta=\beta_{0}$ against $H_{a}:\beta=\beta_{1}$. The Bayesian assurance function is

where $\Delta=\beta_{1}-\beta_{0}$. Given $n_{0}$, we will compute the sample size needed to detect a critical difference of $\Delta$ with probability $1-\beta$ as $n=\arg\min\{n:\delta(\Delta,n)\geq 1-\beta\}$. However, the limiting properties of the function in (8) are not without problems. When the prior is vague, i.e., $n_{0}\to 0$, then

while in the case when the prior is precise, i.e., $n_{0}\to\infty$ we obtain

This is undesirable. Vague priors are customary in Bayesian analysis, but they propagate enough uncertainty that the marginal distribution of the data under the given model will force the assurance to be lower than 0.5. In other words, regardless of how large a sample size we have, we cannot assure the investigator with probability greater than 50% that the null hypothesis will be rejected. At the other extreme, where the prior is fully precise, it fully dominates the data (or the likelihood) and there is no information from the data that is used in the decision. Therefore, the assurance is a function of the prior only and we will always or never reject the null hypothesis depending upon whether $\Delta>0$ or $\Delta<0$.

In order to resolve this issue, we work with two different sets of priors, one at the design stage and another at the analysis stage. Building upon O’Hagan and Stevens (2001), we elucidate with the Bayesian linear regression model in the next section and offer a simulation-based framework for computing the Bayesian assurance curves.

We consider two scenarios that are driven by the amount of information given in a study. We develop the corresponding algorithms based on these assumptions. The first case assumes that the population variance $\sigma^{2}$ is known and the second case assumes $\sigma^{2}$ is unknown, prompting us to consider additional prior distributrions for $\sigma^{2}$ in the design and analysis stage. Our context remains testing the tenability of $H:u^{\top}\beta>C$ given realized data from a study, where $C$ is a known constant.

The first application selects a sample size based on the cost-effectiveness of new treatments undergoing Phase 3 clinical trials (O’Hagan and Stevens,, 2001). As delineated in Section 2, we construct the two-stage paradigm under the context of a conjugate linear model and generalize it to the case where the population variance $\sigma^{2}$ is unknown. We cast the example in O’Hagan and Stevens, (2001) within our framework to assess overall performance and our ability to emulate the analysis in O’Hagan and Stevens, (2001).

Consider designing a randomized clinical trial where $n_{1}$ patients are administered Treatment 1 and $n_{2}$ patients are administered Treatment 2 under some suitable model and study objectives. Let $c_{ij}$ and $e_{ij}$ denote the observed cost and efficacy values, respectively, corresponding to patient $j=1,2,\ldots,n_{i}$ receiving treatment $i$ for $i=1,2$ treatments, where $n_{i}$ is the number of patients in the $i$-th treatment group. Furthermore, the expected population mean cost efficacy under treatment $i$ are set to be $E(c_{ij})=\gamma_{j}$ and $E(c_{ij})=\mu_{i}$, respectively. The variances are taken to be $Var(c_{ij})=\tau_{i}^{2}$ and $Var(e_{ij})=\sigma_{i}^{2}$. For simplicity, we assume equal sample sizes within the treatment groups so that $n=n_{1}=n_{2}$. We also assume equal sample variances for the costs and efficacies such that $\tau^{2}$ = $\tau_{1}^{2}=\tau_{2}^{2}$ and $\sigma^{2}=\sigma_{1}^{2}=\sigma_{2}^{2}$.

O’Hagan and Stevens, (2001) utilize the net monetary benefit measure,

where $\gamma_{2}-\gamma_{1}$ and $\mu_{2}-\mu_{1}$ denote the true differences in costs and efficacies, respectively, between Treatment 1 and Treatment 2, and $K$ represents the maximum price that a health care provider is willing to pay in order to obtain a unit increase in efficacy, also known as the threshold unit cost. The quantity $\xi$ acts as a measure of cost-effectiveness.

Since the net monetary benefit formula expressed in Equation (17) involves assessing the cost and efficacy components conveyed within each of the two treatment groups, we set $\beta=(\mu_{1},\gamma_{1},\mu_{2},\gamma_{2})^{\top}$, where $\mu_{i}$ and $\gamma_{i}$ denote the efficacy and cost for treatments $i=1,2$. Next, we specify $y_{n}$ as a $4n\times 1$ vector consisting of $2\times 1$ vectors $y_{ij}=(c_{ij},e_{ij})^{\top}$, $i=1,2$ and $j=1,2,\ldots,n$. Each individual observation is allotted one row in the linear model. The design matrix $X_{n}$ is a $4n\times 4$ block diagonal with the $n\times 1$ vector of ones, $1_{n}$, as the blocks. With $n=n_{1}=n_{2}$, $\sigma^{2}=\sigma_{1}^{2}=\sigma_{2}^{2}$ and $\tau^{2}=\tau_{1}^{2}=\tau_{2}^{2}$, our variance matrix is $\displaystyle\sigma^{2}\underset{4n\times 4n}{V_{n}}=\sigma^{2}\begin{pmatrix}I_{n}&O&O&O\\ O&\frac{\tau^{2}}{\sigma^{2}}I_{n}&O&O\\ O&O&I_{n}&O\\ O&O&O&\frac{\tau^{2}}{\sigma^{2}}I_{n}\\ \end{pmatrix},$ where $\sigma^{2}$ is factored out to comply with our conjugate linear model formulation expressed in Equation (4).

In the analysis stage, we use the posterior distribution for $\beta$ if $\sigma^{2}$ is fixed or for $\{\beta,\sigma^{2}\}$ if $\sigma^{2}$ needs to be estimated; recall Sections 2.3.1 and 2.3.2. The posterior distribution is needed only for the analysis stage, hence it is computed using the analysis priors. Since there is no data in the design stage, there is no posterior distribution. We use the design priors as specifications for the population from which the data is generated. That is, the design priors yield the sampling distribution $y_{n}\,|\,\sigma^{2(d)}$. O’Hagan and Stevens, (2001) define $\mu_{\beta}^{(d)}=(5,6000,6.5,7200)^{\top}$ and $\displaystyle V_{\beta}^{(d)}=\begin{pmatrix}4&0&3&0\\ 0&10^{7}&0&0\\ 3&0&4&0\\ 0&0&0&10^{7}.\end{pmatrix}$. We factor out $\sigma^{2}$ from $V_{\beta}^{(d)}$ to be consistent with the conjugate Bayesian formulation in Section 2.3.1 so $\displaystyle\sigma^{2}V_{\beta}^{(d)}=\sigma^{2}\begin{pmatrix}4/\sigma^{2}&0&3&0\\ 0&10^{7}/\sigma^{2}&0&0\\ 3&0&4/\sigma^{2}&0\\ 0&0&0&10^{7}/\sigma^{2}\end{pmatrix}$. Lastly, we set the posterior probability of deciding in favor of $H$ to at least $0.975$, which is equivalent to a Type-I error of $\alpha=0.025$ in frequentist two-sided hypothesis tests.

Consider designing a trial to evaluate the cost effectiveness of a new treatment with an original treatment. We seek the tenability of $H:\xi>0$, where $\xi$ is the net monetary benefit defined in Equation (17). Using the inputs specified in Section 3.1 we execute simulations in the known and unknown $\sigma^{2}$ cases to emulate O’Hagan and Stevens, (2001) as a special case.

We now consider a few alternate Bayesian approaches for sample size determination and demonstrate how these methods can be embedded within the two-stage Bayesian framework. We also identify special cases that overlap with the frequentist setting.

Adcock, (1997) constructs rules based on a fixed precision level $d$. In the frequentist setting, if $X_{i}\sim N(\theta,\sigma^{2})$ for $i=1,...,n$ observations and variance $\sigma^{2}$ is known, the precision can be calculated using $d=z_{1-\alpha/2}\frac{\sigma}{\sqrt{n}}$, where $z_{1-\alpha/2}$ is the critical value for the $100(1-\alpha/2)\%$ quartile of the standard normal distribution. Simple rearranging leads to following expression for sample size,

Given a random sample with mean $\bar{x}$, suppose the goal is to estimate population mean $\theta$. The analysis objective entails deciding whether or not the absolute difference between $\bar{x}$ and $\theta$ falls within a margin of error no greater than $d$. Given data $\bar{x}$ and a pre-specified confidence level $\alpha$, the assurance can be formally expressed as