A primer on stable parameter estimation and forecasting in epidemiology by a problem-oriented regularized least squares algorithm.

Public health officials are increasingly recognizing the need to develop disease-forecasting systems to respond to epidemic and pandemic outbreaks. For instance, simple epidemic models relying on a small number of parameters can play an important role in characterizing epidemic growth and generating short-term epidemic forecasts. In the absence of reliable information about transmission mechanisms of emerging infectious diseases, phenomenological models are useful to characterize epidemic growth patterns without the need to explicitly model transmission mechanisms and the natural history of the disease. In this article, our goal is to discuss and illustrate the role of regularization methods for estimating parameters and generating disease forecasts using the generalized Richards model in the context of the 2014-15 Ebola epidemic in West Africa.

Introduction

Developing tools for stable parameter estimation and reliable forecasting of emerging and re-emerging infectious disease epidemics represents a key priority for public health officials and government agencies in their work to prevent and mitigate disease threats. At the early stage of an outbreak, when incidence data are limited and subject to reporting delays, it is often premature to characterize the transition rates of individuals between various disease epidemiological compartments using, for instance, SEIR-type systems of differential equations, which may involve a substantial number of unknown parameters. For the early transmission period, phenomenological models of a logistic type, describing the progression of the epidemic in terms of the cumulative number of reported cases, C, provide a simple alternative. In what follows, we employ the generalized Richards model (Chowell et al., 2016, Smirnova et al., 2017, Turner et al., 1976)to estimate crucial disease parameters, such as the intrinsic growth rate (r), the deceleration of growth (p), the final size of the epidemic (K), the disease turning point (τ), and the extent of deviation from the S-shaped dynamics of the classical logistic-growth curve (a), see Fig. 1. In (1.1), T is the duration of the outbreak.

Fig. 1

Parameter identification process from early incidence data using the generalized Richards model.

For the original Richards model (), the closed-form solution is availablewhich simplifies its theoretical and numerical study. However, the statistical analysis of the fitting accuracy for the original and generalized Richards models indicates that the generalized Richards model outperforms its original version for epidemics characterized by an early sub-exponential growth phase (Chowell et al., 2016). For such outbreaks, the value of parameter p in (1.1) is less than 1 (), and the reduction in the residual sum of squares appears to be statistically significant. Several mechanisms could give rise to initial sub-exponential growth in case incidence including (Chowell et al., 2015, Viboud et al., 2015): (i) spatially clustered contact structures (e.g., high clustering levels) (Chowell et al., 2015, Szendroi and Csanyi, 2004); (ii) early onset of population behavioral changes and control interventions (Chowell et al., 2015, Szendroi and Csanyi, 2004); and (iii) substantial heterogeneity in susceptibility and infectivity of the host population that introduces high local variability in the local reproduction number that fluctuates around the epidemic threshold at 1.0. The two limiting cases of the generalized Richards model are those in which (constant incidence growth) or (exponential growth) (Viboud et al., 2015).

Fig. 2 illustrates representative epidemic profiles, , that the generalized Richards model supports, as the deceleration of growth parameter, p, is varied. Overall, as parameter p increases, the epidemic profile shows faster growth that converges to exponential rate as while keeping the epidemic size, K, fixed. On the other hand, parameter a modulates the epidemic turning point (Fig. 2), which occurs earlier as this parameter increases.

Fig. 2

Simulations from the Generalized Richards Model (A) and the Original Richards Model (B). A) Parameter p is varied from 0.5 to 1.0 while keeping parameter a fixed at 1.0. B) Parameter a is varied from 0.3 to 2.0 while keeping parameter p fixed at 1.0. Parameter r is fixed at 0.3 per day, , and the initial number of cases, .

Thus, in our investigation we utilize the generalized four-parametric model (1.1), and solve the ODE-constrained least squares problem to estimate r, p, K, τ, and a from early data of the 2014–15 Ebola epidemics in Guinea, Sierra Leone, and Liberia. We then use the reconstructed parameter values to forecast future incidence cases by propagating uncertainty in the system forward in time.

The regularized least squares problem

To present a regularized numerical algorithm for stable parameter estimation, we introduceand arrive at the normalized model (Cavallini, 1993)to be fitted to early (noisy) incidence data, , which approximate the vector , and is defined by (1.1). This yields the following least squares problem (LSP)

