Sample size determination is a vital part of clinical studies where cost and safety concerns lead to greater importance of not using more subjects and resources than are required. The Bayesian approach to sample size determination has the advantages of being able to use prior data and expert opinion to possibly reduce the total sample size while also acknowledging all uncertainty at the design stage. We apply a Bayesian decision theoretic approach to the problem of comparing two Poisson rates and find the required sample size to obtain a desired power while controlling the Type I error rate.

Keywords: bayes factor, poisson rate, power, Type I error

Sample size determination continues to be an important research area of statistics. Perhaps nowhere is this truer than pharmaceutical statistics, where cost and time constraints have made finding the appropriate sample size before conducting a study of the utmost importance. The problem is quite simple: too small a sample can lead to under-powered studies, while too large a sample size wastes precious resources. In this article we consider the sample size determination problem as it pertains to the two-sample testing of Poisson rates from a Bayesian perspective subject to operating characteristics constraints.

There are several advantages to the Bayesian perspective when trying to determine a study's requisite sample size, a topic that is expounded in Adcock.¹S Their construction does not depend on asymptotic approximations or bounds. Classical solutions to the sample size determination problem typically hinge on asymptotic arguments that require the researcher to specify one parameter value (perhaps vector valued) as representative for the entire parameter space, a process that is typically done using bounding arguments. This, for example, is what is done when determining the requisite sample size for a confidence interval of a fixed level and given length. The resulting sample sizes are consequently conservative. On the other hand, the Bayesian approach provides the statistician with the ability to model his indecision about the parameter through expert knowledge or previous studies. As noted by Bayarri and Berger,² this can allow the Bayesian approach to have better operating characteristics, such as a smaller required sample sizes or better Type I and II error rates.

Various Bayesian sample size determination methods have been studied for binomial and Poisson data. Stamey et al.³considered one and two sample Poisson rates from the perspective of interval based criteria such as coverage and width. Hand et al.⁴ extend those ideas by considering both interval-based and test-based criteria, albeit without considering power. Katsis and Toman⁵ used more decision theoretic criteria for the two sample binomial case, but only to the extent of controlling the posterior risk with a pre specified bound. Zhao et al.⁶ extend those ideas by using computational methods to consider expected Bayesian power of the test. In this article, we extend these results to the Poisson data model. We also consider the problem the subject to Type I and Type II constraints. This is thus an important extension of,⁶ because it

• Extends the ideas to the Poisson case
• Enables the incorporation of operating characteristics.

A subtle difference between the classical and Bayesian methods of sample size determination merits discussion before proceeding. One of the novel contributions of this article is an algorithmic solution to the sample size determination problem subject to operating characteristics constraints for Poisson data. However, since the entire problem is treated in a Bayesian context, the concept of Type I and Type II error rates is understood in an average, or expected, sense; see.⁷ For example, the power of a test'' retains the interpretation of the probability the decision rule rejects when the null hypothesis when it is false; but rather than being a function defined over the alternative space, here it is averaged over that space and weighted by the prior distribution specified on the alternative hypothesis. To make the distinction more clear, we refer to this as the expected Bayesian power (EBP), as is done in;⁸ alternatively, it may be referred to as the probability of a successful test. These ideas, though apparently not considered in the literature previously, can also apply to the concept of the significance level of a test. While frequentist methods typically report one value for significance level, what they are really doing (in non point null hypotheses) is taking the largest possible significance level; thus, taking an expectation of a significance level curve could be done as well. Consequently, in this article we also consider the expected Bayesian significance level (EBSL), defined as the expected value of the test under the prior distribution given on the null space. In both cases, any particular instance of the actual Type I and Type II error rates can be greater than or less than nominal.

This article proceeds as follows. In Section 2 we introduce the theoretical formulation of the sample size determination problem for two Poisson variates including consideration of operating characteristics. In Section 3, we present an algorithmic solution to the sample size determination problem posed in Section 4. Section 5 contains an application of the method in the area of pharmaceutical statistics. We then conclude with a discussion.

