Wednesday, August 05, 2020

Could R still be less than 1?

It's been suggested that things might all be fine, maybe the increase in case numbers is just due to more/better testing. There certainly could be a grain of truth in the idea, as the number of tests undertaken has risen a little and the proportion of tests that have been positive has actually kept fairly steady over recent weeks at around 1%. On the other hand, you might reasonably expect the proportion of positives to drop with rising test numbers even if the number of ill people was constant, let alone falling as SAGE claim - consider at the extreme, 65 million tests couldn't find 650,000 positives if only a few tens of thousands are actually ill at any one time. Also, the ONS pilot survey is solid independent evidence for a slight increase in cases, albeit not entirely conclusive. But let's ignore that inconvenient result (as the BBC journalist did), and consider the plausibility of R not having increased in recent weeks. 

This is fairly easy to test with my data assimilation system. I can just stop R from varying at some point in time (by setting the prior variance on the daily step to a negligible size). For the first experiment, I replaced the large jump I had allowed on 4th July, with fixing the value of R from that point on. Note however that the estimation is still using data subsequent to that date, ie it is finding the (probabilistic) best fit for the full time series, under the constraint that R cannot change past 4 July. I've also got a time-varying case ascertainment factor which I'll call C, which can continue to vary throughout the full interval.

Here are the results, which are not quite what I expected. Sure, R doesn't vary past the 4th of July, but in order to fit the data, it shoots up to 1 in the few days preceding that date (red plume on 2nd plot). The fit to the death data in the bottom plot looks pretty decent (the scatter of the data is very large, due to artefacts in the counting methodology) and also the case numbers in the top plot are reasonable. See what has happened to the C factor though (red plume on top graph). After being fairly stable through May and most of June, it takes a brief nose-dive to compensate for R rising at the end of June, and then has to bounce back up in July to explain the rise in case numbers.

While this isn't impossible, it looks a bit contrived, and also note that even so, we still have R=1, firmly outside the SAGE range of 0.8-0.9. Which isn't exactly great news with school opening widely expected to raise this value by 0.2-0.5 (link1link2).

So, how about fixing R to a more optimistic level, somewhere below 1? My code isn't actually set up very well for that specific experiment, so instead of holding R down directly, I just put the date back at which R stops varying. In the simulations below it can't change past the 1st June. It still climbs up just prior to that date, but only to 0.9 this time, right at the edge of the range of SAGE values. The fit to the death data is similar, but tis time the swoop down for C on the upper plot is a bit more pronounced (because R is higher through June) and then it has to really ramp up suddenly in July to match the rise in case numbers. You can see that it starts to underestimate the case numbers towards the present day too, C would have to keep on ramping up even more to match that properly.

So R being in the SAGE range isn't completely impossible, but requires some rather contrived behaviour from the rest of the model which doesn't look reasonable to me. I don't believe it and think that unfortunately there is a much simpler explanation for (some of) the rise in case numbers.

More what-ifs

It was pointed out to me that my previous scenarios were roughly comparable to those produced by some experts, specifically this BBC article  referring to this report. And then yesterday another analysis which focussed on schools opening.

The experts, using more sophisticated models, generated these scenarios (the BBC image is simplified and the full report has uncertainties attached):

and for the schools opening report:

The tick marks are not labelled on my screenshot but they are at 3 month intervals with the peaks being Dec on the left hand and March on the right hand panel.

While these are broadly compatible with my analyses, the second peak for both of them is significantly later than my modelling generates. I think one important reason for this is that my model has R a little greater than 1 already at the start of July, whereas they are assuming ongoing suppression right through August until schools reopen. So they are starting from a lower baseline of infection. The reports themselves are mutually inconsistent too, with the first report having a 2nd peak (in the worst case) that is barely any higher than the first peak, and the second report having a markedly worse 2nd peak, despite having a substantially lower R number over the future period that only briefly exceeds 1.5. It's a bit strange that they differ so significantly, now I think about it...I'm probably missing something obvious in the modelling.

Of course in reality policy will react to observations, so all scenarios are liable to being falsified by events one way or another.

Sunday, August 02, 2020

What if?

It's a while since I did any real forecasting, the current system just runs on a bit into the future with the R values gradually spreading out due to the daily random perturbations, and the end result is pretty obvious. Now with the effective R value probably just above 1, and various further relaxation planned (e.g. end of furlough, schools returning) jules thought it would be interesting to see what might possibly happen if R goes up a bit.

Here are two ensembles of simulations, both tuned the same way to historical data, which gives an R_effective of about 1.1 right now. The step up on the 4th July is a modelling choice I made through choice of prior, in allowing a large change on that one day only rather than a gradual ramping up around that time. In the first set of forecasts, I ramp up R by 0.5 over 30 days through September. For the modellers, I'm actually using R as my underlying parameter, calculating R_effective based on the proportion of people who (it is assumed!) have acquired immunity through prior infection. So typically the underlying R value is going up from 1.2 to 1.7 or thereabouts. You can see the resulting ramp up in R_eff on the plot, with the subsequent drop entirely due to the herd immunity factor kicking in as the second wave peaks. The new peak in deaths is...not pretty. I'm disappointed it is so severe, in my head I'd been assuming that a much lower R number (compared to the 3-3.5 at the start) and non-negligible level of current immunity would have helped to keep it lower.

