Adaptively stacking ensembles for influenza forecasting with incomplete data
Thomas McAndrew Nicholas G. Reich
stat.AP cs.LG stat.ML




  

  
    A]Thomas McAndrew[label=e1]mcandrew@umass.edu,
    A]Nick Reich[label=e2]nick@schoolph.umass.edu,

    [A] Department of Biostatistics and Epidemiology,

School of Public Health and Health Sciences

University of Massachusetts Amherst,

Amherst, Massachusetts, United States of America e1

  
    

Seasonal influenza infects between 10 and 50 million people in the United States every year. 
Accurate forecasts of influenza and influenza-like illness (ILI) have been named by the CDC as an important tool to fight the damaging effects of these epidemics. Multi-model ensemble forecasts have shown positive results.
But ensemble forecasts of influenza outbreaks have been static, training on past ILI data, generating a set of optimal component model weights, and keeping the weights constant.
We propose an adaptive ensemble that (i) updates weights week-by-week throughout the influenza season, (ii) only needs the current season's data to make predictions, and (iii) by introducing a prior distribution, shrinks weights toward an equal weighting.
We investigate the prior's ability to impact adaptive ensemble performance and, finding an optimal prior, compare our adaptive ensemble’s performance to equal-weighted and static ensembles.
Applied to forecasts of short-term ILI incidence at the regional and national level, our adaptive model outperforms an equal-weighted ensemble and has similar performance to the static ensemble, needing years of training data.
Adaptive ensembles can quickly train and forecast during epidemics, providing a practical tool to public health officials looking for forecasts that can conform to unique features of a specific season. 



Combination forecasting
Forecast aggregation
Influenza
Statistics
Public health






§ INTRODUCTION


Every year seasonal influenza costs hospitals time and resources, and causes significant loss of life, especially among patients with cardiac disease, previous respiratory illness or allergies, children and the elderly <cit.>.
During the season's peak hospitals admit patients beyond capacity, and in severe cases, for greater lengths of stay, postponing lower priority but still important surgeries <cit.>.
The total economic burden of influenza outbreaks (including medical costs, lost earnings, and loss of life) in the US alone are estimated to exceed $85 billion annually <cit.>.

Accurate forecasts of influenza can help with public health response.
Public health officials at the national, regional, and state level benefit from the advance warning forecasts can provide, by preparing hospitals for increased disease burden <cit.>, shifting vaccines to places in need <cit.>, and alerting the public <cit.>.

Recognizing the value of accurate forecasting, the US Centers for Disease Control and Prevention (CDC) have organized forecasting challenges around seasonal influenza outbreaks <cit.>. 
These challenges focus on generating weekly forecasts of influenza-like illness (ILI) rates at the national and regional level in the US. 
ILI is a syndromic classification and defined as a fever plus an additional symptom such as a cough or sore throat.
Public health officials regard ILI surveillance data as an important indicator of the trends and severity of a given influenza season <cit.>.

The CDC influenza forecasting challenges, called FluSight, have spurred methodological development in the area of influenza forecasting and infectious disease forecasting more broadly.
This effort has led to multiple research groups developing models trained on historical ILI data to forecast future incidence across the US.
These models include statistical and mechanistic models, among other approaches <cit.>.

Decision-makers have requested a single unified forecast, and several groups have developed multi-model ensemble forecasts of influenza. 
These models synthesize predictions from a set of different component models and have shown better predictive ability compared to any single model <cit.>.

Multi-model ensembles aggregate diverse component model predictions into a final `meta-prediction' that is often, but not always, more robust than any single component model <cit.>.
This robustness comes from the ensemble's mixture of component models.
Freed from having to make a single set of assumptions about the data, the ensemble never depends on just one model and associated limitations <cit.>.
In contrast to ensemble modeling approaches like 
bagging <cit.> and boosting <cit.> which create a library of similar component models, multi-model ensembles combine distinct component models together, possibly models designed and implemented in completely different modeling paradigms. Application examples include weather <cit.>, economics<cit.>, energy production<cit.>, education<cit.>, and infectious disease <cit.>. 

The multi-model ensembles addressed in this work are related, but conceptually distinct from Bayesian Model Averaging (BMA) <cit.>.
BMA places a prior over models that are proposed to have generated the observed data <cit.>.
BMA assumes, that from a set of possible models, a single model generated the observed data, and the goal of BMA is to collect data and draw inference about which single model generated the entire data stream.
The multi-model ensembles in this work take a different approach <cit.>.
These models assume that every observation was generated by a probability distribution whose true form is a combination of multiple distributions.
In an ensemble like this, the goal is to determine how the models are combined to from a single data generating distribution—no single model is assumed to generate the entire data stream.

Extending an existing ensemble implementation <cit.>, we developed a new method for combining component models that relies on recently observed in-season data to adaptively estimate a convex combination of models, sometimes referred to as a “linear pool”.
Though past work has studied methods for training ensembles from a fixed data set <cit.>, and continuously adapting ensembles to streaming data <cit.>,
a unique and new challenge in our setting is the reliance on noisy, revision-prone observations like weekly ILI surveillance data.
Recent ILI observations, because they come from a real-time public health surveillance system, are subject to revisions, which are occasionally substantial and can negatively impact forecast accuracy  <cit.>.
Like past work in ensemble forecasting, this work combines a fixed set of probabilistic forecasting models provided by different forecasting teams <cit.> and does not presume to have the capability to refit  component model parameters.

To protect the adaptive multi-model ensemble framework from relying too heavily on recent, revision-prone ILI data, we developed a Bayesian model combination algorithm that uses a prior to regularize ensemble weights.
Previous results in this field show equally-weighted ensembles often perform very well <cit.>, and to align with these results, we chose a uniform Dirichlet prior that shrinks model weights towards equal.
Our prior is also time dependent. 
Unlike a typical model where the prior becomes less influential with increasing data, our prior is designed to exert a constant influence on the final model weights. 
Our model never allows model weights to depend only on revision-prone data.
We show our Bayesian framework can be solved with a Variational Inference algorithm (devi-MM), and demonstrate that an existing Expectation-Maximization approach to ensemble weight estimation <cit.> is a special case.

We compared ensembles that assign equal weights to models (EW) and those with possibly unequal but static and unchanging weights throughout a season (static),  against our adaptive ensemble (adaptive) which updated component model weights every week.
Static ensembles were trained on all cumulative forecasting data before the beginning of each season and were unable to modify their model weights in-season.
In contrast, adaptive ensembles were only trained on within-season data but could modify their weights as each additional week of ILI data was reported.
Comparing static and equal to adaptive ensembles, we highlight: (i) the adaptive model's ability to buffer against sparse and revision-prone data by using a prior over ensemble weights and (ii) similar, and sometimes improved, performance when comparing adaptive to static models, where adaptive models require substantially less data to train.

The manuscript is structured as follows: Section 2 describes the CDC influenza training data and component models included in all ensembles, defines performance metrics, and develops the methodology for our adaptive ensemble. Section 3 presents results, investigating the prior's impact on adaptive ensemble performance and compares the adaptive ensemble to static and equally weighted ensembles. Section 4 relates this work to linear pooling algorithms, and discusses our approach in the broader context of data-driven, evidence-based public health decision making.



§ METHODS




 §.§ Data



  §.§.§ Influenza data

Every week throughout the season, for 10 different reporting regions (HHS1, HHS2, ... , HHS10) and a national average, the CDC publishes surveillance data on the number of patients who visited a sentinel site and the number of patients that were diagnosed with influenza-like illness (ILI).
Percent ILI is defined as the number of patients presenting with a fever (greater than 100F) plus cough or sore throat divided by the number of all patient visits times one hundred.
More than 3,500 outpatient clinics report ILI percentages to the CDC as part of the ILINet surveillance system.
In addition to the newest ILI data, updates to past ILI percentages are also published weekly.

Because of reporting delays, ILI percentages undergo revisions every week after they're first reported, finalized only several weeks after the season ends.
Revisions make predicting future ILI difficult. Component model forecasts need to update week-to-week for two reasons: (i) the newest ILI data is reported, and (ii) all previous ILI values have been revised.
When forecasting models attempt to adjust for this revision process and the uncertainty in reported data they typically outperform models that do not account for data revisions <cit.>.



  §.§.§ Forecasting data


The FluSight Network (FSN) is a collaborative group of influenza forecasters that use historical ILI data to build retrospective forecasts and the performance of models on this retrospective data to build ensemble forecasts.
Teams train their models and generate forecasts as if their model had been applied in real-time across all seasons since 2010/2011.
The resulting collection of forecasts serves as a historical record of how each model would have performed in real-time <cit.>.

Probabilistic forecasts for 21 `component' models were submitted to the FSN for the 2018/2019 season.
The majority of forecasts were submitted retrospectively (models from teams KoT and CU submitted in real-time), as if they were forecasting in real-time and contending with ILI data revisions.

