#  tsbootstrap: Enhancing Time Series Analysis with Advanced Bootstrapping
Techniques

Sankalp Gilda<EMAIL_ADDRESS>
ML Collective
San Francisco, CA 98195-4322, USA Benedikt Heidrich
<EMAIL_ADDRESS>
sktime project Franz Kiraly<EMAIL_ADDRESS>
sktime project

###### Abstract

In time series analysis, traditional bootstrapping methods often fall short
due to their assumption of data independence, a condition rarely met in time-
dependent data. This paper introduces tsbootstrap, a python package designed
specifically to address this challenge. It offers a comprehensive suite of
bootstrapping techniques, including Block, Residual, and advanced methods like
Markov and Sieve Bootstraps, each tailored to respect the temporal
dependencies in time series data. This framework not only enhances the
accuracy of uncertainty estimation in time series analysis but also integrates
seamlessly with the existing python data science ecosystem, making it an
invaluable asset for researchers and practitioners in various fields.

Keywords: Time Series, Bootstrapping, Resampling, Uncertainty Quantification

##  1 Introduction

Time series data, characterized by their sequential observations, are the
backbone of critical decision-making in various scientific and economic
fields. Unlike static or cross-sectional datasets, time series data are
dynamic, evolving over time and exhibiting intricate dependencies that reflect
the complexities of the real world. This dynamic nature poses unique
challenges, particularly when it comes to predicting future events
(forecasting) and assessing the reliability of these predictions (uncertainty
quantification). Traditional statistical methods, designed for independent and
identically distributed (IID) data, falter in this context, as they overlook
the temporal correlations that are fundamental to time series data. This gap
in the statistical toolkit highlights a pressing need for methodologies that
can adeptly navigate the intricacies of time series, preserving their inherent
structure while facilitating rigorous analysis and inference.

Bootstrapping (Efron, 1992; Chernick, 2011) offers a powerful tool for
estimating the variability of a statistic without making stringent assumptions
about the underlying data distribution. However, the classic IID bootstrapping
approach, which resamples data points as if they were independent, is ill-
suited for time series where each observation is a chapter in an ongoing
story, not an isolated event. To capture the essence of time series,
bootstrapping must adapt, preserving the internal structure and dependencies
within the data.

Addressing this need, various time series bootstrapping methods have emerged,
each designed to respect the chronological integrity of the data while
enabling robust uncertainty quantification (Kreiss and Lahiri, 2012). Yet, the
software implementation of these methods presents its own set of challenges.
They require not just statistical acumen but also computational finesse,
demanding algorithms that can handle the complexity of time series data while
being accessible to the broader scientific and industrial community.

tsbootstrap (Gilda, 2024) emerges as a response to the complexities and
specific needs in time series analysis. It is a python package designed to
provide a range of bootstrapping methods tailored for time series data. This
package includes traditional techniques such as Block Bootstrap (Kunsch, 1989)
and extends to more sophisticated methods like Markov and Sieve Bootstrap
(Bühlmann, 1997). Each method within tsbootstrap is implemented to respect and
preserve the chronological order and inherent correlations of time series
data, ensuring more accurate and reliable bootstrapping.

We have designed tsbootstrap as a versatile and user-friendly tool to empower
both academic researchers and industry practitioners, offering them new
capabilities to navigate and interpret the complexities of time series data
with unprecedented precision and insight. Specifically, we fill in two
existing gaps:

  * •

Framework for Enhanced Uncertainty Quantification: tsbootstrap provides a
robust framework that significantly advances the practice of uncertainty
quantification in time series analysis. This is particularly crucial in fields
such as finance, where accurate risk assessment is essential; in meteorology,
where predictive accuracy can have wide-ranging implications; and in
epidemiology, where understanding uncertainty can inform public health
decisions.

  * •

Seamless Integration with the Time Series Analysis Ecosystem: The package is
meticulously designed to complement and enhance the existing python data
science ecosystem. By “time series analysis ecosystem”, we refer to the
collection of tools, libraries, and methodologies that are commonly employed
in the analysis of time series data. This includes integration with popular
libraries like sktime (sktime developers, 2024) and pandas (pandas development
team, 2024), ensuring that tsbootstrap is not only a standalone tool but also
a synergistic component that enhances the capabilities of other time series
analysis tools.

In the rest of this paper, we delve into tsbootstrap in more depth, offering a
comprehensive examination of its theoretical underpinnings, implementation
details, and the breadth of its application. We demonstrate its utility in
elevating time series analysis through case studies and comparative analyses,
emphasizing its role as an essential instrument for researchers and
practitioners. In subsequent sections we provide a structured exploration. In
Section 2 we introduce bootstrapping and its significance. In Section 3, we
focuse on tsbootstrap’s design principles and functionalities. In Section 4 we
detail the supported bootstrap algorithms. In Section 5, we illustrate
practical usage, integration with sktime, and extension capabilities. Finally,
we conclude in Section 6, and outline the future trajectory of tsbootstrap in
Section 7.

Bootstrapping Method |  Class (…Bootstrap) |  Salient Features
---|---|---
Moving Block |  MovingBlock |  Suitable for series with significant intra-block dependencies; avoid for weak dependencies
Circular Block |  CircularBlock |  Ideal for seasonal/cyclical data; not recommended for non-cyclical data
Stationary |  StationaryBlock |  Best for data with varying dependency lengths; less effective for uniform dependency lengths
Non-Overlapping Block |  NonOverlappingBlock |  Effective for data with distinct segments; not suitable for continuous or highly dependent data
Tapered Block |  Bartletts, Hamming, Blackman, Tukey |  Mitigates edge effects using window functions; requires matching window functions to data characteristics
Residual |  WholeResidual, BlockResidual |  Evaluates model-based uncertainty; not intended for non-model analyses
Statistic-Preserving |  WholeStatisticPreserving, BlockStatisticPreserving |  Maintains key statistical properties; not for general purposes without specific preservation needs
Distribution |  WholeDistribution, BlockDistribution |  Leverages known distributions; inappropriate for unknown distributions
Markov |  WholeMarkov, BlockMarkov |  Suitable for Markovian time series; avoid for non-Markovian series
Sieve |  WholeSieve, BlockSieve |  Designed for autoregressive models; not recommended for non-autoregressive series
Table 1: Overview of bootstrapping classes in tsbootstrap.

##  2 Bootstrapping

Bootstrapping is a resampling technique fundamental to statistics, providing a
way to estimate the distribution of a sample statistic (like the mean or
variance) by resampling with replacement from the data. It allows
statisticians to assess the reliability of inferential statistics, such as
confidence intervals and significance tests, without relying on stringent
assumptions about the underlying population distribution.

At its core, bootstrapping involves repeatedly drawing samples, typically
thousands of times, from a single observed dataset. Each sample is drawn with
replacement, meaning the same data point can appear multiple times in a
bootstrap sample. By applying the statistical measure of interest to each
resample, a distribution of these statistics is generated, offering insight
into the variability and bias of the estimate.

To underscore the importance of bootstrapping in various statistical
processes, consider the following key aspects:

  * •

Model-Free Inference: Bootstrapping does not assume a specific parametric form
for the data distribution, making it a non-parametric approach. This
flexibility is crucial when the underlying distribution is unknown or complex,
allowing statisticians to avoid potentially incorrect model assumptions.

  * •

Assessment of Uncertainty: It provides a straightforward method to estimate
the confidence intervals and standard errors of estimates, which are vital for
understanding the precision of statistical inferences.

  * •

Validation of Models: In model building, bootstrapping can be used to validate
the stability and reliability of predictive models. It helps in assessing the
model’s variance and bias, offering insights into its generalizability.

  * •

Complex Estimators: Bootstrapping is particularly valuable for complex
estimators where the theoretical distribution is difficult or impossible to
derive. It allows for the empirical examination of the estimator’s
distribution.

Transitioning to the specific challenges of time series data, it’s crucial to
recognize that the inherent dependencies within these data sets necessitate a
more nuanced approach to bootstrapping. Traditional IID bootstrapping
techniques, which treat each data point as independent, are not suitable for
time series data due to their sequential nature (Lahiri, 2013). This
misalignment highlights the need for specialized time series bootstrapping
methods that respect the data’s temporal structure, ensuring the integrity of
the analysis.

##  3 tsbootstrap

Figure 1: Optional (soft) and required (hard) dependencies for tsbootstrap
v0.1.0. The latter are indicated via lime-colored folders and dashed arrows.

###  3.1 Design Principles

  1. 1.

Modularity and Extensibility: The package’s architecture is highly modular,
allowing users to easily incorporate custom bootstrapping algorithms. This
flexibility makes it an adaptable tool for a wide range of time series data
applications. tsbootstrap’s design embraces the strategy pattern and
composition, ensuring that new bootstrapping strategies can be easily added,
promoting a collaborative environment for the community to contribute and
extend the library.

  2. 2.

User-Centric API: The API is crafted with user experience in mind, ensuring
simplicity in performing complex bootstrapping tasks. This design principle is
evident in the package’s comprehensive documentation and illustrative examples
provided on the GitHub repository. Following scikit-learn-like conventions,
tsbootstrap provides a consistent API , making it accessible and user-friendly
for those accustomed to scikit-learn’s patterns.

  3. 3.

Integration with the python Ecosystem: tsbootstrap is built to seamlessly
integrate with popular python libraries like numpy and pandas, enhancing its
utility in the prevalent data science ecosystem. Moreover, its design
facilitates easy integration with other python libraries, particularly sktime,
allowing tsbootstrap to be used within existing pipelines, enhancing the
robustness and reliability of time series analysis.

  4. 4.

Alignment with scikit-learn’s Cross-Validation: Mirroring scikit-learn’s
cross-validation class, tsbootstrap makes a slight adjustment: it focuses on
bootstrapping for time series data, replacing the split method with a
bootstrap method. This approach is particularly suited for time series data,
respecting and preserving their temporal structure during the resampling
process, and providing more reliable estimates of model performance.

###  3.2 Core Functionalities

  1. 1.

Diverse Bootstrapping Methods: We include various methods such as Block
(Kunsch, 1989; Politis and Romano, 1994; Paparoditis and Politis, 2001),
Residual, Statistic-Preserving, and Distribution Bootstrap. Each method is
fine-tuned for specific types of time series data, offering users a tailored
approach to their analytical needs.

  2. 2.

Customization and Flexibility: Users can customize various parameters like
block size, sampling methods, and random seed settings. This level of
customization is instrumental in conducting precise and controlled
bootstrapping analyses.

  3. 3.

Efficient Data Handling and Processing: We employ advanced data handling
techniques, ensuring efficient processing of large time series datasets, input
data and parameter validation, and fine-grained logging.

###  3.3 Integration into the python ML Universe

tsbootstrap is designed to integrate seamlessly with established and widely
used libraries like sklearn (Pedregosa et al., 2011; Grisel et al., 2024) and
sktime (sktime developers, 2024).

tsbootstrap adopts the task-specific unified interface principle, and
composable specification language. It also provides task independent interface
points common to sklearn and sktime, such as get_params, set_params,
get_fitted_params.

This high-level interoperability between tsbootstrap, sklearn, and adjacent
framworks such as sktime is mediated by the skbase library (RNKuhns et al.,
2024), providing shared interface points and base classes. This API
integration also allows tsbootstrap to leverage sktime’s handling of time
series specific issues such as forecasting and time-dependent splitting.

Interoperability is in-principle possible both ways, and methodologically
advantageous: sktime can leverage tsbootstrap to make, or improve,
probabilistic forecasts, via custom bootstrapping schemas. Conversely,
tsbootstrap can use sktime forecasters or transformations to construct
bootstrap schemas with specific model assumptions, e.g., forecast or model
residual bootstraps.

Finally, this design of tsbootstrap emphasizes extensibility and community
engagement, encouraging contributions that expand its range of bootstrapping
techniques and enhance its integration with the broader python data science
landscape. Using an extensible, sklearn-like interface mediated by skbase,
users can easily create interface compatible bootstrap algorithms, to add to
the core library, or in their own code bases. Details on extension are
presented in Section 4.3.

##  4 Supported Bootstrap Algorithms

tsbootstrap offers an extensive array of bootstrapping techniques, specially
designed for time series data to preserve their inherent temporal structures
during the resampling process. Below we provide a detailed overview of the
supported bootstrap algorithms within the package:

Figure 2: Top left: BaseTimeSeriesBootstrap is the base-class, and has three
main child classes, indicated in Orange. Top right, middle, and bottom: Class
structure for the three main child classes depicted from the top-left figure.
gray blocks/classes are helper classes that are not exposed to the user. All
classes/blocks in green, orange, and blue have Bootstrap as the suffix in the
package.

###  4.1 Block Bootstrap Methods

For time series bootstrapping, Block bootstrap methods stand out as critical
tools. These methods are indispensable because they allow the preservation of
dependency structures within the data, a fundamental aspect that distinguishes
time series from cross-sectional data. By resampling blocks of observations
rather than individual data points, these methods maintain the chronological
order and inherent dependencies, which are crucial for accurate time series
analysis.

Here is a deeper look into the Block Bootstrap methods provided by
tsbootstrap:

  * •

Moving Block Bootstrap (MBB): This method is particularly effective for time
series with significant autocorrelation within contiguous segments. In MBB,
overlapping blocks of data are resampled, ensuring preservation of the
temporal dependence. This technique is especially useful in analyzing
financial time series, where capturing the short-term dependencies is crucial
for accurate forecasting and risk assessment.

  * •

Circular Block Bootstrap (CBB): CBB extends the concept of MBB by treating the
data as if it were circular, meaning that the end of the series connects back
to the beginning. This method is particularly beneficial for handling seasonal
or cyclical time series, as it allows the bootstrapping process to maintain
the inherent cyclical properties, such as those found in meteorological or
economic seasonal patterns.

  * •

Stationary Block Bootstrap (SBB): This approach introduces a random element to
the block length, providing a more dynamic resampling process. Unlike MBB or
CBB, where the block size is fixed, Stationary Bootstrap (Politis and Romano,
1994) randomly varies the length of each block according to a predetermined
probability distribution. This method is adept at capturing both short- and
long-term dependencies within the time series, offering a versatile approach
that can adapt to various structures within the data.

  * •

Non-Overlapping Block Bootstrap (NBB): NBB is a simpler variant where blocks
of data are sampled in a non-overlapping manner. This method ensures that each
data point is included only once in each bootstrap sample, reducing
redundancy. NBB can be particularly effective for time series data with clear
and distinct segments where the independence between blocks can be reasonably
assumed.

  * •

Tapered Block Bootstrap (TBB): TBB (Paparoditis and Politis, 2001) is an
advanced variant that incorporates window functions to taper the blocks,
reducing the potential discontinuities at the block edges. This method can be
applied on top of the existing block bootstrap methods – such as MBB, CBB, SB,
and NBB – as a base. The tapering process involves applying window functions
like Hamming, Hanning, Blackman, or Tukey to the data blocks, which helps
mitigate edge effects and improve the representativeness of the bootstrap
sample. This technique is particularly useful in scenarios where the edge
effects might introduce bias or when the time series exhibits strong seasonal
components or trends. By employing different window functions, users can fine-
tune the tapering to match the specific characteristics and requirements of
their time series data, enhancing the flexibility and effectiveness of the
bootstrapping process.

###  4.2 Advanced Bootstrap Methods

tsbootstrap includes a comprehensive suite of time series bootstrapping
techniques, each designed for specific types of time series data:

  * •

Residual Bootstrap: This method is particularly adept at handling model-based
uncertainties (Kreiss and Lahiri, 2012; Imon and Ali, 2005). After fitting a
time series model to the data, the residuals (the differences between the
observed values and the model’s predictions) are bootstrapped. This approach
allows for the assessment of the variability and uncertainty of the model
predictions, accommodating a range of models from simple linear regressions to
complex ARIMA models. The Residual Bootstrap is invaluable for validating the
stability and accuracy of the predictions made by these time series models.

  * •

Statistic-Preserving Bootstrap: This technique ensures that the bootstrapped
samples retain a user-defined key statistical property of the original
dataset, be it mean, variance, or something else. It is particularly
beneficial in applications where maintaining the original data’s
distributional characteristics is crucial for accurate analysis and inference.
This method is essential for studies where the underlying distributional
properties play a critical role in the analysis outcomes or model validations.

  * •

Distribution Bootstrap: When the underlying distribution of the time series is
known or can be estimated, the Distribution Bootstrap method leverages this
information to generate resampled data. This approach is ideal for datasets
where specific distributional assumptions hold true, allowing for more precise
and theoretically grounded inferences; in situations where the time series is
believed to follow a known distributional pattern, aiding in more accurate
modeling and forecasting.

  * •

Markov Bootstrap: Tailored for time series exhibiting Markovian properties,
this method is suitable for data where the future state is dependent only on
the current state and not on the sequence of events that preceded it. Commonly
used in financial time series analysis, the Markov Bootstrap (Härdle et al.,
2003) can capture the stochastic processes underlying the financial markets,
offering insights into the dynamics and potential future behaviors of these
markets.

  * •

Sieve Bootstrap: This method is specifically designed for time series that can
be modeled using autoregressive structures. By active on the coefficients of
an autoregressive model, the Sieve Bootstrap provides a robust framework for
analyzing and forecasting time-dependent data. This approach is particularly
powerful in scenarios where the data exhibits autoregressive behavior,
enabling analysts to make informed forecasts and understand the underlying
dynamics of the time series (Bühlmann, 1997).

Each of these methods demonstrates tsbootstrap’s capability to address the
unique challenges presented by different types of time series data, making it
a versatile tool in the field of time series analysis.

###  4.3 Extension Template

tsbootstrap is meant to be extensible, designed to make it easy for power
users to add new bootstrap algorithms, to the main code base, or (potentially
closed or differently licensed) third party extensions.

For this, tsbootstrap follows the combined strategy/template design pattern
also seen in sktime or sklearn, using the foundations of skbase for extension
and testing:

  * •

Extension template: implementers of new algorithms are provided with an
extension template, i.e., a fill-in python file with step-by-step instructions
on filling in a new algorithm.

  * •

Extension contract: boilerplate is abstracted away in the extension template -
only a private _bootstrap needs to be implemented, with simplified inputs,
while class inheritance and boilerplate layering ensure public interface
compliance.

  * •

Contract checker: extenders are provided with a compliance checker method,
check_estimator, which allows testing at runtime, and can be used in 3rd party
testing setups for continuous API compliance checks.

Further details, including a 90 minute general-purpose tutorial to extending
skbase-templated packages, can be found in the skbase documentation (RNKuhns
et al., 2024).

##  5 Usage and Integration with sktime

###  5.1 Basic Usage

In Figure 1, we detail the various building blocks and user-facing classes in
tsbootstrap. Below, in Example Code 1 and Figure 1, we demonstrate a simple
usage of MovingBlockBootstrap using simulated data.

[htbp] from tsbootstrap import MovingBlockBootstrap import numpy as np import
matplotlib.pyplot as plt

# Create custom time series data. While below is for univariate data, #
”tsbootstrap” can handle multivariate time series as well. n_samples = 10 X =
np.arange(n_samples)

# Instantiate the bootstrap object n_bootstraps = 3 block_length = 3 rng = 42
mbb = MovingBlockBootstrap( n_bootstraps=n_bootstraps, rng=rng,
block_length=block_length )

# Generate bootstrapped samples return_indices = False bootstrapped_samples =
mbb.bootstrap( X, return_indices=return_indices)

# Collect bootstrap samples X_bootstrapped = [] for data in
bootstrapped_samples: X_bootstrapped.append(data)

X_bootstrapped = np.array(X_bootstrapped)

# Plot the bootstrapped samples and the original time series plt.figure()
plt.plot(X, label=”Original Time Series”, color=”black”) for i in
range(n_bootstraps): _= plt.plot(X_bootstrapped[i].reshape(-1,),
label=f”Bootstrapped Sample i+1”, ls=”–”)

plt.legend() plt.xlabel(”Time”) plt.ylabel(”Value”) plt.show()

Example Code 1:  Example code snippet demonstrating the use of tsbootstrap
Figure 3: Output of Example Code 1.

###  5.2 Integration with sktime

sktime is a python library for time series providing unified interfaces for
various tasks, such as classification, regression, or forecasting. For all of
these tasks, bootrapping can improve the results. Thus, in the following, we
first describe how tsbootstrap is integrated with sktime. We then show the
combined usage of sktime and tsbootstrap with a probabilistic forecasting
example.

####  5.2.1 sktime’s TSBootstrapAdapter

sktime’s TSBootstrapAdapter follows the unified API design as the other
transformers implemented in sktime. In other words, these methods provide fit,
transform, get_params, set_params methods, and __init__ as the constructor.
Below we briefly present these methods.

__init__

is the constructor. It takes as parameters, the tsbootstrap object that should
be used within sktime.

fit

fits the transformer. It is empty for this adapter.

transform

creates the bootstraps. It takes as arguments, the time series that should be
bootstrapped as $X$.

get_params

returns the current parameters of the adapter.

set_params

allows to set the different parameters.

This interface enables the usage of all tsbootstrap methods within sktime
since it requires that an arbitrary tsbootstrap object is passed to the
adapter and that all of these classes have the same interface. This interface
consistency guarantees that the classes are exchangeable.

####  5.2.2 Usage: Probabilistic Forecasting

Point forecasts comprises one value per time step. Thus, such forecasts do not
quantify uncertainty. To quantify the uncertainty in forecasting, often
probabilistic forecasts are generated. However, many forecasting algorithms do
not support probabilistic forecasting directly. Thus, in sktime, there exist
different methods to provide uncertainty quantification for point forecasts.
One of these approaches that smoothly interacts with bootstrapping is the
BaggingForecaster. The BaggingForecaster fits one model per bootstrap. During
inference, each fitted model makes a forecast. Finally, out of this set of
forecasts, the uncertainty can be derived.

Example Code 2 shows an exemplary usage tsbootstrap and sktime for
forecasting. First, the required modules are imported. Afterwards, the time
series is loaded and split into train and test set. Third, the forecaster is
initialised. In this example we use a BaggingForecaster. This forecaster takes
as arguments a bootstrapper and a base forecaster. For the bootstrapper, we
use BlockResidualBootstrap from tsbootstrap, which is wrapped using
TSBootstrapAdapter, and as base forecaster, we use a SARIMA model. Fourth, we
fit the forecaster using the training data. Finally, in the last step, we
forecast the time series and also forecast intervals.

# Example code snippet demonstrating the use of tsbootstrap together with #
sktime from sktime.forecasting.arima import ARIMA from
sktime.transformations.bootstrap import TSBootstrapAdapter from
sktime.forecasting.compose import BaggingForecaster from
sktime.transformations.series.detrend import Detrender, Deseasonalizer from
tsbootstrap import MovingBlockBootstrap from sktime.utils.plotting import
plot_series

from sktime.datasets import load_airline from
sktime.forecasting.model_selection import temporal_train_test_split

import matplotlib.pyplot as plt

# Load time series data with sktime or custom data with pandas data =
load_airline() train, test = temporal_train_test_split(data, test_size=12)

# Initialize the bootstrapper with custom settings adapter =
TSBootstrapAdapter( BlockResidualBootstrap( MovingBlockBootstrap(20) ) )

# Create a BaggingForecaster using the bootstrap adapter and a SARIMA model
forecaster = Deseasonalizer(sp=12, model=”multiplicative”) * Detrender() *
BaggingForecaster(adapter, ARIMA(order=(1, 1, 0), seasonal_order=(0, 1, 0,
12)),)

# Fit the forecaster forecaster.fit(train)

# Perform predictions pred = forecaster.predict(range(1, len(test) + 1))
pred_intervals = forecaster.predict_interval(range(1, len(test) + 1))

# Plot plt.close(”all”) fig, ax = plot_series(data, pred,
pred_interval=pred_intervals, labels=[”y”, ”y_pred”]) ax.set_xlabel(”Time”)
plt.show()

Example Code 2:  Exemplary usage of tsbootstrap together with sktime to
perform a forecast on the airline dataset.

The resulting forecast together with the intervals are shown in Figure 4.

Figure 4: Exemplary usage of tsbootstrap and sktime; output from Example Code
2.

##  6 Conclusion

We have demonstrated that tsbootstrap advances the field of time series
analysis by offering a robust and user-friendly python package designed for
sophisticated bootstrapping techniques. Its comprehensive suite of methods not
only addresses the specific needs of time series data but also enhances their
analysis with a focus on usability and integration. This tool can play a
crucial role in fields where accurate uncertainty quantification is essential,
such as finance, meteorology, and epidemiology. Through tsbootstrap,
researchers and practitioners can better predict and manage risks, ultimately
leading to more informed decision-making.

We invite the community to participate in the ongoing development of
tsbootstrap. We are particularly keen on collaborations that explore new
domains of application or enhance our understanding of time series analysis
challenges. By joining our efforts, you can help shape the future of
bootstrapping techniques in time series analysis. The continuous development
and enhancement of tsbootstrap, as highlighted in its GitHub repository, and
its close integration with sktime, ensures its relevance and utility in the
evolving landscape of time series analysis. We are excited to see how
tsbootstrap will continue to evolve and how you, the community, will be a part
of this journey.

##  7 Future Work

We have a clear yet evolving roadmap for future enhancements, including the
integration of more sophisticated bootstrapping methods and improved
performance optimizations. These are delineated in Issue #144 in the
tsbootstrap repository.

  1. 1.

Performance and Scaling: Improve handling of large datasets, by integrated
with established distributed backends such as Dask, Spark, and Ray. Explore
ways to optimize and profile the code for improvements in both time- and
space-complexity.

  2. 2.

Tuning and Automation: Implement adaptive algorithms that dynamically adjust
block sizes based on autocorrelation, and integrate fractional block lengths
for precise resampling (Bühlmann and Künsch, 1999). Additionally, implement
automated feedback mechanisms to refine the resampling process by analyzing
real-time dataset characteristics and iteratively adjusting to minimize
discrepancies between original and bootstrapped data.

  3. 3.

Real-time and Streaming Data: Process data in real-time by incorporating
event-driven programming or reactive frameworks that handle data streams
efficiently.

  4. 4.

Enhanced integration with sktime: Develop standardized evaluation tools that
utilize sktime’s metrics for detailed performance comparisons between original
and bootstrapped time series data. This includes establishing a repository of
shared datasets for benchmarking across both libraries, and creating extensive
documentation and examples that demonstrate the combined use of tsbootstrap
and sktime. Additionally, introduce functionalities to generate distribution
or sampler-like objects for advanced probabilistic forecasting.

  5. 5.

API Extension: Support pandas DataFrames, and expand API capabilities to
handle panel and hierarchical time series. At the same time, improve handling
of exogenous data within bootstraps for complex models, and refine model state
management to accommodate both fittable and pretrained models.

