Gompertz curve

Our model is based on the sigmoidal growth curve pattern that characterizes epidemics. We start by assuming that the cumulative number of cases follows a Gompertz curve, which is a parsimonious model for the canonical sigmoid shape of cumulative case numbers in a one-wave epidemic. Over the course of a wave, the number of new infected cases increases up to a peak before declining to zero as the pool of susceptible individuals declines. Specifically, if the cumulative number of cases at time t, μ(t), follows a Gompertz curve, we can write μ(t) = μ̄exp {γ₀e^γt}, where μ̄ is the unknown saturation level for the cumulative number of cases, γ₀ < 0 is a parameter related to μ(0) and γ < 0 is the growth rate parameter. Defining μ̇(t) ≡ dμ(t)/dt and g(t) ≡ μ̇(t)/μ(t), it is straightforward to show that

where δ = ln γ₀γ.

The observational model needs to be specified in discrete, rather than continuous, time. This is straightforward. Let Y_t be the observed cumulative number of cases on day t and y_t = Y_t − Y_t − 1 be the number of daily new cases.¹ We can then define the growth rate of Y_t as g_t = y_t/Y_t − 1 and replace ln g(t) with ln y_t − ln Y_t − 1.

Dynamic Gompertz model

The deterministic trend implied by (1) is too inflexible for practical time-series modeling of an epidemic. Replacing it with a stochastic trend allows the model to adapt to changes in dynamics during the course of the epidemic. We call this stochastic-trend counterpart of (1) the dynamic Gompertz model. It is a local linear trend model specified as

where ln g_t = ln y_t − ln Y_t − 1 and where the disturbances ε_t and ζ_t are mutually independent, and NID(0, σ²) denotes normally and independently distributed with mean zero and variance σ². Note that the larger the signal-to-noise ratio, q_ζ = σ_ζ²/σ_ε², the faster the estimate of the slope parameter, γ_t, which can be interpreted as the growth rate of the growth rate of cumulative cases, changes in response to new observations. Conversely, a lower signal-to-noise ratio induces more smoothness to the estimates. When σ_ζ² = 0, the trend is deterministic as in (1).

State space form and estimation

It is convenient to write the dynamic Gompertz model in general state space form: with This model can be estimated using techniques based on the Kalman filter once a prior is specified. The prior is

where a₁ is a 2 × 1 vector of prior means and P₁ a 2 × 2 prior variance matrix. We use a diffuse prior due to the absence of prior information about the epidemic when the model is first estimated: i.e., we set a₁ = (0, 0)′, P₁ = κI, and let κ → ∞. Model estimation, including implementation of the diffuse prior, is carried out using the KFAS package Helske (2017).

The Kalman filter outputs estimates of the state vector (δ_t, γ_t)^′. The estimates at time t conditional on information up to and including time t are denoted (δ̂_t ∣ t, γ̂_t ∣ t)^′ and given by the contemporaneous filter; the predictive filter estimates the state at time t + 1 from the same information set, outputting (δ̂_{t + 1 ∣ t}, γ̂_{t + 1 ∣ t})^′.

It may be useful to review past movements of the state vector (δ_t, γ_t)^′. This can be done using the smoothed estimates (δ̂_t, γ̂_t)^′, which denotes the estimates of the state vector at time t based on all T observations in the series.

Estimation of the unknown variance parameters (σ_ε² and σ_ζ²) is by maximum likelihood (ML) and is carried out using KFAS following the procedure described in Helske (2017). We retain the option of either estimating the signal-to-noise ratio q_ζ, or of fixing it at a plausible value. In practice, for coronavirus applications, we set the value of q_ζ based on experience and judgment, reducing the number of parameters to be estimated by one. Tests for normality and residual serial correlation are based on the standardized innovations, that is one-step ahead prediction errors, v_t = ln g_t − δ_{t ∣ t − 1}, t = 3, ..., T.

Daily effects, which are generally quite pronounced in the coronavirus data, can be included in the model as described in the Appendix.

Forecasts and peak prediction

Forecasts of future observations are obtained from the predictive recursions

so that and Ŷ_{T + ℓ ∣ T} = μ̂_{T + ℓ ∣ T}; the initial value is μ̂_T ∣ T = Y_T.

