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

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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
WITH_MASK_OBSERVATIONS<-2392
WITHOUT_MASK_OBSERVATIONS<-2470

WITH_MASK_POSITIVES<-0
WITHOUT_MASK_POSITIVES<-5

### here assume uniform priors on 
### chances of getting sick
PRIOR_INFECTION_RISK_MASKLESS<-function(){
  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%
PRIOR_CHANGE_RISK_MASK<-
  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)<TARGET_ACCEPTED_RUNS)
{
  ## draw a "real" infection rate from your priors
  infection_maskless<-PRIOR_INFECTION_RISK_MASKLESS()
  mask_multiplier<-PRIOR_CHANGE_RISK_MASK()
  infection_masked<- mask_multiplier * infection_maskless
  ## observe the infection rate!
  observed_infections_maskless<- sum(
rbinom(n=WITHOUT_MASK_OBSERVATIONS,
       size=1,
       prob=infection_maskless))
  observed_infections_mask<-sum(
rbinom(n=WITH_MASK_OBSERVATIONS,
       size=1,
       prob=infection_masked))
  ## if this is EXACTLY what we observe, store the drawn infection rates
  ## they will be part of our posterior
  if(observed_infections_maskless==WITHOUT_MASK_POSITIVES &&
 observed_infections_mask==WITH_MASK_POSITIVES
  ){
pb$tick()
    accepted_runs<-
      append(accepted_runs,
             list(data.frame(maskless = infection_maskless,
                             masked = infection_masked,
                             multiplier = mask_multiplier)))
  }
  
}