Introduction
The aedseo package performs automated and early
detection of seasonal epidemic onsets and estimates the breakpoints for
burden levels from time series data stratified by season. The seasonal
onset (seasonal_onset()) estimates growth rates for
consecutive time intervals of cases and calculates the rolling average
of observations (cases or incidence) in the selected time intervals. The
burden levels (seasonal_burden_levels()) use the previous
seasons to estimate the burden levels of the current season. The
algorithm allows for surveillance of pathogens, by alarming when the
observations have significant growth in the selected time interval and
based on the disease-specific threshold, while also evaluating the
burden of current observations based on previous seasons.
Seasonal data
To apply the aedseo algorithm, data needs to be
transformed into a tsd object. If you have your own data,
the to_time_series() function can be used with the
arguments: cases, incidence,
population, incidence_denominator,
time and time_interval. As default both
seasonal_onset() and seasonal_burden_levels()
use cases, but if incidence is in the
tsd object the output will be in incidence
instead. If population is additionally given as arguments
the function will calculate the incidence per 100.000
(default for incidence_denominator which can be changed in
input). In the following section, the application of the algorithm is
shown with simulated cases created with the
generate_seasonal_data()function. More information about
the function can be found in the
vignette("generate_seasonal_wave").
In the following figure simulated data (solid circles) are visualised as cases over time (weeks). The solid line connects these circles, representing the underlying mean trend over five years of weekly data.
Determining season
Respiratory viruses can circulate in different seasons based on the location. In the nordic hemisphere they mostly circulate in the fall and winter seasons, hence surveillance is intensified from week 40 to week 20 in the following year. To include all data, the season in the example is set from week 21 to week 20 in the following year. In this example burden levels and seasonal onset will be estimated for season 2024/2025.
Determining the disease specific threshold
When cases are low there is a risk that randomness will result in significant growth estimates in isolated periods. To increase the robustness of the method a disease-specific threshold is introduced. It should be set such that subsequent estimates of growth are likely to be significant as well. The disease-specific threshold can be determined by examining continuous periods with sustained significant growth, and determine at what number of cases these events occur.
In this example the disease-specific threshold is determined based on consecutive significant weeks from all available previous seasons. Significant weeks are defined as those with a case count that has a significant positive growth rate.
To capture short-term changes and fluctuations in the cases, a
rolling window of size \(k = 5\) is
used to create subsets of the cases for model fitting, and the
quasipoisson family is used to account for
overdispersion.
The estimate_disease_threshold() function can be used to
automatically estimate the disease-specific threshold. As default it
uses;
- Three time steps that have to be significant for the sequence to be considered in the estimation.
- A max gap time of one time step that is allowed to be non-significant between two significant sequences.
- Three historical seasons.
- The longest sequence closest to the peak.
- A decay factor of 0.8 such that the influence of older seasons diminishes exponentially.
The disease-specific threshold can also be estimated analytically. In
this example it is done with the four available previous seasons. The
seasonal_onset() function can be used for this purpose,
without providing the disease-specific threshold. Then the
consecutive_growth_warnings() function can be used to
create groups with subsequent significant weeks. The data can then be
analysed, else you can use plot/autoplot to visualise the sequences of
significant weeks for each season.
The average_observation_window variable represents the
average of cases over a five-week window, which is used to define the
disease-specific threshold.
tsd_onset <- seasonal_onset(
tsd = tsd_data,
k = 5,
family = "quasipoisson",
na_fraction_allowed = 0.4,
season_start = 21, # Season starts in week 21
season_end = 20, # Season ends in week 20 the following year
only_current_season = FALSE
)
consecutive_gr_warn <- consecutive_growth_warnings(
onset_output = tsd_onset
)
autoplot(
consecutive_gr_warn,
k = 5,
skip_current_season = TRUE
) +
ggplot2::geom_vline(
ggplot2::aes(xintercept = 25, linetype = "Threshold"),
color = "black", linewidth = 0.6
) +
ggplot2::scale_linetype_manual(
name = "",
values = c("Threshold" = "dashed")
)
From the plot above, we observe the length of periods (weeks) with subsequent significant growth rates (y-axis). The season with the longest consecutive period of growth is 2021/2022, lasting 14 weeks. However, since we are determining a threshold specifically for the 2024/2025 season, it’s important to prioritize the most recent seasons. The 2023/2024 season shows two periods of significant growth, with the first being the longest and coinciding closely in timing with the consecutive growth period observed in 2022/2023. We select a disease-specific threshold of 25 to ensure early detection of the seasonal onset while minimizing false positives.
In other words, a season onset is declared when the average case count over five weeks surpasses 25 and is accompanied by a significantly positive growth rate.
Inspect the exact conditions around each detected season start
consecutive_gr_warn |>
dplyr::filter(!is.na(significant_counter)) |>
dplyr::filter(season != max(consecutive_gr_warn$season)) |>
dplyr::group_by(season) |>
dplyr::filter(significant_counter == max(significant_counter)) |>
dplyr::mutate(disease_threshold = average_observations_window,
week = ISOweek::ISOweek(reference_time)) |>
dplyr::select(season, week, disease_threshold)
#> # A tibble: 5 × 3
#> # Groups: season [4]
#> season week disease_threshold
#> <chr> <chr> <dbl>
#> 1 2020/2021 2020-W46 69.8
#> 2 2021/2022 2021-W41 22.6
#> 3 2022/2023 2022-W42 28.2
#> 4 2023/2024 2023-W41 19.6
#> 5 2023/2024 2023-W48 128.By inspecting the output from the above code, the disease-specific
threshold is established at 25 cases.
Run the main algorithm: