From 89007c262d35931d32295dafefff48f3aa325194 Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson Date: Tue, 20 Oct 2026 16:33:43 +0100 Subject: [PATCH 01/11] Sync analysis plan with the current Google Docs draft Update the title, wording fixes, the model-based evaluation section (new covariate list, Gamma GAM equations and mgcv formula) and the Results figure list (WIS instead of CRPS, restructured Figures 3-5) to match the draft. LaTeX typos from the draft are fixed so the equations render (\vu -> \nu, missing closing parentheses, x -> \times). Co-Authored-By: Claude Opus 5.5 --- docs/analysis_plan.qmd | 103 +++++++++++++++++++++++++---------------- 1 file changed, 63 insertions(+), 40 deletions(-) diff --git a/docs/analysis_plan.qmd b/docs/analysis_plan.qmd index 756c939..4fda6c0 100644 --- a/docs/analysis_plan.qmd +++ b/docs/analysis_plan.qmd @@ -1,5 +1,5 @@ --- -title: "Evaluating the influence of incorporating wastewater on forecast performance: a case study applied to COVID-19 Hospitalizations in Germany during the 2024-2025 respiratory virus season" +title: "Evaluating the influence of incorporating wastewater on real-time forecast performance: a case study applied to COVID-19 Hospitalizations in Germany during the 2024-2025 respiratory virus season" authors: "Kaitlyn E. Johnson" format: docx editor: visual @@ -35,7 +35,7 @@ The goal of these analyses will be to assess, in the context of the 2024-2025 re ## What we do here -In this paper, we retrospectively produced 28-day ahead forecasts of incident COVID-19 hospital admissions, with and without incorporating wastewater, for each state in Germany from July 2024 to July 2025 using a semi-mechanistic renewal model, `wwinference` and a baseline ARIMA model. In all cases, we used GitHub's commit history to ensure that we only used the wastewater and hospital admissions data that would have been available as of the forecast date. Retrospective forecasts were evaluated against a final evaluation dataset using proper scoring metrics. We assess the value of incorporating wastewater by comparing the scores from each models overall, across horizons, locations, and forecast dates. To assess the impact of different characteristics of the wastewater surveillance system on forecast performance, we used a generative additive model with wastewater-informed forecast scores as an outcome while including covariate indicators which we hypothesized may have an impact on forecast performance. Based on these analyses, we assess the value of incorporating wastewater for this particular model and discuss areas for future development in wastewater-informed forecasting. +In this paper, we retrospectively produced 28-day ahead forecasts of incident COVID-19 hospital admissions, with and without incorporating wastewater, for each state in Germany from July 2024 to July 2025 using a semi-mechanistic renewal model, `wwinference` and a baseline ARIMA model. In all cases, we used GitHub's commit history to ensure that we only used the wastewater and hospital admissions data that would have been available as of the forecast date. Retrospective forecasts were evaluated against a final evaluation dataset using proper scoring metrics. We assess the value of incorporating wastewater by comparing the scores from each model overall, across horizons, locations, and forecast dates. To assess the impact of different characteristics of the wastewater surveillance system on forecast performance, we used a generalized additive model with wastewater-informed forecast scores as an outcome while including covariate indicators which we hypothesized may have an impact on relative forecast performance. Based on these analyses, we assess the value of incorporating wastewater for this particular model and discuss areas for future development in wastewater-informed forecasting. # Methods @@ -79,7 +79,7 @@ The hospital admissions data provided are indexed by the data of the positive te ## Evaluation -We will use a single evaluation dataset from November X, 2025. It is assume that as of this date, all forecast horizons from the entire period of study (July 1st 2-24 to July 7th 2025) will have reached a stable value. We will evaluate all forecasts using the 7-day rolling sum of hospital admissions. Posterior samples of forecasts will be log-transformed, along with the observations the models are evaluated against. Forecasts will be scored using proper scoring rules, with metrics including the Continuous Ranked Probabilist Score (CRPS), its component breakdown (underprediction, overprediction, dispersion) and interval coverage. +We will use a single evaluation dataset from November X, 2025. It is assume that as of this date, all forecast horizons from the entire period of study (July 1st 2024 to July 7th 2025) will have reached a stable value. We will evaluate all forecasts using the 7-day rolling sum of hospital admissions. Posterior samples of forecasts will be log-transformed, along with the observations the models are evaluated against. Forecasts will be scored using proper scoring rules, with metrics including the Continuous Ranked Probability Score (CRPS), its component breakdown (underprediction, overprediction, dispersion) and interval coverage. Visual comparisons of forecasts overlaid with evaluation data will be used to assess the forecast performance of each model for individual locations. @@ -93,18 +93,57 @@ Results which include the baseline ARIMA model will be presented in the suppleme In addition to stratifying results by forecast date, location, and horizon, we perform an exploratory analysis into the factors that may contribute to the relative performance of the `ww + hosp` model compared to the `hosp only` model. -We will include the following variables with which we would like to understand their impact on relative forecast performance: +To do so, we will model the WIS of the log-transformed predictions and observations as outcome variables, and account for the following variables which might impact relative forecast performance: -1. score of the `hosp only` model (as an offset) -2. number of wastewater treatment plants (sites) -3. population coverage of all the wastewater treatment plants (proportion of the state population on wastewater) -4. average sampling frequency across sites (this is likely very similar so may not be useful) -5. average latency from last collection date to the forecast date -6. minimum latency from last collection date to forecast date -7. variability in wastewater data as measured by the coefficient of variation across sites and time during the calibration period +1. Location +2. Horizon +3. Forecast date +4. Population size of the state +5. Number of wastewater treatment plants (sites) +6. Population coverage of all the wastewater treatment plants (proportion of the state population on wastewater) +7. Average sampling frequency across sites (this is likely very similar so may not be useful) +8. Average latency from last collection date to the forecast date +9. Minimum latency from last collection date to forecast date +10. variability in wastewater data as measured by the coefficient of variation across sites and time during the calibration period Resulting in the following model formulation: +Let $WIS_{h, d, l, m}$ denote the WIS of the forecast at horizon $h$, forecast date $d$, location $l$, and model type $m$. The WIS of the log transformed predictions and observations are non-negative and right-skewed so we treat the WIS as Gamma-distributed with a log link, letting the linear predictor $\eta_{h,d,l,m}$ act multiplicatively on the expected score + +$$ WIS_{h,d,l,m} \sim Gamma(\mu_{h,d,l,m}, \sigma)$$ + +$$ \eta_{h,d,l,m} = \log(\mu_{h,d,l,m}) $$ + +Where the $\eta_{h,d,l,m}$ is modelled as: + +$$ \log(\mu_{h,d,l,m}) = \beta_0 + \beta_{ww}d_m + f(h) + b_l + f(d) + f(pop) + f(n_{sites}) + f(\nu) + f(\rho) + f(\gamma) + f(\phi) + f(\alpha) + d_m \times \left( g(h) + u_l + g(d) + g(pop) + g(n_{sites}) + g(\nu) + g(\rho) + g(\gamma) + g(\phi) + g(\alpha) \right) $$ + +Where $d_m$ is an indicator variable for whether wastewater was used to fit the model, $b_l$, $f(h)$, $f(pop)$, $f(n_{sites})$, $f(\nu)$, $f(\rho)$, and $f(\alpha)$ are the overall effects, respectively, of location, horizon, state population, number of wastewater treatment plants, average wastewater data collection frequency $\nu$, average wastewater latency between latest collection date and forecast date $\rho$, minimum latency between collection date and forecast date $\gamma$, wastewater population coverage $\phi$, and average wastewater data variability $\alpha$. The interaction terms with $d_m$ then reflect the additional effect from incorporating wastewater of each of those variables, after already accounting for the overall impact of that variable on forecast performance across both wastewater and hospital admissions informed models. The formula passed to `mgcv` is provided below. + +``` r +gam_fit <- bam( + wis ~ include_ww + + s(horizon, k = 4) + s(horizon, k = 4, by = include_ww) + + s(location, bs = "re") + s(location, bs = "re", by = include_ww) + + s(forecast_date_num, k = 40) + + s(forecast_date_num, k = 40, by = include_ww) + + # One value per state, so keep this smooth simple + s(state_pop, k = 3) + s(state_pop, k = 3, by = include_ww) + + s(n_sites, k = 5) + s(n_sites, k = 5, by = include_ww) + + s(pop_coverage, k = 5) + s(pop_coverage, k = 5, by = include_ww) + + s(avg_sampling_freq, k = 5) + + s(avg_sampling_freq, k = 5, by = include_ww) + + s(avg_latency, k = 5) + s(avg_latency, k = 5, by = include_ww) + + s(min_latency, k = 5) + s(min_latency, k = 5, by = include_ww) + + s(avg_data_variability, k = 5) + + s(avg_data_variability, k = 5, by = include_ww), + data = data_to_fit, + family = Gamma(link = "log"), + method = "fREML", + discrete = TRUE +) +``` + `CRPS^{ww+hosp}_{h,d,l} ~ \beta + CRPS^{hosponly}_{h,d,l}+ s(numberofsites, k ) + s(popcoverageww, k ) + s(avgfreqsampling, k) + s(avglatencyreporting,k) + s(minlatencyreporting, k) + s(variability_wwdata)` The idea behind this analysis will be that including the CRPS score for the `hosp only` model as an offset accounts for the baseline "difficulty" in that particular location $l$, at forecast date $d$ and horizon $h$. Each of the additional components has some (assumed to be non-linear) positive or negative impact on the relative forecast performance of the wastewater-informed model that we would like to estimate. @@ -123,44 +162,28 @@ Figure 2. Schematic diagram describing the retrospective production of forecasts Q: Do we include baseline in all of the below? Probably not necessary in main text but could make supplement version where we include it -Figure 3. Visual comparison of fit and forecast for a subset of locations and forecast dates (can pick ones where ww improves, wastewater has no impact, wastewater degrades performance) + fit to wastewater data - -3 rows, 4 columns - -row = location - -column 1 = hospital admissions fit with `ww+hosp` vs `hosp_only` - -columns 2:4 = fit to wastewater data by site - -Figure 4. Visual comparison of forecast performance across forecast dates for a subset of locations for horizons 1-7 - -2 rows, 3 columns - -row 1 = forecast eval overlay `ww+hosp` vs `hosp only` - -row 2 = score underlay (bar chart) - -columns = locations +Figure 3. Visual comparison of fit and forecast for a subset of locations and forecast dates for 28 day horizons rows, (row1 is model forecast + data overlay, row 2 are WIS summaries, columns are locations) Supplement: -1. Fig 3 & 4 with baseline ARIMA -2. Fig 4 for additional horizons (7-14, 14-21, 21-28) +1. Fig 3 with baseline ARIMA +2. Add the wastewater fit and forecast and hospital forecast comparison figure to supplement ## Quantitative forecast comparison -Figure 5. "traditional" score comparison: overall, by horizon, by location, by forecast date, by location and forecast date +Statistical significance test results using scoringutils/ mcs R package (same set of forecasts so can use either of these) + +Figure 4. "traditional" score comparison: overall, by horizon, by location, by forecast date, by location and forecast date (make a version with and without ARIMA baseline) -A. CRPS by model summarised across locations and forecast dates, broken down by under/over prediction and dispersion, in the form of a horizontal bar chart for each of the models across all 28 day horizons. +A. WIS by model summarised across locations and forecast dates, broken down by under/over prediction and dispersion, in the form of a horizontal bar chart for each of the models across all 28 day horizons. C. 50th and 90th interval coverage summarised across locations and forecast dates (horizontal bar chart, color = model) B. rWIS by horizon overall -D. CRPS by model stratified by location ordered by CRPS of the hosp-only model (stacked bar chart x-axis = location, y-axis = CRPS breakdown, color = model, shading = CRPS breakdown) +D. WIS by model stratified by location ordered by CRPS of the hosp-only model (stacked bar chart x-axis = location, y-axis = CRPS breakdown, color = model, shading = WIS breakdown) -E. CRPS by model stratified by forecast date (stacked bar chart x-axis = date, y-axis = CRPS breakdown, color = model, shading = CRPS breakdown) +E. WIS by model stratified by forecast date (stacked bar chart x-axis = date, y-axis = WIS breakdown, color = model, shading = WIS breakdown) F. National hospital admissions over time @@ -170,13 +193,13 @@ G. Heatmap of rWIS (relative to hosp only) averaged across horizons for each for Result table from GAM output (what is $\beta$) -Figure 6: Model based eval results +Figure 5: Model based eval results -A. scatterplot of CRPS from wastewater model vs CRPS from hospital admissions only model summarised by location-forecast date +A. scatterplot of WIS from wastewater model vs WIS from hospital admissions only model summarised by location-forecast date -B. scatterplot of CRPS from wastewater vs hosp only summarised by location +B. scatterplot of WIS from wastewater vs hosp only summarised by location -C. scatter plot of CRPS from wastewater vs hosp only summarise by forecast date +C. scatter plot of WIS from wastewater vs hosp only summarise by forecast date D. Partial effect of min latency, avg latency, avg frequency, pop coverage, variability in ww, number of sites From 80e152202a564a6c9922479afc6d3b81d6ca6ec9 Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson Date: Tue, 20 Oct 2026 16:59:40 +0100 Subject: [PATCH 02/11] Remove outdated GAM text, CRPS references; add supplement Remove the CRPS offset formula and text and the brms sentence from the model-based evaluation, and replace remaining CRPS references with WIS. Move all ARIMA comparisons to the supplement so the main text focuses on the head-to-head comparison. Add docs/supplement.qmd with the ARIMA comparison figures and the epidemic phase definition and figures. Co-Authored-By: Claude Opus 5.5 --- docs/analysis_plan.qmd | 21 ++++++++------------- docs/supplement.qmd | 25 +++++++++++++++++++++++++ 2 files changed, 33 insertions(+), 13 deletions(-) create mode 100644 docs/supplement.qmd diff --git a/docs/analysis_plan.qmd b/docs/analysis_plan.qmd index 4fda6c0..659dad1 100644 --- a/docs/analysis_plan.qmd +++ b/docs/analysis_plan.qmd @@ -79,7 +79,7 @@ The hospital admissions data provided are indexed by the data of the positive te ## Evaluation -We will use a single evaluation dataset from November X, 2025. It is assume that as of this date, all forecast horizons from the entire period of study (July 1st 2024 to July 7th 2025) will have reached a stable value. We will evaluate all forecasts using the 7-day rolling sum of hospital admissions. Posterior samples of forecasts will be log-transformed, along with the observations the models are evaluated against. Forecasts will be scored using proper scoring rules, with metrics including the Continuous Ranked Probability Score (CRPS), its component breakdown (underprediction, overprediction, dispersion) and interval coverage. +We will use a single evaluation dataset from November X, 2025. It is assume that as of this date, all forecast horizons from the entire period of study (July 1st 2024 to July 7th 2025) will have reached a stable value. We will evaluate all forecasts using the 7-day rolling sum of hospital admissions. Posterior samples of forecasts will be log-transformed, along with the observations the models are evaluated against. Forecasts will be scored using proper scoring rules, with metrics including the weighted interval score (WIS), its component breakdown (underprediction, overprediction, dispersion) and interval coverage. Visual comparisons of forecasts overlaid with evaluation data will be used to assess the forecast performance of each model for individual locations. @@ -87,7 +87,7 @@ We will present the results by analysing the mean scores by model and then strat When comparing forecast performance with and without wastewater, we will compute the relative WIS compared to the hospital admissions only `wwinference` model. -Results which include the baseline ARIMA model will be presented in the supplement. +The main text will focus on the head-to-head comparison of the `wwinference` models with and without wastewater. All results which include the baseline ARIMA model will be presented in the supplement. ### Model-based evaluation: investigation into drivers of differences in forecast performance with and without wastewater @@ -144,13 +144,9 @@ gam_fit <- bam( ) ``` -`CRPS^{ww+hosp}_{h,d,l} ~ \beta + CRPS^{hosponly}_{h,d,l}+ s(numberofsites, k ) + s(popcoverageww, k ) + s(avgfreqsampling, k) + s(avglatencyreporting,k) + s(minlatencyreporting, k) + s(variability_wwdata)` - -The idea behind this analysis will be that including the CRPS score for the `hosp only` model as an offset accounts for the baseline "difficulty" in that particular location $l$, at forecast date $d$ and horizon $h$. Each of the additional components has some (assumed to be non-linear) positive or negative impact on the relative forecast performance of the wastewater-informed model that we would like to estimate. - We will plot the partial effects as a function of the number of sites, population coverage of wastewater surveillance, the average sampling frequency in the wastewater data and the average latency in reporting of the wastewater data. -The estimated $\beta$ reflects the adjusted additive impact of incorporating wastewater on forecast performance taking into account the other covariates. We will implement this using a linear-link function in the R package `brms`. These exploratory analyses will be used to generate hypotheses into the impact of the characteristics of wastewater surveillance system on the performance of the `wwinference` model. +The estimated $\beta_{ww}$ reflects the adjusted impact of incorporating wastewater on forecast performance taking into account the other covariates. Because of the log link, $\exp(\beta_{ww})$ is a multiplicative effect on the expected WIS. These exploratory analyses will be used to generate hypotheses into the impact of the characteristics of wastewater surveillance system on the performance of the `wwinference` model. # Results @@ -160,20 +156,19 @@ Figure 2. Schematic diagram describing the retrospective production of forecasts ## Visual comparison of forecast performance -Q: Do we include baseline in all of the below? Probably not necessary in main text but could make supplement version where we include it - Figure 3. Visual comparison of fit and forecast for a subset of locations and forecast dates for 28 day horizons rows, (row1 is model forecast + data overlay, row 2 are WIS summaries, columns are locations) Supplement: -1. Fig 3 with baseline ARIMA -2. Add the wastewater fit and forecast and hospital forecast comparison figure to supplement +1. Add the wastewater fit and forecast and hospital forecast comparison figure to supplement + +The version of Figure 3 which includes the baseline ARIMA model is in the supplement (`supplement.qmd`). ## Quantitative forecast comparison Statistical significance test results using scoringutils/ mcs R package (same set of forecasts so can use either of these) -Figure 4. "traditional" score comparison: overall, by horizon, by location, by forecast date, by location and forecast date (make a version with and without ARIMA baseline) +Figure 4. "traditional" score comparison: overall, by horizon, by location, by forecast date, by location and forecast date (the version including the ARIMA baseline is in the supplement) A. WIS by model summarised across locations and forecast dates, broken down by under/over prediction and dispersion, in the form of a horizontal bar chart for each of the models across all 28 day horizons. @@ -181,7 +176,7 @@ C. 50th and 90th interval coverage summarised across locations and forecast date B. rWIS by horizon overall -D. WIS by model stratified by location ordered by CRPS of the hosp-only model (stacked bar chart x-axis = location, y-axis = CRPS breakdown, color = model, shading = WIS breakdown) +D. WIS by model stratified by location ordered by WIS of the hosp-only model (stacked bar chart x-axis = location, y-axis = WIS breakdown, color = model, shading = WIS breakdown) E. WIS by model stratified by forecast date (stacked bar chart x-axis = date, y-axis = WIS breakdown, color = model, shading = WIS breakdown) diff --git a/docs/supplement.qmd b/docs/supplement.qmd new file mode 100644 index 0000000..e40045b --- /dev/null +++ b/docs/supplement.qmd @@ -0,0 +1,25 @@ +--- +title: "Supplement: Evaluating the influence of incorporating wastewater on real-time forecast performance: a case study applied to COVID-19 Hospitalizations in Germany during the 2024-2025 respiratory virus season" +authors: "Kaitlyn E. Johnson" +format: docx +editor: visual +bibliography: references.bib +--- + +# S1. Comparison with the baseline ARIMA model + +Figure S1. Figure 3 with the baseline ARIMA model included: forecast and evaluation data overlay (row 1) and WIS summaries (row 2) for a subset of locations and forecast dates. + +Figure S2. Figure 4 with the baseline ARIMA model included: WIS by model broken down by underprediction, overprediction and dispersion, 50% and 90% interval coverage, and relative WIS by horizon, location and forecast date. + +# S2. Epidemic phase + +We classified each forecast into an epidemic phase based on the change in the final evaluation data over the 28-day forecast period. + +For hospital admissions, we compared the 7-day rolling sum of hospital admissions on the forecast date with the 7-day rolling sum 28 days later, for each location and forecast date. A forecast was classified as increasing if admissions increased by more than 15% and by more than 5 admissions, as decreasing if admissions decreased by more than 15% and by more than 5 admissions, and as stable otherwise. The absolute threshold prevents small fluctuations at low counts from being classified as increasing or decreasing. When there were 0 admissions on the forecast date, the percent change is undefined and the phase was based on the absolute change alone. + +For wastewater, we classified the phase separately for each wastewater treatment plant and lab. We compared the mean log concentration in the 7 days up to the forecast date with the mean log concentration in the 7 days up to 28 days later, and classified the change as increasing (more than 15%), decreasing (less than -15%) or stable. Sites with no observations in either 7-day window were not classified. + +Figure S3. Epidemic phase by forecast date for each location. Top panel: hospital admissions (7-day rolling sum) from the final evaluation dataset, with each forecast date to 7 days later shaded by the hospital admissions phase over the following 28 days. Bottom panels: log wastewater concentrations at each site, shaded by that site's phase. + +# References From 0375f303a3a1fa199cb903f4110d8be24a021197 Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson Date: Tue, 20 Oct 2026 17:03:56 +0100 Subject: [PATCH 03/11] Add epidemic phase definition to the Methods Describe how the epidemic phase is defined for hospital admissions (28-day change beyond 15% and 5 admissions) and for wastewater (28-day change in mean log concentration per site and lab, 15%) in a new Methods subsection, and point the supplement to it instead of repeating the definition. Co-Authored-By: Claude Opus 5.5 --- docs/analysis_plan.qmd | 8 ++++++++ docs/supplement.qmd | 6 +----- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/docs/analysis_plan.qmd b/docs/analysis_plan.qmd index 659dad1..347dbeb 100644 --- a/docs/analysis_plan.qmd +++ b/docs/analysis_plan.qmd @@ -89,6 +89,14 @@ When comparing forecast performance with and without wastewater, we will compute The main text will focus on the head-to-head comparison of the `wwinference` models with and without wastewater. All results which include the baseline ARIMA model will be presented in the supplement. +### Epidemic phase + +To compare forecast performance across epidemic phases, we will classify each forecast into an epidemic phase based on the change in the final evaluation data over the 28-day forecast period. + +For hospital admissions, we will compare the 7-day rolling sum of hospital admissions on the forecast date with the 7-day rolling sum 28 days later, for each location and forecast date. A forecast will be classified as increasing if admissions increased by more than 15% and by more than 5 admissions, as decreasing if admissions decreased by more than 15% and by more than 5 admissions, and as stable otherwise. The absolute threshold prevents small fluctuations at low counts from being classified as increasing or decreasing. When there were 0 admissions on the forecast date, the percent change is undefined and the phase will be based on the absolute change alone. + +For wastewater, we will classify the phase separately for each wastewater treatment plant and lab. We will compare the mean log concentration in the 7 days up to the forecast date with the mean log concentration in the 7 days up to 28 days later, and classify the change as increasing (more than 15%), decreasing (less than -15%) or stable. Sites with no observations in either 7-day window will not be classified. + ### Model-based evaluation: investigation into drivers of differences in forecast performance with and without wastewater In addition to stratifying results by forecast date, location, and horizon, we perform an exploratory analysis into the factors that may contribute to the relative performance of the `ww + hosp` model compared to the `hosp only` model. diff --git a/docs/supplement.qmd b/docs/supplement.qmd index e40045b..24d9e02 100644 --- a/docs/supplement.qmd +++ b/docs/supplement.qmd @@ -14,11 +14,7 @@ Figure S2. Figure 4 with the baseline ARIMA model included: WIS by model broken # S2. Epidemic phase -We classified each forecast into an epidemic phase based on the change in the final evaluation data over the 28-day forecast period. - -For hospital admissions, we compared the 7-day rolling sum of hospital admissions on the forecast date with the 7-day rolling sum 28 days later, for each location and forecast date. A forecast was classified as increasing if admissions increased by more than 15% and by more than 5 admissions, as decreasing if admissions decreased by more than 15% and by more than 5 admissions, and as stable otherwise. The absolute threshold prevents small fluctuations at low counts from being classified as increasing or decreasing. When there were 0 admissions on the forecast date, the percent change is undefined and the phase was based on the absolute change alone. - -For wastewater, we classified the phase separately for each wastewater treatment plant and lab. We compared the mean log concentration in the 7 days up to the forecast date with the mean log concentration in the 7 days up to 28 days later, and classified the change as increasing (more than 15%), decreasing (less than -15%) or stable. Sites with no observations in either 7-day window were not classified. +The epidemic phase of each forecast was defined from the change in the final evaluation data over the 28-day forecast period, for hospital admissions by location and for wastewater by site (see Methods). Figure S3. Epidemic phase by forecast date for each location. Top panel: hospital admissions (7-day rolling sum) from the final evaluation dataset, with each forecast date to 7 days later shaded by the hospital admissions phase over the following 28 days. Bottom panels: log wastewater concentrations at each site, shaded by that site's phase. From 94be4ac5bfce3524e882e1ada9057691a58ae11e Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson Date: Tue, 20 Oct 2026 17:47:28 +0100 Subject: [PATCH 04/11] Incorporate the analyses in #127 into the analysis plan Add a wastewater concentration forecast evaluation aim, methods section and figure; frame the model-based evaluation as exploratory and add hospital epidemic phase interacting with wastewater (plus a less flexible forecast date smooth sensitivity analysis and a latency x forecast date interaction); add bootstrapped intervals, significance tests and MCMC convergence checks to the evaluation. Restructure the figures: a combined schematic (Fig 1), example forecasts by phase with WIS and wastewater fits (Fig 2), the head-to-head summary including relative WIS by phase (Fig 3), the wastewater forecast evaluation (Fig 4) and the exploratory GAM results (Fig 5). Update the supplement to match. Co-Authored-By: Claude Opus 5.5 --- docs/analysis_plan.qmd | 85 ++++++++++++++++++++++++------------------ docs/supplement.qmd | 26 +++++++++++-- 2 files changed, 72 insertions(+), 39 deletions(-) diff --git a/docs/analysis_plan.qmd b/docs/analysis_plan.qmd index 347dbeb..841077f 100644 --- a/docs/analysis_plan.qmd +++ b/docs/analysis_plan.qmd @@ -23,19 +23,21 @@ Wastewater surveillance has shown promise for monitoring for presence/absence, f 1. Larger gap: Lack of evaluation of impact of quantitative levels of pathogen concentration in wastewater on forecast performance: While wastewater has shown promise for detecting presence/absence of vaccine preventable diseases or novel variants, and a number of models have been developed to use quantitative levels of pathogen concentration in wastewater to forecast hospital admissions, few have been evaluated to assess the impact of incorporating wastewater data on forecast performance. 2. Smaller gap: Unknown how a particular model, developed for one setting, generalises to different settings: The `wwinference` R package was developed to produce state-level forecasts of hospital admissions in the United States using state-level hospital admissions data and wastewater data from individual wastewater treatment plants. It is not known how this model generalises to a different context with a similar data structure. The data in Germany is similar to the U.S. in that they collect hospital admissions data at the state-level and measurements of the concentration of SARS-CoV-2 genomes in individual wastewater treatment plants within the state. Unknown how this model will perform when applied to a different setting -3. Even smaller gap: Unknown what factors in the wastewater surveillance system impact forecast performance from this particular model. Heterogeneous characteristics of the wastewater surveillance system across locations in Germany can be used to explore the impact of these factors on relative forecast performance of wastewater. +3. Even smaller gap: Unknown what factors in the wastewater surveillance system impact forecast performance from this particular model. Heterogeneous characteristics of the wastewater surveillance system across locations in Germany can be used to explore, in a hypothesis-generating way, the impact of these factors on relative forecast performance of wastewater. ## Aims The goal of these analyses will be to assess, in the context of the 2024-2025 respiratory virus season in Germany using the `wwinference` model: -1. Quantify the changes in forecast performance of a semi-mechanistic renewal model with and without wastewater data. +1. Quantify the changes in forecast performance of hospital admissions from a semi-mechanistic renewal model with and without wastewater data, overall and by epidemic phase. -2. How do characteristics of the wastewater surveillance system impact the relative forecast performance of the wastewater-informed model? +2. Evaluate the performance of the wastewater-informed model in forecasting wastewater concentrations at individual wastewater treatment plants. + +3. Explore how characteristics of the wastewater surveillance system are associated with the relative forecast performance of the wastewater-informed model (exploratory). ## What we do here -In this paper, we retrospectively produced 28-day ahead forecasts of incident COVID-19 hospital admissions, with and without incorporating wastewater, for each state in Germany from July 2024 to July 2025 using a semi-mechanistic renewal model, `wwinference` and a baseline ARIMA model. In all cases, we used GitHub's commit history to ensure that we only used the wastewater and hospital admissions data that would have been available as of the forecast date. Retrospective forecasts were evaluated against a final evaluation dataset using proper scoring metrics. We assess the value of incorporating wastewater by comparing the scores from each model overall, across horizons, locations, and forecast dates. To assess the impact of different characteristics of the wastewater surveillance system on forecast performance, we used a generalized additive model with wastewater-informed forecast scores as an outcome while including covariate indicators which we hypothesized may have an impact on relative forecast performance. Based on these analyses, we assess the value of incorporating wastewater for this particular model and discuss areas for future development in wastewater-informed forecasting. +In this paper, we retrospectively produced 28-day ahead forecasts of incident COVID-19 hospital admissions, with and without incorporating wastewater, for each state in Germany from July 2024 to July 2025 using a semi-mechanistic renewal model, `wwinference` and a baseline ARIMA model. In all cases, we used GitHub's commit history to ensure that we only used the wastewater and hospital admissions data that would have been available as of the forecast date. Retrospective forecasts were evaluated against a final evaluation dataset using proper scoring metrics. We assess the value of incorporating wastewater by comparing the scores from each model overall, across horizons, locations, forecast dates and epidemic phases. We also evaluate the wastewater-informed model's nowcasts and forecasts of wastewater concentrations at each wastewater treatment plant. As an exploratory analysis of the characteristics of the wastewater surveillance system that may impact forecast performance, we used a generalized additive model with forecast scores as an outcome while including covariate indicators which we hypothesized may have an impact on relative forecast performance. Based on these analyses, we assess the value of incorporating wastewater for this particular model and discuss areas for future development in wastewater-informed forecasting. # Methods @@ -85,7 +87,11 @@ Visual comparisons of forecasts overlaid with evaluation data will be used to as We will present the results by analysing the mean scores by model and then stratifying by location, forecast date, and horizon for the models fit to the vintaged hospitalization datasets. -When comparing forecast performance with and without wastewater, we will compute the relative WIS compared to the hospital admissions only `wwinference` model. +When comparing forecast performance with and without wastewater, we will compute the relative WIS compared to the hospital admissions only `wwinference` model, as the ratio of the mean WIS of the two models. + +To quantify uncertainty in the relative WIS, we will bootstrap the forecasts by resampling location-forecast date pairs with replacement, keeping all horizons of a forecast together since scores across horizons of the same forecast are correlated. We will report percentile 95% intervals from 1,000 bootstrap samples, for the relative WIS overall, by horizon and by epidemic phase. We will also test for differences in forecast performance between the two models using pairwise comparisons in `scoringutils`. + +For each `wwinference` model run we will record MCMC convergence diagnostics (divergent transitions, transitions hitting the maximum tree depth, E-BFMI, R-hat and effective sample size) and flag runs which fail the default diagnostic thresholds in `wwinference`. We will report the proportion of runs flagged, and as a sensitivity analysis repeat the comparison excluding location-forecast dates where either model was flagged. The main text will focus on the head-to-head comparison of the `wwinference` models with and without wastewater. All results which include the baseline ARIMA model will be presented in the supplement. @@ -97,9 +103,13 @@ For hospital admissions, we will compare the 7-day rolling sum of hospital admis For wastewater, we will classify the phase separately for each wastewater treatment plant and lab. We will compare the mean log concentration in the 7 days up to the forecast date with the mean log concentration in the 7 days up to 28 days later, and classify the change as increasing (more than 15%), decreasing (less than -15%) or stable. Sites with no observations in either 7-day window will not be classified. +### Evaluation of wastewater concentration forecasts + +The wastewater-informed model also produces predictions of the wastewater concentrations at each wastewater treatment plant. We will evaluate these against the final flow-normalised wastewater concentration data, separately for each wastewater treatment plant and lab. Predictions for dates on or before the forecast date which were not available at the time of the forecast, either because they were reported after the forecast date or were added to the data later, are treated as nowcasts, and predictions after the forecast date are treated as forecasts. Observations which were available when the model was fit are not scored. Observations below the limit of detection are excluded. Since the predictions are on the log scale, the predicted and observed log concentrations are scored directly, using the WIS, bias and interval coverage. As the hospital admissions only model does not predict wastewater concentrations, this is an absolute evaluation of the wastewater-informed model rather than a comparison between models. We will summarise performance by epidemic phase (using each site's wastewater phase), over time, and by wastewater treatment plant. + ### Model-based evaluation: investigation into drivers of differences in forecast performance with and without wastewater -In addition to stratifying results by forecast date, location, and horizon, we perform an exploratory analysis into the factors that may contribute to the relative performance of the `ww + hosp` model compared to the `hosp only` model. +In addition to stratifying results by forecast date, location, horizon and epidemic phase, we perform an exploratory, hypothesis-generating analysis into the factors that may contribute to the relative performance of the `ww + hosp` model compared to the `hosp only` model. To do so, we will model the WIS of the log-transformed predictions and observations as outcome variables, and account for the following variables which might impact relative forecast performance: @@ -113,10 +123,11 @@ To do so, we will model the WIS of the log-transformed predictions and observati 8. Average latency from last collection date to the forecast date 9. Minimum latency from last collection date to forecast date 10. variability in wastewater data as measured by the coefficient of variation across sites and time during the calibration period +11. Epidemic phase of hospital admissions over the forecast period (increasing, decreasing or stable) Resulting in the following model formulation: -Let $WIS_{h, d, l, m}$ denote the WIS of the forecast at horizon $h$, forecast date $d$, location $l$, and model type $m$. The WIS of the log transformed predictions and observations are non-negative and right-skewed so we treat the WIS as Gamma-distributed with a log link, letting the linear predictor $\eta_{h,d,l,m}$ act multiplicatively on the expected score +Let $WIS_{h, d, l, m}$ denote the WIS of the forecast at horizon $h$, forecast date $d$, location $l$, and model type $m$, and let $p = p(d, l)$ denote the hospital admissions epidemic phase for forecast date $d$ and location $l$. The WIS of the log transformed predictions and observations are non-negative and right-skewed so we treat the WIS as Gamma-distributed with a log link, letting the linear predictor $\eta_{h,d,l,m}$ act multiplicatively on the expected score $$ WIS_{h,d,l,m} \sim Gamma(\mu_{h,d,l,m}, \sigma)$$ @@ -124,13 +135,13 @@ $$ \eta_{h,d,l,m} = \log(\mu_{h,d,l,m}) $$ Where the $\eta_{h,d,l,m}$ is modelled as: -$$ \log(\mu_{h,d,l,m}) = \beta_0 + \beta_{ww}d_m + f(h) + b_l + f(d) + f(pop) + f(n_{sites}) + f(\nu) + f(\rho) + f(\gamma) + f(\phi) + f(\alpha) + d_m \times \left( g(h) + u_l + g(d) + g(pop) + g(n_{sites}) + g(\nu) + g(\rho) + g(\gamma) + g(\phi) + g(\alpha) \right) $$ +$$ \log(\mu_{h,d,l,m}) = \beta_0 + \beta_{ww}d_m + c_p + f(h) + b_l + f(d) + f(pop) + f(n_{sites}) + f(\nu) + f(\rho) + f(\gamma) + f(\phi) + f(\alpha) + d_m \times \left( c^{ww}_p + g(h) + u_l + g(d) + g(pop) + g(n_{sites}) + g(\nu) + g(\rho) + g(\gamma) + g(\phi) + g(\alpha) \right) $$ -Where $d_m$ is an indicator variable for whether wastewater was used to fit the model, $b_l$, $f(h)$, $f(pop)$, $f(n_{sites})$, $f(\nu)$, $f(\rho)$, and $f(\alpha)$ are the overall effects, respectively, of location, horizon, state population, number of wastewater treatment plants, average wastewater data collection frequency $\nu$, average wastewater latency between latest collection date and forecast date $\rho$, minimum latency between collection date and forecast date $\gamma$, wastewater population coverage $\phi$, and average wastewater data variability $\alpha$. The interaction terms with $d_m$ then reflect the additional effect from incorporating wastewater of each of those variables, after already accounting for the overall impact of that variable on forecast performance across both wastewater and hospital admissions informed models. The formula passed to `mgcv` is provided below. +Where $d_m$ is an indicator variable for whether wastewater was used to fit the model, $c_p$ is the overall effect of the epidemic phase (with stable as the reference phase), $b_l$, $f(h)$, $f(pop)$, $f(n_{sites})$, $f(\nu)$, $f(\rho)$, and $f(\alpha)$ are the overall effects, respectively, of location, horizon, state population, number of wastewater treatment plants, average wastewater data collection frequency $\nu$, average wastewater latency between latest collection date and forecast date $\rho$, minimum latency between collection date and forecast date $\gamma$, wastewater population coverage $\phi$, and average wastewater data variability $\alpha$. The interaction terms with $d_m$ then reflect the additional effect from incorporating wastewater of each of those variables, after already accounting for the overall impact of that variable on forecast performance across both wastewater and hospital admissions informed models. In particular, $c^{ww}_p$ is the change in the effect of incorporating wastewater in increasing and decreasing phases relative to the stable phase, so that the relative WIS in phase $p$ (at the reference values of the other covariates) is $\exp(\beta_{ww} + c^{ww}_p)$. This allows us to check whether differences in relative WIS by epidemic phase persist after adjusting for the other covariates. The formula passed to `mgcv` is provided below. ``` r gam_fit <- bam( - wis ~ include_ww + + wis ~ include_ww * hosp_phase + s(horizon, k = 4) + s(horizon, k = 4, by = include_ww) + s(location, bs = "re") + s(location, bs = "re", by = include_ww) + s(forecast_date_num, k = 40) + @@ -152,59 +163,61 @@ gam_fit <- bam( ) ``` -We will plot the partial effects as a function of the number of sites, population coverage of wastewater surveillance, the average sampling frequency in the wastewater data and the average latency in reporting of the wastewater data. +The estimated $\beta_{ww}$ reflects the adjusted impact of incorporating wastewater on forecast performance in the stable phase, taking into account the other covariates. Because of the log link, $\exp(\beta_{ww})$ is a multiplicative effect on the expected WIS. We will report the relative WIS for each epidemic phase, with 95% intervals computed from the model's linear predictor, and plot the multiplicative effects of the wastewater surveillance characteristics on the relative WIS. -The estimated $\beta_{ww}$ reflects the adjusted impact of incorporating wastewater on forecast performance taking into account the other covariates. Because of the log link, $\exp(\beta_{ww})$ is a multiplicative effect on the expected WIS. These exploratory analyses will be used to generate hypotheses into the impact of the characteristics of wastewater surveillance system on the performance of the `wwinference` model. +Because epidemic phases occur in runs of consecutive weeks, the forecast date smooth for the wastewater-informed model, $g(d)$, may absorb some of the epidemic phase effect. We will therefore also fit the model with a less flexible forecast date smooth and compare the estimated epidemic phase effects. Latency was high for most locations at the same forecast dates (around the end of December and beginning of January), so the effect of latency may be confounded by forecast date; we will also consider an interaction between minimum latency and forecast date. As there are only 16 locations, the location-level characteristics (population size, number of sites and population coverage) cannot be fully separated from each other or from the location effects. -# Results +These exploratory analyses will be used to generate hypotheses into the impact of the characteristics of wastewater surveillance system on the performance of the `wwinference` model. -Figure 1. Model schematic? +# Results -Figure 2. Schematic diagram describing the retrospective production of forecasts and the evaluation procedure +Figure 1. Schematic of the model (A) and of the retrospective production of forecasts and the evaluation procedure (B). ## Visual comparison of forecast performance -Figure 3. Visual comparison of fit and forecast for a subset of locations and forecast dates for 28 day horizons rows, (row1 is model forecast + data overlay, row 2 are WIS summaries, columns are locations) +Figure 2. Example forecasts for three location-forecast dates, one each in an increasing, decreasing and stable epidemic phase (columns). Row 1: forecasts of hospital admissions from the `ww + hosp` and `hosp only` models overlaid with the evaluation data, with the background shaded by the hospital admissions epidemic phase. Row 2: WIS of each model by horizon. Row 3: fit and forecast of wastewater concentrations at one wastewater treatment plant in the same location, overlaid with the evaluation data and shaded by that site's epidemic phase. -Supplement: - -1. Add the wastewater fit and forecast and hospital forecast comparison figure to supplement - -The version of Figure 3 which includes the baseline ARIMA model is in the supplement (`supplement.qmd`). +The version of Figure 2 which includes the baseline ARIMA model is in the supplement. ## Quantitative forecast comparison Statistical significance test results using scoringutils/ mcs R package (same set of forecasts so can use either of these) -Figure 4. "traditional" score comparison: overall, by horizon, by location, by forecast date, by location and forecast date (the version including the ARIMA baseline is in the supplement) +Figure 3. Head-to-head comparison of hospital admissions forecasts with and without wastewater, across all locations and forecast dates (the version including the ARIMA baseline is in the supplement). Intervals are bootstrapped 95% intervals. -A. WIS by model summarised across locations and forecast dates, broken down by under/over prediction and dispersion, in the form of a horizontal bar chart for each of the models across all 28 day horizons. +A. WIS by model, broken down by underprediction, overprediction and dispersion (horizontal bar chart) -C. 50th and 90th interval coverage summarised across locations and forecast dates (horizontal bar chart, color = model) +B. Relative WIS by horizon -B. rWIS by horizon overall +C. Relative WIS by epidemic phase -D. WIS by model stratified by location ordered by WIS of the hosp-only model (stacked bar chart x-axis = location, y-axis = WIS breakdown, color = model, shading = WIS breakdown) +D. 50% and 90% interval coverage by horizon for each model -E. WIS by model stratified by forecast date (stacked bar chart x-axis = date, y-axis = WIS breakdown, color = model, shading = WIS breakdown) +E. Bias by epidemic phase for each model -F. National hospital admissions over time +F. Heatmap of relative WIS averaged across horizons for each forecast date and location (x-axis = forecast date, y-axis = location, color = relative WIS) -G. Heatmap of rWIS (relative to hosp only) averaged across horizons for each forecast date and location (x-axis = forecast date, y-axis = location, color = rWIS) +## Evaluation of wastewater concentration forecasts -## Investigation into drivers of differences in forecast performance with and without wastewater +Figure 4. Performance of the wastewater-informed model's nowcasts and forecasts of wastewater concentrations, across all wastewater treatment plants (nowcasts and forecasts shown separately). + +A. Bias by wastewater epidemic phase (increasing, decreasing, stable) -Result table from GAM output (what is $\beta$) +B. Calibration: proportion of observations within each central prediction interval against its nominal coverage -Figure 5: Model based eval results +C. 50% and 90% interval coverage over time, averaged across wastewater treatment plants + +D. 50% and 90% interval coverage for each wastewater treatment plant, grouped by location + +## Investigation into drivers of differences in forecast performance with and without wastewater -A. scatterplot of WIS from wastewater model vs WIS from hospital admissions only model summarised by location-forecast date +Figure 5. Exploratory model-based evaluation results. -B. scatterplot of WIS from wastewater vs hosp only summarised by location +A. Relative WIS for each epidemic phase, adjusted for the other covariates ($\exp(\beta_{ww} + c^{ww}_p)$), with 95% intervals -C. scatter plot of WIS from wastewater vs hosp only summarise by forecast date +B. Multiplicative effect on the relative WIS of the wastewater surveillance characteristics (number of sites, population coverage, average sampling frequency, average and minimum latency, variability in wastewater data) and state population size, on a log scale -D. Partial effect of min latency, avg latency, avg frequency, pop coverage, variability in ww, number of sites +Supplement: location effects, horizon and forecast date smooths, model diagnostics, and the sensitivity analysis with a less flexible forecast date smooth. ## Extensions to still consider: diff --git a/docs/supplement.qmd b/docs/supplement.qmd index 24d9e02..e213bf8 100644 --- a/docs/supplement.qmd +++ b/docs/supplement.qmd @@ -8,14 +8,34 @@ bibliography: references.bib # S1. Comparison with the baseline ARIMA model -Figure S1. Figure 3 with the baseline ARIMA model included: forecast and evaluation data overlay (row 1) and WIS summaries (row 2) for a subset of locations and forecast dates. +Figure S1. Figure 2 with the baseline ARIMA model included: hospital admissions forecasts overlaid with the evaluation data and WIS by horizon, for the same example location-forecast dates. -Figure S2. Figure 4 with the baseline ARIMA model included: WIS by model broken down by underprediction, overprediction and dispersion, 50% and 90% interval coverage, and relative WIS by horizon, location and forecast date. +Figure S2. Figure 3 with the baseline ARIMA model included: WIS by model broken down by underprediction, overprediction and dispersion, relative WIS by horizon and epidemic phase, interval coverage and bias. + +Figure S3. WIS by model stratified by location and by forecast date, broken down by underprediction, overprediction and dispersion, including the baseline ARIMA model, with national hospital admissions over time. # S2. Epidemic phase The epidemic phase of each forecast was defined from the change in the final evaluation data over the 28-day forecast period, for hospital admissions by location and for wastewater by site (see Methods). -Figure S3. Epidemic phase by forecast date for each location. Top panel: hospital admissions (7-day rolling sum) from the final evaluation dataset, with each forecast date to 7 days later shaded by the hospital admissions phase over the following 28 days. Bottom panels: log wastewater concentrations at each site, shaded by that site's phase. +Figure S4. Epidemic phase by forecast date for each location. Top panel: hospital admissions (7-day rolling sum) from the final evaluation dataset, with each forecast date to 7 days later shaded by the hospital admissions phase over the following 28 days. Bottom panels: log wastewater concentrations at each site, shaded by that site's phase. + +# S3. MCMC convergence + +Table S1. Number and proportion of `wwinference` model runs flagged for each convergence diagnostic, by model. + +Figure S5. Relative WIS overall, by horizon and by epidemic phase, excluding location-forecast dates where either model was flagged for convergence issues. + +# S4. Wastewater concentration forecasts by site + +Figure S6. WIS, bias and interval coverage of the wastewater nowcasts and forecasts for each wastewater treatment plant. + +# S5. Model-based evaluation + +Figure S7. Location effects on the relative WIS, and the horizon and forecast date smooths. + +Figure S8. Model diagnostics (residuals, QQ plot, observed against fitted). + +Figure S9. Relative WIS by epidemic phase from the model with a less flexible forecast date smooth. # References From 79770b539d8737a6325d45f53dad8266d60855bc Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson Date: Wed, 21 Oct 2026 10:30:16 +0100 Subject: [PATCH 05/11] small edits --- docs/analysis_plan.qmd | 165 ++++++++++++++++------------------------- 1 file changed, 63 insertions(+), 102 deletions(-) diff --git a/docs/analysis_plan.qmd b/docs/analysis_plan.qmd index 841077f..6144542 100644 --- a/docs/analysis_plan.qmd +++ b/docs/analysis_plan.qmd @@ -1,5 +1,5 @@ --- -title: "Evaluating the influence of incorporating wastewater on real-time forecast performance: a case study applied to COVID-19 Hospitalizations in Germany during the 2024-2025 respiratory virus season" +title: "Evaluating the influence of incorporating wastewater on forecast performance: a case study applied to COVID-19 Hospitalizations in Germany during the 2024-2025 respiratory virus season" authors: "Kaitlyn E. Johnson" format: docx editor: visual @@ -23,43 +23,41 @@ Wastewater surveillance has shown promise for monitoring for presence/absence, f 1. Larger gap: Lack of evaluation of impact of quantitative levels of pathogen concentration in wastewater on forecast performance: While wastewater has shown promise for detecting presence/absence of vaccine preventable diseases or novel variants, and a number of models have been developed to use quantitative levels of pathogen concentration in wastewater to forecast hospital admissions, few have been evaluated to assess the impact of incorporating wastewater data on forecast performance. 2. Smaller gap: Unknown how a particular model, developed for one setting, generalises to different settings: The `wwinference` R package was developed to produce state-level forecasts of hospital admissions in the United States using state-level hospital admissions data and wastewater data from individual wastewater treatment plants. It is not known how this model generalises to a different context with a similar data structure. The data in Germany is similar to the U.S. in that they collect hospital admissions data at the state-level and measurements of the concentration of SARS-CoV-2 genomes in individual wastewater treatment plants within the state. Unknown how this model will perform when applied to a different setting -3. Even smaller gap: Unknown what factors in the wastewater surveillance system impact forecast performance from this particular model. Heterogeneous characteristics of the wastewater surveillance system across locations in Germany can be used to explore, in a hypothesis-generating way, the impact of these factors on relative forecast performance of wastewater. +3. Even smaller gap: Unknown what factors in the wastewater surveillance system impact forecast performance from this particular model. Heterogeneous characteristics of the wastewater surveillance system across locations in Germany can be used to explore the impact of these factors on relative forecast performance of wastewater. ## Aims The goal of these analyses will be to assess, in the context of the 2024-2025 respiratory virus season in Germany using the `wwinference` model: -1. Quantify the changes in forecast performance of hospital admissions from a semi-mechanistic renewal model with and without wastewater data, overall and by epidemic phase. +1. Quantify the effect of incorporating wastewater data on the performance of hospital admissions forecasts from a semi-mechanistic renewal model in Germany. -2. Evaluate the performance of the wastewater-informed model in forecasting wastewater concentrations at individual wastewater treatment plants. +2. Quantify the model's ability to forecast wastewater concentrations. -3. Explore how characteristics of the wastewater surveillance system are associated with the relative forecast performance of the wastewater-informed model (exploratory). +3. Explore how characteristics of the wastewater surveillance system impact the forecast performance of the model. ## What we do here -In this paper, we retrospectively produced 28-day ahead forecasts of incident COVID-19 hospital admissions, with and without incorporating wastewater, for each state in Germany from July 2024 to July 2025 using a semi-mechanistic renewal model, `wwinference` and a baseline ARIMA model. In all cases, we used GitHub's commit history to ensure that we only used the wastewater and hospital admissions data that would have been available as of the forecast date. Retrospective forecasts were evaluated against a final evaluation dataset using proper scoring metrics. We assess the value of incorporating wastewater by comparing the scores from each model overall, across horizons, locations, forecast dates and epidemic phases. We also evaluate the wastewater-informed model's nowcasts and forecasts of wastewater concentrations at each wastewater treatment plant. As an exploratory analysis of the characteristics of the wastewater surveillance system that may impact forecast performance, we used a generalized additive model with forecast scores as an outcome while including covariate indicators which we hypothesized may have an impact on relative forecast performance. Based on these analyses, we assess the value of incorporating wastewater for this particular model and discuss areas for future development in wastewater-informed forecasting. +In this paper, we retrospectively produced 28-day ahead forecasts of the 7-day rolling sum of incident COVID-19 hospital admissions, with and without incorporating wastewater, for each state in Germany from July 2024 to July 2025 using a semi-mechanistic renewal model, `wwinference`. In all cases, we used GitHub's commit history to ensure that we only used the wastewater and hospital admissions data that would have been available as of the forecast date. Hospital admissions and wastewater concentration forecasts were evaluated against a final evaluation dataset using proper scoring metrics. We assess the value of incorporating wastewater by comparing the scores from each model overall, across horizons, locations, forecast dates, and epidemic phase. Forecasts of wastewater concentrations were evaluated using interval coverage and bias metrics. To explore the potential impact of different characteristics of the wastewater surveillance system on forecast performance, we used a generalized additive model (GAM) with all forecast scores as an outcome, with covariates for wastewater presence/absence and other variables we hypothesized may have an impact on forecast performance or the relative forecast performance of incorporating wastewater. We conclude by discussing areas for future development in wastewater-informed forecasting. # Methods ## Data -This analysis will use publicly available data from RKI of [COVID-19 hospital admissions at the state-level](https://github.com/robert-koch-institut/COVID-19-Hospitalisierungen_in_Deutschland/) and [wastewater concentration data](https://github.com/robert-koch-institut/Abwassersurveillance_AMELAG) at the site-level in Germany from July 1st, 2024 to July 7th, 2025. This time period was chosen because starting July 1st, 2024, RKI began regularly committing updated wastewater data to their GitHub. Combined with already regular updates of COVID-19 hospital admissions, this allows us to use the commit history to access exactly the datasets that were available publicly as of each forecast date. To generate retrospective forecasts, we created weekly "vintages", or snapshots of both datasets as they were available as of the forecast date, every Monday throughout this period. Forecasts of hospital admissions are evaluated using a final evaluation dataset accessed on November X, 2025. +This analysis will use publicly available data from RKI of [COVID-19 hospital admissions at the state-level](https://github.com/robert-koch-institut/COVID-19-Hospitalisierungen_in_Deutschland/) and [wastewater concentration data](https://github.com/robert-koch-institut/Abwassersurveillance_AMELAG) at the site-level in Germany from July 1st, 2024 to July 7th, 2025. This time period was chosen because starting July 1st, 2024, RKI began regularly committing updated wastewater data to their GitHub. Combined with already regular updates of COVID-19 hospital admissions, this allows us to use the commit history to access exactly the datasets that were available publicly as of each forecast date. To generate retrospective forecasts, we created weekly "vintages", or snapshots of both datasets as they were available as of the forecast date, every Monday throughout this period. Forecasts of hospital admissions and wastewater concentrations are evaluated using a final evaluation dataset accessed on November X, 2025. ### Hospital admissions data -Hospital admissions data are reported by age group and state as 7-day rolling sums. We fit and generate forecasts for the number of admissions summed across all age groups, at the state-level. We will take the difference between successive days to compute daily hospital admissions incidence, as this format is a required input for the renewal-based model(s) (`wwinference`). These data are indexed by the date of positive test of patients admitted to the hospital, which means the initial reports may have a downward bias due to right-truncation that occurs as a result of the delay between the date of a patient's first positive test and their subsequent hospital admission and reporting of their positive test. Because solving the right-truncation problem is not the main focus of this analysis, but we do want to use the data as it would have been available to us in real-time, we will use the nowcast generated by RKI in the hospital admissions dataset as an input to our model. This nowcast corrects for right-truncation in the initial reported data as of the forecast date. The methods used to produce the nowcast are described in [@Wollfram_et_al]. The forecasts will be evaluated using the 7-day rolling sum of the updated counts of hospital admissions from the dataset available over 3 months following the final forecast date, when we presume the case counts on each date of positive test to be relatively finalized. +Hospital admissions data are reported by age group and state as 7-day rolling sums. We fit and generate forecasts for the number of admissions summed across all age groups, at the state-level. We will take the difference between successive days to compute daily hospital admissions incidence, as this format is a required input for the renewal-based model(s) (`wwinference`). These data are indexed by the date of positive test of patients admitted to the hospital, which means the initial reports may have a downward bias due to right-truncation that occurs as a result of the delay between the date of a patient's first positive test and their subsequent hospital admission and reporting of their positive test. Because solving the right-truncation problem is not the main focus of this analysis, but we do want to use the data as it would have been available to us in real-time, we will use the nowcast generated by RKI in the hospital admissions dataset as an input to our model. This nowcast corrects for right-truncation in the initial reported data as of the forecast date. The methods used to produce the nowcast are described in [@Wollfram_et_al]. The forecasts of hospital admissions will be evaluated using the 7-day rolling sum of the updated counts of hospital admissions from the dataset available over 3 months following the final forecast date, when we presume the case counts on each date of positive test to be relatively finalized. ### Wastewater concentration data -Data on wastewater concentrations of SARS-CoV-2 are available from individual wastewater treatment plants, which we will refer to as sites, indexed by sample collection date. Data are also available aggregated at the state-level. However, we do not plan to use the state-level aggregated data, as these end up getting back-corrected because multiple sites contribute with different reporting delays. The `wwinference` model is designed to be able to fit to the raw observations (concentrations) from multiple wastewater treatment plants, grouping observations by wastewater treatment plant and lab. The wastewater concentration data contains an indicator for a change in lab method, which we will use as a way of grouping sets of observations. Each wastewater treatment plant is contained within a state, and contains the metadata on the catchment area (population size) served by the wastewater treatment plant. We will assume all individuals served by the wastewater treatment plant reside within the state. We will jointly fit to all wastewater treatment plants-lab method combinations within a state. - -To ensure that we are fairly assessing retrospective performance with data that would have been available as of the forecast date, we will use the GitHub commit history to create snapshots of data available as of each forecast date, rather than simply truncating the final dataset. +Data on wastewater concentrations of SARS-CoV-2 are available from individual wastewater treatment plants, which we will refer to as sites, indexed by sample collection date. Data are also available aggregated at the state-level. However, we do not plan to use the state-level aggregated data, as these end up getting back-corrected because multiple sites contribute with different reporting delays. The `wwinference` model is designed to be able to fit to the raw observations (concentrations) from multiple wastewater treatment plants, grouping observations by wastewater treatment plant and lab. The wastewater concentration data contains an indicator for a change in lab method, which we will use as a way of grouping sets of observations, meaning we will essentially treat each change in lab method as a new "lab". Each wastewater treatment plant is contained within a state, and contains the metadata on the catchment area (population size) served by the wastewater treatment plant. We will assume all individuals served by the wastewater treatment plant reside within the state. We will jointly fit to all wastewater treatment plants and lab methods within a state. Each state's model is fit independently. #### Limit of quantification of the wastewater concentrations The limit of quantification (LOQ) of a measured wastewater concentration (genome copies per unit volume) is calculated as the geometric mean of the limit of quantification for each measured gene for measured sample on a particular date and processed in a particular lab. -In the case where half or fewer than half of the measured genes are below the the LOQ for a measured sample, the concentration of that sample is calculated as the geometric mean across the measured genes, using half the LOQ as a replacement for the measured gene below the LOQ. This process is described in detail in [@Abunjela\_et_al, @BMChealth, @RKIGithub] (check this). +In the case where half or fewer than half of the measured genes are below the the LOQ for a measured sample, the concentration of that sample is calculated as the geometric mean across the measured genes, using half the LOQ as a replacement for the measured gene below the LOQ. This process is described in detail in [@Abunjela\_et_al, @BMChealth, @RKIGithub] (check this). In the model fit, values below the LOQ have been flagged and are treated as censored observations, see @wwinference for more details. ## Model specification @@ -69,9 +67,9 @@ We will use the following models to produce forecasts of the 7-day rolling sum o 2. `wwinference` (`hosp only`): This model fits to only to the daily hospital admissions at the state-level. It will be calibrated using the most recent 90 days of data on hospital admissions. -3. Baseline ARIMA (`arima_baseline`): This model will be fit only to the 7-day rolling sum of hospital admissions at the state level, using the `auto.arima()` function in the `forecast` R package, which will perform an optimization of both the model parameters and the autoregressive order $p$, differencing order $d$ and the moving average order $q$. We intend to use this as a reference point for comparing the magnitude of the absolute scores so that we can put the forecast comparison with and without wastewater in context. It will be calibrated using the entire historical admissions time series. +3. Baseline ARIMA (`arima_baseline`): This model will be fit only to the 7-day rolling sum of hospital admissions at the state level, using the `auto.arima()` function in the `forecast` R package, which will perform an optimization of both the model parameters and the autoregressive order $p$, differencing order $d$ and the moving average order $q$. We intend to use this as a reference point for comparing the magnitude of the absolute scores so that we can put the forecast comparison with and without wastewater in context. It will be calibrated using the entire historical admissions time series. Results compared to the ARIMA will be presented in the supplement. -Both semi-mechanistic renewal models, which we refer to as the `wwinference` models, will be implemented using the R package `wwinference`. This model expects the input in terms of daily hospital admissions and returns forecasts of the number of daily hospital admissions. Forecasts will be transformed from daily values to 7-day rolling sums for visualisation and evaluation. Full mathematical details on the model are available in \[\\cite R package\] as well as an abbreviated version provided in S1. Text: Mathematical model +Both semi-mechanistic renewal models, which we refer to as the `wwinference` models, will be implemented using the R package `wwinference`. This model expects the input in terms of daily hospital admissions and returns forecasts of the number of daily hospital admissions and wastewater concentrations. Hospital admissions forecasts will be transformed from daily values to 7-day rolling sums for visualisation and evaluation. Full mathematical details on the model are available in [@wwinference] as well as an abbreviated version provided in S1. Text: Mathematical model. ## Implementation @@ -81,143 +79,106 @@ The hospital admissions data provided are indexed by the data of the positive te ## Evaluation -We will use a single evaluation dataset from November X, 2025. It is assume that as of this date, all forecast horizons from the entire period of study (July 1st 2024 to July 7th 2025) will have reached a stable value. We will evaluate all forecasts using the 7-day rolling sum of hospital admissions. Posterior samples of forecasts will be log-transformed, along with the observations the models are evaluated against. Forecasts will be scored using proper scoring rules, with metrics including the weighted interval score (WIS), its component breakdown (underprediction, overprediction, dispersion) and interval coverage. +We will use a single evaluation dataset from November X, 2025. It is assume that as of this date, all forecast horizons from the entire period of study (July 1st 2-24 to July 7th 2025) will have reached a stable value. We will evaluate all forecasts using the 7-day rolling sum of hospital admissions. Posterior samples of forecasts will be log-transformed, along with the observations the models are evaluated against. Forecasts will be scored using proper scoring rules, with metrics including the Continuous Ranked Probabilist Score (CRPS), its component breakdown (underprediction, overprediction, dispersion) and interval coverage. Visual comparisons of forecasts overlaid with evaluation data will be used to assess the forecast performance of each model for individual locations. We will present the results by analysing the mean scores by model and then stratifying by location, forecast date, and horizon for the models fit to the vintaged hospitalization datasets. -When comparing forecast performance with and without wastewater, we will compute the relative WIS compared to the hospital admissions only `wwinference` model, as the ratio of the mean WIS of the two models. - -To quantify uncertainty in the relative WIS, we will bootstrap the forecasts by resampling location-forecast date pairs with replacement, keeping all horizons of a forecast together since scores across horizons of the same forecast are correlated. We will report percentile 95% intervals from 1,000 bootstrap samples, for the relative WIS overall, by horizon and by epidemic phase. We will also test for differences in forecast performance between the two models using pairwise comparisons in `scoringutils`. - -For each `wwinference` model run we will record MCMC convergence diagnostics (divergent transitions, transitions hitting the maximum tree depth, E-BFMI, R-hat and effective sample size) and flag runs which fail the default diagnostic thresholds in `wwinference`. We will report the proportion of runs flagged, and as a sensitivity analysis repeat the comparison excluding location-forecast dates where either model was flagged. - -The main text will focus on the head-to-head comparison of the `wwinference` models with and without wastewater. All results which include the baseline ARIMA model will be presented in the supplement. - -### Epidemic phase - -To compare forecast performance across epidemic phases, we will classify each forecast into an epidemic phase based on the change in the final evaluation data over the 28-day forecast period. - -For hospital admissions, we will compare the 7-day rolling sum of hospital admissions on the forecast date with the 7-day rolling sum 28 days later, for each location and forecast date. A forecast will be classified as increasing if admissions increased by more than 15% and by more than 5 admissions, as decreasing if admissions decreased by more than 15% and by more than 5 admissions, and as stable otherwise. The absolute threshold prevents small fluctuations at low counts from being classified as increasing or decreasing. When there were 0 admissions on the forecast date, the percent change is undefined and the phase will be based on the absolute change alone. - -For wastewater, we will classify the phase separately for each wastewater treatment plant and lab. We will compare the mean log concentration in the 7 days up to the forecast date with the mean log concentration in the 7 days up to 28 days later, and classify the change as increasing (more than 15%), decreasing (less than -15%) or stable. Sites with no observations in either 7-day window will not be classified. - -### Evaluation of wastewater concentration forecasts +When comparing forecast performance with and without wastewater, we will compute the relative WIS compared to the hospital admissions only `wwinference` model. -The wastewater-informed model also produces predictions of the wastewater concentrations at each wastewater treatment plant. We will evaluate these against the final flow-normalised wastewater concentration data, separately for each wastewater treatment plant and lab. Predictions for dates on or before the forecast date which were not available at the time of the forecast, either because they were reported after the forecast date or were added to the data later, are treated as nowcasts, and predictions after the forecast date are treated as forecasts. Observations which were available when the model was fit are not scored. Observations below the limit of detection are excluded. Since the predictions are on the log scale, the predicted and observed log concentrations are scored directly, using the WIS, bias and interval coverage. As the hospital admissions only model does not predict wastewater concentrations, this is an absolute evaluation of the wastewater-informed model rather than a comparison between models. We will summarise performance by epidemic phase (using each site's wastewater phase), over time, and by wastewater treatment plant. +Results which include the baseline ARIMA model will be presented in the supplement. ### Model-based evaluation: investigation into drivers of differences in forecast performance with and without wastewater -In addition to stratifying results by forecast date, location, horizon and epidemic phase, we perform an exploratory, hypothesis-generating analysis into the factors that may contribute to the relative performance of the `ww + hosp` model compared to the `hosp only` model. +In addition to stratifying results by forecast date, location, and horizon, we perform an exploratory analysis into the factors that may contribute to the relative performance of the `ww + hosp` model compared to the `hosp only` model. -To do so, we will model the WIS of the log-transformed predictions and observations as outcome variables, and account for the following variables which might impact relative forecast performance: +We will include the following variables with which we would like to understand their impact on relative forecast performance: -1. Location -2. Horizon -3. Forecast date -4. Population size of the state -5. Number of wastewater treatment plants (sites) -6. Population coverage of all the wastewater treatment plants (proportion of the state population on wastewater) -7. Average sampling frequency across sites (this is likely very similar so may not be useful) -8. Average latency from last collection date to the forecast date -9. Minimum latency from last collection date to forecast date -10. variability in wastewater data as measured by the coefficient of variation across sites and time during the calibration period -11. Epidemic phase of hospital admissions over the forecast period (increasing, decreasing or stable) +1. score of the `hosp only` model (as an offset) +2. number of wastewater treatment plants (sites) +3. population coverage of all the wastewater treatment plants (proportion of the state population on wastewater) +4. average sampling frequency across sites (this is likely very similar so may not be useful) +5. average latency from last collection date to the forecast date +6. minimum latency from last collection date to forecast date +7. variability in wastewater data as measured by the coefficient of variation across sites and time during the calibration period Resulting in the following model formulation: -Let $WIS_{h, d, l, m}$ denote the WIS of the forecast at horizon $h$, forecast date $d$, location $l$, and model type $m$, and let $p = p(d, l)$ denote the hospital admissions epidemic phase for forecast date $d$ and location $l$. The WIS of the log transformed predictions and observations are non-negative and right-skewed so we treat the WIS as Gamma-distributed with a log link, letting the linear predictor $\eta_{h,d,l,m}$ act multiplicatively on the expected score +`CRPS^{ww+hosp}_{h,d,l} ~ \beta + CRPS^{hosponly}_{h,d,l}+ s(numberofsites, k ) + s(popcoverageww, k ) + s(avgfreqsampling, k) + s(avglatencyreporting,k) + s(minlatencyreporting, k) + s(variability_wwdata)` -$$ WIS_{h,d,l,m} \sim Gamma(\mu_{h,d,l,m}, \sigma)$$ +The idea behind this analysis will be that including the CRPS score for the `hosp only` model as an offset accounts for the baseline "difficulty" in that particular location $l$, at forecast date $d$ and horizon $h$. Each of the additional components has some (assumed to be non-linear) positive or negative impact on the relative forecast performance of the wastewater-informed model that we would like to estimate. -$$ \eta_{h,d,l,m} = \log(\mu_{h,d,l,m}) $$ +We will plot the partial effects as a function of the number of sites, population coverage of wastewater surveillance, the average sampling frequency in the wastewater data and the average latency in reporting of the wastewater data. -Where the $\eta_{h,d,l,m}$ is modelled as: +The estimated $\beta$ reflects the adjusted additive impact of incorporating wastewater on forecast performance taking into account the other covariates. We will implement this using a linear-link function in the R package `brms`. These exploratory analyses will be used to generate hypotheses into the impact of the characteristics of wastewater surveillance system on the performance of the `wwinference` model. -$$ \log(\mu_{h,d,l,m}) = \beta_0 + \beta_{ww}d_m + c_p + f(h) + b_l + f(d) + f(pop) + f(n_{sites}) + f(\nu) + f(\rho) + f(\gamma) + f(\phi) + f(\alpha) + d_m \times \left( c^{ww}_p + g(h) + u_l + g(d) + g(pop) + g(n_{sites}) + g(\nu) + g(\rho) + g(\gamma) + g(\phi) + g(\alpha) \right) $$ +# Results -Where $d_m$ is an indicator variable for whether wastewater was used to fit the model, $c_p$ is the overall effect of the epidemic phase (with stable as the reference phase), $b_l$, $f(h)$, $f(pop)$, $f(n_{sites})$, $f(\nu)$, $f(\rho)$, and $f(\alpha)$ are the overall effects, respectively, of location, horizon, state population, number of wastewater treatment plants, average wastewater data collection frequency $\nu$, average wastewater latency between latest collection date and forecast date $\rho$, minimum latency between collection date and forecast date $\gamma$, wastewater population coverage $\phi$, and average wastewater data variability $\alpha$. The interaction terms with $d_m$ then reflect the additional effect from incorporating wastewater of each of those variables, after already accounting for the overall impact of that variable on forecast performance across both wastewater and hospital admissions informed models. In particular, $c^{ww}_p$ is the change in the effect of incorporating wastewater in increasing and decreasing phases relative to the stable phase, so that the relative WIS in phase $p$ (at the reference values of the other covariates) is $\exp(\beta_{ww} + c^{ww}_p)$. This allows us to check whether differences in relative WIS by epidemic phase persist after adjusting for the other covariates. The formula passed to `mgcv` is provided below. +Figure 1. Model schematic? -``` r -gam_fit <- bam( - wis ~ include_ww * hosp_phase + - s(horizon, k = 4) + s(horizon, k = 4, by = include_ww) + - s(location, bs = "re") + s(location, bs = "re", by = include_ww) + - s(forecast_date_num, k = 40) + - s(forecast_date_num, k = 40, by = include_ww) + - # One value per state, so keep this smooth simple - s(state_pop, k = 3) + s(state_pop, k = 3, by = include_ww) + - s(n_sites, k = 5) + s(n_sites, k = 5, by = include_ww) + - s(pop_coverage, k = 5) + s(pop_coverage, k = 5, by = include_ww) + - s(avg_sampling_freq, k = 5) + - s(avg_sampling_freq, k = 5, by = include_ww) + - s(avg_latency, k = 5) + s(avg_latency, k = 5, by = include_ww) + - s(min_latency, k = 5) + s(min_latency, k = 5, by = include_ww) + - s(avg_data_variability, k = 5) + - s(avg_data_variability, k = 5, by = include_ww), - data = data_to_fit, - family = Gamma(link = "log"), - method = "fREML", - discrete = TRUE -) -``` +Figure 2. Schematic diagram describing the retrospective production of forecasts and the evaluation procedure -The estimated $\beta_{ww}$ reflects the adjusted impact of incorporating wastewater on forecast performance in the stable phase, taking into account the other covariates. Because of the log link, $\exp(\beta_{ww})$ is a multiplicative effect on the expected WIS. We will report the relative WIS for each epidemic phase, with 95% intervals computed from the model's linear predictor, and plot the multiplicative effects of the wastewater surveillance characteristics on the relative WIS. +## Visual comparison of forecast performance -Because epidemic phases occur in runs of consecutive weeks, the forecast date smooth for the wastewater-informed model, $g(d)$, may absorb some of the epidemic phase effect. We will therefore also fit the model with a less flexible forecast date smooth and compare the estimated epidemic phase effects. Latency was high for most locations at the same forecast dates (around the end of December and beginning of January), so the effect of latency may be confounded by forecast date; we will also consider an interaction between minimum latency and forecast date. As there are only 16 locations, the location-level characteristics (population size, number of sites and population coverage) cannot be fully separated from each other or from the location effects. +Q: Do we include baseline in all of the below? Probably not necessary in main text but could make supplement version where we include it -These exploratory analyses will be used to generate hypotheses into the impact of the characteristics of wastewater surveillance system on the performance of the `wwinference` model. +Figure 3. Visual comparison of fit and forecast for a subset of locations and forecast dates (can pick ones where ww improves, wastewater has no impact, wastewater degrades performance) + fit to wastewater data -# Results +3 rows, 4 columns -Figure 1. Schematic of the model (A) and of the retrospective production of forecasts and the evaluation procedure (B). +row = location -## Visual comparison of forecast performance +column 1 = hospital admissions fit with `ww+hosp` vs `hosp_only` -Figure 2. Example forecasts for three location-forecast dates, one each in an increasing, decreasing and stable epidemic phase (columns). Row 1: forecasts of hospital admissions from the `ww + hosp` and `hosp only` models overlaid with the evaluation data, with the background shaded by the hospital admissions epidemic phase. Row 2: WIS of each model by horizon. Row 3: fit and forecast of wastewater concentrations at one wastewater treatment plant in the same location, overlaid with the evaluation data and shaded by that site's epidemic phase. +columns 2:4 = fit to wastewater data by site -The version of Figure 2 which includes the baseline ARIMA model is in the supplement. +Figure 4. Visual comparison of forecast performance across forecast dates for a subset of locations for horizons 1-7 -## Quantitative forecast comparison +2 rows, 3 columns -Statistical significance test results using scoringutils/ mcs R package (same set of forecasts so can use either of these) +row 1 = forecast eval overlay `ww+hosp` vs `hosp only` -Figure 3. Head-to-head comparison of hospital admissions forecasts with and without wastewater, across all locations and forecast dates (the version including the ARIMA baseline is in the supplement). Intervals are bootstrapped 95% intervals. +row 2 = score underlay (bar chart) -A. WIS by model, broken down by underprediction, overprediction and dispersion (horizontal bar chart) +columns = locations -B. Relative WIS by horizon +Supplement: -C. Relative WIS by epidemic phase +1. Fig 3 & 4 with baseline ARIMA +2. Fig 4 for additional horizons (7-14, 14-21, 21-28) -D. 50% and 90% interval coverage by horizon for each model +## Quantitative forecast comparison -E. Bias by epidemic phase for each model +Figure 5. "traditional" score comparison: overall, by horizon, by location, by forecast date, by location and forecast date -F. Heatmap of relative WIS averaged across horizons for each forecast date and location (x-axis = forecast date, y-axis = location, color = relative WIS) +A. CRPS by model summarised across locations and forecast dates, broken down by under/over prediction and dispersion, in the form of a horizontal bar chart for each of the models across all 28 day horizons. -## Evaluation of wastewater concentration forecasts +C. 50th and 90th interval coverage summarised across locations and forecast dates (horizontal bar chart, color = model) -Figure 4. Performance of the wastewater-informed model's nowcasts and forecasts of wastewater concentrations, across all wastewater treatment plants (nowcasts and forecasts shown separately). +B. rWIS by horizon overall -A. Bias by wastewater epidemic phase (increasing, decreasing, stable) +D. CRPS by model stratified by location ordered by CRPS of the hosp-only model (stacked bar chart x-axis = location, y-axis = CRPS breakdown, color = model, shading = CRPS breakdown) -B. Calibration: proportion of observations within each central prediction interval against its nominal coverage +E. CRPS by model stratified by forecast date (stacked bar chart x-axis = date, y-axis = CRPS breakdown, color = model, shading = CRPS breakdown) -C. 50% and 90% interval coverage over time, averaged across wastewater treatment plants +F. National hospital admissions over time -D. 50% and 90% interval coverage for each wastewater treatment plant, grouped by location +G. Heatmap of rWIS (relative to hosp only) averaged across horizons for each forecast date and location (x-axis = forecast date, y-axis = location, color = rWIS) ## Investigation into drivers of differences in forecast performance with and without wastewater -Figure 5. Exploratory model-based evaluation results. +Result table from GAM output (what is $\beta$) + +Figure 6: Model based eval results + +A. scatterplot of CRPS from wastewater model vs CRPS from hospital admissions only model summarised by location-forecast date -A. Relative WIS for each epidemic phase, adjusted for the other covariates ($\exp(\beta_{ww} + c^{ww}_p)$), with 95% intervals +B. scatterplot of CRPS from wastewater vs hosp only summarised by location -B. Multiplicative effect on the relative WIS of the wastewater surveillance characteristics (number of sites, population coverage, average sampling frequency, average and minimum latency, variability in wastewater data) and state population size, on a log scale +C. scatter plot of CRPS from wastewater vs hosp only summarise by forecast date -Supplement: location effects, horizon and forecast date smooths, model diagnostics, and the sensitivity analysis with a less flexible forecast date smooth. +D. Partial effect of min latency, avg latency, avg frequency, pop coverage, variability in ww, number of sites ## Extensions to still consider: From 15861822c0c092ad025fd19c4a465797b08bd713 Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson Date: Wed, 21 Oct 2026 12:20:54 +0100 Subject: [PATCH 06/11] Restore #127 changes to the analysis plan, keeping later edits 79770b5 was made from the original version of the plan, so it reverted the earlier changes. Three-way merge of the original (18fd616), the #127 version (94be4ac) and 79770b5: keeps the edits from 79770b5 (Aims, What we do here, Data, wastewater data, LOQ censoring, ARIMA and model description) and restores the epidemic phase, wastewater evaluation, bootstrap and convergence methods, the updated GAM section and the revised figures. Co-Authored-By: Claude Opus 5.5 --- docs/analysis_plan.qmd | 143 ++++++++++++++++++++++++++--------------- 1 file changed, 90 insertions(+), 53 deletions(-) diff --git a/docs/analysis_plan.qmd b/docs/analysis_plan.qmd index 6144542..003213e 100644 --- a/docs/analysis_plan.qmd +++ b/docs/analysis_plan.qmd @@ -1,5 +1,5 @@ --- -title: "Evaluating the influence of incorporating wastewater on forecast performance: a case study applied to COVID-19 Hospitalizations in Germany during the 2024-2025 respiratory virus season" +title: "Evaluating the influence of incorporating wastewater on real-time forecast performance: a case study applied to COVID-19 Hospitalizations in Germany during the 2024-2025 respiratory virus season" authors: "Kaitlyn E. Johnson" format: docx editor: visual @@ -23,7 +23,7 @@ Wastewater surveillance has shown promise for monitoring for presence/absence, f 1. Larger gap: Lack of evaluation of impact of quantitative levels of pathogen concentration in wastewater on forecast performance: While wastewater has shown promise for detecting presence/absence of vaccine preventable diseases or novel variants, and a number of models have been developed to use quantitative levels of pathogen concentration in wastewater to forecast hospital admissions, few have been evaluated to assess the impact of incorporating wastewater data on forecast performance. 2. Smaller gap: Unknown how a particular model, developed for one setting, generalises to different settings: The `wwinference` R package was developed to produce state-level forecasts of hospital admissions in the United States using state-level hospital admissions data and wastewater data from individual wastewater treatment plants. It is not known how this model generalises to a different context with a similar data structure. The data in Germany is similar to the U.S. in that they collect hospital admissions data at the state-level and measurements of the concentration of SARS-CoV-2 genomes in individual wastewater treatment plants within the state. Unknown how this model will perform when applied to a different setting -3. Even smaller gap: Unknown what factors in the wastewater surveillance system impact forecast performance from this particular model. Heterogeneous characteristics of the wastewater surveillance system across locations in Germany can be used to explore the impact of these factors on relative forecast performance of wastewater. +3. Even smaller gap: Unknown what factors in the wastewater surveillance system impact forecast performance from this particular model. Heterogeneous characteristics of the wastewater surveillance system across locations in Germany can be used to explore, in a hypothesis-generating way, the impact of these factors on relative forecast performance of wastewater. ## Aims @@ -79,106 +79,143 @@ The hospital admissions data provided are indexed by the data of the positive te ## Evaluation -We will use a single evaluation dataset from November X, 2025. It is assume that as of this date, all forecast horizons from the entire period of study (July 1st 2-24 to July 7th 2025) will have reached a stable value. We will evaluate all forecasts using the 7-day rolling sum of hospital admissions. Posterior samples of forecasts will be log-transformed, along with the observations the models are evaluated against. Forecasts will be scored using proper scoring rules, with metrics including the Continuous Ranked Probabilist Score (CRPS), its component breakdown (underprediction, overprediction, dispersion) and interval coverage. +We will use a single evaluation dataset from November X, 2025. It is assume that as of this date, all forecast horizons from the entire period of study (July 1st 2024 to July 7th 2025) will have reached a stable value. We will evaluate all forecasts using the 7-day rolling sum of hospital admissions. Posterior samples of forecasts will be log-transformed, along with the observations the models are evaluated against. Forecasts will be scored using proper scoring rules, with metrics including the weighted interval score (WIS), its component breakdown (underprediction, overprediction, dispersion) and interval coverage. Visual comparisons of forecasts overlaid with evaluation data will be used to assess the forecast performance of each model for individual locations. We will present the results by analysing the mean scores by model and then stratifying by location, forecast date, and horizon for the models fit to the vintaged hospitalization datasets. -When comparing forecast performance with and without wastewater, we will compute the relative WIS compared to the hospital admissions only `wwinference` model. +When comparing forecast performance with and without wastewater, we will compute the relative WIS compared to the hospital admissions only `wwinference` model, as the ratio of the mean WIS of the two models. -Results which include the baseline ARIMA model will be presented in the supplement. +To quantify uncertainty in the relative WIS, we will bootstrap the forecasts by resampling location-forecast date pairs with replacement, keeping all horizons of a forecast together since scores across horizons of the same forecast are correlated. We will report percentile 95% intervals from 1,000 bootstrap samples, for the relative WIS overall, by horizon and by epidemic phase. We will also test for differences in forecast performance between the two models using pairwise comparisons in `scoringutils`. -### Model-based evaluation: investigation into drivers of differences in forecast performance with and without wastewater +For each `wwinference` model run we will record MCMC convergence diagnostics (divergent transitions, transitions hitting the maximum tree depth, E-BFMI, R-hat and effective sample size) and flag runs which fail the default diagnostic thresholds in `wwinference`. We will report the proportion of runs flagged, and as a sensitivity analysis repeat the comparison excluding location-forecast dates where either model was flagged. -In addition to stratifying results by forecast date, location, and horizon, we perform an exploratory analysis into the factors that may contribute to the relative performance of the `ww + hosp` model compared to the `hosp only` model. +The main text will focus on the head-to-head comparison of the `wwinference` models with and without wastewater. All results which include the baseline ARIMA model will be presented in the supplement. -We will include the following variables with which we would like to understand their impact on relative forecast performance: +### Epidemic phase -1. score of the `hosp only` model (as an offset) -2. number of wastewater treatment plants (sites) -3. population coverage of all the wastewater treatment plants (proportion of the state population on wastewater) -4. average sampling frequency across sites (this is likely very similar so may not be useful) -5. average latency from last collection date to the forecast date -6. minimum latency from last collection date to forecast date -7. variability in wastewater data as measured by the coefficient of variation across sites and time during the calibration period +To compare forecast performance across epidemic phases, we will classify each forecast into an epidemic phase based on the change in the final evaluation data over the 28-day forecast period. -Resulting in the following model formulation: +For hospital admissions, we will compare the 7-day rolling sum of hospital admissions on the forecast date with the 7-day rolling sum 28 days later, for each location and forecast date. A forecast will be classified as increasing if admissions increased by more than 15% and by more than 5 admissions, as decreasing if admissions decreased by more than 15% and by more than 5 admissions, and as stable otherwise. The absolute threshold prevents small fluctuations at low counts from being classified as increasing or decreasing. When there were 0 admissions on the forecast date, the percent change is undefined and the phase will be based on the absolute change alone. -`CRPS^{ww+hosp}_{h,d,l} ~ \beta + CRPS^{hosponly}_{h,d,l}+ s(numberofsites, k ) + s(popcoverageww, k ) + s(avgfreqsampling, k) + s(avglatencyreporting,k) + s(minlatencyreporting, k) + s(variability_wwdata)` +For wastewater, we will classify the phase separately for each wastewater treatment plant and lab. We will compare the mean log concentration in the 7 days up to the forecast date with the mean log concentration in the 7 days up to 28 days later, and classify the change as increasing (more than 15%), decreasing (less than -15%) or stable. Sites with no observations in either 7-day window will not be classified. -The idea behind this analysis will be that including the CRPS score for the `hosp only` model as an offset accounts for the baseline "difficulty" in that particular location $l$, at forecast date $d$ and horizon $h$. Each of the additional components has some (assumed to be non-linear) positive or negative impact on the relative forecast performance of the wastewater-informed model that we would like to estimate. +### Evaluation of wastewater concentration forecasts -We will plot the partial effects as a function of the number of sites, population coverage of wastewater surveillance, the average sampling frequency in the wastewater data and the average latency in reporting of the wastewater data. +The wastewater-informed model also produces predictions of the wastewater concentrations at each wastewater treatment plant. We will evaluate these against the final flow-normalised wastewater concentration data, separately for each wastewater treatment plant and lab. Predictions for dates on or before the forecast date which were not available at the time of the forecast, either because they were reported after the forecast date or were added to the data later, are treated as nowcasts, and predictions after the forecast date are treated as forecasts. Observations which were available when the model was fit are not scored. Observations below the limit of detection are excluded. Since the predictions are on the log scale, the predicted and observed log concentrations are scored directly, using the WIS, bias and interval coverage. As the hospital admissions only model does not predict wastewater concentrations, this is an absolute evaluation of the wastewater-informed model rather than a comparison between models. We will summarise performance by epidemic phase (using each site's wastewater phase), over time, and by wastewater treatment plant. -The estimated $\beta$ reflects the adjusted additive impact of incorporating wastewater on forecast performance taking into account the other covariates. We will implement this using a linear-link function in the R package `brms`. These exploratory analyses will be used to generate hypotheses into the impact of the characteristics of wastewater surveillance system on the performance of the `wwinference` model. +### Model-based evaluation: investigation into drivers of differences in forecast performance with and without wastewater -# Results +In addition to stratifying results by forecast date, location, horizon and epidemic phase, we perform an exploratory, hypothesis-generating analysis into the factors that may contribute to the relative performance of the `ww + hosp` model compared to the `hosp only` model. -Figure 1. Model schematic? +To do so, we will model the WIS of the log-transformed predictions and observations as outcome variables, and account for the following variables which might impact relative forecast performance: -Figure 2. Schematic diagram describing the retrospective production of forecasts and the evaluation procedure +1. Location +2. Horizon +3. Forecast date +4. Population size of the state +5. Number of wastewater treatment plants (sites) +6. Population coverage of all the wastewater treatment plants (proportion of the state population on wastewater) +7. Average sampling frequency across sites (this is likely very similar so may not be useful) +8. Average latency from last collection date to the forecast date +9. Minimum latency from last collection date to forecast date +10. variability in wastewater data as measured by the coefficient of variation across sites and time during the calibration period +11. Epidemic phase of hospital admissions over the forecast period (increasing, decreasing or stable) -## Visual comparison of forecast performance +Resulting in the following model formulation: + +Let $WIS_{h, d, l, m}$ denote the WIS of the forecast at horizon $h$, forecast date $d$, location $l$, and model type $m$, and let $p = p(d, l)$ denote the hospital admissions epidemic phase for forecast date $d$ and location $l$. The WIS of the log transformed predictions and observations are non-negative and right-skewed so we treat the WIS as Gamma-distributed with a log link, letting the linear predictor $\eta_{h,d,l,m}$ act multiplicatively on the expected score + +$$ WIS_{h,d,l,m} \sim Gamma(\mu_{h,d,l,m}, \sigma)$$ -Q: Do we include baseline in all of the below? Probably not necessary in main text but could make supplement version where we include it +$$ \eta_{h,d,l,m} = \log(\mu_{h,d,l,m}) $$ -Figure 3. Visual comparison of fit and forecast for a subset of locations and forecast dates (can pick ones where ww improves, wastewater has no impact, wastewater degrades performance) + fit to wastewater data +Where the $\eta_{h,d,l,m}$ is modelled as: -3 rows, 4 columns +$$ \log(\mu_{h,d,l,m}) = \beta_0 + \beta_{ww}d_m + c_p + f(h) + b_l + f(d) + f(pop) + f(n_{sites}) + f(\nu) + f(\rho) + f(\gamma) + f(\phi) + f(\alpha) + d_m \times \left( c^{ww}_p + g(h) + u_l + g(d) + g(pop) + g(n_{sites}) + g(\nu) + g(\rho) + g(\gamma) + g(\phi) + g(\alpha) \right) $$ -row = location +Where $d_m$ is an indicator variable for whether wastewater was used to fit the model, $c_p$ is the overall effect of the epidemic phase (with stable as the reference phase), $b_l$, $f(h)$, $f(pop)$, $f(n_{sites})$, $f(\nu)$, $f(\rho)$, and $f(\alpha)$ are the overall effects, respectively, of location, horizon, state population, number of wastewater treatment plants, average wastewater data collection frequency $\nu$, average wastewater latency between latest collection date and forecast date $\rho$, minimum latency between collection date and forecast date $\gamma$, wastewater population coverage $\phi$, and average wastewater data variability $\alpha$. The interaction terms with $d_m$ then reflect the additional effect from incorporating wastewater of each of those variables, after already accounting for the overall impact of that variable on forecast performance across both wastewater and hospital admissions informed models. In particular, $c^{ww}_p$ is the change in the effect of incorporating wastewater in increasing and decreasing phases relative to the stable phase, so that the relative WIS in phase $p$ (at the reference values of the other covariates) is $\exp(\beta_{ww} + c^{ww}_p)$. This allows us to check whether differences in relative WIS by epidemic phase persist after adjusting for the other covariates. The formula passed to `mgcv` is provided below. -column 1 = hospital admissions fit with `ww+hosp` vs `hosp_only` +``` r +gam_fit <- bam( + wis ~ include_ww * hosp_phase + + s(horizon, k = 4) + s(horizon, k = 4, by = include_ww) + + s(location, bs = "re") + s(location, bs = "re", by = include_ww) + + s(forecast_date_num, k = 40) + + s(forecast_date_num, k = 40, by = include_ww) + + # One value per state, so keep this smooth simple + s(state_pop, k = 3) + s(state_pop, k = 3, by = include_ww) + + s(n_sites, k = 5) + s(n_sites, k = 5, by = include_ww) + + s(pop_coverage, k = 5) + s(pop_coverage, k = 5, by = include_ww) + + s(avg_sampling_freq, k = 5) + + s(avg_sampling_freq, k = 5, by = include_ww) + + s(avg_latency, k = 5) + s(avg_latency, k = 5, by = include_ww) + + s(min_latency, k = 5) + s(min_latency, k = 5, by = include_ww) + + s(avg_data_variability, k = 5) + + s(avg_data_variability, k = 5, by = include_ww), + data = data_to_fit, + family = Gamma(link = "log"), + method = "fREML", + discrete = TRUE +) +``` -columns 2:4 = fit to wastewater data by site +The estimated $\beta_{ww}$ reflects the adjusted impact of incorporating wastewater on forecast performance in the stable phase, taking into account the other covariates. Because of the log link, $\exp(\beta_{ww})$ is a multiplicative effect on the expected WIS. We will report the relative WIS for each epidemic phase, with 95% intervals computed from the model's linear predictor, and plot the multiplicative effects of the wastewater surveillance characteristics on the relative WIS. -Figure 4. Visual comparison of forecast performance across forecast dates for a subset of locations for horizons 1-7 +Because epidemic phases occur in runs of consecutive weeks, the forecast date smooth for the wastewater-informed model, $g(d)$, may absorb some of the epidemic phase effect. We will therefore also fit the model with a less flexible forecast date smooth and compare the estimated epidemic phase effects. Latency was high for most locations at the same forecast dates (around the end of December and beginning of January), so the effect of latency may be confounded by forecast date; we will also consider an interaction between minimum latency and forecast date. As there are only 16 locations, the location-level characteristics (population size, number of sites and population coverage) cannot be fully separated from each other or from the location effects. -2 rows, 3 columns +These exploratory analyses will be used to generate hypotheses into the impact of the characteristics of wastewater surveillance system on the performance of the `wwinference` model. -row 1 = forecast eval overlay `ww+hosp` vs `hosp only` +# Results -row 2 = score underlay (bar chart) +Figure 1. Schematic of the model (A) and of the retrospective production of forecasts and the evaluation procedure (B). -columns = locations +## Visual comparison of forecast performance -Supplement: +Figure 2. Example forecasts for three location-forecast dates, one each in an increasing, decreasing and stable epidemic phase (columns). Row 1: forecasts of hospital admissions from the `ww + hosp` and `hosp only` models overlaid with the evaluation data, with the background shaded by the hospital admissions epidemic phase. Row 2: WIS of each model by horizon. Row 3: fit and forecast of wastewater concentrations at one wastewater treatment plant in the same location, overlaid with the evaluation data and shaded by that site's epidemic phase. -1. Fig 3 & 4 with baseline ARIMA -2. Fig 4 for additional horizons (7-14, 14-21, 21-28) +The version of Figure 2 which includes the baseline ARIMA model is in the supplement. ## Quantitative forecast comparison -Figure 5. "traditional" score comparison: overall, by horizon, by location, by forecast date, by location and forecast date +Statistical significance test results using scoringutils/ mcs R package (same set of forecasts so can use either of these) -A. CRPS by model summarised across locations and forecast dates, broken down by under/over prediction and dispersion, in the form of a horizontal bar chart for each of the models across all 28 day horizons. +Figure 3. Head-to-head comparison of hospital admissions forecasts with and without wastewater, across all locations and forecast dates (the version including the ARIMA baseline is in the supplement). Intervals are bootstrapped 95% intervals. -C. 50th and 90th interval coverage summarised across locations and forecast dates (horizontal bar chart, color = model) +A. WIS by model, broken down by underprediction, overprediction and dispersion (horizontal bar chart) -B. rWIS by horizon overall +B. Relative WIS by horizon -D. CRPS by model stratified by location ordered by CRPS of the hosp-only model (stacked bar chart x-axis = location, y-axis = CRPS breakdown, color = model, shading = CRPS breakdown) +C. Relative WIS by epidemic phase -E. CRPS by model stratified by forecast date (stacked bar chart x-axis = date, y-axis = CRPS breakdown, color = model, shading = CRPS breakdown) +D. 50% and 90% interval coverage by horizon for each model -F. National hospital admissions over time +E. Bias by epidemic phase for each model -G. Heatmap of rWIS (relative to hosp only) averaged across horizons for each forecast date and location (x-axis = forecast date, y-axis = location, color = rWIS) +F. Heatmap of relative WIS averaged across horizons for each forecast date and location (x-axis = forecast date, y-axis = location, color = relative WIS) -## Investigation into drivers of differences in forecast performance with and without wastewater +## Evaluation of wastewater concentration forecasts + +Figure 4. Performance of the wastewater-informed model's nowcasts and forecasts of wastewater concentrations, across all wastewater treatment plants (nowcasts and forecasts shown separately). -Result table from GAM output (what is $\beta$) +A. Bias by wastewater epidemic phase (increasing, decreasing, stable) -Figure 6: Model based eval results +B. Calibration: proportion of observations within each central prediction interval against its nominal coverage + +C. 50% and 90% interval coverage over time, averaged across wastewater treatment plants + +D. 50% and 90% interval coverage for each wastewater treatment plant, grouped by location + +## Investigation into drivers of differences in forecast performance with and without wastewater -A. scatterplot of CRPS from wastewater model vs CRPS from hospital admissions only model summarised by location-forecast date +Figure 5. Exploratory model-based evaluation results. -B. scatterplot of CRPS from wastewater vs hosp only summarised by location +A. Relative WIS for each epidemic phase, adjusted for the other covariates ($\exp(\beta_{ww} + c^{ww}_p)$), with 95% intervals -C. scatter plot of CRPS from wastewater vs hosp only summarise by forecast date +B. Multiplicative effect on the relative WIS of the wastewater surveillance characteristics (number of sites, population coverage, average sampling frequency, average and minimum latency, variability in wastewater data) and state population size, on a log scale -D. Partial effect of min latency, avg latency, avg frequency, pop coverage, variability in ww, number of sites +Supplement: location effects, horizon and forecast date smooths, model diagnostics, and the sensitivity analysis with a less flexible forecast date smooth. ## Extensions to still consider: From 755fdd5f8fd83ca083d35a0ae7de3f4d24918a7a Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson Date: Wed, 21 Oct 2026 12:32:16 +0100 Subject: [PATCH 07/11] Drop "real-time" from the titles; remove repeated ARIMA note Co-Authored-By: Claude Opus 5.5 --- docs/analysis_plan.qmd | 4 ++-- docs/supplement.qmd | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/analysis_plan.qmd b/docs/analysis_plan.qmd index 003213e..88f9bef 100644 --- a/docs/analysis_plan.qmd +++ b/docs/analysis_plan.qmd @@ -1,5 +1,5 @@ --- -title: "Evaluating the influence of incorporating wastewater on real-time forecast performance: a case study applied to COVID-19 Hospitalizations in Germany during the 2024-2025 respiratory virus season" +title: "Evaluating the influence of incorporating wastewater on forecast performance: a case study applied to COVID-19 Hospitalizations in Germany during the 2024-2025 respiratory virus season" authors: "Kaitlyn E. Johnson" format: docx editor: visual @@ -67,7 +67,7 @@ We will use the following models to produce forecasts of the 7-day rolling sum o 2. `wwinference` (`hosp only`): This model fits to only to the daily hospital admissions at the state-level. It will be calibrated using the most recent 90 days of data on hospital admissions. -3. Baseline ARIMA (`arima_baseline`): This model will be fit only to the 7-day rolling sum of hospital admissions at the state level, using the `auto.arima()` function in the `forecast` R package, which will perform an optimization of both the model parameters and the autoregressive order $p$, differencing order $d$ and the moving average order $q$. We intend to use this as a reference point for comparing the magnitude of the absolute scores so that we can put the forecast comparison with and without wastewater in context. It will be calibrated using the entire historical admissions time series. Results compared to the ARIMA will be presented in the supplement. +3. Baseline ARIMA (`arima_baseline`): This model will be fit only to the 7-day rolling sum of hospital admissions at the state level, using the `auto.arima()` function in the `forecast` R package, which will perform an optimization of both the model parameters and the autoregressive order $p$, differencing order $d$ and the moving average order $q$. We intend to use this as a reference point for comparing the magnitude of the absolute scores so that we can put the forecast comparison with and without wastewater in context. It will be calibrated using the entire historical admissions time series. Both semi-mechanistic renewal models, which we refer to as the `wwinference` models, will be implemented using the R package `wwinference`. This model expects the input in terms of daily hospital admissions and returns forecasts of the number of daily hospital admissions and wastewater concentrations. Hospital admissions forecasts will be transformed from daily values to 7-day rolling sums for visualisation and evaluation. Full mathematical details on the model are available in [@wwinference] as well as an abbreviated version provided in S1. Text: Mathematical model. diff --git a/docs/supplement.qmd b/docs/supplement.qmd index e213bf8..095ad67 100644 --- a/docs/supplement.qmd +++ b/docs/supplement.qmd @@ -1,5 +1,5 @@ --- -title: "Supplement: Evaluating the influence of incorporating wastewater on real-time forecast performance: a case study applied to COVID-19 Hospitalizations in Germany during the 2024-2025 respiratory virus season" +title: "Supplement: Evaluating the influence of incorporating wastewater on forecast performance: a case study applied to COVID-19 Hospitalizations in Germany during the 2024-2025 respiratory virus season" authors: "Kaitlyn E. Johnson" format: docx editor: visual From e976b039863936c6a960ef6297a084bd51d4b5a6 Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson Date: Wed, 21 Oct 2026 13:20:18 +0100 Subject: [PATCH 08/11] Add references for the RKI wastewater and hospital admissions data Co-Authored-By: Claude Opus 5.5 --- docs/references.bib | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/docs/references.bib b/docs/references.bib index fb32fef..e1dc034 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -35,3 +35,19 @@ @article{Wolffram2023 number = {8}, } + +@misc{RKI_AMELAG, + author = {{Robert Koch-Institut}}, + title = {Abwassersurveillance {AMELAG}}, + year = {2025}, + howpublished = {GitHub repository}, + url = {https://github.com/robert-koch-institut/Abwassersurveillance_AMELAG} +} + +@misc{RKI_hospitalisations, + author = {{Robert Koch-Institut}}, + title = {{COVID-19}-Hospitalisierungen in Deutschland}, + year = {2025}, + howpublished = {GitHub repository}, + url = {https://github.com/robert-koch-institut/COVID-19-Hospitalisierungen_in_Deutschland} +} From 581e9d4b436bc9cf2d41928a254ef9779510cfdb Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson Date: Wed, 21 Oct 2026 14:10:58 +0100 Subject: [PATCH 09/11] make edits to analysis plan to reflect latest changes --- docs/analysis_plan.qmd | 100 ++++++++++++++++++++++++----------------- docs/references.bib | 49 ++++++++++++++++++++ 2 files changed, 109 insertions(+), 40 deletions(-) diff --git a/docs/analysis_plan.qmd b/docs/analysis_plan.qmd index 88f9bef..7d1fcfb 100644 --- a/docs/analysis_plan.qmd +++ b/docs/analysis_plan.qmd @@ -9,7 +9,7 @@ date: 2025-07-30 # Summary -This document describes the analysis plan to perform a retrospective evaluation of the impact of incorporating wastewater data on the forecast performance of a semi-mechanistic renewal-based model, specifically the model available in the [@wwinference] R package. Forecasts of COVID-19 hospital admissions will be generated weekly and evaluated against the final observed hospital admissions data, summed across all age groups in Germany during the 2024-2025 respiratory virus season. +This document describes the analysis plan to perform a retrospective evaluation of the impact of incorporating wastewater data on the forecast performance of a semi-mechanistic renewal-based model, specifically the model available in the `wwinference` R package [@wwinference; @johnson2026wwinference]. Forecasts of COVID-19 hospital admissions will be generated weekly and evaluated against the final observed hospital admissions data, summed across all age groups in Germany during the 2024-2025 respiratory virus season. # Introduction @@ -47,17 +47,17 @@ This analysis will use publicly available data from RKI of [COVID-19 hospital ad ### Hospital admissions data -Hospital admissions data are reported by age group and state as 7-day rolling sums. We fit and generate forecasts for the number of admissions summed across all age groups, at the state-level. We will take the difference between successive days to compute daily hospital admissions incidence, as this format is a required input for the renewal-based model(s) (`wwinference`). These data are indexed by the date of positive test of patients admitted to the hospital, which means the initial reports may have a downward bias due to right-truncation that occurs as a result of the delay between the date of a patient's first positive test and their subsequent hospital admission and reporting of their positive test. Because solving the right-truncation problem is not the main focus of this analysis, but we do want to use the data as it would have been available to us in real-time, we will use the nowcast generated by RKI in the hospital admissions dataset as an input to our model. This nowcast corrects for right-truncation in the initial reported data as of the forecast date. The methods used to produce the nowcast are described in [@Wollfram_et_al]. The forecasts of hospital admissions will be evaluated using the 7-day rolling sum of the updated counts of hospital admissions from the dataset available over 3 months following the final forecast date, when we presume the case counts on each date of positive test to be relatively finalized. +Hospital admissions data are reported by age group and state as 7-day rolling sums on RKI's GitHub [@RKI_hospitalisations]. We fit and generate forecasts for the number of admissions summed across all age groups, at the state-level. We will take the difference between successive days to compute daily hospital admissions incidence, as this format is a required input for the renewal-based model(s) (`wwinference`). These data are indexed by the date of positive test of patients admitted to the hospital, which means the initial reports may have a downward bias due to right-truncation that occurs as a result of the delay between the date of a patient's first positive test and their subsequent hospital admission and reporting of their positive test. Because solving the right-truncation problem is not the main focus of this analysis, but we do want to use the data as it would have been available to us in real-time, we will use the nowcast generated by RKI in the hospital admissions dataset as an input to our model. This nowcast corrects for right-truncation in the initial reported data as of the forecast date. The methods used to produce the nowcast are described in [@Wolffram2023]. The forecasts of hospital admissions will be evaluated using the 7-day rolling sum of the updated counts of hospital admissions from the dataset available over 3 months following the final forecast date, when we presume the case counts on each date of positive test to be relatively finalized. ### Wastewater concentration data -Data on wastewater concentrations of SARS-CoV-2 are available from individual wastewater treatment plants, which we will refer to as sites, indexed by sample collection date. Data are also available aggregated at the state-level. However, we do not plan to use the state-level aggregated data, as these end up getting back-corrected because multiple sites contribute with different reporting delays. The `wwinference` model is designed to be able to fit to the raw observations (concentrations) from multiple wastewater treatment plants, grouping observations by wastewater treatment plant and lab. The wastewater concentration data contains an indicator for a change in lab method, which we will use as a way of grouping sets of observations, meaning we will essentially treat each change in lab method as a new "lab". Each wastewater treatment plant is contained within a state, and contains the metadata on the catchment area (population size) served by the wastewater treatment plant. We will assume all individuals served by the wastewater treatment plant reside within the state. We will jointly fit to all wastewater treatment plants and lab methods within a state. Each state's model is fit independently. +Data on wastewater concentrations of SARS-CoV-2 are available from individual wastewater treatment plants, which we will refer to as sites, indexed by sample collection date on RKI's GitHub [@RKI_AMELAG]. Data are also available aggregated at the state-level. However, we do not plan to use the state-level aggregated data, as these end up getting back-corrected because multiple sites contribute with different reporting delays. The `wwinference` model is designed to be able to fit to the raw observations (concentrations) from multiple wastewater treatment plants, grouping observations by wastewater treatment plant and lab. The wastewater concentration data contains an indicator for a change in lab method, which we will use as a way of grouping sets of observations, meaning we will essentially treat each change in lab method as a new "lab". Each wastewater treatment plant is contained within a state, and contains the metadata on the catchment area (population size) served by the wastewater treatment plant. We will assume all individuals served by the wastewater treatment plant reside within the state. We will jointly fit to all wastewater treatment plants and lab methods within a state. Each state's model is fit independently. #### Limit of quantification of the wastewater concentrations The limit of quantification (LOQ) of a measured wastewater concentration (genome copies per unit volume) is calculated as the geometric mean of the limit of quantification for each measured gene for measured sample on a particular date and processed in a particular lab. -In the case where half or fewer than half of the measured genes are below the the LOQ for a measured sample, the concentration of that sample is calculated as the geometric mean across the measured genes, using half the LOQ as a replacement for the measured gene below the LOQ. This process is described in detail in [@Abunjela\_et_al, @BMChealth, @RKIGithub] (check this). In the model fit, values below the LOQ have been flagged and are treated as censored observations, see @wwinference for more details. +In the case where half or fewer than half of the measured genes are below the the LOQ for a measured sample, the concentration of that sample is calculated as the geometric mean across the measured genes, using half the LOQ as a replacement for the measured gene below the LOQ. This process is described in detail in [@abunijuela2026; @Radermacher2025, @RKI_AMELAG]. In the model fit, values below the LOQ have been flagged and are treated as censored observations, see [@wwinference; @johnson2026wwinference] for more details. ## Model specification @@ -69,63 +69,83 @@ We will use the following models to produce forecasts of the 7-day rolling sum o 3. Baseline ARIMA (`arima_baseline`): This model will be fit only to the 7-day rolling sum of hospital admissions at the state level, using the `auto.arima()` function in the `forecast` R package, which will perform an optimization of both the model parameters and the autoregressive order $p$, differencing order $d$ and the moving average order $q$. We intend to use this as a reference point for comparing the magnitude of the absolute scores so that we can put the forecast comparison with and without wastewater in context. It will be calibrated using the entire historical admissions time series. +## Implementation + Both semi-mechanistic renewal models, which we refer to as the `wwinference` models, will be implemented using the R package `wwinference`. This model expects the input in terms of daily hospital admissions and returns forecasts of the number of daily hospital admissions and wastewater concentrations. Hospital admissions forecasts will be transformed from daily values to 7-day rolling sums for visualisation and evaluation. Full mathematical details on the model are available in [@wwinference] as well as an abbreviated version provided in S1. Text: Mathematical model. -## Implementation +Each `wwinference` model will be fit using 4 chains with 750 warmup and 500 sampling iterations each, an adaptation target acceptance statistic of 0.95 and a maximum tree depth of 12 (the `wwinference` defaults). For each model run we will record MCMC convergence diagnostics: the number of divergent transitions, the number of transitions hitting the maximum tree depth, the Bayesian fraction of missing information (E-BFMI), the potential scale reduction factor ($\hat{R}$) and the bulk and tail effective sample sizes. ### Addressing right-truncation -The hospital admissions data provided are indexed by the data of the positive test -- which is sometimes well before patients end up hospitalised, leading to a right-truncation problem described in detail in [@Wolffram2023]. Because solving the nowcasting problem is not a focus of this work (see baselinenowcast, wolffram, etc for more details on advancements in nowcasting methods), but we are focused on evaluating retrospectively as realistically as possible what forecasts could have been produced in real-time, we will use the already nowcasted hospital admissions available as of the forecast date as an input into our model. That is, we will treat this nowcast as the observed hospital admissions data. Details of the methodology used to produce this nowcast are available in [@Wolffram2023, RKI data]. +The hospital admissions data provided are indexed by the date of the positive test -- which is sometimes well before patients end up hospitalised, leading to a right-truncation problem described in detail in [@Wolffram2023]. Because solving the nowcasting problem is not a focus of this work (see @johnson2025, @Wolffram2023, for more details on advancements in nowcasting methods), but we are focused on evaluating retrospectively as realistically as possible what forecasts could have been produced in real-time, we will use the already nowcasted hospital admissions available as of the forecast date as an input into our model. That is, we will treat this nowcast as the observed hospital admissions data. Details of the methodology used to produce this nowcast are available in [@Wolffram2023, @RKI_hospitalisations]. ## Evaluation -We will use a single evaluation dataset from November X, 2025. It is assume that as of this date, all forecast horizons from the entire period of study (July 1st 2024 to July 7th 2025) will have reached a stable value. We will evaluate all forecasts using the 7-day rolling sum of hospital admissions. Posterior samples of forecasts will be log-transformed, along with the observations the models are evaluated against. Forecasts will be scored using proper scoring rules, with metrics including the weighted interval score (WIS), its component breakdown (underprediction, overprediction, dispersion) and interval coverage. - -Visual comparisons of forecasts overlaid with evaluation data will be used to assess the forecast performance of each model for individual locations. +We will use a single evaluation dataset from November X, 2025 for both the hospital admissions and wastewater concentration data. It is assume that as of this date, all forecast horizons from the entire period of study (July 1st 2024 to July 7th 2025) will have reached a stable value. We will evaluate all forecasts using the 7-day rolling sum of hospital admissions. Posterior samples of forecasts will be log-transformed, along with the observations the models are evaluated against. Forecasts will be scored using proper scoring rules, with metrics including the weighted interval score (WIS), its component breakdown (underprediction, overprediction, dispersion), interval coverage and bias. -We will present the results by analysing the mean scores by model and then stratifying by location, forecast date, and horizon for the models fit to the vintaged hospitalization datasets. +#### Convergence criteria -When comparing forecast performance with and without wastewater, we will compute the relative WIS compared to the hospital admissions only `wwinference` model, as the ratio of the mean WIS of the two models. +To be included in the analysis, both models will need to pass tests for convergnece diagnostics. Any model fits which will result in exclusion of that location-forecast date from the analysis. We will follow the default thresholds in `wwinference`, a model run will be flagged as failing to converge if any of the following hold: -To quantify uncertainty in the relative WIS, we will bootstrap the forecasts by resampling location-forecast date pairs with replacement, keeping all horizons of a forecast together since scores across horizons of the same forecast are correlated. We will report percentile 95% intervals from 1,000 bootstrap samples, for the relative WIS overall, by horizon and by epidemic phase. We will also test for differences in forecast performance between the two models using pairwise comparisons in `scoringutils`. +1. the number of divergent transitions in any chain is at least 1% of the total number of post-warmup draws across chains, +2. the number of transitions hitting the maximum tree depth in any chain is at least 1% of the total number of post-warmup draws across chains, +3. at least 5% of the model's parameters and generated quantities have $\hat{R} > 1.05$, or +4. the mean E-BFMI across chains is at most 0.2. -For each `wwinference` model run we will record MCMC convergence diagnostics (divergent transitions, transitions hitting the maximum tree depth, E-BFMI, R-hat and effective sample size) and flag runs which fail the default diagnostic thresholds in `wwinference`. We will report the proportion of runs flagged, and as a sensitivity analysis repeat the comparison excluding location-forecast dates where either model was flagged. +We will report the proportion of runs flagged, and as a sensitivity analysis will repeat the comparison including all forecasts even if model convergence criteria were not met. -The main text will focus on the head-to-head comparison of the `wwinference` models with and without wastewater. All results which include the baseline ARIMA model will be presented in the supplement. -### Epidemic phase +#### Epidemic phase classification To compare forecast performance across epidemic phases, we will classify each forecast into an epidemic phase based on the change in the final evaluation data over the 28-day forecast period. For hospital admissions, we will compare the 7-day rolling sum of hospital admissions on the forecast date with the 7-day rolling sum 28 days later, for each location and forecast date. A forecast will be classified as increasing if admissions increased by more than 15% and by more than 5 admissions, as decreasing if admissions decreased by more than 15% and by more than 5 admissions, and as stable otherwise. The absolute threshold prevents small fluctuations at low counts from being classified as increasing or decreasing. When there were 0 admissions on the forecast date, the percent change is undefined and the phase will be based on the absolute change alone. -For wastewater, we will classify the phase separately for each wastewater treatment plant and lab. We will compare the mean log concentration in the 7 days up to the forecast date with the mean log concentration in the 7 days up to 28 days later, and classify the change as increasing (more than 15%), decreasing (less than -15%) or stable. Sites with no observations in either 7-day window will not be classified. +For wastewater, we will classify the phase separately for each wastewater treatment plant. We will compare the mean log concentration in the 7 days up to the forecast date with the mean log concentration in the 7 days up to 28 days later, and classify the change as increasing (more than 15%), decreasing (less than -15%) or stable. Sites with no observations in either 7-day window will not be classified. + + +### Evaluation of hospital admissions forecasts + + +Head-to-head comparison of hospital admissions forecast performance will include all location-forecast dates for which wastewater is present and up-to-date [consider adding definition of this if defined in code] and both models converge (See Convergence criteria above). This means that the forecast set for each model is exactly identical, which allows us to use summaries of the average WIS overall and across strata without introducing bias due to comparing different sets of forecast problems. + +Visual comparisons of forecasts overlaid with evaluation data will be used to assess the forecast performance of each model for individual locations. + +We will present the results by analysing the overall scores by model and then stratifying by location, forecast date, and horizon. + +When comparing hospital admissions forecast performance with and without wastewater, we will compute the relative WIS compared to the hospital admissions only `wwinference` model, as the ratio of the mean WIS of the two models. + +To quantify uncertainty in the relative WIS, we will bootstrap the forecasts by resampling location-forecast date pairs with replacement, keeping all horizons of a forecast together since scores across horizons of the same forecast are correlated. We will report percentile 95% intervals from 1,000 bootstrap samples, for the relative WIS overall, by horizon and by epidemic phase. We will also test for differences in forecast performance between the two models using pairwise comparisons in `scoringutils`. + +The main text will focus on the head-to-head comparison of the `wwinference` models with and without wastewater. All results which include the baseline ARIMA model will be presented in the supplement. ### Evaluation of wastewater concentration forecasts -The wastewater-informed model also produces predictions of the wastewater concentrations at each wastewater treatment plant. We will evaluate these against the final flow-normalised wastewater concentration data, separately for each wastewater treatment plant and lab. Predictions for dates on or before the forecast date which were not available at the time of the forecast, either because they were reported after the forecast date or were added to the data later, are treated as nowcasts, and predictions after the forecast date are treated as forecasts. Observations which were available when the model was fit are not scored. Observations below the limit of detection are excluded. Since the predictions are on the log scale, the predicted and observed log concentrations are scored directly, using the WIS, bias and interval coverage. As the hospital admissions only model does not predict wastewater concentrations, this is an absolute evaluation of the wastewater-informed model rather than a comparison between models. We will summarise performance by epidemic phase (using each site's wastewater phase), over time, and by wastewater treatment plant. +The wastewater-informed model produces predictions of the wastewater concentrations at each wastewater treatment plant. We will evaluate these against the final flow-normalised wastewater concentration data, separately for each wastewater treatment plant and lab. Predictions for dates on or before the forecast date which were not available at the time of the forecast, either because they were reported after the forecast date or were added to the data later, are treated as nowcasts, and predictions after the forecast date are treated as forecasts. Observations which were available when the model was fit are not scored. Observations below the limit of quantification are excluded. Since the predictions are on the log scale, the predicted and observed log concentrations are scored directly, using the bias and interval coverage. As the hospital admissions only model does not predict wastewater concentrations, this is an absolute evaluation of the wastewater-informed model rather than a comparison between models. We will summarise performance by epidemic phase (using each site's wastewater phase), over time, and by state. + -### Model-based evaluation: investigation into drivers of differences in forecast performance with and without wastewater +### Model-based evaluation: investigation into drivers of forecast performance and relative forecast performance with wastewater -In addition to stratifying results by forecast date, location, horizon and epidemic phase, we perform an exploratory, hypothesis-generating analysis into the factors that may contribute to the relative performance of the `ww + hosp` model compared to the `hosp only` model. +In addition to stratifying results by forecast date, location, horizon, and epidemic phase, we perform an exploratory, hypothesis-generating analysis into the factors that may contribute to the forecast performance of both the `ww + hosp` and the `hosp only` model, accounting for variables that impact the forecast performance and including interaction terms which might reveal variables which impact the relative performance of the model with wastewater compared to the model fit without it. -To do so, we will model the WIS of the log-transformed predictions and observations as outcome variables, and account for the following variables which might impact relative forecast performance: +To do so, we will model the WIS of the log-transformed predictions and observations of hospital admissions as outcome variables, and account for the following variables which might impact relative forecast performance: -1. Location -2. Horizon -3. Forecast date -4. Population size of the state -5. Number of wastewater treatment plants (sites) -6. Population coverage of all the wastewater treatment plants (proportion of the state population on wastewater) -7. Average sampling frequency across sites (this is likely very similar so may not be useful) -8. Average latency from last collection date to the forecast date -9. Minimum latency from last collection date to forecast date -10. variability in wastewater data as measured by the coefficient of variation across sites and time during the calibration period -11. Epidemic phase of hospital admissions over the forecast period (increasing, decreasing or stable) +1. Wastewater presence/absence +2. Location +3. Horizon +4. Forecast date +5. Population size of the state +6. Number of wastewater treatment plants (sites) +7. Population coverage of all the wastewater treatment plants (proportion of the state population on wastewater) +8. Average sampling frequency across sites (this is likely very similar so may not be useful) +9. Average latency from last collection date to the forecast date +10. Minimum latency from last collection date to forecast date +11. variability in wastewater data as measured by the coefficient of variation across sites and time during the calibration period +12. Epidemic phase of hospital admissions over the forecast period (increasing, decreasing or stable) Resulting in the following model formulation: -Let $WIS_{h, d, l, m}$ denote the WIS of the forecast at horizon $h$, forecast date $d$, location $l$, and model type $m$, and let $p = p(d, l)$ denote the hospital admissions epidemic phase for forecast date $d$ and location $l$. The WIS of the log transformed predictions and observations are non-negative and right-skewed so we treat the WIS as Gamma-distributed with a log link, letting the linear predictor $\eta_{h,d,l,m}$ act multiplicatively on the expected score +Let $WIS_{h, d, l, m}$ denote the WIS of the forecast at horizon $h$, forecast date $d$, location $l$, and model type $m$, and let $p = p(d, l)$ denote the hospital admissions epidemic phase for forecast date $d$ and location $l$. The WIS of the log transformed predictions and observations are non-negative and right-skewed so we treat the WIS as Gamma-distributed with a log link, letting the linear predictor $\eta_{h,d,l,m}$ act multiplicatively on the expected score. $$ WIS_{h,d,l,m} \sim Gamma(\mu_{h,d,l,m}, \sigma)$$ @@ -183,15 +203,15 @@ Statistical significance test results using scoringutils/ mcs R package (same se Figure 3. Head-to-head comparison of hospital admissions forecasts with and without wastewater, across all locations and forecast dates (the version including the ARIMA baseline is in the supplement). Intervals are bootstrapped 95% intervals. -A. WIS by model, broken down by underprediction, overprediction and dispersion (horizontal bar chart) +A. WIS by model, broken down by underprediction, overprediction and dispersion (horizontal bar chart). Bands indicated 95% intervals on overall WIS. -B. Relative WIS by horizon +B. Relative WIS by horizon, with bootstrapped intervals. -C. Relative WIS by epidemic phase +C. Relative WIS by epidemic phase, with bootstrapped intervals. -D. 50% and 90% interval coverage by horizon for each model +D. 50% and 90% interval coverage by horizon for each model. -E. Bias by epidemic phase for each model +E. Bias by epidemic phase for each model (dot plot centered around 1). F. Heatmap of relative WIS averaged across horizons for each forecast date and location (x-axis = forecast date, y-axis = location, color = relative WIS) @@ -199,13 +219,13 @@ F. Heatmap of relative WIS averaged across horizons for each forecast date and l Figure 4. Performance of the wastewater-informed model's nowcasts and forecasts of wastewater concentrations, across all wastewater treatment plants (nowcasts and forecasts shown separately). -A. Bias by wastewater epidemic phase (increasing, decreasing, stable) +A. Bias by wastewater epidemic phase (increasing, decreasing, stable) (do plot centered around 1). B. Calibration: proportion of observations within each central prediction interval against its nominal coverage -C. 50% and 90% interval coverage over time, averaged across wastewater treatment plants +C. 50% and 90% interval coverage by horizon, averaged across wastewater treatment plants and states. -D. 50% and 90% interval coverage for each wastewater treatment plant, grouped by location +D. 50% and 90% interval coverage over time, averaged across wastewater treatment plants and states. ## Investigation into drivers of differences in forecast performance with and without wastewater diff --git a/docs/references.bib b/docs/references.bib index e1dc034..3112238 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -36,6 +36,46 @@ @article{Wolffram2023 } +@article{abunijuela2026, +title = {Wastewater-based surveillance as a tool for monitoring and estimating COVID-19 incidence and trends: Insights from Germany, 2022–2024}, +journal = {Science of The Total Environment}, +volume = {1018}, +pages = {181290}, +year = {2026}, +issn = {0048-9697}, +doi = {https://doi.org/10.1016/j.scitotenv.2025.181290}, +url = {https://www.sciencedirect.com/science/article/pii/S0048969725029328}, +author = {Susan Abunijela and Peter Pütz and Timo Greiner and Ann-Sophie Lehfeld and Alexander Schattschneider and Udo Buchholz and Jakob Schumacher} +} + +@article{johnson2026wwinference, + title = {Bayesian generative modeling for heterogeneous wastewater data + applied to {COVID-19} forecasting}, + author = {Johnson, Kaitlyn E and Vega Yon, George and Brand, Samuel P C and + Bernal Zelaya, Christian and Bayer, Damon and Volkov, Igor and + Susswein, Zachary and Magee, Andrew and Gostic, Katelyn M and + English, Kayla M and Ghinai, Isaac and Hamlet, Arran and + Olesen, Scott W and Pulliam, Juliet R C and Abbott, Sam and + Morris, Dylan H}, + journal = {medRxiv}, + year = {2026}, + doi = {10.64898/2026.02.23.26346887} +} + +@article{Radermacher2025, + title = {Wastewater as an early indicator for short-term forecasting COVID-19 hospitalization in Germany}, + author = {Radermacher, Jonas and Thiel, Steffen and Kannt, Aimo and Fr{\"o}hlich, Holger}, + journal = {BMC Public Health}, + volume = {25}, + number = {1}, + pages = {2910}, + year = {2025}, + month = aug, + issn = {1471-2458}, + doi = {10.1186/s12889-025-24149-2}, + url = {https://doi.org/10.1186/s12889-025-24149-2} +} + @misc{RKI_AMELAG, author = {{Robert Koch-Institut}}, title = {Abwassersurveillance {AMELAG}}, @@ -51,3 +91,12 @@ @misc{RKI_hospitalisations howpublished = {GitHub repository}, url = {https://github.com/robert-koch-institut/COVID-19-Hospitalisierungen_in_Deutschland} } + +@article{johnson2025, + title = {Baseline nowcasting methods for handling delays in epidemiological data}, + author = {Kaitlyn E. Johnson and Maria Tang and Emily Tyszka and Barbora Nemcova and Daniel Wolffram and Rosa Ergas and Nicholas G. Reich and Sebastian Funk and Jonathon Mellor and Johannes Bracher and Sam Abbott}, + year = {2025}, + journal = {Wellcome Open Research}, + doi = {10.12688/wellcomeopenres.25027.1}, + url = {https://wellcomeopenresearch.org/articles/10-614/v1}, + } From 52b58be68de39ff1e55a4ad4075ccab849206ed5 Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson Date: Wed, 21 Oct 2026 14:52:13 +0100 Subject: [PATCH 10/11] Revise analysis plan and supplement after review Show unadjusted and adjusted relative WIS by epidemic phase side by side (Fig 5A) and explain how they relate; use minimum latency only; move the relative WIS heatmap to the supplement; centre bias on 0; define wastewater horizons relative to the forecast date (including nowcasts); add Table 1; fix headings, the inclusion criterion and the convergence sentence. In the supplement, add the covariate scatterplots and coverage by plant and by state, and show the convergence sensitivity analysis including all forecasts. Co-Authored-By: Claude Opus 5.5 --- docs/analysis_plan.qmd | 44 ++++++++++++++++++++++-------------------- docs/supplement.qmd | 28 +++++++++++++++++---------- 2 files changed, 41 insertions(+), 31 deletions(-) diff --git a/docs/analysis_plan.qmd b/docs/analysis_plan.qmd index 7d1fcfb..ace4c4a 100644 --- a/docs/analysis_plan.qmd +++ b/docs/analysis_plan.qmd @@ -83,9 +83,9 @@ The hospital admissions data provided are indexed by the date of the positive te We will use a single evaluation dataset from November X, 2025 for both the hospital admissions and wastewater concentration data. It is assume that as of this date, all forecast horizons from the entire period of study (July 1st 2024 to July 7th 2025) will have reached a stable value. We will evaluate all forecasts using the 7-day rolling sum of hospital admissions. Posterior samples of forecasts will be log-transformed, along with the observations the models are evaluated against. Forecasts will be scored using proper scoring rules, with metrics including the weighted interval score (WIS), its component breakdown (underprediction, overprediction, dispersion), interval coverage and bias. -#### Convergence criteria +### Convergence criteria -To be included in the analysis, both models will need to pass tests for convergnece diagnostics. Any model fits which will result in exclusion of that location-forecast date from the analysis. We will follow the default thresholds in `wwinference`, a model run will be flagged as failing to converge if any of the following hold: +To be included in the analysis, both models will need to pass tests for convergence diagnostics. Any model fits which fail these criteria will result in exclusion of that location-forecast date from the analysis. We will follow the default thresholds in `wwinference`, a model run will be flagged as failing to converge if any of the following hold: 1. the number of divergent transitions in any chain is at least 1% of the total number of post-warmup draws across chains, 2. the number of transitions hitting the maximum tree depth in any chain is at least 1% of the total number of post-warmup draws across chains, @@ -95,7 +95,7 @@ To be included in the analysis, both models will need to pass tests for convergn We will report the proportion of runs flagged, and as a sensitivity analysis will repeat the comparison including all forecasts even if model convergence criteria were not met. -#### Epidemic phase classification +### Epidemic phase classification To compare forecast performance across epidemic phases, we will classify each forecast into an epidemic phase based on the change in the final evaluation data over the 28-day forecast period. @@ -107,7 +107,7 @@ For wastewater, we will classify the phase separately for each wastewater treatm ### Evaluation of hospital admissions forecasts -Head-to-head comparison of hospital admissions forecast performance will include all location-forecast dates for which wastewater is present and up-to-date [consider adding definition of this if defined in code] and both models converge (See Convergence criteria above). This means that the forecast set for each model is exactly identical, which allows us to use summaries of the average WIS overall and across strata without introducing bias due to comparing different sets of forecast problems. +Head-to-head comparison of hospital admissions forecast performance will include all location-forecast dates for which wastewater data are available and both models converge (See Convergence criteria above). This means that the forecast set for each model is exactly identical, which allows us to use summaries of the average WIS overall and across strata without introducing bias due to comparing different sets of forecast problems. Visual comparisons of forecasts overlaid with evaluation data will be used to assess the forecast performance of each model for individual locations. @@ -121,7 +121,7 @@ The main text will focus on the head-to-head comparison of the `wwinference` mod ### Evaluation of wastewater concentration forecasts -The wastewater-informed model produces predictions of the wastewater concentrations at each wastewater treatment plant. We will evaluate these against the final flow-normalised wastewater concentration data, separately for each wastewater treatment plant and lab. Predictions for dates on or before the forecast date which were not available at the time of the forecast, either because they were reported after the forecast date or were added to the data later, are treated as nowcasts, and predictions after the forecast date are treated as forecasts. Observations which were available when the model was fit are not scored. Observations below the limit of quantification are excluded. Since the predictions are on the log scale, the predicted and observed log concentrations are scored directly, using the bias and interval coverage. As the hospital admissions only model does not predict wastewater concentrations, this is an absolute evaluation of the wastewater-informed model rather than a comparison between models. We will summarise performance by epidemic phase (using each site's wastewater phase), over time, and by state. +The wastewater-informed model produces predictions of the wastewater concentrations at each wastewater treatment plant. We will evaluate these against the final flow-normalised wastewater concentration data, separately for each wastewater treatment plant and lab. Predictions for dates on or before the forecast date which were not available at the time of the forecast, either because they were reported after the forecast date or were added to the data later, are treated as nowcasts, and predictions after the forecast date are treated as forecasts. Observations which were available when the model was fit are not scored. Observations below the limit of quantification are excluded. Since the predictions are on the log scale, the predicted and observed log concentrations are scored directly, using the bias and interval coverage. As the hospital admissions only model does not predict wastewater concentrations, this is an absolute evaluation of the wastewater-informed model rather than a comparison between models. Horizons for the wastewater predictions are defined relative to the forecast date, so nowcasts have horizons of 0 or less and forecasts have positive horizons. In the main text we will summarise performance by epidemic phase (using each site's wastewater phase), by horizon and over forecast dates, averaged across wastewater treatment plants. In the supplement we will summarise interval coverage for each wastewater treatment plant and for each state. ### Model-based evaluation: investigation into drivers of forecast performance and relative forecast performance with wastewater @@ -138,10 +138,9 @@ To do so, we will model the WIS of the log-transformed predictions and observati 6. Number of wastewater treatment plants (sites) 7. Population coverage of all the wastewater treatment plants (proportion of the state population on wastewater) 8. Average sampling frequency across sites (this is likely very similar so may not be useful) -9. Average latency from last collection date to the forecast date -10. Minimum latency from last collection date to forecast date -11. variability in wastewater data as measured by the coefficient of variation across sites and time during the calibration period -12. Epidemic phase of hospital admissions over the forecast period (increasing, decreasing or stable) +9. Minimum latency from last collection date to forecast date across wastewater treatment plants +10. variability in wastewater data as measured by the coefficient of variation across sites and time during the calibration period +11. Epidemic phase of hospital admissions over the forecast period (increasing, decreasing or stable) Resulting in the following model formulation: @@ -153,9 +152,9 @@ $$ \eta_{h,d,l,m} = \log(\mu_{h,d,l,m}) $$ Where the $\eta_{h,d,l,m}$ is modelled as: -$$ \log(\mu_{h,d,l,m}) = \beta_0 + \beta_{ww}d_m + c_p + f(h) + b_l + f(d) + f(pop) + f(n_{sites}) + f(\nu) + f(\rho) + f(\gamma) + f(\phi) + f(\alpha) + d_m \times \left( c^{ww}_p + g(h) + u_l + g(d) + g(pop) + g(n_{sites}) + g(\nu) + g(\rho) + g(\gamma) + g(\phi) + g(\alpha) \right) $$ +$$ \log(\mu_{h,d,l,m}) = \beta_0 + \beta_{ww}d_m + c_p + f(h) + b_l + f(d) + f(pop) + f(n_{sites}) + f(\nu) + f(\gamma) + f(\phi) + f(\alpha) + d_m \times \left( c^{ww}_p + g(h) + u_l + g(d) + g(pop) + g(n_{sites}) + g(\nu) + g(\gamma) + g(\phi) + g(\alpha) \right) $$ -Where $d_m$ is an indicator variable for whether wastewater was used to fit the model, $c_p$ is the overall effect of the epidemic phase (with stable as the reference phase), $b_l$, $f(h)$, $f(pop)$, $f(n_{sites})$, $f(\nu)$, $f(\rho)$, and $f(\alpha)$ are the overall effects, respectively, of location, horizon, state population, number of wastewater treatment plants, average wastewater data collection frequency $\nu$, average wastewater latency between latest collection date and forecast date $\rho$, minimum latency between collection date and forecast date $\gamma$, wastewater population coverage $\phi$, and average wastewater data variability $\alpha$. The interaction terms with $d_m$ then reflect the additional effect from incorporating wastewater of each of those variables, after already accounting for the overall impact of that variable on forecast performance across both wastewater and hospital admissions informed models. In particular, $c^{ww}_p$ is the change in the effect of incorporating wastewater in increasing and decreasing phases relative to the stable phase, so that the relative WIS in phase $p$ (at the reference values of the other covariates) is $\exp(\beta_{ww} + c^{ww}_p)$. This allows us to check whether differences in relative WIS by epidemic phase persist after adjusting for the other covariates. The formula passed to `mgcv` is provided below. +Where $d_m$ is an indicator variable for whether wastewater was used to fit the model, $c_p$ is the overall effect of the epidemic phase (with stable as the reference phase), $b_l$, $f(h)$, $f(pop)$, $f(n_{sites})$, $f(\nu)$, $f(\gamma)$, $f(\phi)$ and $f(\alpha)$ are the overall effects, respectively, of location, horizon, state population, number of wastewater treatment plants, average wastewater data collection frequency $\nu$, minimum latency between the latest collection date and the forecast date $\gamma$, wastewater population coverage $\phi$, and average wastewater data variability $\alpha$. The interaction terms with $d_m$ then reflect the additional effect from incorporating wastewater of each of those variables, after already accounting for the overall impact of that variable on forecast performance across both wastewater and hospital admissions informed models. In particular, $c^{ww}_p$ is the change in the effect of incorporating wastewater in increasing and decreasing phases relative to the stable phase, so that the relative WIS in phase $p$ (at the reference values of the other covariates) is $\exp(\beta_{ww} + c^{ww}_p)$. This allows us to check whether differences in relative WIS by epidemic phase persist after adjusting for the other covariates. The formula passed to `mgcv` is provided below. ``` r gam_fit <- bam( @@ -170,7 +169,6 @@ gam_fit <- bam( s(pop_coverage, k = 5) + s(pop_coverage, k = 5, by = include_ww) + s(avg_sampling_freq, k = 5) + s(avg_sampling_freq, k = 5, by = include_ww) + - s(avg_latency, k = 5) + s(avg_latency, k = 5, by = include_ww) + s(min_latency, k = 5) + s(min_latency, k = 5, by = include_ww) + s(avg_data_variability, k = 5) + s(avg_data_variability, k = 5, by = include_ww), @@ -183,7 +181,9 @@ gam_fit <- bam( The estimated $\beta_{ww}$ reflects the adjusted impact of incorporating wastewater on forecast performance in the stable phase, taking into account the other covariates. Because of the log link, $\exp(\beta_{ww})$ is a multiplicative effect on the expected WIS. We will report the relative WIS for each epidemic phase, with 95% intervals computed from the model's linear predictor, and plot the multiplicative effects of the wastewater surveillance characteristics on the relative WIS. -Because epidemic phases occur in runs of consecutive weeks, the forecast date smooth for the wastewater-informed model, $g(d)$, may absorb some of the epidemic phase effect. We will therefore also fit the model with a less flexible forecast date smooth and compare the estimated epidemic phase effects. Latency was high for most locations at the same forecast dates (around the end of December and beginning of January), so the effect of latency may be confounded by forecast date; we will also consider an interaction between minimum latency and forecast date. As there are only 16 locations, the location-level characteristics (population size, number of sites and population coverage) cannot be fully separated from each other or from the location effects. +Because epidemic phases occur in runs of consecutive weeks, the forecast date smooth for the wastewater-informed model, $g(d)$, may absorb some of the epidemic phase effect. We will therefore also fit the model with a less flexible forecast date smooth and compare the estimated epidemic phase effects. Latency was high for most locations at the same forecast dates (around the end of December and beginning of January). The forecast date smooths account for effects shared by all locations on a forecast date, so the effect of latency is estimated from differences in latency between locations on the same forecast date. We include only the minimum latency, rather than both the minimum and average latency, since the two are correlated. As there are only 16 locations, the location-level characteristics (population size, number of sites and population coverage) cannot be fully separated from each other or from the location effects. + +The relative WIS by epidemic phase is therefore estimated in two ways. The bootstrapped relative WIS by epidemic phase (Figure 3) is the primary, unadjusted result. The model-based estimate is adjusted for location, horizon, forecast date and the wastewater surveillance characteristics, and is exploratory. If the two agree, this strengthens the finding; if they differ, this suggests the unadjusted differences between phases are partly explained by the other covariates (for example, if increasing phases are concentrated in particular locations or periods). We will present the two estimates side by side. These exploratory analyses will be used to generate hypotheses into the impact of the characteristics of wastewater surveillance system on the performance of the `wwinference` model. @@ -191,6 +191,8 @@ These exploratory analyses will be used to generate hypotheses into the impact o Figure 1. Schematic of the model (A) and of the retrospective production of forecasts and the evaluation procedure (B). +Table 1. Summary of the wastewater surveillance system and the forecasts included for each state: number of wastewater treatment plants, population coverage, average sampling frequency, minimum latency (median and range across forecast dates), variability in wastewater data, and the number of location-forecast dates included, excluded because wastewater data were not available, and excluded because either model failed the convergence criteria. + ## Visual comparison of forecast performance Figure 2. Example forecasts for three location-forecast dates, one each in an increasing, decreasing and stable epidemic phase (columns). Row 1: forecasts of hospital admissions from the `ww + hosp` and `hosp only` models overlaid with the evaluation data, with the background shaded by the hospital admissions epidemic phase. Row 2: WIS of each model by horizon. Row 3: fit and forecast of wastewater concentrations at one wastewater treatment plant in the same location, overlaid with the evaluation data and shaded by that site's epidemic phase. @@ -203,7 +205,7 @@ Statistical significance test results using scoringutils/ mcs R package (same se Figure 3. Head-to-head comparison of hospital admissions forecasts with and without wastewater, across all locations and forecast dates (the version including the ARIMA baseline is in the supplement). Intervals are bootstrapped 95% intervals. -A. WIS by model, broken down by underprediction, overprediction and dispersion (horizontal bar chart). Bands indicated 95% intervals on overall WIS. +A. WIS by model, broken down by underprediction, overprediction and dispersion (horizontal bar chart). Bands indicate 95% intervals on overall WIS. B. Relative WIS by horizon, with bootstrapped intervals. @@ -211,19 +213,19 @@ C. Relative WIS by epidemic phase, with bootstrapped intervals. D. 50% and 90% interval coverage by horizon for each model. -E. Bias by epidemic phase for each model (dot plot centered around 1). +E. Bias by epidemic phase for each model (dot plot centered around 0). -F. Heatmap of relative WIS averaged across horizons for each forecast date and location (x-axis = forecast date, y-axis = location, color = relative WIS) +A heatmap of relative WIS for each forecast date and location is in the supplement. ## Evaluation of wastewater concentration forecasts Figure 4. Performance of the wastewater-informed model's nowcasts and forecasts of wastewater concentrations, across all wastewater treatment plants (nowcasts and forecasts shown separately). -A. Bias by wastewater epidemic phase (increasing, decreasing, stable) (do plot centered around 1). +A. Bias by wastewater epidemic phase (increasing, decreasing, stable) (dot plot centered around 0). B. Calibration: proportion of observations within each central prediction interval against its nominal coverage -C. 50% and 90% interval coverage by horizon, averaged across wastewater treatment plants and states. +C. 50% and 90% interval coverage by horizon (including nowcast horizons of 0 or less), averaged across wastewater treatment plants and states. D. 50% and 90% interval coverage over time, averaged across wastewater treatment plants and states. @@ -231,11 +233,11 @@ D. 50% and 90% interval coverage over time, averaged across wastewater treatment Figure 5. Exploratory model-based evaluation results. -A. Relative WIS for each epidemic phase, adjusted for the other covariates ($\exp(\beta_{ww} + c^{ww}_p)$), with 95% intervals +A. Relative WIS for each epidemic phase: the unadjusted estimate with bootstrapped 95% intervals (as in Figure 3) side by side with the estimate adjusted for the other covariates ($\exp(\beta_{ww} + c^{ww}_p)$) with 95% intervals from the model -B. Multiplicative effect on the relative WIS of the wastewater surveillance characteristics (number of sites, population coverage, average sampling frequency, average and minimum latency, variability in wastewater data) and state population size, on a log scale +B. Multiplicative effect on the relative WIS of the wastewater surveillance characteristics (number of sites, population coverage, average sampling frequency, minimum latency, variability in wastewater data) and state population size, on a log scale -Supplement: location effects, horizon and forecast date smooths, model diagnostics, and the sensitivity analysis with a less flexible forecast date smooth. +Supplement: scatterplots of the relative WIS against each covariate, location effects, horizon and forecast date smooths, model diagnostics, and the sensitivity analysis with a less flexible forecast date smooth. ## Extensions to still consider: diff --git a/docs/supplement.qmd b/docs/supplement.qmd index 095ad67..9f16403 100644 --- a/docs/supplement.qmd +++ b/docs/supplement.qmd @@ -14,28 +14,36 @@ Figure S2. Figure 3 with the baseline ARIMA model included: WIS by model broken Figure S3. WIS by model stratified by location and by forecast date, broken down by underprediction, overprediction and dispersion, including the baseline ARIMA model, with national hospital admissions over time. -# S2. Epidemic phase +# S2. Additional comparisons of hospital admissions forecasts + +Figure S4. Heatmap of relative WIS averaged across horizons for each forecast date and location (x-axis = forecast date, y-axis = location, color = relative WIS). + +# S3. Epidemic phase The epidemic phase of each forecast was defined from the change in the final evaluation data over the 28-day forecast period, for hospital admissions by location and for wastewater by site (see Methods). -Figure S4. Epidemic phase by forecast date for each location. Top panel: hospital admissions (7-day rolling sum) from the final evaluation dataset, with each forecast date to 7 days later shaded by the hospital admissions phase over the following 28 days. Bottom panels: log wastewater concentrations at each site, shaded by that site's phase. +Figure S5. Epidemic phase by forecast date for each location. Top panel: hospital admissions (7-day rolling sum) from the final evaluation dataset, with each forecast date to 7 days later shaded by the hospital admissions phase over the following 28 days. Bottom panels: log wastewater concentrations at each site, shaded by that site's phase. -# S3. MCMC convergence +# S4. MCMC convergence Table S1. Number and proportion of `wwinference` model runs flagged for each convergence diagnostic, by model. -Figure S5. Relative WIS overall, by horizon and by epidemic phase, excluding location-forecast dates where either model was flagged for convergence issues. +Figure S6. Relative WIS overall, by horizon and by epidemic phase, including all forecasts regardless of whether the models met the convergence criteria. + +# S5. Wastewater concentration forecasts by site and state + +Figure S7. 50% and 90% interval coverage of the wastewater nowcasts and forecasts for each wastewater treatment plant, averaged across forecast dates and horizons (dot plot grouped by state, with binomial 95% intervals and the number of observations). -# S4. Wastewater concentration forecasts by site +Figure S8. 50% and 90% interval coverage of the wastewater nowcasts and forecasts for each state, averaged across wastewater treatment plants, forecast dates and horizons. -Figure S6. WIS, bias and interval coverage of the wastewater nowcasts and forecasts for each wastewater treatment plant. +# S6. Model-based evaluation -# S5. Model-based evaluation +Figure S9. Scatterplots of the relative WIS against each covariate in the model (state population, number of sites, population coverage, average sampling frequency, minimum latency, variability in wastewater data and forecast date). Each point is a location-forecast date, with the relative WIS computed from the WIS averaged across horizons, shown on a log scale. Points are coloured by epidemic phase, with a smoothed trend. -Figure S7. Location effects on the relative WIS, and the horizon and forecast date smooths. +Figure S10. Location effects on the relative WIS, and the horizon and forecast date smooths. -Figure S8. Model diagnostics (residuals, QQ plot, observed against fitted). +Figure S11. Model diagnostics (residuals, QQ plot, observed against fitted). -Figure S9. Relative WIS by epidemic phase from the model with a less flexible forecast date smooth. +Figure S12. Relative WIS by epidemic phase from the sensitivity analysis with a less flexible forecast date smooth. # References From f03649aff8571ac9e896cab18a752a634f1c05bd Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson Date: Wed, 21 Oct 2026 14:55:11 +0100 Subject: [PATCH 11/11] Stop tracking the rendered analysis plan docx The docx is generated from analysis_plan.qmd; the local copy is kept and docs/ is already gitignored. Co-Authored-By: Claude Opus 5.5 --- docs/analysis_plan.docx | Bin 19077 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 docs/analysis_plan.docx diff --git a/docs/analysis_plan.docx b/docs/analysis_plan.docx deleted file mode 100644 index 37ab5002a10944a1b8c6b427e0334ccd3d950576..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 19077 zcmZ6yV{|4Bw=^8vwryi#+qP}nnb@{%Ua>Q=ZF6F4l6-SN&v(vx@AId7_39tJ)~@RA zUA1c~%7B8Q0f9h40r?17=|EX_&V&F10mVWB0igl`0qKZ1*t?q9yBer^IhwiXF?ib9 zwx&+X4>2K$1&~i*klSk!g}!S7=?RlhK(;@Z#76L?@tq&uWRkH7LnWmFvL$Dxa(w(u zXp7Hh(B>h0$L8Y>jebOsHo&;qYaG7%G0R*Zo}0F;@jxNBZ7{{Hie3vv6s#DUO3nuz zONZ@99+VnDQ!a<=2DRhEU55(G7*ZETAuEV0B#bu|oOzAb<7yGa7mt?JC*h1PbdYEy zwbT8sdUcx0eVSv-3lketOF#rij6nvfQ1&+BP5Tx zWX0=v*sqQ_OGbflfl0-7-lp;AF)otRJMoo#sX-gak9Zhi!#wH*hgL$X29vYp3&E?0}n3o3{dJko1m(zQX%qd z5ugg(uA@6qw2$ACdyb>k=ZoZ%f|fHo&LQ0DJ=_;9jC(kgJtxl>ZK2UTmbBt0fZYgT z6*mLoDA2bI*a*|ax;s8?*9+mL@6~Zi#eZSe=23EU`WLGO5FjAfzh?tyGg}u%hJTK= zNpkX_%y8oMC|_dal@~4-d)5Z_YPu zzcXE&m@P9n=&2AWv@yNmqVb*oanGH@fTAwl|w%JhU@amb)Q;z=?>*pE(Vj})TK(LlZ&ZeJN(`lO=O`I<^-7D z3Nbi(xh_)X8f5PQG5jnB@}7=4(VYLT2i|5mVjdw7P;V6w5c1!8csMwlGMYM=xY_+B z@;?l|(z|fp6i@R{W%yo8q9G%fUi!7+qzp)nW^p>ul+VeXEXxv1zQhuAnJcg~w+7rJ zJlIsZ8yFLs5|Su1o#ylw00*N{!{T{02cSfQeRyA=@<8@|9O-2>;3nEsHc`gE_Y6tM znkGh>arob#b*<*R%@X?E-`$*GHb$9Oj=@08TD%*79nF_em0?0c{>X|*xyZ|kg6qYt zkLzN}-8qz^fs{|n8ZJtTT+}WFgGnG zqjH-z@>9hsdc}=78q*Gc)HO>IxsBXFRR&bx`{bi~(*D~?W+{HBdf;B zYHt>tMM;MZb@m=ni5wrF0KW4!KRInAJV2KG3_i8Z?Us`zzf~+++fh#i7jh;|%W6(H z7-4fTZ(!(#AkQrx5%(0E=85C@ZSx}d6kL}0kTO25-!|XCl=o-QBi+w>Ot$?_H9ul5 z?YB>Qb$7e49p{hyaa~MCubENqe0e~|^^>Ue34Xlh8CAMSAvTi0A+`brPH#R1ugOni zb!wIwXQeSw46#c6#$>Rj5xH*z9Pj*E%4GKU+$z^ zCjpPMC^6+dVxg@1%eeFy-;V3>WAHs^^0CEI9xowfYyKSY%#n^|C^vV>B(tSj!0 z0GY947}nvU$@_JEzbC>E(Feg0WG^~+%*tM3bvUAQ2eN+_q<`0TW^T2{&&v{)>CB#O z&Ffwbl!3g=2EYuULsE?qg<+lA)tOAqAFPtemU8O>|5O+zs9 z*2KYSKpEZ)-7Qdz5Hck?UPKB0xl!t<376aMqdHGS6*CfC(G{3Fob}0^^wY5Kv}!0~ z4L!lD?t9C`A&WDf)@)V5K+hw;`ld^wC6V^9{A0Yzp+%{tjDw8|KJ?-j84Cn{!k+il$#}6!P2vJ(8l_9W0Xj$P)ZVfM#O?tm~ zJvGydam72{>joQ7@4>t{XD$p{=w5#pKlpESmuJWURNeeHT{I=cxC2=B(d?@|HJRMc z_ep`hiG5!eqR_ViUlL0-@N-gCc>fnxv`0@vZn#Z)Q*CPsvLbHO8%6g1WU{AoLZ*yh zA14@CcY?`z$s@_+#0i+YxPzoJ(gI_TTo(Ie#pjoI)VMgqZw&Txz5)F_MsJh@&e;XxG$ zZ;3zqMDn9Ncyn=nFL3uHUt3U{`VA$sYzA2-y`?{FL-ymA10#1II;|riO`vfN`n_3y?=VobOMAT zzBAVeo-R3(afeI71g!9HMaLIyG*3MI3p6s68v^-;rQ~O5zXAQezKfJic(s{JpG3rc z<%Lt!ye`%T#JVVrSu^oT8~U6Zd^)le8Y6}A;qYS6P7J*wYlR=8J97F-t)RW2a-hY! zhYIWsP|C(!Gc=Hs8JO%z_`v1KQZP8q5Pi^@kTGN?-o=ahf!^K1F0h)NhUAe6+AJgD zSMP64g<#jDmgbJX%xDrzETKEcEK&x46T4|f*FHWQdcJVMNkh(cU1tt4aiq7It=ev5 ziS08dF2g#*P>je~V>(n!puRFO)ft7;u%En}TZrQVUc{4MX#Sw%5pb`}jDewYm~cet z2}n=MAnKL*!u$?Arx*~%xONBBqW8x(;ujhI80)M7@j-Cah|?CX$P~Ar4L1^%-zN;% zZ=13kLXM1*2A!BMbGsg^4YWbIDTAhKEMwhcR4sKM5>%l`^lSd1u2yGh;&Wk{+vBFY zp`%;Y2_JGn9n;k%SZrwWlbPJ3eIydq%yZ)6Ca#!?GU&2&mw*fAW*utX2kOI*sWTeu z>|h0VTRH)6TPnd<=6QDv%p}Bf{5eRU(Mk!D_@ZBcNm*$zvMt+&XN}#nh)&%DO-l%C zPhCp|zYf+F+@RC(l`=aj|IqEgLNSUk<#W5^Y@~_@xx-QfrKb=mr+;yX+MPV5j-@=b z3`;O0+S-^@iz}PaT;g;)gApsg|D;NtupN^{oV1Bv56o|5$v7&jN; zn@Ic|_w*!m!Lds)G$$YplHCp6*R4rB89zTYYw84iV+^R5NyH)N4-Z08z*x4O<+2Qx zN4X$KJQh(h0ca5XdL>Kf-R~ET%bvaDei9PM^RKjH*NQhaU)Oibo~rl3PjgnjR*uGA zyuu2?HiLaW8-NdgOY;X`3sb6gyo9|liPe{ma@gKP=4-#ehl++ z=3OY*l&Ay80=pCpEWrq_g*Ndt@kTBDI*}GWG?Tc*5N2uliJc-Pq4?RdMquM5VIK-@ z3KyH`$2FH!gaO~k+Ru!@hFDB{S)Qp{sDwD=!-BwNyyAExZ`<{SePrW_zg-tPhJr~g z^^FMlcd(Cwcsqd{BN zULU}-GlUv5_pUNe69~Li_%qSP?x3_)#lVQx9g9K)^YZ7!oI&N$b15K27Zi~_3J_&< z?i-zIURsE4;_eGSLQopz-;mS7O>)w$8fCAv?X0wIhjx62F>w)%{`_8-a$5Cvj^PXp zz}5%B+9c2s!^m(#3__Y2Sh(BwQ@<jkeFN|16|f$N@1gXVNI_pFh2q*gsPxQPpQ9@;iZzS)Qf>!5+y zJRN5qe3+3%B;xEK^-#oh=g_iCo9=zB#27^@K5-__uW}0!__nB-Ar_%I7xPcab93WY zB{Qbmz#Mn=@z@91bazmGrm^`Yvr#ybP~()nm5f5rWWZ1uD`ruai>7XY@bs~jC`aw{ z)wx%TFmFsGTV7jAB<^w!m>2XC7>l6H;K-21FLX(eVDv=ff{K3?*Xaotsj@hy zTp;Q|LO;4@XK-DI*G8Tx;`dspq~7;|7DJ@B2jxr=neMK1t39ei5gD#R0alUGXv8fI zg(QI`U8h#;PCYhG%54s64A_cMtrXT0BQYz$uo(^%2Rp6MUew_L59@M`i^Q>a<_dA| zT?nk`Ag;w+p?iJg+iJ3+z2z)TZ6>4cc--q(_ZEoB9JPtMXh#PYHWs3AazAjfZ{Di) zh9Rik^DY4lRH95duvd|qji)k@ zUXUI;N2WK&(fno|dRn_O&(zYh?6m0{ zwPPcYhR*h5+hB{I-mt}mUQhu)e@UR?z>~JyI#iXhs?^LfopjIIUWNZ{U+Rn5zB?$y z0trN@e>{$W7V<@A?3(?F&EBH~-{CrkDI(0Vgy(ae%ofxzfr2S(0B7FUx&ku&^rj97 zO|l<0lQ%-OvwpQHOk_H@DVf|>k(>mSK{K(32FSvd;UqFW8`vcs3o1aiuMfU*D>AvdT@Dr)J>Q>Qhv9-_vh|>U(1z)m-q5T z;Ul~1d4NNV0VRHN$fxm)j?P!7($)hr9FtC~aipid_sxCx)dWKnvmV3K2Z9qaU%pWi zZ1E?nBkWqNbl~WV)az|%!F@$mYz#K6nRDckd#8w2D?$G0bG6Xu=n{6jLQO|tok2zz zNp{lj09i=5d&;W#0$GRr=K7|E>s(TI;C2>lBS{CclQ`5V2~ygyke)2AO+~U$TDBPI zbW=6yE2JOOiV`SL&lwyYXt*m{`ZA-P*e%*^C>N(8g0WD}=O%Q^JhW{<0!<9%JvcTH zp7n?W=R}hHP~V8Y4hhN0>kag)4%L2QsK`B~4Z~j}{0oPvT8oF?cc!sM#1phA^_+QXKC5mU zDHL&?GggD;Di01>0cV4;6@fn;B>|eaqDDwSi z=mF3*AbU0mHapNjqR&|U%rTDV-MJRX(B zkMuYkF?wu3{krW38f5BArfiQZwy6XKg-@HKl185R)5 zCWpqr-^6Kb0~A+2{#dWx)#fDL)LyAa6H7CNto zghJ-p4H*Z=tpM<(?4uV|Lh&KabP9E~KXvM8SZnkG%o2^ptE>V+P(3b)k~vVW6+T^M zGSc**cUNt0xyftUM7rBe%4kFn{{Gjm-wz9dK}n^JiFxLa6^SU*n3%GA-DwlFhw*oN zqpIdLKnmP(bNq^vkcF+0aJ;+mrOTjK%o@c271Us^cQ{Elb@y$y3b6Y&gO?6MsB7#> z%`}^*5e?G%?zgWH4~}K$KzJ%lwPO5NKA~8p&yLoh;coYpdAqEER-z$^d1VOZx4;R3 zX;f*l<2{%;O7zF4w5Y2h6L`j344z_)eL4@}=$_b5gjGOOQ1M#w@`aQxzX?9fmS$@4 z!F?_layuu}_wK4hV=>|PnjFKj6nCM4aSH?(d})<;R4JEaTtn)AA1k6u`D#laDDwdd zY51Rtzga=<*HXxEDT|;J<;8!Yplt(08U?}cu{^0sWD#(cz=a@5R4Oz_)N0oqY=mu9 z2EaF6(1|W9bZRM?;f1EFA_H7xy_F*iIi}Ocep4$UAiQX3MfSpjtyy~ZsN&F)HVy`| zgtn}jl@L2>a5jwH!lA8zNMDR=?6K{{RcS!4-7iwO#CNx}cf7OZlMkXxZb~l|NTOKy zATVT3nXe%pq0jCN3Ni|qtq;ys(VCO7~f6y!u70fg&j&I)7@|`z*+xT zba!T{c-CE(E}&Z=@4+Pc@-@0<$U@Sd;M|_j#RLOONA;czZqTCxYF=Bp3TK>*S^@#C z%xWSo1`m40@W;&kq5P|aCZc092T=GMDA&9?K-%y4K(~(90yCVoFwTkrZ^1L!#u%qm z+x-|Ll>pVDuI_4Yk{ztjd50za_EvYlRs@WveM2K{L=RII&HKbZ>m9@M3t(#;BLFis zTef1hiKpa-_bBldl%F$7d3X5Gyh!Tw)h>XO6H5^PI==VnI0gwiFCf4Z2O*(^+lx#7 zVpl0oJ;z{+YsBrI3cZ%YkbdSD{`5VcHNs+(VXUzN!;w|60~|#(Ug!vOq$9_c%0i2z zJUk+_wEbAiLZ9Hf&7O-A>9g4pMTbZC2%Rv2;CVsT@9dHZIdr;ty45B@y&Ad;y8fYe z--f6($URLk0YDKE`D6@>oCB1GNysEOXOFd&*~8FtGNt;(YOVvD*j1No8nL{rI5;?# z%d%iIBSxRP3wzkka)@Kc_J?f*oUN*cwXeFxZWSjbK=x$ z%`9=~6Kg&W+3X?;bLdXpotPnORmw7%35$+Zr=4#;BL-VX60;fm3?LH z%Ez>Mf=NPt+nSz1!8G01c#K9qD-1(<9L0UHCq~b3AvX>J<8d`%%i*7BH?RxTRZ0P* zeq>vHyo$~IQX%PY zoEi5o4$3Ay+}taD56dy?8VqIl^8-6`T?lD%kI+JJt3I>TM}L_Oc&KS{x@)KHWFF(0 zir^LE5nH^GuvKb-g6R*|h8errqfp&a=dvb-CdAp7$r|3C9&=l~nFYvNFqP}2rd5&E zDu2eXzJ+bHigzGiW1T`I8^ikY`3NSF(8Y(>Rr1UDWp&z|=+qmaVjnHL5Jt|dmm5*2 zfb$yOiBJunA=_Qbcm1Wzu8yvkbDfi*@3)(kOis+=VbWpGKZ`LypuKPB5cR5B^C@V! zBG@)FkLe<6W+>4=y)ZNFLaQdpD#E$3efr&qSb5h=zkafFJ%@x`!aulX1DxtGmz+&G zEt77PeR|8&ZBS{0&m4~VAmrek9US6hSoZcI!oDM+3n1vkd%kUNm(PFK&#Vw}K?o(5 z%AvK)x2~MzyC1LkQH-vE7Veq=iddU8$4~N}UTD&;wj}u6uguLFhv_}Ravadr?S}St zFSX=;@zA4@?N=UP)Aa3^G(d{D)_Jjo+YLwhxfp@7|qI3)gDMUywc5TCO((rx{9|)nK{+t9`$$juO z()W@gI`se_cz}yy_{H(lY-UNsTMjsYsgllJJY$q^#TOb7S=2reYB@4?0MzF<8e&lDP{8LDYx`!)ex@`G+NYV%>I_ja#U`(e(z_#c2)`v?6Qup z9>)r4KT!!JLdA`Jf}1-To3AMGq5XLnrYb$M4Wv*}dc_cMg ztw4JS&ztLKhiWzU>eWU>7NX&3{?pp(R`pw06g(PT#JN=fb^f7EwpZdk^#b>fYT-wsjzS@7~&jP7(l4X9$6G9lGdrcVR{-a&R# zBEUy{TvU{NEliGfn##KjO-#7*^w&xxwV)5BMcs|pE%$KS13-Wb)1_0)T^WU@#-Ubl z@qN3hRBhLVVNz5)M#DCnrvV0Yr6&>t#kdFqsIF0H_SRF8FRWXOHllEkU8xZ|= zqzAh5XO}`UJ7rmS5k$z!J))RWmcLBpPu5()EY)I!!fkZ^TDNOtIw_6`ki|cheybQn z5G>*n{=PcoUiLk#HieGIiN4>}?6!X#X7X&NcSm}c5)DnVtn=qQcF|ro&*}m{xF1ZWW&7)tDmF$!dlDPaG)ILA(szA0j zPf>%f``czk<5s>1Z&O*ihxjpRASvH{m`^Tu9UElwKst$v7&zS*AgJ4feCbB{kTW$ojG%Ax?*FM-inM(oMp_9U!ETD z-FX^av%6#m+;zEerlDPI%)rT4+aCQmHZ0Boyq?86Y>VPOJdqn82&D%LL@~yRrj3zd z6}`spQlTKpJ8-ZU#u>>`6Y)yfA8h|5XcL{>GJ&SnB@CN3 zrh7n)BjV4AExZ2d0MJf3wU7QKiO12?T`kIw%Y(X=0fHQj6*5^-X4rmH4vjL|wU=o% z&q0QnS$iRw5&_>*9Kt!2LsWHQj(n~cE%4UIM7w@PWlas45(SNWadRFmko@7hm@Xsu z=Wg?c4u?HZLz!L1fIn7CIu%~NQX3*81oQa$L528Ngl2H00inD`eQ zUx|#z&{5OJPvFFjOYp^@J!;%d$H2$U{1`+&6Qp>0>Y9OAS$5aB_?}|eQFWr)EAVFh zTO{5Q<3Gg>mUIZ_F8sVv<&L&u<`*VQr1IXwH@3DSVyQXm=BtROhBR?Jo*N`>zYFss zNSSkCzZL;AlKAXeRV2u;%zRho?;|IsUELyFU-$V}i+xw6gH>{*KNNd)>1Dc3!$m4& zf#d2tp==I*zZm9E8MN@ocyNUAWUH0DR5Xb7N+K@!c(-Y{L6)|b!~vqSCmEwJJmwcJ zujuLfzBKEFopAgCy-HNA_b|FP(o8yL09wn$lqIVf1It7KpX%8F;D2MCDT&S`_!sMT zP#_@ue`9U#;NWWS;A-abPk!USU?(dn3^E~&_Dw9j)Y1|~7$FTRN|8~Rze{&s19Dxj z7Tg1OcdD2ZEhBod2(ouRuU}?clh+Fq2M}V!YEcG@GQbK#85Gr^y{77hWLUE02?U)_ zX>DU`;eDfZw|b!ooHs+Y0{3Rmu(J3S>nAbZx!BTQEmOtb19;m|4g4wO`?t z(#c^YThDM5uOXRg`Nkzxbx>-eLNnBu+lhomR; zMk9X3JdK9T*0nNpZ{F|i)rNvS3T%JuZqvHDrAr8Nx+XHv^b~*iMO1+uazYv%x*c=3 zzlc+A$IcA=H6O*Uu4+YAxqr;QY-xYv@ic>`Bz=fxzj=Fc_Hjo`4)sZBZwb^cB*YHs z*y*_H7Kt(CRtefE+9{vqv^`3}dXtdwy$3I9{JZogc7(;)fAd3&f9sC&Kc%;Kvokhx zwz9YQPiQ97z^0lBYospNs!Mx_HRR0*76aVb+-ynk+{ zFyE<$tPypHW@{R!t4ePF=-tuCLx^0~lQ7&2uwr>%Ju+FLPTBoULSZG(@N-@2$`pNn zz5B*gQGi8XfW=>crJXRdfH1RxFtdU%^EJS%*Q4HI+h&_b`{jM5XQijujS+C&dW~P} z1A3uCshvPB&Vq(JWs+zGVLqLypIq_!DBkguCx6kU?^A3?&I&$1Xb@f`h2H zM}PwtX`)LY;PP?>m5g9;My=~3eBNGYbQr^#pl-uJBbqb$IsaytukYzB8omdE_R7;? zT+K(O@&mAVH$}e%e1mu?96rvGaS3E5bw4+d4bQR<&+)=@HYY#JB#sYaC%B&6P4C>} zzP90#xRS@*s2ePI=2c_E!(-_HnXu^pov$AKG>xzZKXc$Cr0_n&^Xa0(z$OElYJ9y8tT7dRL%(w_aMa5%!+z!-$mU0lK!?jyIaKBcX;cf02n zJ3wFk7Lxvtpnf6xQbc|S1l)0M@-{UtclI##&o{1By z9j1m26Psoq(aQtLo!M+8@G~0aC$&d7f?39MzkV1!dFgkCl7j4s2(i+-WF9!{so-tA zGs4ljmJ3n3AJhXv*?_1IeyRAMxW027t2a`z!M&bHRSnsKp+vgK@>!#->XMc&DJBZz zhH+)vTV^#2)^_q_*&-({V3>g6TBnqjAEDF5rw$ZHKDEdGa;)uS=C5a3~kRnxSZ(y0gW@o*k(iX;y43vr@p(Sh7;xjf+g}E7X0Xh0_Qp zKj=Yg0U$p={ZaN6AXlrb`C`RU+92YS!PugFhXUU1y}}r2g{Nuhq$F2kv5n|VA*_*) zzYt|quN~Zc3{wjrb3(eIHYf>*64@>4=Dz@uK}Q2l=2_yR)3y)#mDR2v(LSaLH&Se$ zR*Z*-ofqOf#Cjy7kUiMbNxijL(m?56nCzykGlMDhYR|vl)km}c(PKMEX-W>YyGr5P zXFvE7fT}-?(-~TC0Vz(||oI?ta9lIP}`PLP_ zmWTI@Nw%135rZ`>puE?_D~Gn*b!!=!$`HyFW$u*7KYga9rppQajIS^u>J)f_dXc`H z6nnEM7_!2#)8Cin0R^;(w4<1lI!*Bm-H`R`SLm9Qj1h_gm6s-(EXnM>e;i*``FIO# zc*MZcdY|NQUDw|C9RPfUm!HkfE|!}CL64;(QQJA411V#kMO3To6wc%9Wwz>!3uOvz zTjCv`N4?=*k&7^C_n?tJBsh zXDORm>zf!09sKIB())maKnb4L2M zo`y-%yX&$n(x$D@aIT8lto*2`-Q%mWMs`Mt15iZMFnd7cfofF47IlXqNBc^UJJ?^A zo{SWg370dc3Lu|ZPt-3>tfm(YJIQNwNqBXEB5r8@k>1dh_8UGiQiw9RiHXv{g&XY^ zKm`g;GanW0QMboLgk7@wdVHnDnu?sI%yrgdsC|<*e1jOrQEXWpL{XavNjk2G%Wb{C zLF6U?1X0GxBx_+?=#>8$&=pR0xqk3(3&YtJhq7)hvi<=>-QD{E!kHjA8v3(~e8nTO z3f=Ar3|j7&?9PsvP^wgdjGHXdjT19bs7mK6gobW0Nk&d``d)0*uM!_W6_}0R_5|12 z==nF-C-h2`sJXv$u|IPQzADrAlnfJJf1}C3l`D~4psI%7n8?9~QCRrNnBYIDYa7=c znz%_RUzgMi;4qYVQTjoDjf+k?d3q|yu3$&k-k8;&L~HlP@O^3;63NM&OZ~!hA#M78 z4d?7@wihG6Ge^a^1sGebwHPsB+a*jVWnY)|-6?pGat!w)I*!g9a$SE9V(z4)`}6ka zB2gRB&LskE|Bx)@hc&|bd=kVY}69lLWp$3tLc^_Kfk3ytYBI#pz)tPZ|w$eR=Rlt}_Ch$O(SJOTu* zh3>2Oy}Q1Gm13A(o0`X(hArhZ#L0B6D+%=LcVb2Y`*N}VO`@DbstNPohgZMdh6s;N z-m+60sUYv5Vtru~vrCwHsc>?br__?K_*6YsQ5&I$7Vmq`@hq;4V7~UZzhV-_|MoIS z<4kRAlfU(7U=c4Rv#&DvM2JXk_eE~!VJKVngs3)C0C zzz- zqb1{}^|sTiHz>&Hpjl6w(pSBXdWS3{A|*vh{2^crgT(4n=o7&F^`tOiJf1(WBphOK zG5Us5PeG*$N2RrCejsXy2nk|O1A0egkuOCo>sk338h}VN!X`YD zh>iMkx4NKtR!DwC#_|xrb$Is&(u7^dIJ_7ALywf^0^WUQSKI0Hv$F5ko2%>uQQ%Wz z_12A6sX81p1UY&3K=27Z#Dz<(1PJ~^syz#_ThN^5$V z*7;X(@Ngkf3pz2+?~K!v0dNK!fH6ruMiqnw6!wSaDSSnE6COEpfbAbK1V>z(T+Sn)kWz+ZnUx&Y*lA^PN zqYI;ngR|K`vbR)uK>m;k=`&U7y7CZ&He$mu(9EM$j~>NLXp7EgLz@g>rhqfJ$|_~G zCcV`zn&st?x7$5+a9>(!6wa)5la`r-j$N?0quHiq{!*R&hqtQ+<44MdM0YH|Az3uE z=Txb3xEX*PrW8rap;5_d1a5S1V^kbD-H9}KVb(gtSRZXyG^8YkgdAdkS2k+Bt7Jae z)Y60x28+I71OlI&*h1OY5!xh~A0I36J0rG~f3JKAl&jGoDK|oIb-MF!1pKSi^UuSF zy)JtjL!vWIv$c_{RGM{7?~Up<&B0?J0BWW7BJDx+)`r3Jc<+$lT9xDyZ4>M!@KAUf ze~)eC%(Dy(d`m=fDT_cAs#Og1X#e+3pTTNJcVIq0xJ(tGp6<8a0p}hh!XaM<2Z8NW z^qbX8K5N7n=w<&Xfjb8P8l~!JTUeGBZt^Swv>azsWkBOL#_zOPDh5yOCct>sEIG%p zWA-DEH%C{m-oDAYAyFn16SMlDJEe)2#}6luJL?k%Uz06|q!=g|8kYjOe1 zaKQp*?H9v{xDpqyr)DnAQF_15?Cublv2WCmo!A*+9VcuKH`F?69pV`#TQly;&p?N5 ziO{JWt<>DP!^usg+~y`#iLVl@d@s#LH~>?eu}Y4I%Q}b*&A)v8XZA=#77OABEj1D1 zCSU#-s9d3Mr{`VF2KpPTqy70Mkvn`0Z6PAIR+fy+W-Uf)2Qt$OTy&#|U+YwV-kDqn z?uYduH}~4(lrVf3++0Xpf%!pk6?>a-d?zD*(dQzKUiEP#cqvJf#%onLtvu)zP46H66q_ybqb zVg5)VNT&R1X28PAW-)|IVt?)OKz(75-9Ll*8u`_$j;n_&o1mjOzrm&$4&PsFhGJU~ zkixgGI@N&nO?mf!99yzs0cYpGh~Pl{cZKkO5ixOdadoi!2b7+afqzyAMFYt1w_}r0 zHdLjMTOOwh9X)d7OvW2YeN-xtzT>KsglU|BOF4D)r?d6yDY|@YR2v{KO6IvEDPq71 zeo~k293Odg6}mXLP36^C4HSW4zItCdE0Y4EnOy%mcU_SID_^B7fc>+2cg(D+B#;$x zS9nV}lM#E7d2Tj%yIJR>lD4%&19BPCxa+OS=Q-~lvuARDs457v%;XohChdCEw7FMi zz4-_+M7G>A)`^ptUx;_s$vXo62?cx@qhxd4K{>YwHYV}jKDrgp%Y}O8`?D^yZZNw~AQ<`$Pm3E%Bz>9_}h zFq^pwZN9nF1cd%SCXK75nVlKqe~!%mG5oaU95y+ydS5hwfkQ! zsCkOLQV9wUiui-Cz#BxINjHx1al^ft;S?&)z*EB(B*fdQ;xT#BS%E6I$Su-CAcaFf z84+pSAcU`p6kw|0tTF&Io!%h5QKIn@$cXa|B6ORcvMppxQ_@11*nFk%K2SR#L&OnP z;yQl=^IaE83<7~U6I1>U%lw1&3tU2pI1x>>rdW2~EtxQtK`cSivf9;l$>25>PB+)eN)pe?+H9HhVk;%j zY8t88)-Jy~nomi$%volGU$kG5_F|hsiS5cN+pnZra+rWcp#W21%jAx)&x;m=+%n-paXU=gnv<;cdw72uU(}#f z2oUgkzC1dwH^hwhEhJqF+NYAcy)Jk7^135`{b(cH<@F|im|+tR4W0pVssrCOQxPlOR4lXYvA9aKDiF@-H=#GhRL%-(Id4B6 zB3iikDIx&hIX=pwOte-m%t)wa!sVX^2|M&e@B~EaGy~^mK(Xa(5o19ii(4I0#ce7b zhGf-d(;!IA1t0d5x0Vx=qPQs*(Qm>f2({H$3BC|Wm+G(b1QWpmUu3Hs{) zarvQmcl?@J51ma7JDZC@_IfMQfA+^v3wU?8h)+_n*X5ZAqHE)Tv@zM&7Vj3?zlaj6 zDw4t{$7!-_sb>VA9)a`9vaKjNRc=3`pM+6)&|uaprZ#n_FO?f0rC1hUj*1E z+mdL)GjOfI#Q!Z8-YJIvd!xp6KK-PK@qHn>_p$j}W|_d%-GHYrny)Jler00nS3xGSI76#;8OqNVr5|`z?`c(V~Sp12E7YyM>R_H6!_Gc zS`4~QWp(n{%kFu#ipX+XOJr~AAGOofRThEF!B#o6I&H$bTU)Zj{l=~PA<~#Wd$$Qz z+A&w6-wSb!I+UNCSCREu6xL{q zYuQ%1mG)-7hbGv6f?rgPN{hf+ep<4}T+!Q`O++8Lq!h z)pp^k?WDGH3Q5F7a~yra<4k0v5 zn;@`&%&TmFdsFj*B^Q9Fge8vj_1uuR-l)djc4pA6NOCL2yVA&DX07oS30V5{B*fZ@ zb~;y8$raJZ%%E3FJw|CQy0(3z%t_hFf{>lxUa#R@y`Ov-f#0UITUIuouPDqYZ7I~^ zxc&g#Jt#{jqC;C8y4&Nvw{(N0J~{~fwkIRd^loHh0Q(fleBQ8xq7=P7IH&R^!?K+= z@XJmwS|45~7$?nO(&Z$`mM`9F%2gkz`H?H!p~?ZPR<1H+dt0Jn%484s~uR1O}L3)ZOeJ7Ubqa{ z5?;6O`e|LdMAc&3X2)q->CzQT3Be_ep)>8?wdr;lyqR<{VLn97mQ|`gQ0t!MzWO_5 zvrJc)m$jv3aBJA*gKo<2YN7An?I&VBCJ6oC){_eS|Itw{X0EP(r)K}vQL2CcqD*i- zHyU94jE&hPTF7iyg?cL@JZO%cT_lT*isNWhzw$mVBm$^*<_clmo1R|BDESw9dq1rn zOF*1qCgoh2Ws0CvP*~LnIOZqq|G!!9?aIqv_j)Vyd&y$qXVdd?%>p0pjXl0`GONVH zw2XyMBKlHVo-4~Hohep*DEIWD!^_8&lh!2hiCqj}F~2zR)YpK-tSKrbGk&jh+5B#D zyZea+|EH>6I8@xdF5p^gmf!r~b^NxgY&o->gjeynuIit)qK9wKFPXJcnwwjM<|J?K ztXeXsvy{Ljf5e!72V3UK}J*5f|7ca!L0-ggl#0qjxFkI71T+q`|SWW|G-3xWdd zwO{tL*q)o8{w&YVSM7SCg~6*z<^{6%7EAnYF8Vd8O*+HtwrQh!V($CK32&8Rauz*Q zwMzWlaMnFXaKq6KUFK~v8&5aH_TA7*SS@y2J|cLv^F{A>{GT7*Sow5?>9q;ACVN)y zNmxJoh3u72U50Ug`nL0?HpT6}m8o|+Vll)Pg5toeYm;4Z)C94hk9(GckRB zX~_XbpG_*$a`K;S73w?hWZ-WOC*K@6#O)3BQiK@O>nqzBN&1dIz(>#y4}VGHOHa z-B5V#)A4>sA>a1wa=};Oy*CnHo#je@b+hf;Y1m? z$twv-Ov(Xwwck!TowwLPr0shp*XiGyZ4;g?ny|{e%(y%7Hv`{{gp=Y=9=AMud)z|B zdFr$jy26H9oM_tDwTwAW7@8W%z6sd;6vKX%dY~wU#?+ zPkUg?T(en9_qFT7sjN@z!Z#Rj>aATX=Ex9vQl(Wp*6`^bCuz^M)uyu@#97bZXg;@X zlBTpvOALD(XI5;cMkVWx*tP%8WNkI8Of5Oa5VJl==Gx>b+IyrFub6ebYF%{cYi)ES z@4b(4u2b2>tg78zn>=0KaVT!AmU}P3zWipSC$B;EEMKcQM^*9Pi&;(CU+fUw9$9Hp zex+cS*qX=8_btk~vnM5dXKr!XtZAo_7_Q(+@8<6~+YV|7l2?Xs8g3E%x@o>MDN zKk!VR*_)m5^G+vz87;@}|37~fTvx4ieVyl*$+tH3hY7E8-Z5(e@6yuD9UfPB!{0Et zef0dv7gd#QVSKjkO|xlf$1}SS;r=`cot4FD_m+l!5TEw_?cv*>YbEBhzdrpfQF=|i zUU2WHy)zdzSqIfF2%4{4?-9%JYSQ9YZ(UqJemU*;?J-B#V*3DZMqm+%`;Zx6+C~Cc zj+{X^0R6-hgzhcC6oz#E3A#q~qd^cFj{seWaZm`lcJ#$a2<;PHg&gSm(P!)t`W-c(`VsT@=$g?-V-cDc89_Cp z_yAo$`Zy#){|jTNe&kU}bnWQFcL?ntOpwNcV1s$++R+D(5ZX7IA!&yWC86s^ALu~n zcD6v$jTrVoHvqjGjxgYvCDiliopN*|(AysfBMcm2k%v@Wple63Vi4MUoT1v$DjakZ n(2I412~jRk6ViZ3NWlyF0B=^{xDaTpPJn@jg@Ga24a5Tg