In (2.3), the function is evaluated at finitely many points, …, and is the Euclidian norm in . By the first order necessary condition for the minimum, one concludes that , and thereforewhere is the inner product in . Substituting (2.4) into (2.3), one gets the LSP with respect to the parameters b, p, a, and τ:

Minimization problem (2.5) is ill-posed due to the lack of stability with respect to noise in the input data. To overcome the danger of noise propagation, which is evident from the severely ill-conditioned Jacobian, our next step is

to regularize (stabilize) the functional in (2.5) by incorporating a priori information about the unknown parameters;

to develop an efficient problem oriented numerical solver that would further increase our chances of computing a reasonably accurate solution.

Without regularization, the algorithm may diverge due to error accumulation at every step, or the uncertainty could get out of control. The most common a priori information is the one that uses our best guess on the potential values of . Less common, but very helpful for this particular application, is the information on the expected orders of magnitude for b, p, a, and τ. Specifically, since τ has the meaning of the disease turning point, its value should be about one or two orders of magnitude larger than b, p, and a. This suggests that, as we introduce a priori information, parameters b, p, and a must be appropriately weighted to balance the sensitivity of our model to all four unknowns. Combining all of the above, we add the following term to the cost functional (2.5)

Here incorporates the expected (reference) value of , and L is the weighting matrix. Since the information contained in is less reliable than the one contained in the fidelity term, , we penalize by a small regularization parameter, α, and obtain the following auxiliary minimization problem (Tikhonov et al., 1995, Tikhonov et al., 1998, Yagola, 2010, pp. 17–34)which is a particular case of the variational (Tikhonov's) regularization aimed at stabilizing the original ill-posed model. In (2.7),

Numerical optimization method with Hessian reduction

Since defined in (2.8) is a non-linear operator, the regularization needs to be combined with a robust problem-oriented numerical algorithm capable of producing an accurate approximation to the unknown parameters. If one solves (2.7) by a classical Gauss-Newton method and drives α to 0 as iterations progress, then one gets what is known as iteratively regularized Gauss-Newton (IRGN) scheme (Bakushinsky, 1993, Bakushinsky and Kokurin, 2005)where is the adjoint to . In (3.1), and execute regularization and line search, respectively. While this method is generally efficient, for our model one is forced to sacrifice too much accuracy to stabilize the process, and even with a relatively large regularization parameter the uncertainty in the computed solution remains high. Hence our goal is to modify IRGN (3.1) by truncating a “non-essential” part of the Hessian approximation in order to bring down its condition number. To that end, consider the gradient of f, where

One can easily verify that the residual, , is in the kernel of the matrix . Indeed, from definitions of in (2.4) and in (2.8), one concludes that for any

Therefore and . This implies the Hessian approximationwhich we reduce further by dropping and setting This permits to lower the condition number and to enforce the symmetric non-negative definite structure. Dividing both sides of the modified linear system by , which enables us to move all noise from the matrix to the right-hand side, one arrives at what we call Reduced IRGN (RIRGN) (Smirnova et al., 2017)where , and

Simulation results and discussion

In this section we address the important aspects of practical implementation of RIRGN algorithm (3.4). To that end, we investigate the trajectories of the 2014–2015 Ebola epidemic in Guinea, Liberia and Sierra Leone (World Health Organization, 2016). For illustration, we use early incident reported data to estimate the values of five unknown parameters, r, τ, p, a, and K, by solving the regularized least squares problem with the RIRGN method (3.4). Then we substitute the estimated parameters into (1.1) to find the approximate values of future incident cases using the built-in Matlab ode23s solver. The data for this numerical experiment have been retrieved in 2015 from the World Health Organization (WHO) web site. The weekly series of cases were directly obtained from the WHO patient database.

For successful application of any regularized procedure, the right choice of the stabilizing sequence, , is very important. The value of has to be sufficiently large to ensure that we incorporate enough stability into the system to reduce the excessive noise propagation. On the other hand, over-regularization may slow down the convergence of the iterative sequence , and result in unnecessary loss of precision. In our experiment, we take for Sierra Leone and Liberia and for Guinea, where the data are more noise contaminated. We set , since this choice gives us the most aggressive convergence rate for while maintaining stability for every value of k.

The outer step of the RIRGN scheme (3.4) identifies the direction of the decent. In order to optimize the step size in that direction, we implement a version of the Armijo-Goldstein line search strategy (Nocedal & Wright, 2000), i.e., a backtracking with untilwhich is commonly used for Gauss-Newton type algorithms. Iterations are stopped at the moment when the discrepancy for the approximate solution is consistent with the level of noise in our data in order to guarantee convergence and, at the same time, it is not too small to prevent over-fitting that may compromise the accuracy of the estimated parameters, i.e., the stopping index, , is evaluated by a posteriori stopping rule

To balance the sensitivity of the cost functional to all four unknowns in the penalty term, we choose L based on our prior knowledge of the difference in the orders of magnitude between the value of τ and the values of three other parameters. This suggests L in the form (2.6) with ( works well for our code).

The selection of the initial vector, , and the reference vector, , in the penalty term is rather difficult. In the beginning of an emerging outbreak, it is not straightforward to make a sophisticated guess about when the disease is going to peak. Hence, in our experiment, we take in to be slightly greater than , the last week for which the data is available, and we take τ in to be equal to 60, i.e., we assume the incidence curve will turn before 60 weeks of the epidemic have passed. The values of p in and are taken to be 0.7 and 0.1, respectively, to enforce the sub-exponential growth rate. The values of two remaining parameters, b (used to compute r) and a are taken to be 1 in both vectors for the lack of better information.

Fig. 3 illustrates parameter estimation and the epidemic forecasts based on the generalized Richards model with 20 weeks of incident case data for the Ebola epidemic in Sierra Leone. After the five unknown parameters have been recovered from early epidemic data, 100 additional data bootstrap curves are generated by adding a Poisson error structure to the weekly series of reported cases (see (Chowell, Ammon, Hengartner, & Hyman, 2006)) in order to quantify uncertainty in the recovered parameters. An alternative is to employ analytic results from asymptotic theory to estimate parameter uncertainty (Banks, Davidian, Samuels, & Sutton, 2009). The histograms in the upper row show the recovered values of the parameters, while the collection of curves at the bottom of the figure demonstrates the accuracy of the forecasting.

Fig. 3

Sierra Leone: Estimation of r, τ, p, a, and K, and Forecasting. r = 0.52 (95% CI: 0.31, 0.71), τ = 27 (95% CI: 25, 29), p = 0.86 (95% CI: 0.83, 0.91), a = 1.03 (95% CI: 1.02, 1.04), K = 1.7 (95% CI: 1.5, 1.8).

In a similar manner, we estimated parameters and produced forecasts for the Ebola epidemics in Guinea (Fig. 4) and Liberia (Fig. 5). Not surprisingly the model has most difficulties tracking the epidemic in Guinea due to its multi-mode nature (Fig. 4). Since this outbreak is also the longest, we need more data (30 weeks) to forecast future incidence cases. On the other hand, the Liberia outbreak, which is 45 weeks total, can be projected from as little as 12 weeks of early data. This inverse modeling approach could be adapted to analyze outbreaks of other disease system, e.g., avian influenza, Middle East Respiratory Syndrome (MERS), and Zika. One could envision disease forecasts for an emerging outbreak in real time with updates every few weeks or so (depending on the disease system).

Fig. 4

Guinea: Estimation of r, τ, p, a, and K, and Forecasting. r = 1.3 (95% CI: 0.03, 2.7), τ = 41 (95% CI: 33, 49), p = 0.87 (95% CI: 0.7, 0.9), a = 0.96 (95% CI: 0.95, 1), K = 1.1 (95% CI: 1, 5.9).

Fig. 5

Liberia: Estimation of r, τ, p, a, and K, and Forecasting. r = 0.12(95% CI: 0.63, 1.9), τ = 15(95% CI: 13, 17), p = 0.82 (95% CI: 0.7, 0.88), a = 1 (95% CI: 1, 1.04), K = 3.5 (95% CI: 2.2, 4.4).

In conclusion, we have illustrated the potential of rigorous mathematical approaches for generating reliable parameter estimates and useful epidemic forecasts in the context of the generalized Richards equation, a simple phenomenological 4-parameter model. Our results here suggest that carefully linking mathematical models with regularization techniques could lead to stable parameter estimation and associated forward projections.