- Email: [email protected]
Bayesian sample-size determination for two independent Poisson rates
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
journal homepage: www.intl.elsevierhealth.com/journals/cmpb
Bayesian sample-size determination for two independent Poisson rates Austin L. Hand ∗ , James D. Stamey, Dean M. Young Department of Statistics, Baylor University, Waco, TX, USA
a r t i c l e
i n f o
a b s t r a c t
Article history:
Because of the high cost and time constraints for clinical trials, researchers often need to
Received 27 March 2010
determine the smallest sample size that provides accurate inferences for a parameter of
Received in revised form
interest. Although most experimenters have employed frequentist sample-size determina-
7 October 2010
tion methods, the Bayesian paradigm offers a wide variety of sample-size determination
Accepted 14 October 2010
methodologies. Bayesian sample-size determination methods are becoming increasingly more popular in clinical trials because of their flexibility and easy interpretation inferences.
Keywords:
Recently, Bayesian approaches have been used to determine the sample size of a single
Bayesian sample-size determination
Poisson rate parameter in a clinical trial setting. In this paper, we extend these results to
Comparison of two Poisson rates
the comparison of two Poisson rates and develop methods for sample-size determination
Average power
for hypothesis testing in a Bayesian context. We have created functions in R to determine
Average length criterion
the parameters for the conjugate gamma prior and calculate the sample size for the average length criterion and average power methods. We also provide two examples that implement our sample-size determination methods using clinical data. © 2010 Elsevier Ireland Ltd. All rights reserved.
1.
Introduction
In clinical trials using count data of rare events, the Poisson distribution is a natural and commonly used model. Often times, interest is in the comparison of two Poisson rate parameters, as in Gu et al. [1], Ng and Tang [2], Thode [3], and Shiue and Bain [4]. In general, Bayesian methods are becoming increasingly more popular in clinical trials because of their flexibility and their ease of interpretation. Because of the high cost and time constraints for clinical trials, researchers often need to determine the smallest sample size n that provides accurate inferences for a parameter of interest. Although most experimenters have employed frequentist sample-size determination methods, the Bayesian paradigm offers a wider variety of sample-size determination methodologies and has been recommended over frequentist methods by Adcock [5]. The main limitation of the frequentist
∗
approach is that one must assume a “known” parameter value or vector. For example, one must assume that the variance or noncentrality parameters are known to perform a sample-size determination calculation for comparing two means. Obvious concerns are how one obtains this “known” parameter value and the consequences on the determined sample size if this “known” parameter is inserted without some incorporation of the uncertainty or variability. In contrast, Bayesian methods not only incorporate uncertainty about unknown parameter values, but also provide a variety of sample-size determination criteria. Two commonly implemented Bayesian sample-size methods are the average length criterion (ALC) method proposed by Joseph et al. [6] and the average power (AP) method proposed by Wang and Gelfand [7]. Recently, Zaslavsky [8] has used prior data to express uncertainty about parameters and an interval-based criterion to determine the sample size for the estimation of a Poisson rate. In this paper we extend Zaslavsky’s sample-size determina-
Corresponding author. Tel.: +1 254 315 5854. E-mail address: Austin [email protected] (D.M. Young). 0169-2607/$ – see front matter © 2010 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2010.10.010
272
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
tion method to the comparison of two Poisson rates, develop methods for sample-size determination via the AP method, and expand methods described by Stamey et al. [9] for the ALC method. Although a fully Bayesian approach requires the specification of a loss or utility function, the ALC method uses a sample-size criterion without specifying either, thus making the ALC method simpler to apply than other sample-size determination methods. Adcock [5] has noted that one may encounter considerable difficulty in formulating a loss or utility function that applies to all possible sample data. For the ALC method we fix the credible level, ˛, and determine the sample size n such that, on average, the length of the highest posterior density (HPD) interval is equal to a predetermined value l. The AP method uses the average posterior power to determine a suitable sample size n. This method also requires no specification of a loss or utility function. Although the AP method may be computationally intensive, it is very useful for the hierarchical models that dominate the current Bayesian landscape. Because of the potentially large influence the prior may have on the posterior distribution for small samples, choosing a prior is a crucial step in Bayesian modeling. However, the prior should be justified and, therefore, the influence it has on the posterior should be thought of as a valid source of uncertainty incorporated with the parameter. A desirable trait of a prior distribution is that it allows researchers to include available information, or lack thereof, with a degree of uncertainty. One common practice is to base the prior parameters choice on previous studies or trials. This prior selection approach is commonly referred to as objective Bayesian inference. In this paper we formulate an R function that uses prior information from previous clinical trials to define the parameters of the conjugate gamma prior through the implementation of marginal maximum likelihood estimation that uses the R function optim. Using these parameters, we also develop R functions that provide the required sample size for the ALC and AP methods for the ratio of two Poisson rates. The remainder of the paper is organized as follows. In Section 2 we provide background information on the Poisson distribution as well as review the marginal maximum likelihood estimation method for the gamma prior parameters. In Section 3 we propose and derive sample-size determination methods for the ratio of two independent Poisson rate parameters. We continue in Section 4 with two real-life examples for sample-size determination using the conjugate gamma prior. Finally, in Section 5 we conclude with brief comments.
2. Marginal maximum likelihood estimation 2.1.
Xi ∼ Poisson(ni i ) so that
(1)
where xi = 0, 1, 2, . . . and i > 0. Although we can use one of many different priors for i such as the uniform or log-normal, we use the conjugate gamma prior. Therefore,
(i |ai , bi ) =
bai i (ai )
iai −1 e−bi i ,
(2)
where 0 < i < ∞ and ai , bi > 0 with mean ai /bi and variance ai /b2i . Thus, the posterior density of i is a +x
(ni + bi ) i i ai +xi −1 −i (ni +bi ) e . (ai + xi ) i
(i |xi , ai , bi , ni ) =
(3)
Also, the marginal density of Xi is
m(xi |ai , bi ) =
∞
P(Xi = xi |i )(i |ai , bi )di 0
=
(ai + xi ) (ai ) (xi + 1)
b ai n xi i i ni + bi
ni + bi
,
(4)
which is a negative binomial probability function. Given k independent previous Poisson trials with xij events each, where i = 1, 2 and j = 1, 2, . . ., k, the marginal likelihood is
M(ai , bi ) =
k
m(xij |ai , bi ).
(5)
j=1
The marginal maximum likelihood estimators (MMLEs) for the parameters ai and bi , i = 1, 2, that we use in (2) are
(aˆ i , bˆ i ) = arg max M(ai , bi )|xi1 , . . . , xik .
2.2.
(6)
Marginal maximum likelihood estimation
To determine (aˆ i , bˆ i ), i = 1, 2, we apply a simplified minimization method developed by Byrd et al. [10] in the R function optim (see Appendix (A)). Because the method developed by Byrd et al. is a minimization instead of maximization routine, we use the negative of the marginal as our criterion function. Thus, the minimized log -likelihood is
Theoretical background
Let Xi be the number of events in ni patient time, and let i be the event rate for population i, i = 1, 2. We assume the model
x
e−ni i (ni i ) i , xi !
P(Xi = xi |i ) =
L(ai , bi ) =
k j=1
− log m(xij |ai , bi ) ,
(7)
273
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
we find the smallest positive value n such that
where j = 1, 2, . . ., n. Hence, the marginal log -likelihood is
L(ai , bi )
k
=−
log
j=1 k
=−
(ai ) (xij + 1)
∞ ∞ 1
b ai n xij i i
(ai + xij )
ni + bi
ni + bi
x2 =0
log (ai + xij ) − log ( (ai )) − log (xij + 1)
j=1
(8) Therefore, the MMLEs are
(aˆ i , bˆ i ) = arg min L(ai , bi )|xi1 , . . . , xik ,
(9)
where L(ai , bi ) is defined in (8).
3.
Sample-size determination criteria
3.1.
Average length criterion approach
The ALC approach, a method proposed by Joseph et al. [6] and discussed by Adcock [5], proposes that for a fixed posterior credible interval coverage, 1 − ˛, we determine the smallest sample size n such that
l(x, n)m(x)dx ≤ l, where l(x, n) is the length of the 100(1 − ˛)% posterior credible interval for the data x and l is a set length. We now consider the case for comparing two independent Poisson rate parameters using the ALC sample-size determination method. That is, we assume Xi ∼ Poisson(ni i ), i = 1, 2. Although several criteria for comparing the two rate parameters are available, one common methodology is the ratio = 1 /2 . As derived by Price and Bonnet [11], the posterior distribution is k∼Fv1 ,v2 ,
(10)
(a2 + x2 )(b1 + n1 ) (a1 + x1 )(b2 + n2 )
with v1 = a1 + x1 and v2 = a2 + x2 . Therefore the ALC solution, as described by Stamey et al. [9], is a sum since both x1 and x2 are discrete. The sum is over all possible outcomes and averages the (1 − ˛)100% posterior intervals, weighted by the marginal distribution. Thus, we have
∞
Fv1 ,v2
− Fv1 ,v2
˛ 2
m(x1 |aˆ 1 , bˆ 1 )
m(x2 |aˆ 2 , bˆ 2 ) ≤ l,
3.2.
(11)
∞
l(x1 , x2 , n)m(x1 |aˆ 1 , bˆ 1 )m(x2 |aˆ 2 , bˆ 2 ).
x2 =0x1 =0
The length of the interval for a particular value of x1 and x2 is computed by finding the (1 − ˛/2) and the ˛/2 percentiles of the F distribution and multiplying by the constant 1/k. That is
Hypothesis testing
We are often interested in testing the difference of the placebo and treatment groups in clinical trials. Lindley [12] proposed a multiplicative model where denotes the rate for the placebo group and denotes the rate for the treatment group. Then, < 1, = 1, and > 1 correspond to an improvement, no change, or a deterioration, respectively, for the treatment group. In general, we are interested in a one-sided alternative hypothesis that represents an improvement or a deterioration for the treatment group. The likelihood is P(, |Xi = xi ) = e−(n1 +n2 ) (x1 +x2 ) x2 ,
(12)
where xi = 0, 1, 2 . . ., i = 1, 2, and > 0. In practice one is often required to have a balance between the placebo and treatment groups. Thus, we assume n1 = n2 = n. Given the likelihood (12) and the conjugate gamma prior (2), ∼ gamma(a1 , b1 ), the marginal posterior distribution for is
P(|Xi = xi )
=
∞
P(, |Xi = xi )(|a1 , b1 )d 0
∝
x2 x +x +a , (n + n + b1 ) 1 2 1
where xi = 0, 1, 2 . . ., i = 1, 2, and > 0. The marginal density, m(x1 | a1 , b1 ), for the placebo group has a closed form as defined in (4). However, the treatment group marginal, m(x2 | a2 , b2 ), does not have a closed form. Hence, we employ numerical integration to approximate the marginal posterior distribution for .
3.3.
where k=
k
˛ 1− 2
where m(xi |aˆ i , bˆ i ), i = 1, 2, refers to the negative binomial marginal distribution (4).
+ ai log (bi ) − ai log
(ni + bi ) + xij log (ni ) − xij log (ni + bi ) .
x1 =0
Average power method
One can find a general algorithm for computing the average power (AP) in Wang and Gelfand [7]. Here, we determine the smallest sample size n given the parameter of interest such that
E I [P( > 0 |x1 , . . . , xn ) > 1 − ˛]
≥ 1 − ˇ,
(13)
where I( · ) is one if H0 is rejected and zero otherwise. Note that (13) is the expected probability of rejecting the null hypothesis. For the calculation of the average power, we first assign a value for the parameter of interest, while assuming that the alternative hypothesis is true, which we denote by ∗ . We then generate values for our treatment and placebo group, where x1 ∼ Poisson(n) and x2 ∼ Poisson(n ∗ ). Letting x1 and x2 be our data, we use Monte Carlo simulation to approximate the pos-
274
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
terior probability P( > 0 | xi ), where i = 1, 2. Using the posterior probability, for those posteriors with P( > 0 | x1 , . . ., xn ) > 1 − ˛ we assign a value of one to an element of the vector, b, and a value of zero, otherwise. By averaging over the vector b, we have an estimation for the average power. This process can be repeated for several different sample-size values, and then we can use a parametric or spline fit to determine the sample size that satisfies (13). In order to perform a Monte Carlo simulation, we must place priors on all parameters. We place a conjugate gamma prior (2) on the rate parameter , ∼ gamma(a1 , b1 ). For the parameter we use a relatively noninformative gamma prior, ∼ gamma(0.01, 0.01), which has a mean of one and is relatively flat around = 1. Thus, the alternative hypothesis will not be preferred.
4. Example: thromboembolism events for the CarboMedics mitral valves 4.1.
Description of Data
Thromboembolism is the formation of a clot (thrombus) in a blood vessel that breaks loose and is carried by the blood stream to plug another vessel. We are specifically interested in thromboembolism in the mitral valve, which is a dual-flap valve that lies between the left atrium and the left ventricle in the heart. A normally-functioning mitral valve opens secondary to contraction from the chordae tendineae, allowing blood to flow into the left ventricle during diastole, and closes at the end of the atrial contraction to prevent blood from backflowing into the atria during left ventricle systole [13]. Data is given in Zaslavsky [8] on thromboembolism when we use CarboMedics mitral valves to replace dysfunctional mitral valves in Table 1. By optimizing the marginal log -likelihood (9) for the thromboembolism data given in Table 1, we determined that the MMLEs for the first four series are â1 = 0.68 and bˆ 1 = 4.22, and the MMLEs for the latter four series are â2 = 1.75 and bˆ 2 = 8.72. Hence, the priors are 1 ∼ gamma(0.68, 4.22) and 2 ∼ gamma(1.75, 8.72), respectively.
Table 1 – Thromboembolism events for the CarboMedics synthetic mitral valves. First four series Events (x) Times (n)
16 43.1
12 58
0 67.7
6 80
Latter four series 16 109.6
30 139.8
64 157.1
6 181.5
4.2. Example: average length criterion sample-size determination for thromboembolism data We desire to compare the two independent Poisson rate parameters. Thus, we wish to determine the smallest sample size n that yields a relatively high probability of finding a difference in the two rates. Using the average length criterion for the two sample case (11), we set the cutoff level to be l = 1. We have developed a function in R that calculates the smallest sample size n for the comparison of two Poisson rates with gamma priors. Using this function, we found that the smallest sample size is n = 2614 0.1 person-years or 261.4 person-years (see Fig. (1)). Code for the ALC method can be obtained from the authors.
4.3. Example: average power sample-size determination for thromboembolism data For the thromboembolism data, we use only the first four series to estimate the prior for the control group. We are interested in the smallest sample size n that allows us to compare old CarboMedics mitral valves with new mitral valves. If the new mitral valves are better than the old CarboMedics mitral valves, then the rate parameter should be less for the new valves, and, hence < 1. The corresponding null and alternative hypotheses would be H0 : = 1 vs.
HA : < 1.
Using WinBugs, we placed a noninformative prior on the parameter , ∼ gamma(0.01, 0.01), along with the conjugate gamma prior on the first four series, ∼ gamma(0.68, 4.22), and let ∗ = 0.8. Using these values, we simulated the average power of our data for various sample sizes n (see Appendix (B)). Because it is common practice to use a minimum power of 0.8,
Fig. 1 – Average length criterion sample-size determination for CarboMedics synthetic mitral valves.
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
275
Fig. 2 – Simulated average power sample-size determination for CarboMedics synthetic mitral valves with priors ∼ gamma(0.68, 4.22) and ∼ gamma(0.01, 0.01).
we determine the smallest sample size n that gives an average power of 0.8. A plot of the power versus sample size n is given in Fig. (2). From this, we can see that the smallest sample size with an average power of 0.8 is n ≈ 1882 0.1 person-years or 188.2 person-years.
5.
Comments
Zaslavsky [8] has developed methods for sample-size determination for a single Poisson rate, which we expanded to the comparison of two independent Poisson rates by developing methods for sample-size determination using the AP method proposed by Wang and Gelfand [7] and expanding methods described by Stamey et al. [9] for sample-size determination via the ALC method. Although many methods for sample-size determination exist, Bayesian methods provide a wider variety and incorporate uncertainty for the parameter of interest with the specification of the prior distribution. Two commonly implemented Bayesian methods are the average length criterion (ALC) and the average power method (AP). One can implement both methods without the specification of a loss or utility function, thus making them easier to apply than competing Bayesian sample-size determination methods. The AP method is generally computationally more rigorous, but it is very useful for the hierarchical models that dominate the current Bayesian literature. We have created functions in R that determine the parameters for the conjugate gamma prior (see Appendix (A)) and calculate and display the sample size for the ALC and AP methods (see Appendix (B)). We compared the different methods available in the optim function in R and found that the method developed by Byrd et al. [10] consistently performed better than the other methods offered. Another reason for choosing this method is that it allows one to specify limits for the parameters of interest. For our situation this is important because we require ai , bi > 0 for the conjugate gamma prior (2). We have also provided two example implementations using the thromboembolism data in Sections 4.2 and 4.3.
Although the ALC method provides an exact sample size, the AP method must be approximated because the marginal distribution of the treatment group has no closed form. A major advantage of the ALC and AP methods are that they provide fast and easy methods to determine an appropriate sample size for the comparison of two independent Poisson rates.
Acknowledgements The authors would like to thank the Associate Editor and referee for comments that improved the manuscript.
Appendix A. Marginal maximization likelihood estimator function MMLE=function(x,n){ # define the negative log likelihood loglike=function(parameters){-(sum( lgamma(x+parameters[1]))-length(x)* lgamma(parameters[1])-sum(lgamma(x+1)) +length(x)*parameters[1]*log(parameters[2])sum(parameters[1]*log(parameters[2]+n))+ sum(x*log(n))-sum(x*log(n+parameters[2])))} # load the library stats for optim function library(stats) # set the lower limit for the parameters of # interest, if the function returns the error # given below, adjust this value # lowerlim=0.000000000000000001 # optimize the negative log likelihood function # with the function optim; many methods are # available, we use a method developed by # Byrd et al. that determines the minimum opt=optim(c(100,100),loglike,lower=lowerlim, method=”L-BFGS-B”) # check that the lower limit is set small # enough if(opt$par[1]==lowerlim){return(”The lower
276
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
limit in the optim function is not small enough.”)} if(opt$par[2]==lowerlim){return(”The lower limit in the optim function is not small enough.”)} # return the optimized value opt$par}
Appendix B. Average power sample-size determination # function that calculates the power for given # alpha and beta values power=function(a,b,a.bugs,b.bugs,theta,n,m,n. iter,greater){ # load the package R2WinBUGS to implement MCMC # methods library(R2WinBUGS) POW=seq(0,0,length=m) Post.prob=vector() # specifies whether testing greater than or # less than in the alternative hypothesis, # if greater is 1 then we are interested in # testing greater than, and 0 otherwise if(greater==1){parameters=list(”post.prob. great”)} if(greater==0){parameters=list(”post.prob. less”)} # for loop that runs m iterations to calculate the average power # for given alpha and beta values for(i in 1:m){ lambda=rgamma(1,a,rate=b) # generate a random value for our placebo and # treatment groups, respectively y1=rpois(1,n*lambda) y2=rpois(1,n*lambda*theta) data=list(”y1”,”y2”,”n”,”a.bugs”,”b. bugs”) pow=bugs(data,inits=NULL,parameters,”powerbugs. txt”,n.chains=1,n.burnin=1000,n.iter=n.iter) Post.prob[i]=pow$summary[1,1] # assign a value of 1 if the posterior # probability is greater than 0.95 and a # value of 0 otherwise, vector b in text. if(Post.prob[i]>0.95){POW[i]=1}} V # the average power for given alpha and beta values Power=mean(POW) return(Power)} # function that calculates the smallest sample # size for given alpha or beta values POWpoint=function(a,b,a.bugs,b.bugs,theta,m,n. iter,greater,beta0){ n=c(1,2000) POW=c(0,0) # check to see that the power chosen is larger # than the power for n=1, if not displays error # message
POW[1]=power(a,b,a.bugs,b.bugs,theta,n[1],m,n. iter,greater) if(POW[1]>=(1-beta0)){return(“The power chosen is too small. The sample size n=1 has a power greater than one minus beta.”)} POW[2]=power(a,b,a.bugs,b.bugs,theta,n[2],m,n. iter,greater) # increases n if n=2000 is too small for the # specified power value while(POW[2]<(1-beta0)){ n[2]=2*n[2] POW[2]=power(a,b,a.bugs,b.bugs,theta,n[2],m,n. iter,greater)} nmax=n[2] # calls the function Powplot to create a plot # for the power given specified alpha and # beta values, it also returns coefficients # for sample-size calculation coeff=Powplot(a,b,a.bugs,b.bugs,theta,m,n.iter, greater,beta0,nmax) # uses coefficients from a polynomial fit to calculate the # smallest sample size for specified power while((n[2]-n[1])>1){ N=ceiling((n[2]+n[1])/2) POW2=coeff[1]+N*coeff[2]+N ∧2*coeff[3]+N ∧3*coeff[4] if(POW2>=(1-beta0)){n[2]=N}else{n[1]=N}} # returns the smallest sample size calculated return(n[2])} # function that creates a power plot # for specified alpha and beta values POWplot=function(a,b,a.bugs,b.bugs,theta,m,n. iter,greater,beta0,nmax){ n=seq(50,nmax,by=50) powplot=vector() for(i in 1:length(n)){ # calls power function to calculate the power # for specified alpha and beta values powplot[i]=power(a,b,a.bugs,b.bugs,theta,N,m,n. iter,greater)} # fits a third degree polynomial to the power # plot fit=lm(powplot∼n+I(n ∧2)+I(n ∧3)) plot(n,powplot,xlab=”Sample Size n”, ylab=”Power”,main=”Average power sample-size determination for CarboMedics synthetic mitral valves”) lines(n,fit$coeff[1]+n*fit$coeff[2]+n ∧2*fit$coeff[3]+n ∧3*fit$coeff[4]) abline(h=1-beta0) # returns the coefficients from the third # degree polynomial fit return(fit$coeff)} # winbugs model for sample-size calculation model{ # specify the poisson distribution for the # placebo and treatment group, # respectively
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
y1 ∼dpois(lambda1) y2 ∼dpois(lambda2) # specify the noninformative gamma prior theta ∼dgamma(0.01,0.01) # specify the conjugate gamma prior lambda ∼dgamma(a.bugs,b.bugs) # parameters for the placebo and treatment # group, respectively lambda1<-n*lambda lambda2<-n*lambda*theta post.prob.less<-1-step(theta-1) post.prob.great<-1-step(1-theta)} Similar code for the ALC method can be obtained from authors.
references
[1] K. Gu, H.K.T. Ng, M.L. Tang, W.R. Schucany, Testing the ratio of two Poisson rates, Biometrica 50 (2) (2008) 283–298. [2] H.K.T. Ng, M.-L. Tang, Testing the equality of two Poisson means using the rate ratio, Statistics in Medicine 24 (2005) 955–965. [3] H. Thode, Power and sample size requirements for tests of differences between two Poisson rates, The Statistician 46 (1997) 227–230. [4] W.K. Shiue, L.J. Bain, Experiment size and power comparisons for two-sample Poisson tests, Applied Statistics 31 (2) (1982) 130–134.
277
[5] C. Adcock, Sample size determination: a review, The Statistician 46 (1997) 261–283. [6] L. Joseph, D.B. Wolfson, R. Berger, Some comments on Bayesian sample size, The Statistician 44 (2) (1995) 167–171. [7] F. Wang, A.E. Gelfand, A simulation-based approach to Bayesian sample size determination for performance under a given model and for separating models, Statistical Science 17 (2) (2002) 193–208. [8] B.G. Zaslavsky, Empirical Bayes models of Poisson clinical trials and sample size determination, Pharmaceutical Statistics (2009). [9] J.D. Stamey, D.M. Young, T.L. Bratcher, Bayesian sample-size determination for one and two Poisson rate parameters with applications to quality control, Journal of Applied Statistics 33 (6) (2006) 583–594. [10] R. H. Byrd, C. Zhu, P. Lu, J. Nocedal, L-BFGS-B—Fortran Subroutines for Large-Scale Bound Constrained Optimization, Technical Report, ACM Trans. Math. Software, (1994). [11] R. Price, D. Bonnett, Estimating the ratio of two Poisson rates, Computational Statistics and Data Analysis 34 (2000) 345–356. [12] D.V. Lindley, A holistic analysis of Poisson data with application to a trial of selenium and cancer deaths, The Indian Journal of Statistics 64 (3) (2002) 722–739. [13] W. Carey, Current Clinical Medicine 2009: online + print, Saunders/Elsevier, (2009).
journal homepage: www.intl.elsevierhealth.com/journals/cmpb
Bayesian sample-size determination for two independent Poisson rates Austin L. Hand ∗ , James D. Stamey, Dean M. Young Department of Statistics, Baylor University, Waco, TX, USA
a r t i c l e
i n f o
a b s t r a c t
Article history:
Because of the high cost and time constraints for clinical trials, researchers often need to
Received 27 March 2010
determine the smallest sample size that provides accurate inferences for a parameter of
Received in revised form
interest. Although most experimenters have employed frequentist sample-size determina-
7 October 2010
tion methods, the Bayesian paradigm offers a wide variety of sample-size determination
Accepted 14 October 2010
methodologies. Bayesian sample-size determination methods are becoming increasingly more popular in clinical trials because of their flexibility and easy interpretation inferences.
Keywords:
Recently, Bayesian approaches have been used to determine the sample size of a single
Bayesian sample-size determination
Poisson rate parameter in a clinical trial setting. In this paper, we extend these results to
Comparison of two Poisson rates
the comparison of two Poisson rates and develop methods for sample-size determination
Average power
for hypothesis testing in a Bayesian context. We have created functions in R to determine
Average length criterion
the parameters for the conjugate gamma prior and calculate the sample size for the average length criterion and average power methods. We also provide two examples that implement our sample-size determination methods using clinical data. © 2010 Elsevier Ireland Ltd. All rights reserved.
1.
Introduction
In clinical trials using count data of rare events, the Poisson distribution is a natural and commonly used model. Often times, interest is in the comparison of two Poisson rate parameters, as in Gu et al. [1], Ng and Tang [2], Thode [3], and Shiue and Bain [4]. In general, Bayesian methods are becoming increasingly more popular in clinical trials because of their flexibility and their ease of interpretation. Because of the high cost and time constraints for clinical trials, researchers often need to determine the smallest sample size n that provides accurate inferences for a parameter of interest. Although most experimenters have employed frequentist sample-size determination methods, the Bayesian paradigm offers a wider variety of sample-size determination methodologies and has been recommended over frequentist methods by Adcock [5]. The main limitation of the frequentist
∗
approach is that one must assume a “known” parameter value or vector. For example, one must assume that the variance or noncentrality parameters are known to perform a sample-size determination calculation for comparing two means. Obvious concerns are how one obtains this “known” parameter value and the consequences on the determined sample size if this “known” parameter is inserted without some incorporation of the uncertainty or variability. In contrast, Bayesian methods not only incorporate uncertainty about unknown parameter values, but also provide a variety of sample-size determination criteria. Two commonly implemented Bayesian sample-size methods are the average length criterion (ALC) method proposed by Joseph et al. [6] and the average power (AP) method proposed by Wang and Gelfand [7]. Recently, Zaslavsky [8] has used prior data to express uncertainty about parameters and an interval-based criterion to determine the sample size for the estimation of a Poisson rate. In this paper we extend Zaslavsky’s sample-size determina-
Corresponding author. Tel.: +1 254 315 5854. E-mail address: Austin [email protected] (D.M. Young). 0169-2607/$ – see front matter © 2010 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2010.10.010
272
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
tion method to the comparison of two Poisson rates, develop methods for sample-size determination via the AP method, and expand methods described by Stamey et al. [9] for the ALC method. Although a fully Bayesian approach requires the specification of a loss or utility function, the ALC method uses a sample-size criterion without specifying either, thus making the ALC method simpler to apply than other sample-size determination methods. Adcock [5] has noted that one may encounter considerable difficulty in formulating a loss or utility function that applies to all possible sample data. For the ALC method we fix the credible level, ˛, and determine the sample size n such that, on average, the length of the highest posterior density (HPD) interval is equal to a predetermined value l. The AP method uses the average posterior power to determine a suitable sample size n. This method also requires no specification of a loss or utility function. Although the AP method may be computationally intensive, it is very useful for the hierarchical models that dominate the current Bayesian landscape. Because of the potentially large influence the prior may have on the posterior distribution for small samples, choosing a prior is a crucial step in Bayesian modeling. However, the prior should be justified and, therefore, the influence it has on the posterior should be thought of as a valid source of uncertainty incorporated with the parameter. A desirable trait of a prior distribution is that it allows researchers to include available information, or lack thereof, with a degree of uncertainty. One common practice is to base the prior parameters choice on previous studies or trials. This prior selection approach is commonly referred to as objective Bayesian inference. In this paper we formulate an R function that uses prior information from previous clinical trials to define the parameters of the conjugate gamma prior through the implementation of marginal maximum likelihood estimation that uses the R function optim. Using these parameters, we also develop R functions that provide the required sample size for the ALC and AP methods for the ratio of two Poisson rates. The remainder of the paper is organized as follows. In Section 2 we provide background information on the Poisson distribution as well as review the marginal maximum likelihood estimation method for the gamma prior parameters. In Section 3 we propose and derive sample-size determination methods for the ratio of two independent Poisson rate parameters. We continue in Section 4 with two real-life examples for sample-size determination using the conjugate gamma prior. Finally, in Section 5 we conclude with brief comments.
2. Marginal maximum likelihood estimation 2.1.
Xi ∼ Poisson(ni i ) so that
(1)
where xi = 0, 1, 2, . . . and i > 0. Although we can use one of many different priors for i such as the uniform or log-normal, we use the conjugate gamma prior. Therefore,
(i |ai , bi ) =
bai i (ai )
iai −1 e−bi i ,
(2)
where 0 < i < ∞ and ai , bi > 0 with mean ai /bi and variance ai /b2i . Thus, the posterior density of i is a +x
(ni + bi ) i i ai +xi −1 −i (ni +bi ) e . (ai + xi ) i
(i |xi , ai , bi , ni ) =
(3)
Also, the marginal density of Xi is
m(xi |ai , bi ) =
∞
P(Xi = xi |i )(i |ai , bi )di 0
=
(ai + xi ) (ai ) (xi + 1)
b ai n xi i i ni + bi
ni + bi
,
(4)
which is a negative binomial probability function. Given k independent previous Poisson trials with xij events each, where i = 1, 2 and j = 1, 2, . . ., k, the marginal likelihood is
M(ai , bi ) =
k
m(xij |ai , bi ).
(5)
j=1
The marginal maximum likelihood estimators (MMLEs) for the parameters ai and bi , i = 1, 2, that we use in (2) are
(aˆ i , bˆ i ) = arg max M(ai , bi )|xi1 , . . . , xik .
2.2.
(6)
Marginal maximum likelihood estimation
To determine (aˆ i , bˆ i ), i = 1, 2, we apply a simplified minimization method developed by Byrd et al. [10] in the R function optim (see Appendix (A)). Because the method developed by Byrd et al. is a minimization instead of maximization routine, we use the negative of the marginal as our criterion function. Thus, the minimized log -likelihood is
Theoretical background
Let Xi be the number of events in ni patient time, and let i be the event rate for population i, i = 1, 2. We assume the model
x
e−ni i (ni i ) i , xi !
P(Xi = xi |i ) =
L(ai , bi ) =
k j=1
− log m(xij |ai , bi ) ,
(7)
273
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
we find the smallest positive value n such that
where j = 1, 2, . . ., n. Hence, the marginal log -likelihood is
L(ai , bi )
k
=−
log
j=1 k
=−
(ai ) (xij + 1)
∞ ∞ 1
b ai n xij i i
(ai + xij )
ni + bi
ni + bi
x2 =0
log (ai + xij ) − log ( (ai )) − log (xij + 1)
j=1
(8) Therefore, the MMLEs are
(aˆ i , bˆ i ) = arg min L(ai , bi )|xi1 , . . . , xik ,
(9)
where L(ai , bi ) is defined in (8).
3.
Sample-size determination criteria
3.1.
Average length criterion approach
The ALC approach, a method proposed by Joseph et al. [6] and discussed by Adcock [5], proposes that for a fixed posterior credible interval coverage, 1 − ˛, we determine the smallest sample size n such that
l(x, n)m(x)dx ≤ l, where l(x, n) is the length of the 100(1 − ˛)% posterior credible interval for the data x and l is a set length. We now consider the case for comparing two independent Poisson rate parameters using the ALC sample-size determination method. That is, we assume Xi ∼ Poisson(ni i ), i = 1, 2. Although several criteria for comparing the two rate parameters are available, one common methodology is the ratio = 1 /2 . As derived by Price and Bonnet [11], the posterior distribution is k∼Fv1 ,v2 ,
(10)
(a2 + x2 )(b1 + n1 ) (a1 + x1 )(b2 + n2 )
with v1 = a1 + x1 and v2 = a2 + x2 . Therefore the ALC solution, as described by Stamey et al. [9], is a sum since both x1 and x2 are discrete. The sum is over all possible outcomes and averages the (1 − ˛)100% posterior intervals, weighted by the marginal distribution. Thus, we have
∞
Fv1 ,v2
− Fv1 ,v2
˛ 2
m(x1 |aˆ 1 , bˆ 1 )
m(x2 |aˆ 2 , bˆ 2 ) ≤ l,
3.2.
(11)
∞
l(x1 , x2 , n)m(x1 |aˆ 1 , bˆ 1 )m(x2 |aˆ 2 , bˆ 2 ).
x2 =0x1 =0
The length of the interval for a particular value of x1 and x2 is computed by finding the (1 − ˛/2) and the ˛/2 percentiles of the F distribution and multiplying by the constant 1/k. That is
Hypothesis testing
We are often interested in testing the difference of the placebo and treatment groups in clinical trials. Lindley [12] proposed a multiplicative model where denotes the rate for the placebo group and denotes the rate for the treatment group. Then, < 1, = 1, and > 1 correspond to an improvement, no change, or a deterioration, respectively, for the treatment group. In general, we are interested in a one-sided alternative hypothesis that represents an improvement or a deterioration for the treatment group. The likelihood is P(, |Xi = xi ) = e−(n1 +n2 ) (x1 +x2 ) x2 ,
(12)
where xi = 0, 1, 2 . . ., i = 1, 2, and > 0. In practice one is often required to have a balance between the placebo and treatment groups. Thus, we assume n1 = n2 = n. Given the likelihood (12) and the conjugate gamma prior (2), ∼ gamma(a1 , b1 ), the marginal posterior distribution for is
P(|Xi = xi )
=
∞
P(, |Xi = xi )(|a1 , b1 )d 0
∝
x2 x +x +a , (n + n + b1 ) 1 2 1
where xi = 0, 1, 2 . . ., i = 1, 2, and > 0. The marginal density, m(x1 | a1 , b1 ), for the placebo group has a closed form as defined in (4). However, the treatment group marginal, m(x2 | a2 , b2 ), does not have a closed form. Hence, we employ numerical integration to approximate the marginal posterior distribution for .
3.3.
where k=
k
˛ 1− 2
where m(xi |aˆ i , bˆ i ), i = 1, 2, refers to the negative binomial marginal distribution (4).
+ ai log (bi ) − ai log
(ni + bi ) + xij log (ni ) − xij log (ni + bi ) .
x1 =0
Average power method
One can find a general algorithm for computing the average power (AP) in Wang and Gelfand [7]. Here, we determine the smallest sample size n given the parameter of interest such that
E I [P( > 0 |x1 , . . . , xn ) > 1 − ˛]
≥ 1 − ˇ,
(13)
where I( · ) is one if H0 is rejected and zero otherwise. Note that (13) is the expected probability of rejecting the null hypothesis. For the calculation of the average power, we first assign a value for the parameter of interest, while assuming that the alternative hypothesis is true, which we denote by ∗ . We then generate values for our treatment and placebo group, where x1 ∼ Poisson(n) and x2 ∼ Poisson(n ∗ ). Letting x1 and x2 be our data, we use Monte Carlo simulation to approximate the pos-
274
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
terior probability P( > 0 | xi ), where i = 1, 2. Using the posterior probability, for those posteriors with P( > 0 | x1 , . . ., xn ) > 1 − ˛ we assign a value of one to an element of the vector, b, and a value of zero, otherwise. By averaging over the vector b, we have an estimation for the average power. This process can be repeated for several different sample-size values, and then we can use a parametric or spline fit to determine the sample size that satisfies (13). In order to perform a Monte Carlo simulation, we must place priors on all parameters. We place a conjugate gamma prior (2) on the rate parameter , ∼ gamma(a1 , b1 ). For the parameter we use a relatively noninformative gamma prior, ∼ gamma(0.01, 0.01), which has a mean of one and is relatively flat around = 1. Thus, the alternative hypothesis will not be preferred.
4. Example: thromboembolism events for the CarboMedics mitral valves 4.1.
Description of Data
Thromboembolism is the formation of a clot (thrombus) in a blood vessel that breaks loose and is carried by the blood stream to plug another vessel. We are specifically interested in thromboembolism in the mitral valve, which is a dual-flap valve that lies between the left atrium and the left ventricle in the heart. A normally-functioning mitral valve opens secondary to contraction from the chordae tendineae, allowing blood to flow into the left ventricle during diastole, and closes at the end of the atrial contraction to prevent blood from backflowing into the atria during left ventricle systole [13]. Data is given in Zaslavsky [8] on thromboembolism when we use CarboMedics mitral valves to replace dysfunctional mitral valves in Table 1. By optimizing the marginal log -likelihood (9) for the thromboembolism data given in Table 1, we determined that the MMLEs for the first four series are â1 = 0.68 and bˆ 1 = 4.22, and the MMLEs for the latter four series are â2 = 1.75 and bˆ 2 = 8.72. Hence, the priors are 1 ∼ gamma(0.68, 4.22) and 2 ∼ gamma(1.75, 8.72), respectively.
Table 1 – Thromboembolism events for the CarboMedics synthetic mitral valves. First four series Events (x) Times (n)
16 43.1
12 58
0 67.7
6 80
Latter four series 16 109.6
30 139.8
64 157.1
6 181.5
4.2. Example: average length criterion sample-size determination for thromboembolism data We desire to compare the two independent Poisson rate parameters. Thus, we wish to determine the smallest sample size n that yields a relatively high probability of finding a difference in the two rates. Using the average length criterion for the two sample case (11), we set the cutoff level to be l = 1. We have developed a function in R that calculates the smallest sample size n for the comparison of two Poisson rates with gamma priors. Using this function, we found that the smallest sample size is n = 2614 0.1 person-years or 261.4 person-years (see Fig. (1)). Code for the ALC method can be obtained from the authors.
4.3. Example: average power sample-size determination for thromboembolism data For the thromboembolism data, we use only the first four series to estimate the prior for the control group. We are interested in the smallest sample size n that allows us to compare old CarboMedics mitral valves with new mitral valves. If the new mitral valves are better than the old CarboMedics mitral valves, then the rate parameter should be less for the new valves, and, hence < 1. The corresponding null and alternative hypotheses would be H0 : = 1 vs.
HA : < 1.
Using WinBugs, we placed a noninformative prior on the parameter , ∼ gamma(0.01, 0.01), along with the conjugate gamma prior on the first four series, ∼ gamma(0.68, 4.22), and let ∗ = 0.8. Using these values, we simulated the average power of our data for various sample sizes n (see Appendix (B)). Because it is common practice to use a minimum power of 0.8,
Fig. 1 – Average length criterion sample-size determination for CarboMedics synthetic mitral valves.
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
275
Fig. 2 – Simulated average power sample-size determination for CarboMedics synthetic mitral valves with priors ∼ gamma(0.68, 4.22) and ∼ gamma(0.01, 0.01).
we determine the smallest sample size n that gives an average power of 0.8. A plot of the power versus sample size n is given in Fig. (2). From this, we can see that the smallest sample size with an average power of 0.8 is n ≈ 1882 0.1 person-years or 188.2 person-years.
5.
Comments
Zaslavsky [8] has developed methods for sample-size determination for a single Poisson rate, which we expanded to the comparison of two independent Poisson rates by developing methods for sample-size determination using the AP method proposed by Wang and Gelfand [7] and expanding methods described by Stamey et al. [9] for sample-size determination via the ALC method. Although many methods for sample-size determination exist, Bayesian methods provide a wider variety and incorporate uncertainty for the parameter of interest with the specification of the prior distribution. Two commonly implemented Bayesian methods are the average length criterion (ALC) and the average power method (AP). One can implement both methods without the specification of a loss or utility function, thus making them easier to apply than competing Bayesian sample-size determination methods. The AP method is generally computationally more rigorous, but it is very useful for the hierarchical models that dominate the current Bayesian literature. We have created functions in R that determine the parameters for the conjugate gamma prior (see Appendix (A)) and calculate and display the sample size for the ALC and AP methods (see Appendix (B)). We compared the different methods available in the optim function in R and found that the method developed by Byrd et al. [10] consistently performed better than the other methods offered. Another reason for choosing this method is that it allows one to specify limits for the parameters of interest. For our situation this is important because we require ai , bi > 0 for the conjugate gamma prior (2). We have also provided two example implementations using the thromboembolism data in Sections 4.2 and 4.3.
Although the ALC method provides an exact sample size, the AP method must be approximated because the marginal distribution of the treatment group has no closed form. A major advantage of the ALC and AP methods are that they provide fast and easy methods to determine an appropriate sample size for the comparison of two independent Poisson rates.
Acknowledgements The authors would like to thank the Associate Editor and referee for comments that improved the manuscript.
Appendix A. Marginal maximization likelihood estimator function MMLE=function(x,n){ # define the negative log likelihood loglike=function(parameters){-(sum( lgamma(x+parameters[1]))-length(x)* lgamma(parameters[1])-sum(lgamma(x+1)) +length(x)*parameters[1]*log(parameters[2])sum(parameters[1]*log(parameters[2]+n))+ sum(x*log(n))-sum(x*log(n+parameters[2])))} # load the library stats for optim function library(stats) # set the lower limit for the parameters of # interest, if the function returns the error # given below, adjust this value # lowerlim=0.000000000000000001 # optimize the negative log likelihood function # with the function optim; many methods are # available, we use a method developed by # Byrd et al. that determines the minimum opt=optim(c(100,100),loglike,lower=lowerlim, method=”L-BFGS-B”) # check that the lower limit is set small # enough if(opt$par[1]==lowerlim){return(”The lower
276
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
limit in the optim function is not small enough.”)} if(opt$par[2]==lowerlim){return(”The lower limit in the optim function is not small enough.”)} # return the optimized value opt$par}
Appendix B. Average power sample-size determination # function that calculates the power for given # alpha and beta values power=function(a,b,a.bugs,b.bugs,theta,n,m,n. iter,greater){ # load the package R2WinBUGS to implement MCMC # methods library(R2WinBUGS) POW=seq(0,0,length=m) Post.prob=vector() # specifies whether testing greater than or # less than in the alternative hypothesis, # if greater is 1 then we are interested in # testing greater than, and 0 otherwise if(greater==1){parameters=list(”post.prob. great”)} if(greater==0){parameters=list(”post.prob. less”)} # for loop that runs m iterations to calculate the average power # for given alpha and beta values for(i in 1:m){ lambda=rgamma(1,a,rate=b) # generate a random value for our placebo and # treatment groups, respectively y1=rpois(1,n*lambda) y2=rpois(1,n*lambda*theta) data=list(”y1”,”y2”,”n”,”a.bugs”,”b. bugs”) pow=bugs(data,inits=NULL,parameters,”powerbugs. txt”,n.chains=1,n.burnin=1000,n.iter=n.iter) Post.prob[i]=pow$summary[1,1] # assign a value of 1 if the posterior # probability is greater than 0.95 and a # value of 0 otherwise, vector b in text. if(Post.prob[i]>0.95){POW[i]=1}} V # the average power for given alpha and beta values Power=mean(POW) return(Power)} # function that calculates the smallest sample # size for given alpha or beta values POWpoint=function(a,b,a.bugs,b.bugs,theta,m,n. iter,greater,beta0){ n=c(1,2000) POW=c(0,0) # check to see that the power chosen is larger # than the power for n=1, if not displays error # message
POW[1]=power(a,b,a.bugs,b.bugs,theta,n[1],m,n. iter,greater) if(POW[1]>=(1-beta0)){return(“The power chosen is too small. The sample size n=1 has a power greater than one minus beta.”)} POW[2]=power(a,b,a.bugs,b.bugs,theta,n[2],m,n. iter,greater) # increases n if n=2000 is too small for the # specified power value while(POW[2]<(1-beta0)){ n[2]=2*n[2] POW[2]=power(a,b,a.bugs,b.bugs,theta,n[2],m,n. iter,greater)} nmax=n[2] # calls the function Powplot to create a plot # for the power given specified alpha and # beta values, it also returns coefficients # for sample-size calculation coeff=Powplot(a,b,a.bugs,b.bugs,theta,m,n.iter, greater,beta0,nmax) # uses coefficients from a polynomial fit to calculate the # smallest sample size for specified power while((n[2]-n[1])>1){ N=ceiling((n[2]+n[1])/2) POW2=coeff[1]+N*coeff[2]+N ∧2*coeff[3]+N ∧3*coeff[4] if(POW2>=(1-beta0)){n[2]=N}else{n[1]=N}} # returns the smallest sample size calculated return(n[2])} # function that creates a power plot # for specified alpha and beta values POWplot=function(a,b,a.bugs,b.bugs,theta,m,n. iter,greater,beta0,nmax){ n=seq(50,nmax,by=50) powplot=vector() for(i in 1:length(n)){ # calls power function to calculate the power # for specified alpha and beta values powplot[i]=power(a,b,a.bugs,b.bugs,theta,N,m,n. iter,greater)} # fits a third degree polynomial to the power # plot fit=lm(powplot∼n+I(n ∧2)+I(n ∧3)) plot(n,powplot,xlab=”Sample Size n”, ylab=”Power”,main=”Average power sample-size determination for CarboMedics synthetic mitral valves”) lines(n,fit$coeff[1]+n*fit$coeff[2]+n ∧2*fit$coeff[3]+n ∧3*fit$coeff[4]) abline(h=1-beta0) # returns the coefficients from the third # degree polynomial fit return(fit$coeff)} # winbugs model for sample-size calculation model{ # specify the poisson distribution for the # placebo and treatment group, # respectively
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 4 ( 2 0 1 1 ) 271–277
y1 ∼dpois(lambda1) y2 ∼dpois(lambda2) # specify the noninformative gamma prior theta ∼dgamma(0.01,0.01) # specify the conjugate gamma prior lambda ∼dgamma(a.bugs,b.bugs) # parameters for the placebo and treatment # group, respectively lambda1<-n*lambda lambda2<-n*lambda*theta post.prob.less<-1-step(theta-1) post.prob.great<-1-step(1-theta)} Similar code for the ALC method can be obtained from authors.
references
[1] K. Gu, H.K.T. Ng, M.L. Tang, W.R. Schucany, Testing the ratio of two Poisson rates, Biometrica 50 (2) (2008) 283–298. [2] H.K.T. Ng, M.-L. Tang, Testing the equality of two Poisson means using the rate ratio, Statistics in Medicine 24 (2005) 955–965. [3] H. Thode, Power and sample size requirements for tests of differences between two Poisson rates, The Statistician 46 (1997) 227–230. [4] W.K. Shiue, L.J. Bain, Experiment size and power comparisons for two-sample Poisson tests, Applied Statistics 31 (2) (1982) 130–134.
277
[5] C. Adcock, Sample size determination: a review, The Statistician 46 (1997) 261–283. [6] L. Joseph, D.B. Wolfson, R. Berger, Some comments on Bayesian sample size, The Statistician 44 (2) (1995) 167–171. [7] F. Wang, A.E. Gelfand, A simulation-based approach to Bayesian sample size determination for performance under a given model and for separating models, Statistical Science 17 (2) (2002) 193–208. [8] B.G. Zaslavsky, Empirical Bayes models of Poisson clinical trials and sample size determination, Pharmaceutical Statistics (2009). [9] J.D. Stamey, D.M. Young, T.L. Bratcher, Bayesian sample-size determination for one and two Poisson rate parameters with applications to quality control, Journal of Applied Statistics 33 (6) (2006) 583–594. [10] R. H. Byrd, C. Zhu, P. Lu, J. Nocedal, L-BFGS-B—Fortran Subroutines for Large-Scale Bound Constrained Optimization, Technical Report, ACM Trans. Math. Software, (1994). [11] R. Price, D. Bonnett, Estimating the ratio of two Poisson rates, Computational Statistics and Data Analysis 34 (2000) 345–356. [12] D.V. Lindley, A holistic analysis of Poisson data with application to a trial of selenium and cancer deaths, The Indian Journal of Statistics 64 (3) (2002) 722–739. [13] W. Carey, Current Clinical Medicine 2009: online + print, Saunders/Elsevier, (2009).