{'fraction_non_alphanumeric': 0.03826586234464714, 'fraction_numerical': 0.022603386874046078, 'mean_word_length': 4.832768122085629, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Forecasting Based on Surveillance Data

Forecasting the future course of epidemics has always been one of the main goals of epidemic modelling. This chapter reviews statistical methods to quantify the accuracy of epidemic forecasts. We distinguish point and probabilistic forecasts and describe different methods to evaluate and compare the predictive performance across models. Two case studies demonstrate how to apply the different techniques to uni- and multivariate forecasts. We focus on forecasting count time series from routine public health surveillance: weekly counts of influenza-like illness in Switzerland, and age-stratified counts of norovirus gastroenteritis in Berlin, Germany. Data and code for all analyses are available in a supplementary R package.


Introduction
Epidemic modelling has at least three distinct aims: Understanding the spread of infectious diseases, identifying suitable measures to control the spread of an epidemic, for example through isolation or vaccination (Daley and Gani, 1999, Section 1.5), and predicting the future course of an epidemic. Mathematical models are often used to better understand the dynamics of infectious disease spread and the effects of control measures (Keeling and Rohani, 2008, Section 1.5), but are less oriented towards predictions. In recent years, more emphasis has been placed on the development of predictive models and methods. The goal of this chapter is to review the literature in this area and to describe how general statistical principles from the forecasting literature can be applied to evaluate the quality of epidemic forecasts. The described methods will also be illustrated in two case studies, where we assess competing approaches to forecast time series of infectious disease counts.

The history of forecasting epidemics
Predicting the future course of epidemics has been a desire of mankind for a long time. Scientific forecasting based on mathematical models dates back to the pioneering work by Baroyan et al. (1971), who predicted the course of an influenza epidemic for the main cities in the Soviet Union. The rise of new infectious diseases has been a major trigger of novel forecasting methods, such as for AIDS (Daley and Gani, 1999, Section 6.2) and SARS (Hufnagel et al., 2004). Furthermore, meteorologic forecasting methods have been adopted in epidemiological research, including modelling (Viboud et al., 2003;Shaman and Karspeck, 2012) and assessment (Paul and Held, 2011;Held et al., 2017) techniques. Recent developments in influenza forecasting have focussed on the integration of search logs from Google (Dukic et al., 2012;Shaman and Karspeck, 2012;Yang et al., 2015), social media data from Twitter (Paul et al., 2014), combinations thereof , Wikipedia article views (Generous et al., 2014;Hickmann et al., 2015), or human mobility data (Pei et al., 2018). Compared to weather forecasting, epidemic forecasting is still in its infancy, and the human component makes it particularly challenging (World Health Organization, 2014;Moran et al., 2016).

Forecasting competitions
In order to develop better epidemic forecasting methods, the World Health Organization has joined forces in an informal consultation with more than 130 global experts (World Health Organization, 2016), and the Centers for Disease Control and Prevention in the USA have organized several seasonal forecasting competitions (Biggerstaff et al., 2016(Biggerstaff et al., , 2018. Realtime influenza forecasts are now provided by Nicholas Reich and co-workers at the ReichLab (Tushar et al., 2017), see http://reichlab.io/flusight/. Other recent competitions include the DAPRA challenge on forecasting chikungunya (Del Valle et al., 2018), the White House/NOAA challenge on forecasting dengue (Buczak et al., 2018), and the RAPIDD ebola forecasting challenge (Viboud et al., 2018).

Forecasting targets
Several quantities are of interest in epidemic forecasting, such as timing of and incidence in the peak week , onset week (Pei et al., 2018), cumulative incidence (Lega and Brown, 2016), weekly incidence (Paul and Held, 2011;Reich et al., 2016a), outbreak size and duration (Farrington et al., 2003), and the epidemic curve (Jiang et al., 2009). Also of public health relevance are stratified forecasts of the above quantities, for example by region or by age group . Forecasting targets for seasonal influenza epidemics in particular have been reviewed previously (Chretien et al., 2014;Nsoesie et al., 2014). More recently, the weekly proportion of doctor visits due to influenza-like illness is becoming a popular forecasting target (Biggerstaff et al., 2018).

Evaluating point forecasts
A point forecast is usually made for a continuous or integer (often count) outcome, for example the disease incidence in the next week, the number of newly confirmed cases in a certain time interval, or the timing of the peak week of an epidemic. Several measures have been proposed and used to evaluate the quality of a point forecast. Gneiting (2011) gives a comprehensive discussion of suitable measures for point predictions.
To introduce some notation, letŷ i denote the point forecasts of the observations y i , i = 1, . . . , n. Given a suitable non-negative scoring function S(ŷ i , y i ), the overall predictive performance can be assessed with the mean scorē Scoring functions are usually negatively oriented, so the smaller a score, the better the forecast. Commonly used scoring functions are the absolute error (AE) S(ŷ, y) = |ŷ − y|, the squared error (SE) S(ŷ, y) = (ŷ − y) 2 , the absolute percentage error S(ŷ, y) = |ŷ − y|/y and the relative error S(ŷ, y) = |ŷ − y|/ŷ. Note that the latter two scoring functions require y andŷ to be positive, respectively. Other summary measures used in the infectious disease literature may be based on one of the scoring functions mentioned above, for example the root mean squared error (RMSE) or the relative mean absolute error (Hyndman and Koehler, 2006;Reich et al., 2016b). The latter is defined as the ratio of the mean absolute errors of two competing forecasts and is not to be confused with the mean relative absolute error of one forecasting method. The small simulation study reported by Gneiting (2011) reveals the key differences of the four scoring functions listed above. It is based on a time series model commonly used in econometrics, but the results equally apply to other fields. A simple conditionally heteroscedastic Gaussian time series model is used to generate a non-negative time series y 1 , . . . , y n . Model-based one-step-ahead point forecasts (based on the mean of the forecast distribution) are then compared to naive approaches, such as the so-called "fence-sitter forecast", a forecast with constant predictions. Quite surprisingly, one of the naive methods outperforms the model-based one-step-ahead forecasts both under the absolute error and the absolute percentage error scoring functions. Gneiting (2011) points out that the model-based forecast is only optimal under squared error loss but not necessarily for other loss functions, where other functions of the predictive distribution will be optimal. A scoring function is called consistent for a loss function, if the expected score is minimized when following this loss function. The squared error scoring function is consistent for the mean, and the absolute error scoring function for the median, both standard loss functions. However, the mean absolute percentage error commonly used in influenza forecasting (Chretien et al., 2014, Table 4), is consistent for a nonstandard and rather exotic functional. To quote Gneiting (2011), "it thus seems prudent that authors using this functional consider the intended or unintended consequences and reassess its suitability as a scoring function." Some scoring functions may be problematic per se, such as the commonly used correlation coefficient between predictions and observations. For example, point predictions always twice as large as the observations are obviously inappropriate, but the correlation coefficient between predictions and observations will be one.

Evaluating probabilistic forecasts
In many areas of science, researchers argue that forecasts should be probabilistic (Gneiting and Katzfuss, 2014) and this plea has recently been taken up in the infectious disease literature Funk et al., 2018).
Evaluating probabilistic forecasts requires suitable scores to quantify the "distance" between a cumulative distribution function F , the forecast, and the observation y that later realizes. Proper scoring rules play a key role, as they allow for a fair comparison of different probabilistic forecasting methods. In addition, the notions of calibration and sharpness are important. Calibration is a property of both the forecast and the observations, whereas sharpness is a property of the forecasts only. The paradigm of probabilistic forecasting is to "maximize sharpness subject to calibration" .
Much of the literature on probabilistic forecasting focusses on continuous forecast distributions. However, in infectious disease epidemiology, the quantities to be predicted are often (incidence or prevalence) counts, which require suitable extensions of the assessment techniques (Czado et al., 2009;Held et al., 2017).

Calibration and sharpness
Calibration is defined as the statistical consistency of probabilistic forecasts and the observations . For continuous outcomes, calibration is usually assessed with the probability integral transform (PIT). Specifically, the PIT is F (y), so will be uniformly distributed if the observation y is a realisation from the forecast F . In practice we will compute PIT values for a series of forecasts and visually assess their distribution in a histogram. Modifications for count data exist (Czado et al., 2009).
For binary data, the calibration slope studies the association between the (logittransformed) predicted probabilities with the binary outcomes in a logistic regression (Cox, 1958;Steyerberg, 2009). Calibration can also be visually assessed in a calibration curve (Gneiting and Katzfuss, 2014).
Statistical calibration tests are commonly employed to assess the evidence for miscalibration. Different methods exist for continuous outcomes (Mason et al., 2007;Held et al., 2010), count outcomes (Wei and Held, 2014) and binary outcomes (Spiegelhalter, 1986), but they may require certain distributional assumptions on the forecasts.
For continuous forecasts, sharpness is usually defined as the width of the associated prediction intervals. For multivariate forecasts, sharpness is defined based on the predictive covariance matrix, for example its determinant (Gneiting et al., 2008).

Proper scoring rules
Proper scoring rules are summary measures of the predictive performance of probabilistic predictions Y i ∼ F i , i = 1, . . . , n, allowing for a joint assessment of calibration and sharpness (Gneiting and Katzfuss, 2014). As in Section 6.2 for point forecasts, the mean score across forecasts is usually reported and compared across different forecasting methods. For ease of presentation we drop the index i in this section and compare the probabilistic forecast Y ∼ F with the actual observation y.
The definition of propriety is mathematically somewhat demanding and we refer the interested reader to the review paper by . Put simply, a proper scoring rule ensures that a forecaster reports his actual forecast and does not obtain a better score in expectation by digressing from his belief. Here we take scores to be negatively oriented penalties that forecasters wish to minimize, so the smaller a score, the better the forecast.
For forecasts of binary data y ∈ {0, 1}, commonly used scores are the logarithmic score (LS), i. e. the negative log-likelihood of the observation under the forecast distribution, and the Brier score (BS), also known as probability score. Specifically, let π denote the predicted probability of the observation y = 1, then Both scores are proper, whereas the absolute score AS(F, y) = |π − y| can be shown to be improper (Held and Sabanés Bové, 2014, Chapter 9), so should not be used.
A probabilistic forecast for count data y ∈ {0, 1, 2, . . .} is represented by probabilities π k = Pr(y = k), k ∈ {0, 1, 2, . . .}, where ∞ k=1 π k = 1. The logarithmic score for the observation y = j now simply reads LS = log(π j ) where π j is the probability of the observation y = j. An extension of the Brier score to forecasts of count data is the ranked probability score where P j = j k=0 π k and 1 denotes the indicator function. Probabilistic count forecasts are often summarized with the mean µ and the variance σ 2 . A proper score that only uses these first two moments is the Dawid-Sebastiani score This score can easily be generalized to multivariate predictions, which depends only on the mean vector µ and the covariance matrix Σ of the predictive distribution. The first term log|Σ| in (6.1) is called the log determinant sharpness (logDS) and is recommended as a (multivariate) measure of sharpness (Gneiting et al., 2008). To avoid unnecessarily large numbers it is common practice to report scaled versions of logDS and mDSS, obtained through division by 2d where d is the dimension of the observation y.
A possible alternative to mDSS is the energy score, a multivariate extension of the ranked probability score (Gneiting et al., 2008;Held et al., 2017). The question arises which scoring rule should be used in practice. The logarithmic score is known to be sensitive to outliers if π j = 0 is close to zero. The ranked probability score is reported to be less sensitive to extreme observations  and provides an attractive alternative. The DSS is particularly useful if the first two predictive moments are available but not the whole predictive distribution. It has also been argued that the choice of a scoring rule should take into account the costs of bad forecasts (Merkle and Steyvers, 2013). In practice we recommend to evaluate several scores to obtain a more robust comparison of predictive performance.
The computation of the different scores depends on whether the forecast distributions are known analytically or derived from simulations. Care has to be taken with simulationbased forecasts as scores can become numerically instable, for example the logarithmic score for count forecasts will be infinite if an observation y = j has empirical frequency equal to zero in the forecast samples. A possible remedy is to apply a Rao-Blackwell/importance sampling scheme (Gelfand and Smith, 1990) and to average the conditional predictive distributions, see also Ray et al. (2017, Supp Mat Section 3.2). Alternatively one may use kernel density estimation as implemented in the R package scoringRules (Jordan et al., 2018). The Dawid-Sebastiani score for multivariate forecasts can be numerically unstable if the predictive covariance matrix Σ is estimated from the empirical covariance matrix of Monte Carlo samples (Scheuerer and Hamill, 2015).
For independent forecasts, a simple paired t-test or a permutation test can then be used to assess the evidence for differences in predictive performance between two competing forecasting methods. If the forecasts are dependent as in a series of sequential k-step ahead forecasts, the Diebold-Mariano test can be used to account for the correlation between scores (Diebold and Mariano, 1995;Gneiting and Katzfuss, 2014).

Applications
We now describe two applications, in which we compare the quality of forecasts provided by different models and methods. The data and code to reproduce these analyses are available in the R package HIDDA.forecasting (Meyer, 2018) at https://HIDDA.github.io/ forecasting. The code is presented in several package vignettes, which also give some additional results.

Univariate forecasting of influenza incidence
We consider a time series of weekly incidence counts on influenza-like illness (ILI) in Switzerland from 2000 to 2016 (Figure 6.1). The last 213 weeks (starting from 2012-12-04) are used to assess one-week-ahead forecasts. Four sets of 30-weeks-ahead long-term forecasts (shaded time periods in Figure 6.1) are computed for each of the last four seasons, conditional upon all data prior to the respective first week of December.
We have used five different methods to produce forecasts of influenza activity in Switzerland: an ARMA(2,2) time series model as estimated by auto.arima() from the forecast package (Hyndman and Khandakar, 2008), an observation-driven ARMA model for negative-binomial counts as implemented in the glarma package (Dunsmuir and Scott, 2015), the endemic-epidemic negative-binomial time-series model hhh4 from the surveillance package , Facebook's forecasting tool prophet (Taylor and Letham, 2018), as well as a recently proposed method based on kernel conditional density estimation , which we implemented following the code provided at https://github.com/reichlab/article-disease-pred-with-kcde. Note that forecasts from arima and prophet are generated on a log-scale. All models were configured to account for yearly seasonality. In the first three models, seasonality was represented by parametric sine-cosine regressors (Held and Paul, 2012) with frequency 2π/52.1775 (derived from the average number of calendar weeks per year). For prophet and kcde, we followed the documented examples (see the supplementary vignettes for the exact model formulations). We also included a separate Christmas effect for calendar week 52 (not with kcde). For reference, we computed naive historical forecasts for the 213 weeks of the test period, based on fitting log-normal distributions to the observed counts in the same calendar week of previous years.
We computed the RMSE and the mean DSS and LS for the one-week-ahead forecasts and the long-term forecasts. The results are summarized in Table 6.1. As expected, the one-weekahead foreacasts have smaller (i. e. better) mean scores than the long-term forecasts for all methods considered. Note that we have not attempted to compute long-term forecasts with kcde due to excessively long runtimes. The time taken to compute a single one-week-ahead forecast varied between 0.001 (naive) and 3 (prophet) seconds.
The best one-week-ahead forecasts in terms of all scores are obtained with hhh4 followed by glarma, whose point forecasts are slightly worse with an average error of 2450 cases. The arima and kcde methods come next and the worst one-week-ahead predictions are obtained from prophet, which achieves similar scores as the naive approach. Figure 6.2 shows the weekly forecasts and associated scores as well as the overall PIT histograms, which are computed based on the method for count data (Czado et al., 2009). There is no clear evidence for miscalibration of any of the one-week-ahead forecasts, but the first three and the naive (not shown) forecasting methods have a distinct peak of the histogram in the first bin, indicating some evidence for biased or underdispersed predictions. The prophet and kcde (not shown) forecasts don't have this problem. PIT histograms for the kcde and naive forecasts can be found in the corresponding supplementary vignettes.
An interesting result is obtained from the long-term forecasts, where the best scores are achieved by taking simple historical averages, except for RMSE where glarma is better. For the more sophisticated forecasts shown in Figure 6.3, the ranking is less clear. While hhh4 is the best model in terms of DSS, arima is the best model in terms of LS. Quite surprising is the large DSS value of the glarma method. Closer inspection of Figure 6.3 suggests that this is caused by the excessively large uncertainty of the glarma long-term predictions, compared to the other three methods. From Figure 6.4 we can see that the prophet method has a particularly large DSS score in the 2015/2016 season, where the DSS score is roughly twice as large as for the other methods. Closer inspection of Figure 6.3 reveals that this is due to under-prediction of the incidence combined with an insufficient amount of uncertainty at the end of the epidemic season.
It may also be of interest to predict the peak week in each of the four years based on the long-term forecasts. For illustration, we have computed a probabilistic prediction of the peak week based on the samples from the hhh4 model only. The median peak week forecasts is always calendar week 5. This is quite close to the observed peak weeks 6, 8, 6, and 6, respectively, shown in Figure 6.1. However, considerable uncertainty is attached to these predictions with 2.5% quantile the second week and 97.5% quantile week 16 of each season.

Multivariate forecasting of norovirus incidence
To illustrate statistical methods for the assessment of multivariate forecasts we use agestratified surveillance data on norovirus gastroenteritis from Berlin, Germany. We use the well-established hhh4 modelling framework from the surveillance package to generate forecasts, building on our previous analyses of these data Held et al., 2017). Here we only consider models for spatially aggregated counts Y gt , g = 1, . . . , G, i. e. without additional stratification by city district. The surveillance data are available from the package hhh4contacts (Meyer, 2017), cover five norovirus seasons from 2011-W27 to 2016-W26, and are stratified into G = 6 age groups: 0-4, 5-14, 15-24, 25-44, 45-64, and 65+ years of age. These groups reflect distinct social mixing of pre-school vs. school children, and intergenerational mixing. Figure 6.5 shows the age-specific incidence time series. Over the five years, the reported numbers of cases by age group were 2783, 326, 553, 1909, 2530, and 8335, respectively, corresponding  FIGURE 6.5: Age-stratified incidence of norovirus gastroenteritis in Berlin, Germany. The calendar weeks 52 and 1 (Christmas break period) are highlighted. The last 52 weeks (to the right of the vertical dashed line) are used to assess forecasts. Note that the y-axis employs a square-root scale.
Given the counts from the previous week, Y g ,t−1 , g = 1, . . . , G, we assume Y gt to follow a negative binomial distribution with a group-specific overdispersion parameter and mean The endemic log-linear predictor ν gt contains group-specific intercepts, a Christmas break effect (via a simple indicator for the calendar weeks 52 and 1), and group-specific seasonal effects via sin(ωt) and cos(ωt) terms, ω = 2π/52. The epidemic log-linear predictor φ gt also contains group-specific intercepts, but shared seasonality and no Christmas break effect. The contact weights c g g are estimated from the German subset of the POLYMOD survey (Mossong et al., 2008), taking into account the reciprocal nature of contacts (Wallinga  et al., 2006). The resulting contact matrix (Figure 6.6) is subsequently row-normalized to a transition matrix, removing differences in group-specific overall contact rates. The model described above is the same as reference model 6 in Held et al. (2017, Section 3.3). As alternative models we consider homogeneous mixing between age groups (c g g = 1), no mixing between age groups (c g g = 1{g = g}) and a model where we estimate a power transformation C κ of the contact matrix as described in . See the vignette("BNV") in the supplementary package HIDDA.forecasting for how to implement these models in R.
We have fitted the models to the first four seasons (2011-W27 to 2015-W26) and evaluated the quality of forecasts in the subsequent season (2015-W27 to 2016-W26). Figure 6.7 shows the fitted mean components from the AIC-optimal model with power-adjusted contact matrix. The proportion of the incidence that can be explained by the epidemic component varies between the age groups. The endemic-epidemic decomposition actually resembles Figure 4 in Meyer and Held (2017) quite well, which is based on a different model, additionally stratified by city district. Within the epidemic part, the incorporated social contact matrix results in predominant within-group reproduction of the disease for the 0-4 and 65+ age groups, whereas the 25-44 and 45-64 year-old persons inherit a substantial proportion of cases from other age groups.
Turning to the predictive performance, Figure 6.8 shows the one-week-ahead forecasts and associated scores during the last season for the power-adjusted contact model. Average scores of the different models are compared in Table 6.2. Although the power-adjusted model gives the best AIC score on the training data, it is outperformed by the "no mixing" model on the test season, both in terms of DSS and LS. This can be explained by Table 6.3, which gives the log-likelihood contributions of the different models from a separate fit to the whole time period (training and test). We can see that the "no mixing" model has indeed a larger log-likelihood contribution in 2015/16 than all other three models. However, the poweradjusted model is best in 2011/12, 2013/14 and 2014/15 and second best in 2012/13 (after "reciprocal"). This explains why the power-adjusted model performs best in the training, but not in the test period.   Multivariate long-term predictions from the four models are evaluated in Table 6.4 for the test period. The first two columns give mDSS and logDS for multivariate predictions (of dimension 6×52) across weeks and age groups. The "no mixing" model is again best in terms of mDSS, although both the "reciprocal" and the power-adjusted model provide sharper predictions in terms of logDS. Predictions of the age-stratified total number of cases (now of dimension 6) are displayed in Figure 6.9 and the corresponding scores listed in the last two columns of Table 6.4. More uncertainty of the predictions (in age group 65+) is visible for the "no mixing" and power-adjusted models, resulting in larger values of logDS. However, the increased uncertainty seems to be benificial for the quality of the forecasts, yealding the smallest mDSS scores for these two formulations. In contrast, the "homogeneous" prediction barely covers the observed count in age group 65+, resulting in a poor mDSS score. For comparison, a naive forecaster who independently fits negative binomial distributions to the age-stratified sizes of the past four seasons would reach mDSS = 4.084 (with logDS = 3.820), which is superior to the predictions from the homogeneous contact model.

Discussion
We have provided a non-comprehensive review of the literature on forecasting epidemics. We have focussed on statistical methods to quantify the accuracy of predictions, distinguishing between point and probabilistic forecasts. Two applications show how the different techniques can be applied to uni-and multivariate forecasts. Inspired by similar techniques in weather forecasting and other areas of science, recent work in this area has focussed on model averaging and stacking in order to improve the predictive quality of single model forecasts (Ray and Reich, 2018). Perhaps the biggest challenge to epidemic forecasting is the incorporation of reporting delays and underreporting, as described in Chapter V.3 of this handbook. A rigorous analysis requires surveillance and internet search data archived in a way where you can actually see what was available at a given time (McIver and Brownstein, 2014). Such real-time epidemiological data should become the standard in surveillance systems, facilitating the development of novel forecasting techniques.