## References

  * Bühlmann (1997) Peter Bühlmann.  Sieve bootstrap for time series.  _Bernoulli_ , pages 123–148, 1997.
  * Bühlmann and Künsch (1999) Peter Bühlmann and Hans R Künsch.  Block length selection in the bootstrap for time series.  _Computational Statistics & Data Analysis_, 31(3):295–310, 1999.
  * Chernick (2011) Michael R Chernick.  _Bootstrap methods: A guide for practitioners and researchers_.  John Wiley & Sons, 2011.
  * Efron (1992) Bradley Efron.  Bootstrap methods: another look at the jackknife.  In _Breakthroughs in statistics: Methodology and distribution_ , pages 569–593. Springer, 1992.
  * Gilda (2024) Sankalp Gilda.  tsbootstrap, March 2024.  URL https://doi.org/10.5281/zenodo.10866090.
  * Grisel et al. (2024) Olivier Grisel, Andreas Mueller, Lars, Alexandre Gramfort, Gilles Louppe, Thomas J. Fan, Peter Prettenhofer, Mathieu Blondel, Vlad Niculae, Joel Nothman, Guillaume Lemaitre, Arnaud Joly, Loïc Estève, Jérémie du Boisberranger, Jake Vanderplas, manoj kumar, Hanmin Qin, Nicolas Hug, Nelle Varoquaux, Robert Layton, Adrin Jalali, Rajagopalan (Venkat) Raghav, Johannes Schönberger, Julien Jerphanion, Roman Yurchak, Lucy Liu, Christian Lorentzen, Tom Dupré la Tour, Wei Li, and Chiara Marmo.  scikit-learn/scikit-learn: Scikit-learn 1.4.2, April 2024.  URL https://doi.org/10.5281/zenodo.10951361.
  * Härdle et al. (2003) Wolfgang Härdle, Joel Horowitz, and Jens-Peter Kreiss.  Bootstrap methods for time series.  _International Statistical Review_ , 71(2):435–459, 2003.
  * Imon and Ali (2005) AHM Imon and M Masoom Ali.  Bootstrapping regression residuals.  _Journal of the Korean Data and Information Science Society_ , 16(3):665–682, 2005.
  * Kreiss and Lahiri (2012) Jens-Peter Kreiss and Soumendra Nath Lahiri.  Bootstrap methods for time series.  In _Handbook of statistics_ , volume 30, pages 3–26. Elsevier, 2012.
  * Kunsch (1989) Hans R Kunsch.  The jackknife and the bootstrap for general stationary observations.  _The annals of Statistics_ , pages 1217–1241, 1989.
  * Lahiri (2013) Soumendra Nath Lahiri.  _Resampling methods for dependent data_.  Springer Science & Business Media, 2013.
  * pandas development team (2024) The pandas development team.  pandas-dev/pandas: Pandas, April 2024.  URL https://doi.org/10.5281/zenodo.10957263.
  * Paparoditis and Politis (2001) Efstathios Paparoditis and Dimitris N Politis.  Tapered block bootstrap.  _Biometrika_ , 88(4):1105–1119, 2001.
  * Pedregosa et al. (2011) Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  _the Journal of machine Learning research_ , 12:2825–2830, 2011.
  * Politis and Romano (1994) Dimitris N Politis and Joseph P Romano.  The stationary bootstrap.  _Journal of the American Statistical association_ , 89(428):1303–1313, 1994.
  * RNKuhns et al. (2024) RNKuhns, Franz Király, Tomas Pereira de Vasconcelos, Anirban Ray, LGTM Migrator, MBristle, and Sagar Mishra.  sktime/skbase: v0.7.7, April 2024.  URL https://doi.org/10.5281/zenodo.10984314.
  * sktime developers (2024) sktime developers.  sktime/sktime: v0.28.1, April 2024.  URL https://doi.org/10.5281/zenodo.7117735.


# A general method for estimating the prevalence of Influenza-Like-Symptoms
with Wikipedia data

Giovanni De Toni  DISI, University of Trento<EMAIL_ADDRESS>,
Cristian Consonni  Eurecat - Centre Tecnòlogic of Catalunya
<EMAIL_ADDRESS>and  Alberto Montresor  DISI, University of
Trento<EMAIL_ADDRESS>
††copyright: none

## Abstract

Influenza is an acute respiratory seasonal disease that affects millions of
people worldwide and causes thousands of deaths in Europe alone. Being able to
estimate in a fast and reliable way the impact of an illness on a given
country is essential to plan and organize effective countermeasures, which is
now possible by leveraging unconventional data sources like web searches and
visits. In this study, we show the feasibility of exploiting information about
Wikipedia’s page views of a selected group of articles and machine learning
models to obtain accurate estimates of influenza-like illnesses incidence in
four European countries: Italy, Germany, Belgium, and the Netherlands. We
propose a novel language-agnostic method, based on two algorithms,
Personalized PageRank and CycleRank, to automatically select the most relevant
Wikipedia pages to be monitored without the need for expert supervision. We
then show how our model is able to reach state-of-the-art results by comparing
it with previous solutions.

##  1\. Introduction

Influenza is a widespread acute respiratory infection that occurs with
seasonal frequency, mainly during the autumn-winter period. The European
Centre for Disease and Control (ECDC) estimates that influenza causes up to 50
millions cases and up to $70,000$ deaths annually in Europe (ECDC2018, ;
thompson2003mortality, ). Globally, the number of deaths caused by this
infection ranges from $290,000$ to $650,000$ (WHO2019, ). The most affected
categories are children in the pediatric age range and seniors. Severe cases
are recorded in people above 65 years old and with previous health conditions,
for example, respiratory or cardiovascular illnesses (thompson2009estimates, ;
nair2011global, ).

The ECDC monitors the influenza virus activities and provides weekly bulletins
that aggregate the open data coming from all the European countries. The local
data are gathered from the national networks of volunteer physicians
(ECDC2017, ). The Centre supplies only data about the current situation, which
means that they do not provide any valuable information about which will be
the future impact of the influenza illness in a European country, namely the
number of future infected people. Moreover, ECDC data are available with a
delay ranging between one and two weeks, which prevents taking preventive
measures to mitigate the impact of influenza. In recent years, thanks to
machine learning techniques, new methods have been developed to estimate the
activity level of diseases by harnessing the power of the internet and big
data. Researchers tried to use unconventional sources of data to make
predictions, for example, tweets from Twitter (signorini2011use, ;
santos2014analysing, ; kim2013use, ), Google’s search keywords
(yang2015accurate, ) or social media (woo2016estimating, ).

Wikipedia111https://www.wikipedia.org is a multilingual encyclopedia, written
collaboratively by volunteers online. As of this writing, the English version
alone has more than 6M articles and 50M pages and it is edited on average by
more than 65k active users every month (wikipedia2016statistics, ).

Research has shown that Wikipedia is a first-stop source for information of
all kinds, including information about science (yasseri2014most, ;
spoerri2007popular, ). This is of particular significance in fields such as
medicine, where it has been shown that Wikipedia is a prominent result for
over $70\%$ of search results for health-related keywords (laurent2009seeking,
). Wikipedia is also the most viewed medical resources globally
(heilman2015wikipedia, ), and it is used as an information source by $50\%$ to
$70\%$ of practicing physicians (heilman2011wikipedia, ).

In this study, we investigate how to exploit Wikipedia as a source of data for
doing rapid now-casting of the influenza activity levels in several European
countries. McIver and Brownstein (mciver2014wikipedia, ) developed a method
based on Wikipedia’s pageviews, which behaves quite well when predicting
influenza levels in the United States provided by the Centers for Disease
Control and Prevention (CDC). They use pageviews of a predefined set of pages
from the English language edition of Wikipedia chosen by a panel of experts.
Generous et al. (generous2014global, ) also developed linear machine learning
models and demonstrated the ability to forecast incidence levels for 14
location-disease combinations.

We extend upon the state of the art by proposing a novel method that can be
applied automatically to other languages and countries. We focused our
attention on predicting influenza levels in Italy, Germany, Belgium and the
Netherlands. In extending this work, we have developed a method that is more
flexible with respect of the characteristics of the edition of Wikipedia and
that can be applied to any language thus enabling us to build and deploy
machine learning models for Influenza-like Illnesses (ILI) estimation without
the need of expert supervision.

The paper is structured as follows: in Section 2 we review the relevant
datasets that we have used in our study. In Section 3 and 4 we discuss our
models and how we selected the relevant features to use for predicting
influenza activity levels from Wikipedia data. We present the results of our
study in Section 5 and 6, and we discuss them in Section 7. Finally, we
conclude the paper in Section 8.

##  2\. Data sources

##### Wikipedia page views

Data about Wikipedia page views are freely available to download. For our
study, we have used two datasets: the pagecounts dataset222Pagecounts-raw on
wikitech.wikimedia.org, https://w.wiki/gnq, visited on 2020-10-14 (December
2007-May 2016, for a total of 461 weeks) and the pageviews dataset (October
2016 to April 2019, for a total of 134 weeks). Both of them count the number
of page hits to Wikipedia articles. The pagecounts dataset contains non-
filtered page views, including automatically bot-generated ones; moreover, it
does not contain data of page hits from mobile devices. The pageviews dataset
has been recently developed, superseding pagecounts in October 2016, and
includes only human traffic from both desktop and mobile devices. By
considering both these datasets, we can examine a longer time range while
estimating the impact of different models on our results.

From these datasets, we extracted the page views of three different versions
of Wikipedia: Italian, German, and Dutch. We selected these languages because
the majority of page views for each of them come almost exclusively from one
single area.333Report: Breakdown of page views by language on
stats.wikimedia.org, https://w.wiki/gnr visited on 2020-10-14 For instance,
just $41.5\%$ of the page views of the English Wikipedia come from US and UK,
while all the rest comes from all over the world. On the contrary, $87.0\%$ of
the page views of the Italian version of Wikipedia come from Italy, $77.0\%$
of the page views of the German version of Wikipedia come from Germany, and
almost $89\%$ of the Dutch Wikipedia’s page views come from the Netherlands
and Belgium ($20.3\%$ coming from Belgium and $69.1\%$ coming from the
Netherlands alone). This enables us to build datasets which contain less noise
and which are only affected by the internet behaviour of the population of the
country or area we are examining.

Each of row is composed of four columns, namely: the project’s language, the
page title, the number of requests and the size in byte of the page retrieved.
Data are available with hourly granularity; to use them in our project, we
filtered them to select only specific entries and then we computed an
aggregated weekly value of the page views for each of the entries selected.

We generated weekly aggregates for page views for the period ranging from the
50th week of 2007 to the 15th week of 2019, for a total of 591 weeks.
pageviews data were used to cover the period from September 2016 to April 2019
period, while for the previous period pagecounts data were used.

We processed weekly page views by standardizing them, i. e. removing the mean
and then dividing by the standard deviation. In case of missing data, we used
a padding technique to propagate the last valid value to fill the empty
entries. Since Wikipedia articles are created collaboratively by the community
of editors, their creation date may vary: we assigned the value of zero to the
page view count to the period when they did not exist. We included the week
number as a feature of our model. In order to do so, we one-hot encoded the
week value into binary vectors $w_{i}\in\\{0,1\\}^{52}$ where the bit
corresponding to the given week is set to $1$.

We then built two different datasets named PV and PC+PV. PV contains only the
pageviews data, thus covering only four influenza seasons (from 2015-2016 to
2018-2019). PC+PV was obtained by merging pagecounts and pageviews data, thus
covering all the influenza seasons (from 2007-2008 to 2018-2019).

