Thursday, April 09, 2020

Updated method and forecasts

I've tweaked the method just a little since my earlier forecasts. 

You may have noticed that the initialisation of my short-term forecasts was fairly poor. This is an inevitable consequence of naively fitting parameters in a simple imperfect deterministic model using a long time series. It gives good overall parameter estimates but the inevitable limitations of the model mean that the final state (from which we look ahead in a forecast) isn't truly a realistic estimate. This can easily be seen in my very first forecast for the UK for example. The spread of model trajectories is quite wide and biased high compared to the last few data points:

A consequence is that even if the trend is perfect, the one-day forecast will be quite poor. We weren't really going to jump from the observed 621 to 1900 deaths in a single day for example, but the day-ahead forecast spread reached that high a value. Various fudges can be used to address this problem, using anomalies is one but that's problematic during exponential growth and won't satisfy the nonlinear model dynamics. The route I have chosen is to weight the more recent data more highly in the fit. The justification for this is that it is a reasonable approximation to allowing the model to forget older data as would happen automatically if I used a stochastic model with sequential updating, or otherwise represented dynamical model error in a more explicit manner. It clearly works at a heuristic level and I'm working on how best to calibrate the effect and implement it in a more technically justifiable manner. I suppose I may end up publishing this somewhere, so I ought to make sure it is defensible.

Last week I looked at Lombardy, and argued that the lockdown was probably working. Here is a re-run of my forecast, withholding the last 7 data points, using a version of the updated method.


You see those cyan flashing circles? They are the 7 latest data points, which were not used in any way in this estimation. Promise!

For giggles, have a look at what IC forecast for Italy last week, and how well that turned out. Each of their forecasts run from Monday to Sunday, the last 4 weeks of experimental forecasts have obs plotted alongside (under!) them and the pink plume is what they forecast just a few days ago for this week. That last one might be just about ok.


 Here is what my method predicted for all of Italy starting from the same Sunday as their most recent failed forecast, again with validation data plotted:


I must admit I am impressed as to how well it has worked, though it should be noted that at the heart of it, it's basically a 3 parameter fit for a model that generates two straight lines (in log-space) with a bend joining them so not a terribly difficult problem. If you know the shape of the bend and how it is formulated, there's a good chance you can predict roughly where and when it's going to stop bending. The remarkable thing is that the experts at IC haven't managed to achieve this. I find it baffling that they are quite so hopeless. There are about 50 or 60 of them co-authoring that report, and they do it for a living. 

Here is another example, my forecast for Spain, initialised a week last Sunday, versus the IC group forecast at the same date (the rightmost greenish plume with dots below was initialised at the same date).





And here is the current UK forecast...before today's figure comes out.

This is the IC forecast for the UK for this week again (pink plume again, below). The data were already outside their range by yesterday. What on earth were they thinking?


Wednesday, April 08, 2020

Dumb and dumber

All these people exhorting amateurs to "stay in their lane" and not muddy the waters by providing analyses and articles about the COVID-19 pandemic would have an easier job of it if it wasn't for the supposed experts churning out dross on an industrial scale.

Two new forecasts for the UK caught my eye yesterday. First from IHME, it predicts a rising death toll to a daily value of over 200 by the end of the week, peaking at about 3000 later in the month (with much uncertainty). It seems to be some sort of curve-fitting exercise. They are regarded as "world-leading disease data experts" by the Guardian at least. Here is their forecast for the UK as reported in the Guardian:


and again from their own report with uncertainties (lines are at 1k intervals). So you can see they expect a strong rise and guarantee that it will go over 1k deaths later in the month.


The article describing their method is here, it's some sort of fancy curve fitting that doesn't seem to make much use of what is known about disease dynamics. I may be misrepresenting them somewhat but we'll see below what a simple disease model predicts.

Next, the famous IC group. And these are the real ones with Ferguson as last author, not the jokers who generated the silly 5700 number that got a brief bit of attention. This is the MRC-funded, WHO-collaborating lot. And what have they done here? They've also basically fitted a curve through historical death data,  and extrapolated it into the future. Their prediction is for a whopping 14k deaths with a range of 8-16k - just in the next week! A minimum of 1500 on the 7th day in that forecast too.




So why do I think these are bonkers? Because we KNOW that a lot changed two weeks ago. And we KNOW that this feeds though into deaths in about 3 weeks. Typical time from infection to death may be a little over 3 weeks but a bunch of people die significantly quicker and this means the effect of the social distancing controls should be feeding though into death rate, oh, about now.

This below is what I get from fitting my simple SEIR model to the available death data, backdated to Sunday to form a fair comparison with the IC group. 

Now, before looking at it too carefully, I do need to give a "health warning" on its validity. In particular, I don't expect the uncertainty range to be very well calibrated. I'm really focussed mainly on the best fit, and the parametric uncertainty included here is a rather subjective component. I suspect the uncertainty intervals to be rather too broad, especially at the top end (which essentially allows for continued strong growth of the epidemic, contrary to what has been seen in other countries that clamped down sooner). The corollary being that I'm pretty confident of reality lying inside the bounds, and optimistic of my median being a much better estimate of the truth than the central forecasts provided by the experts above :-)

