# Sample size calculation in COVID-19 study

From A Randomized Trial of Hydroxychloroquine as Postexposure Prophylaxis for Covid-19 by Boulware et al. in the New England Journal of Medicine (https://www.nejm.org/doi/full/10.1056/NEJMoa2016638), I'm curious about the following sentence when it comes to appropriate sample size calculations:

Using Fisher’s exact method with a 50% relative effect size to reduce new symptomatic infections, a two-sided alpha of 0.05, and 90% power, we estimated that 621 persons would need to be enrolled in each group.

I'm interested in how this calculation is performed. I've never heard of "effect size" being used in the context of Fisher's exact test (I'm familiar with Coehn's $$d$$), and I'm not exactly sure how power calculations would work in this case (what's an appropriate alternative hypothesis?).

Keep in mind that I have zero expertise in clinical trials. I'm very comfortable with statistics at the level of Casella and Berger's text.

Textbooks and journal articles would be very helpful to have for further study.

• "relative effect size" doesn't mean anything. is it relative risk? risk difference? odds ratio? Jun 7 '20 at 17:34
• anyway, the alternative hypothesis is having that 50% "relative effect size", wathever it means. Jun 7 '20 at 17:36

A glib answer is that they probably just plugged their numbers into a power calculator. I've attached a screenshot re-creating this power analysis in G*Power 3.1, a freely available power calculator. Note to match their result of 621 I had to go to "Options" and select "Maximize Alpha".

The paper says "We anticipated that illness compatible with Covid-19 would develop in 10% of close contacts exposed to Covid-19" as well as "50% relative effect size". I interpret the second part to mean they assume that the effect of treatment will reduce the illness rate from 10% to 5%.

This leads to the values of $$0.05$$ and $$0.1$$ for Proportions p1 and p2 respectively.

Sadly I don't know how G*Power makes this calculation, but I can attempt to explain the idea at least.

We are given our proportions of 0.1 and 0.05. For a given sample size $$n$$, we can randomly sample a 2x2 contingency table by sampling from two binomial random variables. The power calculation asks, "how often will Fischer's Exact Test reject the null hypothesis for a contingency table created using this process?".

In particular, we want to find the smallest $$n$$ such that Fischer's test will reject the null hypothesis at least 90% of the time.

One way to approximate this is with simulation. For a given $$n$$, sample say 10,000 contingency tables, run Fischer's test, and see how often the p-value is below 0.05. Keep increasing $$n$$ until the p-value is below 0.05 90% of the time or more... I know I am several months late, but just want to respond to the other answers. All answers use simulations and/or claim the exact Fisher calculation is too computationally intensive. If you code this efficiently, you can get an exact computation very quickly. Below is a comparison time of the sample code fisherpower() function vs. the power.exact.test() function in the Exact R package:

    > system.time(power1 <- fisherpower(0.1,0.05,621))
user  system elapsed
698.23    0.93  700.23
> system.time(power2 <- Exact::power.exact.test(n1=621,
n2=621, p1=0.1, p2=0.05,
method="Fisher")$power) user system elapsed 0.32 0.00 0.33 > power1  0.9076656 > power2  0.9076656  The calculation only takes 0.33s using power.exact.test() function compared with 700s using the fisherpower() function. Note the power.exact.test() function computes the exact power without simulations, so there's no uncertainty and it's faster than simulating. I also strongly recommend using Barnard's exact test over Fisher's exact test for comparing two proportions. Below is the power calculation as the group sample size increases:  nGroup <- 570:630 powerFisher <- vapply(nGroup, FUN = function(xn) { Exact::power.exact.test(n1=xn, n2=xn, p1=0.1, p2=0.05, method="Fisher")$$power }, numeric(1) ) powerBarnard <- vapply(nGroup, FUN = function(xn) { Exact::power.exact.test(n1=xn, n2=xn, p1=0.1, p2=0.05, method="Z-pooled")$$power }, numeric(1) ) plot(NA, xlim=range(nGroup), ylim = c(0.85,0.95), xlab="Sample Size per Group", ylab = "Power") lines(nGroup, powerFisher, col='red', lwd=2) points(nGroup, powerFisher, pch = 21, col = 'red', bg = "red", cex = 0.8) lines(nGroup, powerBarnard, col='blue', lwd=2) points(nGroup, powerBarnard, pch = 21, col = 'blue', bg = "blue", cex = 0.8) abline(h=0.9, lty=2) abline(v=c(579, 606), col=c('blue', 'red')) legend(610, 0.875, c("Barnard", "Fisher"), col = c('blue', 'red'), lty = 1, pch=21, pt.bg=c('blue', 'red'), cex=1.2) @heropup is correct that the group sample size should be 606 (not 621) as shown in the figure. However, Barnard's test is more powerful and only requires 579 participants in each group using the "Z-pooled" test statistic. Since this is a rare event, one may want to use the Berger and Boos (1994) interval approach, which brings the sample size down to 573 participants (code not shown, takes some time). Importantly, these alternatives still control for the type 1 error rate and are simply superior to Fisher's exact test for 2x2 tables. For analyzing the dataset, I would recommend using Exact::exact.test() which takes only 0.3s for the example dataset that @SextusEmpiricus provided instead of Barnard::barnard.test() which takes 47s. However, both yield the same results and I am the maintainer of the Exact R package so may be biased. • This fisherpower() is a function written for the didactical purpose to demonstrate how such calculation is performed (maybe I should have written it as a formula instead of sloppy code). But how does power.exact.test() do it? Is it just smarter execution of the computations, like not repeating the calculations of the binomial distribution. Or is it something more? Sep 29 '20 at 5:39 • There are three pieces to implement this efficiently. 1) Implement Fisher's exact test quickly. For example, simply adding conf.int = FALSE to your fisherpower() function cuts the computation time down to 336s (I created my own faster test for 2x2 tables). 2) Utilize symmetry of the rejection region; this halves the number of tables considered. 3) Utilize convexity property whereby knowing one test is rejected can imply other more extreme tables are rejected without needing to calculate p-value. Altogether, I ended up calculating 1237 p-values to determine the entire rejection region. Sep 29 '20 at 22:21 They used Fisher's exact test, which relates to sampling without replacement. But in reality this is not exactly like that, and it is more like binomial distributed data. For that case you get the following: • For the null hypothesis it is sampling where you have equal probabilities that the people get covid-19, no matter whether they are in the placebo group or the effect group. • The alternative hypothesis, for which they computed the power, is that you have 10% probability for the placebo group to get covid-19 and 5% for the treatment group (so the treatment reduces probability by 50%). ### Exact computation of power You can compute the probability of rejecting the null hypothesis given certain sample sizes and probabilities simply by trying all possibilities and see which ones result in a negative/positive fisher's test. Then you sum over the probabilities to get the cases where you reject the test. $$P(\text{reject})= \sum_{\substack{\text{over all i,j} \\ \text{where Fisher test is rejected}}} P(\text{i placebo cases and j treatment cases})$$ Below is a code example fisherpower <- function(p1, p2, n) { pf <- 0 for (i in 1:n) { for (j in 1:n) { M <- matrix(c(i,n-i,j,n-j),2) if (fisher.test(M)$p.value <= 0.05) {
pf <- pf + dbinom(i,n,p1)*dbinom(j,n,p2)
}
}
}
pf
}