We construct forecast intervals for y_t based on the prediction intervals for δ_t. The conditional distribution of future values of δ̂_t is Gaussian. We replace δ̂_T + j|T in (5) with the upper bound of a prediction interval for δ̂_T + j|T to compute the upper bound of our forecast interval for y_T + ℓ and likewise for the lower bound.² In effect, the forecast intervals are based on inference on the log cumulative growth rate, δ_t.

The filtered growth rate ĝ_{y, t ∣ t} of new cases y_t, can be extracted from the continuous-time incidence curve: μ^′(t) = g(t)μ(t), where μ(t) is the growth curve and g(t) is its growth rate. Taking logarithms and differentiating we get

where ĝ_t ∣ t = exp δ̂_t ∣ t. The sampling variability of ĝ_t ∣ t is dominated by that of γ̂_t ∣ t (see Harvey and Kattuman (2021)). Therefore when constructing confidence intervals for ĝ_t ∣ t we treat ĝ_y, t as if it has a normal distribution centered on ĝ_{y, t ∣ t} with variance Var(γ̂_t ∣ t).

Even when the nowcast ĝ_{y, T ∣ T} is positive and daily cases are growing, there will be a saturation level for the cumulative total, Y_t, so long as γ̂_T ∣ T is negative. The nowcasts of y_t peak when ĝ_{y, t ∣ t} = 0, which requires γ̂_t ∣ t to be sufficiently negative to outweigh ĝ_t ∣ t, which is, of course, always positive. This can be seen from the expression for the growth rate of daily cases:

When γ̂_T ∣ T is negative, there is a flattening of the curve and a signaling of an upcoming peak in the trend of y_t. As shown in [HK, p10], the peak in the trend is predicted to be ℓ_T days ahead where³

$$\ell_{T}=\frac{\ln (-\hat\gamma_{T\mid T})-\hat\delta_{T\mid T}}{\hat\gamma_{T\mid T}}=\frac{\ln (-\hat\gamma_{T\mid T}/\hat{g}_{T\mid T})}{\hat\gamma_{T\mid T}},\;\;\;\; -\hat{g}_{T\mid T}<\hat\gamma_{T\mid T}<0. $$

The generation of forecasts is demonstrated in Section 4.

Reproduction Number R_t

The path of the epidemic is best tracked by nowcasts and forecasts of g_y, t, the growth rate of y_t, which are constructed by HK from the filtered estimates in the state space model, (2), (3) and (4). Wallinga and Lipsitch (2007) describe how the estimates of g_y, t can be translated into estimates of the instantaneous reproduction number R_t. Harvey and Kattuman (2021) propose where τ is the generation interval – the typical number of days between an infected person becoming infected and them transmitting the disease to someone else. We construct credible intervals for R̃_t, τ and R̃_τ, t^e by substituting the upper and lower bounds of the confidence intervals for g_y, t into (8) to get the upper and lower bounds of the credible intervals. See Harvey, Kattuman, and Thamotheram (2021) for an application.

The estimates of R_t can be used for tracking and forecasting the epidemic. The nowcasts of y_t peak when ĝ_{y, t ∣ t} = 0, corresponding to R̃_t, τ = R̃_τ, t^e = 1. Based on (7), predictions of g_y, t are given by

We can then obtain predictions of R_t, as in (8). If γ̂_T ∣ T is zero, the estimated growth of y_t is exponential and it is helpful to characterize it by the doubling time, ln 2/ĝ_{y, T ∣ T} = 0.693exp (−δ̂_T ∣ T).

When exp δ̂_T ∣ T + γ̂_T ∣ T > 0, the nowcast ĝ_{y, T ∣ T} is positive and the estimate of R_t given by (8) is greater than one. So long as γ̂_T ∣ T is negative, then as T → ∞, R̃_{τ, T + ℓ ∣ T}^e → exp (τγ̂_T ∣ T) < 1, and a saturation level for Y appears on the horizon.

We now turn to case where γ_t potentially turns positive in a typically short-lived phase, as a new wave emerges.

Reinitalizing the data series

Let t = r denote the re-initialization date and let r₀ denote the date at which the cumulative series is set to 0. Then:

where Y_t^r denotes the cumulative cases after re-initialization. We set Y_r − 1^r = 0, so that the growth rate of cumulative cases is available from t = r + 1 onwards. Note that Y_t^r = Y_t − Y_r0.

The gap between the two series becomes apparent by writing

