Authors: Nigel Stallard, Frank Miller, Simon Day, Siew Wan Hee, Jason Madan, Sarah Zohar, Martin Posch
Categories: Special Topic: ISCB2015, Bayesian, Clinical trial design, Decision theory, Exponential family form, Optimal sample size, Research Paper
Source: Biometrical Journal. Biometrische Zeitschrift
The problem of choosing a sample size for a clinical trial is a very common one. In some settings, such as rare diseases or other small populations, the large sample sizes usually associated with the standard frequentist approach may be infeasible, suggesting that the sample size chosen should reflect the size of the population under consideration. Incorporation of the population size is possible in a decision‐theoretic approach either explicitly by assuming that the population size is fixed and known, or implicitly through geometric discounting of the gain from future patients reflecting the expected population size. This paper develops such approaches. Building on previous work, an asymptotic expression is derived for the sample size for single and two‐arm clinical trials in the general case of a clinical trial with a primary endpoint with a distribution of one parameter exponential family form that optimizes a utility function that quantifies the cost and gain per patient as a continuous function of this parameter. It is shown that as the size of the population, N, or expected size, N∗ in the case of geometric discounting, becomes large, the optimal trial size is O(N1/2) or O(N∗1/2). The sample size obtained from the asymptotic expression is also compared with the exact optimal sample size in examples with responses with Bernoulli and Poisson distributions, showing that the asymptotic approximations can also be reasonable in relatively small sample sizes.
Keywords: Bayesian, Clinical trial design, Decision theory, Exponential family form, Optimal sample size
The problem of determining the sample size for a clinical trial is a very common one. For large‐scale definitive phase III clinical trials, a frequentist approach is usually adopted, with the sample size chosen so as to control the type I error rate at a specified level, α, and to give specified power 1−β, to detect some appropriately chosen size of treatment effect (see, e.g. Pocock, 1983, for details). Choice of α=0.05 and β=0.1 or 0.2 is typical.
The sample sizes obtained using the frequentist approach do not always seem appropriate. In particular they do not reflect the size of the population to which the results of the trial apply. The population size is relevant, however, when considering decisions made on the basis of trial results. This is particularly true for clinical trials conducted in rare diseases or other small populations, when the population size means that a large trial would be infeasible, or even impossible.
One way in which the size of the population can influence the sample size is to use a decision theoretic approach in which the benefits to future patients in the population, sometimes called the “patient horizon”, are explicitly considered so that future benefit depends on the size of this population. Such an approach has been proposed and discussed by numerous authors over the last 50 years (see, e.g. Anscombe, 1963; Colton, 1963; Sylvester, 1988; Berry et al., 1994; Cheng et al., 2003; Kikuchi and Gittins, 2009 and reviews by Pezeshk et al., 2013; Hee et al., 2016). Although this approach has very rarely been implemented in practice, it can nevertheless provide important insight into an appropriate choice of sample size for a clinical trial.
The role of the population size in determination of the optimal sample size for a clinical trial has been considered by Cheng et al. (2003). They considered single‐arm and two‐arm clinical trials with a primary endpoint following a Bernoulli distribution indicating either success or failure of treatment for each patient. Adopting a decision‐theoretic approach, they obtained designs that maximize the total number of successes. Denoting the population size by N, they show that asymptotically as N→∞, the optimal sample size for a clinical trial is O(N1/2) and give an expression for the asymptotically optimal sample size that depends on the prior distributions for the unknown probability of success for each trial arm.
In this paper, we extend the work of Cheng et al. (2003), both to more general distributional forms for the primary endpoint and to situations in which the aim is to maximize some general utility expressed as a function of a parameter of these distributions. We show that the result that the optimal sample size is O(N1/2) applies for any continuous utility function and for responses with a distribution of any one‐parameter exponential family form assuming a conjugate prior distribution. We also consider the case where no finite patient horizon is assumed, but gains from future patients are geometrically discounted, that is the gain from the jth patient is multiplied by λj−1 for some discounting parameter, λ<1 (Berry and Fristedt, 1985). As considered in the discussion section, the size of λ can also reflect the size of the patient population in this setting, via an effective number of patients, (1−λ)−1, which will be denoted N∗. We show that in this case as λ→1 so that N∗ is large, the optimal sample size is O(N∗1/2). We also investigate through exact calculation the small‐sample accuracy of the the large sample approximations. Although the results obtained depend on asymptotics, we show that, depending on the exact form of the utility function chosen, these may be reasonable even for extremely rare diseases, for example for patient populations of 1000 or less when the optimal sample size can be less than 50.
Suppose that a clinical trial is to be conducted to choose between two treatments with n
1 and n
2 patients receiving treatment 1 and treatment 2, respectively, where treatment 2 may be the current standard treatment included as a control. Note that taking n2=0 corresponds to a single‐arm trial, though the decision to be taken at the end of such a trial remains comparative, with a choice being made regarding treatment of future patients.
It is assumed that the gain associated with treatment of patients in the trial or patients outside the trial who receive each treatment can be specified as a function of a parameter of the distribution for the response for patients receiving that treatment. It is noted that here “gain” is to be interpreted widely to include any kind of costs, losses, gains, or benefits associated with treatment. Following the trial, the most preferable treatment, that is the treatment for which the posterior expected gain given the observed data, is highest, will be selected. The remaining patients will then receive this treatment. We wish to determine the optimal values, n1∗ and n2∗, of the sample sizes n
1 and n
2 and to determine how n1∗ and n2∗ depend on the population size.
We will assume that responses for patients follow some distribution of natural one‐parameter exponential family form. In detail, let Yij denote the response for patient j receiving treatment i and assume Yi1,…,Yini are i.i.d. with density fi(y∣ψi)=ai(y)exp(yψi−bi(ψi)) for some ai(y) and bi(ψi). Typically responses from patients in the two treatment groups will follow distributions of the same form, that is functions ai and bi will not depend on i, with ψ1 and ψ2 differing, though this is not assumed. We will assume that ψ1 and ψ2 are taken to have independent prior distributions of conjugate form, that is with ψi having density πi(ψi∣n0i,y0i)=ci(n0i,y0i)exp(n0iy0iψi−n0ibi(ψi)) for some y0i and n0i and normalising constant ci(n0i,y0i),i=1,2. The values y0i and n0i can be interpreted respectively as the prior mean of ξi=bi′(ψi)=E(Yij) and the number of observations to which the prior information is equivalent, so that following a trial with ni patients receiving treatment and observation of Y=(Y11,…,Y1n1,Y21,…,Y2n2), the posterior mean for ξi given Y is equal to
with Y¯i=∑j=1niYij/ni (see Bernado and Smith, 2000).
Suppose that the expected gain from a patient receiving treatment i in the trial is hi(ξi), and that the expected gain to a future patient receiving treatment i is gi(ξi) where hi and gi,i=1,2 are such that the expected values E0(hi(ξi)), E0(gi(ξi)) and E0(maxi=1,2(gi(ξi))) where E
0 denotes the expected value taken over the prior distribution of ξ, exist and hi and gi are assumed to be differentiable with gi strictly increasing and with finite derivative. Assume further that
This ensures that the gain from treating patients in the trial cannot exceed that from treating them outside the trial. This is considered further in the discussion section below.
We will consider two cases. In the first, the population is considered to be finite with known size, N. The number of patients treated following the trial is thus N−n1−n2. In the second case, no finite population size is assumed, but the gain from future patients is geometrically discounted, so that the gain from patient j if they receive treatment i is λj−1hi(ξi) if they are included in the trial and λj−1gi(ξi) if they are treated following the trial, for some λ<1.
The geometric discounting of gains from future patients can be interpreted in a number of ways. One interpretation is that gains further in the future are reduced to reflect either opportunity loss or loss of financial interest on an investment (see, e.g. Fergusson, 2008). With this interpretation it might be appropriate to take λ<1 and N finite. An alternative interpretation is to imagine that gain from each future patient is of constant value, as is assumed in the finite horizon case, but that the size of the population, N, rather than being fixed in advance, is random, following a geometric distribution. This could be the case if, for example, the population is limited by some new treatment becoming available at which point the trial, or use of the recommended treatment following the trial, will be terminated, with the probability of this event constant over time (see Berry and Fristedt, 1985). In this interpretation N∗=(1−λ)−1 is the expected population size prior to this new treatment becoming available. The size of λ and the resulting N∗ can thus reflect the population size, and in some ways a smaller value of λ, corresponding to fewer patients being available prior to the new treatment becoming available, more reasonably models a small population than assuming the number of patients has some fixed and known value, N.
Consider first the setting of a finite patient horizon of size N. Following observation of data Y= (Y11,…,Y1n1,Y21,…,Y2n2), the total expected gain if treatment i is recommended for all further N−n1−n2 patients is n1Eξ∣Y(h1(ξ1)∣Y)+n2(Eξ∣Y(h2(ξ2)∣Y)+(N−n1−n2)Eξ∣Y(gi(ξi)∣Y) where Eξ∣Y(.∣Y) denotes the expected value taken over the posterior distribution of ξ given Y.
The optimal action at the end of the trial is thus to select the treatment with the largest value of Eξ∣Y(gi(ξi)∣Y) and the expected gain assuming this action is taken is equal to n1Eξ∣Y(h1(ξ1)∣Y)+n2Eξ∣Y(h2(ξ2)∣Y)+(N−n1−n2)maxi=1,2Eξ∣Y(gi(ξi)∣Y).
Since, prior to the commencement of the trial, Y is unknown, the expected gain from the trial is equal to E0(G) where the expectation is taken over the prior distribution for ξ1 and ξ2 and G is the function of ξ1, ξ2, N,n1, and n
2 given by
where EY denotes the expectation taken over Y for a given value of ξ1 and ξ2 so that expectations are taken first over the posterior distributions of ξi given Y and then over Y given ξ1 and ξ2.
Since E0(EY(Eξ∣Y(hi(ξi)∣Y))) is equal to the prior expectation E(hi(ξi)) for any function hi for which the expectations exist, we get
We wish to find the optimal values of n
1 and n
2, that is the values for which E0(G) is maximized. For small N, when the optimal sample sizes, n1∗ and n2∗ will also be small, it may be feasible to evaluate the prior expected gain given by (4) directly, taking the expectation over the prior predictive distribution for Y and to find n1∗ and n2∗ by a numerical search. For larger N, such an approach may be infeasible. In this case asymptotic expressions for the optimal sample sizes, n1∗ and n2∗, as the population size, N, becomes large are more useful.
For finite n
01 and n
02, the expectation E0(EY(maxi=1,2Eξ∣Y(gi(ξi)∣Y))) is increasing in n
1 and n
2. Thus for N>n1−n2, (N−n1−n2)E0(EY(maxi=1,2Eξ∣Y(gi(ξi)∣Y))) and hence E0(G) is also increasing in n
1 and n
2. Thus as N→∞ the optimal trial design has both (n01+n1)→∞ and (n02+n2)→∞. The case in which both n
01 and n
02 are both infinite corresponds to both ξ1 and ξ2 being known a priori. We will therefore consider the case in which, without loss of generality, n
01 is finite, and derive the optimal value n1∗.
Note that as N→∞, so that the optimal n01+n1 and n02+n2 also approach infinity, E0(EY(maxi=1,2Eξ∣Y(gi(ξi)∣Y)))→maxi=1,2E0(gi(ξi)) so that the optimal sample sizes are such that ni/N→0,i=1,2.
By the central limit theorem we have Y¯i→dN(ξi,vi(ξi)/ni) where vi(ξi) denotes the variance of Yij. Thus from (1), applying the delta method since gi is assumed to be differentiable and strictly increasing so that the derivative, gi′(ξi), is non‐zero, we get
Using an expression for the expected value of the maximum of two normally distributed random variables given by Clark (1961), we have
where ξ∼i=gi((n0iy0i+niξi)/(n0i+ni)), a2=∑i=12ni2(gi′(ξi))2vi(ξi)/ni(n0i+ni)2, δ=(ξ∼1−ξ∼2)/a and ϕ and Φ denote standard normal density and distribution functions. Thus
In order to find n1∗, we obtain the derivative ∂E0(G)/∂n1=E0(∂G/∂n1). From (3),
We will thus find a large‐sample approximation for this derivative.
Equation (6) gives an approximation for EY(maxi=1,2Eξ∣Y(gi(ξi)∣Y)). The derivative of the right hand side of (6) is equal to
(see Web Appendix A for details), so that ∂G/∂n1 can be approximated by
The limit of this as N, n
1 and n02+n2 all approach infinity with n1/N and n2/N approaching 0 (again, see Web Appendix A for details), has expected value
where π(ξ1,ξ2) denotes the joint prior distribution of ξ1 and ξ2.
Setting this to zero and solving for n
1, the maximum expected gain is found to be at
Note that the fact that g
1 is increasing and the requirement (2) ensures that both numerator and denominator are positive so that the square root exists. It is also interesting to note that the asymptotic optimal sample size for arm 1, n
1, does not depend on n
2 since we have assumed either that n
02 is infinite or that n
2 also approaches infinity.
When n
02 is finite, by symmetry, the optimal value of n
2 is given by
When n02=∞, the prior distribution has mass at ξ2=y02 only, and so may be written as a univariate density π(ξ1), so that the optimal value of n
1 becomes
As n02→∞, ∂EY(maxi=1,2Eξ∣Y(gi(ξi)∣Y))/∂n2→0. Thus from an expression similar to (8) giving the derivative with respect to n
2 and (2), the derivative is negative and n2∗=0.
When g
1 and g
2 are identical, (10) becomes
and when gi(ξi)=hi(ξi)=ξi,i=1,2, this becomes
showing that this is a generalization of the expression obtained by Cheng et al. (2003) for the case in which Yij has a Bernoulli distribution with parameter ξi and vi(ξi)=ξi(1−ξi).
Consider next the second setting introduced above; that of an infinite population with geometric discounting.
In a two‐arm trial it is assumed that n
1 and n
2 are sufficiently large and randomisation to treatments 1 and 2 sufficiently balanced that the gain to patients receiving treatment i in the trial can be taken to be ni∑j=1n1+n2λj−1/(n1+n2)hi(ξi). The total expected gain if treatment i is recommended for all further N−n1−n2 patients is then
The optimal action at the end of the trial is thus again to treat all future patients with the treatment with the largest value of Eξ∣Y(gi(ξi)∣Y) and the expected gain assuming this action is taken is equal to
The prior predicted expected utility for a trial with ni patients receiving treatment i is thus
Optimal sample sizes, n1∗ and n2∗, can be found directly using this expression and a numerical search in cases when this is computationally feasible.
As above, it is of interest to seek asymptotic approximations to n1∗ and n2∗, in this case as the geometric discounting parameter, λ, approaches 1 from below. We will again assume that n
01 is finite and obtain first an approximation for n1∗.
The derivatives of ∑j=1n1+n2λj−1=(1−λn1+n2)/(1−λ) and ∑j=n1+n2+1∞λj−1=λn1+n2/(1−λ) with respect to n
1 are respectively −λn1+n2logλ/(1−λ) and λn1+n2+1logλ/(1−λ) which, by L'Hospital's rule, tend to 1 and −1, respectively as λ→1, so that, since the limit as λ→1 of ∑j=1n1+n2λj−1 is n1+n2, the derivative of the term ni∑j=1n1+n2λj−1/(n1+n2) with respect to n
1 tends to 1 if i=1 and 0 if i=2.
As λ→1, the derivative of (13) with respect to n
1 thus approaches
The argument above gives an approximation to this derivative of
which, since ∑j=n1+n2+1∞λj−1=∑j=1∞λj−1−∑j=1n1+n2λj−1, with the former term, which is equal to (1−λ)−1, dominating, gives the approximation
Writing N∗=∑j=1∞λj−1=(1−λ)−1, (14) can be written as
directly analogous to (10) with N∗ replacing N.
For n
02 finite, n2∗ is again given by symmetry by an expression analogous to (11). For n02=∞, n2∗=0 and n1∗ is given by an expression analogous to (12) with N replaced by N∗.
We consider first the case of Bernoulli data. In this case the distribution of Yij, the responses in treatment group i, can be parameterised with ξi equal to the probability of treatment success. We will take ξ1 to have a conjugate beta prior with parameters a
1 and b
1 so that y01=a1/(a1+b1) and n01=a1+b1, and assume ξ2 is known with value y
02, that is with n02=∞, so that n2∗=0 and a single‐arm trial is optimal.
Given observation of data Y1j=y1j,j=1,…,n1, ξ1 has a Beta(a1+∑j=1n1y1j,b1+n1−∑j=1n1y1j) posterior distribution, and the prior predictive distribution of ∑j=1n1Y1j is Betabinomial(n
1, a
1, b
1).
We will assume that the gain from patients receiving treatment i will be determined by whether or not the treatment is successful, so that gi(ξi)=hi(ξi)=ξi,i=1,2. From (4) with n2=0, the prior expected gain is thus n1E0(ξ1)+(N−n1)E0(EY(maxi=1,2Eξ∣Y(ξi)∣Y)).
Following one of the examples considered by Cheng et al. (2003), we take a1=1 and b1=1, the known value of ξ2 to be 0.5 and N=100. Figure 1 shows the prior expected gain for a range of values of n
1, here plotted on a logarithmic scale. As given by Cheng et al. (2003), the optimal value of n
1 is equal to 9, which is marked with a plus sign. The approximation to the prior expected gain given by (7) is also shown on the figure as a dashed line, showing that in this case the approximate and true values are close even for small n
1. The approximately optimal value of n
1 given by (12) is 10. This is shown by the circle on the figure, showing that this is close to the true optimum. In this case, this is close to the value of n
1 maximizing the approximation given by (7), though in general the additional approximation leading to (12) means that this need not be the case.
Figure 1 Prior expected gain (solid line) and approximate prior expected gain (dashed line) for a range of n
1for the first single‐arm Bernoulli example with a1=b1=1, ξ2=0.5, and N=100. The optimal and approximately optimal values of n1are marked by + and ○, respectively
Figure 2 gives values of n1∗ from (12) along with exact optima for a range of values of N, with both n1∗ and N plotted on logarithmic scales so that the square root relationship between n1∗ and N given by (12) corresponds to a straight line with slope 1/2. The approximation given by (12) is close to the true value, and approaches it as N increases as would be expected.
Figure 2 Optimal (solid lines with points) and approximately optimal (dashed lines) values for n
1for a range of N values for the first single‐arm Bernoulli example with a1=b1=1 and ξ2=0.5.
The effect of varying the prior distribution for ξ1 was also investigated. Web Figure 1 in Web Appendix B shows the optimal sample size n1∗ with N=100 and ξ2=0.5 for a range of n01=a1+b1 and y01=a1/(a1+b1) values. When the prior mean value of ξ1, that is y
01, is equal or greater than the fixed value of ξ2, the optimal sample size increases with prior weight n
01. This is reasonable since as n
01 increases we are increasingly confident that patients are not harmed by being in the trial and a larger sample size gives more information for the final decision. When y
01 is less than ξ2, the optimal sample size increases with n
01 for small n
01, so that more information can be collected for the final decision and decreases for larger n
01 when there is strong prior belief that treatment 2 is superior to treatment 1. It is interesting to note that when the prior weight, n
01, is kept fixed, the optimal sample size increases as the prior mean y
01 of treatment 1 increases. This is in contrast to frequentist sample size that would decrease as the assumed success rate for treatment 1 is increased (and is larger than the success rate of treatment 2).
We next consider an example based on Stallard (1998), who also considered a single‐arm phase II trial with a Bernoulli outcome. In this case the gain function was chosen to reflect the financial costs and rewards associated with the conduct of the trial assuming that, if successful, it would be followed by a further trial with a frequentist design and analysis. Assuming the probability of success for treatment 2, here taken to be the current standard treatment, to be known, costs and rewards were taken relative to continuing to give all patients the current treatment. Thus if this treatment is recommended, the gain to patients outside the trial is taken to be zero, that is g2=0. The gain per patient outside the trial if the experimental treatment is recommended was taken to be of the form g1(ξ1)=l(1−Φ(zα/2−(zα/2+zβ)θ/θ1))−m where zα=Φ−1(1−α), θ=logit(ξ1)−logit(ξ0), θ1=logit(ξ0+δ0)−logit(ξ0), and logit(ξ)=log(ξ/(1−ξ)) for some l, α, β, and δ0. This form reflects a fixed cost of m per patient with a gain of l per patient if the treatment is shown to be effective in the subsequent trial where that trial has frequentist (one‐sided) type I error rate of α/2 and power 1−β to detect a log‐odds ratio of θ1. Stallard assumed linear discounting for j less than some constant, n
0, with geometric discounting for j>n0 whereas we will assume geometric discounting for all j. Patients in the trial were taken to have constant (discounted) cost, that is h1(ξ1)=−k, for some k. Since we know n2=0, it is not necessary to specify h
2. The parameter ξ1 was again taken to have a Beta(a
1, b
1) prior distribution. As the gain function g2(ξ2) does not depend on ξ2, it is not necessary to specify a value for the point‐prior for this parameter in this case.
Since g2(ξ2)=0 for all ξ2, the expression (12) becomes
where π(ξ1) is the (univariate) prior density for ξ1. The value of E0(max{g1(ξ1),0}) can be evaluated using numerical integration and, as h1=−k, we have E0(h1(ξ1))=−k.
Although the form of utility function proposed by Stallard (1998) was not exactly that proposed here, based on values given, we took λ=exp(−0.00173)=0.99827, l=12.79, m=0.346, α=0.05, β=0.1, ξ0=0.2, δ0=0.15, a1=0.845 and b1=9.155. Figure 3 shows the prior distribution for ξ1 along with the form of g1(ξ1) in this case. Figure 4 shows the prior expected gain calculated exactly using the betabinomal prior predictive distribution for ∑j=1n1Y1j for a range of values of n
1 values, plotted on a logarithmic scale, along with the approximation given by (7). It can be seen that in this case the approximation (7) to the expected gain is rather poor, particularly for smaller n
1. The value of n1∗ obtained using (15) in this case is 102. This value and the associated prior expected gain is marked on the figure by a circle. Note again that n1∗ does not maximise the approximate gain given by (7) that is shown by the dashed line on the plot. In this case n
1 the value of n1∗ is quite far from the value of n
1 maximizing the approximate gain, though it is closer to the true optimal value of 95, which marked on the figure by a plus sign.
Figure 3 Prior distribution for ξ
1(upper panel) and gain function g1(ξ1) giving gain from treating each future patient with treatment 1 (lower panel). The gain function g2(ξ2)=0 is shown as a dashed line on the right hand panel for comparison (see text for details).
Figure 4 Prior expected gain (solid line) and approximate prior expected gain (dashed line) for a range of n
1for the second single‐arm Bernoulli example with N=5000 (see text for details of gain function and prior distribution parameter values). The optimal and approximately optimal values of n1are marked by + and ○ respectively
Figure 5 gives values of n1∗ from (15) along with exact optima for a range of values of N, again with both plotted on a logarithmic scale. The approximation given by (15) is again close to the true value, and approaches it as N increases as would be expected.
Figure 5 Optimal (solid lines with points) and approximately optimal (dashed lines) values for n
1for a range of N values for the second single‐arm Bernoulli example (see text for details of gain function and prior distribution parameter values).
The effect of varying the prior distribution for ξ1 was again investigated, and again illustrated in Web Appendix B. Web Figure 2 shows the optimal sample size n1∗ with N=5000 for a range of n01=a1+b1 and y01=a1/(a1+b1) values. In this case as y
01 increases from 0.04 to 0.0845 the optimal sample size is increased in a similar way to that noted for the example above. In this case as y
01 increases further, however, the optimal sample size is reduced. Here, since hi=−k, there is a cost associated with experimentation so that when there is strong prior belief that treatment 1 is superior to treatment 2, rather than giving many patients treatment 1 in a trial, a smaller trial is optimal.
At the end of the trial, treatment 1 will be recommended if Eξ∣Y(g1(ξ1)∣Y)>0. For large n
1, that is approximately if Eξ∣Y(ξ1∣Y)>g1−1(0), which, using expression (1), is true when n1Y¯1>g1−1(0)(n01+n1)−n01y01. Considering recommendation of treatment 1 to correspond to rejection of the null hypothesis that ξ1=g1−1(0)=0.201, frequentist error rates attained for the optimal designs obtained can be derived.
With a prior distribution with a1=0.845 and b1=9.155, taking n1=102 gives a type I error rate of 0.39. The form of g1(ξ1) shown in Fig. 3 suggests that we might require a test with high power when ξ=0.35, since this value of ξ1 is associated with a high gain value. The power of the optimal test in this case is 0.999.
Type I and type II error rate values for the optimal designs and prior distributions considered in Web Fig. 2 are shown in Web Fig. 3. As the prior weight becomes small, the optimal decision at the end of the trial is to select treatment 1 whenever it has observed mean exceeding g1−1(0), so the type I error rate approaches 0.5. As the prior weight increases, since in this case the optimal value of n
1 is relatively small, prior information comes to dominate the final decision and the type I error rate approaches zero or one, with the type II error approaching one or zero, depending on whether the prior mean is less than or greater than g1−1(0).
The third example is based on an example given by Berry et al. (1994), who describe a trial of an HIB vaccine in Navajo children aged 2–18 months. The number of HIB cases is assumed to follow a Poisson distribution. Rather than expressing n
1, n
2, and N in terms of child‐months, we will assume that all children are followed up for the entire 16‐month period, and refer to the number of children in the trial and population. The observed number of cases per child j in group i, will be denoted by Yij,j=1,…,ni,i=1,2. The distribution of Yij can be parameterised such that ξ1 and ξ2 are the expected numbers of cases per child for treatments 1 (the new vaccine) and 2 (the placebo), respectively. Thus p(yij∣ξi)=ξiyijexp(−ξi)/yij!,j=1,…,ni,i=1,2, so that yij has mean ξi, with ξi following independent prior gamma(αi, βi) distributions, that is with density π(ξ1,ξ2)=∏i=12ξiαi−1exp(−ξiβi)βiαi/Γ(αi) for some αi,βi,i=1,2. Note that ξi has prior mean αi/βi and prior variance αi/βi2. The posterior distribution of ξi given Y is a gamma (αi+∑j=1niyij, βi+ni) distribution and the prior predictive distribution of ∑j=1niyij is NegBin (αi,(1+ni/βi)−1).
Berry et al. (1994) include in their gain function a term that depends on the observed data that reflects the probability of obtaining regulatory approval for the vaccine. Here, we assume the gain from a child receiving treatment i depends on ξi alone and, since ξi gives the rate of HIB cases, which we would like to minimize, we take gain functions hi(ξi)=gi(ξi)=−ξi, i=1,2. The case of gain functions that depend on the observed data is considered briefly in the discussion section below.
The optimal values may be approximated using expressions (10) and (11). In this case vi(ξi)=ξi, gi′(ξi)=−1, and g2−1(g1(ξ))=g1−1(g2(ξ))=ξ so that, for example, n1∗ is
the integral in the numerator being equal to
Following Berry et al. (1994), we take (α1,β1)=(1,200) and (α2,β2)=(5,667), the latter corresponding to the placebo (note that the βi values given by Berry et al. are 16 times those used here as they take ξi to give the rate of cases per child‐month). Berry et al. report that approximately 5400 Navajo are born each year so that minimization of HIB cases over a 20‐year period would correspond to N=108,000. Figure 6 shows a contour plot giving the prior expected gain for this N for a range of n
1 and n
2 values (plotted on logarithmic scales) together with the approximation given by (7) (dashed lines). It can be seen that even for small sample sizes, (7) gives a close approximation to the true prior expected gain. The optimal design has n1=3162 and n2=1585, and is marked by the plus sign. The approximately optimal design given by (10) and (11) has n1∗=3524 and n2∗=2089, and is marked by a circle. The prior expected gain, in this case corresponding to minus one times the prior expected number of HIB cases in the population over the 20 year period, is −416.9 using the optimal design and −417.4 using the approximately optimal design.
Figure 6 Contour plot of prior expected gain (solid lines) and approximate prior expected gain (dashed lines) for N=108,000 for a range of n
1and n2values assuming gamma (1, 200) and gamma (5, 667) prior distributions. The optimal and approximately optimal values of n1and n2are marked by + and ○, respectively.
Figure 7 shows the values of n1∗ and n2∗ along with the approximations from (10) and (11) (dashed lines) for a range of N values, again with both plotted on logarithmic scales. It can again be seen how the approximations become increasingly accurate as N increases.
Figure 7 Optimal (solid lines with points) and approximately optimal (dashed lines) values for n
1(upper lines) and n2(lower lines) for a range of N values assuming gamma (1, 200) and gamma (5, 667) prior distributions.
The effect of varying the prior distributions for ξ1 and ξ2 was again investigated. Web Fig. 4 in Web Appendix B shows the optimal sample sizes, n1∗ and n2∗, for a range of prior means and prior weights, here equal for the two priors, when N=108,000. When the prior means are equal, the optimal sample size increases with prior weight. For unequal prior means, more patients are assigned to the arm considered a priori to be superior, in this case corresponding to a lower prior mean since hi(ξi)=gi(ξi)=−ξi, with the number assigned to the inferior arm increasing with prior weight when this is small and decreasing for larger prior weight values. Web Fig. 5 shows the effect of changing the prior weight for ξ1 alone when the prior mean is equal to, greater than or less than that for ξ2. In this case increasing the prior weight leads to an increase in the optimal sample size for both arms, with the arm with the lower prior mean having a smaller optimal sample size.
The examples above compare the large sample approximation for the optimal sample size for arm i, ni∗, given by expression (10), with that obtained by exact numerical optimization in two examples, enabling assessment of the approximation for smaller values of N. The derivation of (10) relies on large sample approximations in two ways; first the distribution of the posterior expected utility E(gi(ξi)∣Y¯) is approximated by its asymptotic normal form in (5) using the central limit theorem and the delta method, and second, the derivative of the expected gain, given by (8) is approximated by (9). The first approximation is exact when Yij are normally distributed and gi(ξi) are linear. The first and second examples suggest that both approximations are sufficiently accurate for Bernoulli data even for quite small N when gi(ξi) are linear, but less so for nonlinear gi(ξi), when the first approximation may be poor, as noted by Bernado and Smith (2000). For nonnormal data the accuracy of the first approximation also improves as the prior weight n0i increases, though, as illustrated in the third example above, this also leads to smaller ni∗, so that the overall accuracy of the asymptotic approximation to the optimal sample size might be poorer.
The work reported above leads to expressions for the optimal sample size in a clinical trial to compare two treatments or to compare a single experimental treatment with a historical control the properties of which are assumed known. The observed data for patients receiving treatment i are assumed to follow a distribution of one parameter exponential family form, with mean ξi assumed to have a conjugate prior distribution. Optimization is based on consideration of the costs and benefits both from patients in the trial, given by some utility function, hi(ξi) for patients receiving treatment i, and from subsequent patients who will receive treatment i based on the results of the trial, given by some function gi(ξi) if they receive treatment i. Although the expressions obtained could be directly used to design a trial, it is also of more general interest to see how the optimal sample size depends on the size of the population under investigation. We have shown that if the population is assumed to be of some known size, N, for any hi and gi satisfying sufficient regularity conditions for expectations to exist, differentiable with gi strictly increasing, and satisfying the condition given by (2); that is E0(hi(ξi))≤E0(maxi=1,2(gi(ξi))), the optimal sample size is O(N1/2) as N→∞. If it is assumed that there is an infinite population with geometric discounting with discounting factor λ, under the same conditions the optimal sample size is O(N∗1/2) as N∗→∞ where N∗=(1−λ)−1. This extends previous work by Cheng et al. (2003).
Although we have considered general functions hi(ξi) and gi(ξi), giving the gain to patients inside and outside the trial who receive treatment i, i=1,2, we have assumed that these are functions of ξi only. This is a common assumption and it seems reasonable that the benefit to a patient from taking a given treatment will depend only on the properties of that treatment (see, e.g. Lindley, 1997, who cites the seminal work by Raiffa and Schlaifer, 1961). Noting, however, that the gain functions gi correspond to gain from future patients if the trial indicates that treatment i is superior, some authors have proposed gain functions for patients treated following the trial that depend, in addition to ξi, on the observed trial data, Y. In particular, gain functions have been proposed that reflect the fact that use of a novel treatment following a trial may depend on regulatory decisions that in turn depends on whether trial results are sufficiently compelling (see, e.g. Posch and Bauer, 2013). In both of the more realistic examples described above, the gain functions given by Stallard (1998) and Berry et al. (1994) depended on Y, and we have simplified the gain functions when discussing these examples above.
The forms of the utility functions hi and gi,i=1,2 given above were motivated by consideration of the gain to each patient from participation in the trial or from being treated with treatment i following the trial, suggesting that the trial sample size is optimised from the patient's perspective. The general form of the expected gain given by (3), however, can express any gain so long as this can be specified on a per‐patient basis. The results obtained could thus also apply to financial gains from a commercial perspective or to societal gains from development of a novel therapy. In the latter cases it may be more appropriate for hi and gi to have a more complex form or to depend on trial data as discussed above.
The condition (2) ensures that the gain per patient in the trial does not exceed that per patient outside the trial if patients were to receive optimal treatment. If this does not hold the optimal design will be to continue with trial forever, giving all patients the treatment for which the prior expected gain E0(hi(ξi)) is the largest. This restriction on hi and gi seems reasonable if hi reflects not only the benefit to patients in the trial receiving treatment i, but also the cost of the trial, either in financial terms for the trial sponsor or funder or in terms of commitment by the patient, both of which may be considerable.
It is interesting to compare the optimal sample sizes obtained above with sample sizes typical for clinical trials. In particular, it might be of interest to consider the size of population for which conventional sample sizes would correspond to that of the optimal design. A method for frequentist sample size calculations for a single arm trial with a Bernoulli response is given by Fleming (1982), who shows that the sample size required for a trial with (one‐sided) type I error rate α and power 1−β to detect an improvement to a success probability of p
1 from a control success probability of p
0 as p1→p0 is (p0(1−p0)Φ−1(1−α)+p1(1−p1)Φ−1(1−β))2/(p1−p0)2. As discussed above, the form of g1(ξ1) given in the second single arm Bernoulli data example above and shown in Fig. 3 suggests that an appropriate value for p
1 might be about 0.35, since this value of ξ1 is associated with a high gain value. For α=0.025 and β=0.9, this would give a sample size of 111. Using the gain function described above, this would be optimal for a population of size, or expected population size in the case of geometric discounting, of about 3000. The prior distribution used in the example above, and also shown in Figure 3 is such that a value of ξ1 as large as 0.35 is highly unlikely, suggesting that a smaller value could be used for p
1 in the frequentist sample size calculation. The 95th percentile of the prior distribution is 0.256. To give power of 0.9 to detect a treatment effect corresponding to this value of p
1 would require a sample size of 704. This would be optimal for an expected population of size of a little over 100,000. It is important to note that even when the sample sizes are similar, the optimal designs obtained above may be very different from those obtained using the usual frequentist approach as, following the trial, a treatment is recommended depending on the posterior expected gains rather than on the basis of type I error rate control. As seen above, depending on the prior distribution, this can lead to type I error rates considerably higher than those conventionally used in large‐scale confirmatory studies. In this regard, the designs obtained are more similar to those sometimes used for early‐phase clinical trials or pilot studies (Schoenfeld, 1980; Stallard, 2012). Further comparison of frequentist and decision‐theoretic approaches is an area where further research would be of interest.
The authors have declared no conflict of interest.
The authors are grateful to two anonymous referees for their helpful comments. The work was conducted as part of the InSPiRe (Innovative methodology for small populations research) project funded by the European Union's Seventh Framework Programme for research, technological development and demonstration under grant agreement number FP HEALTH 2013‐602144.