The second set of results is a more optimistic assumption where R only goes up by 0.2, this time in a single step when the schools go back near the start of Sept (don't quote me on the date, it was just a guess).'s still not great I'm afraid. The lower R gives a more spread out peak and there is a chance of things turning out not too badly but a lot of the trajectories still go up pretty high, with most of them exceeding the April max in daily deaths, and sustaining this for quite a while.

So...that's all a bit of a shame. There are however reasons why this may be a bit too pessimistic: it is well-known that this simple model will overestimate the total penetration of the disease as it doesn't account for heterogeneity in the population, which could make a significant difference. Also, I've kept the fatality rate at 0.75% despite advances in treatment which have definitely nudged it lower than it was at the start. On the other hand, the model does not account for loss of immunity here among people who have had the virus. Not clear if that simplification is truly valid over this time scale.

Anyway, these are not predictions, I just put in some reasonable-sounding (to me!) numbers to see what would happen. It does look like any further significant increase in R will have serious consequences.

Thursday, July 30, 2020

The price of freedom

The Govt changed the lockdown rules substantially from the 4th July, with pubs, restaurants reopening and a new “1m plus” rule to replace the previous 2m distancing requirement. Predictably, the tabloids announced a new free-for-all which they labelled “Independence day”.

Up to this time, the R number had been fairly stable at around 0.8, meaning that each infected person would pass the disease onto less than one person on average and the rate of illness (and death) was dropping fairly steadily at about 20% reduction per week.

Below is how my model fits to the data up to 4th July (red circles in both plots). You can click on the plots for bigger and clearer versions. The left hand plot is daily reported cases, and the right hand plot is daily deaths. The green plume shows the model fit to each of these, with a few lines from the ensemble drawn on (dark blue) and the median prediction in magenta. The thin blue plume is the total modelled number of new infections each day, which is much higher. There is also a red plume on this plot representing the “case ascertainment factor”, ie the proportion of infections that is actually observed. This uses the scale on the right hand side of the plot, and so rises from about 1% at the start of the epidemic, to around 10% now. The blue circles represent data that had not been observed by the 4th July, and you can see in the LH plot that they tend to drift above the model forecast.

On the right hand plot, the red plume is the R number (which again uses the axis on the right hand side of that plot). It starts off around 3ish, then drops sharply when the lockdown controls were imposed, and wobbles around a bit after that point. The “current” number quoted there (mean and range) is the estimate as of the 4th July. The data observed subsequent to that date agree better with the model than was the case in the LH plot, but still look to be more above than below the forecast.

Redoing the analysis as of yesterday's data (i.e. including all data points in the estimation), and we get the following:

Now the rise in cases is reflected in the LH graph, and the corresponding rise in R is shown at the bottom of the RH plot. R is probably greater than 1, meaning that the epidemic is starting to take off again. It seems that something happened around the 4th of July to increase the rate of infection. I wonder what that could have been?

So, emboldened by these results and Peter's comment below we can try adding in a step change on the 4th July - this is just a high variance step in the prior, I'm not imposing a rise specifically, just allowing a large change. This generates the result below and it looks like a rather better fit especially to cases. However I'm not really that confident about what is going to happen and especially wouldn't be surprised if there is a bit of a decoupling between cases and deaths due to differences in the age range of people infected (eg mostly younger working age with a much lower fatality rate).

Thursday, July 23, 2020 Back to the future

Way back in the mists of time (ie, 2006), jules and I saw what was going on with people estimating climate sensitivity, and in particular how this literature was interpreted by the authors of the IPCC AR4. And we didn’t like it. We thought that any reasonable synthesis should consider the multiple lines of evidence in a coherent fashion in order to form a credible overall view. This resulted in the paper "Using multiple observationally‐based constraints to estimate climate sensitivity" described in this blog post (paper here), which people unfamiliar with the story might like to glance at before progressing further…

It’s fair to say that our intervention was not met by universal approval at the time, with the established researchers mostly finding excuses as to why our result might not be entirely trustworthy. Fine, do your own calculations, we said. And they didn’t.

Time passed, and a new generation of people with different backgrounds became interested in estimating climate sensitivity. The World Climate Research Program (WCRP) made it a central theme in one of their Grand Challenges in climate science. There were a couple of meetings in Ringberg that jules and then I attended sequentially.

In 2016, several of leaders of this WCRP steering group wrote a paper which kicked off a project to perform a new synthesis of the evidence on climate sensitivity. Their idea was to form an overall synthesis of the multiple lines of evidence, roughly along the lines that we had originally proposed, but in a far more comprehensive and thorough fashion. This is something that the IPCC isn’t really equipped to do, as it just assesses and summarises the literature. The project leaders considered three main strands of evidence: that arising from process studies (ie the behaviour of clouds, including simulations from GCMs), the transient warming over the historical record, and paleoclimate. Jules was one of the lead authors for the paleo chapter, but I wasn’t involved at the outset. However when invited to join the group I was of course happy to contribute to it, having thought about the problem off and on for the past decade.

Writing it was a lengthy and at times frustrating process, due to the huge range of ideas, topics, backgrounds and knowledge of the author team. That is also what gives this review its strength, of course, as we have genuine experts in multiple areas of modelling and data analysis, covering a huge range of time scales and techniques, and the different perspectives meant we gave each other quite a workout in testing the robustness of our approaches and ideas. During the 4 year process we had regular videoconferences, typically 9pm UK time, being 6am for Japan, 10am in Australia and afternoon for the continental USA. Luckily we had an 8-9h gap in the global spread so no-one actually had to get up in the middle of the night each time! We also had a single major writing meeting in Edinburgh in summer 2018 which almost all the main authors were able to attend in person, and a handful of "meet-ups of opportunity" when subsets happened to go to other conferences. In all, it was good practice for the new normal that we are enjoying due to COVID.

The peer review was probably the most extensive I’ve experienced, with something like 10 sets of comments – this was something we were all keen on, as we suspected it would be beyond the compass of just the usual 2-3 people. Comments were basically encouraging but gave us quite a lot to work on and in fact we reorganised the paper substantially for the better resulting in the 2nd set of reviews being very positive. Finally got it done a couple of months ago and it was accepted subject to very minor corrections (which were mostly things we had spotted ourselves, in fact).

The new paper has now been published, actually I’m not entirely sure it is up yet (minor snafu on the embargo timing) but anyone who needs an urgent look can find it here. I may write more on the details if pressed, but for now here is a quick peek at the main results:

The "baseline" calculation is what we get from putting together all the evidence, with a resulting 2.6-3.9C "likely" range. The coloured curves are various sensitivity tests, with the purple line at the top defined as the range from the lowest 17th percentile, and the highest 83rd percentile, across these tests. This isn’t really a probability range and doesn’t correspond to any particular calculation.

Tuesday, July 21, 2020

That Russian Report, in full, in brief

We hear no evidence of Russian interference. We discuss no evidence of Russian interference. We see no evidence of Russian interference.

Sunday, July 19, 2020

Patrick Vallance's faulty memory

On reflection, perhaps it shouldn't be surprising. We expect the Chief Scientist to be a genius with a brain the size of a planet who is perpetually on top of their game, but in fact they are a human frequently operating under great stress, and fallible like the rest of us. Nevertheless, his first responsibility - and ours - is to the truth, and it is therefore my task to explain that he unfortunately misled the House of Commons Science and Technology Committee when he appeared before it on Thursday 16th July.

The topic under consideration is SAGE's recommendations around mid-March, when the various restrictions were being introduced - some have argued (and I'm among them) that this happened rather too late, with the result that the country suffered many more deaths, and far greater economic damage, than would have been the case with prompt action.