In the next section, where we illustrate the working of the program, it can be seen that in contrast to the original ln g_t series, which continues to increase, the reinitialized ln g_t series begins to decrease from the reinitialization date. The reinitialization enforces the canonical Gompertz curve with the log of growth rate of cumulative cases sloping down.

Reinitalizing the model

We reinitialize the model by specifying the appropriate prior distribution for the initial states: α₁^r ∼ N(a₁^r, P₁^r) with a₁^r = (a_δ, 1^r, a_γ, 1^r)′ and P₁^r defined as follows. Let ℱ_t = ln g_t, ln g_t − 1, …, ln g₁ and define a_t = E(α_t|ℱ_t − 1), a_t|t = E(α_t|ℱ_t), P_t = Var(α_t|ℱ_t − 1), P_t|t = Var(α_t|ℱ_t). Then,

where a_{δ, r + 1} and P_r + 1 are obtained from the non-reinitialized model estimated over t = 1, …, r via the usual Kalman filter recursions. Adding ln (Y_r/y_r) to a_{δ, r + 1} corrects for the shift down in the log cumulative cases caused by reinitializing the cumulative case series. Setting a_γ, 1^r = 0 ensures the model starts off with exponential, rather than super-exponential, growth.

We reinitialize the model through the priors in this way rather than simply re-estimating the model from scratch for two reasons. First, it allows us to impose a proper (rather than diffuse) prior centered on zero for γ, so that the starting point is exponential growth. Second, it enables us to make use of data from before the reinitialization date. One needs a reasonable sample size for the estimated model and forecasts to be reliable, but if a new wave is taking off, forecasts need to be generated quickly. This was particularly true with the emergence of the Omicron variant, which caused an explosive increase in infection over a short period of time.

We do not re-estimate the σ_ε² or σ_ζ² parameter in the reinitialized model. Rather, we use the values estimated in the original model over t = 1, …, r. The one-step-ahead prediction error at t = r is the same in both the initialized and reinitialized models, but after t = r, the prediction errors diverge.

The reinitialization procedure is very similar in the case where we have seasonal terms. If we let α_s, t be the vector of seasonal states and maintain an analogous notation to that above, the prior mean of the seasonal components in the reinitialized model is a_s, 1^r = a_{s, r + 1}. The prior variance of α₁^r remains P_r + 1 where P_r + 1 is appropriately re-defined to include the seasonal term, as described in the Appendix.

Setting up the forecasting exercise

We begin by specifying a number of options for the forecasting exercise, as defined below.

• Y is the data, in the form of time series of cumulative confirmed cases. In this example the object holding this series is called gauteng.

• estimation.date.start is the date of the first observation in the sample to be used for estimating the model. By default, it is the first date in the xts object Y.

• estimation.date.end is the date of the last observation in the sample to be used for estimating the model. By default, it is the last date in the xts object Y.

• n.forecasts is the number of days or periods for which forecasts are to be made. E.g., if n.forecasts = 14, forecasts will be generated for up to 14 days following estimation.date.end.

• q is the signal-to-noise ratio, which controls the smoothness of the estimated trend. A lower value will lead to more smoothness. By default, we use q = 0.005, which in our experience ensures a good balance between the smoothness of the trend and the speed with which changes in estimates respond to new observations. Alongside, q can be estimated and compared with the default value.

• confidence.level sets the coverage of the confidence intervals for ln (g_t) which is then used to generate the prediction intervals for forecasts. Here, we use 0.68, corresponding to the probability that the forecast lies within one standard deviation of the point forecast.

• plt.length sets a truncation date to enhance the clarity of plots, e.g. showing only the last 30 days of the estimation sample. The date range for plotting can be set as plt.length days upto estimation.date.end.

In this example the data is the cumulative confirmed cases time series for Gauteng. The start and end dates (estimation.date.start and estimation.date.end) that define the sample used for estimation are chosen as appropriate for the exercise. In this example, we begin with the sample period set from 1 February to 19 April 2021. This marks the beginning of the third wave in Gauteng as can be seen in Figure 1. The options are specified as below.

file.path = here()
res.dir <- here::here(file.path, 'results')
date.format <- "%Y-%m-%d"
Y = gauteng
estimation.date.start = as.Date("2021-02-01")
estimation.date.end = as.Date("2021-04-19")
n.forecasts = 14
q = 0.005
confidence.level = 0.68
plt.length = 30