The forecast is necessarily vague at this point, because although I believe R changed at the breakpoint, I don't yet know what it changed to. "About 1" is a plausible estimate, it might be higher but we can certainly hope that it's lower as it proved to be in Wuhan and also Italy after they clamped down.





My model predicts a total of 8k deaths next week, with a 5-95% range of 4-19k. Yes it's a wide uncertainty range, I think my prior on Rt is probably still too broad as I don't really expect to see a value lower than 0.5 or higher than 1.5 (and these are just the 1sd spread limits in the above). But I am very optimistic that the median estimate generated by this method is better than the experts have provided, and they don't seem to believe that anything in the lower half of my range is possible at all. Ominously for them, the recent data are already consistently lower than the median, albeit marginally.

I'll come back to these in a week to see how we fared.

Sunday, April 05, 2020

A time to die?

Edit: To save the reader from having to plough through to the bottom, it is apparent that people in some countries are only being tested once seriously ill and indeed close to death. Hence, the time from "case" to "death" is very small, even though the time from infected/infectious may be much longer.
---

A cheery title for a cheery topic.

One issue that causes significant difficulties in assessing the epidemic is not knowing the fatality rate or case numbers with any real certainty. We do at least have a ballpark estimate for fatality rate of around 1%, probably higher if you only count symptomatic cases, but maybe a bit lower if you assume a significant number of asymptomatic ones too. I've generally been assuming 0.75% for the total fatality rate, based on 1.5% for symptomatic cases and an equal number asymptomatic. This number doesn't really affect many of my shorter-term calculations anyway though it matters for the total death toll under full epidemics.

However there is another issue that is also extremely important in understanding and interpreting observations, which that has not received anywhere so much attention. This is the length of time that it takes people who die, to die. Since death data are the only semi-reliable numbers that we have (not perfect, but better than anything else in the UK for sure), even if we assume the death rate is known, and so we know how many cases there must have been at some point in the past, if we don't know when these cases occurred, we are little the wiser. Given a 3 day doubling time, 100,000 cases 5 days ago and 100,000 cases 15 days ago mean very different things for today!

The time-to-die calculation in the IC work seems to be very strongly reliant on a very small data set from early in the Wuhan outbreak. The analysis of these data is presented in this paper here. I think it's a nice piece of work that does the best job available with the limited data. But it only considers 24 deaths, to which it fits a parametric curve. Here is the data set, and the fitted distribution.





The red curve is a direct fit to the data, the blue one is after adjusting for the epidemic growth (which will have biased the observed sample towards shorter times, simply because when looking at deaths on a given day, there were more new cases 10 days ago than there were 20 days ago etc). Ticks are every 10 days, ie the longest observed time to death is 21 days and there isn't a single death in the upper 30% end of the estimated (blue) distribution which falls above this value. This long tail has been created entirely by a parametric extrapolation. I have no idea what theoretical basis there may be for using a Gamma distribution as they did. It may just be convenience - I have no problem with that in principle, but it does mean that this result is rather speculative. 10% of the deaths are estimated to take fully 30 days from symptoms appearing.

During the exponentially-rising part of the epidemic, this time to death can't be easily deconvolved from the death rate anyway, and it might be concealed from view by uncertainties in this parameter. However, we now have epidemics that have played out at least in part down the exponential decline of a suppressed outbreak. This exponential decline focusses attention on the long tail of the time to death distribution, as in this situation there were many more new cases 20 and even 30 days ago than there were 10 days ago.

It is a simple matter to convolve the observed case statistics (faulty as they may be) with a hypothesised time to death distribution to see what this predicts for the shape of the death time series, and look at how this compares to what was actually seen.

In all the plots below, the red case data is expected to lead the other two curves, but the timing of the predicted deaths (black) really should match the observed deaths (blue). I have normalised the curves to have the same magnitude in order to aid visual comparison. Cumulative curves are probably easier to assess by eye but I'll present the daily values too. Click to see larger, I deliberately kept them a bit small as there are so many.

Here are results from Wuhan (Hubei province actually):


 and South Korea

 and Italy

 and Spain




There is a very clear tendency towards too big a lag in the predicted death data (black line) according to the time-to-death curve. South /Korea is pretty good but the others are quite badly wrong. Deaths should still be shooting up in Italy and Spain!

As a simple test, I just arbitrarily shortened the time-to-death by using a gamma distribution with the same overall shape but 15 days shorter mean time (yes, 3.8 rather than 18.8 - I said, just an arbitrary test, and a rather extreme one). And look...













It's a startling improvement in all cases except South Korea. SK may be a bit unrepresentative with the early spread focussed among young church-goers (probably a low death rate and I'm guessing a long time to die) but I shouldn't throw out all inconvenient results just because it suits me.

