Análisis de los datos del estudio de la máscara danesa de Nassim Nicholas Taleb (GLM binomial con separación completa)

Question

Hace poco, Nassim Nicholas Taleb hecho ce poste sobre la reciente Estudio sobre la máscara danesa Un ensayo controlado aleatorizado que concluyó que la proporción de infecciones por coronavirus diagnosticadas recientemente no era significativamente diferente entre el grupo con mascarillas y el grupo de control sin mascarillas (42/2392=1,8% frente a 53/2470=2,1%), basándose en un valor p de 2 colas a partir de una regresión logística.

Taleb señala que si nos centramos sólo en los casos en los que la infección se ha confirmado mediante una prueba qRT PCR obtenemos 0/2392 frente a 5/2470. Escribe "Consideremos ahora el error más obvio. ¿Cuáles son las probabilidades de obtener 0 PCR frente a 5 al azar?

La probabilidad de tener 0 realizaciones en 2392 si la media es 5/2470 es de 0,0078518, es decir, 1 en 127. Podemos volver a expresarlo en valores p, que serían muy inferiores a 0,05, es decir, muy por encima de los valores p en el intervalo de 0,21 a 0,33 que muestra el artículo. No entiendo cómo estos investigadores pasaron por alto este punto".

Para mí, este cálculo no parece tener mucho sentido, ya que ese cálculo basado en la binomial PMF ignora la incertidumbre de muestreo en 5/2470. A esto, Taleb respondió "Yo no hago valores P: https://arxiv.org/pdf/1603.07532.pdf ", sólo para añadir más tarde este cálculo basado en el doble Bernoulli de Monte Carlo, en notación Mathematica, en el que obtiene un valor p de 1 cola de

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

¿Podría alguien arrojar luz sobre lo que Taleb estaba calculando exactamente aquí? Si él está tratando de lograr una prueba binomial de 2 muestras utilizando Monte Carlo aquí no entiendo muy bien de dónde viene el 2400 y por qué no hay Bernoulli allí con una expectativa 0/2392 (que en sí mismo sería problemático como uno tendría varianza cero entonces). Para una prueba binomial MC de 2 muestras hubiera esperado algo como (en R, y usando un ajuste +1 de todos los recuentos para evitar la expectativa binomial p=0 en un grupo) :

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

pero parece que en su lugar intenta algo como (he corregido el 2400, que presumiblemente era un error tipográfico, y también he cambiado el > por un >=):

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

lo cual creo que está mal, ¿no? En todo caso, me habría parecido más lógico :

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

aunque no estoy seguro de si ésta sería una forma aceptada de llevar a cabo una prueba binomial de 2 muestras de este tipo (parecería una especie de versión Monte Carlo de Prueba de la tabla de contingencia 2x2 de Liddell ).

El propio Taleb no fue de mucha ayuda, señalando que sólo estábamos calculando una "distribución conjunta de tabla de doble columna a la Fisher" y luego comentando en su típico estilo brusco "Pareces ser muy muy ignorante, repitiendo fórmulas como un loro, sin entender de qué va la probabilidad. He dejado de entablar conversación con usted".

Sí le dije que, dado que las máscaras a priori también podrían haber empeorado las cosas (por ejemplo, cuando se utilizan mal), era más seguro utilizar pruebas de 2 colas. Y que debido a la separación completa no se puede hacer una regresión logística normal (GLM binomial) (que es lo que los autores utilizaron en su artículo), por ejemplo, en 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

