key: cord-0459971-kmrp7y11
authors: Thomas, Maud; Rootz'en, Holger
title: Real-time prediction of severe influenza epidemics using Extreme Value Statistics
date: 2019-10-23
journal: nan
DOI: nan
sha: ebbe88012f95c2f364c6464e2d55237d84b2a102
doc_id: 459971
cord_uid: kmrp7y11

Each year, seasonal influenza epidemics cause hundreds of thousands of deaths worldwide and put high loads on health care systems. A main concern for resource planning is the risk of exceptionally severe epidemics. Taking advantage of recent results on multivariate Generalized Pareto models in Extreme Value Statistics we develop methods for real-time prediction of the risk that an ongoing influenza epidemic will be exceptionally severe and for real-time detection of anomalous epidemics and use them for prediction and detection of anomalies for influenza epidemics in France. Quality of predictions is assessed on observed and simulated data.

Every year, seasonal influenza epidemics cause 250,000-500,000 deaths worldwide [Rambaut et al., 2008] and put high strain on public health care systems in a short time frame, due essentially to increased doctor visits, and overcrowded emergency departments and intensive care units. Predicting the likelihood of an exceptionally severe epidemic in the future is of paramount importance for health resource planning [Bresee and Hayden, 2013, Khan and Lurie, 2014] . The three following questions are thus of central interest to public health policy makers.

(Q1) estimation of risks of occurrence of a very severe epidemic during the following years, (Q2) real-time prediction of the risk that an ongoing epidemic will be exceptionally severe, and (Q3) real-time detection of unusual, thus potentially dangerous, epidemics.

