Two Wrongs

Forecasting Covid-19 Variants

Forecasting Covid-19 Variants

A recent question on Metaculus asks whether the who will name BA.2.86 as a variant of interest before October 1, 2023.

At the time of writing, the community has this question at 35 %. Apparently Ars Technica reports that it has

34 mutations in its spike gene relative to BA.2, the omicron sublineage from which it descended. This number of spike mutations between BA.2.86 and BA.2 is chillingly similar to the number of mutations seen between the original omicron (BA.1) and the ancestral Wuhan strain.

Does this mean 35 % is appropriate? Heck if I know! I’m no epidemiologist. I don’t know how likely it is that a variant is upgraded from vum (variant under monitoring) to voi (variant of interest) when its “spike mutations are chillingly similar to the number of mutations seen between the ancestral Wuhan strain and the oroginal omicron.”

Hunting for base rate data

Let’s back up and take the outside view.

I dug up a somewhat sloppily produced pdf published by the who that purports to be a changelog of variant classifications. The list seems incomplete, because there are variants that change classification before being assigned a classification, and some variants are assigned a classification but then don’t get mentioned in the rest of the document, which spans more than two years.

Here is a brief summary of the relevant variant changes listed in the document.

  • There are 28 variants that entered the list as variants under monitoring.
  • Of these, 4 were upgraded to variants of interest.
  • Another 20 were de-escalated, meaning they ceased being under monitoring without entering a higher severity class.

We can sense some complications in this, like how the fate of four variants is unknown. In survival analysis, we say that the data is censored.1 Additionally, there is one variant that was upgraded to voi from vum without first entering as vum, and four variants that were de-escalated from vum without entering as vum. This is known as truncation. We could use tools from survival analysis to deal with this, but we will – at least for now – focus on the specific question asked on Metaculus, assume the effects of censoring and truncation are relatively small and look just at the cases for which we know the full timeline.

  • 24 variants enter under monitoring, of which
  • 4 are eventually upgraded to variants of interest.

Making sense of data

We can already provide a very naïve estimate for the base rate here. When a variant enters under monitoring, we have – thanks to Laplace’s law of succession – a posterior of

\[\frac{4+1}{24+2} \approx 19 \%\]

probability that it eventually becomes a variant of interest. Already, this is lower than the community 35 %.

But it doesn’t account for the time horizon of the question! Some of those 19 % of variants will be upgraded to a variant of interest only after the question time horizon is over, which means they don’t count for our forecast. What’s the timeline of upgrades? The variants under monitoring that became variants of interest stayed under monitoring for

7 days 22 days 49 days 70 days

This means that even assuming that a variant becomes a variant of interest eventually, there’s only about a 70 % chance that it happens within 44 days of listing.2 Linear interpolation between 50 % at 22 days and 75 % at 49 days. You might find it odd that I dare linearly interpolate an empirical distribution of just four values, but as we’ll see later, the distribution of monitoring times also for de-escalated variants is fairly linear, which lends credibility to the method.

So in total, the probability that a variant is upgraded to a variant of concern within 44 days of being listed as a variant of concern is 19 % × 70 % = 14 %. This probably suffices as a first forecast, but it doesn’t tell us anything about how to update our forecast as time passes.

De-escalated variants linger longer

Note that the variants that were upgraded were under monitoring for an average of 37 days. The variants that are under monitoring but get de-escalated stay under monitoring for an average of 119 days! For every day that passes and a variant is not upgraded, that should slightly increase our belief that this variant belongs to the group that will eventually be de-escalated.

We can compute this using Bayes’ rule in its odds form:

\[o(u \mid T > t) = \frac{P(T > t \mid u)}{P(T > t \mid \neg u)} o(u)\]

where \(u\) is the hypothesis “BA.2.86 will eventually be upgraded” and \(T > t\) is the observation “BA.2.86 has been under monitoring for at least \(t\) days”.

It helps here to know that

\[o(E) = \frac{P(E)}{1 - P(E)}\]

and conversely

\[P(E) = \frac{o(E)}{1+o(E)}\]

Applying this, we’ll learn the odds that BA.2.86 will eventually be upgraded given that it has been under monitoring for \(t\) days, and then we will recompute this for each day that passes.

We already know our prior odds \(o(u)\) – that’s the 19 % probability we computed before, an odds of 0.23.

Survival function gives us the odds ratio

To get the posterior out of Bayes’ rule, we still have to compute the odds ratio, which had the numerator and denominator

  • \(P(T > t \mid u)\), the probability that a variant that will eventually be ugraded sticks around under monitoring for at least \(t\) days; and
  • \(P(T > t \mid \neg u)\), the probability that a variant that will eventually be de-escalated sticks around under monitoring for at least \(t\) days.