Problem specification and the bayes rule
We now follow the general framework of⁶ in the development of this problem, adapting the binomial case to fit the Poisson data model. Suppose $Y_1~Pois(t{\lambda }_1)$  and $Y_2~Pois(t{\lambda }_2)$ , independently, where ${\lambda }_1$  and ${\lambda }_2$  represent the rate parameters of interest and t represents a common sample (or “opportunity”) size. Together, we write these $Y=(Y_1,Y_2{)}^{\prime }$  with observations $y=(y_1,y_2{)}^{\prime }$ . The sample size determination problem is to calculate the necessary sample size required to test the hypotheses
$H_0:{\lambda }_1={\lambda }_2$ vs   $H_1:{\lambda }_1\ne {\lambda }_2$                                                            
using a given decision rule; here we use the Bayes rule. Denoting the parameter pair $\lambda =({\lambda }_1,{\lambda }_2{)}^{\prime }$ , the associated null and alternative spaces are therefore ${\Lambda }_0=\{\lambda \in R_{+}^2:{\lambda }_1={\lambda }_2\},$  which we identify with $R_{+}$  with elements generically denoted $\lambda$ , and ${\Lambda }_1=\{\lambda \in R_{+}^2:{\lambda }_1\ne {\lambda }_2\}$ .
As this problem is being considered from the Bayesian perspective, we place prior probabilities of ${\pi }_0$  and ${\pi }_1=1-{\pi }_0$  on $H_0$  and $H_1$ , respectively. Conditional on the null being true, we represent the expert opinion regarding $\lambda$  as $p_0(\lambda )$ , defined over the set ${\Lambda }_0$ . Alternatively, conditionally on $H_1$  being true, we represent the belief concerning ${\lambda }_1\ne {\lambda }_2$  with a joint prior $p({\lambda }_1,{\lambda }_2)$  with support ${\Lambda }_1$ . Marginalizing over the hypotheses, we have the unconditional prior
$p({\lambda }_1,{\lambda }_2)=p_0(\lambda ){\pi }_0{I}_{H_0}+p_1({\lambda }_1,{\lambda }_2){\pi }_1{I}_{H_1},$ Where $I_{H_0}$  and $I_{H_1}$  are the indicator function of $H_0$  and $H_1$ , respectively.
In practice, the prior distributions on $\lambda$ , ${\lambda }_1$ , and ${\lambda }_2$  summarize expert opinion concerning the parameters in each of the two scenarios $H_0$  and $H_1$ . This information can be obtained in a number of ways including past data, prior elicitation of expert opinion (see especially⁹), or based on uninformative criteria. For simplicity, we consider conjugate priors for all three $\lambda$ ’s, so that under $H_0$ , $\lambda \sim \text{Gamma}(\alpha ,\beta )$ , and under $H_1$ , ${\lambda }_i\sim \text{Gamma}({\alpha }_i,{\beta }_i)$  independently. This assumption is not very restrictive and allows us to specify parameters of prior distribution sin stead of distributions themselves.
We now derive the optimal Bayes (decision) rule in deciding between the hypotheses presented in (1). In solving the related problem between two binomial proportions,⁶ use the classical decision theoretic setup using the 0-1 loss function.¹⁰ Here we use the more general unequal loss function
$L(H,a)=\{\begin{array}{c}{c}_1\text{if }{H}_0\text{ is true and }a\text{ = 1,}\\ c_2\text{if }{H}_1\text{ is true and }a\text{ = 0,}\\ 0,\text{otherwise }\end{array}$
Where $a=\delta (y)$  is the decision rule with a = 0 representing selection of $H_0$  and a = 1 selection of $H_1$ . Thus, c1 represents the loss associated with a Type I error, and c2 that of a Type II error.
The Bayes action is simply the one that minimizes posterior expected loss.¹¹ Since the posterior expected loss associated with an action $a\in \{0,1\}$  is
$$\begin{array}{l}\rho (a)=E_{\lambda |y}[L(\phi ,a)]\\ =\{\begin{array}{c}{c}_{1}P(H_0|Y=y),\text{ if }a=1,\\ c_{2}P(H_1|Y=y),\text{ if }a=0,\end{array}\end{array}$$ Setting $c=c_1/c_2$  we can express the optimal decision rule, $a^{*}$ , as
$a^{*}(y)=\{\begin{array}{c}0,\text{ if }P(H_1|Y=y)<cP(H_0|Y=y),\\ 1,\text{ if }P(H_1|Y=y)\ge cP(H_0|Y=y).\end{array}$

The rejection region, $W$ , is therefore
$W=\{y:P(H_1|Y=y)\ge cP(H_0|Y=y)\}.$ The optimal rule in (2) can be nicely represented in terms of the Bayes factor. The Bayes factor is defined as the ratio of the posterior odds of $H_1$  to the prior odds of $H$ , so that a large Bayes factor is evidence for rejecting the null hypothesis.¹² Specifically, the Bayes factor is defined

This ratio is useful in Bayesian inference because it is often interpreted as partially eliminating the influence of the prior on the posterior, instead emphasizing the role of the data. Moreover, the decision rule is a function of a Bayes factor:
$\begin{array}{l}W=\left\{y:P(H_1|Y=y)\ge cP(H_0|Y=y)\right\}\\ =\left\{y:P(H_1|Y=y)\frac{{\pi }_0}{{\pi }_1}\ge cP(H_0|Y=y)\frac{{\pi }_0}{{\pi }_1}\right\}\end{array}$ = $$=\left\{y:B\ge c\frac{{\pi }_0}{{\pi }_1}\right\}.$$

This is particularly useful because it allows for the interpretation of the Bayes factor B as the test statistic for the decision rule in (2); this is the condition in (4).

We now derive closed-form expression for the posterior probabilities of the null and alternative hypotheses. Using Bayes’ theorem, we have

Note that the probability of the data under the null hypothesis is the product of two independent negative binomial likelihoods.
Combining (3) with (5) and (6) allows us to represent W in terms of the null and alternative likelihoods as follows:

Consequently, (7) and (8) give explicit conditions for the rejection of the optimal decision rule:
$W=\left\{y:B=\frac{\Gamma (\alpha ){(2t+\beta )}^{y_1+y_2+\alpha }}{{\beta }^{\alpha }\Gamma (y_1+y_2+\alpha )}{\displaystyle \prod _{i=1}^2\frac{{\beta }_i{}^{{\alpha }_i}\Gamma (y_i+{\alpha }_i)}{\Gamma ({\alpha }_i){(t+{\beta }_i)}^{y_i+{\alpha }_i}}}\ge c\frac{{\pi }_0}{{\pi }_1}\right\}$
Note that the left side of (9) is our test statistic and Bayes factor, B, so that (9) is an explicit formulation of the condition presented in (4).

Sample Size Determination for the Bayes Rule
The explicit description of the decision rule in (9) allows us to compute all sorts of quantities of interest. For given prior parameters ${\pi }_0$ , ${\pi }_1$ , α, β, ${\alpha }_1$ , ${\beta }_1$ , ${\alpha }_2$ , and ${\beta }_2$  and loss penalties $c_1$ and $c_2$ (or simply c), the Expected Bayesian Power (EBP) ${\omega }_t$ is defined
${\omega }_t=P(Y\in W|H_1)={\displaystyle \sum _{y\in W}P(Y=y|H_1})={\displaystyle \sum _{y\in W}{\displaystyle \prod _{i=1}^2\frac{t^{y_i}{\beta }_i{}^{{\alpha }_i}\Gamma (y_i+{\alpha }_i)}{y_i!\Gamma ({\alpha }_i){(t+{\beta }_i)}^{y_i+{\alpha }_i}}}},$
and the Expected Bayesian Significance Level (EBSL)  ${\alpha }_t$ is
${\alpha }_t=P(Y\in W|H_0)={\displaystyle \sum _{y\in W}P(Y=y|H_0})={\displaystyle \sum _{y\in W}\frac{t^{y_1+y_2}{\beta }^{\alpha }\Gamma (y_1+y_2+\alpha )}{y_1!y_2!\Gamma (\alpha ){(2t+\beta )}^{y_1+y_2+\alpha }}}.$

Note three things. First the inclusion of the t subscripts highlights the fact that these quantities depend on t. Second, both ${\omega }_t$  and ${\alpha }_t$ marginalize over the corresponding alternative and null spaces ${\Lambda }_1$  and ${\Lambda }_0$ , respectively, this is the sense in which the power and significance level are expectations. Third, the constant c (or c1 and c2) is represented in the expressions through Wt, which is itself dependent on t.

In their articles,^5,12 demonstrate that as the sample size tends to infinity, the Bayes factor  converges to either 0 or 1. As a consequence, in the current context as the sample size t tends to infinity the Bayes factor $B$  converges to either 0 or 1, so that¹² implies that ${\omega }_t\stackrel{a.s.}{\to }1$  and ${\alpha }_t\stackrel{a.s.}{\to }0$ . Thus, for any pre-specified power $\omega$  and significance level $\alpha$ , there exists a t such that for all $t^{\prime }\ge t^{*}$ , ${\omega }_{t^{\prime }}\ge \omega$  and ${\alpha }_{t^{\prime }}\ge \alpha$ . We define $t_{\alpha ,\omega }^{*}$  to be the in fimum of this collection of lower bounds, i.e.

$t_{\alpha ,\omega }^{*}$ is said to be the optimal sample size for $\omega$ and $\alpha$ , and computing $t^{*}$  is called the sample size determination problem. If only a power is specified, or if only a significance level is specified, the other quantity is left off of the subscript and out of the definition. We often write simply $t^{*}$  for $t_{\alpha ,\omega }^{*}$ .

Practically speaking, in our experience the appearance of these jumps is monotonically decreasing in magnitude and dissipate quite quickly in t so that, while there are jumps, they become relatively minor even for quite small t. Thus, while out-of-the-box numerical routines are insufficient for the task, a straight-forward heuristic algorithm suffices to certify $t^{*}$  to a reasonable level of accuracy.

Initialize $t=t_0$ with a value small enough to satisfy ${\omega }_t<\omega$  and ${\alpha }_t>\alpha$ , typically a value like .1 will suffice. Then, double t until both conditions are met. Once both conditions are met, decrement t by 1 until the conditions are no longer satisfied; then decrement by .1, and so on to achieve the desired precision. Alternatively, one may move in a binary search manner; this method is faster but loses the certificate of the solution up to the highest evaluated point. By contrast, the first heuristic certifies that the sample size achieved is optimal up to the largest t observed, which is of the form $2^k{t}_0$ .

One computational detail is relevant for implementing this procedure. By definition, since $Y_1$  and $Y_2$  are Poisson variates, their sample spaces are infinitely large, and thus computing ${\omega }_t$  and ${\alpha }_t$  is not numerically possible - one cannot check every $y$  to determine whether or not it is included in Wt. There is, however, a very reasonable work-around for this problem using prior predictive distributions. Under H1, the prior predictive distributions on $Y_1$  and $Y_2$  are both negative binomial. Specifically,
$Y_1~\text{NegBin}\left({\alpha }_1,\frac{{\beta }_1}{t+{\beta }_1}\right)$ and
$Y_2~\text{NegBin}\left({\alpha }_2,\frac{{\beta }_2}{t+{\beta }_2}\right).$

While one cannot enumerate the entire sample space of , one can be confident they have a satisfactory approximation by computing small (e.g. .00001) and large (e.g. .99999) quintiles of these distributions and then simply taking every combination of the ranges from low to high. This is the approach we take to computing the sums listed in (10) and (11) above.

An Example from Cancer Therapeutics

1. An example of the proposed methodology is readily available from the field of cancer therapeutics. Suppose (1) Drug A is an industry standard therapy for a certain type of cancer that is known to have the common side effect of mild seizures every hour of infusion
2. (2) Drug B is a novel compound believed to have the same side effect but at a lessened rate
3. (3) The goal is to design a design a clinical trial that compares the two using the minimum resources required to meet 5% significance level (EBSL) and 80% power (EBP). Moreover, suppose that the losses associated with Type I and Type II errors have a ratio of 1:1.

From past studies, it is known that the uncertainty in the rate of seizures with Drug A (per hour) is well-represented by a Gamma (4, 4) distribution so that ${\lambda }_A~\text{Gamma}(4,4)$ . Drug B, by contrast, is believed to be a bit worse, with perhaps a rate that is double that of Drug A with experts 90% sure the value is less than about 3. This translates into roughly ${\lambda }_B~\text{Gamma}(4,8)$ . Figure 1 shows both of these priors graphically.

Figure 1 Prior structures used in Poisson sample size determination example.

Assuming that the rates are the same, the past indicators of Drug A supersede the lack of evidence for Drug B, so that under the null hypothesis $\lambda ={\lambda }_1={\lambda }_2~\text{Gamma}(4,4)$  is most appropriate. Assuming that the null and alternative hypotheses are given the same belief (50%), the rejection rule from (9) is therefore

$$W=\left\{y:\frac{\Gamma (4){(2t+4)}^{y_1+y_2+4}}{4^4\Gamma (y_1+y_2+4)}\frac{4^4\Gamma (y_1+4)}{\Gamma (4){(t+4)}^{y_1+4}}\frac{4^8\Gamma (y_2+8)}{\Gamma (8){(t+4)}^{y_1+8}}\ge 3\right\}$$ (12)
$$=\left\{y:\frac{{(2t+4)}^{y_1+y_2+4}}{{(t+4)}^{y_1+4}}\ge \frac{15120}{65536}\frac{(y_1+y_2+3)!}{(y_1+3)!(y_2+7)!}\right\}$$(13)

To design a test with 80% power and at the 5% significance level, we need but to run the algorithm described in Section 3. To illustrate the scenario, we plot the significance level and power for every t from 1 to 80; these are included in Figures 2 and Figure 3, respectively. Note how quickly (in t) the functions become smooth, nullifying any concern about jumps. The 80% level is achieved at t = 37, with a power of 80.1%, and the 5% significance level is achieved at t = 57, where the significance level is 4.9%. Thus, to achieve both, we select the higher sample size, t = 57.

Figure 2 EBP curve.

Figure 3 EBSL curve.

We can also verify these results via simulation by generating one million random values of ${\lambda }_1$ and ${\lambda }_2$  from the priors given above. To validate the significance level result, we simply (1) sample from the prior, (2) generate two Poisson observations with mean equal to the sample in (1) times 57, and (3) determine whether the test rejects (1) or not (0). Averaging the million test results, the value is confirmed to Monte Carlo error (code available upon request). To validate the power result, we (1) sample one number from a Gamma(4, 4)and one number from a Gamma(8, 4), (2) sample two Poisson observations from distributions with mean 37 times the two variates generated in (1), and (3) determine whether the test rejects or not. Averaging the million results validates the theoretical result.

None.

Authors declare that there are no conflicts of interests.

1. Adcock CJ. Sample Size Determination: A Review. The Statistician. 1997;46(2):261–283.
2. Bayarri MJ, Berger JO. The Interplay of Bayesian and Frequentist Analysis. Statistical Science. 2004;19(1):58–80.
3. Stamey JD, Young DM, Bratcher TL. Bayesian Sample-Size Determination for One and Two Poisson Rate Parameters with Applications to Quality Control. Journal of Applied Statistics. 2006;33(6):583–594.
4. Hand A, Stamey JD, Young DM. Bayesian Sample-Size Determination for Two Independent Poisson Rates. Computer Methods and Programs in Biomedicine. 2011;104(2):271–277.
5. Katsis A, Toman B. Bayesian Sample Size Calculations for Binomial Experiments. Journal of Statistical Planning and Inference. 1992;81(2):349–362.
6. Zhao Z, Tang N, Li Y. Sample-Size Determination for Two Independent Binomial Experiments. Journal of Systems Science and Complexity. 2011;24(5):981–990.
7. O'Hagan A, Stevens J, Campbell M. Assurance in Clinical Design. Pharmaceutical Statistics. 2005;4(3):187–201.
8. Speighelhalter DJ, Abrams K R, Myles JP. Bayesian Approaches to Clinical Trials and Health Care Evaluation. Wiley. 2004.
9. Garthwaite PH, Kadane JB, O'Hagan A. Statistical methods for eliciting probability distributions. Journal of the American Statistical Association. 2005;100:680–701.
10. Shao J. Mathematical statistics. 2nd ed. Springer; 2003.
11. Berger J. Statistical decision theory and Bayesian analysis. Springer; 1985.
12. Kass R, Raftery A. Bayes Factors. Journal of the American Statistical Association. 1995;90(430):773–795.