In France, countrywide data on seasonal influenza morbidity have been available since 1985 (see Figure 1 below). The Sentinelles network monitors cases of Influenza-like Illness (ILI) defined by the presence of fever in excess of 39 • C, respiratory symptoms, and muscle pains [Réseau Sentinelles, 2019] . Though only a fraction of the ILI reported visits are in fact caused by influenza, their total number reflects the burden on the health care system. Figure 1 suggests that influenza epidemics, or at least a proportion of them, may be conceptualised as extreme episodes that occur within the ILI time series. This led us to lean on Extreme Value Statistics (EVS) to address the questions above. EVS have been developed to handle extreme events, such as extreme floods, heat waves or episodes with huge financial losses [e.g. Katz et al., 2002 , Embrechts et al., 1997 . EVS allow for prediction of risks of episodes outside of the observed range. The contributions of this paper are to develop methods for real-time prediction of severe epidemics through estimation of conditional probabilities of exceedances of very high levels; methods for detection of anomalous epidemics; and to use the new methods for prediction of influenza epidemics in France.

Question (Q1) may be answered by standard and well established EVS methods such as the univariate Peaks over Threshold method described in Section 2.1 (for a more thorough presentation see Coles [2001] ). For instance, they were applied by Chen et al. [2015] to avian influenza, and by Thomas et al. [2016] to influenza mortality, and emergency department visits.

Previous approaches to Question (Q2) include high dimensional times series prediction [Davis et al., 2016] and the method of analogues [Réseau Sentinelles, 2019] . The US Centre for Disease Control initiated a data challenge to predict the 2013-2014 US influenza epidemic. Their conclusion was "Forecasting has become technically feasible, but further efforts are needed to improve forecast accuracy so that policy makers can reliably use these predictions" [Biggerstaff et al., 2016] . Building on recently developed EVS results and methods based on multivariate Generalized Pareto (GP) distributions [Rootzén et al., 2018 , Kiriliouk et al., 2019 , we extend the toolbox of available techniques and develop methods to tackle Question (Q2). Question (Q3) is an anomaly detection problem. Recent publications have used EVS in this context [Guillou et al., 2014 , Thomas et al., 2017 , Goix, 2016 , Chiapino et al., 2019 . In the present paper, the multivariate GP models estimated in the course of solving Question (Q2) are used to detect anomalous epidemics.

For predictions to be useful in practice, it is necessary to have an understanding of their reliability. Methods for evaluating the quality of prediction of extremes have only very recently started to be approached in the literature. Here we use a strategy that uses standardized Brier scores, Precision-Recall curves and Average Precision scores [Steyerberg et al., 2010 , Brownlee, 2020 , Saito and Rehmsmeier, 2015 .

Section 2 describes the methods used in this study and in particular the multivariate GP distributions. Section 3 presents the Sentinelles network data and the definition of epidemics. In Section 4 the EVS methods are applied to the Sentinelles ILI data. Finally, Section 5 contains the conclusion.

The main interest of EVS methods is to allow prediction outside of the range of the observations. In this section, we describe and develop the EVS tools required to answer the questions raised in this paper.

Section 2.1 presents the univariate Peaks over Thresholds (PoT) method which is used for Question(Q1). Our approach to Question (Q2) is described in Section 2.2: it is based on the multivariate PoT method, building on recent results on multivariate EVS models [Rootzén et al., 2018 , Kiriliouk et al., 2019 . Section 2.3 deals with Question (Q3) and combines a general framework for anomaly detection with Generalized Pareto modelling. Finally, a strategy for assessing real-time prediction of exceedances of very high levels is proposed in Section 2.4.

The one-dimensional PoT method, introduced in the hydrological literature and combined with the univariate GP distribution in Smith [1984] consists in selecting a suitably high threshold u and defining excesses above u as the differences between the observations and u. Under general conditions, the conditional distribution of the excesses, given that they are positive, is asymptotically as u → ∞ a GP distribution with cumulative distribution function (cdf)

Here σ > 0 is a scale parameter, and γ ∈ R is a shape parameter. For a ∈ R, (a) + = a if a ≥ 0 and 0 if a < 0. The parametrisation is chosen so that the cdf for γ = 0 is the limit as γ → 0 of the cdf with non-zero γ. A useful account of this method is given in [Coles, 2001] .

Assuming the observations are independent and identically distributed with common cdf F , the parameters of the cdf (1) are estimated from the observed excesses. For y > u, F (y) is estimated, by

where H is the GP cdf (1) with parameters replaced by their estimated values, and p u is the empirical frequency of observations which exceed the threshold u [Coles, 2001] .

The level that the maximum of n observations will exceed with probability 1 − α is estimated by the α-th quantile y α of F (y) n . For example, if γ = 0, H is an exponential distribution then

for α such that p u ≥ 1 − α 1/n .

The multivariate PoT method was introduced by [Rootzén and Tajvidi, 2006 , Michel, 2009 , Brodin and Rootzén, 2009 . A suitably high threshold is chosen for each component of a d-dimensional random vector Y = (Y 1 , . . . , Y d ). Let denote u = (u 1 , . . . , u d ) the d-dimensional vector of thresholds. Excess vectors X = (X 1 , . . . , X d ) are defined as the component-wise differences between observations and thresholds and considered as positive if at least one of the components exceeds its threshold. Under general conditions, the joint distribution of the positive excess vectors is asymptotically a multivariate GP distribution as the thresholds tend to ∞. For the sake of completeness, we briefly recall the definitions and properties of multivariate GP distributions that will be needed in this paper [for further details, see Rootzén et al., 2018 , Kiriliouk et al., 2019 . For simplicity, we describe the method for d = 3, but it can easily extended to any dimension d.

Unlike in the univariate case, the family of multivariate GP distributions cannot be described as a parametric family. Rootzén et al. [2018] have developed four representations of multivariate GP distributions for which closed formulas of densities are available. In this paper, we use the U -representation since it presents nicer properties across the dimension, which are useful for prediction. Moreover, at the price of standardization, the marginals of the distributions can be assumed to be standard exponential distributions [see Section 3 in Kiriliouk et al., 2019] .

Let U = (U 1 , U 2 , U 3 ) be a 3-dimensional random vector such that E[e max U ] < ∞, where max U = max{U 1 , U 2 , U 3 }, and let f U be the probability density function of U . For i = 1, 2, 3, let f i and F i denote the density and distribution functions of U i , respectively. The vector U or its distribution is referred to as the generator.

According to Equation (3.4) in [Kiriliouk et al., 2019] , the density function h U of the 3-dimensional GP distribution with standard exponential margins generated by U is

where x = (x 1 , x 2 , x 3 ) and x + log t = (x 1 + log t, x 2 + log t, x 3 + log t). The indicator function 1 {x 0} equals one if at least one of the components of x are positive, and is zero otherwise. Different choices of distributions for U yield different GP models.

Our answer to Question (Q2) is to estimate the conditional probability that Y 3 will exceed some level > u 3 given Y 1 and Y 2 , or equivalently, for X = (X 1 , X 2 , X 3 ) = (Y 1 − u 1 , Y 2 − u 2 , Y 3 − u 3 ), the conditional probability that X 3 > − u 3 given X 1 and X 2 .

For this we assume that X is positive and follows a multivariate GP distribution with generator U,

(4)

Formulas for general U distributions are given in the Supplementary Material (Proposition 1, Section 3).

Question (Q3) belongs to the field of anomaly detection. In the present context, a natural approach is to use the estimated GP model from Section 2.2 to detect new observations that exhibit a significantly different pattern from the data used to fit the model. A statistical test for anomalous observations may be based on the GP negative log-likelihood, with a very large value suggesting that the new observation might be anomalous. The quantiles of the GP negative log-likelihood distribution were estimated by simulation in order to define the decision region of the test [e.g. Section 2 of Root et al., 2015] . However, it must be stressed that "anomalous" has to be understood with respect to the fitted GP model.

There is a substantial literature on assessment of forecasting [see Lerch et al., 2017, and references therein] . This literature provides metrics aimed at comparing predictions with observations. However, standard assessment metrics are not appropriate when the frequency of exceedances in the data is small and sometimes zero. This issue was tackled in [Renard et al., 2013 ], but they aimed at a rather different context, spatial environmental data, where prediction is done at not just one, but many stations. Here, we wish to compare estimated prediction probabilities and the observed outcome for just one data set, not many.

For this purpose, we use (i) standardized Brier scores; (ii) Precision-Recall Curves; and (iii) Average Precision scores, together with simulations from estimated models. We also illustrate prediction quality visually with stratified plots of predicted probabilities of exceedance.

(i) The standardized Brier score is defined as

where N is the number of predictions, p i is the prediction probability of exceedance, o i = 1 if an exceedance was observed, and 0 otherwise, and p = 1 [Steyerberg et al., 2010] . The score is bounded by 1, with larger values indicating better prediction. Using the predictors p i = p gives the value 0.

(ii) A question such as "will the maximum of the past observations be exceeded in the future?" may be formulated as a binary classification problem. The data are divided into two classes: Positives (exceedances) and Negatives (no exceedances). The strategy consists first in computing p i and then in comparing this estimate to some cut-off probability value p c . If p i ≥ p c then the observation is assigned to the Positives class, and to the Negatives class otherwise. The "true positives" ("false positives") correspond to the observations that are correctly (incorrectly) assigned to the Positives class, and the "true negatives" ("false negatives") to the observations that are correctly (incorrectly) assigned to the Negatives class.

Common methods to assess the performance of binary classifiers include true positive and true negative rates, and ROC (Receiver Operating Characteristics) curves. These methods, however, are uninformative when the classes are severely imbalanced, which is the case when predicting very high level exceedances which are rare by nature. In this context, Saito and Rehmsmeier [2015] have argued that Precision-Recall curves are more informative. These curves as the cut-off probability p c varies from 0 to 1. Precision quantifies the number of correct positive predictions out all positive predictions made; and Recall (often also called Sensitivity) quantifies the number of correct positive predictions out of all positive predictions that could have been made. Both focus on the Positives class (the minority class) and are unconcerned with the Negatives (the majority class). The Precision-Recall curve of a skillful model bows towards the point with coordinates (1, 1). The curve of a no-skill classifier will be a horizontal line on the plot with a y-coordinate proportional to the number of Positives in the dataset. For a balanced dataset this will be 0.5 [Brownlee, 2020] .

(iii) The Average Precision score is an approximation to the area under the Precision-Recall curve [Su et al., 2015] . A perfect prediction model would have an Average Precision score equal to 1, and the closer the score is to 1, the better the prediction performance of the model.

The French nationwide Sentinelles network consists of approximately 1,500 general practitioners in France who participate on a voluntary basis in the ILI surveillance, and report new cases of ILI observed in their practice. Based on these data, nationwide weekly ILI incidence rates-i.e. numbers of new cases in France per week per 100,000 individualsare estimated from a Poisson model which also provides confidence intervals [for further details, see Réseau Sentinelles, 2019, Souty et al., 2014] . Estimating incidence is difficult, especially since there is no available gold standard for comparison. The quality of the data produced by the Sentinelles network has been assessed in comparison with ILI incidence estimated from independent sources [see e.g. Kalimeri et al., 2019 , Carling et al., 2013 , Pelat et al., 2017 . In the Sentinelles network, epidemics are identified using the Serfling method [Serfling, 1963] . It consists in fitting a cyclic regression model to the weekly ILI rates and setting the start of the epidemic at the first of the first two consecutive weeks during which the ILI incidence rates exceed the upper bound of the 90% prediction interval [for further details see Réseau Sentinelles, 2019].

Weekly ILI incidence rates from January 1985 to February 2019 were downloaded for analysis. Figure 1 shows the time series of weekly ILI incidence rates, which includes 35 epidemics. The epidemic with the highest peak corresponds to the 1989-epidemic, with a value of 1,793. The lowest peak was observed in 2014, with a value of 325. The durations of the Serfling epidemics range from 5 to 16 weeks, in 1991/2014 and 2010, respectively. In this study, the 1985-2018 data were used for estimation and the 2019 data were kept as a test sample. Since we were interested in very high ILI incidence rates and relied on EVS methods, we focused on the most active part of the epidemics. The Serfling method was thus adapted as follows. The start of the epidemic was set at the first of the first two consecutive weeks during which the ILI incidence rates exceeded the 0.88-quantile (=272) of the 1985-2018 data. The end was set at the end of the Serfling epidemic. This quantile was chosen to synchronise the peaks of the epidemics as much as possible, see Figure 2: a) Weekly ILI incidence rates for epidemics identified with the Serfling method, and b) Weekly ILI incidence rates for epidemics obtained from the definition in this paper.

3 to 12 weeks in 2014 and 1985/2018, respectively. The size of an epidemic was defined as the sum of the weekly ILI incidence rates from the start to the end of the epidemic. The smallest epidemic size was 847 in 2014 and the largest was 8,062 in 1989.

The size of an influenza epidemic depends on multiple factors including immunity and vaccination prevalence in the population. Some of these factors may likely be impacted by the sizes of past epidemics. Nevertheless, there is little indication that this translates into statistical dependence between epidemics-as shown by the correlation and distance correlation plots in the Supplementary Material (Section 1)-and we thus assumed that the ILI incidence rates of different epidemics were mutually independent.

4 Results: prediction of very high ILI loads on the French health care system

In this section, the methods described in Section 2 are applied to the Sentinelles ILI data.

Recall that the 1985-2018 data were used for estimation and the 2019 data were kept as a test sample. For each epidemic, we shall refer to the first week of the epidemic as Week 1, the second week as Week 2 and so on until the end of the epidemic. For j = 1, 2, 3 and k = 1985, . . . , 2019, we let Y k j denote the ILI incidence rate of Week j in year k, and S k denote the size of the epidemic of year k. We omit the index k when we refer to a generic epidemic.

To address Questions (Q1), (Q2) and (Q3), we choose to focus on predictions for the Week 3 ILI incidence rate Y 3 and the epidemic size S given observation of the Week 1 and 2 incidence rates. This choice was driven by both epidemiological and technical reasons. Policy makers will be predominantly interested in early predictions of the ongoing epidemic. Thus, basing predictions on observation of the first two weeks Table 1 : Estimated levels that the Week 3 ILI incidence rate and the epidemic size will exceed with either probability 10% or 1% during the following year and the following 10 years.

one year one year 10 years 10 years 1 − α 10% 1% 10% 1% Week 3 ILI incidence rates 1,192 2,094 2,076 2,994 Epidemic sizes 6,165 9,452 9,385 12,733 was a reasonable compromise between accuracy of prediction and early prediction. If a very high ILI incidence rate in Week 3 was predicted, this would suggested a rapid intensification of the ongoing epidemic, and prediction of very large epidemic size would mean that the ongoing epidemic will be severe. Technically the number of parameters in models of dimension above 3 was too large for reliable computation, given the amount of available observations.

The univariate PoT method was applied to the Week 3 ILI incidence rates (Y 1985 3 , . . . , Y 2018

3 ) and to the epidemic sizes (S 1985 , . . . , S 2018 ).

As explained in Section 2.1, the first step is to choose a suitably high threshold for each series of observations. The threshold must be high enough to ensure that the asymptotic model is valid, but low enough to yield a sufficient number of positive excesses. For the ILI rates, the threshold u I was chosen as the 0.9-quantile (=339) of the whole 1985-2018 series; for epidemic sizes, the threshold u S as the 0.6-quantile (=4,144) of the series of epidemic sizes from 1985 to 2018. These thresholds yield 30 positive excesses for Week 3 ILI rates and 14 for epidemic sizes.

One-dimensional GP distributions (Equation (1)) were fitted to the ILI rate and size positive excesses. Likelihood ratio tests showed that the hypothesis γ = 0 was not rejected (p = 0.64, and p = 0.98 respectively), so that cdf's were assumed exponential. QQ-plots and estimates of scale parameters are shown in the Supplementary Material (Section 2, Figure 2 and Table 1) . Table 1 presents estimates of the levels y α such that during the following year and during the following 10 years the Week 3 ILI incidence rate and the epidemic size will exceed them with given probability 1 − α. These estimates were computed using Equation (2). For Week 3 ILI incidence rates, u I = 339 yielded p u I = 0.88 and for epidemic sizes u S = 4, 441 yielded p u S = 0.41.

For the 2019 epidemic, the Week 3 ILI incidence rate was 599 and the epidemic size was 1,505. The observations were well below the 90% estimated quantile for the Week 3 ILI incidence rate and the epidemic size given in Table 1 . For example, it was estimated that there is a 10% probability that the epidemic size will exceed 9,385 at least once during the next 10 years. 

. . , 2018. The thresholds u I of 339 for the ILI incidence rates Y 1 , Y 2 and Y 3 and u S of 4, 144 for the epidemic size S were as in the previous section. For both (Y 1985 To define a multivariate GP model, the distribution of the generator U = (U 1 , U 2 , U 3 ) must be specified. Assuming that the three components U 1 , U 2 , U 3 of U were mutually independent, three families of marginal distributions were considered: Gumbel, reverse exponential, and reverse Gumbel distributions. The corresponding models were fitted to the 32 positive excess vectors. Formulas for the densities h U for the three GP models are given in the Supplementary Material (Section 4). Table 2 shows that both in terms of AIC and BIC, the best fit was that of the GP family with Gumbel generator, for both Y and S. According to Equation (3), the corresponding density is

with α 1 , α 2 , α 3 > 1 and β 1 , β 2 , β 3 ∈ R (for further details see Supplementary Material, Section 4 and [Kiriliouk et al., 2019] ). In the sequel, we refer to this GP model as the Gumbel model. The estimated parameters of the Gumbel model are given in the Supplementary Material (Table 2 ). More parsimonious sub-models were considered and consistently rejected on the basis of AIC, BIC and log-likelihood ratio tests (Supplementary Material, Table 3 ). The estimated model was used to provide estimates of the probability that Y 3 and S will exceed a specified level given that Y 1 and Y 2 have been observed. Since Equation (4) is valid for positive excess vectors only, two situations must be considered i) If at least one of Y 1 or Y 2 exceeds its threshold, Equation (4) is valid whatever Y 3 (or S).

ii) If neither Y 1 nor Y 2 exceeds its threshold, then whether the excess vector will be positive or not is unknown. In this case, the right-hand side of Equation (4) must be multiplied by the probability that Y 3 (or S) exceeds its threshold given that Y 1 and Y 2 do not exceed theirs. The corresponding empirical probability was 0.33 for both Y 3 and S.

The procedure is illustrated in Table 3 which presents predictions for the 2019 epidemic. The largest observed Week 3 ILI incidence rate between 1985 and 2018 was 1,729. The table shows the estimated probabilities that the 2019 Week 3 ILI incidence rate exceeds a fraction κ of the 1985-2018 maximum ILI incidence rate 1,729 for κ = 0.5, 0.75, 0.95, 1. Estimates for epidemic sizes are presented similarly with a largest observed epidemic size of 8,062. The prediction probabilities, even for the lowest level were quite small, and in fact this level was not exceeded in 2019. For the 2019 epidemic the Weeks 1, 2 and 3 ILI incidence rates were 366, 540, and 599, respectively, and the epidemic size was 1,505.

The quantiles of the GP negative log-likelihood were obtained as follows: the estimated Gumbel model was used to generate 1,500 datasets, each consisting of 33 threedimensional positive excess vectors. For each simulated dataset, a Gumbel model was fitted to the 32 first vectors and the estimated negative log-likelihood was computed at the 33rd vector. The significance levels obtained as the empirical quantiles of these negative log-likelihoods are shown in Table 4 .

To illustrate the meaning of "anomalous", we simulated two extra positive excess vectors: (i) a positive excess vector with a very high third component equal to the 0.99-quantile of the third component of simulated positive excess vectors and (ii) an anomalous positive excess vector obtained from (i) by multiplying the first and the third components by 1.5 and the second by 0.5. Figure 3 shows that only the anomalous point (ii) exceeds the 0.999-quantile. Interestingly, the 2009-10 swine flu pandemic was not unusual compared to the other epidemics.

Quality of real-time predictions was assessed on both the Sentinelles 1985-2019 series and on simulated data. For purpose of comparison, prediction probabilities were also estimated using a standard logistic regression model, with Y 1 and Y 2 as covariates and response variable coded as 1 in the presence of an exceedance and 0 otherwise. Levels used in this section are those defined in Table 3 (Section 4.2).

Assessment on the 1985-2019 ILI data A leave-one-out procedure was performed on the 1985-2019 epidemics yielding 35 estimates of prediction probabilities for the two lower levels (κ = 0.5, 0.75). Prediction probabilities for the two higher levels could not be estimated since there was only one exceedance for κ = 0.95 and none for κ = 1. Figure 4 shows the prediction probabilities of level exceedances stratified according to whether an exceedance was observed or not. Contrary to GP prediction, logistic regression was never able to discriminate between the two outcomes. Table 5 presents the corresponding standardized Brier scores and confirms that GP prediction performs much better than the logistic regression.

Precision-Recall curves and Average Precision scores are not shown since they were uninformative due to insufficient data. Table 5 : Standardized Brier scores derived from a leave-one-out procedure on observed data for the GP prediction and the logistic regression for predictions of exceedances of 1,729×κ for Week 3 ILI incidence rates and 8,062×κ for epidemic sizes.

Week 3 ILI incidence rates Epidemic sizes κ 0.5 0.75 0.5 0. Assessment on simulated data 1,500 datasets consisting of 33 three-dimensional vectors were simulated from the estimated Gumbel models for Week 3 ILI incidence rates and epidemic sizes, respectively. The simulations were carried out following Section 7 of [Rootzén et al., 2018] . A Gumbel model was fitted to the first 32 vectors of each dataset and the estimated model was then used to predict the third component of the 33rd vector, conditionally on the first two components. Figure 5 shows boxplots of the prediction probabilities stratified according to whether an exceedance was observed or not, for the GP prediction for both Week 3 ILI incidence Logistic regression q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Outcome Prediction probability GP prediction q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q rates and epidemic sizes for the four levels (κ = 0.5, 0.75, 0.95, 1). A table is associated with each boxplot. These tables first give the proportion of 0-outcomes and of 1-outcomes, which shows that the classes are indeed very imbalanced, and then the quartiles of the prediction probabilities. One can see that the boxes for the 0-outcomes are very concentrated and close to 0 while the boxes for the 1-outcomes are larger. There is no overlap between the two boxes. Moreover, the widths of the boxes indicate that the performance of the GP prediction is better for incidence rates than for epidemic sizes. The boxplots for predictions obtained from the logistic regression are shown in Figure 6 for the two lower levels. Predictions could not be made for the higher levels since, for these levels, there were too few exceedances in the simulated data to allow estimation of the parameters. The quality of the GP prediction is much better than that of the logistic regression.

Standardized Brier scores and Average Precision scores for predictions with GP and logistic predictions are presented in Table 6 , and Precision-Recall curves are shown in Figure 7 . For the simulated data, the true parameters for the Gumbel model are known. The prediction probabilities obtained from the true model, perhaps surprisingly, give results which are not noticeably better than the results from fitted Gumbel model (Figure 3 

Extreme events that can put high stress on the health care system are a major concern for resource planning in public health. Both real-time prediction of the development of an ongoing epidemic and prediction of sizes of future epidemics are important for policy makers.

In Section 4.1 we used existing univariate EVS to estimate risks of very high ILI incidence rates and of extreme epidemic sizes in France during the following year, and during the following 10 year, as input for long term resource planning.

A main result of this paper is the development of a methodology for real-time prediction of risks as early as at the start of an epidemic. The prediction method is based on estimates of conditional probabilities of exceedance of a very high level, using recent results on multivariate GP distributions [Rootzén et al., 2018 , Kiriliouk et al., 2019 . Sections 4.2 and 4.4 indicate that our predictions work well for the French Sentinelles data. We chose to focus on prediction of very high levels for the third week of the epidemic and of epidemic size, given observation of the first two weeks of the epidemic. Policy makers may be most interested predicting the overall severity of epidemic as measured by epidemic size. However, according to Section 4.4, our method gives better predictions for the Week 3 incidence rates than for epidemic sizes. This is as expected since the epidemic size is the sum of the incidence rates from the beginning to the end of the epidemic and predicting it hence requires prediction further into the future than just to the following week .

The choice of 3-dimensional GP models was motivated both by epidemiological and by technical considerations. The methods can be used to predict any other characteristic of the epidemic, e.g. the fourth week given the first two weeks, and in other dimension if the number of available data allows it, and provided dependence is strong enough to make prediction is possible. We also developed a method to use observation of the first three weeks of an epidemic to warn decision makers if something anomalous and potentially dangerous was happening.

Real-time prediction provides an important support for health care resource planning. However, our methodology is only useful for diseases, like influenza, for which historical data are available but does not apply to emerging diseases such as Covid-19.

To a large extent, assessment of predictions of extreme events remains an open issue since by definition there are always few observations of very extreme events, which does not fit well into standard theory. Here we used an assessment strategy which uses standardized Brier scores, Precision-Recall curves and Average Precision scores [Steyerberg et al., 2010 , Brownlee, 2020 , Saito and Rehmsmeier, 2015 .

Web-based supporting materials for "Real-time prediction of severe influenza epidemics using Extreme Value Statistics"

by Maud Thomas and Holger Rootzén

We made autocorrelation, cross-correlation, auto-distance-correlation, and cross-distance-correlation plots of all pairs of the ILI incidence rates of Weeks 1, 2 and 3, and epidemic sizes. None of the plots suggested dependence. For example, see Figure 1 . Table 1 gives the estimates of the scale parameters for the exponential fit of the Week 1, 2 and 3 ILI incidence rates and epidemic sizes and the corresponding QQ-plots are shown in Figure 2 . The p-values for the likelihood ratio tests for the hypothesis γ = 0 were p = 0.66, p = 0.57 and p = 0.64, p = 0.98 for the Week 1, 2 and 3 ILI incidence rates and epidemic sizes, respectively. Observed Theorical c) d) Figure 2 : Quantile-quantile plots of the exponential fit of the excesses of the 0.9-quantile of the ILI incidence rates for a) ILI incidence rates for Week 1, b) ILI incidence rates for Week 2, c) ILI incidence rates for Week 3 and of the excesses of the 0.6-quantile of the sizes for d) the epidemic sizes.

Generalizing the notation in Section 3.2, write x i:j = (x i , . . . , x j ) and x = (x 1 , . . . , x d ). Then, using the same conventions as in Section 3.2, the U-representation of the density of d-dimensional GP distributions is

where σ = (σ 1 , . . . , σ d ) and γ = (γ 1 , . . . , γ d ). Note that in the general case h U depends both on the generator U and on the marginal parameters σ, γ. For σ = 1, γ = 0 this expression reduces to Equation (3) in Section 3.2. The next proposition gives the conditional prediction formulas for prediction of the exceedances of the dth component given the first c components of a d-dimensional vector.

Proposition 1. Let X = (x 1 , . . . , x d ) a d-dimensional GP vector with density f U and let 1 ≤ c < d. The probability that x d will exceed a threshold v d given x 1:c = x 1:c is given by

ii) for x 1:c ≤ 0 and v d > 0

iii) for x 1:c ≤ 0 and v d ≤ 0