Para solucionarlo, podríamos añadir 1/2 a nuestras observaciones como corrección de continuidad (equivalente a utilizar un Jeffrey's prior en un GLM binomial bayesiano, creo):

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

Entonces señalé que sería mejor hacer un regresión logística exacta por ejemplo, en 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

Y que esto también estaría muy cerca del resultado de una prueba de 2 colas Prueba exacta de Fisher basado en la distribución hipergeométrica, que también da un valor p de 2 colas de 0,06:

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

o un valor p de 1 cola de 0,03 :

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

Aunque una prueba exacta de Fisher supondría que tanto los márgenes de las filas como los de las columnas son fijos, lo que en realidad no es del todo correcto en este caso, ya que sólo los márgenes de las filas son fijos, y esto haría que una regresión logística / binomial de 2 muestras fuera más apropiada.

Otra alternativa que señalé fue una Regresión logística de Firth lo que daría un valor p de 2 colas de 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

A lo que respondió: "Por favor, no me dé bibliotecas. Por favor, dame derivaciones". (sin tener en cuenta que no existe una solución de forma cerrada ni siquiera para la solución de máxima verosimilitud de un MLG binomial).

De todos modos, ¿podría alguien aquí dar alguna opinión sobre toda esta discusión, preferiblemente desde un ángulo estadístico formal, de modo que potencialmente pudiera complacer a Taleb? Y específicamente también sobre la cuestión de la separación completa y la mejor manera de tratarla en pruebas binomiales de 2 muestras o regresiones logísticas, y cuáles serían las mejores opciones para obtener valores p exactos.

EDIT: Pensando un poco más en las posibles opciones, una prueba exacta incondicional para comparar dos proporciones binomiales independientes sería probablemente lo más correcto. Por ejemplo, utilizando la prueba de Boschloo ( 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

Aunque eso exact.test parece tener bastantes métodos diferentes y yo mismo no estoy seguro de cuál sería el mejor (en particular para el caso de recuentos bajos y un grupo con una expectativa binomial de p=0), ya que no he profundizado en los detalles de todos esos métodos. Ej. method="Z-pooled" da valores p más optimistas, más cercanos al valor p que obtengo mediante el método MC de Liddell para probar las proporciones frente a un común p=5/(2392+2470) arriba:

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

Del mismo modo, el uso de library(exact2x2) y utilizando method="FisherAdj" Estoy obteniendo valores p más optimistas :

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

Se agradecerán las opiniones sobre cuál de estas pruebas sería la más adecuada en este caso.

Por otra parte, si se tuvieran en cuenta los falsos negativos de las pruebas PCR (del mismo modo que a Taleb le gusta tener en cuenta los falsos positivos de las pruebas de anticuerpos), las conclusiones podrían cambiar bastante... También estoy bastante seguro de que habría que saber qué individuos fueron sometidos a cada tipo de prueba, y cuáles fueron los recuentos de todos los demás virus respiratorios a los que se les hicieron pruebas.

Answer 1

Su prueba de 2 caras asigna implícitamente exactamente la mitad de su nivel de significación del 5% a "las máscaras son perjudiciales" ( $M_-$ ) y la otra mitad a "las máscaras son beneficiosas" ( $M_+$ ). Para un bayesiano como Taleb eso podría sugerir que no estás pensando adecuadamente en tu prior, porque implica que la cantidad de evidencia que te llevaría a aceptar $M_-$ es exactamente igual a la cantidad que le costaría aceptar $M_+$ aunque $M_+$ es intuitivamente más probable (al menos para mí: si hace un año me hubieran preguntado si llevar mascarilla aumentaría o disminuiría el riesgo de contraer un virus respiratorio, habría dicho que disminuiría).

Y me parece que, dado que la única variable independiente era binaria, utilizar la regresión logística de Firth en lugar de la prueba exacta de Fisher no añade mucho valor.

Pero su argumento principal era que la mayor parte de la duda de que el uso de mascarillas redujera las infecciones por PCR en el ensayo proviene de la presencia de incertidumbre en la estimación $p=\frac{5}{2470}$ . Esto parece indiscutible.

Taleb dice: "No lo entiendes. BTW I added a double column table joint distribution a la Fisher". Quizá sea su forma de salvar las apariencias admitiendo que has planteado una crítica válida.

Answer 2

Esta simulación es un intento de estimar el resultado de una prueba exacta en la que la probabilidad de observar una disparidad tan extrema sea

$$\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}$$

Este cálculo se justifica por la presunción de que todos los $2470+2392$ Los sujetos se distribuyeron de forma aleatoria e independiente en los grupos de prueba y de control. Los términos del numerador cuentan las formas en que los cinco resultados positivos podrían haber acabado aleatoriamente en el primer grupo o en el segundo. El denominador cuenta todos los posibles subconjuntos de cinco elementos del grupo de todos los sujetos, todos los cuales son igualmente probables bajo este supuesto de aleatorización.

Para los que sostienen (incorrectamente, en mi opinión) que aquí procede una prueba unilateral, he mostrado los valores de cada una de las dos fracciones por separado.

Creo que este sencillo resultado es equivalente tanto a la prueba de Fisher como a la de Boschloo.

También merece la pena señalar que el error estándar en el cálculo de Monte Carlo de Taleb es

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

teniendo en cuenta que parece haberse aproximado $2392$ por $2400$ (lo que supone poca diferencia). Esto sitúa su resultado de $0.03483$ aproximadamente $1.9$ errores estándar por debajo del valor real, dentro de lo que cabría esperar de una simulación de este tipo.

Answer 3

Si quieres consentir
(i) evitar los métodos que hacen p-valores
(ii) elaborar una prueba más "a medida" para este problema
(objetivos mencionados por Taleb, no necesariamente algo con lo que esté de acuerdo)

Una solución sería simular y hallar los parámetros mediante muestreo de rechazo (en realidad, cálculo bayesiano aproximado).

Así que digamos que hay dos parámetros clave: (a) la tasa de infección cuando sin máscara $r$ y (b) el % de cambio en la tasa de infección con mascarilla $m$ .
Pongamos dos priores uniformes $r \sim U [.0001,.005]$ y $m \sim U[.5,1.5]$ a priori crees que la máscara puede reducir el riesgo en un 50% o aumentarlo en un 50% (dos colas como sugerías).

Podemos realizar una simulación en la que extraigamos esos dos parámetros, los utilicemos para muestrear 2392 observaciones con máscara y 2470 sin ella. Sin embargo, sólo aceptamos la simulación si arroja exactamente 5 resultados positivos sin máscara y 0 resultados positivos con máscara. Seguimos haciendo esto hasta reunir, digamos, 5000 simulaciones.

Cuando observamos los valores posteriores, la distribución del multiplicador se ha desplazado claramente hacia la izquierda y se aleja del centro original en 1. Sin embargo, sigue siendo cierto que la probabilidad de que el multiplicador sea superior a 1 es de aproximadamente el 17% (muy por encima del umbral del 5% que nos habíamos fijado). ¡Todo el mundo gana!

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)))
  }

}