##### Influenza activity levels

Data about the incidence of influenza in the population are available through
different online official websites and repositories. We collected influenza
data for these countries: Italy, Germany, Austria and Netherlands. In the
following we describe briefly the collection procedure we followed for each
country.

Italian data were downloaded from the InfluNet system444InfluNet on
epicentro.iss.it, https://www.epicentro.iss.it/influenza/influnet, visited on
2020-10-15, which is the Italian flu surveillance program managed by the
Istituto Superiore di Sanità; German data were obtained from the Robert Koch
Institute555Robert Koch Institute, on www.rki.de,
https://www.rki.de/DE/Home/homepage_node.html, visited on 2020-10-20; Dutch
and Belgian data were take from the GISRS (Global Influenza Surveillance and
Response System) online tool which is available from the WHO website666WHO
Influenza Surveillance on www.who.it/influenza/en,
https://www.who.int/influenza/en/, visited on 2020-10-20, since at the best of
our knowledge, the two countries do not offer open data from national health
institute websites.

For Germany and Italy, we cover 12 influenza seasons (from 2007-2008 to
2018-2019); for the Netherlands and Belgium, we have data that cover 9
influenza seasons (from 2010-2011 to 2018-2019). Each of the datasets records
the influenza incidence aggregated weekly over 100 000 patients. For each
influenza seasons, each of the dataset entries represents the influenza level
on a specific week.

For each influenza seasons, data cover 26 weeks, from the 42nd week of the
previous year to the 15th week of the following year; for instance, the
2009-2010 influenza season comprises data ranging from the 42nd week of 2009
to the 15th week of 2010. The complete dataset has the following size: 312
weeks of data for Germany and Italy, 234 for the Netherlands and Belgium.

Table 1. Selected Wikipedia categories. We first selected the Italian categories and then we chose the corresponding German and Dutch translation. The English categories are reported for translation purposes. English | Italian | German | Dutch
---|---|---|---
Viral diseases | Malattie infettive virali | Virusziekte | Virale Infektionskrankheit
Infectious diseases | Malattie infettive | Infectieziekte | Infektionskrankheit
Epidemics | Epidemie | Epidemie | Epidemie
Viruses | Virus | Virus | Viren, Viroide und Prionen
Vaccines | Vaccini | Vaccin | Impfstoff
Medical signs | Segni clinici | Symptoom | Krankheitssymptom

##  3\. Feature Selection

The selection of the pages to monitor was based on the following assumption:
during the influenza seasons, we should see an increment in the page view
number of some specific Wikipedia’s pages, which should be related to flu
topics (e.g., _“Influenza”_ , _“Headache”_ , _“Flu vaccine”_ , etc.). That
occurs because most people tend to search online the symptoms of their disease
when they are sick, trying to get more information about their illness.

We use three different methodologies in order to avoid to make unnecessary
assumptions about which pages to choose.

###  3.1. Exploiting Wikipedia’s Categories

The first method consists of determining a series of Wikipedia’s categories
related to the medical sector, from which to extract the single entries. We
started by choosing the categories from the Italian version of Wikipedia.
After doing that, we replicated the procedure for all the other languages
(German, Dutch) by selecting the corresponding translation of the given
category in the Italian Wikipedia. We ended up having a dedicated set of
Wikipedia’s pages for each of the different Wikipedia’s versions. In the end,
we produced three different lists of Wikipedia’s pages which were used to
filter the Wikipedia’s page view dumps.

The categories selected by this method can be seen in Table 1.

###  3.2. Graph-based methods

The strategies that we propose in the following section can choose
automatically the set of relevant pages to monitor. In this way, we could keep
our page list always up-to-date with minimum effort, yet retaining its
effectiveness.

Wikipedia can be represented as a graph in which each article is a node and
links between them are directed edges. Wikipedia editors write articles and
insert link to other pages which are contextually relevant; links are also
used by readers to navigate between pages. When choosing which pages are
relevant we should consider the entire graph structure; for example, we may
consider the pages linked by the _“Influenza”_ article to be relevant, but
also the ones that link to it. Previous work has showed that Wikipedia links
provide a way to explore the context of an article (west2015, ) and we can
exploit this knowledge to extract the best pages to monitor. Ideally, we would
like to be able to give to each Wikipedia page a score that measures its
relative importance to a given “source” page, in this case the _“Influenza”_
page for each language.

There are several algorithms that can be used to given scores to graph nodes
based on specific metrics. We adopt here two approaches: the well-known
Personalized PageRank (page1999pagerank, ) and a more recent proposal called
CycleRank (consonni2020, ).

#### Personalized PageRank

The Personalized PageRank algorithm (page1999pagerank, ) (PPageRank for short)
ranks web pages in order of importance with respect to a set of source pages.
PPageRank models the relevance of nodes around the selected nodes as the
probability of reaching each of them, when following random walks starting
from one of the sources nodes.

This algorithm has been applied successfully to a broad range of graph
structures in order to find relevant items (doi:10.1137/140976649, ;
spark2019, ) and we include it as a baseline in our analysis. Its performance,
however, is hindered by pages with high in-degree that function as hubs and
obtain high scores regardless of the starting point.

#### CycleRank

The CycleRank algorithm (consonni2020, ) is a novel approach based on cycles,
i.e. circular walks that starts and end in a given node, exploiting both
incoming and outgoing links to reduce the importance of in-degree hubs. The
idea behind this method is that by using circular walks, it is possible to
identify pages that are more pertaining to the context of the topic of
interest. Cycles guarantee that we rank higher pages that are, at least
indirectly, both linked from and linking to the reference article, thus
avoiding the pitfalls of PPageRank where random walks may easily lead far away
from the original topic, into pages that are not relevant at all.

Table 2. Percentages of features shared by each feature set given the method used to extract them. Since the various datasets have different size, the percentages are not symmetrical. Each row shows the result for one of the methods. For instance, the CycleRank row shows the fraction of features it shares with the other feature sets. | CycleRank | PPageRank | Categories
---|---|---|---
| Italian | German | Dutch | Italian | German | Dutch | Italian | German | Dutch
CycleRank | 100% | 100% | 100% | 53,75% | 49,36% | 83,72% | 18,75% | 10,97% | 41,86%
PPageRank | 35,10% | 42,23% | 26,27% | 100% | 100% | 100% | 10,20% | 8,30% | 15,32%
Categories | 6,38% | 6,82% | 9,42% | 5,31% | 6,04% | 10,99% | 100% | 100% | 100%

The CycleRank score can be seen as the time spent on a given node when
following random loops from the reference node, assuming a fixed overall time
for each loop, equally split among all the nodes encountered. The CycleRank
algorithm is similar to PPageRank, because it can be explained as the
probability of ending up in a node when following a random path that starts
and ends in the source nodes.

###  3.3. Discussion

We then analyze the set of features extracted by using these different
methodologies. For the Dutch, German and Italian Wikipedia versions,
respectively, the Categories feature sets contain 382, 381 and 470 Wikipedia’s
pages; the PPageRank feature sets contain 274, 277 and 245 pages; and the
CycleRank feature sets contain 86, 237 and 160 pages.

The feature sets built from the categories contain more pages than the other
two, this is expected since we do not impose any constraint on them. The
CycleRank and PPageRank feature sets are smaller since we enforce that all
those pages must be connected to the _“Influenza”_ page; CycleRank sets are
the smallest ones because they require connections in both directions. In
general, the Categories feature set encompass a broad spectrum of pages
related to medical topics, which may or may not have any relation with the
_“Influenza”_ concept.

Table 2 shows how much the features sets overlaps across the three methods
under consideration. CycleRank and PPageRank feature lists share many common
pages. For instance, $42\%$ of the pages selected by CycleRank are also
selected by the PPageRank, while Categories feature sets are more diverse,
with many pages appearing only with this method. However, for the Dutch
Wikipedia, CycleRank shares $42\%$ of its features with the Categories Dutch
feature set.

At a first glance, Categories method outlined above is difficult and prone to
errors. Choosing the correct pages to monitor is a time-consuming task which
requires some expertise about the medical and epidemiological domain.
Moreover, creating and maintaining this list of pages is a process which needs
to be done manually and frequently. Wikipedia’s structure can undergo rapid
changes and each of its pages can be deleted, moved or renamed, thus forcing
us to keep our list of features constantly updated . On the other hand,
CycleRank and PPageRank do not require human intervention and they are surely
faster than the manual Categories.

##  4\. Models

We trained a simple linear model regularized with the least absolute shrinkage
and selection (LASSO) method (tibshirani1996regression, ). The LASSO
regularization diminishes the possible over-fitting and automatically reduces
the number of features used by the models. The models were trained by
stochastic gradient descent (SGD) and by trying to minimize the Mean Squared
Error (MSE). During training, we used cross-validation to estimate the perfect
balance for the LASSO regularization term. We tried to use the simplest
technique possible, linear regression, because we assumed that the correlation
between the influenza incidence and the page views of Wikipedia was linear.
Moreover, we can use the linear model’s weights to visualize the relative
importance of each page towards the final prediction.

We tested and evaluated three different linear models: Categories, PPageRank
and CycleRank, each reflecting the set of features used. Each model was
trained using the PV and PC+PV datasets mentioned above.

The models were trained on all the influenza seasons by leaving out the one we
wanted to estimate. Basically, we trained a different model for each influenza
season. We computed an estimate of the final performances as an average over
all the trained models. Apart from the last year, every model has been
evaluated on a given influenza season (e.g., 2015-2016), but the model has
been trained by using also future influenza season data (e.g., 2016-2017,
2017-2018, etc.). We argue that this is acceptable since we assume the various
influenza period to be conditionally independent. Despite many factors could
make this assumption less precise, we claim that this experimental setting is
still appropriate to give us an idea about the performance of the final
complete models, which will be used to estimate ILI case on new data.