For the CDC FluSight challenges, component forecast models and ensembles built from them, forecast 7 key targets (Figure <ref>). These targets are 1, 2, 3, 4 week-ahead ILI percentages, onset week (defined as the week where 3 consecutive weeks are above the CDC-specified baseline ILI percentage), the season peak week (the epidemic week with highest ILI percentage), and the season peak percentage (the highest ILI percentage throughout the season).
The CDC cannot determine the start of the influenza season and when the peak occurs until mid-summer due to continued revisions to data.
This means final performance data on forecasts—the data used for building ensembles—for the 3 seasonal targets (baseline onset, season peak week and percentage) can not be compiled until after the season closes.
Because our adaptive ensemble is trained on within-season data, only component model forecasts of the 4 week-ahead targets will be used for training from the FluSight Network.

The difficulty estimating ensemble weights comes from ILI data revisions.
An adaptive ensemble must compute weights based on component model performance—the probability component model forecasts place on true ILI values.
But past ILI data is revised weekly and changes component model performance in past epidemic weeks.
Unlike an ensemble that trains on component models scored on finalized ILI data, an adaptive ensemble must account for ILI revisions because they change past component model performance week to week.

Seasons from 2011/2012 up to 2017/2018 will be used to compare equal, static, and adaptive ensembles. The 2010/2011 season will be used as hold-out data to build a prior for the adaptive ensemble.






  §.§.§ Performance metrics
 
A useful predictive model must be well-calibrated and accurate, consistently placing high probability on potential ILI percentages, that after all revisions, become the final ILI percentage.
A proper scoring rule <cit.> widely accepted in the forecasting community is the logscore <cit.>, and defined as the predicted probability placed on a future true (or finally revised) ILI value.
The CDC adopted a modified logscore as their standard scoring metric.
Given a model with probability density function f(z) and true ILI value i*, both the logscore and modified logscore can be written

    logscoreł(f)̊ = logł( ∫_z = i*-w^z=i*+wf(z)  dz )̊,

where w is a small positive number.
The traditional logscore sets w=0% and this is the metric we use to compare ensembles.





In practice, we bin predictive probabilities of ILI percentage from 0 to 13% by 0.1% plus a single bin from 13% to 100%.
The proper logscore reduces to computing the log of the probability given to the single true ILI bin

Models that place high probabilities on the bins containing the observed final ILIs are scored closer to 0 (perfect score), while models that place low probabilities on observed ILI bins are scored closer to -∞ (the worst score).
Logscores are truncated to -10 to follow convention from CDC FluSight scoring <cit.>.



 §.§ Ensemble Framework



  §.§.§ Ensemble model specification

Static and adaptive ensemble algorithms use a mixture model to combine (or average) component models forecasts. 
We developed two approaches to finding ensemble weights: a maximum likelihood method and a maximum aposteriori Bayesian approach. 
A Bayesian approach lets us regularize component model weights by modifying our log likelihood with a prior.
The choice of prior for our adaptive ensemble is meant to discourage erratic shifts when assigning weights to component models as we move, week-by-week through an influenza season.

Our ensemble approach assumes a generative model for our data—ILI percentage at time t (i_t)—by first sampling a component model m with probability π_m, and second, sampling the ILI percentage from the component model m with probability mass function f_m(i).
(In practice, the submitted probabilistic forecasts f_m(i) are probability mass functions defined over discrete bins of possible ILI incidence.
These are often discretized versions of probability densities.)

    z(t) ∼Catł(π_1, π_2, ⋯, π_M)̊
    
      i_t | z(t) ∼ f_z(t)

where z(t) is a vector of length M, each entry corresponding to one of M possible component models, and f_z(t) is the probability density the variable z(t) selected to generate i_t, the ILI percentage at epidemic week t.
The variable z(t) is an indicator variable and assigns the value 1 when this model generated the ILI percentage i_t and 0 otherwise.
While a mixture model is likely not the true generating process for ILI data, this framework conveniently expresses the probability distribution over possible ILI values as a convex combination of component models by averaging over the random variable z.
This also makes clear the difference between combination forecasting and BMA.
Our combination forecast assumes there is a probabilistic distribution over which models produce ILI values at every time step (z(t)).
BMA would place a single probabilistic distribution over which model produced all ILI values.
We define the probability of observing an ILI percentage at time t,

    p(i_t | π) = ∑_z(t) pł[z(t)]̊× pł[i_t | z(t)]̊ = ∑_m=1^Mπ_m f_m(i_t)

and require the πs are non-negative and sum to one.
This mixture model framework for probabilistic model averaging is the standard for the field <cit.>.













  §.§.§ Expectation-Maximization method for weight estimation

Though a convenient way to specify the model, directly optimizing the loglikehood is difficult.
The loglikelihood over all T time points equals

    log p(𝒟 | π ) = ∑_t=1^Tlogł[∑_m=1^Mπ_m f_m(i_t)]̊,

and is a summation over log-sums, where 𝒟 = [i_1, i_2, ⋯, i_T] is the observed T ILI percentages.
An alternative method to optimize the loglikelihood <cit.> considers the loglikelihood over both ILI percentages and the set of hidden indicator variables z(t)—one for each week of data—that decide which component model (from M possible models) was selected to generate data point i_t.
The hidden variables Z simplify the likelihood to

    p(𝒟,Z | π )    = ∏_t=1^T∏_m=1^Mł[π_m f_m(i_t)]̊^z(m,t)

where z(m,t) equals 1 when the m^th model generated the t^th ILI observation and 0 otherwise, and Z is a M × T matrix of all hidden variables.
The Expectation-Maximization (EM) algorithm <cit.> iteratively optimizes a lowerbound on the logliklihood, the expectation over p(𝒟,Z | π ) with respect to p(Z | 𝒟, π^t-1) or E_p(Z | 𝒟, π^t-1) p(𝒟,Z | π ) to find the component model weights that maximize the loglikelihood (deEM-MM Fig. <ref> left).
The prefix de (degeneracy) refers to this algorithm's inability to change parameters inside component models and the suffix MM (mixture model) references the assumed data generating process.

 


  §.§.§ Variational Inference for weight estimation


While the deEM-MM algorithm find an optimal set of weights, it considers π a fixed vector of parameters.
Assuming π is a random variable, we can restructure our problem as Bayesian and infer a posterior probability over weights.
Our Bayesian framework extends (<ref>) by modeling the posterior probability of π

    p(π | Z, 𝒟) ∝ p(π) × p(𝒟,Z | π).

While many different choices are available for a prior over ensemble weights, p(π), we chose a Dirichlet distribution because it is conjugate to the likelihood.
We further decided on a Dirichlet prior that has a single parameter (α) shared across all M models. 
The scale of the shared parameter (α) governs the influence of the prior distribution on the posterior.
A larger value of α shrinks ensemble weights closer towards an equally weighted ensemble.
Equal weights are a reasonable first choice for many ensemble settings.
Smaller values allow the model to rely more heavily on the observed data. 

Typically the prior is fixed during training. 
However, in our setting, where the amount of data is growing and recent past ILI percentages will be revised as we move through a season, we allow the prior to be time-dependent.
A time-dependent prior can act to continually weaken the influence of the (likely to be revised) data on the ensemble weights and consistently pull component model weights towards equal over the course of the season.

We specify the prior as 

    π_t   ∼Dirł[α(t)]̊
    
        p[π_t | α(t)]    = Γł[∑α(t) ]̊∏Γł[α(t)]̊∏_m=1^Mπ_m^α(t)-1

where α(t) is the parameter (the same across all M component models) that defines the Dirichlet distribution at time t.
We chose α(t) to be a constant fraction (ρ) of the number of training data points for fitting at time t divided equally amongst all M component models

    α(t) = ρN(t)M,

