Data

CDC [16] compiled data from three countries of West Africa, based on weekly Ebola reports by the WHO (Fig. 1) that counted the confirmed, probable and suspected infections and fatalities in each country. CDC [16] added the numbers of infections since the start of the outbreak to the number of total cases.

Fig. 1. Total weekly and monthly count of Ebola cases in West Africa; blue dots connected with a black line are the values of Table 1 and brown rings are the counts from CDC [16]; plotted using MS Excel.

Unfortunately, there were some flaws in the data. First, there were false-positive and false-negative diagnoses: initially Ebola was diagnosed as fatal diarrhoea (December 2013 and January 2014) and at the peak of the outbreak villages reported Ebola cases that were later falsified, but the total counts in the reports were not corrected retrospectively. Second, the published data contained typos, incorrect arithmetic and an obvious outlier at 24 June of 2015. (We suppose that it resulted from a typo: 24 472 instead of 27 472.) We therefore aggregated the CDC list to Table 1, which eliminated data uncertainty as much as possible without altering the data. It also reduced random fluctuations. Nevertheless, it still was representative of the raw data (Fig. 1).

Table 1. Total monthly count of Ebola cases in West Africa

^a For each month, the table records the maximum of the total number of cases reported from Guinea, Liberia and Sierra Leone since the beginning of the outbreak.

^b Month 0 is December 2013 (index patient).

Source: Adapted from [16].

In addition, we also fitted the BP models to the full set of data from CDC [16] about West Africa. As amongst the 265 data-points (brown rings in Fig. 1) there were contradictory values for 19 April, 3 May and 30 June of 2015, we took the larger counts (reducing the size of the dataset to 262 data-points) and left the data otherwise unaltered (not removing the outlier).

The BP model class

The BP differential equation (1) describes the cumulative count of infected individuals, y(t), as a function of time, t, since the start of an infection. It uses five model parameters that are determined from fitting the model to cumulative epidemiological data; non-negative exponent pair a < b, non-negative scaling constants p and q and the initial value y(0) = c > 0:

(1)$$y^{\prime}\,\lpar t\rpar = p\cdot y\,\lpar t\rpar ^a\,-\,q\cdot y\,\lpar t\rpar ^b$$

Equation (1) appeared originally in Pütter [Reference Pütter17]; von Bertalanffy [Reference Bertalanffy18, Reference Bertalanffy19] popularised it as a model for animal growth. In epidemiology, and for a ≤ 1, this model has been proposed as the generalised Richards model [20, 21]. It can be solved analytically, although in general not by means of elementary functions [Reference Ohnishi, Yamakawa and Akamine22].

We interpret equation (1) as a model class, where each exponent pair (a, b) defines a unique model BP(a, b) in the BP class of models; BP(a, b) has three free parameters (c, p, q). Equation (1) then includes several trend models that have been used to describe epidemic trajectories [Reference Bürger, Chowell and Lara-Diaz23–Reference Zhao25], such as the Brody [Reference Brody26] model of bounded exponential growth BP(0, 1), [Reference Verhulst27] logistic growth BP(1, 2), the model BP(2/3, 1) of von Bertalanffy [Reference Bertalanffy18] or the model BP(3/4, 1) of West et al. [Reference West, Brown and Enquist28]. Also, the Gompertz [Reference Gompertz29] model fits into this scheme: it is the limit case BP(1, 1), with a different differential equation, where b converges to a = 1 from above [Reference Marusic and Bajzer30]. Equation (1) also includes several classes of models, such as the generalised Bertalanffy model (b = 1 with a variable) and the Richards [Reference Richards31] model (a = 1 with b variable). However (Fig. 2), when compared to the range of models searched for this paper (yellow area), these named models appear as rather exceptional.

Fig. 2. Named models (blue) and initial search region (yellow) of BP models (plot using Mathematica 12.0): 0 ≤ a ≤ 1.3, a < b ≤ a + 3, step size in both directions 0.01. When needed, the grid was extended.

Some authors added further assumptions about the parameter values in order to simplify the model (e.g. setting c = 1 in advance for the indicator case). We do not consider such simplifications.

Model calibration

