Data
We obtained weekly updates of the daily confirmed monkeypox cases by reporting date from publicly available sources by the CDC and the Our World in Data (OWID) GitHub repository [24, 25]. At the global level and for countries that have reported the great majority of the cases, including Brazil, Canada, England, France, Germany, Spain, and the USA, we retrieved daily case series from the GitHub Our World in Data (OWID) repository [8, 25]. We reported forecasts based on CDC and OWID team data for the USA. The CDC and OWID data sources define a confirmed case as a person with a laboratory-confirmed case of monkeypox [26, 27]. Data were downloaded every Wednesday evening from the CDC and every Friday afternoon from the GitHub Our World in Data (OWID) repository from the week of July 28th, 2022, through the week of October 13th, 2022. For the week of July 28th, 2022, data posted by the OWID team on August 9th, 2022, was used to produce the forecast as it was the earliest data available.
The n-sub-epidemic modeling framework
A detailed description of our modeling framework is given in ref. [19]. In this n-sub-epidemic modeling framework, epidemic trajectories are modeled as the aggregation of overlapping and asynchronous sub-epidemics. A sub-epidemic follows the 3-parameter generalized-logistic growth model (GLM), which has displayed competitive performance [28,29,30]. This model is given by the following differential equation:
$$\frac{dC(t)}{dt}={C}^{\prime }(t)=r{C}^p(t)\left(1-\frac{C(t)}{K_0}\right),$$
where C(t) denotes the cumulative number of cases at time t and \(\frac{dC(t)}{dt}\) describes the curve of daily cases over time t. The parameter r is positive, denoting the growth rate per unit of time, K0 is the final outbreak size, and p ∈ [0, 1] is the “scaling of growth” parameter which allows the model to capture early sub-exponential and exponential growth patterns. If =0 , this equation describes a constant number of new cases over time, while p = 1 indicates that the early growth phase is exponential. Intermediate values of p (0 <p < 1) describe early sub-exponential (e.g., polynomial) growth dynamics.
An n-sub-epidemic trajectory comprises n overlapping sub-epidemics and is given by the following system of coupled differential equations:
$$\frac{d{C}_i(t)}{dt}={C_i}^{\prime }(t)={A}_i(t){r}_i{C_i}^{p_i}(t)\left(1-\frac{C_i(t)}{{K_0}_i}\right),$$
where Ci(t) tracks the cumulative number of cases for sub-epidemic i, and the parameters that characterize the shape of the ith sub-epidemic are given by (ri, pi, K0i), for i = 1, …, n. Thus, the 1-sub-epidemic model is equivalent to the generalized growth model described above. When n > 1, we model the onset timing of the (i + 1)th sub-epidemic, where (i + 1) ≤ n, by employing an indicator variable given by Ai(t) so that the (i + 1)th sub-epidemic is triggered when the cumulative curve of the ith sub-epidemic exceeds Cthr.
The (i + 1)th sub-epidemic is only triggered when Cthr ≤ K0i. Then, we have:
$${A}_i(t)=\left\{\begin{array}{c}1,{C}_{i-1}(t)>{C}_{thr}\\ {}\ \\ {}0, Otherwise\end{array}\right.i=2,\dots n,$$
where A1(t) = 1 for the first sub-epidemic. Hence, the total number of parameters that are needed to model an n-sub-epidemic trajectory is given by 3n + 1. The initial number of cases is given by C1(0) = I0, where I0is the initial number of cases in the observed data. The cumulative curve of the n-sub-epidemic trajectory is given by:
$$C_{tot}(t)=\sum_{i=1}^nC_i(t).$$
Hence, this modeling framework is suitable for diverse epidemic patterns including those characterized by multiple peaks.
Parameter estimation for the n-sub-epidemic model
The time series of new weekly monkeypox cases are denoted by:
\({y}_{t_j=}{y}_{t_1,}{y}_{t_2},\dots, {y}_{t_{n_d}}\)where j = 1, 2, …, nd
Here, tj are the time points for the time series data, nd is the number of observations. Using these case series, we estimate a total of 3n + 1 model parameters, namely Θ = (Cthr, r1, p1, K01, …, rn, pn, K0n). Let f(t, Θ) denote the expected curve of new monkeypox cases of the epidemic’s trajectory. We can estimate model parameters by fitting the model solution to the observed data via nonlinear least squares [31] or via maximum likelihood estimation assuming a specific error structure [32]. For nonlinear least squares, this is achieved by searching for the set of parameters \(\hat{\varTheta}\)that minimizes the sum of squared differences between the observed data \({y}_{t_j=}{y}_{t_1,}{y}_{t_2}\dots ..{y}_{t_{n_d}}\)and the model mean f(t, Θ). That is, Θ = (Cthr, r1, p1, K01, …, rn, pn, K0n) is estimated by \(\hat{\varTheta}=\mathit{\arg}\mathit{\min }\ {\sum}_{j=1}^{n_d}{\left(f\left({t}_j,\varTheta \right)-{y}_{t_j}\right)}^2\).
We quantify parameter uncertainty using a bootstrapping approach described in [33], which allows the computation of standard errors and related statistics in the absence of closed-form solutions. To that end, we use the best-fit model \(f\left(t,\hat{\varTheta}\right)\) to generate B-times replicated simulated datasets of size nd, where the observation at time tjis sampled from a normal distribution with mean \(f\left({t}_j,\hat{\varTheta}\right)\) and variance \(\frac{\sum_{j=1}^{n_d}{\left(f\left({t}_j,\hat{\varTheta}\right)-{y}_{t_j}\right)}^2}{n_d-\left(3n+1\right)}\). Then, we refit the model to each B simulated dataset to re-estimate each parameter. The new parameter estimates for each realization are denoted by \({\hat{\varTheta}}_b\)where b = 1, 2, …, B. Using the sets of re-estimated parameters \(\left({\hat{\varTheta}}_b\right),\) the empirical distribution of each estimate can be characterized, and the resulting uncertainty around the model fit can similarly be obtained from \(f\left(t,{\hat{\varTheta}}_1\right),\)\(f\left(t,{\hat{\varTheta}}_2\right),\dots, f\left(t,{\hat{\varTheta}}_B\right)\). We run the calibrated model forward in time to generate short-term forecasts with quantified uncertainty.
Selecting the top-ranked sub-epidemic models
We used the AICc values of the set of best fit models based on one and two subepidemics to select the top-ranked sub-epidemic models. We ranked the models from best to worst according to their AICc values, which is given by [34, 35]:
$${AIC}_c={n}_d\mathit{\log}(SSE)+2m+\frac{2m\left(m+1\right)}{n_d-m-1}$$
where \(SSE={\sum}_{j=1}^{n_d}{\left(f\left({t}_j,\hat{\varTheta}\right)-{y}_{t_j}\right)}^2\), m = 3n + 1 is the number of model parameters, and nd is the number of data points. Parameters from the above formula for AICc are estimated from the nonlinear least-squares fit, which implicitly assumes normal distribution for error.
Constructing ensemble n-sub-epidemic models
We generate ensemble models from the weighted combination of the highest-ranking sub-epidemic models as deemed by the \({AIC}_{c_i}\) for the ith ranked model where \({AIC}_{c_1}\le \dots \le {AIC}_{c_I}\) and i = 1, …, I. An ensemble derived from the top-ranking I models is denoted by Ensemble(I). Thus, Ensemble (2) refers to the ensemble model generated from the weighted combination of the top-ranking 2 sub-epidemic models. We compute the weight wi for the ith model, i = 1, … , I, where ∑wi = 1 as follows:
\(w_i=\frac{l_i}{l_1+l_2+\dots+l_I}\;for\;all\;i=1,2,\dots,I,\)
where li is the relative likelihood of model i, which is given by \({l}_i={e}^{\left(\left({AIC}_{min}-{AIC}_i\right)/2\right)}\) [36], and hence wI ≤ … ≤ w1 . The prediction intervals based on the ensemble model can be obtained using a bootstrap approach similar as before. We employed the first-ranked and the second-ranked models to derive the ensemble forecasts. AICc values of the top models for the most recent forecast can be found in figure 1s (Additional file 1) [24, 25].
Forecasting strategy
Using a 10-week calibration period for each model, we have conducted 324 real-time weekly sequential 4-week ahead forecasts across studied areas and models (week of July 28th–week of October 13th, 2022) thus far. At the national and global levels, we also report forecasting performance metrics for 8 sequential forecasting periods covering the weeks of July 28th, 2022, through September 15, 2022, for which data was available to assess the 4-week ahead forecasts. We also compare the predicted cumulative cases for the 4-week forecasts across models for a given setting. Cumulative cases for a given model were calculated as the sum of median number of new cases predicted during the 4-week forecast. Forecasts were evaluated using data reported during the week of October 13th, 2022.
Performance metrics
Across geographic areas, we assessed the quality of our model fit and performance of the short-term forecasts for each model by using four standard performance metrics: the mean squared error (MSE) [37], the mean absolute error (MAE) [38], the coverage of the 95% prediction interval (PI) [37], and the weighted interval score (WIS) [20, 39]. While MSE and MAE assess the average deviations of the mean model fit to the observed data, the coverage of the 95% PI and the weighted interval score consider the uncertainty of the forecasts.