Example. For d = 3 and c = 2 we get

1 if x 1:2 ≤ 0 and v 3 ≤ 0.

4 Probability density functions for the 3-dimensional GP distributions

In this section, we recall the density functions for the 3-dimensional GP distributions that we use in Section 4.2 (for further details see Section 4.1. in [Kiriliouk et al. 2019] ). Let U = (U 1 , U 2 , U 3 ) be a 3-dimensional random vector such that E[e max U ] < ∞, let f U be the probability density function of U and f i the density function of its i-th component U i , for i = 1, 2, 3.

If U = (U 1 , U 2 , U 3 ) with U i , i = 1, 2, 3 independent and distributed according to a Gumbel distribution with parameters α i and β i . The marginal density f i is given by

where the requirement E[e max U ] < ∞ imposes the restriction α i > 1 for i = 1, 2, 3. The density of a 3-dimensional GP distribution with independent Gumbel generators is then given by

Generators with independent reverse exponential components If U = (U 1 , U 2 , U 3 ) with U i , i = 1, 2, 3 independent and distributed according to a reverse exponential distribution with parameters α i and β i . The marginal density f i is given by

The density of a 3-dimensional GP distribution with independent reverse exponential generators is then given by

where β = (β 1 , β 2 , β 3 ).

If U = (U 1 , U 2 , U 3 ) with U i , i = 1, 2, 3 independent and distributed according to a reverse Gumbel distribution with parameters α i and β i . The marginal density f i is given by