Most of the interesting action during his appearance was under questioning from Graham Stringer MP, from about 50 minutes in to the video, or Q1041 on the transcript. Stringer is pressing him on the promtness (or otherwise) of introducing the lockdown, and particularly the speed of response to the data showing more rapid doubling than they had originally assumed:
Q1041 Graham Stringer: As a scientist, I was always taught to forget hypotheses, theories and ideas and look at the data, because having preconceived ideas can distort the way you look at things. When we went into this, scientists in this country were looking at data from China that showed a doubling of the infection every six or seven days. When you looked at our data closely, the infection death rates were doubling every 30 to 36 hours. Why didn’t you and SAGE advise the Government to change their attitude because, if you had looked at that and given that advice, the lockdown might have happened earlier?
To start with, to avoid the usual tedious ducking and weaving from the usual tedious suspects, it's important to be clear about the terms. When Stringer and Vallance are talking about “lockdown”, they mean the strict policies from the 23rd March onwards, when we were told to stay at home, all non-essential shopping and travel was forbidden, etc. As Vallance puts it:
there was a series of steps in the run-up to lockdown, which started with the isolation of people who had come from China, but the main ones were: case isolation; household isolation; and recommendations not to go to pubs, theatres and so on.
So, “lockdown” here means policies of the 23rd March, as also confirmed by Hancock in Hansard:
the level of daily deaths is lower than at any time since lockdown began on 23 March.
Sorry for this tedious pedantry, but experience has shown some people will, having lost the argument about timing, duck and weave about what "lockdown" means in the first place.

So, back to the timing. Vallance's main claim, which I will argue is incorrect, is contained in the following sentences:
When the SAGE sub-group on modelling, SPI-M, saw that the doubling time had gone down to three days, which was in the middle of March, that was when the advice SAGE issued was that the remainder of the measures should be introduced as soon as possible. I think that advice was given on 16 or 18 March, and that was when those data became available.
Note how clear he is that this advice to introduce the remainder of the measures - ie implementation of the full lockdown - was based on the realisation that the doubling time was as short as 3 days. I'll let him off with his use of “had gone down to” - in reality the doubling time had not changed at all, it was just SAGE's realisation that had gone down, but I will be generous and attribute this to sloppy language. He emphasises this reliance on the new data repeatedly:
Sir Patrick Vallance: Knowledge of the three-day doubling rate became evident during the week before. 
Q1042 Graham Stringer: Did it immediately affect the recommendations on what to do?  
Sir Patrick Vallance: It absolutely affected the recommendations on what to do, which was that the remaining measures should be implemented as soon as possible. I think that was the advice given.
and again:
Sir Patrick Vallance: The advice changed because the doubling rate of the epidemic was seen to be down to three days instead of six or seven days. We did not explicitly say how many weeks we were behind Italy as a reason to change; it was the doubling time, and the realisation that, on the basis of the data, we were further ahead in the epidemic than had been thought by the modelling groups up until that time.
So he is absolutely certain that the advice to proceed full steam ahead on the lockdown was predicated on the new 3 day doubling time.

However, he also claimed that this advice was given “on 16 or 18 March.” This is the critical error in his statements, that prompted this blog. Some people have jumped on this claim (and to be fair to Vallance, he was obviously unsure of the exact date in his response) to argue that the Govt was slow to react to SAGE's recommendation, and that this was the cause of the late lockdown and large death toll.

