# Analysis of Danish mask study data by Nassim Nicholas Taleb (binomial GLM with complete separation)

Recently, Nassim Nicholas Taleb made this post about the recent Danish mask study, a randomized controlled trial which concluded that the proportions of newly diagnosed coronavirus infections was not significantly different among the group with masks and the non-mask wearing control group (42/2392=1.8% vs. 53/2470=2.1%), based on a 2-tailed p value from a logistic regression.

Taleb points out that if you focus just on the cases where infection has been confirmed by a qRT PCR test you get 0/2392 vs. 5/2470. He writes "Now consider the more obvious error. What are the odds of getting 0 PCRs vs 5 from random?

The probability of having 0 realizations in 2392 if the mean is 5/2470 is 0.0078518, that is 1 in 127. We can reexpress it in p values, which would be well <.05, that is far in excess of the ps in the range .21-.33 that the paper shows. How these researchers missed this point is beyond me."

To me, this calculation doesn't seem to make much sense, since that binomial PMF-based calculation ignores the sampling uncertainty on 5/2470. To this, Taleb responded "I don’t do P values: https://arxiv.org/pdf/1603.07532.pdf", only to later add this Monte Carlo double Bernoulli based calculation, in Mathematica notation, in which he obtains a 1-tailed p value of

ta=Table[data1=RandomVariate[BernoulliDistribution[5/2470],2470]//Total;
data2=RandomVariate[BernoulliDistribution[5/2470],2400]//Total;
data1-data2,{10^5}];
[[Select[ta,#<-5&]//Length] / 10^5]//N
0.03483


Could somebody shine a light though on what exactly Taleb was calculating here? If he is trying to accomplish a 2-sample binomial test using Monte Carlo here I don't quite understand where the 2400 is coming from and why there is no Bernoulli in there with a 0/2392 expectation (which in itself would be problematic as one would have zero variance then). For a 2-sample MC binomial test I would rather have expected something like (in R, and using a +1 adjustment of all counts to avoid the p=0 binomial expectation in one group) :

p1=rbinom(1E8, 2470, (5+1)/(2470+1))/(2470+1)
p2=rbinom(1E8, 2392, (0+1)/(2392+1))/(2392+1)
mean(p1<=p2) # 1-tailed p = 0.03401019
2*mean(p1<=p2) # 2-tailed p = 0.06802038


but it seems he instead tries something like (I corrected the 2400, which presumably was a typo, and also changed the > to a >=):

mean((rbinom(1E7, 2470, 5/2470)-rbinom(1E7, 2392, 5/2470))>=5)
# 1-tailed p = 0.0811906


which I believe is just wrong, right? If anything, I would have found this more logical :

mean((rbinom(1E8, 2470, 5/(2470+2392))/2470-rbinom(1E8, 2392, 5/(2470+2392))/2392)>=(5/2470-0/2392))
# 1-tailed p = 0.01446185
mean(abs((rbinom(1E8, 2470, 5/(2470+2392))/2470-rbinom(1E8, 2392, 5/(2470+2392))/2392))>=(5/2470-0/2392))
# 2-tailed p = 0.03479425


though I am not sure if this would be an accepted way to carry out such a 2-sample binomial test (it would seem sort of the Monte Carlo version of Liddell's 2x2 contingency table test).

Taleb himself was not very helpful, noting that we was just calculating a "double column table joint distribution a la Fisher" and then remarking in his typical blunt style "You seem to be very very ignorant, repeating formulas like a parrot, not understanding what probability is about. I’ve stopped engaging with you.".

I did tell him that because a priori masks could also have made things worse (e.g. when badly used) it was safer to use 2-tailed tests. And that because of complete separation you can't do a regular logistic regression (binomial GLM) (which is what the authors used in their paper), e.g. in R :

summary(glm(cbind(pcrpos, pcrneg) ~ treatment, family=binomial, data=data.frame(treatment=factor(c("masks","control")),pcrpos=c(0,5), pcrneg=c(2392,2470-5))))
# 2-tailed p = 1, obviously not correct


To solve this, we could add 1/2 to our observations as a continuity correction (equivalent to using a Jeffrey's prior in a Bayesian binomial GLM I believe):

summary(glm(cbind(pcrpos+1/2, pcrneg+1/2) ~ treatment, family=binomial, data=data.frame(treatment=factor(c("masks","control")),pcrpos=c(0,5), pcrneg=c(2392,2470-5))))
# 2-tailed p = 0.11


I then pointed out that it would be better to do an exact-like logistic regression, e.g. in R:

library(elrm)
fit = elrm(pcrpos/n ~ treatment, ~ treatment, r=2, iter=400000, burnIn=1000, dataset=data.frame(treatment=factor(c("masks", "control")), pcrpos=c(0, 5), n=c(2392, 2470)) )
fit$$p.values # 2-tailed p value = 0.06 fit$$p.values.se # standard error on p value = 0.0003


And that this would also be very close to the result of a 2-tailed Fisher exact test, based on the hypergeometric distribution, which also gives a 2-tailed p value of 0.06:

fisher.test(rbind(c(0,2392), c(5,2470-5))) # 2-tailed p value = 0.06


or a 1-tailed p-value of 0.03 :

fisher.test(rbind(c(0,2392), c(5,2470-5)), alternative="less") # 1-tailed p value = 0.03


Even though a Fisher exact test would assume both the row and column margins to be fixed, which is in fact not quite correct here, as only the row margins are fixed, and this would make a logistic regression / 2-sample binomial more appropriate.

Another alternative that I pointed out was a Firth's logistic regression, which would give a 2-tailed p value of 0.11 :

library(brglm)
summary(brglm(cbind(pcrpos, pcrneg) ~ treatment, family=binomial, data=data.frame(treatment=factor(c("masks","control")), pcrpos=c(0,5), pcrneg=c(2392,2470-5))))
# 2-tailed p = 0.11


To that he responded "Please don’t give me libraries. Please provide derivations." (never mind that there is no closed-form solution even for the maximum likelihood solution of a binomial GLM).

Anyway, would somebody here be able to give some feedback on this whole discussion, preferably from a formal statistical angle, so that it could potentially please Taleb? And specifically also on the issue of complete separation and how to best deal with it in 2-sample binomial tests or logistic regressions, and what the best options would be there to obtain exact p values.

EDIT: Thinking about the possible options a bit more, an exact unconditional test to compare two independent binomial proportions would probably be most correct. E.g. using Boschloo's test (https://en.wikipedia.org/wiki/Boschloo%27s_test):

library(Exact)
exact.test(rbind(c(0,2392), c(5,2470-5)), method="Boschloo", alternative="two.sided", model="Binomial")
# Boschloo's test, 2-tailed p = 0.06223
exact.test(rbind(c(0,2392), c(5,2470-5)), method="Boschloo", alternative="less", model="Binomial")
# Boschloo's test, 1-tailed p = 0.03196


Though that exact.test function seems to have quite a lot of different methods, and I am unsure myself which would be best (in particular for the case with low counts and a group with a binomial expectation of p=0), as I haven't dug into the details of all those methods. E.g. method="Z-pooled" gives more optimistic p values, closer to the p-value I get via the Liddell-like MC method to test the proportions against a common p=5/(2392+2470) above:

exact.test(rbind(c(0,2392), c(5,2470-5)), method="Z-pooled", alternative="two.sided", model="Binomial")
# 2-tailed p-value = 0.02809
exact.test(rbind(c(0,2392), c(5,2470-5)), method="Z-pooled", alternative="less", model="Binomial")
# 1-tailed p-value = 0.01425


Likewise, using library(exact2x2) and using method="FisherAdj" I am getting such more optimistic p values :

uncondExact2x2(0, 2392, 5, 2470, alternative="two.sided", method="FisherAdj")
# 2-tailed p = 0.03417
uncondExact2x2(0, 2392, 5, 2470, alternative="greater", method="FisherAdj")
# 1-tailed p = 0.01709


Thoughts on which of these tests would be most appropriate here would be welcome.

On another note, if one would take into account false-negative PCR tests (in the same way that Taleb likes to take into account false-positives in the antibody tests) this might change conclusions quite a bit... Also pretty sure that one would need to know which individuals were subjected to each type of test, and what the counts were for all the other respiratory viruses they were tested for.

• Pleasing Taleb is probably too tall an order; I hope you will consider answers that meet more reasonable demands. Dec 4, 2020 at 11:54
• @Forgottenscience Yes, very hard to engage in a discussion with him. Seems he just likes to throw random insults at you, even when you make sensible comments. Is there a name for that condition? :-) Dec 4, 2020 at 12:11
• Taleb Syndrome?
– hd1
Dec 5, 2020 at 23:06
• Just a note for others - the study also tested for several other respiratory viruses via PCR finding 9 in the mask group and 11 in the non-mask group. So even if this study suggests mask protect against covid, it also suggests they do not protect against other respiratory viruses. Dec 14, 2020 at 13:50

Your 2-sided test implicitly allots exactly half of your 5% significance level to "masks are harmful" ($$M_-$$) and the other half to "masks are beneficial" ($$M_+$$). To a Bayesian like Taleb that might suggest that you aren't properly thinking about your prior, because it implies that the amount of evidence it would take you to accept $$M_-$$ is exactly equal to the amount it would take you to accept $$M_+$$, even though $$M_+$$ is intuitively more likely (to me, at least - if you had asked me a year ago whether wearing a mask would be more likely to increase or decrease the risk of catching a respiratory virus, I would have said decrease).

And it seems to me that since the single independent variable was binary, using Firth's logistic regression over Fisher's exact test doesn't add much value.

But your core point was that most of the doubt that mask wearing reduced PCR infections in the trial comes from the presence of uncertainty in the estimate $$p=\frac{5}{2470}$$. This seems unarguable.

Taleb says: "You don’t get it. BTW I added a double column table joint distribution a la Fisher." Maybe that's his way of saving face while admitting that you raised a valid criticism.

• Thanks for this - that makes sense! Injecting priors in medical trials in general is quite tricky, because of the inherent subjectivity. Though I still quite didn't get what he was doing in that double Bernoulli MC simulation - how come the 0/2392 expectation didn't feature there, and where did his 2400 come from? Any thoughts on that? And was he attempting to do a 2-sample binomial test there? Dec 4, 2020 at 16:56
• I think he was estimating the probability of a difference of at least 5 between 2 samples, assuming a population mean of 5/2470 in both. (So not quite a traditional 1-sample or 2-sample binomial test.) It could be intended as a back-of-the-envelope estimation, with the implication that the estimation's result is significant enough that it should be obvious that accounting for the sampling error in the usual way wouldn't be enough to push the p-value above 5%. Dec 4, 2020 at 17:40
• And my guess is that the 2400 is him rounding 2392 to 2 significant figures, on the assumption that this wouldn't change the result too much. Dec 4, 2020 at 17:41
• Ha thanks - of course, yes, I get it now... Still very "back of the envelope" though indeed. And that rounding if 2392 to 2400 is very weird. He also has another "back of the envelope" approximate 1-tailed p value there, calculated as 0.5^5 = 0.03125 (ignoring the difference in sample size between both groups). A bit strange these kind of heuristic p values in a post critcising another study claiming they didn't get their hypothesis tests & p values right... Dec 4, 2020 at 20:46
• Also, you would expect that if Taleb really wanted to do it "a la Fisher" that he would just have done a = 0; b = 2392; c = 5; d = 2470 - 5; n = a + b + c + d; Binomial[a + b, a]*Binomial[c + d, c]/Binomial[n, b + d] // N, 1-tailed p p=0.033771! Dec 4, 2020 at 20:46

This simulation is an attempt to estimate the result of an exact test in which the chance of observing such an extreme disparity is

\begin{aligned} \Pr(\text{all positive results in one group}) &= \frac{\binom{2470}{5} + \binom{2392}{5}}{\binom{2470+2392}{5}} \\&= 0.03377 + 0.02876 = 0.06253. \end{aligned}

This calculation is justified by the presumption that all $$2470+2392$$ subjects were independently randomized into the test and control groups. The terms in the numerator count the ways in which all five positive results could have ended up randomly in the first group or the second group. The denominator counts all possible five-element subsets of the group of all subjects, all of which are equally likely under this randomization assumption.

For those who maintain (incorrectly, IMHO) that a one-sided test is appropriate here, I have shown the values of each of the two separate fractions.

I believe this simple result is equivalent to both Fisher's Test and Boschloo's test.

It might also be worth noting that the standard error in Taleb's Monte Carlo calculation is

$$\operatorname{se} = \sqrt{0.03377(1-0.03377)/10^5} = 0.0005712,$$

bearing in mind he appears to have approximated $$2392$$ by $$2400$$ (which makes little difference). This places his result of $$0.03483$$ approximately $$1.9$$ standard errors below the true value, well within what one would expect of such a simulation.

• When you say "I believe this simple result is equivalent to both Fisher's Test and Boschloo's test" - that can't be quite correct right, since both tests return different p values, among others because Fisher's exact test is based on the hypergeometric distribution and assumes both column and row margins are fixed, while Boschloo's test only assumes that row margins are fixed and is derived on the assumptions that both groups conform to independent binomial distributions. Dec 5, 2020 at 22:18
• You're right: I shouldn't have said "equivalent." These tests are all made in a similar spirit, but with slightly different models. The model I use here is explicitly the model most people use in controlled studies: they randomly and independently put people into treatment and control groups.
– whuber
Dec 6, 2020 at 2:20
• Though what you write makes sense, I feel that this cannot be what he tried to calculate, since he also presents a similar calculation for a comparison 5/2393 vs 15/2470. I think what he really meant was maybe some MC version of Liddell's 2x2 contingency table test, rss.onlinelibrary.wiley.com/doi/abs/10.2307/2988087, i.e. mean((rbinom(1E8, 2470, 5/(2470+2392))/2470-rbinom(1E8, 2392, 5/(2470+2392))/2392)>=(5/2470-0/2392)) # 1-tailed p=0.01446 or mean(abs((rbinom(1E8, 2470, 5/(2470+2392))/2470-rbinom(1E8, 2392, 5/(2470+2392))/2392))>=(5/2470-0/2392)) # 2-tailed p=0.03479 Dec 6, 2020 at 21:30

If you want to acquiesce to
(i) avoiding methods that do p-values
(ii) producing a more "bespoke" test for this problem
(objectives mentioned by Taleb, not necessarilly something you agree with)

One solution would be to simulate and find parameters by rejection sampling (actually approximate Bayesian computation).

So let's say there are two key parameter: (a) the infection rate when maskless $$r$$ and (b) the % change in infection rate wearing a mask $$m$$.
Let's put two uniform priors $$r \sim U [.0001,.005]$$ and $$m \sim U[.5,1.5]$$: a priori you think mask can cut risk by 50% or increase it by 50% (two-tailed as you suggested).

We can run a simulation where we draw those two parameters, use them to sample 2392 observations with mask and 2470 without it. However we accept the simulation only where we it outputs precisely 5 positives when wearing no masks and 0 positives when wearing it. We keep doing this until we collect, say, 5000 simulations.

When we look at the posteriors, the multiplier distribution has clearly shifted left and away from the original center at 1. However it is still true that the probability that the multiplier is above 1 is about 17% (far above the threshold of 5% we set ourselves). Everybody wins! library(tidyverse)
library(progress)

### hard data from the paper

### here assume uniform priors on
### chances of getting sick
runif(n=1,min=.0001,max=.005)
}
### let's assume we have a uniform prior that masks can do anything
### between cutting risks by 50% and increasing it by 50%
function(){
runif(n=1,min=.5,max=1.5)
}

### run simple ABC where we look for parameters
### where we get exactly the observations above
accepted_runs<-list()
TARGET_ACCEPTED_RUNS<-5000
pb <- progress_bar$$new(total = TARGET_ACCEPTED_RUNS) ## to watch the time go by while(length(accepted_runs)tick()
accepted_runs<-
append(accepted_runs,