Calculations very similar to the derivation of Equation (7.2) in [Kiriliouk et al., 2019] show that the resulting density of a 3-dimensional GP distribution with independent reverse exponential generators is then given by

(4.5)

To ensure identifiability, the first component of the location parameter β 1 was fixed, equal to 0. Table 2 gives the estimated parameters for the selected Gumbel model. Table 3 gives AIC, BIC and log-likelihood tests for simplification of the Gumbel model and indicates that the full model (M1) should not be simplified further. For the simulated data presented in Section 5.2, the true parameters for the Gumbel model are known. Figure 3 shows the boxplots of the prediction probabilites obtained from the model with the true parameters. Figure 3 : Boxplots of prediction probabilities of level exceedances for simulated data obtained from the Gumbel model with the true parameters for a) Week 3 ILI incidence rates and epidemic sizes. Outcome is equal to 0 if there was no exceedance and 1 otherwise. The tables give the proportions of 0-outcomes and of 1-outcomes, and the quartiles of the prediction probabilities.

Results from the centers for disease control and prevention's predict the 2013-2014 influenza season challenge

Epidemic influenza-responding to the expected but unpredictable

Univariate and bivariate GPD methods for predicting extreme wind storm losses

Imbalanced Classification with Python: Better Metrics, Balance Skewed Classes, Cost-Sensitive Learning