Unfortunately, Vallance was mistaken with his dates. In fact, SAGE actually still thought the doubling time was 5-6 days on the 16th March (minutes):
UK cases may be doubling in number every 5-6 days.
and by the 18th March their estimate was even slightly longer (minutes):
Assuming a doubling time of around 5-7 days continues to be reasonable.
It is therefore not at all surprising that the minutes of these two meetings do not contain any recommendation, or even a hint of a suggestion of a recommendation, that we should proceed with haste to a full lockdown. In fact the minutes of the 18th March make the very specific and detailed recommendation that schools should be shut, with the clear statement that further action would only be necessary “if compliance rates are low” (NB compliance with all measures has been consistently higher than in the modelling assumptions):
2. SAGE advises that available evidence now supports implementing school closures on a national level as soon as practicable to prevent NHS intensive care capacity being exceeded.
3. SAGE advises that the measures already announced should have a significant effect, provided compliance rates are good and in line with the assumptions. Additional measures will be needed if compliance rates are low.
Incidentally, this is why we have to be precise about what “lockdown” means, so that certain people don't pivot to “Aha! They said we should shut something! Vallance was right all along!” SAGE here is not recommending “lockdown” in the sense used by Vallance, Stringer, Hancock, or anyone else. They are only recommending school closures, which the Govt did implement promptly at that time.

Now let's go back to this from Vallance:
When the SAGE sub-group on modelling, SPI-M, saw that the doubling time had gone down to three days, which was in the middle of March, that was when the advice SAGE issued was that the remainder of the measures should be introduced as soon as possible.
The relevant SPI-M meeting at which they reduced their estimate of doubling time was actually on the 20th March (minutes). At this meeting, they abruptly realised:
Nowcasting and forecasting of the current situation in the UK suggests that the doubling time of cases accruing in ICU is short, ranging from 3 to 5 days.
The observed rapid increase in ICU admissions is consistent with a higher reproduction number than 2.4 previously estimated and modelled; we cannot rule out it being higher than 3.
All well and good, but a week late.

The nest SAGE meeting was on the 23rd (21st-22nd was a weekend) and at this point they conclude (minutes):
The accumulation of cases over the previous two weeks suggests the reproduction number is slightly higher than previously reported. The science suggests this is now around 2.6-2.8. The doubling time for ICU patients is estimated to be 3-4 days.
(NB doubling time is in principle the same for all measures of the outbreak, ignoring transient effects as the epidemic gets established. That's why it is such a useful concept and measure.)

SAGE also state at this meeting on the 23rd:
Case numbers in London could exceed NHS capacity within the next 10 days on the current trajectory.
They don't explicitly minute the need for a tight lockdown, but certainly provide statements that point in that direction, such as:
High rates of compliance for social distancing will be needed to bring the reproduction number below one and to bring cases within NHS capacity.
It seems reasonable to conclude that the message taken from this meeting was that London at least was on the verge of exceeding capacity and that strong measures needed to be urgently taken to slow transmission. As Vallance had put it:
the remaining measures should be implemented as soon as possible.
So it seems that Vallance has described the narrative arc precisely as the minutes of all the meetings around this time describe, but for the important fact that he got the date of this final recommendation wrong. He appears to have created a false memory of a world where the heroes of SAGE worked it all out in the nick of time, and told the government....who then sat on this information and delayed lockdown for a week. It's a nice story, but it's not actually what happened. The data were certainly clear to many by mid-March (ie the 14th, prior to the famously uncalibrated runs of the Imperial College model) but SAGE resolutely ignored and rejected this evidence for a further week, and this delay caused huge unnecessary harm to the country.

This would be a minor tale of a small slip of memory, were it not for the unfortunate fact that various factions have glommed onto Vallance's statement as proof that the scientists were blameless and the Govt guilty. Most egregiously, SAGE member Jeremy Farrar tweeted:
To make the mistake that Vallance did, under pressure of live questioning, is forgivable. To double down on the error from the comfort of your own computer, when the documentation is freely available, is not. The minutes prove that SAGE did not accept the evidence of the short doubling time on the 16th and 18th March. It is quite possible that some SAGE members - perhaps including Farrar - had tried to sound the alarm about the rapid doubling at an earlier time. However, they did not carry the day and I find no evidence that they spoke up in public either. SAGE did not recommend lockdown prior to the 23rd March, however much it suits various peoples' agendas to claim so.

Saturday, July 18, 2020

Mountains and molehills

A couple of weeks ago, I heard about an issue with the way COVID deaths are counted in England. It seems that PHE are going through the lists of people who had died every day, and checking to see if they had previously had a positive COVID test. If so, they are added to the total number of COVID deaths for the day, even if they had long since made a full recovery and were run over by a bus (or died of some other illness).

Clearly this is wrong, and will tend to overstate the number of people killed by the disease. Equally clearly, there aren't many deaths falling into this situation. Take the total number of 300k positive tests, assume this means 300k people (which it doesn't, as many people are tested more than once) and that they have an average remaining life span of 40 years. Then we'd expect to see 20 of them die every day from all causes, implying about this many of the daily "COVID deaths" in the PHE stats are bogus. That took me under 5 mins to work out, so I shrugged and ignored the issue. My number might not be quite right, the 40m remaining years of life thing will depend on the precise age/gender distribution of those who have tested positive but it's hard to see it being too far wrong. In the face of 100 deaths per day, about 20 of them being erroneously counted is not a huge issue though it would become more of a problem as/when the daily death toll shrinks further. It certainly has little bearing on any retrospective analysis of the size of the outbreak so far.

Two weeks later, and Loke (who I now note is who I first learnt about this issue from) and Heneghan write an article covering this issue, and promote it all round the press. I'm sure it is just an unfortunate accident that they make it sound like it's a really big issue that is major factor in explaining why the death toll in England has remained so high, as they are surely competent enough to have reproduced the calculation I presented above. Unsurprisingly, it's been picked up by the denialist wing of the media which is desperately trying to pretend that the response in England has been anything other than absolutely terrible. It's probably worth mentioning that Heneghan has form for minimising the dangers of COVID: in this piece he argues that the fatality rate is down around 0.2% which is far below all credible estimates I've seen and implies that a very large proportion of UK population has had the disease, which is robustly refuted by a hefty pile of evidence. 

