Estimating and Predicting Epidemic Behavior for the 2014 West African Ebola Outbreak

A Quick Stochastic Spatial SEIR Modeling Approach

Grant Brown
Jacob Oleson
Last Updated: 9/9/2014

Table of Contents

  1. Introduction
  2. Analysis 1: Polynomial Basis
  3. Analysis 2: Natural Spline Basis
  4. Conclusions

Introduction

About this document.

This is a living document, as the Ebola epidemic is rapidly evolving. Presented below is a preliminary modeling approach to the crisis, originally conceived as an illustration of the spatial SEIR model family generally and the capabilities of the rapidly developing (and totally unfinished) libspatialSEIR software library particularly. This is not yet peer reviewed research, or even a particularly complete analysis. Nevertheless, we hope that the initial exploration given below is instructive and useful.


Past versions of this analysis are cached in a repository on Github, and remain available:


Summary of Recent Changes
  • Sep 4 The situation continues to worsen, which is especially worrying given that the data used are almost certainly undercounts. Removed the polynomial analysis in favor of the spline basis analysis to save on computation time.
  • Aug 31 Basic reproductive number clearly very heterogenous between countries. Other interpretations TBD.
  • Aug 28 With the addition of Nigeria and new data, the situation appears dire across the board, especially in Liberia. The basic reproductive number calculations no longer look reasonable, so more must be done to estimate them separately for each country. The changes may also be due to recent changes to the R0 estimation code, which must now be reevaluated. This work is underway.
  • Aug 12 Liberia continues to worsen, though Sierra Leone appears to have leveled off. Guinea might be worsening slightly.
  • Aug 10 Predictions remain much the same, though perhaps not so immediately catestrophic as the August 5 predictions.
  • Aug 5 Models begin to show catastrophic predictions. Epidemic is changing too quickly for these simple models.
    • Shorten prediction window
    • Increase intensity process flexibility (quartic rather than cubic)
    • Recall that these models can't account for changing inteventions, of which there are many. The predictions are therefore made under the assumption that the epidemic will continue to evolve according to the same process it has so far.
  • Aug 3 Predicted epidmic curves start to strongly favor a uniformly worsening situation.
  • July 30 Previous predictions of mid-fall resolution of epidemic begin to shift.

The Outbreak

The 2014 Ebola outbreak in West Africa is an ongoing public health crisis, which has killed thousands of people so far. The cross-border nature of this epidemic, which emerged in Guinea, Liberia and Sierra Leone has complicated mitigation efforts, as has the poor health infrastructure in the region. This document explores a simple spatial SEIR model to make some initial predictions.

The Data

A summary of the WHO case reports is very helpfully compiled on wikipedia. It can be easily read into R with the xml library:

With data in hand, let's begin where every analysis should begin: graphs.

plot of chunk unnamed-chunk-2

These represent cumulative counts, but because case reports can be revised downward due to non-Ebola illnesses the graphs are not monotone. A quick, but effective solution to this problem is to simply "un-cumulate"* the data and bound it at zero to get a rough estimate of new case counts over time.

For better graphical representation, the "un-cumulated" counts are scaled to represent average number of infections per day, and linearly interpolated. The process is a bit noisier from this perspective when compared to the original cumulative counts.

plot of chunk unnamed-chunk-4

One can also represent this data geographically to get an idea of the spatial epidemic pattern, and to place the problem in a more relatable context.

Average Number of Infections Per Day:

Day:


Compartmental Models