where N(t) is the number of training data points, M is the number of component models, and ρ is a value between 0 and 1.
Throughout the season the prior parameter will grow linearly as the ensemble trains on an increasing number of observations and regularize component model weights at a constant rate.

We can plug our specific prior in for time t

    p(π_t | Z, 𝒟)    ∝ pł(π_t)̊× pł(𝒟,Z | π_t)̊
       = Dirł[π_t | α(t)]̊× pł(𝒟,Z | π_t)̊.

and find an optimal set of mixture weights by optimizing the following function of π and Z

    log pł(π_t|Z,𝒟)̊∝∑_t=1^T∑_m=1^M1ł[z(m,t)]̊logł[π_mtf(m,i_t)]̊ + ∑_m=1^Mł[α(t)-1]̊logπ_mt.

We note that the first term on the right-hand side of equation (<ref>) is the same as taking the log of the frequentist expression for the likelihood shown in equation (<ref>).
Intuitively, we expect a stronger (increasing ρ) even prior to encourage weights across components to move towards 1/M, uniformly distributing probability mass to every component model.


To estimate the posterior probability over ensemble weights (<ref>),  We will use a variational inference algorithm (deVI-MM).
Variational inference is similar to the EM algorithm and both can be motivated by introducing a distribution over hidden variables and decomposing the loglikehood into two parts: a lowerbound ℒ on the log likelihood log p(𝒟|π) and an error term that is the difference between the lowerbound and loglikehood ł( ℒ - log pł(𝒟|π)̊)̊.

Given weights π, hidden variables Z, and fixed data 𝒟, we first rewrite the posterior over Z and π

    p(Z,π | 𝒟) =  p(Z, π, 𝒟) p(𝒟)

Reordering terms and taking the log,

    logł[ p(𝒟) ]̊ = logł[  p(𝒟, Z,  π)p(Z, π | 𝒟)]̊ = logł[ p(𝒟 | Z,π) p(Z,π)  p(Z,π | 𝒟)]̊.

We can express the marginal loglikelihood in terms of the complete loglikelihood (numerator) and log of the posterior distribution. 
The last step introduces a distribution over hidden variables Z and π

    logł[ p(𝒟 ) ]̊ = logł[ p(𝒟 | Z,π) p(Z,π) q(Z,π)×q(Z,π) p(Z,π | 𝒟)]̊

and integrates over q(z)

    logł[ p(𝒟) ]̊   = E_qł{logł[ p(𝒟 | Z,π) p(Z,π)q(Z,π)]̊}̊ + KLł[q(Z,π)  p(Z,π| 𝒟) ]̊

where KL is the Kullback-Leibler divergence <cit.>, a non-negative function of q.
Since KL is a non-negative function, the first term on the right hand side is a lower bound to the marginal loglikeihood and the second term, the KL term, is the difference between the marginal loglikelihood and the lower bound.
Iteratively optimizing this lower bound can be shown to monotonically improve the loglikelihood and converge on the best possible representation (q) of the complete loglikelihood. 

The EM and VI algorithms diverge on how to choose q: the EM algorithm considers π a fixed set of parameters and chooses q = p(Z|π,𝒟), zeroing out the above Kullback-Leibler divergence but also assuming p(Z | π, 𝒟) can be computed (see <cit.> for theoretical development and <cit.> for an application to infectious disease forecasting).
The VI algorithm allows any choice for q.
The mean-field approximation, the most common choice for q, factors q into discrete pieces q = ∏_jq_j(z) <cit.> and is the approach we take.
We chose to factor our distribution over hidden variables into a distribution over Z and separate distribution over π, q(π,Z) = q(Z) · q(π), dividing indicator variables from mixture weights.

This choice yields (see details in Suppl. <ref>) a q(π) that follows a Dirichlet distribution and indicator variable z(m,t) that follows a Bernoulli distribution:

    q(π)     ∼Dirł( α(t) + ∑_t=1^T  r(m,t)  )̊
    
       qł[z(m,t)]̊   ∼Bernł[ r(m,t) ]̊,

where the responsibility r(m,t), or the probability that model m generated datapoint i_t, equals

    r(m,t) = expł{E_πlogł[ π_m ]̊ + logł[ f_m(i_t) ]̊}̊∑_m=1^Mexpł{E_πlogł[ π_m ]̊ + logł[ f_m(i_t) ]̊}̊

and the responsibilities summed over component models (m) equals one.
The distribution qł[z(m,t)]̊ computes the probability ILI value i_t was generated by model m.

The Variational approach here can also be interpreted as a generalization of the EM solution.
If we fix π to be a non-random parameter, the responsibilities reduce to

    r(m,t) = π_mf_m(i_t)∑_m=1^Mπ_mf_m(i_t) ,

the same responsibility function for the EM algorithm <cit.>.
Both VI and EM algorithms will find a unique global maximum of the log-likelihood (<ref>) due to the convexity of the log-likelihood function (see analysis of convexity in Suppl. <ref>).



  §.§.§ How prior impacts ensemble weights

The prior can be recognized as a regularizer.
The maximum aposteriori estimate (MAP) for the m^th component model can be computed  computed directly by dividing the weight given to model m by the sum over all weights (sum over m)

    MAPł[π_m]̊ = α(t) + ∑_t=1^Tr(m,t)∑_mα(t) + N(t)

where ∑_m,tr(m,t)=N(t) is the total number of ILI values used for training at time t.
Though useful for computation, an alternative characterization expresses the MAP as a convex combination of the prior and the percent responsability assigned to component model m.
Define the prior percent weight assigned to component model m as α_m(t)/∑_mα_m(t) and the percent responsibility assigned to component model m due to the data as ∑_tr(m,t)/N(t).
The MAP of π can be reexpressed as a convex combination of the prior responsibility plus responsibility learned form the data










    MAPł[π_m]̊ = ł[ α(t)∑_mα(t)]̊ł(ρ1+ρ)̊ + ł(11+ρ)̊ł[∑_t=1^Tr(m,t)N(t)]̊.

Strengthening the prior, in our case increasing ρ, shifts the MAP estimate away from the data-estimated responsibility and towards the prior.
The constant shift, despite increased training data, is a result of our time-dependent prior.
By setting our prior percentage to equal across all models in the ensemble (1/M), the prior weight ρ interpolates between an equally-weighted ensemble (ignoring past model performance) and a completely data-driven weighting.




























  §.§.§ deEM-MM and deVI-MM Algorithms

A comparison of the EM and VI algorithms is shown below in Fig.<ref>.
The EM and VI algorithms are similar to one another: both rely on adding hidden variables to the loglikelihood and iteratively maximizing a lower bound.
When computing Z, both algorithms need the probability that component model m generated data point i for models 1 to M and data points 1 to T. 
A key difference is how the EM and VI algorithms approximate the distribution over hidden variables Z.
The EM algorithm requires the previous estimate of ensemble weights, but the VI algorithm considers the weights a random variable and needs to compute the ensemble weight's expectation (see Suppl. <ref> for details on calculating E[π]).
These differences are present in step 7, updating Z, and step 10 updating π.
Another difference is in evaluating model fit.
The EM algorithm computes the loglikelihood, log of (<ref>), while the VI algorithm computes a related quantity called the Evidence Lower Bound (ELBO).
The ELBO is defined as 

    ELBO = loglikelihood - ł{log[q(π)] + log[q(Z)]}̊,

or the difference between the loglikelihood and our approximation over hidden variables.
When updating the ensemble weights, both algorithms sum over the probability i_t belongs to model m, but the VI algorithm adds an additional term to the ensemble weight, the prior.


[th]8cm
  


[th]8cm
  


figure
  An Expectation maximization (left) and Variational Inference (right) algorithm estimating ensemble weights for a Frequentist and Bayesian adaptive algorithm.
  Although both algorithms have different underlying probability models, they follow similar steps.
  The major differences between EM and VI algorithms are: how Z is computed (step 7), the prior over π (step 10), and how weights are evaluated (steps 11 and 12).
   



 §.§ Experimental design