Now Hancock has called for an "urgent inquiry" into this and is using it as a excuse to halt publication of the daily statistics. Even though he's a bit dim, it's hard to believe he doesn't have any numerate advisors who could tell him why it's not that big a deal. Indeed PHE quickly put out a rebuttal which supports my analysis - they pointed out that 90% of the COVID deaths occurred within 28 days of diagnosis, and of the remaining 4000, half of them were directly attributed to COVID on the death certificate anyway. Leaving perhaps 2000 bogus deaths which should not have been added. Over a ~100 day period that's pretty much the same as the 20 per day figure I came up with above.

Compare and contrast with the known under-reporting which is clear from the total death statistics and perhaps most stark in care homes, where the total "non-COVID" deaths have a massive bump coincident with the epidemic. We know that patients were pushed out of hospitals into care homes, without any testing or facilities for safe care and treatment, and it's clear that many thousands of these people died without being counted in the COVID statistics. See the huge yellow bump in the official ONS statistics below: 

This miscoding of unrelated deaths is small beer in comparison.

One way of getting the "correct" answer would be to use excess deaths, but that involves a certain amount of statistical work (excess over what, and how is this calculated?) and is not so quick and easy to come by. So I don't know what they will come up with as a solution. Using a cut-off date might be a reasonable solution, perhaps in conjunction with death certificate where it didn't specify COVID as the cause. Ie, cut out those 2k deaths where they both (a) took more than 28d from diagnosis to death and (b) were not directly attributed to COVID by a doctor. That would seem to minimise any errors in a straightforward manner, so probably they will do something more complicated...

Friday, July 03, 2020

Ho hum

Haven't posted for a while, so how about a few minutes of James O'Brien to pass the time.

Friday, June 26, 2020 Like a phoenix redux

Even odder than finding that our old EnKF approach for parameter estimation was particularly well suited to the epidemiological problem, was finding that someone else had independently invented the same approach more recently…and had started using it for COVID-19 too!

In particular, this blogpost and the related paper, leads me to this 2013 paper wherein the authors develop a method for parameter estimation based on iterating the Kalman equations, which (as we had discovered back in in 2003) works much better than doing a single update step in many cases where the posterior is very small compared to the prior and the model is not quite perfectly linear – which is often the case in reality.

The basic idea behind it is the simple insight that if you have two observations of an unknown variable with independent Gaussian errors of magnitude e, this is formally equivalent to a single observation which takes the average value of the two obs, with an error of magnitude e/sqrt(2). This is easily shown by just multiplying the Gaussian likelihoods by hand. So conversely, you can split up a precise observation, with its associated narrow likelihood, into a pair of less precise observations, which have exactly the same joint likelihood but which can be assimilated sequentially in which case you use a broader likelihood, twice. In between the two assimilation steps you can integrate the model so as to bring the state back into balance with the parameters. It works better in practice because the smaller steps are more consistent with the linear assumptions that underpin the entire assimilation methodology.

This multiple data assimilation idea generalises to replacing one obs N(xo,e) with n obs of the form N(xo,e*sqrt(n)). And similarly for a whole vector of observations, with associated covariance matrix (typically just diagonal, but it doesn’t have to be). We can sequentially assimilate a lot of sets of imprecise obs in place of one precise set, and the true posterior is identical, but the multiple obs version often works better in practice due to generating smaller increments to the model/parameter samples and the ability to rebalance the model between each assimilation step.

Even back in 2003 we went one step further than this and realised that if you performed an ensemble inflation step between the assimilation steps, then by choosing the inflation and error scaling appropriately, you could create an algorithm that converged iteratively to the correct posterior and you could just keep going until it stopped wobbling about. This is particularly advantageous for small ensembles where a poor initial sample with bad covariances may give you no chance of reaching the true posterior under the simpler multiple data assimilation scheme.

I vaguely remembered seeing someone else had reinvented the same basic idea a few years ago and searching the deep recesses of my mind finds this paper here. It is a bit disappointing to not be cited by any of it, perhaps because we’d stopped using the method before they started….such is life. Also, the fields and applications were sufficiently different they might not have realised the methodological similarities. I suppose it’s such an obvious idea that it’s hardly surprising that others came up with it too.

Anyhow, back to this new paper. This figure below is a set of results they have generated for England (they preferred to use these data than accumulate for the whole of the UK, for reasons of consistency) where they assimilate different sets of data: first deaths, then deaths and hospitalised, and finally adding in case data on top (with some adjustments for consistency).

Screenshot 2020-06-24 21.25.10

The results are broadly similar to mine, though their R trajectories seem very noisy with extremely high temporal variability – I think their prior may use independently sampled values on each day, which to my mind doesn’t seem right. I am treating R as taking a random walk with small daily increments except on lockdown day. In practice this means my fundamental parameters to be estimated are the increments themselves, with R on any particular day calculated as the cumulative sum of increments up to that time. I’ve include a few trajectories for R on my plot below to show what it looks like.


Monday, June 15, 2020 Like a phoenix…

So, the fortnightly chunks in the last post were doing ok, but it’s still a bit clunky. I quickly found that the MCMC method I was using couldn’t really cope with shorter intervals (meaning more R values to estimate). So, after a bit of humming and hawing, I dusted off the iterative Ensemble Kalman Filter method that we developed 15 years ago for parameter estimation in climate models I must put a copy up on our web site, it looks like there’s a free version here. For those who are interested in the method, the equations are basically the same as in the standard EnKF used in all sorts of data assimilation applications, but with a couple of tweaks to make it work for a parameter estimation scenario. It had a few notable successes back in the day, though people always sneered at the level of assumptions that it seemed to rely on (to be fair, I was also surprised myself at how well it worked, but found it hard to argue with the results).

