From 3ff2d1ea7ef84264a6cb1075c2b17d67998c65c1 Mon Sep 17 00:00:00 2001 From: Kaitlyn Johnson <94390107+kaitejohnson@users.noreply.github.com> Date: Thu, 31 Oct 2024 09:38:35 -0400 Subject: [PATCH] Issue 184: Add outputs to `generate_simulated_data()` fxn and package data (#220) --- NEWS.md | 4 +- R/data.R | 116 ++++++++++++++++++++++++++++- R/generate_simulated_data.R | 131 +++++++++++++++++++++++++++++++-- R/model_component_fwd_sim.R | 46 ++++++++++++ data-raw/vignette_data.R | 17 ++--- data/hosp_data.rda | Bin 607 -> 612 bytes data/hosp_data_eval.rda | Bin 655 -> 644 bytes data/subpop_hosp_data.rda | Bin 0 -> 1122 bytes data/subpop_hosp_data_eval.rda | Bin 0 -> 1349 bytes data/true_global_rt.rda | Bin 2169 -> 2182 bytes data/ww_data.rda | Bin 1657 -> 1577 bytes data/ww_data_eval.rda | Bin 0 -> 1925 bytes man/format_subpop_hosp_data.Rd | 31 ++++++++ man/generate_simulated_data.Rd | 5 ++ man/hosp_data.Rd | 4 +- man/subpop_hosp_data.Rd | 46 ++++++++++++ man/subpop_hosp_data_eval.Rd | 48 ++++++++++++ man/ww_data_eval.Rd | 55 ++++++++++++++ scratch/sim_data_script.R | 1 + 19 files changed, 482 insertions(+), 22 deletions(-) create mode 100644 data/subpop_hosp_data.rda create mode 100644 data/subpop_hosp_data_eval.rda create mode 100644 data/ww_data_eval.rda create mode 100644 man/format_subpop_hosp_data.Rd create mode 100644 man/subpop_hosp_data.Rd create mode 100644 man/subpop_hosp_data_eval.Rd create mode 100644 man/ww_data_eval.Rd diff --git a/NEWS.md b/NEWS.md index 5adf1ca4..eec35d5c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,11 +1,11 @@ # wwinference 0.1.0.99 (dev) ## User-visible changes - +- Add wastewater data into the forecast period to output in `generate_simulated_data()` function and as package data. Also adds subpopulation-level +hospital admissions to output of function and package data. ([#184](https://github.com/CDCgov/ww-inference-model/issues/184)) - `wwinference` now checks whether `site_pop` is fixed per site (see issue [#223](https://github.com/CDCgov/ww-inference-model/issues/226) reported by [@akeyel](https://github.com/akeyel)). ## Internal changes - - Updated the workflow for posting the pages artifact to PRs (issue [#229](https://github.com/CDCgov/ww-inference-model/issues/229)). - Modify `plot_forecasted_counts()` so that it does not require an evaluation dataset ([#218](https://github.com/CDCgov/ww-inference-model/pull/218)) diff --git a/R/data.R b/R/data.R index b196dcf6..3baeeff7 100644 --- a/R/data.R +++ b/R/data.R @@ -39,6 +39,47 @@ #' @source vignette_data.R "ww_data" +#' Example evaluation wastewater dataset. +#' +#' A dataset containing the simulated retrospective wastewater concentrations +#' (labeled here as `log_genome_copies_per_ml_eval`) by sample collection date +#' (`date`), the site where the sample was collected (`site`) and the lab +#' where the samples were processed (`lab`). Additional columns that are +#' required attributes needed for the model are the limit of detection for +#' that lab on each day (labeled here as `log_lod`) and the population size of +#' the wastewater catchment area represented by the wastewater concentrations +#' in each `site`. +#' +#' This data is generated via the default values in the +#' `generate_simulated_data()` function. They represent the bare minumum +#' required fields needed to pass to the model, and we recommend that users +#' try to format their own data to match this format. +#' +#' The variables are as follows: +#' +#' @format ## ww_data_eval +#' A tibble with 126 rows and 6 columns +#' \describe{ +#' \item{date}{Sample collection date, formatted in ISO8601 standards as +#' YYYY-MM-DD} +#' \item{site}{The wastewater treatment plant where the sample was collected} +#' \item{lab}{The lab where the sample was processed} +#' \item{log_genome_copies_per_ml_eval}{The natural log of the wastewater +#' concentration measured on the date specified, collected in the site +#' specified, and processed in the lab specified. The package expects +#' this quantity in units of log estimated genome copies per mL.} +#' \item{log_lod}{The log of the limit of detection in the site and lab on a +#' particular day of the quantification device (e.g. PCR). This should be in +#' units of log estimated genome copies per mL.} +#' \item{site_pop}{The population size of the wastewater catchment area +#' represented by the site variable} +#' \item{location}{ A string indicating the location that all of the +#' data is coming from. This is not a necessary column, but instead is +#' included to more realistically mirror a typical workflow} +#' } +#' @source vignette_data.R +"ww_data_eval" + @@ -57,9 +98,9 @@ #' to match this format. #' #' This data is generated via the default values in the -#' `generate_simulated_data()` function. They represent the bare minumum +#' `generate_simulated_data()` function. They represent the bare minimum #' required fields needed to pass to the model, and we recommend that users -#' try to format their own data to match this formate. +#' try to format their own data to match this format. #' #' The variables are as follows: #' \describe{ @@ -132,6 +173,77 @@ #' @source vignette_data.R "hosp_data_eval" + + + +#' Example subpopulation level hospital admissions dataset +#' +#' A dataset containing the simulated daily hospital admissions +#' (labeled here as `daily_hosp_admits`) by date of admission (`date`) in +#' each subpopulation. +#' Additional columns that are the population size of the +#' population contributing to the hospital admissions. In this instance, +#' the subpopulations here are each of the wastewater catchment areas plus +#' an additional subpopulation for the portion of the population not captured +#' by wastewater surveillance. The data generated are daily hospital +#' admissions but they could be any other epidemiological count dataset e.g. +#' cases. This data should only contain hospital admissions that would have +#' been available as of the date that the forecast was made. +#' +#' This data is generated via the default values in the +#' `generate_simulated_data()` function. +#' +#' The variables are as follows: +#' \describe{ +#' \item{date}{Date the hospital admissions occurred, formatted in ISO8601 +#' standards as YYYY-MM-DD} +#' \item{subpop_name}{A string indicating the subpopulation the hospital +#' admissiosn corresponds to. This is either a wastewater site, or the +#' remainder of the population} +#' \item{daily_hosp_admits}{The number of individuals admitted to the +#' hospital on that date, available as of the forecast date} +#' \item{subpop_pop}{The number of people contributing to the daily hospital +#' admissions in each subpopulation} +#' } +#' @source vignette_data.R +"subpop_hosp_data" + + +#' Example subpopulation level retrospective hospital admissions dataset +#' +#' A dataset containing the simulated daily hospital admissions +#' (labeled here as `daily_hosp_admits`) by date of admission (`date`) in +#' each subpopulation observed retrospectively. +#' Additional columns that are required are the population size of the +#' population contributing to the hospital admissions. In this instance, +#' the subpopulations here are each of the wastewater catchment areas plus +#' an additional subpopulation for the portion of the population not captured +#' by wastewater surveillance. The data generated are daily hospital +#' admissions but they could be any other epidemiological count dataset e.g. +#' cases.This data should contain hospital admissions retrospectively beyond +#' the forecast date in order to evaluate the forecasts. +#' +#' This data is generated via the default values in the +#' `generate_simulated_data()` function. They represent the bare minimumum +#' required fields needed to pass to the model, and we recommend that users +#' try to format their own data to match this format. +#' +#' The variables are as follows: +#' \describe{ +#' \item{date}{Date the hospital admissions occurred, formatted in ISO8601 +#' standards as YYYY-MM-DD} +#' \item{subpop_name}{A string indicating the subpopulation the hospital +#' admissions corresponds to. This is either a wastewater site, or the +#' remainder of the population} +#' \item{daily_hosp_admits_for_eval}{The number of individuals admitted to the +#' hospital on that date, available as of the forecast date} +#' \item{subpop_pop}{The number of people contributing to the daily hospital +#' admissions in each subpopulation} +#' } +#' @source vignette_data.R +"subpop_hosp_data_eval" + + #' COVID-19 post-Omicron generation interval probability mass function #' #' \describe{ diff --git a/R/generate_simulated_data.R b/R/generate_simulated_data.R index c3582bd5..df74cfd5 100644 --- a/R/generate_simulated_data.R +++ b/R/generate_simulated_data.R @@ -59,6 +59,9 @@ #' infection feedback into the infection process, default is `FALSE`, which #' sets the strength of the infection feedback to 0. #' If `TRUE`, this will apply an infection feedback drawn from the prior. +#' @param subpop_phi Vector of numeric values indicating the overdispersion +#' parameter phi in the hospital admissions observation process in each +#' subpopulation #' @param input_params_path path to the toml file with the parameters to use #' to generate the simulated data #' @@ -121,6 +124,7 @@ generate_simulated_data <- function(r_in_weeks = # nolint sigma_eps = 0.05, sd_i0_over_n = 0.5, if_feedback = FALSE, + subpop_phi = c(25, 50, 70, 40, 100), input_params_path = fs::path_package("extdata", "example_params.toml", @@ -322,12 +326,35 @@ generate_simulated_data <- function(r_in_weeks = # nolint ) ## Latent per capita admissions-------------------------------------------- + # This won't be used other than for the unit test model_hosp_over_n <- model$functions$convolve_dot_product( p_hosp_days * new_i_over_n, # individuals who will be hospitalized rev(inf_to_hosp), uot + ot + ht )[(uot + 1):(uot + ot + ht)] + # Also compute per capita hosps for each subpopulation + model_hosp_subpop_over_n <- matrix( + nrow = n_subpops, + ncol = (ot + ht) + ) + for (i in 1:n_subpops) { + model_hosp_subpop_over_n[i, ] <- model$functions$convolve_dot_product( + p_hosp_days * new_i_over_n_site[i, ], + rev(inf_to_hosp), + uot + ot + ht + )[(uot + 1):(uot + ot + ht)] + } + + # unit test to make sure these are equivalent + if (!all.equal( + colSums(model_hosp_subpop_over_n * pop_fraction), + model_hosp_over_n, + tolerance = 1e-8 + )) { + cli::cli_abort("Sum of convolutions not equal to convolution of sums") + } + ## Add weekday effect on hospital admissions------------------------------- pred_hosp <- pop_size * model$functions$day_of_week_effect( @@ -335,12 +362,36 @@ generate_simulated_data <- function(r_in_weeks = # nolint day_of_week_vector, hosp_wday_effect ) + + pred_hosp_subpop <- matrix( + nrow = n_subpops, + ncol = (ot + ht) + ) + for (i in 1:n_subpops) { + pred_hosp_subpop[i, ] <- pop_fraction[i] * pop_size * + model$functions$day_of_week_effect( + model_hosp_subpop_over_n[i, ], + day_of_week_vector, + hosp_wday_effect + ) + } + + ## Add observation error--------------------------------------------------- - # This is negative binomial but could swap out for a different obs error - pred_obs_hosp <- rnbinom( - n = length(pred_hosp), mu = pred_hosp, - size = 1 / ((params$inv_sqrt_phi_prior_mean)^2) + # Use negative binomial but could swap out for a different obs error. + # Each subpopulation has its own dispersion parameter, then we sum + # the observations to get the population total + pred_obs_hosp_subpop <- matrix( + nrow = n_subpops, + ncol = (ot + ht) ) + for (i in 1:n_subpops) { + pred_obs_hosp_subpop[i, ] <- rnbinom( + n = length(pred_hosp_subpop[i, ]), mu = pred_hosp_subpop[i, ], + size = subpop_phi[i] + ) + } + pred_obs_hosp <- colSums(pred_obs_hosp_subpop) @@ -381,6 +432,18 @@ generate_simulated_data <- function(r_in_weeks = # nolint lab_site_reporting_latency = lab_site_reporting_latency ) + # Create evaluation data with same reporting freq but go through the entire + # time period + log_obs_conc_lab_site_eval <- downsample_ww_obs( + log_conc_lab_site = log_conc_lab_site, + n_lab_sites = n_lab_sites, + ot = ot + ht, + ht = 0, + nt = 0, + lab_site_reporting_freq = lab_site_reporting_freq, + lab_site_reporting_latency = rep(0, n_lab_sites) + ) + # Global adjusted R(t) -------------------------------------------------- @@ -406,6 +469,18 @@ generate_simulated_data <- function(r_in_weeks = # nolint lod_lab_site = lod_lab_site ) + ww_data_eval <- format_ww_data( + log_obs_conc_lab_site = log_obs_conc_lab_site_eval, + ot = ot + ht, + ht = 0, + date_df = date_df, + site_lab_map = site_lab_map, + lod_lab_site = lod_lab_site + ) |> + dplyr::rename( + "log_genome_copies_per_ml_eval" = "log_genome_copies_per_ml" + ) + # Artificially add values below the LOD---------------------------------- # Replace it with an NA, will be used as an example of how to format data # properly. @@ -419,16 +494,27 @@ generate_simulated_data <- function(r_in_weeks = # nolint TRUE ~ .data$log_genome_copies_per_ml ) ) + ww_data_eval <- ww_data_eval |> + dplyr::mutate( + "log_genome_copies_per_ml_eval" = + dplyr::case_when( + .data$log_genome_copies_per_ml_eval == + !!min_ww_val ~ 0.5 * .data$log_lod, + TRUE ~ .data$log_genome_copies_per_ml_eval + ) + ) # Make a hospital admissions dataframe for model calibration - hosp_data <- format_hosp_data(pred_obs_hosp, + hosp_data <- format_hosp_data( + pred_obs_hosp = pred_obs_hosp, dur_obs = ot, pop_size = pop_size, date_df = date_df ) - hosp_data_eval <- format_hosp_data(pred_obs_hosp, + hosp_data_eval <- format_hosp_data( + pred_obs_hosp = pred_obs_hosp, dur_obs = (ot + ht), pop_size = pop_size, date_df = date_df @@ -437,6 +523,36 @@ generate_simulated_data <- function(r_in_weeks = # nolint "daily_hosp_admits_for_eval" = "daily_hosp_admits" ) + # Make a subpopulation level hospital admissions data + # For now this will only be used for evaluation, eventually, can add + # feature to use this in calibration + subpop_map <- tibble::tibble( + subpop_index = as.character(1:n_subpops), + subpop_pop = pop_size * pop_fraction, + subpop_name = c(1:n_sites, NA) + ) |> + dplyr::mutate(subpop_name = ifelse(!is.na(subpop_name), + glue::glue("Site: {subpop_name}"), + "remainder of population" + )) + + subpop_hosp_data <- format_subpop_hosp_data( + pred_obs_hosp_subpop = pred_obs_hosp_subpop, + dur_obs = ot, + subpop_map = subpop_map, + date_df = date_df + ) + + subpop_hosp_data_eval <- format_subpop_hosp_data( + pred_obs_hosp_subpop = pred_obs_hosp_subpop, + dur_obs = (ot + ht), + subpop_map = subpop_map, + date_df = date_df + ) |> + dplyr::rename( + "daily_hosp_admits_for_eval" = "daily_hosp_admits" + ) + # Global R(t) true_rt <- tibble::tibble( unadj_rt_daily = as.numeric(unadj_r_daily), @@ -453,8 +569,11 @@ generate_simulated_data <- function(r_in_weeks = # nolint example_data <- list( ww_data = ww_data, + ww_data_eval = ww_data_eval, hosp_data = hosp_data, hosp_data_eval = hosp_data_eval, + subpop_hosp_data = subpop_hosp_data, + subpop_hosp_data_eval = subpop_hosp_data_eval, true_global_rt = true_rt ) diff --git a/R/model_component_fwd_sim.R b/R/model_component_fwd_sim.R index b5449646..956e574d 100644 --- a/R/model_component_fwd_sim.R +++ b/R/model_component_fwd_sim.R @@ -422,6 +422,52 @@ format_hosp_data <- function(pred_obs_hosp, return(hosp_data) } + +#' Format the subpopulation-level hospital admissions data into a tidy +#' dataframe +#' +#' @param pred_obs_hosp_subpop matrix of non-negative integers indicating the +#' number of hospital admissions on each day in each subpopulation. Rows are +#' subpopulations, columns are time points +#' @param dur_obs integer indicating the number of days we want the +#' observations for +#' @param subpop_map tibble mapping the numbered subpopulations to the +#' wastewater sites, must contain columns "subpop_index" and "subpop_name" +#' @param date_df tibble of columns `date` and `t` that map time in days to +#' dates +#' +#' @return a tidy dataframe containing counts of admissions by date alongside +#' population size for each subpopulation +format_subpop_hosp_data <- function(pred_obs_hosp_subpop, + dur_obs, + subpop_map, + date_df) { + subpop_hosp_data <- as.data.frame(t(pred_obs_hosp_subpop)) |> + dplyr::mutate(t = seq_len(ncol(pred_obs_hosp_subpop))) |> + dplyr::filter(t <= dur_obs) |> + tidyr::pivot_longer(!t, + names_to = "subpop_index", + names_prefix = "V", + values_to = "daily_hosp_admits" + ) |> + dplyr::left_join( + date_df, + by = "t" + ) |> + dplyr::left_join( + subpop_map, + by = "subpop_index" + ) |> + dplyr::select( + "date", + "subpop_name", + "daily_hosp_admits", + "subpop_pop" + ) + return(subpop_hosp_data) +} + + #' Back- calculate R(t) from incident infections and the generation interval #' #' @description diff --git a/data-raw/vignette_data.R b/data-raw/vignette_data.R index 38c61081..8f4c2dbd 100644 --- a/data-raw/vignette_data.R +++ b/data-raw/vignette_data.R @@ -1,22 +1,19 @@ set.seed(1) simulated_data <- wwinference::generate_simulated_data() hosp_data_from_sim <- simulated_data$hosp_data -ww_data_from_sim <- simulated_data$ww_data -# Add some columns and reorder sites to ensure package works as expected -# even if sites are not in order -ww_data <- ww_data_from_sim |> - dplyr::mutate( - "location" = "example state", - "site" = .data$site + 1 - ) |> - dplyr::ungroup() |> - dplyr::arrange(desc(.data$site)) +ww_data <- simulated_data$ww_data +ww_data_eval <- simulated_data$ww_data_eval hosp_data <- hosp_data_from_sim |> dplyr::mutate("location" = "example state") hosp_data_eval <- simulated_data$hosp_data_eval +subpop_hosp_data <- simulated_data$subpop_hosp_data +subpop_hosp_data_eval <- simulated_data$subpop_hosp_data_eval true_global_rt <- simulated_data$true_global_rt usethis::use_data(hosp_data, overwrite = TRUE) usethis::use_data(hosp_data_eval, overwrite = TRUE) usethis::use_data(ww_data, overwrite = TRUE) +usethis::use_data(ww_data_eval, overwrite = TRUE) +usethis::use_data(subpop_hosp_data, overwrite = TRUE) +usethis::use_data(subpop_hosp_data_eval, overwrite = TRUE) usethis::use_data(true_global_rt, overwrite = TRUE) diff --git a/data/hosp_data.rda b/data/hosp_data.rda index 7595c3bb63da37ec0c78e995b9ea137bb3e08967..83e0eeb833d8928bc8825143775be9bd48780ed4 100644 GIT binary patch delta 601 zcmcc5@`OdsDJsL#&@oaiIC5`Y&2I+Qef8J>7iTcI{|AGKauW56R=!HsxGW%La)GV? zZ>Wjf0+S5pzai`)M)$ax*^d{WcJyJK z#1d?>wo8-MX}$AT^{z*Yp1WE6bPlQ8v1Aozh}ItI@PzXZC(fVrsh+EL(T2xP*E9@P zaB{}9L@9N!N;)Zh{=Lc7d-kSOK|_%2H6Qre}6CMJbm%vv4>wFgP(TaA@q17HDHo;bLd2Zb&%T%*G>O zfJBMZi+HlCu(M}5?opb+knLoa=aqQhv^Z&ttFfg0%r&}#UN4Q8DNX9x5S`|Iqta+o z*`cm@L63>Ml%{|5&wKOyNWvvSKk2l4jUtEZ`MB2I%JjcdcqG_4pWo_9#+idtT(u?! zb>*cl2yrwEab!}Ma|Q^`y3PcFRov6|=M`M}XVm=Ym2}X;?YV2#&e*LWlv7$Y$x{W4 dOmy^16&eMCIhJp{az*vOc#*3`zywe#0svwC2Ppsm delta 596 zcmV-a0;~Px1m6T6LRx4!F+o`-Q(1%es|^4Jy#Lq#b7%lZ|Nr}c|NW637=NjyH9V$( z(?da#paK7?p!EiTZ}n7Yey9xqWYE!}NBXLhKmY(VXaF(*0MVcT00E{Ms;T;+p@asS zWXJ#-7@7dc28JPs05k}qp|wp<6!e~;&@=-<Bnp1y)Sae%4cZo4K?4@a5SBL(AU@KnGFhCPc085~Fb0Sb zO*33rvI>P+U#U!1l@|1E$I2R`^Ti^Fn~(zVO(50%Nvorbf}8PdFMmSY6c_;jh!hA& zP-YQEW>FV^yC|u1dRd5|ifIF403ra00tmtWl`KFSYruN|7%&VN02#~~!@dAJ{9$!E*;L zaniU^OQIY%CIzaj8qGbJhO3=Dg*24`s4qNJp|bQgTaIU>30`AU~pjI073@F15BI!_?nXx zmFBOwB*nnXz-uD^YgGobLgtIFzVAa7Y!#$D<)7z`)C_>^6(dfkDdT z;wz4s0U5juj0+ebfX$I>iV~P9WpV;&^a94kL2n;GCDACpDXN-=Ee9B!m}WE_QgCq4 zmD-l~@bvP0xdlJ}a2PB+(k*4vX|(jr?Cy=5gVpscMNNv_{ggPTYDaz9wDa_|(4#yd z&t#sS=zhBO#r?Llb=@29FF6r$)UrbG#3|p^(QKRjxV~nDKXlq9^%>_Z*ShIbIX$~?mW+ks9y4Ez3Sf1S^n9A+x9Pv@tBmCw`r26%H%noDnM#df2LNRG3w!_o literal 655 zcmZ>Y%CIzaj8qGbG{4@whkcHn7T?7*YIxPU>@B*~zG!DXsi@K=ruwuX}nC;PG~U*gH&Two;8 z)9_a_>+1&wtblJzvcxfMi5|yAS{p;VIGQdr$#j{lys~<7wAS8k%gUKMZ~l_93fuD4 zG&?){_QtfHTSZwHW^K4W-K%){nJHT|{gvN1G%YcbJ-Aiz!pHa9!rq?ki8A_Ix;Zv` z2JfxS>Dl?!vzK-C?ftMY6xGT)DZ|G3bJ^0!xDji~8-~UNMb23swlR zC^#{=ajXk*1;L_4TDuSKoH={r8q=_Bser4Sxs9E#i|$x-=IhJqy~gu*%)I<-^O{L4 z5=k>Ax;9Pa6jNbwP;zMG=m=<3ahb7UVq57!ub?SQHBpEq$6MAIp3dE|@m5EtlGfgR zUhcA?37s!GUe{>-e17p?*6o$8p}Cg!5#G9<`YK!d&Pi`x;;=ckGca()uiJ6uZXx?- z7@GF)?{{A;yO@<=qIAW`}M(l1^_{uAeR6D diff --git a/data/subpop_hosp_data.rda b/data/subpop_hosp_data.rda new file mode 100644 index 0000000000000000000000000000000000000000..29de916840b52c20232fad3acc5bd125460e1fad GIT binary patch literal 1122 zcmZ>Y%CIzaj8qGb6f`V7#~`g+fBgUNnk(V|f#CcfhyCUE7C1NrFmNz1I5;?7VBKT$ z@STO$SvCW?+{0SSS8j2eG0Wsc=464(3~fHMQZ6o=$^SRhkeP!=z(i@DbFf43SC+ql zvlJIQFi7&+Dm%Vn00O4qOPwB74g3Le8YYYjfC9D-KtfQ|6-@sRWnj1f0|zhb#HtI~n3;o_@|zrj}+j7&1l`UlUNpYli5YpPbr%AlpBQfvBI zg2jB4ZG3lLoHympT+NqfgO^)PUTt-Gy{+es*_*G5!$ zuU{wfz5ml{j;?mn>CH1%AW$9Bfy4-rzEzjC!bHjnY#XLoa@xwMuVM+1{z4|2FP#)%2GZSiJj7#^LNN%|~~=N-oVUS$4dB z;ra8@oxbK1TfH(v!n6dAT%MwJUh>fjNrSGrBC}umtm^u^!EFw&eAuBir)?x%Tb7A# zm~!gGydP6Os9L`Cd9KW-ppwAJ!Q{p{!&PBzlIEs=VUiDK?B?|8oODuixly{YSAe%m zP*Vh_f{F*%*$9rM8jUk~3)cj@x+!V6D4Vxz>P(xW;_2=5*fY^X;Kjvxy;{Z>uPwHf z=`uavl|Ad(&b?NXi%Qj29<)$DZF;tNwe=mz-UF#Vy;>$qGCpsSoE^tp8CtfcUrsRS zh1--hQXEQs0y^hiWvewO_S!`-((!y0MCPokmIfQ0zm~!4?*{fia<%-t3e=@J_ zKjD@5_`~n@`Q6$_-v6-C^snBU^GRk-%9IO_cE`R8o4YP3*zVf7Imt}HUcu9*ob6Lt zp6=|bxlr>Tr>{(=FA(~fEblye-mCD>Nx#Um&z4%wxv%XMylu*%;B%oCcX{VM7FVB{ zv((rt)2ke4_;%YlmAiS@$-dold!Njlhng+REW<4|Prg(VY%CIzaj8qGboWI{!ok5MG{`mjjaaY3s1Ht(}4*SdRuW)j3VBlb2aBy&bz%C(v z?}P6L#(540URWrKddw77Yi2aLxO}nV{!k7R6Qc|h>0th#`d2F#vq&Pz!Dnr|8Hmj<6;M?Bwhw1HZ>-o90NlJo8*#?B}ohnJPiyc zyhaQR7o>oqzFw19X1jb^?91>KD010Bis51huT(}=xhxfEXD@ zm*33lRGH*SIbvI|>fDxV9H9%XH7CtovSFs@v}r+?($uojRkyU9_bEByGv$=emNPzA z&iQ<~;OVJr67*0_ZIatsli9g@-=*vKu08QLLU(@fl8kt_^^524ycRlZx1qp>35WbS zB+n+OcB&uoaQqaawmvMBDPVz7TUN@1B|fdCYI}W}>zx*yvk?&1+%c>ZX-$(UH>ID#5o!YZIe@Zc6epfs`*QeiOwPYX$^HMFdF; zCt67<3aopYS;@hDGEi3C!9`51bm5^@(;b%Y_bzMX3RFomT(hlnn;w%V+u78e!Kyzt z$lG1m!_dgUC_ZnB;{xRlvp`|hWtVQfG5R+9teb|B>dR;O>pi2M>fQ@o_e5D>MhoAS zo@FQ2rm;>o*!S3S*^O3_1s#GeSDG|bUK%R3+^rNmp>lv>)sbxJKlQoq@{8xKp0{ev zEzY>Br#}{nhL+p9%l`dYV-u?Xc~2EX@rxxNe$BOd`?m4%7lxk)Y89-kukU&Mxbb-O z-ld_jNh_0>POop7^mJ#cU}46p^{iT793Cz8zICMZRmrctEIZwGE!tMC_IknG^&q~j z;oE`l3yLhaK3SQ%DsolSTD=Ani>Rfm_4)go)&&Qh`jt@V-CDq2JykDwf=tlE12RF6 zu54YIx-vKD?i6;B9apw3-Kx7X>+sp7)3{AUI-CSlJgz9&?A*3>smlVxfP5+;lr zr50)ha=Dm_`>s{n$s06p$x@R^Q$DC_KJ*ISGi}nLC~wWz{z=)(4h79zxoKX>E#DijLBQzZ4l`cuPSxj{ z9NfoM@(xdx_H+z5)y3K3BFd??!c*XaA}14zg7b#A0TWkpYFyS#%)WIevQI2T%jM9M zW4GMY%x7>htuzqoJEGgG9hI_GM)c9%fEh7f3~K~gJQ6INUhj@9uNAzu;WQ8NXiEm4 z_pAqec^R?`iq|C=J->8Tcj;OU0gv+29Sn2%CWJDq?^H58e*1das&y0W-!fi)%6tO|PqF*= z=c~_QpMH&b-fnA#&(${nCfF7j$n&h)pFi#M@A*>ajW^z5o*?_dJiLQHefH`ax$f;p zzwVJ)?|d}iAZzHt(2TtWpBl2_P`_k`NXk{AB&CeNDeiU%ks@dLEKsyFc~QkH>e#i(D-N ICV+}f0BMI-kN^Mx literal 0 HcmV?d00001 diff --git a/data/true_global_rt.rda b/data/true_global_rt.rda index 399520388d6873a1278e61be94909519651c969b..c1a6d882047826236a0403955ff3e97b28ea2134 100644 GIT binary patch literal 2182 zcmah@c{J1u8~%+M%P_+rOp7&ks6>{KYD6;Du}=)fzJ^Mu zkgUtf*008I* z-~!=z>b?@70f06caUiG29E}DmLz5|LJlsJHrXAn`@@%Rs4i9@=f&t@Yf19arp#OIV zEx`dX8=&eyBJ`y|xG`HG!8VsI3k za5(1RK@4b#x*snKcrFQ+tPSt4;(mL7KSV`S2iPGrs!&j`x|vGS|0!Iy^%d0*PA4es zsowRyuLsRZ<(|Z!&j-Jw)Z;JeN0;v(d=ASeDA53oc(S7a%Rd~=$V)_~!M*0@Ad&C3 zG=(!}6C`!^T0Kcg$&p^;wha|EZ|IMYXFYS{R7DW}6U_jqjX*k4x4f<6927aAZF}>J zvJX;DT2#-_1giXq)Rw|3C}Unjj|Shdlf=ease9d$iS@;Pn_Zst*1IdeB6FCo^J6y{ zCZw$XP|t^rba!!6kpSGGFucJzLHd~oj>c@b4UNsGAS5yic}NgO2=(LwbpOAeBV#sY z7X4!`r5T)f6c%gx5Z1#pm*eBS<1zD-fxhOr{yu|(YAE-IqEiZLx1?ZP;-AAlO>+T5 zF6}-q`12m1a-dwd!9x2ZCqd98ZD!c|Wm~^1PqRS8!u3Liufi(_`k3w-Eqz@#{KZgd zVjSOK-8$f}cT9n?S-TmEFh}sGy+jMWwurw+b^tHbrO>W7Go5+st>_J-kWk`HLEvVl zQeqI9J)NDjD#tHQ?uHresOChjU9N{rCmm;HlJuhU&{<31wDT=-H!$SgKh?~=?%%lc zu7l7)IK{Yp>1#u!5o6(Lwa9Id>!VgLa4&P#{?J9nY?_e3i{a&r(T5w#XKtuTy)HMY zmpYcRp2tr{9B?vxBMFI+WXsp0E``_XZW+g*(JjkO<-|{wCVq%ewWDCC?jQSIN0t-No+n7h#!xow%Y(2mS$VRHsy$i$fhW}cClXU3^EjE%inqO0NpUG zbmBG^5)8?R(5li8t>F~oab>7RU@24nDpL7FnSQpv<;VayH#C8U(0S%#2+b_|@dvM8 z?kmvjqiwWUTzG*M;>SUYcD2nVk6pH7Mh&A_5}tL5+e6l78kxj^yYrV_^ZL?LYAudM z^gLau_^Ir7%Y~i~@T#?^UZVmTlKEcXEy4|}9l*ox*II+9^Q2Q9pXkw0v5V$H36i|m zlPrqx$}MsHO^**$a zpu90Dk5!tZ>=PA7O9iRbzsi(kcgij?{ThY?qkWFtljgt>0FrvPL%U4MaC<^#4m++N zo6dpJO}E^s|a$qN~9<_xZkE7sQ#KG;<*l z--``OHznecoUaq(G`T4oE%{jphw%gg>E@=KS^#xjwztR%R#$cCo~}R;GWY;$3mwuf z@hR`Uk}i@7ZHSy^{m|O|=x71E@LkE2d$4tCh1aMTAd%HCPIC&=?#in8jG<}+Hj6p* RKX;GMM*HL{FdfP2{{pa};P?Ol literal 2169 zcmaKnc{tRI8pnUL7<*x?W7n8LcHnr?NXXNx z3~{p*j4_#SjKRR@1Xs-11JVqc9K{F5*t70VoC@{kmRZ|5vXOCGE0f8|Zoz~xBz}^j zhQo+tB0wij8yQp54OQ3bDI7!v?--h0fkNejFhdey1oUpp7!3u$!C@!{#xMh5h#LZM zwW+zP3a6YxZ-*8N2UV$vs>1M$#GzCihBwBc@epKy2P&~8&=3j$nFnehFcctoHY~#x zMk8YoWJ=L3X-?_yG@dvB%g2G(00IGsBLD;j-~j-P1{`v@1R%CJ2#zccC?JuQD(vk_ z7@fw4WgNbsaPZeSGyui`5U`9nqeg&?Fh0D&02B@bgMbidkQg$E&6{b;d*4{>g#7<1 z`4cL(o7VmL{Hymn3rDB6Bn*iGCK;Jo6!D{RIsfW4?*jm1pVUseG@uf^z*_$>bumr5 z$tYaxlE$;9$HJpKS_5Of7zE7FhFkdEUr6@fDBa43!)7*j=-D~JA)&Cp%lrrE(4h%1 zg={1}0LW4MyvC8)sMFDBH0Q{a7}Ub^my1jDU$9vIxVd?s#dqb}H!30h1_^}&&6`W- z_Cb>{kt0&F2-SW`nc}?B^7IpL00}#1p}V;tx^un_UJbVD7*W?%Wab>%R&5T#*??(6 zG&MEQx(yEM<&jmC%>_Ug`~@@=#;7-UEz`y|AsFn+t`7}N8_pqLz3(Q_dfdFbWWN0* z&7bKAk781Wk~~lX2q$>8k+r&8icq6=9^O%Hqo8Ri~ z>(+$PMTL%c*4*j-a@Mj}t0)@dys4mF-IXLZHnQvxI>Ey$kuHCHm=91nI^l^a-hnyL ze5Hi#SK-SBxs&&3H3Hc|c|Y&T3s$&@MxA>)Q052(S=dY}B*4}eei2ocwO<<|Qj)Ap zKFfoXhAX&HqHBJmneS5I-LdL<1-5IkIP*>ujZMX9NC;O;E43jcVTl&2g6Ex5*6@K} zX>YE8jy#&nU+7-yHvU%Z4gRpkP`_0cT|azERPn;cp7jUS*S#kzgzrMR!IvVcGk*a~ zy3PCVi}KXGv+|L+^a0ll->xxy%4;t`HpQiWWMiInBL!_MaZ29(T$pT+4qg#;Z$(Td zCa^?see%v})|C4V)n0c)18i->R}p!6_e^Ca__X6j#lC4J5d^Ve&1r|5eby5B^4cL~ znbq}G!mg$Lt&!DP*@>U}$4+ge#x=@4b7|mD)$gxQ$R)-Zu0IT6+AcX+hnrx_hS>UF zeGQ;TL_BY|=+`OS^Z(8#-8vBv=*%xRfXZA*fC}b_jki(Nth+;ezT&C0bB#sDhB1(5>P+S1$=^sp}sp8tUv^>`kVL1 zIS>VFJ>HWP^lPM|wK6;VP#vwI2u|P-d`B}i)zq>%fbY#6A-@+6d3~kJprknwbO-AF z6KUx{j`b)z@km^~)Gki%h(9Fv8CuUTBt>y#km+)bKcj=QM}FP&_9_l5Y;^Boe|2c~ z0L%BPisc;pYxb0#+7y`6l(OZwD|*F7ag-iS z^q#cB$(uf&Pw+mp+-Ch8nYPL?_G05}>csVbE<);A7OY5F@z3%UIEZMMr{q|t%kol% zN(n<^^)?OQ`SGiB+FgD}oy_%SHuln38qd_+CcOsVDa?R9VZpMQ*&am|=T1I)(%UiE zV0JQSOU@cAaI^i!@tbA^^#;D+G27`-`YWtbGIzYP6+Ts5-JwJ9{K(QRU`phP>qIQ_anA=mgj`1Lp z*_j^`51wZSUv=@L*Y5_+y3;PQV3Lc&_84_5SM$BjP2{BqvA-5trEiu4YnG3G{fCpp zk$%6MmVPKUbS?4tcwEg;FX&-gbi4q5e058zyVqaCUogTJ$pX)6Ab_|vN10NAgz9k~ za`!b!8gM<3`}U)gL6Q0`wW#3Z{?q9LoLws=a^}9tZ|t(h+0hr_gg_|nysNgdJ~*1~ zWS_nQvDBj`1;#6{An$laR}p4RedF^6bWOICc_QDKK7f`Rgz`j&tRjK{pWsdO$Lh8Z%#m463QcXuwtTZD9@C;S&4#?IE^3WL)yCwA3jOc;QGZUPqzra z#4sQUoe|ra#V*r14qQcK)u;#sSBT z?0Y>JRgLem9OTJX&oxfV+oJ>nR{92e>jD(jC+eP&qfb+69wgJ$H1v9q%|lHaP|#_hdW|v26UDUb|}JrmPV#TpF&00001 zpa1{>&;S|%rkZ9`3|nrQ7QeJBNp?6DCT`A|eYv@el^RL@WenqzFoD z;3!J8BmldtAfbPrNQL;XhC~<{1-gT#A&7t@)MbYuPIFbVp2XMxi=N%GU>Cn$Yn27Plvaf|4obdp&Rk7)8hCeuY8@Yu5zfoU2AWe`kND5>0yLmA2Cwg z00kk;7j_vNW`dZ3n7R7bbTfefEFNH?@*9wzIH4Tu$a?Typ$QPw5TRrsgrJDYQ4<=S zC6Pdw@?#1ryRz%VkKT<`bS^@f11#AMh*%Y2DQ1byQqRnno>SbRko}Tbt z4&>mnS9SUO*jJ_yx1#@L?70WcN}FhtxP`UC((W`?JO8oe**AYEkX<{M=2oFiDgpwEAwohQO7Lc_Hynd+Gzr96q`jDe@7lhpD?CV}FE(maq#N-62006kAg z1JpG1ntGW4+6_VKXk^gZntFi32AVxXKpFv{0009)pn8A?K+pgH4Lv|KGy$eVOqvWz zQ~(c9&@?un4IZOF(0Yc{007VcX`lcA0MkZ;AO?T{F#|vu05Sof00006fDixx0iXte z0000D0000oXbgY=0MGyc00000000000004#Kq8V6Ae$8SsBI>L)Y0kUpm{(Cr~?qw zKn(x}gG_+*o})uR$N&RFO#^Bg13&-(4F{+IWXNa_Q$RgPo>CSdk|-d>1d2?_nSw=- zW+9U^CLtunGG;-UFlZWqc=8IU6;`na2KKjIYV;w*f&ea`{vd)P5Djhx2vn-BrhpWn z)d50#P>>4iFbD=Jc_0O`t4c`(Rh~n63Iq6m-jw(eMSv(%t~{AxQiT8u7_F3&fKrRj z>Od&(FoLz<6dE~ctR|6C@>zSiO+bL%-QdE;6}b>t1O^~)pao=*CmGpV_-pO zo9AB2teXs2S4;k1 zvx*ul^xFl*xp?uoiYCjbV`LW)j^Lm$bMT+M0j@^+t+a;^0GLQu*8zuBg8r1Em1&>;s+v#m$-MFc7|d}ozK|De&pN+{DFp!Vj@tSBh*91H)s?` zURj5c9ttr}V%#rd7+Vn0QmS~Odlb2Dzg)0uJOV5UH6oOrhx$YU1&eu1tY;<(q=cR^ zLnpIuof&Pw+ZmEGM_-!b^Cx4A$>FTB3K0UOc>sqoQb026g}}xLOT<*5S||XUp;V=N zYf9?FN2S7W`CII#qd4gOB&-poxgsN6R{j%OGo+ws*7H$TW^1I*<>kZPr)t_|{gTV3 zkm(AeGlH*xpmm92Rjk28Mm?<`<mN_IAnH!l6tX{p~X<_n$0j)zTE9pcny2zyk~< zf4$@BUPn&s-L$g_#rBm|i5o|zbL=!%^TK+I zCJpyWJUt#SH;MX;9t%Q!Y6CaAQl%&m9oG?{n8&T}$C2*lm2KKtH+8ADjSxUCi2xR7 z`DX_kAV>h}rf<^px82Cwd+?u4A6aD`T2()_Z~$*mHWd?3uUcqv9~EfM*0E3$%!JGJ7ORr=!atwC51zg6YAt23$p8RS4nP1~2sm%l3Ds2e-F^gE?MaZ;$^nKyYp z3j4Eid;Y0mr$|R}B!$^W-yEqsU3tUN+rE^G#j%)jh!$53mq$qQ*fRe3=(7`CKE~MD zJK|fkz#|H~r;_sU)9n8UfC{DY_o-m44HlPU0V2Ot&M%B008u7#WV&4HhR27~7>Bv- z-l9YyKH+!EPkLQ5pnPuI0OBf>HUXMo0#UtSw|2#~&-0`}+`iY>>D}4vAQ65S_<5<1 zm;`+*3tmoVI{X4q;TtlgcX&WLaTYhoIm7}cXTVsMA*+gxtVFmBf5qI9P81{;g%vo! DlbG>! diff --git a/data/ww_data_eval.rda b/data/ww_data_eval.rda new file mode 100644 index 0000000000000000000000000000000000000000..176a52b471dbadb3a38f39b5d9c98f9beb7ba6ef GIT binary patch literal 1925 zcmV;02YUEIT4*^jL0KkKS%_c%%m4`zfB*mg{r~^}|NsC0|NsC0|G)qL|NsC0|G)qL z|Nr0r|LxEPT<`*vd%4-#&t}(1LXk9@Jx@uKPbfT-YBcp8lha7VdY+@j0py!V4@v3` zGzNgvN2#MmY79oy4@e%Ts1H-r^#B0%9+O5*JfXEcA*PuPH1!Wt(F~hOqfDNmqd>}O zp)d)i27+MJW`jVO4I=@dnqZg#7$Zyo01^6t011FK8ek>>00hZ^00h7Uc>pv9p%P6L z^&3xA{U)Q-!8WEN)HKubiIXR(nunr1Q$|6M#A0Y*Lm00000000000000000Te(00000 z00|_Br|MJnPs%n$F{$KF6xvObBP40)nr%-g(0YxjrbnnfN2nTU>Ulxx49Om$>S*;g zQRIx&@{J81pfmvTjRuCA4H^vrka~kn4LpgXKo3eP!Y&v9>BkV@4DJcg>Co=(ojZv- zcOBiG$7gqU2LN%?x(7p(1u>aS)6G=caYWQPP6XhLslny)&6weHE~KTVu0ly-_d-BT zhy^^556h&G=cW*n9gPPeBomAY0?0n#1&=1#IAVeK3q$|&|I^q_3IrsKn7R@P*T51+ z9f}T_$OWQnMnDK5GC%$ly5`_)H3HL$qxYSHIjNOxQA`>b}r(Ewz0JixgB>(3b!dLnx6{EpKBFI?XodQ7>{;65xUw*G}_%gdqe*KwnVrdk4fQ z%m^NudE_Oa&}#9Ok8bLC?WFJ|PhUyvV&d@`v;}*%bkdm*!Ng1))5-{rs+ZP0U$LdE z^e5k?YSuQ)bc~}gZkY3a0SUqr>?D|G229qg)3h+5AiQk=BvHq1n4m_@U+B;{s3m2m zI!iVeG;T|PZw^V%?(d$J+#xR3O7>o_oq!@TmTAhE`nSxbgNP{|#WUKM1sxPUq5<(# z_b$xCm&HGN#yN{^os-MqWu(Knh(~hMvaMQQF{Oitia{uNS5xgnw1YFPfbtN7Dv1yw zVh;vw$Zzu`y-*6IS%@3TiTXUV_HI_Jwp&+D_#C~3kb#&$sqZIh?(WO8ylJWW*H&G( zB6qacc9)&IN0jOWGb9BL3!VT5WFRqF%Zzew(3?Cf9YZB2{ZAknFuku4e(gt7cEitd^FDXXS0yh3^$s;Zd+S=NNr5 z3ep;0p+4WC# zNedpAHh0;ARdFb&2oNMNnFv6V1cU%22muHHW*&w+(e53F@~)0;2U1L6q`^=htAP2_Gn> z%ZSD2P{TWcLqN8${6lb&KDz^7ID;iTvFcd%t)Wn=HSa^LsVn%x8sH@@Wu%y^nkd8J zjh&ghz>usM(!R@1BetuRXHPCv(=%p?NItw*fezK93)Len{ha)M}{+?x~g?sE34M;tY-P*(5{| zqYYbc4+(muG}eaLQtRF|i!isrJcIAB0{|c-0D09U0RenZMjVaw5K*}ZB*0AzaXmGo z7?{5KC1A|zX907CZlT)A*H6S|96e;gAec5u?GPye5Jd|z9V}~SRQic01|$4vHO%9J z`KXnJh71uXY|!tFh#x1~7_>91>w`(>-@I6bO7oNZ3)Pzatu6xTRE&m33uC3+fPh4r z-5;A!WD5Wg2Z0OP)yI=)inHvQX~W{M?|JpVkPrhi5kwyFUEr#ONx?HA# zm*YpKOy}!mA1`b1pAZ0HCkXNX3lH4bt4sm|lbc?R^Q46-2H8;>GHECPfLs2|$wb*; zITSVu5IuYde-tnZpb^8=e8anK0sR1AAc+7VphX0F1Q4Zx#tT<>Bp?8TqeiyHce+Lr z*3WNQPecA_#Jov`ifil_r+v{RI067w9VI58+E_s-pR^aUSlrY8Mt00e3YCKWfnR>U zBwa*Y8d;bv@Yr0Z{qkHPDx{H2Qb2!nH`!a$YJ=!8C062?ZY%@@eE=X3_Sia|`a7+0 zOugnhD~b6gAS2c?eXc^d1OlUTntqLD$g3nk0svvb6}3J?h??9Xn7QM@%mp%R{F1VE zGI6k^?Q+t>Kv!)+1|}@EKDz4LOlg4z6?Da{omAoiMlx#3%>jSP&?r-0o