which gives

> fisherpower(0.1,0.05,621)
 0.9076656


However, this method costs a lot of computation power. You need to try 621 times 621 possibilities. The implementation above can be improved a lot (you do not need to compute all 621 by 621 cases), but it will remain slow so that is why the standard implementation in R uses simulations. A fast implementation of the above is in Peter Calhoun's R package Exact which he explains in his answer here.

### Computing with simulations

You compute multiple times a hypothetical outcome and for that outcome you determine whether a 5% hypothesis test will fail or not.

As a function of the sample sizes you get:

• If the null hypothesis is true then you will always get 5% probability to reject.

Actually this is not entirely true and the Fisher's exact test is slightly conservative when the conditioning is not correct. The Fisher's exact test will reject less often than 5%, even when the null hypothesis is true (in the case that we do not sampling with replacement). In the example graph below we compute the rejection probability when $$p_1 = p_2 = 0.1$$ (in which case the null is true).

• If the null hypothesis is false, and the probabilities are not equal. Then you get a larger probability to reject the null hypothesis when the sample size is larger. ### computing
set.seed(1)
n <- seq(100,1000,20)
power <- sapply(n,
FUN = function(xn) {
statmod::power.fisher.test(0.1,0.05,xn,xn, nsim = 10000)
} )
type1 <- sapply(n,
FUN = function(xn) {
statmod::power.fisher.test(0.1,0.1,xn,xn, nsim = 10000)
} )