Estimation

We begin by selecting the data series (Y) for the defined sample period.

idx.est =
    (zoo::index(Y) >= estimation.date.start) &
    (zoo::index(Y) <= estimation.date.end)
y = Y[idx.est]

We then estimate the model using a diffuse prior distribution for the initial state vector. The signal-to-noise ratio can be left as a free parameter to be estimated, as in the code below.

model_q <- SSModelDynamicGompertz$new(Y = y)
res_q <- model_q$estimate()

In the rest of this example we estimate the model setting the signal-to-noise ratio at 0.005. As mentioned, in our experience this value strikes a useful balance between the smoothness of the estimate of the slope parameter γ, and the speed with which it adapts to new observations.

model = SSModelDynamicGompertz$new(Y = y, q = q)
res = model$estimate()

Results

We can now plot the forecast of ln g_t - the log of the growth rate of Y, the cumulative cases – which is the transformation of the data series that is taken to the model, and we can compare these forecasts to the actual ln g_t series. We do this by passing the output (res) of the estimation step along with an evaluation sample to a plotting function. We specify the evaluation sample by converting the cumulative cases series to the log of the growth rate of the cumulative cases.

First, we create the evaluation sample.

y.eval = Y %>%
subset(index(.) > tail(res$index,1)) %>%
tsgc::df2ldl()

tsgc::plot_forecast then creates and plots forecasts of ln (g_t).

tsgc::plot_forecast(
 res=res,
 y.eval = y.eval, n.ahead = n.forecasts,
 plt.start.date =  tail(res$index, 1) - plt.length
)

Figure 2: Fourteen-day forecast of ln(g_t) from 20 April 2021, for Gauteng province in South Africa.

From these results we can recover the forecasts of new cases from 20 April 2021, with their prediction intervals.

tsgc::plot_new_cases(
  res, Y=y,
  n.ahead=n.forecasts,
  confidence.level=confidence.level,
  date_format="%Y-%m-%d",
  plt.start.date = tail(res$index, 1) - plt.length
)

Figure 3: Fourteen-day forecast of new cases from 20 April 2021 for Gauteng province in South Africa

To assess accuracy, we plot these forecasts against the actual new cases, that have been held back from the estimation sample, using the plot_holdout function. The model forecasts are compared with the the first differences of Y.eval, the cumulative series for the forecast window.

tsgc::plot_holdout(res, Y=y,
                   Y.eval = Y[(tail(res$index,1)+0:n.forecasts)],
                   confidence.level = 0.68,
                   date_format = "%Y-%m-%d")

Figure 4: Accuracy of the fourteen-day forecast of new cases from 20 April 2021 for Gauteng province in South Africa.

Figure 4 shows that the forecasts were accurate over the first seven days, with a mean absolute percentage error (MAPE) of 13.9%. Note that reported cases were unusually low on 27 April due to the fact that it is Freedom Day and a public holiday in South Africa. That day aside, over the six days from 28 April the MAPE was 12.8%. Over the full 14 days of the forecasts, the MAPE was 27%.

As discussed earlier, the reproduction numbers R_t and their 68% credible intervals can be calculated. The plot in Figure 5 reveals that R_t remains above one through the period, indicating that the (third) wave in Gauteng had launched by this time.

r.t <- tail(exp(res$get_gy_ci()*gen_int),7)

## # A tibble: 7 × 5
##   Date          Rt lower upper name   
##   <date>     <dbl> <dbl> <dbl> <chr>  
## 1 2021-04-13  1.18 1.05   1.34 Gauteng
## 2 2021-04-14  1.14 1.01   1.29 Gauteng
## 3 2021-04-15  1.13 1.01   1.28 Gauteng
## 4 2021-04-16  1.12 0.996  1.27 Gauteng
## 5 2021-04-17  1.06 0.941  1.20 Gauteng
## 6 2021-04-18  1.02 0.906  1.15 Gauteng
## 7 2021-04-19  1.06 0.944  1.20 Gauteng

Figure 5: Reproduction numbers for the 7-day period to 19 April 2021, for Gauteng province in South Africa.

A CSV file (named ) is written to the directory specified. The forecast options specified earlier are retained.

tsgc::write_results(
 res=res,
 res.dir = res.dir,
 Y=Y,
 n.ahead = n.forecasts,
 confidence.level =  confidence.level
)