Now that the data is read in (and now that we have several plots to suggest that we haven't done anything too terribly stupid with it) , let's do some compartmental epidemic modeling. Not only has Ebola been well modeled in the past using compartmental modeling techniques, but this author happens to be working on a software library designed to fit compartmental models in the spatial SEIRS family. What a strange coincidence! Specifically, we'll be using heirarchical Bayesian estimation methods to fit a spatial SEIR model to the data.



While a full treatment of this field of epidemic modeling is (far) beyond the scope of this writing, the basic idea is pretty intuitive. In order to come up with a simplified model of a disease process, discrete disease states (aka, compartments) are defined. The most common of these are S, E, I, and R which stand for:

  • Susceptible to a particular disease
  • Exposed and infected, but not yet infectious
  • Infectious and capable of transmitting the disease
  • Removed or recovered

This sequence, traversed by members of a population (S to E to I to R), forms what we might call the temporal process model of our analysis. This analysis belongs to the stochastic branch of the compartmental modeling family, which has its roots in deterministic systems of ordinary and partial differential equations. In the stochastic framework, transitions between the compartments occur according to unknown probabilities. It is the S to E probability, which captures infection activity, into which we introduce spatial structure. Some details of this are given as comments to the code below, and more information than you probably want on the statistical particulars is available in this pdf document. For now, suffice it to say that we'll place a simple spatial structure on the epidemic process which simply allows disease to spread between the three nations involved, and we'll try to estimate the strength of that relationship. Many other potential structures are possible, limited primarily by the amount of additional research and data compilation one is willing to do.



For the purposes of this analysis, we will not do anything fancy with demographic information or public health intervention dates. Demographic parameters are relatively difficult to estimate here, as there are only four spatial units which are all from the same region. Intervention dates are more promising, but their inclusion requires much more background research than we have time for here. In the interest of simplicity and estimability, we'll just fit a different disease intensity parameter for each of the three countries to capture aggregate differences in Ebola susceptibility in addition to using a set of basis functions to capture the temporal trend.

Analysis 1
A polynomial basis analysis is available in previous versions of this document. It was removed (9/4/2014) in favor of the second, spline based, analysis to save on computation time, as polynomials generally extrapolate poorly. This makes them generally unreliable for prediction purposes.


Analysis 2: Spline Basis

Set Up

There are some things we need to define before we can start fitting models and making predictions.

  1. The population sizes need to be determined.
  2. Initial values for the four compartments must be determined.
  3. The time points are not evenly spaced, so we need to define appropriate offset values to capture the amount of aggregation performed (time between reports).
  4. We must define the spatial correlation structure.
  5. A set of basis functions needs to be chosen to capture the temporal trend.
  6. Prior parameters and parameter staring values must be specified for each chain.
  7. A whole bunch of bookkeeping stuff for which I haven't yet programmed sensible default behavior needs to be set up.

Compartment starting values follow the usual convention of letting the entire initial population be divided into susceptibles and infectious individuals. The starting value for the number of infectious individuals was 86, the first infection count available in the data. Offsets are actually calculated in the first code block (above) as the differences between the report times. For temporal basis functions, natural splines were used. Prior parameters for the E to I and I to R transitions were chosen based on well documented values for the average latent and infectious times, and the rest of the prior parameters were left vague. These decisions are addressed in more detail as comments to the code below.



With the set up out of the way, we can finally build and run the models the models. The code presented below has recently (8/27/2014) been set up to use the "parallel" library which is included with the R statistical analysis software. While this allows us to spend considerably less time waiting for the document to compile, the code may be more difficult to understand for those unfamiliar with this useful parallelization library. The previous analyses, linked above, may be a better guide to the basic use of the spatialSEIR library for such users.




Convergence

As this is a Bayesian analysis in which the posterior distribution is sampled using MCMC techniques, we really need some indication that the samplers have indeed converged to the posterior distribution in order to make any inferences about the problem at hand. In the code below, we'll read in the MCMC output files created so far, plot the three chains for each of several important parameters, and take a look at the Gelman and Rubin convergence diagnostic (which should be close to 1 if the chains have converged.)

plot of chunk unnamed-chunk-9
plot of chunk unnamed-chunk-10

Basic Reproductive Number Calculation

A common tool for describing the evolution of an epidemic is a quantity known as the basic reproductive numer, the basic reproductive ratio, or one of several other variants on that theme. The basic idea is to quantify how many secondary infections a single infectious individual is expected to cause in a large, fully susceptible population. Naturally, when this ratio exceeds one we expect the epidemic to spread. Conversely, a basic reproductive number less than one indicates that a pathogen is more likely to die out.



plot of chunk unnamed-chunk-11
While the basic reproductive number is a useful quantity to know, it does not directly make any predictions about future epidemic behavior. In order to do that, we need to simulate epidemics based on the MCMC samples we have obtained and summarize their variability over time.

Epidemic Prediction




Below, we will attempt to predict the course of the epidemic through early fall. We must be cautious when making predictions about a chaotic process this far into the future. The intensity function which drives the exposure process is based on simple smooth functions of time, rather than any external information. While this is perfectly adequate for estimation, it may or may not provide good prediction performance.

plot of chunk unnamed-chunk-12

It can also be helpful to take a look at the data in tabular form.


Estimated and Predicted number of Infectious Individuals


   Guinea         Liberia         Sierra Leone         Nigeria    
 2014-03-24  0 0 0 0
 2014-03-26  36 0 0 0
 2014-03-28  51 8 0 0
 2014-03-31  56 7 0 0
 2014-04-07  62 12 0 0
 2014-04-10  71 21 0 0
 2014-04-14  67 25 0 0
 2014-04-17  94 24 0 0
 2014-04-21  83 22 0 0
 2014-04-24  72 27 0 0
 2014-05-01  41 25 0 0
 2014-05-03  39 25 0 0
 2014-05-10  39 24 0 0
 2014-05-18  42 18 0 0
 2014-05-27  30 12 0 0
 2014-05-29  56 10 16 0
 2014-06-03  80 7 68 0
 2014-06-06  86 5 55 0
 2014-06-15  86 4 50 0
 2014-06-17  87 23 42 0
 2014-06-22  71 18 97 0
 2014-07-02  65 79 138 0
 2014-07-08  47 79 163 0
 2014-07-14  35 96 193 0
 2014-07-20  31 97 202 0
 2014-07-27  39 119 228 0
 2014-08-01  67 211 204 3
 2014-08-06  72 276 263 8
 2014-08-11  67 288 236 10
 2014-08-16  66 412 265 8
 2014-08-20  108 497 295 9
 2014-08-31  132 702 311 7
 2014-09-05  224 857 428 10
 2014-09-07  246 953 442 17
 2014-09-09  241 951 429 30
 2014-09-11  223 916 397 54
 2014-09-13  214 942 381 95
 2014-09-15  209 997 369 140
 2014-09-17  207 1089 367 183
 2014-09-19  210 1210 368 222
 2014-09-21  215 1375 372 262
 2014-09-23  223 1580 384 305
 2014-09-25  236 1848 398 356
 2014-09-27  253 2183 420 418
 2014-09-29  276 2615 447 494
 2014-10-01  303 3165 482 590
 2014-10-03  342 3874 526 715
 2014-10-05  390 4790 583 876
 2014-10-07  454 5979 652 1088
 2014-10-09  535 7542 745 1363
 2014-10-11  642 9598 861 1734
 2014-10-13  787 12322 1011 2229
 2014-10-15  979 15902 1211 2894
 2014-10-17  1243 20581 1473 3802
 2014-10-19  1598 26514 1831 5029
 2014-10-21  2089 33731 2304 6674
 2014-10-23  2758 41926 2947 8831
 2014-10-25  3669 50503 3798 11494
 2014-10-27  4875 58733 4944 14612
 2014-10-29  6460 65864 6386 17913
 2014-10-31  8441 71380 8215 21134
 2014-11-02  10901 74835 10337 23950
 2014-11-04  13751 76240 12756 26244
 2014-11-06  16976 75675 15225 27910
 2014-11-08  20294 73340 17721 29025
 2014-11-10  23596 69543 19993 29615
 2014-11-12  26579 64542 22113 29774
 2014-11-14  29215 58829 23814 29564
 2014-11-16  31341 52656 25073 29016
 2014-11-18  32869 46488 25685 28222
 2014-11-20  33662 40465 25766 27185
 2014-11-22  33694 34831 25278 26010
 2014-11-24  33064 29642 24408 24699
 2014-11-26  31864 24988 23208 23337
 2014-11-28  30279 20850 21813 21882
 2014-11-30  28391 17253 20254 20422
 2014-12-02  26375 14133 18658 18895

Such data can also be visualized in map form, though the recent surge in predicted cases has the effect of swamping the earlier dynamics:


Total Infection Size - Estimated and Predicted:

Day:



Conclusions

This epidemic is evolving extremely rapidly. As of 8/12, it looked like the situation in Liberia was set to continue worsening, and that Guinea is at risk of the same (though not to nearly the same degree). On the other hand, the epidemic in Sierra Leone appeared to be leveling off (though not disappearing). As of 8/28, these look like reasonable predictions, however the models appear to have resumed predicting a fairly catastrophic continued spread, especially in liberia. In particular, the models predict that the epidemic will take off in Nigeria, as the countries are assumed in this case to share several intensity parameters. We may hope that this particular simplifying assumption is invalid, however it is not a hopeful sign that WHO predictions are also becoming catastrophic. These models can not anticipate public health interventions and sudden changes in governmental policy and individual behavior, but recent news from the region gives little reason to hope for a swift end to the epidemic. It is more important now than ever to support the efforts of involved governmental and non-governmental organizations like the WHO and MSF


That wraps up the analyses for now. This document will continue to be updated as the epidemic progresses, reflecting new data and perhaps additional analysis techniques. As the document is tracked via source control it will be easy to see how well past predictions held up and how they change in response to new information. Questions and comments can be shared here