And….rather to my surprise….it works brilliantly! I have a separate R value for each day, a sensible prior on this being Brownian motion (small independent random perturbation each day) apart from a large jump on lockdown day. I’ve got 150 parameters in total and everything is sufficiently close to Gaussian and linear that it worked at the first time of asking with no additional tweaks required. One minor detail in the application is that the likelihood calculation is slightly approximate as the algorithm requires this to be approximated by a (multivariate) Gaussian. No big deal really – I’m working in log space for the number of deaths, so the uncertainty is just a multiplicative factor. It means you can’t do the “proper” Poisson/negative binomial thing for death numbers if you care about that, but the reporting process is so much more noisy that I never cared about that anyway and even if I had, model error swamps that level of detail.

The main thing to tweak is how big a daily step to put into the Brownian motion. My first guess was 0.05 and that worked well enough. 0.2 is horrible, generating hugely noisy time series for R, and 0.01 is probably inadequate. I think 0.03 is probably about ok. It’s vulnerable to large policy changes of course but the changes we have seen so far don’t seem to have had much effect. I haven’t done lots of validation but a few experiments suggest it’s about right.

Here are a few examples where (top left) I managed to get a validation failure with a daily step of 0.01 (top right) used 0.2 per day but no explicit lockdown, just to see how it would cope (bottom left) same as top left but with a broader step of 0.03 per day (bottom right) the latest forecast.

I’m feeling a bit smug at how well it’s worked. I’m not sure what other parameter estimation method would work this well, this easily. I’ve had it working with an ensemble of 50, doing 10 iterations = 500 simulations in total though I’ve mostly been using an ensemble of 1000 for 20 iterations just because I can and it’s a bit smoother. That’s for 150 parameters as I mentioned above. The widely-used MCMC method could only do about a dozen parameters and convergence wasn’t perfect with chains of 10000 simulations. I’m sure some statisticians will be able to tell me how I should have been doing it much better…

Friday, June 12, 2020

Neoliberalism Kills?

However you "solve" the problem, the pandemic was always going to be very expensive. Government mandated lockdowns might be framed as the point at which the government disrupts market forces, decides to transfer some costs away from individuals, and to strongly shape the future trajectory.

The costs of lockdown vs no-lockdown in terms of both lives and money were probably not that well understood by the decision makers. While epidemiological models can easily predict 100s of thousands of deaths they assume no changes in behaviour by citizens. Rich countries with low inequality may reasonably hope for auto-lockdown by its citizens without such a massive interference by the government.  This has probably helped Japan, and perhaps to a lesser extent Sweden (Sweden has plenty of deaths but has also possibly kept more of its economy going..?). On the other hand, in many countries most people cannot afford to lockdown, and people not at such high risk (in this case, younger people) may feel much lower motivation to do so.

But instead of really thinking any of this through it seems that our dimwitted politicians simply applied their rubbish political ideological theories. They aren't scientists and do not know that theories have to make testable predictions in order to be worthwhile.

And what happened?! In this case socialism wins while neoliberalism both kills lots and lots of people and crashes the economy (because early lockdown => shorter lockdown). Ooops!

I would think it isn't always the case; socialist governments have surely fucked up very badly in the past when faced with other problems. That's the problem with random theories based on no evidence - they only work every so often. But I am still a bit worried that perhaps neoliberalism always kills and that this is just first chance to do a proper job of it.

The big caveat in all this is that we do not know what the endgame is. If herd immunity remains the inevitable consequence, then lockdown might be viewed in terms of the effect on quality of life in terms of months rather than in total lives saved. But those months are still pretty valuable, aren't they?

Tuesday, June 09, 2020 More COVID-19 parameter estimation

The 2 and now 3-segment piecewise constant approach seems to have worked fairly well but is a bit limited. I’m not really convinced that keeping R fixed for such long period and then allowing a sudden jump is really entirely justifiable, especially now we are talking about a more subtle and piecemeal relaxing of controls.

Ideally I’d like to use a continuous time series of R (eg one value per day), but that would be technically challenging with a naive approach involving a whole lot of parameters to fit. Approaches like epi-estim manage to generate an answer of sorts but that approach is based on a windowed local fit to case numbers, and I don’t trust the case data to be reliable. Also, this approach seems pretty bad when there is a sudden change as at lockdown, with the windowed estimation method generating a slow decline in R instead. Death numbers are hugely smoothed compared to infection numbers (due to the long and variable time from infection to death) so I don't think that approach is really viable.

So what I’m trying is a piecewise constant approach, with a chunk length to be determined. I’ll start here with 14 day chunks in which R is held constant, giving us say 12 R values for a 24 week period covering the epidemic (including a bit of a forecast). I choose the starting date to fit the lockdown date into the break between two chunks, giving 4 chunks before and 8 after in this instance.

I’ve got a few choices to make over the prior here, so I’ll show a few different results. The model fit looks ok in all cases so I’m not going to present all of them. This is what we get for the first experiment:

The R values however do depend quite a lot on the details and I’m presenting results from several slightly different approaches in the following 4 plots.

Top right left is the simplest version where each chunk has an independent identically distributed prior for R of N(2,12). This is an example of the MCMC algorithm at the point of failure, in fact a little way past that point as the 12 parameters aren’t really very well identified by the data. The results are noisy and unreliable and it hasn’t converged very well. The last few values of R here should just sample the prior as there is no constraint at all on them. That they do such a poor job of that is an indication of what a dodgy sample it is. However it is notable that there is a huge drop in R at the right time when the lockdown is imposed, and the values before and after are roughly in the right ballpark. Not good enough, but enough to be worth pressing on with….

