We don't need fraud to get a replication crisis

(A longer version of this is going in the technical notes for my book, so there are references here to it.)

We can use an analysis like the one in [book chapter] to understand the so-called ``replication crisis,'' of many research papers failing to replicate.

We'll sidestep the controversies about suspicion of fraud in specific results and show how the same random clusters that happen in car crashes and diseases can lead to this replication crisis (i.e. they create positive results with data that doesn't show them).

(Some authors of the papers that failed to replicate have admitted to various levels of fraud and several papers have been retracted due to strong indications of fraud. We're not disputing that fraud exists. But here we show that large numbers of well-intentioned researchers in a field, combined with a publication bias towards positive results, can generate false results just by probability alone.)

As we saw in the Statistics Are Weird Interlude, when using a sample to make judgments about the population, a traditional (frequentist) way to express confidence in the results is by testing them and reporting a confidence level. For historical reasons many fields accept 95% as that confidence level.

Reminder: The general public assumes that 95% confidence level means that the phenomenon is true with probability $0.95$, given the test results. That's not how it works: the 95% confidence means that the probability that we get the test result when the phenomenon is true is $0.95$. In other words, the public perception is a probability of the truth of the phenomenon conditional on the result of the test, the real meaning is a probability of the result conditioned on the truth of the phenomenon, and in general

\[ \text{[perception] } \Pr(\text{phenomenon}|\text{result}) \neq \Pr(\text{result}|\text{phenomenon}) \text{ [reality]}\]

There's a schism in statistics about whether and how to use uninformative priors and chained Bayes rule ---as described in the Bayesian Interlude--- to deal with multiple experiments (the frequentists' approach being ``not to use them at all''). We'll sidestep those discussions here, but point out their implications further down. End of reminder.

So, the probability that we get a false positive from testing a data set for a non-existent phenomenon at the 95% confidence level is $p \le 0.05$ and we'll use the equal value for our modeling, $p = 0.05$.

Let's say we have a team of scientists investigating the very important phenomenon of whether liking pineapple on pizza makes a person more likely to put real cream, rather than non-dairy creamer, in their coffee.

To investigate this important matter, the team chooses a sample of the population and runs some controlled experiment. Maybe the result is positive, or maybe it's negative; but since they have funding, the scientists can run more experiments, controlling for different variables (is the pizza deep-dish? do people put honey or sugar in the coffee with their cream/creamer? these are called covariates or, sometimes, controls). Say the team runs 100 experiments in total.

Now we make our important assumption, which we can do because we always play in God-mode: in reality there's no relationship between pineapple on pizza and the choice of cream or creamer in coffee.*

Because of that assumption, all the positives in the data will be false positives, and we know that those happen with $p = 0.05$, by our setup. So, in 100 experiments there should be an average of 5 false positives, and the number of false positives itself is a random variable, say $N$, distributed Poisson with parameter 5. The next figure shows that distribution.

We can also ask some simple questions: how likely is the team to observe at least, say 7 positive results? How about 10? These are easily calculated:

\[ \Pr(N\ge 7) = 1-F_N(6) = 0.24 \quad\text{and}\quad \Pr(N\ge 10) = 1-F_N(9) = 0.03.\]

So, in almost one-quarter of the parallel universes that people use to interpret probabilities as frequencies there'll be seven or more positive results, but only in 3% there'll be ten or more.

This may look like cherry-picking of experiments, but it's very easy for a well-meaning team to come up with reasons why an experiment that failed to produce the desired result had a design flaw. Whereas the usual form of cherry-picking, selecting the data points within an experiment, is generally recognized by all researchers as fraud.**

But wait, there's more…

Important scientific questions like the relationship of pineapple on pizza to real cream in coffee attract large numbers of people to scientific research; so there are other teams investigating this important phenomenon. When we consider that, we can ask how likely it is for us to see a paper with 7, 8, 9, or even 10 positive results, given a total number of teams $T$ investigating this question.

The next figure shows the probability of at least one team in $T$ finding at least $N$ positive results at a 95% certainty level when each team runs 100 experiments; remember that this is given that there's no phenomenon.

