diff --git a/docs/analysis_plan.docx b/docs/analysis_plan.docx deleted file mode 100644 index 37ab500..0000000 Binary files a/docs/analysis_plan.docx and /dev/null differ diff --git a/docs/analysis_plan.qmd b/docs/analysis_plan.qmd index 756c939..ace4c4a 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 @@ -23,41 +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 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 effect of incorporating wastewater data on the performance of hospital admissions forecasts from a semi-mechanistic renewal model in Germany. -2. How do characteristics of the wastewater surveillance system impact the relative forecast performance of the wastewater-informed model? +2. Quantify the model's ability to forecast wastewater concentrations. + +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 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 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 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. 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 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 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,116 +69,175 @@ 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. -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 - ## 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. + +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 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 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 + +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, +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. + +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 + +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. 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 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. -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. +We will present the results by analysing the overall scores by model and then stratifying by location, forecast date, and horizon. -When comparing forecast performance with and without wastewater, we will compute the relative WIS compared to the hospital admissions only `wwinference` model. +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. -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 +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. -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. +### Evaluation of wastewater concentration forecasts -We will include the following variables with which we would like to understand their impact on relative forecast performance: +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. -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: +### Model-based evaluation: investigation into drivers of forecast performance and relative forecast performance with wastewater -`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)` +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. -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. +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: -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. +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. 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) -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. +Resulting in the following model formulation: -# Results +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. -Figure 1. Model schematic? +$$ WIS_{h,d,l,m} \sim Gamma(\mu_{h,d,l,m}, \sigma)$$ -Figure 2. Schematic diagram describing the retrospective production of forecasts and the evaluation procedure +$$ \eta_{h,d,l,m} = \log(\mu_{h,d,l,m}) $$ -## Visual comparison of forecast performance +Where the $\eta_{h,d,l,m}$ is modelled as: -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 +$$ \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) $$ -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 $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. -3 rows, 4 columns +``` 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(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 +) +``` -row = location +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. -column 1 = hospital admissions fit with `ww+hosp` vs `hosp_only` +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. -columns 2:4 = fit to wastewater data by site +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. -Figure 4. Visual comparison of forecast performance across forecast dates for a subset of locations for horizons 1-7 +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. -2 rows, 3 columns +# Results -row 1 = forecast eval overlay `ww+hosp` vs `hosp only` +Figure 1. Schematic of the model (A) and of the retrospective production of forecasts and the evaluation procedure (B). -row 2 = score underlay (bar chart) +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. -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). Bands indicate 95% intervals on overall WIS. -B. rWIS by horizon overall +B. Relative WIS by horizon, with bootstrapped intervals. -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, with bootstrapped intervals. -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 (dot plot centered around 0). -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) +A heatmap of relative WIS for each forecast date and location is in the supplement. -## 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). + +A. Bias by wastewater epidemic phase (increasing, decreasing, stable) (dot plot centered around 0). -Result table from GAM output (what is $\beta$) +B. Calibration: proportion of observations within each central prediction interval against its nominal coverage -Figure 6: Model based eval results +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. + +## 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: 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 -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, 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: 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/references.bib b/docs/references.bib index fb32fef..3112238 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -35,3 +35,68 @@ @article{Wolffram2023 number = {8}, } + +@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}}, + 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} +} + +@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}, + } diff --git a/docs/supplement.qmd b/docs/supplement.qmd new file mode 100644 index 0000000..9f16403 --- /dev/null +++ b/docs/supplement.qmd @@ -0,0 +1,49 @@ +--- +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 +bibliography: references.bib +--- + +# S1. Comparison with the baseline ARIMA model + +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 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. 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 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. + +# S4. MCMC convergence + +Table S1. Number and proportion of `wwinference` model runs flagged for each convergence diagnostic, by model. + +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). + +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. + +# S6. 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 S10. Location effects on the relative WIS, and the horizon and forecast date smooths. + +Figure S11. Model diagnostics (residuals, QQ plot, observed against fitted). + +Figure S12. Relative WIS by epidemic phase from the sensitivity analysis with a less flexible forecast date smooth. + +# References