Reinitialization

In all countries the coronavirus pandemic was characterized by a series of recurring waves due to a combination of biological, behavioral, and environmental reasons. In an epidemic, the onset of a new wave is signalled when the slope parameter γ, which measures of the growth rate of the growth rate of new cases, rises above zero for a sustained period. Such a super-exponential phase of the epidemic in which the growth rate of new cases is itself increasing over time is typically short.

This section illustrates the reinitialization procedure which allows us to apply the model to the new wave as it begins, without jettisoning information from the wave that has just ended. We extend the estimation window to 25 June 2021, by which date the third wave is well on course with its peak within sight (see Figure 1). All other options remain the same.

estimation.date.end = as.Date("2021-06-25")

Triggering reinitialization

It is not obvious, a priori, when precisely to reinitialize the model. Based on experiments a reasonable option is to reinitialize when the estimate of the slope parameter, γ_t, breaches a threshold of two standard errors above zero, and at that point backdate reinitialization to when the estimate of γ_t first turned positive. In applying the above heuristic there is a choice between the filtered slope estimate and the smoothed slope estimate. Experiments suggests that the greater noisiness of the filtered estimate of γ_t often triggers reinitialization too early. The smoothed estimate is more reliable.

Figure 6 shows that for the third wave in Gauteng, the smoothed slope estimate exceeded twice its standard error on 1 May 2021, having risen above zero on 21 April 2021.

Figure 6: Trigger for reinitialization: when filtered slope exceeds twice its standard error above zero, reinitialization is triggered to the date when the filtered slope crossed zero.

The date for reinitialization is set accordingly.

reinit.dates = "2021-04-29"

Estimating the reinitialized model

SSModelDynGompertzReinit takes the same arguments as SSModelDynGompertz, with the addition of the reinit.dates argument.

model <- SSModelDynGompertzReinit$new(
  Y = y, q = q,
  reinit.date = as.Date(reinit.dates,
  format = date.format)
)
res.reinit <- model$estimate()

We generate the reinitialized data frame by setting cumulative cases to 0 at the appropriate point, as discussed in the Theory section and extract the evaluation sample from the reinitialized data frame as below.

y.eval.reinit <- Y %>%
  reinitialise_dataframe(., reinit.dates) %>%
  df2ldl() %>%
  subset(index(.) > tail(res.reinit$index,1))

Estimating the model with the reinitialized series, the actual and forecast ln (g_t) can be plotted as in Figure 7, which is analogous to Figure 2 for the non-reinitialized series.

tsgc::plot_forecast(
    res=res.reinit,
    y.eval = y.eval.reinit,
    n.ahead = n.forecasts,
    plt.start.date =  tail(res.reinit$index, 1) - plt.length,
    title='Forecast of $\ln(g_t)$  after reinitialization.'
)

Figure 7: Forecast of ln(g_t) after reinitialization.

The plot of the forecasts of new cases, and that of these forecasts against the actual number of new cases can be produced as before. See Figures 8 and 9. The trend has begun to turn down in model forecasts with with the reinitialized series.

tsgc::plot_new_cases(
  res.reinit, Y=Y.reinit, n.ahead=n.forecasts,
  confidence.level=confidence.level,
  date_format="%Y-%m-%d",
  plt.start.date = tail(res.reinit$index, 1) - plt.length
)

Figure 8: Forecast of new cases after reinitialization.

tsgc::plot_holdout(res = res.reinit, Y=Y.reinit[index(y)],
                   Y.eval = Y[(tail(res.reinit$index,1)+0:n.forecasts)],
                   confidence.level = 0.68,
                   date_format = "%Y-%m-%d")

Figure 9: Forecast accuracy of the reinitialized model over the hold-out sample period: 14 days past 25 June 2021.

Comparing forecasts, without reinitialization the 14-day forecast MAPE is 41.9% (see Figure 10). With reinitialization, it falls to 20.2%. The corresponding MAPE figures for 7-day forecasts is 15.2% without reinitialization and 9.5% with reinitialization.

Figure 10: Forecast accuracy of the model without reinitialization over the hold-out sample period: 14 days from 25 June 2021.

Just as for the standard (non-reinitialized) model, the returned estimation results contained in res.reinit are a FilterResults object, and can be written to CSV using the write_results method.