Risks to patient safety associated with implementation of electronic applications for medication management in ambulatory care-a systematic review

Using extreme value theory approaches to forecast the probability of outbreak of highly pathogenic influenza in Zhejiang

A multivariate extreme value theory approach to anomaly clustering and visualization

An introduction to statistical modeling of extreme values

Sparse vector autoregressive modeling

Modelling extremal events for insurance and finance

Machine learning and extremes for anomaly detection

An extreme value theory approach for the early detection of time clusters. A simulation-based assessment and an illustration to the surveillance of Salmonella

Unsupervised extraction of epidemic syndromes from participatory influenza surveillance self-reported symptoms

Statistics of extremes in hydrology

Health security in 2014: building on preparedness knowledge for emerging health threats

Peaks over thresholds modeling with multivariate generalized Pareto distributions

Forecaster's dilemma: Extreme events and forecast evaluation

Parametric estimation procedures in multivariate generalized Pareto models

Improving regional influenza surveillance through a combination of automated outbreak detection methods: the 2015/16 season in france

The genomic and epidemiological dynamics of human influenza A virus

Data-based comparison of frequency analysis methods: A general framework

Réseau Sentinelles. Inserm/Sorbonne Université

Learning efficient anomaly detectors from k-nn graphs

Multivariate peaks over thresholds models

Multivariate generalized pareto distributions

The precision-recall plot is more informative than the ROC plot when evaluating binary classifiers on imbalanced datasets

Methods for current statistical analysis of excess pneumonia-influenza deaths

Threshold methods for sample extremes

Improving disease incidence estimates in primary care surveillance systems

Assessing the performance of prediction models: a framework for some traditional and novel measures

A relationship between the average precision and the area under the roc curve

Anomaly Detection in Extreme Regions via Empirical MV-sets on the Sphere

Applications of extreme value theory in public health

We thank Tom Britton, Anna Kiriliouk, Andreas Pettersson, Thordis Torainsdottir and Jenny Wadsworth for help and comments, and two referees and an associate editor for comments which have lead to important improvements of the paper.