Five ensemble models will be analyzed in detail (Table. <ref>).
The equally-weighted (EW) ensemble assigns the same weight (1/M) to all component models.
The Static ensemble trains on all past component model forecasts of ILI and assigns weights to component models before the start of the next season.
Weights are kept fixed throughout the next season.
The adaptive model with three types of regularization, that correspond to three values of ρ, will be studied: close to 0% regularization (Adaptive_non), `optimal' regularization (Adaptive_opt), and `over' regularization (Adaptive_over).
We will consider our adaptive model `optimally' regularized by computing logscore statistics for the held-out 2010/2011 season and choosing the prior that has the highest average logscore.





  §.§.§ Training and scoring component models for ensembles


Component models were scored on unrevised ILI data observed throughout each season.
We created a record of the component model scores, as if they had been scored in real time, every week throughout each season.
Adaptive models were trained on component model logscores throughout a season and this impacted how the adaptive model was trained in two ways: (i) every epidemic week new ILI data was observed and generated new component model logscores, and (ii) past ILI data was revised, which in turn changes past component model logscores.
Our static ensemble was trained on component model log scores from past seasons, not the current season.
At the end of a season ILI data is finalized and component model log scores are final, and so training data was constant for the static ensemble.
Ensemble log scores were calculated on ILI percentages reported on EW28, the “final” week used by the CDC FluSight challenges.
ILI values at or after EW28 for each season are considered definitive.
Based on previous work <cit.>, and due to the relatively smaller amount of data available to the adaptive ensemble, we chose to fit a  `constant weight' ensemble model for both the static and adaptive approaches.
The constant weight model assigns one weight to each component model, the weight does not vary by region or target.



  §.§.§ Fitted ensemble models


We computed adaptive ensembles for prior values from 1% to 100% by 1% increments for the 2010/2011 season.
To determine a prespecified `optimal prior', we chose the prior corresponding to the highest average logscore.
The 2010/2011 season was removed from all formal comparisons between ensembles.
The adaptive ensemble corresponding to the optimal prior was used for formal comparisons.

For each season in 2011/2012 through 2017/2018, we computed , , ,  ensembles.
Static ensemble and equally-weighted ensembles were also fit at the beginning of each season, using only data for prior seasons.




  §.§.§ Formal Comparisons


To compare ensemble models, we computed the difference between the logscore assigned to one ensemble versus another.
A random intercept model will describe the difference between logscores, averaged over epidemic weeks and paired by season, region, and target.

    D_w,s,r,t   ∼𝒩ł( β_0 + a_s + b_r + c_t, σ^2)̊
    
        a_s   ∼𝒩ł(0,σ^2_a)̊
    
        b_r   ∼𝒩ł(0,σ^2_b)̊
    
        c_t   ∼𝒩ł(0,σ^2_c)̊

where D is the difference between logscores assigned to two ensembles: for epidemic week w, in season s, region r, and for target t.
The fixed intercept (average difference in logscore between ensemble models) is denoted β_0, and season (a_s), region (b_r), and target (c_t) effects are Normally distributed (𝒩) with corresponding variances.

We will fit two random effects models: one comparing the adaptive vs. equally-weighted ensemble, and the second comparing the adaptive vs. static ensemble.
The conditional mean, 95% confidence interval, and corresponding p-value will be reported.
An additional bootstrapped p-value is also reported.
Random samples (with replacement) are selected from the set of all season-region-target-epidemic week tuples.
For every random sample, a random effects models is fit and conditional means collected for: seasons, regions, and targets. 
The set of random samples is centered to create a null distribution and compared to our original dataset's conditional mean.
We report the empirical probability a random sample from the centered null distribution exceeds our original dataset's conditional mean.



§ RESULTS




 §.§ Choosing prior to maximize performance


We chose an optimal prior for the adaptive model by fitting our adaptive model to the 2010/2011 season for priors from 0% to 100% by 1% and selecting the prior with the highest average logscore (Suppl. <ref> includes adaptive fits and average logscores for all seasons from 2010/2011 to 2017/2018).
The average logscore for the 2010/2011 season peaks at a prior of 8% (Fig. <ref>), and we will consider an adaptive model with 8% prior (ρ=0.08) our ensemble.
After a prior of 8%, the logscore sharply decreases.
This decrease in performance suggest ensemble weights are over regularized, and we chose a 20% prior as our ensemble.
Finally, a 0% prior was chosen as our ensemble.






 §.§ Prior successfully regularizes ensemble weights
 

We compared the , , and ensembles (ρ = 0.01, 0.08, and 0.20 respectively) to investigate how the ensemble weights change with the prior.
Adding a prior regularizes component model weights (Fig. <ref>).
Smaller priors yield higher variability in ensemble weights throughout any given season. 
For example, in 2017/2018 (Fig. <ref>), this is especially evident for the ensemble, when ensemble weights vacillate between assigning almost all weight to one component model or another.
The and weights, in comparison, do not show as high variability over time.
Component model weights for all three ensembles do track one another, with specific model weights moving up and down in unison, albeit in a more muted fashion for the stronger 8% prior () and 20% prior ().
The patterns shown in Fig. <ref> persist in other seasons as well (see Suppl. <ref>).






 §.§ Comparing adaptive vs. equally-weighted and static ensembles





Adaptive ensembles consistently outperform equally-weighted ensembles, and show comparable performance to static ensembles (Fig. <ref>).
The model has higher logscores than either and models, and the EW model.
The , unlike the model, always outperforms the EW model, indicating it is better to over- than under-regularize, at least to some degree.
The static model—trained on multiple seasons of data—performs the best. 
Adaptivity improves over assigning equal weights to component models, and performs similar to, but cannot outperform, the data-rich static model. 




Formal comparisons (see Fig. <ref> and regression table in Suppl. <ref>) demonstrate the ensemble has statistically higher logscores compared to the EW ensembles, and perform similarly to static ensembles.
The ensemble does perform worse against the static models for any particular choice of season, region, and target.
EW and ensembles perform similar for the 2012/2013 season.
The 2015/2016, 2016/2017, and 2017/2018 season does show better (but not significant) performance for the static ensemble.
This suggests that the static model's performance comes from training on multiple seasons of data.
Differences in performance are less variable when stratified by region (see Fig. <ref>B).
Difference in logscores between all ensembles decrease as forecasting horizons increase.
This reflects the difficulty in forecasting further into the future rather than ensemble performance.

In select strata the adaptive ensemble shows a statistically significant higher logscore compared to the EW and close to signficiant difference compared to the static ensemble. 
This significance is not enough to conclude the adaptive ensemble should outperform the static ensemble in most cases, however the data suggest, unsurprisingly, that the static ensemble performs relatively better the more data is has to train on.



 §.§ The adaptive ensemble's ability to learn seasonal variation


The adaptive model has the opportunity to outperform the static model by tailoring component model weights within season.
Training weights based on component model performance data throughout the season is useful, as is shown by the model outperforming the EW model.
But the adaptive ensemble must accrue several weeks of training data before it can perform similar, and in some cases better, than the static ensemble (Fig. <ref>).
When this increase in adaptive performance occurs highlights two key points: (i) the adaptive ensemble is learning seasonal variation and (ii) the static model's better performance early in the season suggests some continuity between component model weights at the end of one season and beginning of the next.




However, the adaptive models also are able to, without any prior data, learn how to create an optimal weighting within a season. 
While the adaptive models may be penalized early on in a season by not possessing the historical knowledge of which models have performed well, they adjust for this by learning which models are performing well in a particular season. 
This is illustrated clearly in Fig. <ref> where the heaviest model at the end of the season (with 46.2% of the weight for the ensemble) was assigned only 4.82% of weight at the beginning of the season. 
This ability to adapt to within-season dynamics appears to provide the model with a slight advantage over the static model in the middle and end of the season.



§ DISCUSSION


We developed a new algorithm for assigning weights to component models.
Our model extends a previous ensemble algorithm <cit.> by assuming component model weights are random variables, rather than fixed quantities.
This model reassigns component model weights every week, but influenza data is prone to change within a season.
For settings with sparse and possibly revisable data, we introduced a time-dependent uniform Dirichlet distribution that regularizes ensemble weights. 
This adaptive ensemble outperforms an equally-weighted ensemble, and performs similarly to a static ensemble that requires multiple years of data to perform well.

It is expected that learning from training data, and the amount available, would contribute to better predictive performance.
The equally-weighted model ignores any training data, and similar to an adaptive ensemble with no training data, assigns equal weights.
The equally-weighted ensemble performs worst.
The next best model, the adaptive ensemble, trains on component model performance within the same season.
Even though the static model has multiple seasons of training data available, and does have the best performance, the adaptive ensemble is not far behind.
It is difficult to tell whether the similar performance is a shortcoming of the static or adaptive ensemble.
Similar performance between ensembles could be because component model performance from past seasons (available to the static ensemble) does not generalize to future seasons.
Alternatively, the adaptive ensemble may not be using the within-season training data efficiently and could be leaving behind hidden patterns in the component model forecasting data. 
But we did find static weights outperforming the adaptive ensemble early on in the season.
These results suggest a model could perhaps use the static ensemble weights, when available, as a prior for the adaptive ensemble.

Regularizing ensemble weights increased adaptive ensemble performance.
Not only can regularization improve adaptive model performance, but could also improve the performance of the static model (Suppl. <ref>).
Improving performance by regularizing ensemble weights with a time-dependent uniform Dirichlet prior may be applicable to any ensemble weighting scheme. 

Though our empirical work suggests an optimal prior near 8%, we expect the optimal prior to vary by season, region, and target.
Different priors will likely be optimal for different biological phenomena too.

A more flexible prior could improve adaptive ensemble performance:
changing the prior percent throughout the season, allowing the prior to depend on performance data from past seasons, or modeling the prior on regional and target information.
Instead of limiting the prior as a regularizer, we can include outside information from past seasons, or patterns in influenza data the ensemble cannot learn from training data.
Future research will explore different methods of modeling prior ensemble weights.

Our work here connects to Gneiting and other's work on linear pooling <cit.>.
Linear pools assume, unlike Bayesian model averaging, that the data is generated by a combination of component models.
Similar to Gneiting's work, our weights are optimized with respect to the loglikelihood of a generative model, but the loglikelihood is really a different representation of a logscore—the log of the probability corresponding to weights π.
Our model could be fit by minimizing a generic loss function based on logscore, but we found the EM (and VI) algorithms fit our model fast and consistently found a global optimum. 

Unlike Gneiting, we did not recalibrate the combined predictive distributions, and methods like the beta transformed linear pool <cit.> or shifted linear pool <cit.> could improve predictive performance even further.

This paper has many limitations future work can address: (i) better choice of prior, (ii) accounting for correlated component models, (iii) post-processing or `recalibrating' our combined forecast.
Our adaptive ensemble examined how a prior can impact ensemble performance, but we only explored a uniform prior.
Future work could study informative priors and how to manipulate priors during the course of a season to improve ensemble performance.
Our model also assumed component models made independent forecasts.
A more advanced ensemble would examine the correlation structure of the data (region, target), and the component model forecasts.
Finally, our ensemble model focused on an optimal set of weights for forecasting, and made no efforts to recalibrate the combined forecast.<cit.>