Next plot on top right is when I impose a smoothness constraint. R can still vary from block to block, but deviations between neighbouring values are penalised. The prior is still N(2,12) for each value, so the last values of R trend up towards this range but don’t get there due to smoothness constraint. The result looks much more plausible to me and the MCMC algorithm is performing better too. However, the smoothness constraint shouldn’t apply across the lockdown as there was a large and deliberate change in policy and behaviour at that point.

So the bottom left plot has smoothness constraints applied before and after the lockdown but not across it. Note that the pre-lockdown values are more consistent now and the jump at lockdown date is even more pronounced.

Finally, I don’t really think a prior of N(2,12) is suitable at all times. The last plot uses a prior of N(3,12) before the lockdown and N(1,0.52) after it. This is probably a reasonable representation of what I really think and the algorithm is working nicely.

Here is what it generates in terms of daily death numbers:

There is still a bit of of tweaking to be done but I think this is going to be a better approach than the simple 3-chunk version I’ve been using up to now.

Thursday, June 04, 2020

Doubling times and the curious case of the dog that didn't bark

I'm sure many of you spend your mornings glued to Parliament TV but for anyone who missed it, there was a House of Lords Science and Technology Select Committee hearing on Tuesday, which discussed the science of Covid-19 with a particular focus on modelling. There are transcripts (two parts) here. and there is a video of the whole event here.

Two things particularly caught my attention. First bit is Question 31 in the transcript about 10:47ish when they are asked about the 20,000 death estimate that Vallance had hoped to undershoot. Matt Keeling first excused this by saying that the lockdown could have been stricter. However, Ferguson had already confirmed that the lockdown was in fact adhered to more strictly than the modellers had anticipated. (Keeling also went on to talk about the failure to protect care homes which is probably fair enough though the 20,000 target would have been badly overshot in any case). It would have surely been appropriate at this point to mention that the 20k target was already missed at the point that the lockdown was imposed on the 23rd March. By this point we had well over a million infected people, which already guaranteed well over 10k deaths, and there was no plausible way of stopping the outbreak dead in its tracks at that point. The decay was always going to be slower than the growth,  meaning more cases and deaths.

Then in the follow-up question, the witnesses were asked what had most surprised them compared to their earlier predictions. Ferguson answered with the biggest factor being that more infection had come from Spain and Italy and this is why they were further ahead on the curve than anticipated! The doubling time has nothing at all to do with how quickly people came with the infection, of course, and getting this wrong made a far greater difference to how things turned out.

It's curious that in two hours of discussion about the modelling, it never once came up that the modellers' estimate of the doubling time was 5-7 days up to the 18th march, and abruptly changed to 3-5 days on the 20th (reaching SAGE on the 23rd, according to the documentation)

Monday, June 01, 2020

The Dim and Dom show

I feel like I should be blogging about something related to the ongoing epidemic, but I can't bring myself to do it. The utterly vacuous, self-destructive, hopelessly incompetent nature of our government is beyond my ability to put into words. I am surprised at the scientists who are still prepared to work with them over the epidemic.

That aside, it's been an interesting couple of weeks. I'd been doing more of the same modelling and forecasting of the epidemic (and have updated our paper and submitted to a real journal), and then suddenly the media got hold of the delayed lockdown story. This is a very simple calculation, initially I thought too trivial to even write into a blog post, but it is of course very eye-catching. After mentions in the Guardian, Telegraph, More or Less, some requests for interviews came in. Initially I ducked them as I didn't really think it was appropriate for a non-expert to be pushing his own research especially as no-one else had backed it up at that point, and ATTP had tried to get results out of the IC model but initially came up with some significantly different answers (after a few more tries at getting the code to do the right things it worked very nicely though). Kit did a very good job on Sky I thought:

and then I found this manuscript (also written by an outsider, mathematical modeller like me) and the research showing essentially the same results for USA (manuscript here) (I think the smaller effect is mostly because they looked at a shorter interval) and also the Sunday Times article which managed to claim it was all new research from the IC team so I relented and did an interview for Vanessa Feltz on Radio 2 (which was live):

and also for the German channel ZDF which was recorded on Friday. Whether it will/did make the cut remains to be seen...the said they would send a link to the final version so I wait with bated breath.

Thursday, May 21, 2020 The EGU review

Well.. that was a very different EGU!

We were supposed to be in Vienna, but that was all cancelled a while back of course. I might have felt sorry for my AirBnB host but despite Austria banning everything they didn’t reply to my communication and refused a refund so when AirBnB eventually (after a lot of ducking and weaving) stepped in and over-ruled them and gave me my money back I didn’t have much sympathy. They weren’t our usual host, who was already full when I booked a bit late this year.

Rather than the easy option of just cancelling the meeting, the EGU decided to put everything on-line. They didn’t arrange videoconferencing sessions – I think this was probably partly due to the short notice, and also to make everything as simple and accessible as possible to people who might not have had great home broadband or the ability to use streaming software – but instead we had on-line chat (typing) sessions with presentation material previously uploaded by authors, that we could refer to as we liked. There was no formal division into posters and oral presentations. Authors could put up whatever they wanted (50MB max) onto the website beforehand and people were free to download and browse through at will. It is all still up there and available to all permanently, and you can comment on individual presentations up to the end of the month (assuming the authors have allowed this, which most seem to). The EGU has posted this blog with statistics of attendance which shows it to have been an impressive success.

Some people put up huge presentations, far more than they would have managed in a 15 minute slot, but most were more reasonable and presented a short summary. We did poster format for ours as we felt that this allowed more space for text explanation and an easier browsing experience than a sequence of slides with bullet points. Unfortunately my personal program of sessions I had decided to attend has been deleted from the system so I can’t review what I saw in much detail. I usually take notes but this time was too busy with computer screens.