The (ordinary) least-squares method is common in epidemiology for large case counts, as for the present data [Reference Smirnova and Chowell32]. It measures the goodness of the fit to the data by means of SSE, the sum of squared errors (fit residuals). In the present context, we sought parameters a, b, c, p and q so that for the solution y(t) of equation (1) the following sum SSE was minimised, whereby (t_i, y_i) were the first n total case counts y_i at time t_i (here 10 ≤ n ≤ 28):

(2)$$\;SSE = \mathop \sum \limits_{i = 1}^n \lpar {y_i-y\,\lpar t_i\rpar } \rpar ^2$$

An implicit assumption of the least-squares method is a normal distribution of the errors: in formal terms, SSE is equivalent to the maximum-likelihood estimate under the assumption of normally distributed data with expected value y(t) and time-independent variance s > 0. If this assumption is not satisfied, as for the animal size-at-age data or the outbreak data, other methods of calibration need to be considered. For the animal data we could work well under the assumption of a log-normal distribution, where the variance increases with the size. However, for the Ebola data (total case counts) this assumption was problematic, too, as the variance was maximal at an intermediate stage (epidemic peak).

We therefore developed another measure for the goodness of fit, weighted least-squares that aimed at finding parameters to minimise the following sum of weighted squared errors SWSE:

(3)$$SWSE = \sum\limits_{i = 1}^n {} \displaystyle{{(y_i-y{{\rm (}t_i{\rm ))}}^2} \over {\vert y^{\prime}{\rm (}t_i{\rm )}\vert}}$$

This weight function tolerates a higher variability in the total count during the epidemic peak. In formal terms, SWSE identifies the maximum-likelihood estimate under the assumptions that the data are normally distributed with expected value y(t) and a variance s⋅y′(t) with s a constant. (Thus, the variance of the total counts is assumed to be proportional to the absolute value of the derivative of y. This derivative corresponds to the number of new cases.)

In the following subsections, we will explain our notation using SSE; the same definitions will also apply to SWSE.

Optimisation

For model class (1) standard optimisation tools (e.g. Mathematica: NonLinearModelFit) may not identify the best-fit parameters (numerical instability). To overcome this difficulty, we used a custom-made variant of the method of simulated annealing [Reference Vidal33] which solves the optimisation problem to a prescribed accuracy [Reference Kühleitner12]. Note that using simulated annealing to optimise all five parameters of equation (1) at once may not always come close to the optimum parameters, because the region of nearly optimal parameters may have a peculiar pancake-like shape ([Reference Renner-Martin15]; extremely small in one direction and extremely large in other ones).

We confined optimisation to a grid (Fig. 2). For each exponent pair on the grid 0 ≤ a ≤ 1.3 and a < b ≤ 3 with distance 0.01 in both directions we searched for the best-fit model BP(a, b); i.e. we minimised SSE for the model BP(a, b). Thus, we defined the following function SSE_opt on the grid:

(4)$$SSE_{opt}\lpar {a\comma \;b} \rpar = \mathop {\min }\limits_{c\comma p\comma q} \lpar {SSE} \rpar \comma \;\;\;\quad SSE_{min} = \mathop {\min }\limits_{a\comma b} SSE_{opt}\lpar {a\comma \;b} \rpar $$

Summarising this notation, for each exponent pair (a, b) on the grid we identified the best-fit model BP(a, b) with certain parameters c, p, q and the corresponding least sum of squared errors SSE_opt(a, b). These values were recorded in a spreadsheet. Finally, the best-fit model had the overall least sum of squared errors (SSE_min). It minimised the function SSE_opt at the optimal exponent pair (a_min, b_min). From the above-mentioned spreadsheet we could then read off the other parameters p_min, q_min, c_min that minimised SSE_opt(a_min, b_min) of the model BP(a_min, b_min). However, if that optimal exponent pair was on the edge of the search grid, then we extended the search grid and continued the optimisation. (Otherwise, if the optimal exponent pair was surrounded by suboptimal ones, then we stopped the search for an optimum.)

Model comparison