From a public health perspective, officials have been moving away from stand-alone subjective decision making in favor of incorporating data-driven analyses into decisions.
This trend has become particularly apparent in infectious disease outbreak management.<cit.>
Adaptive ensembles, providing real-time forecasts of infectious disease, supports this trend towards evidence based decision making for public health.
Combining their expertise with statistical models like our adaptive ensemble, public health officials can work towards impeding outbreaks and changing public health policy.
Public health officials often find themselves making decisions in the middle of crises, times when fast effective forecasting is necessary. 
The adaptive ensemble's flexibility, reliance only on near-term data, and capacity to track unusual disease trends can support accelerated public health decision making.

Adaptive ensembles—reweighting component models, week by week, throughout a season—can handle sparse data scenarios, admit new component models season by season, and shows similar prescience compared to the more data-heavy static ensemble.
This added flexibility makes them an important tool for real-time, accurate forecasts.



§ ACKNOWLEDGEMENTS

This work has been supported by the National Institutes of General Medical Sciences (R35GM119582) and the Defense Advanced Research Projects Agency.
The content is solely the responsibility of the authors and does not necessarily represent the official views of NIGMS, the National Institutes of Health, or the Defense Advanced Research Projects Agency.




The data and code used to train equal, static, and adaptive ensembles can be found at <https://github.com/tomcm39/adaptively_stacking_ensembles_for_infleunza_forecasting_with_incomplete_data>.
Influenza-like illness, component model forecasts, and ensemble performance metrics are published on Harvard's Dataverse.
Links to the data are provided on the above GitHub page.






§ COMPUTING Q(Π) AND Q(Z) FOR THE DEGENERATE VARIATIONAL MIXTURE MODEL (DEVI-MM)
 

The functional forms for q(Z) and q(π) are computed.
Readers interested in more details, and theoretical background, should consult <cit.>.
In particular, <cit.> gives a brief introduction to variational inference focused on the applied statistician, while <cit.> provide more theoretical details. 
The most detailed material can be found in <cit.>.

We find the q(Z) and q(π) that maximize the lower bound ℒ(q) by taking advantage of our factored q and using iterated expectations.

    ℒ(q) = E_π,Zł{logł[p(𝒟,Z,π)]̊ - log[q(π)] - log[q(Z)]}̊

Maximizing q(π), we can take the iterated expectation

    max_q(π)ℒ(q) =     max_q(π) E_π|Z E_Zł{logł[p(𝒟,Z,π)]̊ - log[q(π)] - log[q(Z)]}̊
       =max_q(π) E_π|Zł{ E_Zlogł[p(𝒟,Z,π)]̊ - log[q(π)]}̊
    
         =min_q(π)KLł{log q(π) ||  E_Zlogł[p(𝒟,Z,π)]̊}̊

The Kullback-Leibler divergence, taking values from 0 to ∞, is minimized when

    log q(π) = E_Zlogł[p(𝒟,Z,π)]̊

or when 

    q(π) ∝expł{E_Zlogł[p(𝒟,Z,π)]̊}̊.


Maximizing Z follows a similar pattern.
The optimal hidden distributions of q(π) and q(Z) can then be computed

    q(π)    ∝expł{E_Zlogł[p(D,Z,π)]̊}̊
    
      q(Z)    ∝expł{E_πlogł[p(D,Z,π)]̊}̊.


We first compute q(π), expanding the complete loglikelihood, taking the expectation over Z, and recognizing this expectation as a specific distribution.
Here we don't explicitly describe π's dependence on t for convenience.

    log q(π)    ∝E_Zlogł[ p(D,Z,π)]̊
       = ∑_t=1^T∑_m=1^MEł[z(m,t)]̊logł[ π_m f_m(i_t) ]̊ + ∑_m=1^Mł[α(t)-1]̊log(π_m) - logł{ηł[α(t)]̊}̊
       = ∑_t=1^T∑_m=1^MEł[z(m,t)]̊logł( π_m )̊ + ∑_m=1^Mł[α(t)-1]̊logł(π_m)̊
       = ∑_m=1^Mlogł( π_m )̊ł{α(t) + ∑_t=1^TEł[z(m,t)]̊ - 1 }̊
       = ∑_m=1^Mlogł( π_m )̊ł{α(t) + ∑_t=1^T  r(m,t) - 1 }̊,

where η is the normalizing constant for the Dirichlet distribution and r(m,t) the expected value of the indicator variable z_m,t, the probability model m generated the ILI value at time t.
Studying the form of log q(π),we recognize π is Dirichlet distributed

    q(π)    ∼Dirł( γ)̊
    γł[m]̊   = α(t) + ∑_t=1^T  r(m,t)


The same procedure can also be applied to compute q(Z):

    log q(Z)    ∝E_πlogł[ p(D,Z,π)]̊
       =  ∑_t=1^T∑_m=1^M  z(m,t) ł{E_πlogł( π_m )̊ + logł[ f_m(i_t) ]̊}̊ +  E_π∑_m=1^Mł[α(t)-1]̊log(π_m) - logł[η(α)]̊
       ∝∑_t=1^T∑_m=1^M  z(m,t) ł{E_πlogł[ π_m ]̊ + logł[ f_m(i_t) ]̊}̊
    
      q(Z_t,m)    ∝expł{E_πlogł( π_m )̊ + logł[ f_m(i_t) ]̊}̊,

and we recognize q(Z_t,m) is Bernoulli distributed, and with the additional constraint that all indicators must sum to one for every time period (t), we see Z_t is multinomial for every t. 

    q(Z_t,m)    ∼Bernł[ r(m,t) ]̊
    
      r(m,t)    = expł{E_πlogł( π_m )̊ + logł[ f_m(i_t) ]̊}̊∑_m=1^Mexpł{E_πlogł( π_m )̊ + logł[ f_m(i_t) ]̊}̊.