The code777https://github.com/fluTN/influenza and the dataset
(detoni2020fluTN, ) used for the project are released as an open-source
project to ensure correct peer-review and to assure reproducibility of our
results.

##  5\. Results

We compared the performance of the CycleRank and PPageRank models against the
one trained with the user-defined categories (Categories). We trained the
models to predict the last four influenza seasons (2015-2016, 2016-2017,
2017-2018 and 2018-2019). This analysis was done for all the chosen European
countries: Italy (IT), Germany (DE), Belgium (BE) and the Netherlands (NL). We
measured the Pearson Correlation Coefficient (PCC) to determine which of the
models held the best performances by taking the mean of the various PCC
obtained for each of the four seasons.

We also evaluated the accuracy in estimating the week in which the influenza
peak occurs (highest incidence value over the entire season). We considered a
peak estimation to be acceptable if the model was able to estimate the peak
with a deviation of $\pm 2$ weeks with respect to the ground truth. The models
were trained by using the PV and PC+PV datasets. Namely, we wanted to analyze
the effects of providing additional data during training, even if it was
recorded with a different methodology. In this work we also tried to evaluate
the goodness of the most important features selected by the various models. We
tried to understand which Wikipedia’s pages are the most valuable when
estimating the influenza incidence over a population. Moreover, we wanted also
to assess the effectiveness of our automatic method with respect to the
traditional, manual one (Categories).

###  5.1. Estimation Accuracy

Table 3 shows the complete results for each model and each country. All models
present a high PCC with the influenza incidence data. Generally, CycleRank and
PPageRank yield better results or provide similar performances to the
Categories model. For instance, the latter shows lower performances on Germany
and the Netherlands (both with the PV and PC+PV datasets). The PPageRank model
scores the highest correlation coefficient on 3 out of 4 of the examined
countries (IT, DE, BE) with the best model for NL being CycleRank. In general,
PPageRank and CycleRank show a lower variance in PCC than the Categories
model.

Examining the effect of using different datasets to train the various models,
in 8 over 12 experiments, the PC+PV dataset improves the PCC of about $10\%$
on average. This effect is visible for CycleRank and Categories, while
PPageRank seems to not obtain the same gain when employing a larger page view
dataset. In the remaining four experiments, we noticed that using a larger
training dataset produced an average accuracy loss of around $5\%$.

Figures 1 and 2 show the predicted influenza trends for each of the given
countries when using PV and PC+PV, respectively. Sometimes the models produce
visible “spikes”. Those outliers are related to a sudden page views increment
(or decrement) of one or more of the feature pages. These sudden changes in
page views counts could be caused by many factors. For instance, unexpected
news coverage of a certain topic may artificially increase the total views of
a certain page, thus leading to overestimation of the results.

Ideally, we want our model to be robust against these kinds of situation and
we prefer models which generate less “spikes”. This objective reflects also on
the set of pages we are using to monitor ILI incidence. We would like our
model to select highly informative pages, which however do not present sudden
changes in the page views count. We can see from the graphs how the Categories
model (green line) generates many of these unwanted values. On the contrary,
CycleRank and PPageRank are more resilient and produces better predictions.
Again, by using the PC+PV dataset as training set we can increase the
robustness of all the models thus making them less susceptible to sudden page
views changes.

