Coronavirus Modeling

The Metric We Need to Manage COVID-19

It’s Easter. A couple weeks ago, this was the target for returning to normal. Hospitals are now full of patients, cities bloom as new hotspots, and politicians wrestle with the balance of human and economic costs. We’re left to wonder if we are equipped with the right metrics to guide our path forward. Add to the confusion that metrics are based on noisy data that changes daily. There’s one metric, however, that has the most promise. It’s called Rt – the effective reproduction number. We can estimate it, and it’s the key to getting us through the next few months.

Most people are more familiar with R0. R0 is the basic reproduction number of an epidemic. It’s defined as the number of secondary infections produced by a single infection. If R0 is greater than one, the epidemic spreads quickly **. If R0 is less than one, the epidemic spreads, but limps along and disappears before everyone becomes infected. The flu has an R0 between one and two while measles sits in the high teens. While R0 is a useful measure, it is flawed in an important way: it’s static.

We’ve all witnessed that humans are adaptable. Our behavior changes, whether mandated or self-prescribed, and that changes the effective R value at any point in time. As we socially distance and isolate, R plummets. Because the value changes so rapidly, Epidemiologists have argued that the only true way to combat COVID19 is to understand and manage by Rt.

I agree, and I’d go further: we not only need to know Rt, we need to know local Rt. New York’s epidemic is vastly different than California’s and using a single number to describe them both is not useful. Knowing the local Rt allows us to manage the pandemic effectively.

States have had a variety of lockdown strategies, but there’s very little understanding of which have worked and which need to go further. Some states like California have been locked down for weeks, while others like Iowa and Nebraska continue to balk at taking action as cases rise. Being able to compare local Rts between different areas and/or watch how Rt changes in one place can help us measure how effective local policies are at slowing the spread of the virus.

Tracking Rt also lets us know when we might loosen restrictions. Any suggestion that we loosen restrictions when Rt > 1.0 is an explicit decision to let the virus proliferate. At the same time, if we are able to reduce Rt to below 1.0, and we can reduce the number of cases overall, the virus becomes manageable. Life can begin to return to ‘normal.’ But without knowing Rt we are simply flying blind.

How We Can Calculate Rₜ

It’s impossible to measure Rt directly, so we have to estimate it. Fortunately, there are many ways to this. One particular method Bettencourt & Ribeiro described in their 2008 paper, “Real Time Bayesian Estimation of the Epidemic Potential of Emerging Infectious Diseases.” This solution caught my attention because it focuses on the same principles from my first post, Predicting Coronavirus Cases. It uses Bayesian statistics to estimate the most likely value of Rt and also return a credible interval for the true value of Rt.

What follows is an application of Bettencourt & Ribeiro’s process (with an important modification) to US State COVID19 data. Note that while this post focuses on the high level concepts, those who want to dig in further can find the details in this Jupyter notebook.

Bettencourt & Ribeiro’s original algorithm to estimate Rt is a function of how many new cases appear each day. The relationship between the number of cases yesterday and the number of cases today give us a hint of what Rt might be. However, we can’t rely on any one day too much in trying to guess Rt, as daily case counts are imperfect due to changing testing capacity, lags in data reporting, and random chance. However, using Bayes’ Theorem, we can take the new information we get from each day’s case count to adjust our expectation of what Rt is, getting closer to the true value as more daily data becomes available.

I applied this algorithm to the data to produce a model for each state’s Rt, and how it changes over time. But I noticed something strange. Over time, all states trended asymptotically to Rt = 1.0, refusing to descend below that value. Somehow, the algorithm wasn’t reflecting the reality that Rt could be < 1.0 as well.

To fix this, I made one significant change to their algorithm that maintains the integrity of Bettencourt & Ribeiro’s original work while allowing us to see the real-time picture clearly. The change was simple. Instead of considering every previous day of data we have to estimate Rt, I only use the last seven days. Doing so is mathematically sound and produces more accurate results when the model is compared to actual data, but I admit is not reviewed by anyone. While I invite feedback, I’m sharing these results with that disclaimer well in advance.

The Results

The algorithm produces a most-likely value for Rt over time for each locale. Below you can see select results for Michigan, Washington and New York.

The graphs above show two important things. First, the ‘most likely value’ of Rt for each day, represented by the dots (the more red a dot is, the higher and more dangerous the value of Rt is). Second, the interval of values that Rt might actually be is represented by the gray bands (if there’s a ‘most likely value’, then there are other less likely but still possible values as well).