With 10 teams there's a 50:50 chance of a team finding 9 or more positive results. That's usually considered a very high level of cross validation (each experiment validating and being validated by at least 8 others). And it's all based on a non-existent phenomenon.

Because scientific publications (and the mass media) tend to publish only positive results (results that show something different from what was expected), this publication bias together with the large numbers of teams $T$ and repeated experiments (here we used 100 experiments per team) can create the illusion of phenomena where none exist. Later, when these phenomena are assumed to be the default position, other researchers find them to be non-existent, and the replication crisis that we're seeing now ensues.

It's not a replication crisis, it's a ``too many phenomena were artifacts'' crisis.***

- - - -

* As this is a blog post and not the book or technical notes, here are my opinions on these weighty matters: pineapple on pizza, yes (but no pizza in the P:E diet, so it's a moot point); creamer, cream, honey, sugar, syrups, or any other adulterants in coffee, Heck no! 

** Note how the Bayesian approach deals with chains of experiments and non/low-informative results: the probability of the phenomenon given the result, $\Pr(P|R)$, given by Bayes's rule is

\[\Pr(P|R) = \frac{\Pr(R|P) \, \Pr(P)}{\Pr(R)}.\]

Bayesians can easily integrate new results with previous experiments by using the $\Pr(P)$, therefore chaining inference instead of disposing of experiments that "didn't work." And a low-informative result, where $\Pr(R|P) \simeq \Pr(R)$, i.e. the result happens with almost the same probability whether or not the phenomenon is true, will be automatically accounted as such by making $\Pr(P|R) \simeq \Pr(P)$, in other words negating its effect without disposing of its data as cherry-picking would.

*** Granted, there was some fraud as well.

A thought about the DANMASK study and presenting quantitative results

This post is about the analysis and presentation of the results, not about the substantive question of whether to wear masks. Link to the study. 

The main point here is that the way the results are presented, without a comparison counterfactual, makes it difficult to understand what the study really means:

So, without further ado, the results themselves. Cutting through a lot of important but uninteresting stuff, there are four numbers that matter:

Size of the sub-sample with masks: 2392

Number of infected in the mask sub-sample: 43

Size of the sub-sample with no masks: 2470

Number of infected in the no-mask sub-sample: 52

From these four numbers we can compute the incidence of infection given masks (1.8%) and no masks (2.1%). We can also test these numbers in a variety of ways, including using the disaggregate data to calibrate a logit model (no, I won't call it ``logistic regression''), but for now let's look at those two incidences only.

A likelihood ratio "test"

Here's a simple test we like: the likelihood ratio test between two hypotheses: that both samples are drawn from a common incidence (1.95%) or that each sample is drawn from its own incidence. In other words, we want

\[ LR = \frac{\Pr(52 \text{ pos out of } 2470|p = 0.021)\, \Pr(43 \text{ pos out of } 2392|p = 0.018)}{\Pr(52 \text{ pos out of } 2470|p = 0.0195)\, \Pr(43 \text{ pos out of } 2392|p = 0.0195)}\]

Using log-space computations to get around precision problems, we get $LR = 1.083$.

In other words, it's only 8.3% more likely that the data comes from two groups with different incidences than from a group with a common incidence. In order to be minimally convinced we'd like that likelihood ratio to be 20 or so, at least, so LR = 1.083 supports the frequentist analysis that these numbers seem to come from the same population.

(Yes, this is the same test we apply to Rotten Tomatoes ratings.)

Presenting the results: make comparisons!

A problem with the paper is that the academese of the results is hard for many people to understand. One way to make the lack of effect of masks more obvious is to create a comparison with an alternative. We choose a simple 2:1 ratio of protection, which is weak protection (a person wearing a mask has half the likelihood of infection of that of someone with no mask), but is enough to make the point.

Since we want to make fair comparison, we need to use the same size population for both the mask and no-mask conditions (we'll choose 2400 as it's in the middle of those sample sizes) and infection ratios similar to those of the test (we choose 1.25% and 2.5% for mask and no-mask, respectively). Now all we need to do is plot and compare:

(The more eagle-eyed readers will notice that these are Poisson distributions.)

The comparison with an hypothetical, even one as basic as a 2:1 protection ratio makes the point that the distributions on the left overlap a lot and therefore there's a fair chance that they come from the same population (in other words that there's no difference in incidence of infections between the use and non-use of masks).

Bayesians (and more attentive frequentists) might note at this point that having non-significant differences isn't the same thing as having a zero effect size; and that a richer model (including the distributions of the estimates themselves, which are random variables) might be useful to drive policy.

But for now, the point is that those four lines in the figure are much easier to interpret than the word-and-number salad under the subheading "results" in the paper itself.

Why is it so hard for people to change their minds?

The usual explanation is that people discount information that contradicts their beliefs. Let's build a model to explore this idea. 

We analyze the evolution of beliefs of a decision-maker receiving outside information, and to start we'll assume that the decision-maker is rational in the Bayesian sense, i.e. uses Bayes's rule with the correct conditional probabilities to update beliefs.

We call the variable of interest $X$ and it's binary: true/false, 1/0, red/blue. For the purposes of this discussion $X$ could be the existence of water on Pluto or whether the New England Patriots deflate their footballs. 

There's a true $X \in \{0,1\}$ out there, but the decision-maker doesn't know it. We will denote the probability of $X=1$ by $p$, and since there's information coming in, we index it: $p[n]$ is the probability that $X=1$ given $n$ pieces of information $X_1 \ldots X_n$.

We start with $p[0] = 1/2$, as the decision-maker doesn't know anything. A piece of information $X_i$ comes in, with reliability $r_i$, defined as $\Pr(X_i = X) = r_i$; in other words false positives and false negatives have the same probability, $(1-r_i)$.

Using Bayes's rule to update information, we have

$p[i] = \frac{r_i \, p[i-1]}{r_i \, p[i-1] + (1-r_i)(1-p[i-1])}$, if $X_i = 1$ and

$p[i] = \frac{(1-r_i) \, p[i-1]}{r_i \, (1-p[i-1]) + (1-r_i) \, p[i-1]}$, if $X_i = 0$.

For illustration, let's have the true $X=1$ (so there's indeed water in Pluto and/or the Patriots do deflate their balls), and $r_i = r$, fixed for all $i$; with these definitions, $\Pr(X_i = 1) = r$. We can now iterate $p[i]$ using some random draws for $X_i$ consistent with the $r$; here are some simulations of the path of $p[i]$, three each for $r = 0.6, 0.7, 0.8$.*