I don't have an answer to this, but I think it's an interesting question. A 3.8 day mean time to death is surely unrealistically short as an average (from onset of symptoms). It may be the case that the time to get test results accounts for some of this difference - ie, the case numbers reported in these data sets are already somewhat delayed relative to symptoms developing, due to the time it takes to process and report cases. But is there really 15 days added delay here? It seems unlikely.



Saturday, April 04, 2020

Are we achieving suppression in the UK?

Some people are talking about a peak in deaths over the next week or two. And suppressing the disease such that we might be able to lift the lockdown eventually (end of May has been mentioned). But when will we tell how well the controls are working? I've already shown that it takes a while to really show up convincingly in the daily death data (and I haven't even taken account of additional delays in reporting deaths once they've happened). Clearly the reported case numbers for illness are wholly inaccurate and reflect testing policy more than anything else.

Here are a few hypotheticals, again based on the SEIR model fit to daily death data. This time, I have tried a few scenarios in which the effectiveness of the lockdown is measured as a scaling on the R0 parameter. In order from top to bottom, in the plot below R0 stays the same, then is reduced by 50%, 75% and 95% of its value, at the date of the lockdown. The two end cases are IMO wholly implausible and I've put them in just to show the bounds of what the model could do.

First the daily death figures, with recent observations shown. The data don't distinguish much between the scenarios, though you could say it's slightly optimistic that the most recent couple of days are towards the lower end of the spread, but still compatible with even the blue plume at this stage. The plumes do separate out a bit and in another week we will have a much clearer picture of which is closest to the truth. The original R0 in this fit is 3.3 so the green plume with Rt = 1.6 is still trending up towards a widespread epidemic whereas the red one with Rt = 0.8 is achieving control at a much lower level, albeit taking a bit of time to do so. 



The yellow plume shows that people keep on dying at much the same rate for a while even if the disease is quickly reduced with Rt = 0.16. This is a consequence of the long time that some people take to die - at least according to the fitted distribution of time to death that I've taken from the Ferguson et al work. I think it may be worthwhile to revisit this time to death as it was based (as so much of the modelling has been) on a very small set of Chinese data from early on in their outbreak. Rather amazingly, the entire shape and size of the time-to-death distribution beyond 20 days is based on the extrapolation of a parametric fit into a region which contains not a single data point. I'm not saying it's wrong and they probably did as good a job as was possible with the limited data available at that time, but it could surely be improved (revalidated) now.

I've said I don't believe the data on reported cases, but the level of illness is still something that the healthcare system experiences and it's possible that this will provide an earlier indicator. If fewer people are getting ill then there will be fewer 999 calls and lowered demand on beds etc. The symptomatic cases respond much more rapidly and clearly to controls, as shown below. 



So in the red case of even marginal control, where Rt is just over 0.8 on average, the number of new cases falling ill each day should be dropping off quite quickly (consistently by about 5% per day) and by now be well down from their peak. This may not be immediately felt if the beds are still full of course, but it could be observable soon even if time from symptomatic to seriously ill is another week (some people die faster than this). It will certainly be clear soon enough if the cases are really dropping.

Friday, April 03, 2020

Lombardy Lockdown is working?!?

A few small tweaks to the model I was using yesterday and a bit more data have resulted in quite exciting results. Firstly, I improved the death calculation which I had implemented a bit naively. I didn't think it mattered, and perhaps it doesn't much, but every little helps. Now also using the Gamma infection-to-death distribution in the IC work.

I've also changed the prior a bit, I reckon with hindsight that the mean of Rt = 1.5 was a bit high. I'm really interested in whether the data point to a value lower than or higher than 1 with highest likelihood, and if the prior mean is 1.5, a posterior result of 1.3 (say) doesn't really tell you. So I have changed to N(1,1) for the prior of Rt.

And we have two more days of data, since yesterday's file wasn't quite up to date. And the result is....

Looks pretty impressive to me. The lockdown seems to be working! I wonder how believable it is? With deaths starting to trend slightly down, it does seem hard to imagine that the disease is growing.

Uncertainties in all this are a bit dodgy and open to question. But it certainly seems more likely than not that the lockdown in Lombardy has worked.

Thursday, April 02, 2020

When was the Lombardy Lockdown?

Two posts in one day! 

Wuhan has a nice long time series which gives a good test of the model and the fitting. Which it seemed to pass pretty well. Next on the hit-list was Lombardy, which has had the second-longest lockdown. Jules and I thought (yes this was a prior decision) that it was more logical to consider Lombardy as a region rather than Italy as a whole, as this one area was so far ahead of the rest of the country and also had substantially different policies in place with restrictions coming in at an earlier date, Actually, defining the date of lockdown wasn't that easy in this case. According to the Wikipedia page, the restrictions started on 21 Feb, were extended around the end of Feb and early March, and then went national around 8 March. While it may seem inappropriate to consider highly localised restrictions, if they covered the majority of the outbreak, then this will dominate the behaviour of the total exponential growth until the rest catches up. So it really isn't that obvious (without doing a lot more detailed investigation) what date to use. I'll try a range.