Because we’re using Bayesian methods, the model produces something called a highest-density interval (HDI). In this case, the model claims that 95% of the possible true values lie in the gray band, that is, in the HDI. This is helpful when looking at states where conclusions might not be as certain. Consider Maine.

Comparing Maine to the previous states is instructive. Maine looks to have reduced Rt below 1.0, but the gray band shows us that we cannot safely conclude we are below the safety threshold, since it’s possible that Rt is still above 1.0. In this case, restrictions should continue. The reason Maine’s HDI is wide has to do with the lack of information. Maine has fewer cases than the previous examples, so the margins are wider. The more data we get and the more consistent that data is, the more confident we can be in our conclusions. The explicit HDI helps make sure we’re not misled by the data and prematurely conclude something that might be a grave error. This is useful when it comes to policy decisions: wait until the high end of gray band is well below 1.0 before resuming normal activity.

You can review each state in the main graphic at the top, but a denser view might be helpful, too.

How do States Compare?

In this image, I’ve taken the most recent values for Rt and plotted them along with the highest-density interval (HDI) bands (as of the evening of 4/11). Doing this allows us to see all sorts of interesting things.

Notice how many of the states without a shelter-in-place order—Iowa, Nebraska and South Dakota— have some of the highest Rt values. These states should be aware that their high Rt will lead to exponential growth in cases. As I argued in my previous post, when this happens, we need to lock down, so why not lock down now before cases grow? Doing so will save countless lives.

Now, you might point out that North Dakota looks well below 1.0 and does not have a lockdown. But here it’s important to consider North Dakota’s HDI. North Dakota’s interval is so wide that it includes many possible values above 1.0—therefore we can’t safely conclude that Rt is truly below 1.0. If the actual value lies in the highest portion of the HDI, North Dakota might have one of the highest Rt values. In fact, if we plot the same data, but sort by the high end of the HDI (worst-case scenario) we see this:

Now it’s clear that all non-lockdown states cluster on the right-hand side. While we cannot be sure the true Rt value is that high, this graph should give any policymaker in a non-lockdown state pause.

Instead of worst-case scenarios, let’s look at which states are almost certainly below the 1.0 threshold. These are the states that are clearly doing well.

Here, I’ve plotted states where the high end of their HDI is less than 1.1. This means that even in the worst case scenario, the state likely has the epidemic “under control.” You can see that many of the states in the headlines for the past few weeks are in this group. Surprisingly, Louisiana seems to have turned the corner decisively, even with one of the worst outbreaks. Note, however, that “under control” does not mean good. Louisiana still has far more new cases than Washington per day, but Louisiana’s outbreak will not be growing as nearly quickly as Washington and will soon see cases decline even more rapidly. Whatever they’re doing is working. But even in these states—the most “under control” in the nation— Rt is still far from the comfortable value of ~0.5. Even Idaho has a wide HDI, so it’s too early to tell.

Which states have the epidemic least “under control?” To answer this, I plotted states where the best case (eg. low end of HDI) is above 1.0, indicating the true value of Rt is almost certainly above 1.0. Surprisingly, Rhode Island, Maryland and Massachusetts sit at the top. Part of this may be that Rhode Island is earlier in its infection curve, but seeing large states like Massachusetts and Texas above 1.0 is worrisome —especially because none of these states have hit the headlines as being trouble spots.

While the measure of Rt can be abstract, it shows up in the case data. Compare Massachusetts’s new cases per day to Louisiana’s. Louisiana’s cases per day peaked a week ago while Massachusetts’s continue to ascend. Remember that we’re looking at new cases per day, so in both, total cases continue to climb but Louisiana’s rate of increase is declining rapidly. Rt is a measure of direction, not absolute size.

Looking Forward

Without the use of a clear metric on our ability to contain the coronavirus pandemic, it’s difficult to imagine that we’ll manage a return to normalcy anytime soon. I hope by sharing this work with you, you will consider Rt as the metric that can guide our analysis and decision-making. I’ve kept most of the math and theory in the notebook, but I’d highly suggest diving in if you’re so inclined. And as always, you can get in touch by email or reaching out on Twitter or Instagram @kevin

** Note: In a previous version of this post I wrote that R0 > 1 meant that the virus spread to everyone in the population. This is not true. While R0 values in the range of 3-4 as shown above mean it will spread to nearly all people in the population, a value of 1.5 might only spread to ~60% (depending on your parameter estimates). In any case, our goal should be to reduce R0 to as small a value as possible. Thank you to Jim Fisher for pointing this out.