Essentially, truth wins out eventually. The more reliable the information, the faster the convergence. So, that whole "it's easier to fool someone than to get them to realize they were fooled" was wrong, wasn't it?

Only if people are Bayesian updaters with accurate perception of the reliability. In particular, when they don't let their beliefs bias that perception.

Huh-Oh!

Let us consider the case of biased perception. The simplest approach is to consider that the decision-maker's perception of reliability depends on whether the $X_i$ is in support or against current beliefs.

For simplicity the true reliability of information will still be a constant, denoted $r$; but the decision maker uses a $r_i$ that is dependent on the $X_i$ and the $p[i-1]$: if they agree (for example $p[i-1]>1/2$ and $X_i = 1$),  then $r_i = r$; if they don't (for example $p[i-1]<1/2$ and $X_i = 1$), then $r_i = (1-r)$.

Note that the $X_i$ are still generated by a process that has $\Pr(X_i = 1) = r$, but now the decision-maker's beliefs are updated using the $r_i$, which are only correct  ($r_i = r$) for draws of $X_i$ that are consistent with the beliefs $p[i-1]$, and are precisely opposite ($r_i = 1-r$) otherwise.

To illustrate this biased behavior, in the following charts we force $X_1 = 0$ (recall that $X=1$), so that the decision-maker starts with the wrong information.

There are just a few of the many paths, but they illustrate three elements that tended to be common across most simulations:

1. There's a lot more volatility in the beliefs, and much slower convergence. Sometimes, like the middle case with $r=0.8$, there's a complete flip from $p[i] \simeq 0$ to a quick convergence to the true $p[i] \simeq 1$; this was rare but worth showing one example in the image.
2. There are many cases when the decision-maker stays very close to the wrong $p[i]\simeq 0$ for very long periods (sometimes for the total length of the simulation, 1000 steps; the graphs are for the first 60 because that was enough for illustration).
3. The higher the reliability the more volatile the results can be, unlike in the case with fixed $r_i$. In general increasing reliability $r$ didn't help much with convergence or stability.