Table 3. Mean Pearson Correlation Coefficient. The bold values indicate the best results over all combinations of models and datasets. The green color indicates the best results by either using the PV or the PC+PV dataset, by keeping fixed the country and the model (no color means no changes). | CycleRank | PPageRank | Categories
---|---|---|---
Country | PV | PC+PV | PV | PC+PV | PV | PC+PV
IT |  [HTML]FFCCC90.898 |  [HTML]9AFF990.906 |  [HTML]FFCCC90.879 |  [HTML]9AFF990.923 |  [HTML]FFCCC90.790 |  [HTML]9AFF990.913
DE |  [HTML]FFCCC90.755 |  [HTML]9AFF990.865 |  [HTML]9AFF990.914 |  [HTML]FFCCC90.874 |  [HTML]FFCCC90.482 |  [HTML]9AFF990.548
NL |  [HTML]9AFF990.757 |  [HTML]FFCCC90.681 |  [HTML]9AFF990.746 |  [HTML]FFCCC90.703 |  [HTML]FFCCC90.423 |  [HTML]9AFF990.657
BE |  [HTML]FFCCC90.487 |  [HTML]9AFF990.793 |  [HTML]FFCCC90.625 |  [HTML]9AFF990.809 |  [HTML]9AFF990.804 |  [HTML]FFCCC90.733
Table 4. Number of influenza peaks predicted correctly ($\\#peaks=4$). The bold values indicate the best results over all combinations of models and datasets. The first value in each cell indicates how many peaks were predicted correctly, the second value indicates how many peeks where predicted in the $\pm 2$ weeks range. The green color indicates the best results by either using the PV or the PC+PV dataset, by keeping fixed the country and the model (no color means no changes). | CycleRank | PPageRank | Categories
---|---|---|---
Country | PV | PC+PV | PV | PC+PV | PV | PC+PV
IT | 0 (4) | 0 (4) |  [HTML]FFCCC91 (3) |  [HTML]9AFF991 (4) |  [HTML]FFCCC90 (3) |  [HTML]9AFF991 (4)
DE |  [HTML]FFCCC90 (3) |  [HTML]9AFF992 (4) |  [HTML]FFCCC91 (4) |  [HTML]9AFF992 (4) | 1 (3) | 1 (3)
NL |  [HTML]FFFFFF0 (3) |  [HTML]FFFFFF0 (3) |  [HTML]9AFF990 (3) |  [HTML]FFCCC90 (2) |  [HTML]FFCCC90 (2) |  [HTML]9AFF990 (3)
BE |  [HTML]FFCCC90 (2) |  [HTML]9AFF990 (3) |  [HTML]FFCCC90 (2) |  [HTML]9AFF990 (3) |  [HTML]FFCCC90 (4) |  [HTML]9AFF991 (3)

###  5.2. Peak Accuracy

The peak prediction performance can be seen in Table 4. The models are able to
estimate the highest influenza peak within a 2-weeks range with respect to the
ground-truth value. Unfortunately, they do not reach a level of accuracy high
enough to estimate correctly the peak. For instance, on Dutch and Belgian
data, the models do not estimate correctly any peak week, but they still
position it in the two-weeks range. In this case, by examining the models
trained with the extended dataset (PC+PV), in 8 models over 12 we can see a
small increment in the correctly predicted peaks.

(a)

(b)

(c)

(d)

Figure 1. CycleRank, PPageRank and Categories models predictions using the PV
dataset on the Italian, German, Belgian and Dutch influenza incidence. The
dashed lines delimits the various influenza seasons.

(a)

(b)

(c)

(d)

Figure 2. CycleRank, PPageRank and Categories models predictions using the
PC+PV dataset on the Italian, German, Belgian and Dutch influenza incidence.
The dashed lines delimits the various influenza seasons.

##  6\. Feature Analysis

In this section we analyze which features are selected by each model with the
aim of characterizing how the chosen features affect the final estimate.
Namely, which are the best predictors for estimating ILI incidence with
Wikipedia? Are the selected predictors (Wikipedia’s pages) related in some way
to the _“Influenza”_ topic? Is there a difference between using hand-picked
features or automatically discovered ones?

Table 5. Analysis of the number of features selected by the models for each
influenza season. We record the minimum, maximum and mean number of features
for all the countries and for the different dataset used.

Country | CycleRank | PPageRank | Categories
---|---|---|---
PV | PC+PV | PV | PC+PV | PV | PC+PV
Min | Max | Mean | Min | Max | Mean | Min | Max | Mean | Min | Max | Mean | Min | Max | Mean | Min | Max | Mean
IT | 30 | 65 | 48.75 | 37 | 135 | 64.25 | 29 | 71 | 49.25 | 44 | 165 | 82.25 | 57 | 81 | 74.00 | 17 | 93 | 58.25
DE | 27 | 71 | 55.25 | 84 | 96 | 89.25 | 22 | 51 | 37.50 | 69 | 123 | 101.00 | 34 | 75 | 61.00 | 75 | 89 | 81.75
NL | 7 | 30 | 17.75 | 30 | 40 | 34.75 | 17 | 67 | 34.25 | 19 | 82 | 52.00 | 20 | 67 | 51.75 | 47 | 82 | 66.00
BE | 8 | 39 | 23.75 | 29 | 34 | 31.75 | 20 | 67 | 42.50 | 37 | 56 | 48.50 | 16 | 70 | 38.50 | 7 | 12 | 9.25

###  6.1. Top-5 Wikipedia’s Predictors

First, we focus on the estimated relative importance of each predictor. We
wanted to find which pages are the most useful to estimate the ILI incidence
of a given week. We focused on the weights assigned to the various Wikipedia’s
pages by the models. We directed our attention to those features which
received a positive weight. A positive weight means that there is a positive
correlation between the predictor and the influenza incidence. To detect the
most successful subset of features for the final model, we ranked the features
based on their mean weight over the different models. If a given Wikipedia’s
page is selected multiple times by separate models for different influenza
seasons, it means that we found a feature which can be used to predict the ILI
incidence over many seasons.

We collected the top-5 features found by each model type (CycleRank, PPageRank
and Categories) and for each dataset, PV and PC+PV. We also computed the
shortest path distance $D_{I}$ between the _“Influenza”_ page to the found
predictor and the Pearson Correlation Coefficient (PCC) for each chosen page
against the influenza incidence. The complete results can be seen in Table 6d.
From the results, we can notice how the top-5 predictors of many models are
Wikipedia’s pages related to the influenza topic. If we perform a rapid
semantic analysis, some of them refer directly to the symptoms (e.g.,
_“Fever”_ , _“Bronchitis”_ , _“Chills”_), the pathogens (e.g., _“Influenza A
virus subtype H1N1”_) or to similar concepts related to ILI (e.g.,
_“Peramivir”_ , _“Swine Influenza”_). We can also identify some super-
predictors. Namely, some features which are chosen multiple times by several
models (e.g., _“Febbre”_ for the Italian models or _“Grippe”_ for both the
Belgian and Dutch models). These super-predictors tend to have a very high
PCC. Thus they are very valid in predicting the variation of the ILI
incidence. Ultimately, 9 out the best 12 models considered in Table 6d show to
have at least half of their top-5 predictors with a distance $D_{I}$ less or
equal than 2 from the _“Influenza”_ page. More interestingly, the Categories
models show the same behaviour, thus hinting that there may be some
relationship between the optimal predictors and their position in Wikipedia’s
graph.

###  6.2. Size of the feature subsets

Table 5 shows some statistics about how many features were selected by each
model during training. We record the minimum, maximum and mean value of the
total features chosen by all the models. Generally, we can see how the final
models always produce subsets of the original feature sets. Not all the
Wikipedia’s pages available are good predictors and selecting too many
features would cause loss of generalization. Models trained with the PC+PV
dataset show a higher mean than their PV counterparts. This can be explained
by the fact that by employing a larger dataset we may need more predictors to
accurately fit the training set. Ultimately, we can see how the CycleRank
models tend to choose a smaller feature subset for all the countries still
retaining the same accuracy. In conclusion, we argue that we need to monitor
just a few pages to be able to obtain good estimates of the ILI incidence.

###  6.3. Features shared between models

We were also interested in investigating internal predictor variations which
may occur when training models by changing the initial feature set. Figure 3
shows how many features are shared by the final trained models. First of all,
if we compare the models trained with both the available datasets, PV and
PC+PV, but with fixed feature set (e.g., using only CycleRank), we can notice
how they will both select identical features. The CycleRank models seem to be
less affected (more than $60\%$ shared features between PV and PC+PV), while
PPageRank and Categories show a higher degree of variation. This fact is
interesting because it tells us that providing a better dataset with more
influenza seasons, such as PC+PV, affects the final predictors and it can push
the models toward choosing different Wikipedia’s pages. These changes, in
turn, could lead to an improvement in the performances.

In all the other cases, we can notice how the various models overlap just in
some features. The models end up selecting almost non-intersecting sets of
Wikipedia’s voices. For instance, Categories shares only up to $13\%$ of the
final predictors with CycleRank and PPageRank, regardless of the training
dataset used. Despite this, we can retain comparable performances. Again, this
behaviour suggests that there are just a few Wikipedia’s pages worth to
monitor to estimate the ILI incidence.

Table 6. Top-5 Wikipedia’s pages selected by the models for the influenza
seasons from 2015 to 2019 for all the examined countries. Since each model
could have been trained with either the PV or PV+PC dataset, we report only
the models which performed better (see Table 3 for the best models). For each
model, we report the page name, the shortest-path distance $D_{I}$ between the
page and the corresponding _“Influenza”_ page and the Pearson Correlation
Coefficient (PCC) measured against the influenza incidence. We also report the
corresponding page in the Engligh Wikipedia between parentheses. The value NE
is used to specify when a page has no english equivalent. The value $D_{I}>3$
indicates that the page is more than three hops away from the _“Influenza”_
page.

Cyclerank | PPageRank | Categories
---|---|---
Page Name | PCC | $D_{I}$ | Page Name | PCC | $D_{I}$ | Page Name | PCC | $D_{I}$
Febbre (Fever) | 0.55 | 1 |  Febbre (Fever) | 0.55 | 1 |  Febbre (Fever) | 0.55 | 1
Oseltamivir (Oseltamivir) | 0.36 | 1 |  Influenzavirus A sottotipo H1N1 (Influenza A virus subtype H1N1) | 0.76 | ¿ 3 |  Influenzavirus A sottotipo H1N1 (Influenza A virus subtype H1N1) | 0.76 | ¿ 3
Bronchite (Bronchitis) | 0.39 | 1 |  Oseltamivir (Oseltamivir) | 0.36 | 1 |  Bronchite (Bronchitis) | 0.39 | 1
Tosse post-virale (Post-viral cough) | 0.78 | ¿ 3 |  Bronchite (Bronchitis) | 0.39 | 1 |  Influenza asiatica (Influenza A virus subtype H2N2#Asian flu) | 0.67 | ¿ 3
Vaccino influenza stagionale (Influenza vaccine) | 0.03 | ¿ 3 |  Polmonite (Pneumonia) | 0.45 | 1 |  Fago lambda (Lambda phage) | 0.56 | ¿ 3

(a) Italy

Cyclerank | PPageRank | Categories
---|---|---
Page Name | PCC | $D_{I}$ | Page Name | PCC | $D_{I}$ | Page Name | PCC | $D_{I}$
Influenza-Schnelltest (NE) | 0.94 | 1 |  Influenza-Schnelltest (NE) | 0.94 | 1 |  Erkältung (Common cold) | 0.59 | 1
Rifttalfieber (Rift Valley fever) | 0.61 | 2 |  Pathogen (Pathogenicity) | -0.13 | 2 |  Impfstoff (Vaccine) | 0.09 | 2
Schweineinfluenza (Swine influenza) | 0.60 | 2 |  Totenschein (Death certificate) | 0.24 | ¿ 3 |  Bradykardie (Bradycardia) | 0.66 | 2
Peramivir (Peramivir) | 0.25 | 2 |  Asiatische Grippe (1957–1958 influenza pandemic) | 0.79 | ¿ 3 |  Schüttelfrost (Chills) | 0.70 | 1
Impfung (Vaccination) | 0.17 | 2 |  Pferdeinfluenza (Equine influenza) | 0.31 | 2 |  Rhinosinusitis (NE) | 0.23 | ¿ 3

(b) Germany

Cyclerank | PPageRank | Categories
---|---|---
Page Name | PCC | $D_{I}$ | Page Name | PCC | $D_{I}$ | Page Name | PCC | $D_{I}$
Koorts (Fever) | 0.45 | 1 |  Koorts (Fever) | 0.45 | 1 |  Griep (Influenza) | 0.59 | 1
Griep (Influenza) | 0.59 | 1 |  Griep (Influenza) | 0.59 | 1 |  Cytokinestorm (Cytokine storm) | 0.35 | 2
Mexicaanse griep (2009 swine flu pandemic) | 0.56 | ¿ 3 |  Chordadieren (Chordate) | 0.33 | ¿ 3 |  Flavivirus (Flaviviridae) | 0.33 | ¿ 3
Azerbeidzjan (Azerbaijan) | 0.29 | 2 |  1950-1959 (1950s) | 0.18 | ¿ 3 |  Salkvaccin (NE) | 0.04 | ¿ 3
Virulentie (Virulence) | 0.03 | 2 |  Antigene drift (Antigenic drift) | 0.33 | ¿ 3 |  Lassakoorts (Lassa fever) | 0.39 | ¿ 3

(c) Belgium

Cyclerank | PPageRank | Categories
---|---|---
Page Name | PCC | $D_{I}$ | Page Name | PCC | $D_{I}$ | Page Name | PCC | $D_{I}$
Griep (Influenza) | 0.75 | 1 |  Griep (Influenza) | 0.75 | 1 |  Griep (Influenza) | 0.75 | 1
Hongkonggriep (Hong Kong flu) | 0.51 | 1 |  Hongkonggriep (Hong Kong flu) | 0.51 | 1 |  Koorts (Fever) | 0.57 | 1
Aziatische griep (1957–1958 influenza pandemic) | 0.61 | ¿ 3 |  Influenzavirus A (NE) | 0.38 | ¿ 3 |  Mexicaanse griep (2009 swine flu pandemic) | 0.43 | ¿ 3
Influenzavirus A (NE) | 0.38 | ¿ 3 |  Chirurgie (Surgery) | 0.25 | 2 |  Cytokinestorm (Cytokine storm) | 0.48 | 2
Amantadine (Amantadine) | -0.17 | 1 |  Immuniteit (NE) | 0.14 | ¿ 3 |  Respiratoir syncytieel virus (Respiratory syncytial virus) | 0.07 | ¿ 3

(d) Netherlands

(a)

(b)

(c)

(d)

Figure 3. Percentage of features shared by each model. We compare how many
features are shared by each model (CRank=CycleRank, PPRank=PPageRank,
Cat=Categories). We also report the percentage of feature shared by models
trained with the two different dataset, PV and PC+PV. We recorded all the
features selected by each model for each influenza season (a feature was
included if and only its weight in at least one model was different than
zero). Then, we computed the intersection between each set of features for
each combination model/dataset.

##  7\. Discussion

The methods outlined in the previous sections complement more traditional
techniques, especially if these conventional approaches present a lag between
the collection of the data and the publications of the results. For instance,
some illnesses require to conduct some laboratory tests such to have a
positive/negative response. In such scenarios, it could be beneficial to
employ these Wikipedia-based solutions to gain rapidly additional (albeit
maybe less precise) information. More importantly, we have shown the
feasibility of using automatic methods to extract relevant predictors from
Wikipedia itself by exploiting Wikipedia’s links structure. Wikipedia’s
structure is highly flexible, and its pages could be removed or changed very
frequently, thus making some pages less relevant for the ILI estimation.
Therefore, we may need to periodically retrain the models by incorporating new
influenza data and additional Wikipedia’s pages to ensure proper stability. We
have seen how just a small subset of pages is responsible for a correct ILI
estimation. Therefore we need to be able to find them correctly. CycleRank and
PPageRank allow us to select new subsets of pages automatically, thus reducing
the chances to have a low-performing model. Besides, we can employ these
methods to effortlessly deploy machine learning models in other languages,
thus making it easy to predict ILI incidence in different countries, as long
as there is a Wikipedia’s version with the same target language. Previously,
this would have required to gather a team of expert in the target language in
order to create a new feature dataset. Ultimately, we argue that this
procedure is suitable to be applied to other illnesses, but further studies
and analysis are required.

We discuss here some of the main shortcomings of our proposed approach that we
leave in future work. Firstly, we employed simple linear regression, since it
is easier to train, and it provides interpretable results. However, by using
linear models, we are implying that we can obtain the response variable, the
ILI incidence, from the linear combination of the predictors. This assumption
does not hold very well in practice, and more complex regression models are
needed, such as Poisson regression or Negative Binomial regression. Secondly,
one of the main issues of Wikipedia-based models is that they are sensitive to
the media coverage of the target topics. This excessive coverage could lead to
overestimation of flu incidence because of the rising public attention. In our
opinion, this is one of the primary concern and limit of this class of methods
which relies on human-generated web traffic. This issue was also the reason
behind the failure of several other similar solutions, like Google Flu Trends
(lazer2014parable, ; santillana2014can, ). We can account for these concerns
by extending our model to take into account the media coverage for the
_“Influenza”_ topic. We could perform anomaly detection on the page views, and
we could instruct the model to normalize the various features to account for
increased traffic. However, it is not trivial to monitor media and the news.
We would need to identify reliable sources, and they should provide open
dataset as Wikipedia does. Moreover, we may need to employ natural language
processing analysis to check if the news subject is semantically related to
ILI. This will require further tests and investigation.

##  8\. Conclusion

The recent pandemics have shown the importance of being able to estimate the
spread of diseases in a fast and reliable way. In our case, we focused on the
estimation of Influenza-Like Illnesses in European countries by using machine
learning models, and the page views coming from Wikipedia, the online
encyclopedia. We extended upon the state of the art by tackling the problem of
using different language editions of Wikipedia that are visited by users in
different European countries. Wikipedia is the first stop for many of health-
related searches, when people are sick and look for information on the
internet, they are likely to end up reading Wikipedia. In this work, we have
shown that we can exploit this information to build linear models to do
nowcasting of the incidence of ILI in a given country. The incidence
estimation concerns both the total number of cases and influenza peak
detection, namely, the week in which we will see the highest number of
infected people in all the influenza season. We can use these results to
direct possible efforts of the public authorities and to devise safety
measures in concert with more traditional approaches (e.g., SIR models,
collection of laboratory tests, etc.). We tested our method on the previous
four influenza seasons (from 2015 to 2019). More importantly, we have also
shown how it is possible to use two automated methods, CycleRank and
PPageRank, to extract automatically and rapidly highly-relevant features from
Wikipedia by looking at Wikipedia’s links structure. In previous works, this
time-consuming task relied upon a group of experts physicians who manually
chose every single relevant Wikipedia’s page. The main advantage of CycleRank
and PPageRank methods relies on their being able to generalize over multiple
languages (and countries) without the need for expert supervision. We compared
the performance of these novel methods against a more traditional one, called
Categories, which uses hand-picked pages. We showed that the models trained
with these newly found features are equal or even outperform the models
trained with hand-picked Wikipedia’s pages.

## Acknowledgments

The authors would like to thank Michele Bortolotti and the team of “Gestione
Sistemi” at the University of Trento for their support with the HPC cluster.

## References

  * (1)
  * (2) Generous N, Fairchild G, Deshpande A, Del Valle SY, Priedhorsky R.  Global disease monitoring and forecasting with Wikipedia.  PLoS computational biology. 2014;10(11):e1003892.
  * (3) Heilman JM, Kemmann E, Bonert M, Chatterjee A, Ragar B, Beards GM, et al.  Wikipedia: a key tool for global public health promotion.  Journal of medical Internet research. 2011;13(1).
  * (4) Heilman JM, West AG.  Wikipedia and medicine: quantifying readership, editors, and the significance of natural language.  Journal of medical Internet research. 2015;17(3).
  * (5) Kim EK, Seok JH, Oh JS, Lee HW, Kim KH.  Use of hangeul twitter to track and predict human influenza infection.  PloS one. 2013;8(7):e69305.
  * (6) Laurent MR, Vickers TJ.  Seeking health information online: does Wikipedia matter?  Journal of the American Medical Informatics Association. 2009;16(4):471–479.
  * (7) Lazer D, Kennedy R, King G, Vespignani A.  The parable of Google Flu: traps in big data analysis.  Science. 2014;343(6176):1203–1205.
  * (8) McIver DJ, Brownstein JS.  Wikipedia usage estimates prevalence of influenza-like illness in the United States in near real-time.  PLoS computational biology. 2014;10(4):e1003581.
  * (9) Nair H, Brooks WA, Katz M, Roca A, Berkley JA, Madhi SA, et al.  Global burden of respiratory infections due to seasonal influenza in young children: a systematic review and meta-analysis.  The Lancet. 2011;378(9807):1917–1930. research agenda.
  * (10) Santillana M, Zhang DW, Althouse BM, Ayers JW.  What can digital disease detection learn from (an external revision to) Google Flu Trends?  American journal of preventive medicine. 2014;47(3):341–347.
  * (11) Santos JC, Matos S.  Analysing Twitter and web queries for flu trend prediction.  Theoretical Biology and Medical Modelling. 2014;11(1):S6.
  * (12) Signorini A, Segre AM, Polgreen PM.  The use of Twitter to track levels of disease activity and public concern in the US during the influenza A H1N1 pandemic.  PloS one. 2011;6(5):e19467.
  * (13) Spoerri A.  What is popular on Wikipedia and why?  First Monday. 2007;12(4).
  * (14) Thompson WW, Shay DK, Weintraub E, Brammer L, Cox N, Anderson LJ, et al.  Mortality associated with influenza and respiratory syncytial virus in the United States.  Jama. 2003;289(2):179–186.
  * (15) Thompson WW, Weintraub E, Dhankhar P, Cheng PY, Brammer L, Meltzer MI, et al.  Estimates of US influenza-associated deaths made using four different methods.  Influenza and other respiratory viruses. 2009;3(1):37–49.
  * (16) Foundation W. Wikimedia Foundation reportcard; 2015.  https://reportcard.wmflabs.org/.
  * (17) Foundation W. English Wikipedia, Special:Statistics; 2016.  https://en.wikipedia.org/wiki/Special:Statistics.
  * (18) Woo H, Cho Y, Shim E, Lee JK, Lee CG, Kim SH.  Estimating Influenza Outbreaks Using Both Search Engine Query Data and Social Media Data in South Korea.  J Med Internet Res. 2016;18(7):e177.  doi:10.2196/jmir.4955.
  * (19) Yang S, Santillana M, Kou SC.  Accurate estimation of influenza epidemics using Google search data via ARGO.  Proceedings of the National Academy of Sciences. 2015;112(47):14473–14478.
  * (20) Yasseri T, Spoerri A, Graham M, Kertesz J.  The Most Controversial Topics in Wikipedia: A Multilingual and Geographical Analysis (May 23, 2013).  Global Wikipedia: International and Cross-Cultural Issues in Online Collaboration, edited by Pnina Fichman, and Noriko Hara Lanham: Rowman & Littlefield. 2014; p. 25–47.
  * (21) Factsheet about seasonal influenza; 2018.  Available from: https://www.ecdc.europa.eu/en/seasonal-influenza/facts/factsheet.
  * (22) ; 2019.  Available from: https://www.who.int/news-room/fact-sheets/detail/influenza-(seasonal).
  * (23) European Influenza Surveillance Network (EISN); 2017.  Available from: https://www.ecdc.europa.eu/en/about-us/partnerships-and-networks/disease-and-laboratory-networks/eisn.
  * (24) Page L, Brin S, Motwani R, Winograd T.  The pagerank citation ranking: Bringing order to the web.  Stanford InfoLab; 1999.
  * (25) De Toni G., Consonni C.  fluTN/influenza-wikipedia-dataset: influenza- wikipedia-dataset-v1.1.0; 2020.  Available from: https://doi.org/10.5281/zenodo.3786254.
  * (26) Consonni C., Laniado D., Montresor A.  CycleRank, or there and back again: personalized relevance scores from cyclic paths on directed graphs.  , Proc. R. Soc. A.47620190740, 2020 http://doi.org/10.1098/rspa.2019.0740
  * (27) West R, Paranjape A, Leskovec J.  Mining Missing Hyperlinks from Human Navigation Traces: A Case Study of Wikipedia.  In: Proceedings of the 24th International Conference on World Wide Web. WWW ’15. Republic and Canton of Geneva, CHE: International World Wide Web Conferences Steering Committee; 2015. p. 1242–1252.  Available from: https://doi.org/10.1145/2736277.2741666.
  * (28) Tibshirani R.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological). 1996;58(1):267–288.
  * (29) Gleich, David F.,  PageRank Beyond the Web,  SIAM Review, 2015;57(3):321-363, doi: 10.1137/140976649
  * (30) Park S., Lee W., Choe B., Lee S.,  A Survey on Personalized PageRank Computation Algorithms.,  IEEE Access, 2019, doi: 10.1109/ACCESS.2019.2952653