### plotting of results
plot(n,power, type = "l", ylim = c(0,1),
ylab = "reject probability")
lines(n,type1, col =2)
points(n,power, pch = 21, col = 1, bg = "white", cex = 0.7)
points(n,type1, pch = 21, col = 2, bg = "white", cex = 0.7)

# lines at 0.05 and 0.9
lines(c(0,2000),c(0.05,0.05), col = 2, lty = 2)
lines(c(0,2000),c(0.9,0.9), col = 1, lty = 2)

# legend
legend(1000,0.6,c("if p1 = p2 = 0.1",
"if p1 = 0.1, p2 = 0.05"), title = "reject probability",
col = c(2,1), lty = 1, cex = 0.7, xjust = 1
)


### Alternative tests

There are many other ways to look at it. We can also perform a Barnards test

> Barnard::barnard.test(49,58,414-49,407-58)

Barnard's Unconditional Test

Treatment I Treatment II
Outcome I           49           58
Outcome II         365          349

Null hypothesis: Treatments have no effect on the outcomes
Score statistic = 1.02759
Nuisance parameter = 0.012 (One sided), 0.986 (Two sided)
P-value = 0.16485 (One sided), 0.320387 (Two sided)


or use a GLM model

> summary(glm(cbind(c(49,58),c(414-49, 407-58)) ~ 1+c("chloroquine", "placebo"), family = binomial(link="identity")))

Call:
glm(formula = cbind(c(49, 58), c(414 - 49, 407 - 58)) ~ 1 + c("chloroquine",
"placebo"), family = binomial(link = "identity"))

Deviance Residuals:
  0  0

Coefficients:
Estimate Std. Error
(Intercept)                         0.11836    0.01588
c("chloroquine", "placebo")placebo  0.02415    0.02350
z value Pr(>|z|)
(Intercept)                          7.455 8.98e-14 ***
c("chloroquine", "placebo")placebo   1.028    0.304
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1.0568e+00  on 1  degrees of freedom
Residual deviance: 2.4780e-13  on 0  degrees of freedom
AIC: 15.355

Number of Fisher Scoring iterations: 2

>


Each of these methods shows more or less the same thing, the result 58 vs 49 is not an anomaly (but also, the effect would need to be 50% or more in order for us to have at least 90% probability to detect an anomaly with this test).

• Great answer! I added a response that improves the computation time Sep 28 '20 at 22:52

You are missing a critical piece of information that the article cited immediately prior to your quote:

We anticipated that illness compatible with Covid-19 would develop in 10% of close contacts exposed to Covid-19.

This is the assumed incidence in the control group under the alternative hypothesis; i.e., $$\pi_c = 0.1$$. The 50% relative effect size refers to a reduction in the incidence of Covid-19 infection in the treatment group, i.e. $$\pi_t/\pi_c = 0.5$$ from which it follows that $$\pi_t = 0.05$$, under the alternative hypothesis.

However, when I input these (along with $$\alpha$$ and $$\beta$$) into EAST 6, I don't get $$n = 621$$ per arm. I get $$n = 606$$ per arm, and based on my simulations, I believe the latter value is correct.