So, when people start biasing their perspectives (which might come from changing the reliability, what was simulated here, or from ignoring information that contradicts their beliefs, which is similar in effect), to counter the effect of bad early information (the $X_1 = 0$) it takes a lot of counteracting new information.

Lucky for us, people in the real world don't have these biases, they're all perfect Bayesians. Otherwise, things could get ugly.

😉

- - - - 

* In case it's not obvious, the effect of $r_i$ is symmetric around 0.5: because of how the information is integrated, the patterns for $r=0.6$ and $r=0.4$ are identical. As a consequence, when $r=0.5$ there's no learning at all and the decision-maker never moves away from $p = 0.5$.

OBSERVATION: There are also strategic and signaling reasons why people don't publicly change their minds, because that can be used against them by other players in competitive situations; but that's a more complicated — and to some extent trivial, because obvious in hindsight — situation, since it involves the incentives of many decision-makers and raises questions of mechanism design.

Of martingales and election forecasts

(This post started its life as a response to a video, but during its development I decided that there's enough negativity in the world, so it's now a stand-alone post.)

What are these martingales?

Originally a gambling strategy, martingales are discrete-time stochastic processes... hold on, I sound like the person in that video: pompous, jargon-spewing, and unhelpful.

Let's say we have some metric that evolves over time, like the advantage candidate A (for Aiden) has over candidate B (for Brenna) in an election in the fictional country of Zambonia, and that we get measures of this metric at some discrete points (every time we take a poll, for example). Note that these are a sequence of points, ordered, but not necessarily equidistant. That's what discrete-time means, that the "independent variable" (time) is ordinal but not cardinal.