8 March is when the lockdown was imposed firmly across Lombardy (there's a later date of 11 March that applies nationally) and perhaps this is the biggest effect on R. Though some things took place earlier in March, like some schools and universities being shut down. Here are results when I assume a date of 8 March for the change in R:


So...it sort of looks like there's a bend, and the model makes some effort to fit it, but tends to overshoot a bit. Marginally narrowed compared to the prior which was N(1.5,1) remember. And nudged down a bit to lower values but still likely greater than 1. Though some people will eye up those data points and the median of the result and say the model is being unreasonably pessimistic. Time will tell. I'll come back to this in a minute.

But first, back on 22 Feb there was a local shutdown around the area of the initial outbreak. If I try this date, on the assumption that even though it's small, control in this area might dominate the overall growth, then I get the following: 

Amusing but nonsense. Rt is actually pushed up compared to my prior guess, and while that is not impossible, what really appears to be happening is that the model is trying to fit the initial uncontrolled fast growth phase with this post-intervention value. It really seems to be saying that the 22 Feb interventions did not have a strong effect.

So, how about some time in between? Choosing the mid-point of 1 March gives the following:



 It's a decent fit, Rt is pushed down and constrained a bit. Of course post-hoc playing around with the intervention date isn't really in the spirit of objective scientific testing and research, but I never claimed to be doing objective scientific research so that's ok. There is obviously a risk that I've overfitted the intervention date but it looks to me like this simulation provides a pretty good explanation of what we can see in the data.

As with Wuhan, 3-4 weeks of data after the intervention date is showing some evidence of a reduction in R0, probably below even the prior assumed mean of 1.5, though maybe not below 1. This seems reasonably encouraging news at least in terms of avoiding the worst extreme of a completely uncontrolled explosive epidemic, though it may mean the controls have to be in place for a while.

What can we learn from Wuhan?

Epidemiology Thursday again, so my boss is allowing me to spend some more time doing calculations. I have set up a proper MCMC sampler to do the parameter estimation more efficiently which is also helpful for churning through a lot of different cases.

Wuhan (Hubei) is a great test case for this sort of modelling and parameter estimation as we have a good long time series, albeit I know people have some doubts about the reliability of their data in current situation but that doesn't really affect what I'm doing. There are doubts about the reliability of all data in this sort of situation!

The main aims in this post are to see:

(A) Can we model the epidemic curve, including the effect of the lockdown?
(B) What parameters do we obtain from fitting the data?
(C) How soon after the lockdown can we estimate its effectiveness, and how do rapidly we learn over time?

C is really the fun one. But starting at the beginning...

The model is the same old SEIR, with uncertainties primarily in R0 during the initial growth, and Rt after the lockdown (that's still the reproductive ratio R but it applies after a given time t). I have to also allow the initial state to vary (it is strongly related to R0) in order that we get the right size epidemic at the right time. Including uncertainty in the time constants of the model doesn't actually make much real difference here, changes in R compensate for changes in the reproductive time period anyway. I've allowed for half a day of uncertainty (at one standard deviation) in the latent period. The other parameters that I could have allowed to vary are held fixed for now, there's no point having too many more degrees of freedom than constraints as it just means the estimation process is a bit less well behaved.

For priors in R0 and Rt, I used N(mean, sd) = N(3,2) and N(1.5,1) respectively - the prior for R0 is very broad and although I do expect the lockdown to have an effect, I don't know how much. Rt could easily be less than 1 in this prior, but it's rather more likely to be greater than this value. I could have implemented this change as a scaling factor on R0 which might have been better, but I'm not altering it now.

I should also mention for completeness that I massaged the Hubei data just a little to remove obvious errors. I could only find daily cumulative data starting on the 23rd Jan, so I had to add a couple of weeks prior to that as the first death was reported on the 9th and there were 17 by the 23rd. I also eliminated some obviously erroneous zeros in the daily death figures by shifting 1/3rd of the deaths from the two immediately adjacent days. I'm working primarily in log-space so zeros are particularly irksome. And anyway they are just missing data, not real days with no deaths (except at the very start).

And without further ado,..the fit to the Hubei data. With lockdown date (when R changes) indicated. And the before/after values shown on the plot. As previously, the plume is 5-95% range, I've put the median as the magenta line and the thin blue ones are just a few random samples from the ensemble. The red circles are the obs after the massaging I just described.



That's pretty amazing IMO for such a simple model to do such a good job. The parameter values seem very reasonable too, compared to the literature. Yes I know there's a bit of extra whitespace at the top of the plot, the reason is to put all of them on the same axes for easy comparison.

So.....how quickly could we learn about Rt, the value after the lockdown? One week after, perhaps? Here's the result when we use the first 7 days of data after that date. 


It's a big fat nope. The Rt value shown there is simply the prior choice I made (slight variation due to sampling/rounding error, the exact mean in the ensemble here is 1.551!). Therefore this plot shows the range of futures I assumed under my prior assumption, with some trajectories heading up to herd immunity without much hesitation, and others bending down towards suppression very rapidly. 

How about 2 weeks of data? Does that give us a hint? 



No, not really. If you know what is coming, you could claim there's a slight reduction in the growth rate over the last few days of data, and this is reflected in a slightly lower estimate of R0, but it doesn't constrain Rt. In this model. In fact the time scale for Rt to affect deaths in this model does theoretically allow for it to influence the last few data points and thus for these data to help in the estimation of Rt, but they aren't sufficiently strong evidence to really tell us much.

Three weeks post lockdown, however, and there appears to be something going on. Rt is nudged lower than the prior estimate, and its uncertainty is starting to reduce a touch. Even in the worst case of an ongoing epidemic, this is slower and smaller than in the previous results.




4 weeks, and it's clearer, but even then there is still a substantial chance of Rt > 1. 

5 weeks down the line, it is much better bounded:

That was a very unlucky sample line that shoots off upwards. I checked what was going on in case of a bug or lack of convergence and it just happened to lie at the 99.8th percentile of the posterior 5000-member ensemble. But even so, it's not getting to the top of the graph by the end of March.

Usual caveats apply: it's a simple box model, and I've made a few assumptions that might not be fully justified. But overall I'm really impressed with how well it works.

Tuesday, March 31, 2020

The new study from IC

I haven't had time to examine this new research in great detail but it looks pretty good to me (maybe one or two minor caveats). They have fitted a fairly simple mechanistic statistical model to time series data for deaths in a large number of European countries, using a Bayesian hierarchical modelling approach to estimate the effects of various "non-pharmaceutical interventions" like shutting schools and banning large meetings which have been widely adopted. I suspect (could be wrong tho) that the reason for using a primarily statistical model is that it's quicker and easier for their method than a fully dynamical model like the one I'm using. They say it is fundamentally similar to SIR though, so I don't think their results should be any the worse for it (though perhaps the E box makes a difference?).

First and foremost, by far the most significant result IMO is that they find the initial R0 to be much larger than in the complex mechanistic modelling of Ferguson et al. This confirms what I suggested a few days ago, which was based on the fact that their value of R0 = 2.4 (ish) was simply not compatible with the rapid exponential growth rate in the UK and just about everywhere else in the western world.

Here are pics of their values for a number of countries, for their preferred reproductive time scale of 6.5 days. Left is prior to interventions, right is current post-intervention estimates.

 


You may spot that all the values on the right imply a substantial probability of R0 > 1, as I also suggested in my recent post. If the initial R0 is high, it's hard to reduce it enough to get it below 1. If true, that would mean ongoing growth of the epidemic, albeit at a lower rate and with a longer lower peak than would be the case without the interventions. I will show over the next few days what the implications of this could be over the longer term.

It's also important to realise that these values on the right depend on the data for all countries - this is the "hierarchical" bit - there is no evidence from the UK data itself of any significant drop in R, as you can work out from the below plot which gives the fit and a one-week forecast. There simply hasn't been long enough since the restrictions were put in place, for anything to feed through into the deaths. Though if it bends in the near future, that will be a strong indication that something is going on. They appear to be 3 days behind reality here - it looks like they haven't got the two small drops and most recent large rise.




Despite my best intentions, I don't really see much to criticise, except that they should have done this a week or two ago before releasing the complex modelling results in which they made such poor parameter choices. Obviously, they would not have been able to estimate the effect of interventions using this method, but they'd have got the initial growth rate right.

It's a little like presenting forecasts from one of the new CMIP6 models, without checking that it has roughly the right warming rate over the last few decades. You'd be better off using an energy balance that was fitted to recent history, especially if you were primarily interested in the broad details such as large-scale temperature trend. In fact that's a pretty direct analogy, except that the epidemic model is much simpler and more easily tuned than a state of the art climate model. Though they also had a bit less time to do it all :-)

As for the caveats however - I had to laugh at this, especially the description above:



I remember seeing a similarly impressive "large correspondence" plotted in a climate science paper a while back. I pressed the authors for the correlation between the two sets of data, and they admitted there wasn't one.

But that's a minor quibble. It's a really interesting piece of work.

I am now just about in a position to show results from my attempts at doing proper MCMC fits to time series data with the SEIR model so it will be interesting to see how the results compare. This can also give longer-term forecasts (under various scenario assumptions of course).

Sunday, March 29, 2020

Is R0 larger than people think?

EDIT: as a quick update to this, I'd just like to point to the new IC report out here, which vindicates what I've written below. Yes, R0 is higher than they thought, and this does mean that control is far harder to achieve (and quite possibly will not be achieved with current policies).

You read it here first...

--------


Is R0 really big?  R0 here being the average number of people that each infected person passes the disease on to, assuming a population with no immunity. The reason I wonder, is that the doubling rate in the death data is well under 3 days - under 2.5 in fact over the full time series from the first death to yesterday's data. There hasn't been a single period of 6 days in which the total number of deaths in the UK has not at least quadrupled, apart from when the govt sneaked in a short 8 hour “day” into the official time series recently. I see this has been adjusted for in some data sets, but not the one I have to hand - it won't matter in the long run of course. 

Incidentally, it should not be forgotten that the official death toll now only counts hospitalised cases. Those who don't get that far, and aren't tested, aren't counted. I'll assume this is a minor factor for now and ignore it completely in this analysis.

Doubling time depends on 3 parameters in the SEIR model - R0, and the latent and infectious periods which I'll abbreviate to L_p and I_p. I_p is the time interval during which ill people are passing on the disease, which can hardly be less than a day. I'll fix it at 2, which is also much shorter than most estimates I've seen. It's important to realise that L_p and I_p are basically biological parameters of the disease. They may vary a bit from place to place, but it's hard to imagine wild variation. R0, on the other hand, is substantially a social construct. It varies with how many people we meet, and how we behave with hygiene and cleanliness. It is R0 we aim to change through social distancing and disease control measures such as quarantining etc. A value of R0 calculated from another country and culture may be rather unsuitable for our situation. The value of R0=2.4 that Ferguson et al preferred for their main calculations (with a range of 2.0-2.6) was from the first few people infected in China, and may be unrepresentative for a number of reasons.

Even with I_p fixed at a very small value, in order to get a short doubling time, you still need either to have a small value for L_p, or a large R0. Most studies put L_p around 4-5 days which would mean R0 has to be pretty huge. You do have to be a little careful in interpreting the epidemiological analyses, as people often define L_p as the time taken to develop symptoms, whereas the model actually only cares about when you start to infect others - and there is evidence that this is a slightly quicker process. But still, a latent period of 3 days would be very short compared to what most people are saying.

If the latent period is 3.5, and the infectious period is 2 - both value which seem to be below the range of widely accepted values - then R0 has to be 2.7 even to get a doubling rate of 3 days, which is a bit slower than the data suggest. R0 = 2.7 is outside the range of values that were even tested by Ferguson et al. If you take a more accurate but pessimistic view that the doubling rate is actually 2.5 (it's actually comfortably quicker even than this in the data!) then R0 has to be 3.2.

Why does this matter? The control, mostly. Ferguson et al assume we can reduce R0 substantially by the extreme control measures that are put in place. Their “suppression” strategy takes a basic value of R0=2.2 and assumes it will be reduced by about 70% though social distancing and all the other measures which are put in place.

Reducing 2.2 by 70% gets it well below 1 and control is rapidly achieved. Reducing 3.2 by 70% just gets to below 1, barely. If R0 is any bigger, then a 70% reduction still leaves it greater than 1 meaning we progress through a full epidemic, albeit a smaller one than we'd have had without the controls.

In the graphs below I've tried a range of combinations with different pairs of R0 and L_p, all of which have exactly the same doubling time of 3 days - at the slow end of what might be compatible with the data, though the quick end of what the Govt is currently planning for. These following simulations are also all tuned to have 1000 deaths at the present date.


With no intervention, the epidemics are identical at first and broadly similar in their entirety, only differing in the point at which herd immunity is reached. We see about 600,000 deaths at a 1% death rate.


If we introduce controls in the middle of March to reduce R0 by 70%, then the second plot results. What a difference! For the low R0 values, the epidemic is rapidly controlled, though the deaths lag this a bit. If R0 is large, even though the growth is completely indistinguishable up to this point, the futures are very different. The growth remains positive for a while, before the epidemic levels out and decays. Note this assumes we keep the control measures in place indefinitely.

Here's a blow-up of the last plot just focussing on the period around the present so you can see the detail. And I will also show daily deaths as this gives a clearer indication of changes.



I suspect this model reacts too sharply to changes in parameters, compared to reality where there are numerous sources of variability that this model does not capture. That tends to smooth things out rather. Nevertheless, it suggests we'll be a few days into April that we could possibly expect to recognise a levelling off in the death rate. A key question will be whether we keep below 1000 deaths per day, which I think would be an optimistic hope at this point. I emphasise that the lower lines depend on really rather extreme values for the time scale parameters of the virus, which in principle I'd think should be reasonably well understood by virologists and the like. The lower two lines have latent periods of just under 2 and 3 days respectively!

The Govt's talk of 20,000 deaths looks to be a bit hopeful to me.


Friday, March 27, 2020

More on parameters and rates

So it seems the Govt got the memo (I'll add a link when I see a write-up of Gove's comments). Though their estimate of 3-4 days for the doubling period is still longer than the data indicate - it was probably under 3 in my estimation, and that was way back in the mists of time when the latest figure for deaths was only 422. Two days on (plus a bit, shifting the reporting deadline from 9am to 5pm so really 56 hours) and it's 759. That's a doubling per 2.75 days unless my maths skillz are all gone to pot.

I hope that they have also realised the consequences for the timing and magnitude of the peak, because it really makes a bigly huge difference in bringing it all upward and forward as I showed yesterday. An awful lot hangs on whether this social distancing thing really works, because if not, it's going to be pretty dreadful. It's already pretty bad, and there are several doublings to go.

One important detail to remember in the parameter estimation/model fitting process is that with 3 parameters all playing a role in determining the net growth rate, we have a highly underdetermined problem when we look at population-level time series data alone, as this is just a simple exponential curve with a particular magnitude and growth rate. If I tell you that AxBxC = 35.4, you are still none the wiser as to what A, B, and C are. The converse of this is that this lack of determinacy isn't a huge problem in the simple simulation of the initial exponential growth curve - all that matters is the growth rate, and the data do constrain that. It does start to matter when you consider the later stages where herd immunity comes into play, or the response to a changed social situation such as controls and shutdowns. These are primarily linked only to R0.

So anyway, my first assumption when wondering what had lead Ferguson et al astray was that they might have used R0 from one epidemiological study and the latent period from another, picking values that were separately plausible but jointly wrong. I was a bit surprised to see the same ref listed for both. It reports a doubling time of 7.4 days for the early stages of the Wuhan outbreak. So, shrug. Either the data there are unreliable or the social situation is sufficiently different as to be inapplicable. In any case, it was a clear mistake to rely so strongly on one small data set. We have an abundance of empirical evidence that the intrinsic doubling time in multiple western societies (prior to any large-scale behavioural changes) is about 3 days, quite possibly even shorter.

A consequence of the parametric indeterminacy is that we don't know if the latent period is short and R0 smallish, or whether both are rather larger. This doesn't matter for the historical curve (which is of course precisely why they are undetermined) but it does when you start to consider interventions. If we assume (ok, it's an assumption, but we have little choice here with this sort of model) that a given social intervention acts so as to halve the amount of contact people have and consequently halve the value of R0, then a change of R0 from 8 to 4 (assuming reproductive time scale of 9 days) is much less effective than a change from 4 to 2 (and a 6 day time scale) or 2 to 1 (with a 3 day generation). The doubling time was initially 3 days in each case, and changes to 4.5 days in the first case, 6 days in the second and the epidemic stops in its tracks for the third. I should probably make a pic of this but I can't be bothered. Use your imagination.

It seems to me that the only way these parameters can be reasonably estimated is from detailed case studies, as in the cited Wuhan study which actually looked at time from infection to symptoms. Though, if the disease is infectious sooner than that (as appears to be the case) this isn't really the right time scale to use anyway. Regardless of the reasons, it's doubling much more rapidly than every 7.4 days.



Thursday, March 26, 2020

Fitting the COVID-19 SEIR model to data, part 2

TL;DR version: when the model is fitted to the data, the doubling time is around 3 days. This is far more rapid than the Imperial College (IC) modelling, and (if correct) could mean that the peak would be upon us faster than they predicted.

A continuation of my investigations with the SEIR model which I have presented in previous posts especially here and here. Use the label for the full set. I'm a little surprised there hasn't been more competent statistical modelling appearing yet. I know there are plenty of extremely talented statisticians around who can do this sort of thing in their sleep, and yet the modelling that has appeared recently seems a bit...limited. The Oxford study was plainly silly (at least in its presentation) and even the IC study which is the direct foundation of Govt policy didn't do any sort of meaningful uncertainty analysis. I'm a bit less embarrassed to show my stuff given the level of what the pros have done.

I am however a little embarrassed that it took me so long to discover a significant issue with the IC study which I will now discuss. A while back, Vallance mentioned a 5 day doubling time (on this youtube video, starting from 27 mins in). And people were puzzled. Where did this come from? The data show a much more rapid growth, eg on this excellent analysis which is regularly updated in the FT:



Most countries are close to doubling every 3 days, which is a much more alarming position.

Well, it turns out that 5 days is baked in to the parameter values of the IC model! The doubling time is a function of R0 and the latent and infectious periods, which they set according to their interpretation of the literature and research. Importantly, through their use of a single data point to calibrate their model, they made no attempt at all to fit the observed growth rate, merely setting the start point of the epidemic in their model at the level that ensures their simulation had the correct number of total deaths by mid-March). So, I am trying to tune my model in more detail in this post to take account of the information we have regarding the growth rate.

I've taken a fairly subjective view, which I believe is necessary given the limitations of both the data and the model. The simple model I'm using cannot be expected match the data precisely but I'm trying to get the general gist of it, and allow a reasonable range of plausible model fits. I have changed my priors from the previous attempt, in order to provide more support for shorter doubling times as this is clearly justified by the international data shown above. I have also done an explicit calculation of deaths (based on a proportion of infected cases dying after a significant time delay) and am using both the death data and the reported number of cases in my model fitting, with the latter adjusted as reasonably as I can do for under-reporting.

With no further ado, here are the results. I'm just using a simple rejection sampling (likelihood weighting) approach, with uncertainties in R0, latent period, death rate and starting size. No real need for anything more sophisticated with only 4 parameters and a quick model. Here is my prior predictive distribution, and posterior, over the initial segment where we have observations, presented in log-space as this changes the exponential growth to linear and clarifies what the model is and isn't doing....

Total (UK) cases first. All ranges presented are 5-95%, so there is a chance of being outside, but not a great one:

I have just used the start and end of each observational time series, there's no point in using all the intermediate data points as they are very close to a straight line anyway. Although, one could probably argue that the growth rate in the very very earliest stage might be different to that generally applicable (eg due to contact tracing and the atypical nature of incoming infections), and I might come back to that in future work. For the number of reported cases, I've assumed that the early numbers were fairly reliable and the later one substantially under-reported along the lines of Kucharski, as I mentioned earlier. I should also mention that I make a further assumption of 50% of cases being symptomless and a death rate of about 1% (also an uncertain parameter in this work) for the symptomatic cases, ie 0.5% of all infections. I do have some uncertainty in this death rate parameter so it could be a bit higher or lower. In the above graph I have plotted the official data as black circles, the red H marks indicate the estimates I'm using in the fit, together with their uncertainties (ie, the uncertainties that I have subjectively attached to them, these are not truly uncertainties of the data as such). There is lots of room to quibble with my judgements!

Total deaths next, and I set the axes to give the same aspect ratio, which is why the labels are a bit odd (ok, laziness played a rĂ´le too):

For these death data, I'm assuming the reported values are themselves reliable (up to the day before yesterday!) so not adjusting them directly. But still allowing for significant mismatch in the fit. 

The fitted parameters suggest a slightly higher R0 of about 2.6 and a latent period of only 2.5 days. These, together with the infectious period which I fixed at 4 days, combine to give a doubling interval of just under 3 days which is even shorter than in my adjusted prior let alone the 5 days of the IC report. Clearly these data (at least, my interpretation of them) point towards a rapid growth rate. You can see in the above graph I'm still not matching the growth rate in deaths which is actually doubling every 2.2 days over this period. These death data seem to be going up rather faster than the reported cases even after I have adjusted the latter for estimated under-reporting, and the model does not achieve this. I don't think this can be a model limitation per se as simple theory implies that all the growth rates of all variables should be the same in this early stage of exponential growth - cases, symptomatic cases, serious illness, death, recoveries. All of them should double at the same rate, and therefore be linked by simple proportionality rules. But in fact the first few deaths were at a rate of 1 per hundred current cases, now it's more like 5. Maybe the under-reporting is worse than I think? Food for thought surely.


Now some projections for the coming summer. Here are the epidemic curves under a business-as-usual course of taking no action. I'm showing first the total number of cases, then the total deaths. In the first of these, I also showed a handful of the individual curves to show the typical shape - as I mentioned previously, it is not the case that we might follow along the lower bound of the 5-95% range, this is instead the envelope within which most of the curves sit. All of them give a large rapid epidemic.


It's over by the end of May in terms of infections - even the start of May in some cases - but the bodies will probably continue to pile up for another few weeks after that:






Remember that this is just assuming a 0.5% death rate, which is not particularly pessimistic. I think the model probably over-predicts the proportion of us that get the disease, but not by a huge margin. And 0.5% might well be optimistic in such a compressed time frame when we don't have capacity to treat people. If you don't like the notation for the axis, it's about 250,000-500,000 deaths.

Now, the obvious question is, what should we do about this? The results of Ferguson et al suggest that changes in social behaviour can reduce R0 significantly, and this is easy to simulate, though the reduction itself is just something I have to impose (based loosely on fits to Ferguson et al).

Here are two possible futures. I've assumed we either achieve (a) a 50% reduction in R0, which is pretty optimistic, or (b) a 70% reduction in R0, which is incredibly optimistic. Both of these are applied for a 3 month period (the shaded region) and then released. The results are.....not really all that great. 





If we suppress strongly now (and then relax the controls), the epidemic comes back for a second go later, rather as Thomas showed. The red curve is slightly better than the others, but the total deaths aren't that different. Of course this analysis completely ignores the possible development of treatment or vaccines that could completely change the picture. Let's hope so cos these pictures are not really very appealing!



I'm sure I could think of more things to do, but this is enough for now. Comments welcome. It's a shockingly compressed time scale compared to what the IC group have predicted - they had critical care demands peaking in June! - so I hope their judgement over parameter choices proves to be more appropriate than what this model fitting achieved.