How to correctly calculate the type I error in a two-step clinical trial when stopping after the first step?

Question

I just finished reading the following article by Berger & Berry (1988) in which they explain how subjectivity enters statistical analyses. One of their examples concerns a clinical trial with either one or two stages (pp. 164 in the article). The scenario looks as follows:

• we have n matched pairs of subjects.
• each member of a pair receives vitamin C and a placebo, respectively
• we want to assess whether the subject receiving vitamin C or the subject receiving the placebo exhibits greater relief from cold symptoms (it's a purely fictional experiment, so don't hang me up on the medical aspects)
• our null hypothesis is that there is no difference between vitamin C and the placebo (i.e. p(vitamin C helps) = p(placebo helps) = 1/2; hence, we're dealing with a binomial distribution

1) In the first case, there is only one stage including n=17 pairs of patients. If the number of pairs with preference for vitamin C (i.e. the number of pairs for which vitamin C seemed to help) is an element of [0,1,2,3,4,13,14,15,16,17], i.e. if it is smaller than 5 or bigger than 12, the null-hypothesis shall be rejected. The summed probability of getting a number smaller than 5 or bigger than 12 under the null hypothesis is 0.049

2) In the second case, we are dealing with two stages. In the first stage, we have n=17 pairs of patient. Again, we reject the null hypothesis if the number of pairs with a preference for vitamin C is smaller than 5 or bigger than 12. However, if we are not able to reject the null hypothesis based on this criterion, we add a second stage to our trial by looking at an additional 27 pairs, for a total of 44 pairs, concluding that there is sufficient evidence against H if the total number of preferences for vitamin C is less than 16 or more than 28. The summed probability of these events is again 0.049.

Now, Berger & Berry argue that in the case that we have 13 preferences for vitamin C out of 17 pairs the p-value is different in the first case with only one planned and conducted stage compared to the second stage with two planned, but only one conducted stage.

Their explanation: "To see this, recall the basic process for arriving at a P-value. One assumes that H is true, calculates the probability of the set of possible data which would cast as much or more doubt on H than the observed data, and claims significant evidence against H if this probability is small enough. The set R* of more extreme observations in the two-stage design equals the set R of more extreme observations for the one-stage design (n=17) plus the more extreme observations at the second stage (n=44). Since R is contained in R*, it is clear that R* has a larger probability and hence is less "significant". The probability of hitting a region outside the significance threshold after 17 observations or, failing that, after 44 observations, turn out to be 0.085"

I understand the argument that R* comprises R, so R* cannot be smaller than R. I couldn't figure out, however, how the authors got the number 0.085. My approach would have been to add up the probability of getting a significant result in the first stage (0.049) and the probability of getting a significant result in the second stage (0.049) times the probability of actually getting to the second stage (1-0.049). But this gives me: 0.049 + ((1-0.049)*0.049) = 0.095599 ≠ 0.085.

Another approach was to take into consideration that in order to get to the second stage, the number of pairs preferring vitamin C in the first stage had do lie between 5 and 12 (otherwise we would have rejected the null hypothesis and stopped the experiment). Hence, the number of pairs preferring vitamin C can only lie between 5 and 39. But this information didn't help me either to get to the right result.

Could somebody explain how I can get the type I error of 0.089 for the above-described scenario?

Answer 1

The problem here is that given the fact that the first stage has not reached significance, the second stage has less chance than its theoretical 0.049 to reach significance. Intuitively one can think that preferences are more likely to be centered on half-half if stage 1 has not reached significance. If it was not the case, your computation of your real alpha would have been 1-(1-0.049)*(1-0.049). By the way this number is the type I error, not a p value (anyway, I guess you wrote p-value by inattention) : $$1-(1-P(X1))*((1-P(X2|\bar{X1})$$ where

• X1 represents the event "statistical significance reached as step 1"
• X2 represents the event "statistical significance reached as step 2"

Adding probabilities as you did is not correct, X1 and X2 are not disjoincted events. They can be true together. The tricky part, here, is to compute $(1-P(X2|\bar{X1})$. I do not know if there is a direct formula giving it. However, a simulation, or a crude calculation both gives the alpha you are looking for.

I computed both on python. I hope the code is not too hard to grasp.

Simulation

import random success=0 random.seed(44) trials=1000000 for i in range(trials): preferences=[bool(random.getrandbits(1)) for j in range(44)] countTo_17=sum(preferences[:17]) countTo_44=sum(preferences) if((countTo_17)<=4 or (countTo_17)>=13)or((countTo_44)<16 or(countTo_44)>28): success+=1 print(success/float(trials))

This gave 0.084878

Crude computation

We want to compute $(1-P(X2|\bar{X1})$. The idea is to enumerate all the possible cases giving significance and associating their probabilities. The probability of each number of preferences got at stage 1 is given by the binomial law. But one must not forget the fact that significance has not been reached at stage 1, so the probabilities have to be "normalized" by a coefficient allowing that the sum of the different oucomes reach 1. This is what coeff in the code below do. After the computation have been done, one must do the same for the number of preferences accumulated at stage 2 and multiply both probabilities.

from scipy.stats import binom probaReachingStage1=0.049 coeff=(1/(1-probaReachingStage1)) probaReachingStage2=0 for sumOfPref_stage1 in range(5,13): p1=coeff*binom.pmf(sumOfPref_stage1, 17,0.5) for sumOfPref_stage2 in range(27): TotalsumOfPreferences=sumOfPref_stage1+sumOfPref_stage2 if(TotalsumOfPreferences>28)or(TotalsumOfPreferences<16): p2=binom.pmf(sumOfPref_stage2, 27,0.5) probaReachingStage2+=p1*p2 print(probaReachingStage2) alphaFinal=1-(1-probaReachingStage1)*(1-probaReachingStage2) print(alphaFinal)

Which also gives 0.08488