Although burdensome upfront, this approximate procedure drastically reduces computational time, compared to more intense monte carlo sampling techniques, and gives a good approximate answer, only assuming Z and π independent from one another.
Also note this mixture model algorithm (both EM and VI), unlike a typical Gaussian mixture model, cannot manipulate the paramters control the component model distributions.
But this inability to access the component model parameters opens up greater opportunities for other forecasters to submit models, requiring every forecast model to supply just 1, 2 ,3 , and 4 week ahead forecasts.



§ COMPUTING E[LOG(Π)]



Variational inference over hidden variables (Z,π) requires computing the expected value of the log of π's, a natural consequence of adapting the Dirichlet distribution to exponential form.
Exponential form rewrites a probability distribution in the form

    pł[π|η(α)]̊ = h(π)e^∑_kη_k(α) T_k(π) - c(η).

where c(π) is a normalizing constant.
Two facts about exponential form lead to an analytic formula for E(logπ): taking the gradient of c(π) and relating the set of sufficient statistics T(π) with this expectation.
Starting with the first fact.
If c(π) is a normalizing constant, then it must equal

    c(η) = logł( ∫  h(π)e^η^'(α)T(π)  d π)̊,

and it's gradient must equal

    ∇_ηł[c(η)]̊   = ∫ T_k(π) h(π)e^η^'(α)T(π)∫  h(π)e^η^'(α)T(π)   d π
       = ∫ T(π) p(π) = E[T(π)] .

A powerful consequence of exponential form, the gradient of the normalizing constant equals the expected value of the distribution's sufficient statistics.
The Second fact.
Working with the log likelihood of a distribution in exponential form, 

    logł[ ∏_n=1^N p(π|α)]̊ = N log h(π) +  N ∑_kη_k(α_k)T_k(π) - N c(π),

taking the gradient and setting equal to 0,

    ∇_ηlogł[ ∏_n=1^N p(π|α)]̊   =  ∑ T(π_n) - N ∇ c(η) = 0
    ∇ c(η)    = Eł[ T(η) ]̊ = 1N T(π).

The expected value of the distribution's sufficient statistics are equal to the gradient of c(η).
If the Dirichlet's sufficient statistics take the form log(π), then we only need to take the gradient of the normalizing constant to find an analytic expression.

Looking at the Dirichlet distribution, the loglikelihood equals

    logł [ ∏ p(π_n|α) ]̊ = N logΓł(∑_αα_k)̊ - N ∑_klogΓł(α_k)̊ + N ∑_k (α_k-1) log(π_k)

the sufficient statistics for π_k take the form logł( π_k)̊.
Taking the gradient then will provide an analytic expression for computing the logπ's expected value.

    ∇_ηlogł [ ∏ p(π_n|α) ]̊ = N ψł(∑_αα_k)̊ - N ∑_kψł(α_k)̊ + N log(π_k),

where ψ is the digamma function.
Then the expected value of logł(π)̊ equals a difference of digamma functions

    Eł[π_k]̊ = ψł(α_k)̊ -  ψł(∑_αα_k)̊.

This formula can be used to compute the responsibility function (see Algorithm 2 step 7).



§ CONVEXITY ANALYSIS



The problem of finding an optimal set of mixture weights can be recast as a constrained optimization problem.
Attempting to optimize a convex function over a convex set will prove a global optima exists, and showing strict convexity will prove the a unique vector obtains this optimum.

The original problem searches for weights that maximize a loglikelihood whose sum is constrained to equal one,

    max∑_t=1^Tlogł[ ∑_m=1^Mπ_mf_m(i_t) ]̊
       subject to
       ∑_m=1^Mπ_m = 1 ,

and can be converted to a Langragian with a single constraint (λ)

    ℒ(π,λ) = ∑_t=1^Tlogł[ ∑_m=1^Mπ_mf_m(i_t) ]̊  +λł( ∑_m=1^Mπ_m-1 )̊.


Knowing the solution π needs to lie in the M dimensional simplex, a convex set, we only need to show the above function (<ref>) is convex.
Given  (<ref>) is differentiable at least twice, we will appeal to a second-order condition for convexity—proving the Hessian is positive semidefinite.
After proving the Hessian is positive semidefinite, going further and showing the Hessian is in fact positive definite will prove a unique vector π is a the global optimum of our objective function (<ref>).

First, we compute the m^th element of the gradient for (<ref>),

    ∇_π_m f(π) = ∑_t=1^Tf_m(i_i)∑_m=1^Mπ_mf_m(i_t) .

Then the (m,n) entry of the Hessian is

    H(m,n) = - ∑_t=1^T f_m(i_t) f_n(i_t)ł[∑_m=1^Mπ_mf_m(i_t)  ]̊^2.

The Hessian is always negative, and if we minimized the negative loglikelihood (instead of maximizing the positive logelikihood) would see the Hessian is a positive semidefinite matrix.
But we can prove more by rewriting H

    H = F'F,

where the [m,t] element

    F[m,t] = f_m(i_t)∑_m=1^Mπ_mf_m(i_t).

Positive definite matrices can always be written as a transposed matrix times the matrix itself, and so H must also be positive definite.
Our convex optimization problem must then have a global optimum, attained by a unique vector π.
The inability to change forecasting models is a limitation, but allows us to guarantee a global optima.
When the component model parameters can also be updated, no such guarantee exists.



§ SIMPLEX PLOTS OF ALL SEASONS



Regularizing the adaptive ensemble suppresses large swings in weight assignment, independent of season.
The ensemble stays closer to equal weighting throughout the season compared to the and ensembles.
If one exists, there is no strong relationship between the adaptive and static ensemble weighting.
Some seasons (2014/2015 and 2016/2017) show higher variability in weight assignment than others.
The optimal weight assignments given data from all previous seasons (Static MLE) does not carry forward to the optimal weight assignment for the adaptive ensemble at the end of the current season.  







§ TABULATED COMPARISONS OF ADAPTIVE, STATIC, AND EQUAL LOGSCORES



The ensemble consistently outperforms the EW ensemble, and performs similarly to the static ensemble.
When compared to the EW model, the adaptive cannot show statistical significance in the 2012/2013 season, but all differences are in the positive direction, favoring the adaptive model.
The adaptive and static comparisons are more even.
Some season, regions, and targets favor the static ensemble, others the adaptive model, but absolute differences are small.
Despite markedly less data, the adaptive ensemble bests the EW ensemble and shows similar performance to the static.

































§ REGULARIZATION IMPROVES STATIC ENSEMBLE


Regularizing ensemble weights improves both adaptive and static ensemble performance (Fig. <ref>). 
The static ensemble achieves peak performance with a smaller prior than the adaptive.
This smaller prior reflects the data the static model is trained on. 
Past seasons of finalized ILI data. 
The adaptive model, finding peak performance for a larger prior, needs to account for the lack of training data and the revision-prone state of the data mid-season.
























§ SUPPLEMENTARY FIGURE






imsart-number

57




adejo2018predicting
[author]
Adejo, Olugbenga WilsonO. W. Connolly, ThomasT.
(2018).
Predicting student academic performance using multi-model heterogeneous
  ensemble approach.
Journal of Applied Research in Higher Education
10
61–75.



benkeser2018online
[author]
Benkeser, DavidD.,
  Ju, ChengC.,
  Lendle, SamS. van der Laan, MarkM.
(2018).
Online cross-validation-based ensemble learning.
Statistics in medicine
37
249–260.



biggerstaff2016results
[author]
Biggerstaff, MatthewM.,
  Alper, DavidD.,
  Dredze, MarkM.,
  Fox, SpencerS.,
  Fung, Isaac Chun-HaiI. C.-H.,
  Hickmann, Kyle SK. S.,
  Lewis, BryanB.,
  Rosenfeld, RoniR.,
  Shaman, JeffreyJ.,
  Tsou, Ming-HsiangM.-H. et al.
(2016).
Results from the centers for disease control and prevention’s predict
  the 2013–2014 Influenza Season Challenge.
BMC infectious diseases
16
357.



biggerstaff2018results
[author]
Biggerstaff, MatthewM.,
  Johansson, MichaelM.,
  Alper, DavidD.,
  Brooks, Logan CL. C.,
  Chakraborty, PrithwishP.,
  Farrow, David CD. C.,
  Hyun, SangwonS.,
  Kandula, SasikiranS.,
  McGowan, CraigC.,
  Ramakrishnan, NarenN. et al.
(2018).
Results from the second year of a collaborative effort to forecast
  influenza seasons in the United States.
Epidemics
24
26–33.



biggerstaff2017systematic
[author]
Biggerstaff, MatthewM.,
  Kniss, KristaK.,
  Jernigan, Daniel BD. B.,
  Brammer, LynnetteL.,
  Bresee, JosephJ.,
  Garg, ShikhaS.,
  Burns, ErinE. Reed, CarrieC.
(2017).
Systematic assessment of multiple routine and near real-time indicators
  to classify the severity of influenza seasons and pandemics in the United
  States, 2003–2004 through 2015–2016.
American journal of epidemiology
187
1040–1050.



bishop2006pattern
[author]
Bishop, Christopher MC. M.
(2006).
Pattern recognition and machine learning.
springer.



blei2017variational
[author]
Blei, David MD. M.,
  Kucukelbir, AlpA. McAuliffe, Jon DJ. D.
(2017).
Variational inference: A review for statisticians.
Journal of the American Statistical Association
112
859–877.



breiman1996bagging
[author]
Breiman, LeoL.
(1996).
Bagging predictors.
Machine learning
24
123–140.



brooks2015flexible
[author]
Brooks, Logan CL. C.,
  Farrow, David CD. C.,
  Hyun, SangwonS.,
  Tibshirani, Ryan JR. J. Rosenfeld, RoniR.
(2015).
Flexible modeling of epidemics with an empirical Bayes framework.
PLoS computational biology
11
e1004382.



brooks2018nonmechanistic
[author]
Brooks, Logan CL. C.,
  Farrow, David CD. C.,
  Hyun, SangwonS.,
  Tibshirani, Ryan JR. J. Rosenfeld, RoniR.
(2018).
Nonmechanistic forecasts of seasonal influenza with iterative
  one-week-ahead distributions.
PLoS computational biology
14
e1006134.



cover2012elements
[author]
Cover, Thomas MT. M. Thomas, Joy AJ. A.
(2012).
Elements of information theory.
John Wiley & Sons.



dawid2014theory
[author]
Dawid, Alexander PhilipA. P. Musio, MonicaM.
(2014).
Theory and applications of proper scoring rules.
Metron
72
169–183.



dempster1977maximum
[author]
Dempster, Arthur PA. P.,
  Laird, Nan MN. M. Rubin, Donald BD. B.
(1977).
Maximum likelihood from incomplete data via the EM algorithm.
Journal of the royal statistical society. Series B (methodological)
1–38.



dietterich2000ensemble
[author]
Dietterich, Thomas GT. G.
(2000).
Ensemble methods in machine learning.
In International workshop on multiple classifier systems
1–15.
Springer.



fern2003online
[author]
Fern, AlanA. Givan, RobertR.
(2003).
Online ensemble learning: An empirical study.
Machine Learning
53
71–109.



foley2012current
[author]
Foley, Aoife MA. M.,
  Leahy, Paul GP. G.,
  Marvuglia, AntoninoA. McKeogh, Eamon JE. J.
(2012).
Current methods and advances in forecasting of wind power generation.
Renewable Energy
37
1–8.



centers2009prevention
[author]
for Disease Control, CentersC.,
  Prevention et al.
(2009).
Prevention and control of seasonal influenza with vaccines.
  Recommendations of the Advisory Committee on Immunization Practices (ACIP),
  2009.
MMWR Early release
58
1–54.



garratt2011real
[author]
Garratt, AnthonyA.,
  Mitchell, JamesJ.,
  Vahey, Shaun PS. P. Wakerly, Elizabeth CE. C.
(2011).
Real-time inflation forecast densities from ensemble Phillips curves.
The North American Journal of Economics and Finance
22
77–87.



garten2018update
[author]
Garten, RebeccaR.,
  Blanton, LeneeL.,
  Elal, Anwar Isa AbdA. I. A.,
  Alabi, NoreenN.,
  Barnes, JohnJ.,
  Biggerstaff, MatthewM.,
  Brammer, LynnetteL.,
  Budd, Alicia PA. P.,
  Burns, ErinE.,
  Cummings, Charisse NC. N. et al.
(2018).
Update: Influenza Activity in the United States During the 2017–18
  Season and Composition of the 2018–19 Influenza Vaccine.
Morbidity and Mortality Weekly Report
67
634.



germann2006mitigation
[author]
Germann, Timothy CT. C.,
  Kadau, KaiK.,
  Longini, Ira MI. M. Macken, Catherine AC. A.
(2006).
Mitigation strategies for pandemic influenza in the United States.
Proceedings of the National Academy of Sciences
103
5935–5940.



geweke2011optimal
[author]
Geweke, JohnJ. Amisano, GianniG.
(2011).
Optimal prediction pools.
Journal of Econometrics
164
130–141.



gneiting2005weather
[author]
Gneiting, TilmannT. Raftery, Adrian EA. E.
(2005).
Weather forecasting with ensemble methods.
Science
310
248–249.



gneiting2013combining
[author]
Gneiting, TilmannT.,
  Ranjan, RoopeshR. et al.
(2013).
Combining predictive distributions.
Electronic Journal of Statistics
7
1747–1782.



harper2009seasonal
[author]
Harper, Scott AS. A.,
  Bradley, John SJ. S.,
  Englund, Janet AJ. A.,
  File, Thomas MT. M.,
  Gravenstein, StefanS.,
  Hayden, Frederick GF. G.,
  McGeer, Allison JA. J.,
  Neuzil, Kathleen MK. M.,
  Pavia, Andrew TA. T.,
  Tapper, Michael LM. L. et al.
(2009).
Seasonal influenza in adults and children—diagnosis, treatment,
  chemoprophylaxis, and institutional outbreak management: clinical practice
  guidelines of the Infectious Diseases Society of America.
Clinical infectious diseases
1003–1032.



kleiber2011locally
[author]
Kleiber, WilliamW.,
  Raftery, Adrian EA. E.,
  Baars, JeffreyJ.,
  Gneiting, TilmannT.,
  Mass, Clifford FC. F. Grimit, EricE.
(2011).
Locally calibrated probabilistic temperature forecasting using
  geostatistical model averaging and local Bayesian model averaging.
Monthly Weather Review
139
2630–2649.



krishnamurti1999improved
[author]
Krishnamurti, TNT.,
  Kishtawal, CMC.,
  LaRow, Timothy ET. E.,
  Bachiochi, David RD. R.,
  Zhang, ZhanZ.,
  Williford, C EricC. E.,
  Gadgil, SulochanaS. Surendran, SajaniS.
(1999).
Improved weather and seasonal climate forecasts from multimodel
  superensemble.
Science
285
1548–1550.



krumkamp2011health
[author]
Krumkamp, RR.,
  Kretzschmar, MM.,
  Rudge, JWJ.,
  Ahmad, AA.,
  Hanvoravongchai, PP.,
  Westenhoefer, JJ.,
  Stein, MM.,
  Putthasri, WW. Coker, RR.
(2011).
Health service resource needs for pandemic influenza in developing
  countries: a linked transmission dynamics, interventions and resource demand
  model.
Epidemiology & Infection
139
59–67.



madigan1996bayesian
[author]
Madigan, DavidD.,
  Raftery, Adrian EA. E.,
  Volinsky, CC. Hoeting, JJ.
(1996).
Bayesian model averaging.
In Proceedings of the AAAI Workshop on Integrating Multiple Learned
  Models, Portland, OR
77–83.



mcgowan2019collaborative
[author]
McGowan, Craig JC. J.,
  Biggerstaff, MatthewM.,
  Johansson, MichaelM.,
  Apfeldorf, Karyn MK. M.,
  Ben-Nun, MichalM.,
  Brooks, LoganL.,
  Convertino, MatteoM.,
  Erraguntla, MadhavM.,
  Farrow, David CD. C.,
  Freeze, JohnJ. et al.
(2019).
Collaborative efforts to forecast seasonal influenza in the United
  States, 2015–2016.
Scientific reports
9
683.



minka2000bayesian
[author]
Minka, Thomas PT. P.
(2000).
Bayesian model averaging is not model combination.
Available electronically at http://www. stat. cmu.
  edu/minka/papers/bma. html
1–2.



molinari2007annual
[author]
Molinari, Noelle-Angelique MN.-A. M.,
  Ortega-Sanchez, Ismael RI. R.,
  Messonnier, Mark LM. L.,
  Thompson, William WW. W.,
  Wortley, Pascale MP. M.,
  Weintraub, EricE. Bridges, Carolyn BC. B.
(2007).
The annual impact of seasonal influenza in the US: measuring disease
  burden and costs.
Vaccine
25
5086–5096.



montgomery2010bayesian
[author]
Montgomery, Jacob MJ. M. Nyhan, BrendanB.
(2010).
Bayesian model averaging: Theoretical developments and practical
  applications.
Political Analysis
18
245–270.



murphy2012machine
[author]
Murphy, KK.
(2012).
Machine learning: a probabilistic approach.
Massachusetts Institute of Technology
1–21.



osthus2019even
[author]
Osthus, DaveD.,
  Daughton, Ashlynn RA. R. Priedhorsky, ReidR.
(2019).
Even a good influenza forecasting model can benefit from internet-based
  nowcasts, but those benefits are limited.
PLoS computational biology
15
e1006599.



osthus2018dynamic
[author]
Osthus, DaveD.,
  Gattiker, JamesJ.,
  Priedhorsky, ReidR.,
  Del Valle, Sara YS. Y. et al.
(2018).
Dynamic Bayesian Influenza Forecasting in the United States with
  Hierarchical Discrepancy.
Bayesian Analysis.



pari2018multi
[author]
Pari, RR.,
  Sandhya, MM. Sankar, SharmilaS.
(2018).
A Multi-Tier Stacked Ensemble Algorithm to Reduce the Regret of
  Incremental Learning for Streaming Data.
IEEE Access
6
48726–48739.



pierro2016multi
[author]
Pierro, MarcoM.,
  Bucci, FrancescoF.,
  De Felice, MatteoM.,
  Maggioni, EnricoE.,
  Moser, DavidD.,
  Perotto, AlessandroA.,
  Spada, FrancescoF. Cornaro, CristinaC.
(2016).
Multi-Model Ensemble for day ahead prediction of photovoltaic power
  generation.
Solar energy
134
132–146.



raftery2005using
[author]
Raftery, Adrian EA. E.,
  Gneiting, TilmannT.,
  Balabdaoui, FadouaF. Polakowski, MichaelM.
(2005).
Using Bayesian model averaging to calibrate forecast ensembles.
Monthly weather review
133
1155–1174.



ranjan2010combining
[author]
Ranjan, RoopeshR. Gneiting, TilmannT.
(2010).
Combining probability forecasts.
Journal of the Royal Statistical Society: Series B (Statistical
  Methodology)
72
71–91.



ray2018prediction
[author]
Ray, Evan LE. L. Reich, Nicholas GN. G.
(2018).
Prediction of infectious disease epidemics via weighted density
  ensembles.
PLoS computational biology
14
e1005910.



ray2017infectious
[author]
Ray, Evan LE. L.,
  Sakrejda, KrzysztofK.,
  Lauer, Stephen AS. A.,
  Johansson, Michael AM. A. Reich, Nicholas GN. G.
(2017).
Infectious disease prediction with kernel conditional density
  estimation.
Statistics in medicine
36
4908–4929.



reed2009infection
[author]
Reed, DeoineD. Kemmerly, Sandra AS. A.
(2009).
Infection control and prevention: a review of hospital-acquired
  infections and the economic implications.
The Ochsner Journal
9
27–31.



reich2019collaborativepnas
[author]
Reich, Nicholas GN. G.,
  Brooks, Logan CL. C.,
  Fox, Spencer JS. J.,
  Kandula, SasikiranS.,
  McGowan, Craig JC. J.,
  Moore, EvanE.,
  Osthus, DaveD.,
  Ray, Evan LE. L.,
  Tushar, AbhinavA.,
  Yamana, Teresa KT. K. et al.
(2019).
A collaborative multiyear, multimodel assessment of seasonal influenza
  forecasting in the United States.
Proceedings of the National Academy of Sciences
116
3146–3154.



reich2019collaborative
[author]
Reich, Nicholas GN. G.,
  McGowan, Craig JC. J.,
  Yamana, Teresa KT. K.,
  Tushar, AbhinavA.,
  Ray, Evan LE. L.,
  Osthus, DaveD.,
  Kandula, SasikiranS.,
  Brooks, Logan CL. C.,
  Crawford-Crudell, WillowW.,
  Gibson, Graham CaseyG. C. et al.
(2019).
A Collaborative Multi-Model Ensemble for Real-Time Influenza Season
  Forecasting in the US.
bioRxiv
566604.



Rivers2019
[author]
Rivers, CaitlinC.,
  Chretien, Jean-PaulJ.-P.,
  Riley, StevenS.,
  Pavlin, Julie AJ. A.,
  Woodward, AlexandraA.,
  Brett-Major, DavidD.,
  Maljkovic Berry, IrinaI.,
  Morton, LindsayL.,
  Jarman, Richard GR. G.,
  Biggerstaff, MatthewM.,
  Johansson, Michael AM. A.,
  Reich, Nicholas GN. G.,
  Meyer, DianeD.,
  Snyder, Michael RM. R. Pollett, SimonS.
(2019).
Using “outbreak science” to strengthen the use of models during
  epidemics.
Nature Communications
10
3102.
10.1038/s41467-019-11067-2



russell2008global
[author]
Russell, Colin AC. A.,
  Jones, Terry CT. C.,
  Barr, Ian GI. G.,
  Cox, Nancy JN. J.,
  Garten, Rebecca JR. J.,
  Gregory, VickyV.,
  Gust, Ian DI. D.,
  Hampson, Alan WA. W.,
  Hay, Alan JA. J.,
  Hurt, Aeron CA. C. et al.
(2008).
The global circulation of seasonal influenza A (H3N2) viruses.
Science
320
340–346.



rustagi1976variational
[author]
Rustagi, Jagdish SJ. S.
(1976).
Variational methods in statistics.
In Variational methods in statistics
Academic Press.



schanzer2013impact
[author]
Schanzer, Dena LD. L. Schwartz, BrianB.
(2013).
Impact of seasonal and pandemic influenza on emergency department
  visits, 2003–2010, Ontario, Canada.
Academic Emergency Medicine
20
388–397.



schapire2003boosting
[author]
Schapire, Robert ER. E.
(2003).
The boosting approach to machine learning: An overview.
In Nonlinear estimation and classification
149–171.
Springer.



sewell2008ensemble
[author]
Sewell, MartinM.
(2008).
Ensemble learning.
RN
11.



shaman2012forecasting
[author]
Shaman, JeffreyJ. Karspeck, AliciaA.
(2012).
Forecasting seasonal outbreaks of influenza.
Proceedings of the National Academy of Sciences
109
20425–20430.



steel2011bayesian
[author]
Steel, Mark FJM. F.
(2011).
Bayesian model averaging and forecasting.
Bulletin of EU and US Inflation and Macroeconomic Analysis
200
30–41.



van2007super
[author]
Van der Laan, Mark JM. J.,
  Polley, Eric CE. C. Hubbard, Alan EA. E.
(2007).
Super learner.
Statistical applications in genetics and molecular biology
6.



vaughan2009effective
[author]
Vaughan, ElaineE. Tinker, TimothyT.
(2009).
Effective health risk communication about pandemic influenza for
  vulnerable populations.
American Journal of Public Health
99
S324–S332.



wallis2011combining
[author]
Wallis, Kenneth FK. F.
(2011).
Combining forecasts–forty years later.
Applied Financial Economics
21
33–41.



yamana2017individual
[author]
Yamana, Teresa KT. K.,
  Kandula, SasikiranS. Shaman, JeffreyJ.
(2017).
Individual versus superensemble forecasts of seasonal influenza
  outbreaks in the United States.
PLoS computational biology
13
e1005801.



zhou2012ensemble
[author]
Zhou, Zhi-HuaZ.-H.
(2012).
Ensemble methods: foundations and algorithms.
Chapman and Hall/CRC.