(This makes a difference for many models; in actual electoral metrics it's not very important since most campaigns run daily tracking polls.)

So, we have a metric, say $A_i$, the point advantage of Aiden in poll number $i$. This is just a sequence of numbers. If they come from an underlying process which includes some unobservable or random parts we say that the $A_i$ follow a stochastic process. (Stochastic is a [insert Harvford tuition here] word for random.)

A discrete-time stochastic process is a martingale if the best estimate we have for the metric in the future is the current value, in other words,

\[ E[A_{i+1}] = A_i. \]

In some sense, we already sort-of assume that the elections are some sort of martingale: we treat the daily poll as the best estimate of the future results. Well, we used to. Some people still do, and add a lot of unsupported assumptions to develop option pricing models for... oh, bother, almost got into that negativity again.

Martingales and forecasting

A simple example of a martingale is a symmetric random walk,

\[A_{i+1} = \left\{ \begin{array}{ll}   A_i + a & \text{ with prob.  1/2} \\  A_i - a & \text{ with prob. 1/2} \end{array}\right.\]

Here are two examples, with different $a$, to show how that parameter influences the dispersion.

We can see from that figure that despite the current value being the best estimate of future values, we can make serious errors if we don't consider that dispersion. Consider the red process and note how bad the values for $A_{13}$ (POINT A) and $A_{41}$ (POINT B) are as estimates of the final value. Note also that $A_{13}$ is closer to the final value than $A_{41}$  despite $A_{41}$ being much farther along in the process (and therefore its $i=41$ is closer to the final $i=66$ than $i=13$).

Another example of a martingale is $A_{i+1} = A_i + \epsilon$ where $\epsilon$ is a Normal random variable with mean 0 and standard deviation $\sigma$. Using a standard Normal, $\sigma = 1$, here are two examples of this process:

Note how despite the same parameters and starting point, the processes' evolution is quite different. This becomes more obvious when the processes have different standard deviations:

The main point here is that even though martingales appear very simple, in that the best estimate for the future is the current value of the metric, the actual realizations of the future may be very different from the current metric.

That alone would be a good reason to try to find better ways to model elections. However this is not the only, or even the best argument against models of elections using martingales. As Ron Popeil used to say:

But wait, there's more!

The real argument here is that the process of interest (who people will vote for) and the process being measured (who the people who are willing to answer poll questions say they'll vote for) are not the same.

What's primarily wrong is that the information being used to create the $A_i$ at any point isn't an unbiased measure of the probability of Aiden winning. And that's not on the math, that's on (a) polling technique and (b) political use of polls.

Polling technique depends on people's answers, usually corrected with some measures of demographics and representativeness. For example, if Zambonia has 20% senior citizens and the polling sample only has 10%, that has to be accounted for with some statistical corrections.

Another correction comes from noticing, for example, that in previous elections the model was off by some percentage and dealing with that: if the polls for Zamboni City had Clarisse winning by 10% in the last elections but Hannibal won Zamboni City by 5%, that response bias needs to be corrected, somehow, in newer models.

Political use of polls happens when results that are known to be biased are released for political reasons. For example Aiden may release what their campaign knows to be wrong numbers to discourage Brenna donors, volunteers, and voters.

So, the problem with using martingales as a model of the election is that the information being used to generate the metrics being tracked is not an unbiased representation of the underlying reality. It's possible that the dynamics of the metric are a martingale, but what the metric is measuring is not the electoral vote but a mix of socially acceptable answers (who wants to say they're voting Hannibal rather than Clarisse, even when they are?) and push-poll results designed to influence the electoral process

Many professional political forecasters deal with this mismatch using field-specific knowledge and heuristics. Certain others criticize them for the heuristics and field-specific knowledge while missing the problems implicit in using martingale-based models.

Good, no Taleb references at all. 🤓

Recommendation: readers interested in political (and other) forecasting might want to read Superforecasting, by Phil Tetlock and Dan Gardner.

Fun with COVID-19 Numbers for March 1, 2020

NOTA BENE: The Coronavirus COVID-2019 is a serious matter and we should be taking all reasonable precautions to minimize contagion and stay healthy. But there's a lot of bad quantitative thinking that's muddling the issue, so I'm collecting some of it here.

Death Rate I: We can't tell, there's no good data yet.

This was inspired by a tweet by Ted Naiman, MD, whose Protein-to-Energy ratio analysis of food I credit for at least half of my weight loss (the other half I credit P. D. Mangan, for the clearest argument for intermittent fasting, which convinced me); so this is not about Dr Naiman's tweet, just that his was the tweet I saw with a variation of this proposition:

"COVID-19 is 'like the flu,' except the death rate is 30 to 50 times higher."

But here's the problem with that proposition: we don't have reliable data to determine that. Here are two simple arguments that cast some doubt on the proposition:

⬆︎ How the death rate could be higher: government officials and health organizations under-report the number of deaths in order to contain panic or to minimize criticism of government and health organizations; also possible that some deaths from COVID-19 are attributed to conditions that were aggravated by COVID-19, for example being reported as deaths from pneumonia.

⬇︎ How the death rate could be lower: people with mild cases of COVID-19 don't report them and treat themselves with over-the-counter medication (to avoid getting taken into forced quarantine, for example), hence there's a bias in the cases known to the health organizations, towards more serious cases, which are more likely to die.

How much we believe the first argument applies depends on how much we trust the institutions of the countries reporting, and... you can draw your own conclusions!

To illustrate the second argument, consider the incentives of someone with flu-like symptoms and let's rate their seriousness or aversiveness, $a$, as a continuous variable ranging from zero (no symptoms) to infinity (death). We'll assume that the distribution of $a$ is an exponential, to capture thin tails, and to be simple let's make its parameter $\lambda =1$.

Each sick patient will have to decide whether to seek treatment other than over-the-counter medicine, but depending on the health system that might come with a cost (being quarantined at home, being quarantined in "sick wards," for example); let's call that cost, in the same scale of aversiveness, $c$.

What we care about is how the average aversiveness that is reported changes with $c$. Note that if everyone reported their $a$, that average would be $1/\lambda = 1$, but what we observe is a self-selected subset, so we need $E[a | a > c]$, which we can compute easily, given the exponential distribution, as

\[
E[a | a > c]
=
\frac{\int_{c}^{\infty} a \, f_A(a) da }{1 - F_A(c)}
=
\frac{\left[ - \exp(-a)(a+1)\right]^{\infty}_{c}}{\exp(-c)}
= c + 1
\]
Note that the probability of being reported is $\Pr(a>c) = \exp(-c)$, so as the cost of reporting goes up, a vanishingly small percentage of cases are reported, but their severity increases [linearly, but that's an artifact of the simple exponential] with the cost. That's the self-selection bias in the second argument above.

A plot for $c$ between zero (everyone reports their problems) and 5 (the cost of reporting is so high that only the sickest 0.67% risk reporting their symptoms to the authorities):

Remember that for all cases in this plot the average aversiveness/seriousness doesn't change: it's fixed at 1, and everyone has the disease, with around 63% of the population having less than the average aversiveness/seriousness. But, if the cost of reporting is, for example, equal to twice the aversiveness of the average (in other words, people dislike being put in involuntary quarantine twice as much as they dislike the symptoms of the average seriousness of the disease), only the sickest 13.5% of people will look for help from the authorities/health organizations, who will report a seriousness of 3 (three times the average seriousness of the disease in the general population).*

With mixed incentives for all parties involved, it's difficult to trust the current reported numbers.

Death Rate II: Using the data from the Diamond Princess cruise ship.

A second endemic problem is arguing about small differences in the death rate, based on small data sets. Many of these differences are indistinguishable statistically, and to be nice to all flavors of statistical testing we're going to compute likelihood ratios, not rely on simple point estimate tests.

The Diamond Princess cruise ship is as close as one gets to a laboratory experiment in COVID-19, but there's a small numbers problem. In other words we'll get good estimates when we have large scale, high-quality data. Thanks to @Clarksterh on Twitter for the idea.

Using data from Wikipedia for Feb 20, there were 634 confirmed infections (328 asymptomatic) aboard the Diamond Princess and as of Feb 28 there were 6 deaths among those infections. The death rate is 6/634 = 0.0095.

(The ship's population isn't representative of the general population, being older and richer, but that's not what's at stake here. This is about fixating on the point estimates and small differences thereof. There's also a delay between the diagnosis and the death, so these numbers might be off by a factor of two or three.)

What we're doing now: using $d$ as the death rate, $d = 0.0095$ is the maximum likelihood estimate, so it will give the highest probability for the data, $\Pr(\text{6 dead out of 634} | d = 0.0095)$. Below, we calculate and plot the likelihood ratio between that probability and the computed probability of the data for other candidate death rates, $d_i$.**

\[LR(d_i) = \frac{\Pr(\text{6 dead out of 634} | d = 0.0095)}{\Pr(\text{6 dead out of 634} | d = d_i)}\]

We can't reject any rates between 0.5% and 1.5% with any confidence (okay, some people using single-sided point tests with marginal significance might narrow that a bit, but let's not rehash old fights here), and that's a three-fold range. And there are still a lot of issues with the data.

On the other hand...

It's easy to see that the COVID-19 death rate is much higher than that of the seasonal flu (0.1%): using the data from the Diamond Princess, the $LR(0.001) =  3434.22$, which should satisfy both the most strong-headed frequentists and Bayesians that these two rates are different. Note that $LR(0.03) = 510.01$, which also shows that with the data above the Diamond Princess invalidates the 3% death rate. (Again, noting that the numbers might be off by a factor of two or three in either direction due to the delay in diagnosing the infection and between diagnosis and recovery or death.)

As with most of these analyses, disaggregate clinical data will be necessary to establish these rates, which we're estimating from much less reliable [aggregate] epidemiological data.

Stay safe: wash hands, don't touch your face, avoid unnecessary contact with other people. 

- - - - - 

* A friend pointed out that there are some countries or subcultures where hypochondria is endemic and that would lead to underestimation of the seriousness of the disease; this model ignores that, but anecdotally I've met people who get doctor's appointments because they have DOMS and want the doctor to reassure them that it's normal, prescribe painkillers and anti-inflammatories, and other borderline psychotic behavior...

** We're just computing the binomial here, no assumptions beyond that:

$\Pr(\text{6 dead out of 634} | d = d_i) = C(634,6) \, d_i^6 (1-d_i)^{628}$,

and since we use a ratio the big annoying combinatorials cancel out.