Of course, being in Vienna in spirit, I had to have a schnitzel. I might have to have some more in the future, they were rather good and quite easy to make. Pork fillet, not veal.

The 2nd portion at the end of the week was better as I made my own breadcrumbs rather than using up some ancient panko that was skulking in the back of the cupboard. But we ate them too quickly to take pictures! Figlmüller eat your heart out!

The chat sessions were a bit frenetic. Mostly, the convenors invited each author in turn to post a few sentences in summary, following which there was a short Q-and-A free-for all. This only allowed for about 5 mins per presentation, which meant maybe 2 or 3 questions. But this wasn’t quite as bad as it seems since it was easy to scroll through the uploaded material ahead of time and pick out the interesting ones. Questioning could also run over subsequent presentations, it wasn’t too hard to keep track of who was asking what if you made the effort. As usual, there were only handful of interesting presentations per session for me (at most) so it was easy enough to focus on these. It was also possible to be in several different chat sessions at once, which you can’t do so easily with physical presentations! The structure made it more feasible to focus on whatever piqued our interest, and jules in particular spent more time at those sessions she does not usually get around to attending because they are outside of her main focus. Some convenors grouped presentations into themes and discussed 3-5 of them at a time, for longer. Some naughty convenors thought they would be clever and organise videoconferencing sessions outside of the EGU system, which actually worked pretty well in practice for those (probably a large majority to be honest) who could access it, but not so good for those who had access blocked for a number of reasons. Which is probably why the EGU didn’t organise this themselves. Whether it is actually preferable to the on-line chat is a matter of taste.

Jules was co-convening a couple of sessions and the convenors set up a small zoom session on the side to help coordinate, which added to the fun. A bit of personal chat with colleagues is an important aspect of these conferences. Her presentation is here and outlines some early steps in some work we are currently doing – an update to our previous estimate of the LGM climate, which is now getting on for 10 years (and two PMIP/CMIP cycles) old. I think we should probably find it encouraging that the new models don’t seem very different, though it may just mean that they share the same faults! There is some new data, perhaps not as much as we had hoped. And the method itself could do with a little bit of improvement.

I had actually found it a bit difficult to find the right session for my work when originally submitting it. It didn’t seem to quite fit anywhere, but in the end it turned out fine where I put it. The data assimilation stuff was a little less interesting methodologically speaking, perhaps because it’s a sufficiently mature field that everyone is just getting on with the nuts and bolts of doing it rather than inventing new approaches. I did get one idea out of it that I may end up using though, and this from the Japanese looks absolutely incredible from a technological point of view – nowcasting cloudbursts over Tokyo with a 30 second update cycle! With the extra year they’ve now got, it will probably be operational for the Olympics.

Jules and I also co-authored Martin’s work with us on emergent paleoconstraints which we were originally going to present for him as he wasn’t planning to attend. But, with the remote attendance he ended up able to do it himself which was a small bonus.

Best of all – no coffee queues! Well that and not needing to schlep out at 8pm looking for dinner each night…which is fun but gets pretty tiring by the end of the week. On the downside, we had to buy our own lunches rather than gatecrashing freebies all week like we usually (try to) do.

As for the future…well it seems pretty embarrassing that it took current events into forcing the EGU into moving on-line. Some of us have been pushing them on this for years and it’s always been met with “it’s too complicated” by the powers that be. I suspect they mostly like the idea of being in charge of a huge event and enjoy hobnobbing at all the free dinners (don’t we all!) but that doesn’t justify forcing everyone to fly over there and spend at least €2k minimum – probably rather more for most – to take part. It’s a huge amount of time, money, and carbon and we really ought to do better. If one good thing is to come out of the current mess, it might be that people finally wake up to the idea that working remotely really is fully feasible these days with the level of communication technology that is available. Blue Skies Research has been living your future life for more than 5 years now, and it’s great! Roll on next year. I know that turning up has added benefits, and don’t expect all travel to stop. But with remote access, people can easily “go” to both of the AGU and EGU each year, drop in to the bits that interest them, without having to devote a full week and more to each, with huge costs, jet-lag, the carbon budget of a small country, etc.

I expect that the AGU will want to put on a better show this December. Even if travel is opened up by then (which I wouldn’t be confident about at this point) I doubt this will happen quickly enough for the event to be organised in the usual manner. It will be good to have a bit of friendly rivalry to spur things on. In recent years, the AGU has generally been ahead of the EGU in terms of streaming and remote access – last December we watched a couple of live sessions and even asked a question (via text chat) though we were lucky that the small selection of streamed sessions included stuff of interest to us. The EGU has tended to put up streams of just a few of the public debate sessions rather than the science, and this only after the event with no opportunity for direct interaction. Bandwidth is a problem for streaming multiple sessions from the same location, but maybe even an audio stream with downloadable material would work? One thing is for sure, back to “business as usual” is not going to be acceptable now that they’ve shown it can be done differently.

Here’s Karlskirche which I hope to see again in the flesh some time.


Coincidentally, just a few days after the EGU I took part in this one-day webinar. It had a bit of the same sort of stuff – I presented the same work again, anyway! This was a zoom session which worked pretty well, there were one or two technical problems but you usually get in a real conference anyway with people plugging their laptops into the projector. It was great to have people from a range of countries attend and present at what would normally have been a local UK meeting of climathnet people. I have never quite managed to attend any of these before because they always seemed like a long way to travel for a short meeting that mostly isn’t directly relevant to our research. I expect to see a rapid expansion of remote meetings of various types in the future.