Thread by @gro_tsen

There's something really fishy going on with “effective reproduction numbers” in various analyses of this pandemic. And I'm not talking about predicting the future, I'm talking about analysing the past. Here we see a claim of R~3 for Belgium. •1/21

https://twitter.com/vdwnico/status/1320267782622162945

On the other hand, researchers at the Centre for Mathematical Modelling of Infectious Diseases have a web site computing R(t) for many countries, and their estimates for Belgium never went above 1.5 in the computed time frame: epiforecasts.io/covid/posts/na… •2/21

Now briefly speaking, computing R(t) depends on two things: the exponential rate of growth r (=logarithmic slope, =logarithmic derivative) of the number of new cases, and temporal data on infections, notably the “serial interval” and its distribution. •3/21

If the serial interval d is constant (each infected person infects R new always exactly d days after infection), we can compute R very simply as the number of new cases at date D divided by the same number at date D−d. This is because R=exp(r·d). Simple. BUT! •4/21

First, the data aren't perfect (e.g., there's a strong weekly bias in the number of cases, and reporting isn't instantaneous); second, the timing is variable (btw infection and spreading, btw infection and test, and btw test and reporting): serial intervals vary. •5/21

AFAIU, many people, including France's own health agency, simply compute a “reproduction number” by the simplistic formula above with d=7, i.e., their “R” is: number of cases at date D ÷ number of cases at date D−7 (both numbers averaged over 7 days, I guess?). •6/21

This kind of “reproduction number” is just a rate of exponential growth converted by assuming a constant serial interval. It's useful for simply viewing whether cases are on an upward or downward trend, and how fast, comparing places; but it has many problems. •7/21

The first is that the value will be very noisy. Maybe the real R is also very noisy, but its short-term variations probably bear absolutely no relation to the rapid variations of the value thusly computed. See

https://twitter.com/covid_stat/status/1320072079291654146

for France to get a feel of wobbles: … •8/21

… I'm pretty sure all these rapid variations are just pure noise, due to bad data to start with, and difficulty in removing underlying noise and random timing variations. •9/21

The second problem is that, if the serial interval isn't constant (people don't take a constant time d to infect others), the relation between R and the exponential slope r of infection becomes more delicate than simply R = exp(r·d) … •10/21

… which underlies the computation of R for constant d as a simple ratio of new cases d days apart. Instead, we must average over the distribution of the serial interval. For example, imagine ½ of infections occur after d₁ days and ½ after d₂ (with d₂>d₁ say): … •11/21

… then R will be 1/[½(exp(−r·d₁)+exp(−r·d₂))] (compare Euler-Lotka equation), so even in a constant and steady exponential growth, the d we should use to compute R is −log(½(exp(−r·d₁)+exp(−r·d₂)))/r, which depends on r: … •12/21

… essentially, when r (or R) is large, fast infections (d₁) dominate slow infections (d₂) whereas if r (or R) is small, it's the other way around. So we can't simply use some “average” d value and divide values d days apart, … •13/21

… because doing so will overestimate R when it is large, or underestimate when R is small (or both, depending on what kind of average d has been chosen to be). So in effect, it artificially amplifies the variations of R, even with perfect data. •14/21

Now the “sophisticated” way to compute R(t) is more delicate: we use an a priori distribution on it, hypotheses on the distribution of the serial interval (not assumed constant, that's the point), predict values and condition by the observations made, … •15/21

… and then we get an a posteriori distribution on R, of which we generally retain the maximal likelihood value (plus uncertainty bounds). This is well explained here (in French):

https://twitter.com/gro_tsen/status/1319394804443324418

•16/21

Of course the real nitty-gritty details are, well, nitty and gritty. There are subtly different ways to do the computation, implemented in various packages for the (abominable) “R” stats program, such as EpiEstim cran.r-project.org/web/packages/E…‌; and they are difficult to use. •17/21

(I thank @jeuasommenulle for giving me a start on how to actually use these tools, which I still need to investigate more fully, both in theory and in actual practice.) And of course we need hypotheses on the distribution of the serial interval, which isn't well known! •18/21

All this is to say, computing R(t) is a more delicate matter than dividing two numbers, even though many people are content with this. How it's done doesn't matter to extrapolate a simple exponential growth, of course, but it matters if we try to look further than that. •19/21

So, bottom line, if you see an estimate of an “effective reproduction number” you should question how it was computed: is it a simple quotient or a maximum likelihood? And crucially: what hypotheses have been made on the serial interval? •20/21

When in doubt, I would tend to consider as more accurate the values computed by the CMMID at epiforecasts.io/covid/ — their inputs won't be magically better, uncertainties are huge, and their graphs aren't updated as often as one might wish, but I trust their expertise. •21/21

PS1: Correction: it seems that the French health agency doesn't use a simple quotient but a more sophisticated method, running the computation under Excel (😭). Details, however, are still very unclear on what exactly they do. •22/(21+2)

https://twitter.com/gro_tsen/status/1320334918782734336

PS2: I also didn't discuss the issue that R(t) can be computed from various input sources (lab positives, symptomatic cases, hospital admissions, serious cases, or even deaths), all of which will have their own specific issues. •23/(21+2)

https://twitter.com/gro_tsen/status/1320336012757901313