For each exponent pair (a, b) of the search grid we identified a best-fitting model BP(a, b), whose fit was given by SSE_opt(a, b) of formula (4). Amongst these models we selected the model with the least SSE_opt. In the literature, there are various alternative criteria for the comparison of models, amongst them the root mean squared error RMSE and the Akaike's [Reference Akaike34] information criterion (AIC). RMSE is the square-root of SSE/n and AIC = n ⋅ ln (SSE/n) + 2K, where n is the number of data-points (between 10 and 28 for the monthly truncated data), K = 4 is the number of optimised parameters of the model BP(a, b) with a given exponent pair (counting c, p, q and SSE), and SSE = SSE_opt(a, b). As an additional measure to compare the goodness of fit, we used formula (5) for the relative Akaike weight:

(5)$$\wp = \displaystyle{{{\rm exp}\lpar {-\lpar \lpar AIC-AIC_{min}\rpar /2\rpar } \rpar } \over {1 + {\rm exp}\lpar {-\lpar \lpar AIC-AIC_{min}\rpar /2\rpar } \rpar }}$$

This formula has been interpreted as probability, $\wp$, that the (worse) model with higher AIC would be ‘true’, when compared to the model with the least AIC [Reference Burnham and Anderson35, Reference Motulsky and Christopoulos36]. Thereby (citation from [Reference Burnham and Anderson37], on p. 272, referring to information theory of [Reference Kullback and Leibler38]) ‘true’ means that the model ‘is, in fact, the K-L best model for the data… This latter inference about model selection uncertainty is conditional on both the data and the full set of a priori models [in formula 5: the models with AIC and with AIC_min] considered’. Furthermore, as does AIC, the Akaike weight assumes a normal distribution of the (weighted) fit residuals.

Using SSE or AIC is only meaningful with reference to a fixed dataset (different models are fitted to the same data), whereas RMSE and $\wp$ make also sense, if we consider datasets with different sizes, as for the forecasts, where we successively use truncated initial data with 10, 11, …, 28 months. Thereby, we prefer the probabilities for their more intuitive appearance. Note that $\wp \equals 50\percnt$ is the maximal value that can be attained by $\wp$ (comparison of the best-fit model with its copy).

Multi-model approach

For each exponent pair (a, b) of the search grid we identified the best-fit growth curve y_a,b(t) for the model BP(a, b); its parameters p, q and c were optimised according to formula (4). Thereby, we could observe a high variability in the best-fit exponents: even for exponent pairs (a, b) with notable differences from the optimal exponent pair (a_min, b_min) the best-fit growth curve y_a,b(t) barely differed from the data. To describe this phenomenon, we used the following terminology: a model BP(a, b) and its exponent pair (a, b) are u-near-optimal with model-uncertainty u, if SSE _opt(a, b) ≤ (1 + u) ⋅ SSE _min. Figure 3 illustrates this concept by plotting near-optimal exponent pairs for different levels u = 0.17 and u = 0.5 of model-uncertainty (red and blue dots).

Fig. 3. Optimal exponent-pair (black) for fitting the data of Table 1 with respect to SSE; red and black dots for 245 and 1330 exponent pairs with up to 17% and 50% higher SSE; exponent-pairs of the Bertalanffy, Gompertz and Verhulst models (cyan) and part of the search grid (yellow).

Using formula (5), we translated the level u of near optimality into a probability $\wp = 1/ \lpar 1 + \sqrt {{\lpar {1 + u} \rpar }^n} \rpar$. For example, if we wish to consider models (exponent pairs) with a probability of at least 10%, this corresponds to a model uncertainty of at most u = 0.55 for the 10-month data (n = 10) and u = 0.17 for the 28-month data.

For each dataset, we then may ask about the forecasts that would be supported by models with e.g. $\wp \equals 10\percnt$ or higher, discarding other possible forecasts as unlikely, given the data. Thus, given a model probability, $\wp$, we studied the forecasts that came from the ensemble of the growth curves y_a,b(t) corresponding to the near-optimal models BP(a, b) with that probability or higher. For example, for each of these y_a,b(t) we estimated its final disease size by the asymptotic limit. The upper and lower bounds for these estimates defined a prediction interval that informed about the likely disease size. Note that this is an analysis of model uncertainty, not of data uncertainty, whence the prediction interval is not a confidence interval. (Confidence intervals assess data uncertainty by means of simulations that add random errors to the data and use best-fit functions to compute bounds for the confidence interval. For prediction intervals, different models are fitted to the same data.)