---
title: "Statistical comment on the revision of Bendavid et al."
author: "Will Fithian"
date: "5/1/2020"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


This is a comment based on my Twitter thread, about why I still think the (unweighted) data analysis does not support the conclusion that the prevalence of COVID-19 in Santa Clara County is above zero. Of course, we know based on information from outside of the study that the prevalence is higher than zero -- the point is that if the data in the study cannot even establish that, then they also cannot establish that the prevalence is large.

You can [download the Rmd file here](https://www.stat.berkeley.edu/~wfithian/overdispersionSimple.Rmd)

## New data on specificity

The main new feature in the new study is the addition of 11 more negative control samples in addition to the two from the first draft, bringing the number of samples to 13 and providing more insight into the test's *specificity*:
\[
s = P(\text{test negative} \mid \text{true negative}).
\]

I find it a bit more intuitive to discuss things in terms of the test's *false positive rate* (FPR), which is $1-s$.
\[
\text{FPR} = P(\text{test positive} \mid \text{true negative}) = 1 - s.
\]
This is a population quantity, not a sample quantity: the percent of positive tests in any particular negative control sample is a noisy estimate of that sample's true FPR.

Taken together, the 13 negative control samples have 3324 total specimens. We can copy the data into R, adding in the Santa Clara County study with 50 positives out of 3330 tests:

```{r, echo = TRUE}
dat <- data.frame(
  size     = c(371, 30, 70, 1102, 300, 311, 500, 200, 99, 31, 150, 108, 52, 3330),
  test.neg = c(368, 30, 70, 1102, 300, 311, 500, 198, 99, 29, 146, 105, 50, 3280)
)

dat$test.pos <- dat$size - dat$test.neg
dat$pct.pos <- round(dat$test.pos / dat$size * 100, 2)
dat$sample.type <- as.factor(c(rep("control", 13), "scc"))

dat
```

Right away we see something odd: four of the 13 negative control samples have a *higher* fraction of false positives than the Santa Clara County sample. Looking at the data this way, it would be easy to believe that the Santa Clara County dataset was just one more negative control sample with a false positive rate somewhere in the middle of the pack.

## A heroic assumption?

But Bendavid et al. don't look at the data this way: they just pool the 13 negative control samples into one giant sample with 0.5% false positives overall. This makes sense only if we are assuming that all 14 of the samples have the same false positive rate.[^1]

[^1]: The authors' language consistently discusses specificity as if it's a single number, but their bootstrap algorithm arguably is not assuming the FPR is the same in every sample, only that the FPR in the SCC sample is exactly the weighted average of the negative control samples. If so, that assumption wouldn't make any more sense: if the FPR is variable across samples, there's no reason SCC would be exactly the average of the others; we should just expect it to have its own FPR that's variable in the same way the other samples' FPRs are.

But that assumption doesn't fit the data. In particular, we can reject it using Pearson's $\chi^2$ test on the 13 x 2 table of negative control samples:
```{r, echo = TRUE}
(ctrl.tbl <- as.matrix(dat[dat$sample.type == "control", c("test.pos", "test.neg")]))

chisq.test(ctrl.tbl)
```

Because of the small counts we really should do a permutation test. This is computationally intensive.
```{r, echo = TRUE, cache = TRUE}
perm.groups <- as.factor(rep(1:13, dat$size[1:13]))
perm.labels <- rep(c(TRUE,FALSE), c(sum(dat$test.pos[1:13]),sum(dat$test.neg[1:13])))

perm.dist <- replicate(1E4, {
  perm.pos <- table(perm.groups[sample(perm.labels)])
  perm.tbl <- cbind(perm.pos, dat$size[1:13] - perm.pos)
  suppressWarnings(chisq.test(perm.tbl)$statistic)
})

observed.stat <- suppressWarnings(chisq.test(ctrl.tbl)$statistic)
```

We can plot the permutation distribution:
```{r, echo = TRUE}
hist(perm.dist, breaks=100, xlim=c(0, max(perm.dist, observed.stat)), freq=FALSE, 
     main = "Permutation Chi Squared Distribution",
     xlab="Chi Squared Statistic")
arrows(observed.stat, .02, observed.stat, 0, col="red")
text(observed.stat, .02, "observed", pos = 3, col = "red")
```

The permutation $p$-value:
```{r}
1 - (rank(c(observed.stat, perm.dist))[1] - 1) / (1+length(perm.dist))
```

## Modeling with overdispersed counts

The 13 samples are overdispersed relative to the binomial. One simple analysis is to fit a beta-binomial model to the negative control samples, using the [VGAM library](https://cran.r-project.org/web/packages/VGAM/VGAM.pdf):

```{r, echo = TRUE}
library(VGAM)

ctrl.bbfit <- vglm(ctrl.tbl ~ 1, family=betabinomial)

# back out
mu <- 1/(1+exp(-coef(ctrl.bbfit)[1]))
rho <- 1/(1+exp(-coef(ctrl.bbfit)[2]))
a.plus.b <- 1/rho - 1
a <- a.plus.b * mu
b <- a.plus.b * (1-mu)

# simulate FPRs
sim.fpr <- rbeta(10000, shape1 = a, shape2 = b)
sim.count <- rbinom(10000, sim.fpr, size = 3330)
hist(sim.count, breaks=100, xlim=c(0,250), freq = FALSE,
     main = "Number of false positives expected in Santa Clara County sample\n(Beta-binomial model)", xlab="False positives")
arrows(50, .02, 50, 0.005, col="red")
text(50, .02, "total observed\npositive tests", pos = 3, col = "red")


```

50 is the 77th percentile of this distribution, so there is again no way to be confident that any of the positive test kits represented true positives.

A bit more formally, there are a few simple ways to estimate binomial regression models with overdispersed counts. We can add an indicator variable of the sample type being Santa Clara County instead of control, and see if the coefficient is significant.

```{r, echo = TRUE}
full.tbl <- as.matrix(dat[, c("test.pos", "test.neg")])

full.bbfit <- vglm(full.tbl ~ sample.type, data = dat, family=betabinomial)
summary(full.bbfit)

full.qbfit <- glm(full.tbl ~ sample.type, data = dat, family=quasibinomial)
summary(full.qbfit)

library(lme4)
dat$id <- as.factor(1:14)

full.lmefit <- glmer(full.tbl ~ sample.type + (1 | id), data = dat, family = binomial)
summary(full.lmefit)
```


Note that all of this analysis still underestimates the uncertainty because it doesn't take into account the fact that some of the participants in the SCC study were parent-child pairs.

Finally, this is not an endorsement of the rest of the statistics in the paper. I am not addressing the weighted analysis at all. I would be surprised if the weighted analysis could overturn this conclusion, but without seeing the data I'm reluctant to guess.

## Pre-COVID samples

Some of the negative control samples were based on pre-COVID blood and others were based on clinically excluded patients during the pandemic. For good measure, I redo the permutation test for the Pre-COVID negative control samples only. Nothing much changes.

```{r, echo = TRUE, cache = TRUE}
precovid <- (1:13)[-c(6, 11, 13)]
precovid.perm.groups <- as.factor(rep(precovid, dat$size[precovid]))
precovid.perm.labels <- rep(c(TRUE,FALSE), c(sum(dat$test.pos[precovid]),sum(dat$test.neg[precovid])))

precovid.perm.dist <- replicate(1E4, {
  precovid.perm.pos <- table(precovid.perm.groups[sample(precovid.perm.labels)])
  precovid.perm.tbl <- cbind(precovid.perm.pos, dat$size[precovid] - precovid.perm.pos)
  suppressWarnings(chisq.test(precovid.perm.tbl)$statistic)
})

precovid.observed.stat <- suppressWarnings(chisq.test(ctrl.tbl[precovid,])$statistic)
```


```{r, echo = TRUE}
hist(precovid.perm.dist, breaks=100, xlim=c(0, max(precovid.perm.dist, precovid.observed.stat)), freq=FALSE, 
     main = "Permutation Chi Squared Distribution (Pre-COVID samples only)",
     xlab="Chi Squared Statistic")
arrows(precovid.observed.stat, .02, precovid.observed.stat, 0, col="red")
text(precovid.observed.stat, .02, "observed", pos = 3, col = "red")
```


```{r, echo = TRUE}
library(VGAM)

precovid.ctrl.bbfit <- vglm(ctrl.tbl[precovid,] ~ 1, family=betabinomial)

# back out
mu <- 1/(1+exp(-coef(precovid.ctrl.bbfit)[1]))
rho <- 1/(1+exp(-coef(precovid.ctrl.bbfit)[2]))
a.plus.b <- 1/rho - 1
a <- a.plus.b * mu
b <- a.plus.b * (1-mu)

# simulate FPRs
sim.fpr <- rbeta(10000, shape1 = a, shape2 = b)
sim.count <- rbinom(10000, sim.fpr, size = 3330)
hist(sim.count, breaks=100, xlim=c(0,250), freq = FALSE,
     main = "Number of false positives expected in Santa Clara County sample\n(Beta-binomial model)", xlab="False positives")
arrows(50, .02, 50, 0.005, col="red")
text(50, .02, "total observed\npositive tests", pos = 3, col = "red")
```
```{r, echo = TRUE}
full.tbl <- as.matrix(dat[c(precovid,14), c("test.pos", "test.neg")])

full.bbfit <- vglm(full.tbl ~ sample.type, data = dat[c(precovid,14),], family=betabinomial)
summary(full.bbfit)

full.qbfit <- glm(full.tbl ~ sample.type, data = dat[c(precovid,14),], family=quasibinomial)
summary(full.qbfit)

library(lme4)
dat$id <- as.factor(1:14)

full.lmefit <- glmer(full.tbl ~ sample.type + (1 | id), data = dat[c(precovid,14),], family = binomial)
summary(full.lmefit)
```