Fortunately for us, we don’t have to compute these for more than 44 days out, which means we can just draw directly from the empirical distribution. To re-iterate, this was, for variants that are upgraded,

7 days 22 days 49 days 70 days

and for variants that are de-escalated (also in days),

27 55 56 63 63 63 84 91 91 105
120 120 133 140 141 175 175 203 231 238

These are two empirical survival functions, i.e. they indicate the probability that a variant has survived on the list of variants under monitoring for so many days.

who-vum-01.svg

These two curves express exactly \(P(T > t \mid H)\) where \(H\) is \(u\) for the red curve, and \(\neg u\) for the blue one. If we divide the lower curve with the upper curve for the first 44 days, we get the odds ratio we’re looking for. I did it in a spreadsheet with linear interpolation:

Day P(T > t ⎸ u) P(T > t ⎸ ~u) Odds ratio
0 1 1 1
5 0.82 0.99 0.82
10 0.7 0.98 0.71
15 0.61 0.97 0.63
20 0.53 0.96 0.55
25 0.47 0.95 0.49
30 0.42 0.94 0.44
35 0.37 0.94 0.40
40 0.33 0.93 0.35
44 0.29 0.93 0.31

The last column tells us that, for example, on day 5, the odds that we’re dealing with a (future) variant of interest is already 0.82 times the prior odds, which was 0.23. If we multiply this into the odds ratio, and convert to probabilities, we get our posterior probability that, for \(t\) days passed, we are looking at a variant that will eventually become upgraded. Presented in graphical form:

who-vum-02.svg

Probability of upgrade over time

Now we know how our probability that BA.2.86 eventually becomes a variant of interest changes as time passes. The remaining piece of the puzzle is how the conditional probability, i.e. if we assume that BA.2.86 will be upgraded, that it happens within the first 44 days, changes over time.

Mathematically, we can express this as

\[P(T \le 44 \mid u, T > t)\]

If we think of the probability that a variant under monitoring is upgraded to variant of interest as mass distributed over the real line between 0 and 70, then the conditional probability above expresses “the proportion of the remaining area that comes before 44.”

The survival function \(P(T > t \mid u)\) we we’ve seen already represents the full distribution, so we can subtract values from it cleverly to get that relative area of probability mass. Specifically,

\[P(T \le 44 \mid u, T > t) = \frac{P(T \le 44 \mid u) - P(T \le t \mid u)}{1 - P(T \le t \mid u )}\]

where the survival function

\[P(T > t \mid u) = 1 - P(T \le t \mid u)\]

This will be a decreasing function – first slowly, then more quickly, as the opportunities for upgrade fade away until on day 44 on the list there are no more opportunities for upgrading before the time horizon of the question has passed.

Multiply together for final prediction

The question asks for \(P(T \le 44, u | T > t)\) – what is the probability that this variant will stop being under monitoring within 44 days, and that it is then upgraded to variant of interest, given that it has been \(t\) days under monitoring – and we have implicitly decomposed this as

\[P(T \le 44, u | T > t) = P(T \le 44 | u, T > t) P(u | T > t)\]

We have further used Bayes’ rule to determine the second factor, namely

\[o(u \mid T > t) = \frac{P(T > t \mid u)}{P(T > t \mid \neg u)} o(u)\]

and we just figured out how to evaluate \(P(T \le 44 | u, T > t)\). What remains is to multiply everything together! The result of that multiplication is our forecast, presented graphically:

who-vum-03.svg

As I’m writing this, 8 days have passed since BA.2.86 was first entered as a variant under monitoring, which would put the forecast at 8.7 % – certainly a far way off from the community 35 %. Remains to be seen whether I’ve missed something the community has not, or if the high community forecast merely reflects uncertainty in the subject area.

Appendix A: Tabulation of forecast values

Tabulation of forecast for my own convenience.

Days under monitoring Probability of upgrade within 44 days
0 13.1 %
1 12.5 %
2 12.0 %
3 11.4 %
4 10.8 %
5 10.2 %
6 9.6 %
7 8.9 %
8 8.7 %
9 8.4 %
10 8.1 %
11 7.8 %
12 7.5 %
13 7.2 %
14 6.9 %
15 6.6 %
16 6.3 %
17 5.9 %
18 5.6 %
19 5.3 %
20 5.0 %
21 4.6 %
22 4.3 %
23 4.1 %
24 3.9 %
25 3.8 %
26 3.6 %
27 3.4 %
28 3.2 %
29 3.0 %
30 2.8 %
31 2.6 %
32 2.4 %
33 2.2 %
34 2.0 %
35 1.8 %
36 1.6 %
37 1.4 %
38 1.2 %
39 1.0 %
40 0.8 %
41 0.6 %
42 0.4 %
43 0.2 %
44 0.0 %