From 18734bb184452fd784363a14165765e76b173fd0 Mon Sep 17 00:00:00 2001 From: tristan-myles Date: Tue, 20 Jan 2026 13:31:12 +0000 Subject: [PATCH 01/29] Combines evol and sub plot & sets sort to true - Seperate sub rate plot is commented out --- shared/nipah_functions.R | 4 +- src/nipah_transmission/nipah_transmission.R | 52 ++++++++++++--------- 2 files changed, 32 insertions(+), 24 deletions(-) diff --git a/shared/nipah_functions.R b/shared/nipah_functions.R index c840230d..c1f6e98e 100644 --- a/shared/nipah_functions.R +++ b/shared/nipah_functions.R @@ -132,8 +132,8 @@ forest_plot <- function(df, label, color_column, lims, text_size = 11, df$segment_alpha <- 1 if(qa_alpha!=1){ - df[df$qa_score<0.5, ]$plot_alpha <- qa_alpha - df[df$qa_score<0.5, ]$segment_alpha <- 0.65 * qa_alpha + df[df$qa_score<=0.5, ]$plot_alpha <- qa_alpha + df[df$qa_score<=0.5, ]$segment_alpha <- 0.65 * qa_alpha } cats <- length(unique(df[[color_column]])) diff --git a/src/nipah_transmission/nipah_transmission.R b/src/nipah_transmission/nipah_transmission.R index 3b098b0d..7900f3a4 100644 --- a/src/nipah_transmission/nipah_transmission.R +++ b/src/nipah_transmission/nipah_transmission.R @@ -133,36 +133,40 @@ for (i in seq_along(qa_thresh_vec)){ qa_threshold <- qa_thresh_vec[i] qa_alpha <- qa_alpha_vec[i] - p1 <- forest_plot(d1 |> filter(qa_score>qa_threshold), - expression(Evolutionary~Rate~(s/s/y ~10^{-4})), + p1 <- forest_plot(d1 |> + rbind(d2) |> + filter(qa_score>qa_threshold), + expression(Evolutionary~Substitution~Rate~(s/s/y ~10^{-4})), "genome_site", - c(-0.01,15), custom_colours = custom_colour_genome_groups, + c(-0.01,16), custom_colours = custom_colour_genome_groups, segment_show.legend=c(color=TRUE, shape=FALSE), - text_size=text_size, qa_alpha=qa_alpha) + + text_size=text_size, qa_alpha=qa_alpha, + sort=TRUE) + scale_color_manual(values=custom_colour_genome_groups, limits=names(custom_colour_genome_groups)) + scale_fill_manual(values=custom_colour_genome_groups, limits=names(custom_colour_genome_groups)) + guides(fill = guide_none(), color = guide_legend(title = "Genome type", order = 1, - override.aes = list(fill = custom_colour_genome_groups)), - shape=guide_none()) - - # Should this be segment? - p2 <- forest_plot(d2 |> filter(qa_score>qa_threshold), - expression(Substitution~Rate~(s/s/y ~10^{-4})), - "genome_site", - c(0,16), custom_colours = custom_colour_genome_groups, - text_size=text_size, qa_alpha=qa_alpha) + - guides(shape = guide_legend(title = "Parameter type", order=1), - fill = guide_none(), - color = guide_none()) + override.aes = list( + fill = custom_colour_genome_groups)), + shape=guide_legend(title = "Parameter type", order=1)) + + # p2 <- forest_plot(d2 |> filter(qa_score>qa_threshold), + # expression(Substitution~Rate~(s/s/y ~10^{-4})), + # "genome_site", + # c(0,16), custom_colours = custom_colour_genome_groups, + # text_size=text_size, qa_alpha=qa_alpha, sort=TRUE) + + # guides(shape = guide_legend(title = "Parameter type", order=1), + # fill = guide_none(), + # color = guide_none()) p3 <- forest_plot(d3 |> filter(qa_score>qa_threshold), "Overdispersion (max nr. of cases related to a case)", "population_group", c(0,35), custom_colours = custom_colour_pop_groups, - text_size=text_size, qa_alpha=qa_alpha) + + text_size=text_size, qa_alpha=qa_alpha, + sort=TRUE) + guides(color = guide_none(), shape = guide_none()) @@ -170,7 +174,7 @@ for (i in seq_along(qa_thresh_vec)){ "Primary Attack Rate (%)", "population_group", c(-0.01,3), custom_colours = custom_colour_pop_groups, - text_size=text_size, qa_alpha=qa_alpha) + + text_size=text_size, qa_alpha=qa_alpha, sort=TRUE) + guides(color = guide_none(), shape = guide_none()) @@ -178,7 +182,7 @@ for (i in seq_along(qa_thresh_vec)){ "Proportion of Symptomatic Cases (%)", "population_group", c(-5, 110), custom_colours = custom_colour_pop_groups, - text_size=text_size, qa_alpha=qa_alpha) + + text_size=text_size, qa_alpha=qa_alpha, sort=TRUE) + guides(color = guide_none(), shape = guide_none()) @@ -187,7 +191,7 @@ for (i in seq_along(qa_thresh_vec)){ "population_group", c(0, 1.5), custom_colours = custom_colour_pop_groups, segment_show.legend=c(color=TRUE, shape=FALSE), - text_size=text_size, qa_alpha=qa_alpha) + + text_size=text_size, qa_alpha=qa_alpha, sort=TRUE) + scale_color_manual(values=custom_colour_pop_groups, limits=names(custom_colour_pop_groups)) + scale_fill_manual(values=custom_colour_pop_groups, @@ -198,10 +202,14 @@ for (i in seq_along(qa_thresh_vec)){ shape=guide_none()) # Save plots - patchwork <- (p6 + p3 + p4 + p5 + p1 + p2) + + patchwork <- (p6 + p3 + p4 + p5 + p1 + guide_area()) + plot_layout(ncol = 2, widths = c(1,1), guides = "collect") - patchwork <- patchwork + plot_annotation(tag_levels = "A") + patchwork <- patchwork + + plot_annotation(tag_levels = "A", + theme = theme(legend.position = "bottom", + legend.direction = "horizontal")) + ggsave(paste0("figure_",label,"trans.png"), plot = patchwork, width = 14, height = 10) ggsave(paste0("figure_",label,"trans.pdf"), From aafb3530a2827774d4d355d658395e26bff4688d Mon Sep 17 00:00:00 2001 From: tristan-myles Date: Tue, 20 Jan 2026 14:22:59 +0000 Subject: [PATCH 02/29] Updates Nipah to use dplyr::select - Prevents issues in Nipah_bsl_data_synthesis where MASS::select() masks dplyr::select() --- shared/nipah_functions.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/shared/nipah_functions.R b/shared/nipah_functions.R index c1f6e98e..13764225 100644 --- a/shared/nipah_functions.R +++ b/shared/nipah_functions.R @@ -11,7 +11,7 @@ data_curation <- function(articles, outbreaks, models, parameters, plotting,swit mutate(new_refs = ifelse(refs %in% refs[duplicated(refs)], paste0(sub("\\)$", "", refs),letters[counter],")"), refs)) |> - select(-counter,-refs) |> + dplyr::select(-counter,-refs) |> rename(refs = new_refs) |> mutate(refs = str_to_title(refs)) @@ -70,7 +70,7 @@ data_curation <- function(articles, outbreaks, models, parameters, plotting,swit mutate(central = coalesce(parameter_value, 100*cfr_ifr_numerator/cfr_ifr_denominator, 0.5*(parameter_lower_bound+parameter_upper_bound))) |> - select(-c(no_unc)) + dplyr::select(-c(no_unc)) if (plotting) { parameters <- param4plot From a637bf0930b5619baa89485163f0c2897a3707f8 Mon Sep 17 00:00:00 2001 From: tristan-myles Date: Tue, 20 Jan 2026 14:27:05 +0000 Subject: [PATCH 03/29] Sets population group to be a factor - Ordering in the plot now matches the legend --- src/nipah_transmission/nipah_transmission.R | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/nipah_transmission/nipah_transmission.R b/src/nipah_transmission/nipah_transmission.R index 7900f3a4..334577c6 100644 --- a/src/nipah_transmission/nipah_transmission.R +++ b/src/nipah_transmission/nipah_transmission.R @@ -35,6 +35,17 @@ parameters <- dfs$parameters |> left_join(qa_scores) # *----------------------------- Data preparation -----------------------------* +# Sort population sample type to match legend where Other, Unspecified, or NA +# are included. A neater solution would be to do this during cleaning or for the +# relevant subset of population groups included in the params considered in +# this script. +parameters <- parameters |> + mutate(population_group = factor( + population_group, + levels = c(sort(setdiff(unique(population_group), + c("Other", "Unspecified"))), + "Other", "Unspecified"))) + d1 <- parameters |> filter(parameter_type == "Mutations - evolutionary rate") d2 <- parameters |> filter(parameter_type == "Mutations - substitution rate") d3 <- parameters |> filter(parameter_class == "Overdispersion") From a518cd5ceaea7521e7d1e7fcb9fc7741cf735ba4 Mon Sep 17 00:00:00 2001 From: tristan-myles Date: Tue, 20 Jan 2026 14:27:53 +0000 Subject: [PATCH 04/29] Adds a ^* to Lo Presti to flag for caption --- src/nipah_transmission/nipah_transmission.R | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/nipah_transmission/nipah_transmission.R b/src/nipah_transmission/nipah_transmission.R index 334577c6..22c5419d 100644 --- a/src/nipah_transmission/nipah_transmission.R +++ b/src/nipah_transmission/nipah_transmission.R @@ -69,6 +69,13 @@ d4 <- d4 |> mutate(across(all_of(variables_to_mutate), mutate(parameter_unit = ifelse(parameter_unit == "No units", "Percentage (%)", parameter_unit)) +# Presti reports this metric as an evolutionary rate +# Rahman refers to the metric as both a substitution and evolutionary rate in +# the same text +d1 <- d1 |> + mutate(refs=str_replace(refs, "Lo Presti \\(2016\\)", + "'Lo Presti (2016\\)'^'*'")) + d1 <- d1 |> arrange(genome_site,-central) # TODO: check Upper and lower bound - zero out for now @@ -153,6 +160,7 @@ for (i in seq_along(qa_thresh_vec)){ segment_show.legend=c(color=TRUE, shape=FALSE), text_size=text_size, qa_alpha=qa_alpha, sort=TRUE) + + scale_y_discrete(labels = function(x) parse(text = x)) + scale_color_manual(values=custom_colour_genome_groups, limits=names(custom_colour_genome_groups)) + scale_fill_manual(values=custom_colour_genome_groups, From 2a80c47efd392cbafc49d6184cfe1474ed81a2e3 Mon Sep 17 00:00:00 2001 From: tristan-myles Date: Mon, 26 Jan 2026 21:53:29 +0000 Subject: [PATCH 05/29] Adds bsl plot and linetype to legend --- src/nipah_delays/nipah_delays.R | 52 +++++++++++++++++++++++++-------- 1 file changed, 40 insertions(+), 12 deletions(-) diff --git a/src/nipah_delays/nipah_delays.R b/src/nipah_delays/nipah_delays.R index 26378465..67198b79 100644 --- a/src/nipah_delays/nipah_delays.R +++ b/src/nipah_delays/nipah_delays.R @@ -13,6 +13,10 @@ orderly_parameters(pathogen = NULL) orderly_dependency("db_cleaning", "latest(parameter:pathogen == this:pathogen)", c("articles.csv", "outbreaks.csv", "models.csv", "params.csv")) +orderly_dependency("nipah_bsl_data_synthesis", + "latest(parameter:pathogen == this:pathogen)", + c("bsl_model_plot.RDS")) + orderly_shared_resource("nipah_functions.R" = "nipah_functions.R") source("nipah_functions.R") @@ -243,6 +247,18 @@ p6_oa_o_plots <- list("all"=list(), "qa"=list()) p7_o_a_plots <- list("all"=list()) p7_oo_reduced_plots <- list("all"=list()) +# Read in BSL plot - incubation period +bsl_model_plot <- readRDS("bsl_model_plot.RDS") + +bsl_model_plot <- bsl_model_plot + + labs(title="", x="Incubation period (days)", color="Model", fill="Model") + + theme_minimal() + + theme(panel.border = element_rect(color = "black", linewidth = 1.25, + fill = NA), + text = element_text(size = text_size), + legend.position=c(0.45, 0.95), + legend.direction = "horizontal") + for (i in seq_along(qa_thresh_vec)){ label <- labels[i] qa_threshold <- qa_thresh_vec[i] @@ -417,6 +433,7 @@ for (i in seq_along(qa_thresh_vec)){ p6_oa_o_plots[[plot_type]][[colour_col]] <- p6_oa_o_plots[[plot_type]][[colour_col]] + guides(shape = guide_none(), + linetype = guide_none(), color = guide_legend(title = "Outcome")) } @@ -448,8 +465,9 @@ for (i in seq_along(qa_thresh_vec)){ p7_oo_reduced_plots[[plot_type]][[colour_col]] <- p7_oo_reduced_plots[[plot_type]][[colour_col]] + - guides(shape = guide_none(), - color = guide_legend(title = "Outcome")) + guides(shape = guide_legend(title = "Parameter type", order=1), + color = guide_legend(title = "Outcome"), + linetype = guide_legend(title = "Variation type")) } # Update legends for final plot @@ -458,7 +476,8 @@ for (i in seq_along(qa_thresh_vec)){ p1_incb_plots[[plot_type]][["population_country"]] + guides(shape = guide_none(), fill = guide_none(), - color = guide_none()) + color=guide_legend(title="Country"), + linetype = guide_none()) p4_oo_plots[[plot_type]][["population_country"]] <- p4_oo_plots[[plot_type]][["population_country"]] + @@ -466,7 +485,8 @@ for (i in seq_along(qa_thresh_vec)){ limits=all_groups) + scale_fill_manual(values=custom_colours, limits=all_groups) + - guides(shape = guide_legend(title = "Parameter type", order=1), + guides(shape = guide_none(), + linetype = guide_none(), color=guide_legend(title="Country")) } } @@ -494,8 +514,9 @@ for (i in seq_along(qa_thresh_vec)){ # Alternative is to use guides="collect" (legends) in plot_layout delays_plot <- (p1_incb_plots[[plot_type]][["population_country"]] + p4_oo_plots[[plot_type]][["population_country"]])/( - p6_oa_o_plots[[plot_type]][["parameter_type"]]) + - plot_layout(heights = c(2, 1), widths = c(2, 1)) + + bsl_model_plot + + p6_oa_o_plots[[plot_type]][["parameter_type"]]) + + plot_layout(heights = c(1, 1), widths = c(1, 1)) + plot_annotation(tag_levels = 'A') ggsave(paste0("figure_5", label,"_delays.pdf"), plot = delays_plot, @@ -508,6 +529,7 @@ for (i in seq_along(qa_thresh_vec)){ p3_ao_plots[[plot_type]][["parameter_type"]] <- p3_ao_plots[[plot_type]][["parameter_type"]] + guides(shape = guide_none(), + linetype = guide_none(), color = guide_legend(title = "Outcome")) + common_left_legend # p7_oo_reduced_plots[[plot_type]][["parameter_type"]] <- @@ -519,13 +541,18 @@ for (i in seq_along(qa_thresh_vec)){ # p6_oa_o_plots[[plot_type]][["parameter_type"]] + # plot_annotation(tag_levels = 'A') - delays_plot <- (p1_incb_plots[[plot_type]][["population_country"]] + - p4_oo_plots[[plot_type]][["population_country"]]) / ( - p7_oo_reduced_plots[[plot_type]][["parameter_type"]] + - p3_ao_plots[[plot_type]][["parameter_type"]]) + - plot_layout(heights = c(1.25, 1), widths = c(1.25, 1)) + - plot_annotation(tag_levels = 'A') + left_col <- p1_incb_plots[[plot_type]][["population_country"]] / + bsl_model_plot + + plot_layout(heights = c(17, 20)) + + right_col <- p7_oo_reduced_plots[[plot_type]][["parameter_type"]] / + p4_oo_plots[[plot_type]][["population_country"]]/ + p3_ao_plots[[plot_type]][["parameter_type"]] + + plot_layout(heights = c(14, 17, 6)) + delays_plot <- (left_col | right_col) + + plot_layout(widths = c(1, 1)) + + plot_annotation(tag_levels = "A") ggsave(paste0("figure_5", label,"_delays.pdf"), plot = delays_plot, width = 28, height = 15) @@ -535,6 +562,7 @@ for (i in seq_along(qa_thresh_vec)){ } # ============================================================================== +# *--------------------------------- Not used ---------------------------------* # Incubation facet: incubation_pc_facet <- forest_plot(d1, "Incubation period (days)", From 692bcd552a5640312ff30fc18b4e55fcedd4e8e4 Mon Sep 17 00:00:00 2001 From: tristan-myles Date: Mon, 2 Feb 2026 09:29:04 +0000 Subject: [PATCH 06/29] Updates shared bsl_synthesis & nipah bsl task - bsl sythensis: adds make posterior summary function which returns cdf draws - nipah bsl task: wraps chunks of the code so that the model doesn't need to be rerun, and adds new plots --- shared/bsl_data_synthesis.R | 17 +- .../nipah_bsl_data_synthesis.R | 454 ++++++++++++++---- 2 files changed, 384 insertions(+), 87 deletions(-) diff --git a/shared/bsl_data_synthesis.R b/shared/bsl_data_synthesis.R index 53987dab..7ed4b9eb 100644 --- a/shared/bsl_data_synthesis.R +++ b/shared/bsl_data_synthesis.R @@ -251,17 +251,28 @@ bsl_make_density_summary <- function(dist_name, post, scale <- exp(loc_d); dweibull(x_seq, shape = phi, scale = scale) } }) + dens_mat[i, ] <- rowMeans(dens_l) + # dens_mat[i, ] <- t(apply(dens_l, 1, median)) + } - data.frame( + + if (posterior_cdf){ + dx <- c(x_seq[1], diff(x_seq)) + dens_mat <- dens_mat * dx + dens_mat <- t(apply(dens_mat, 1, cumsum)) + } + + summary_df <- data.frame( x = x_seq, - mean = apply(dens_mat, 2, mean, na.rm = TRUE), + mean = apply(dens_mat, 2, median, na.rm = TRUE), low = apply(dens_mat, 2, quantile, 0.025, na.rm = TRUE), high = apply(dens_mat, 2, quantile, 0.975, na.rm = TRUE), model = dist_name ) -} + return (summary_df) +} # =============================== # AUTOMATED BSL DIAGNOSTIC REPORT diff --git a/src/nipah_bsl_data_synthesis/nipah_bsl_data_synthesis.R b/src/nipah_bsl_data_synthesis/nipah_bsl_data_synthesis.R index 517f0298..0c4778a4 100644 --- a/src/nipah_bsl_data_synthesis/nipah_bsl_data_synthesis.R +++ b/src/nipah_bsl_data_synthesis/nipah_bsl_data_synthesis.R @@ -9,18 +9,35 @@ library(ggsci) library(coda) library(bridgesampling) library(tidybayes) +library(orderly2) +library(patchwork) # *--------------------------------- Orderly ----------------------------------* -orderly_parameters(pathogen = NULL) +# Results presented in the Nipah manuscript used 10000 runs +orderly_parameters(pathogen = NULL, runs=0, actual_dataset=TRUE, + generate_diagnostics=FALSE) orderly_dependency("db_cleaning", "latest(parameter:pathogen == this:pathogen)", c("articles.csv", "outbreaks.csv", "models.csv", "params.csv")) -orderly_shared_resource("nipah_functions.R" = "nipah_functions.R") +orderly_shared_resource("nipah_functions.R" = "nipah_functions.R", + "bsl_data_synthesis.R" = "bsl_data_synthesis.R") + +artefacts <- c("nipah_incubation_period_bsl_model.pdf", + "bsl_model_plot.RDS") + +if (runs==0){ + orderly_shared_resource( + "bsl_model_fits.RDS" = file.path("nipah","bsl_model_fits.RDS")) +}else{ + artefacts <- c(artefacts, "bsl_model_fits.RDS") +} + source("nipah_functions.R") +source("bsl_data_synthesis.R") orderly_artefact(description="Nipah data synthesis figures", - c()) + artefacts) # *------------------------------ Data curation -------------------------------* articles <- read_csv("articles.csv") @@ -32,39 +49,35 @@ dfs <- curation(articles, outbreaks, models, parameters, plotting = TRUE) parameters <- dfs$parameters # *--------------------------------- BSL Hierarchical Model ----------------------------------* - -TESTING <- TRUE - # -------------------------- # Bayesian Synthetic Likelihood for Data Synthesis based on this paper # https://www.tandfonline.com/doi/full/10.1080/10618600.2017.1302882?scroll=top&needAccess=true#d1e213 -# -------------------------- - -orderly_strict_mode() -orderly_parameters(pathogen = NULL) -orderly_dependency() -orderly_shared_resource("nipah_functions.R" = "nipah_functions.R", - "bsl_data_synthesis.R" = "bsl_data_synthesis.R") -source("nipah_functions.R") -source("bsl_data_synthesis.R") - -orderly_artefact() - # -------------------------- # Incubation Period datasets # -------------------------- -if(TRUE){ - datasets <- list( - d1 = list(median = 10.0, min = 9.0, max = 12.0, n = 4), - d2 = list(median = 4.0, min = 2.0, max = 7.0, n = 6), - d3 = list(median = 9.0, min = 6.0, max = 14.0, n = 11), - d4 = list(median = 9.5, min = 4.0, max = 14.0, n = 22), - d5 = list(median = 8.0, min = 3.0, max = 20.0, n = 15), - d6 = list(median = 9.0, min = 6.0, max = 11.0, n = 11), - d7 = list(median = 10.0, min = 8.0, max = 15.0, n = 14), - d8 = list(median = 9.0, min = 6.0, max = 11.0, n = 11) +# TODO: Get the source for this in a reproducible way +if(actual_dataset){ + datasets_with_refs <- list( + d1 = list(median = 10.0, min = 9.0, max = 12.0, n = 4, + access_param_id="113_003"), + d2 = list(median = 4.0, min = 2.0, max = 7.0, n = 6, + access_param_id="113_004"), + d3 = list(median = 9.0, min = 6.0, max = 14.0, n = 11, + access_param_id="138_020"), + d4 = list(median = 9.5, min = 4.0, max = 14.0, n = 22, + access_param_id="121_002"), + d5 = list(median = 8.0, min = 3.0, max = 20.0, n = 15, + access_param_id="081_001"), + d6 = list(median = 9.0, min = 6.0, max = 11.0, n = 11, + access_param_id="030_002"), + d7 = list(median = 10.0, min = 8.0, max = 15.0, n = 14, + access_param_id="151_001"), + d8 = list(median = 9.0, min = 6.0, max = 11.0, n = 11, + access_param_id="033_006") ) + + datasets <- lapply(datasets_with_refs, function(x) x[!names(x) %in% "access_param_id"]) } else { #testing only datasets <- list( d1 = list(mean = 5.2, sd = 2.1, median = 5.0, q25 = 4.3, q75 = 9.0, n = 100), @@ -73,6 +86,7 @@ if(TRUE){ ) } + if( "min" %in% names(datasets[[1]]) ) { summary_names_used <- c("median", "min", "max") } else { @@ -93,8 +107,7 @@ stopifnot(all(is.finite(S_obs))) # -------------------------- # Safe summary function # -------------------------- -if("min" %in% names(datasets$d1)) -{ +if("min" %in% names(datasets$d1)){ summaries_fun <- function(x) { x <- x[is.finite(x)] if (length(x) < 1) x <- 1e-6 @@ -139,8 +152,7 @@ model_list <- list( weibull = bsl_create_model("weibull", theta0_weibull) ) -if( interactive() ) -{ +if(runs!=0){ # -------------------------- # Parallel setup # -------------------------- @@ -173,7 +185,7 @@ if( interactive() ) replicate(chains_per_model, { bsl( S_obs, - n = if(TESTING) 1000 else 10000, + n = runs, M = 2000, model = model_list[[mdl]], covRandWalk = diag(0.125, 3), @@ -211,13 +223,14 @@ ess_results <- lapply(mcmc_lists, effectiveSize) print(ess_results) # --- Trace plots for visual check -par(mfrow = c(3, 3)) -for (mdl in names(mcmc_lists)) { - for (p in 1:ncol(as.matrix(mcmc_lists[[mdl]][[1]]))) { - traceplot(mcmc_lists[[mdl]], varname = colnames(as.matrix(mcmc_lists[[mdl]][[1]]))[p], - main = paste("Trace:", mdl, "-", colnames(as.matrix(mcmc_lists[[mdl]][[1]]))[p])) - } -} +# Commenting out for orderly run +# par(mfrow = c(3, 3)) +# for (mdl in names(mcmc_lists)) { +# for (p in 1:ncol(as.matrix(mcmc_lists[[mdl]][[1]]))) { +# traceplot(mcmc_lists[[mdl]], varname = colnames(as.matrix(mcmc_lists[[mdl]][[1]]))[p], +# main = paste("Trace:", mdl, "-", colnames(as.matrix(mcmc_lists[[mdl]][[1]]))[p])) +# } +# } # -------------------------- # Inspect results @@ -257,13 +270,19 @@ set.seed(2025) posterior_summaries <- bsl_summarise_posteriors(posterior_samples, L = 50) print(posterior_summaries) - - # Combine for all models +# Makes a denisty summary by default density_summary <- do.call(rbind, lapply(names(posterior_samples), function(nm) - bsl_make_density_summary(nm, posterior_samples[[nm]], n_draws = 200, L = 20) + bsl_make_posterior_summary(nm, posterior_samples[[nm]], + x_seq=seq(0,30, length.out=400), + n_draws = 200, L = 20) )) +cdf_summary <- do.call(rbind, lapply(names(posterior_samples), function(nm) + bsl_make_posterior_summary(nm, posterior_samples[[nm]], + x_seq=seq(0, 40, length.out=2000), + n_draws = 200, L = 50, posterior_cdf=TRUE) +)) # -------------------------- # Posterior summaries for plotting (use predictive mean & CI) # -------------------------- @@ -294,96 +313,363 @@ posterior_summaries_90th_long <- posterior_summaries %>% # -------------------------- # Observed summaries from datasets (extract summary-list values) # -------------------------- -observed_summary <- do.call(rbind, lapply(names(datasets), function(nm) { - d <- datasets[[nm]] +observed_summary <- do.call(rbind, lapply(names(datasets_with_refs), function(nm) { + d <- datasets_with_refs[[nm]] data.frame( dataset = nm, + access_param_id = d$access_param_id, median_obs = ifelse(!is.null(d$median), d$median, NA_real_), min_obs = ifelse(!is.null(d$min), d$min, NA_real_), max_obs = ifelse(!is.null(d$max), d$max, NA_real_) ) })) # Note: observed_summary$dataset are d1,d2,... (not model names). We'll map them visually below. - +observed_summary <- observed_summary |> + left_join(parameters[, c("access_param_id", "refs")]) # -------------------------- # Plot: density + predictive mean CI + observed medians/ranges # -------------------------- +# *=================================== pdf ====================================* obs_plot_df <- observed_summary -posterior_summaries_long$y_jitter <- jitter(rep(0.02, nrow(posterior_summaries_long)), amount = 0.02) -posterior_summaries_median_long$y_jitter <- jitter(rep(0.02, nrow(posterior_summaries_median_long)), amount = 0.02) -posterior_summaries_90th_long$y_jitter <- jitter(rep(0.02, nrow(posterior_summaries_90th_long)), amount = 0.02) -obs_plot_df$y_jitter <- jitter(rep(0, nrow(obs_plot_df)), amount = 0.015) +# Observed summary row +obs_strip <- obs_plot_df %>% + transmute( + panel = "Summary statistics", + model = "Summary statistics", + y = "Summary statistics", + x = median_obs, + xmin = min_obs, + xmax = max_obs + ) -bsl_model_plot <- ggplot(density_summary, aes(x = x, y = mean, color = model, fill = model)) + - # Posterior densities - geom_ribbon(aes(ymin = low, ymax = high), alpha = 0.15, colour = NA) + - geom_line(size = 1) + +# Posterior medians (with interval) +med_strip <- posterior_summaries_median_long %>% + filter(model != "lognormal") %>% + transmute( + panel = "Median", + model = model, + y = model, + x = median_value, + xmin = low_value, + xmax = high_value + ) + +# Posterior 90th (point only; add xmin/xmax if you have them) +p90_strip <- posterior_summaries_90th_long %>% + filter(model != "lognormal") %>% + transmute( + panel = "90th", + model = model, + y = model, + x = median_value, + xmin = NA_real_, + xmax = NA_real_ + ) + +strip_df <- bind_rows(obs_strip, med_strip, p90_strip) %>% + mutate( + panel = factor(panel, + levels = c("Summary statistics", "Median", "90th"), + labels=c("Summary statistics", "Median", + "90th percentile")), + y = row_number() + ) + +# lanonc_colours <- ggsci::pal_lancet("lanonc")(9) - # Posterior medians + 95% CrIs (as points and horizontal bars above the curves) +summary_strip_plot <- + ggplot(strip_df, aes(x = x, y = y)) + + # interval rows where available geom_errorbarh( - data = posterior_summaries_median_long, - aes(y = y_jitter + 0.15, xmin = low_value, xmax = high_value, color = model), - height = 0.005, size = 0.8, inherit.aes = FALSE + aes(xmin = xmin, xmax = xmax, color = model), + height = 0.15, linewidth = 0.5, + na.rm = TRUE ) + + # points (median / observed median / 90th point) geom_point( - data = posterior_summaries_median_long, - aes(x = median_value, y = y_jitter + 0.15, fill = model, color = model), - shape = 8, size = 2, inherit.aes = FALSE + aes(color = model), + size = 2 + ) + + # facet_grid(panel ~ ., scales = "free_y", space = "free_y") + + ggforce::facet_col(facets = vars(panel), + scales = "free_y", + space = "free") + + labs(title = "", x = "Incubation period (days)", y = "", colour="Model") + + guides(colour=guide_none())+ + theme_bw(base_size = 10) + + theme( + axis.text.y = element_blank(), + axis.ticks.y = element_blank(), + panel.grid = element_blank(), + plot.margin = margin(-14, 3, 0, -10), + strip.background = element_blank(), + strip.text = element_text(color = "black"), + ) + + scale_color_manual( + values = c( + "Summary statistics" = "black", + gamma = "#ED0000FF", + weibull = "#00468BFF" + ) ) + + scale_y_discrete(expand = expansion(add = 0.75)) - geom_point( - data = posterior_summaries_90th_long, - aes(x = median_value, y = y_jitter + 0.11, fill = model, color = model), - shape = 4, size = 2, inherit.aes = FALSE +bsl_model_plot <- ggplot(density_summary |> filter(model!="lognormal"), + aes(x = x, y = mean, color = model, fill = model)) + + # Posterior densities + geom_ribbon(aes(ymin = low, ymax = high), alpha = 0.15, colour = NA) + + geom_line(linewidth = 1) + + labs( + title = "", + colour="Model", + fill="Model", + x = "Incubation period (days)", + y = "Density" ) + + theme_bw(base_size = 14) + + theme(legend.position = "right") + + scale_color_lancet() + + scale_fill_lancet() + +bsl_model_plot <- bsl_model_plot + theme(text=element_text(size=28)) + +bsl_model_plot <- bsl_model_plot + theme(legend.position = "top") + + inset_element( + summary_strip_plot, + left = 0.65, + bottom = 0.25, + right = 0.99, + top = 0.975 + ) + +saveRDS(bsl_model_plot, "bsl_pdf_model_plot.RDS") +ggsave('nipah_incubation_period_bsl_pdf_model.pdf', plot=bsl_model_plot, width = 11, height = 7) + +# *=================================== cdf ====================================* +# *-------------------------------- Main text ---------------------------------* +obs_plot_df$y_jitter <- -0.3 + seq_len(nrow(obs_plot_df)) * 0.03 - # Observed data medians + ranges (at bottom, in black) +observation_plot <- ggplot() + geom_errorbarh( data = obs_plot_df, aes(y = y_jitter - 0.01, xmin = min_obs, xmax = max_obs), height = 0.005, color = "black", size = 0.8, - inherit.aes = FALSE + inherit.aes = FALSE, linetype="dashed" ) + geom_point( data = obs_plot_df, aes(x = median_obs, y = y_jitter - 0.01), - shape = 18, fill = "white", color = "black", size = 2, + shape = 22, fill="black", color = "black", size = 3, inherit.aes = FALSE ) + + labs( + title = "", + x = "Incubation period (days)", + y = "" + ) + + theme_bw(base_size = 14) + + scale_y_continuous( + breaks = obs_plot_df$y_jitter, + labels = obs_plot_df$refs + ) + + xlim(c(0,30)) + +cdf_quantiles <- cdf_summary |> + group_by(model) |> + arrange(x, .by_group = TRUE) |> + summarise( + x_50 = approx(mean, x, xout = 0.5, ties = base::mean, rule = 2)$y, + x_90 = approx(mean, x, xout = 0.9, ties = base::mean, rule = 2)$y, + q_high_50 = approx(low, x, xout = 0.5, ties = base::mean, rule = 2)$y, + q_high_90 = approx(low, x, xout = 0.9, ties = base::mean, rule = 2)$y, + q_low_50 = approx(high, x, xout = 0.5, ties = base::mean, rule = 2)$y, + q_low_90 = approx(high, x, xout = 0.9, ties = base::mean, rule = 2)$y, + .groups = "drop" + ) +cdf_quantiles_long <- cdf_quantiles |> + pivot_longer( + cols = -model, + names_to = "quantile", + values_to = "x_value" + ) |> + mutate( + y = as.numeric(str_extract(quantile, "(?<=_)\\d+")) / 100 + ) + +cdf_quantiles_wide <- cdf_quantiles_long |> + mutate(level = sub(".*_(\\d+)$", "\\1", quantile)) |> + mutate(type = sub("_(\\d+)$", "", quantile)) |> + dplyr::select(-quantile) |> + pivot_wider( + names_from = type, + values_from = x_value) + +median_labels <- cdf_quantiles_wide |> + filter(level==50) |> + mutate(q_low=q_low-2.75) + +nineth_labels <- cdf_quantiles_wide |> + filter(level==90) |> + mutate(q_low=q_low-4) + +bsl_cdf_model_plot <- ggplot(cdf_summary |> + filter(model=="weibull"), + aes(x = x, y = mean, color = model, + fill = model)) + + # Posterior densities + geom_ribbon(aes(ymin = low, ymax = high), + alpha = 0.15, colour = NA) + + geom_line(size = 1) + + geom_errorbarh( + data=cdf_quantiles_wide |> filter(model=="weibull"), + aes(xmin = q_low, xmax = q_high, y=y), + height = 0.035) + + geom_point(data=cdf_quantiles_wide |> + filter(model=="weibull"), + aes(x = x, y=y), size = 2) + + # Vertical lines + # geom_segment(data = cdf_quantiles_long |> filter(model=="weibull"), + # aes(x = 0, xend = x_value, + # y = y, yend = y), + # color = "black", + # alpha=0.6, + # linetype = "solid", + # linewidth = 0.6 + # ) + + # geom_segment( + # data = cdf_quantiles_long |> filter(model=="weibull"), + # aes( + # x = x_value, xend = x_value, + # y = 0, yend = 0, + # ), + # color = "black", + # alpha=0.6, + # linetype = "solid", + # linewidth = 0.6 + # ) + + # geom_text( + # data = data.frame(x = 7.5, y = 0.905), + # aes(x = x, y = y), + # label = expression(90^th~percentile), + # color = "black", + # size = 5, + # inherit.aes = FALSE + # ) + + # geom_text( + # data = data.frame(x = 3.5, y = 0.505, label = "Median"), + # aes(x = x, y = y, label = label), + # color = "black", + # size = 5, + # inherit.aes = FALSE + # ) + geom_text( - data = data.frame(x = 19, y = 0.135, label = "90th\npercentiles"), - aes(x = x, y = y, label = label), + data = nineth_labels |> filter(model=="weibull"), + aes(x = q_low, y = y+0.01, label = label), + label = expression(90^th~percentile), color = "black", - size = 3.5, + size = 8, inherit.aes = FALSE ) + - geom_text( - data = data.frame(x = 10, y = 0.19, label = "Medians"), - aes(x = x, y = y, label = label), + data = median_labels |> filter(model=="weibull"), + aes(x = q_low, y = y+0.01, label = label), + label = "Median", color = "black", - size = 3.5, + size = 8, inherit.aes = FALSE ) + - scale_color_lancet() + scale_fill_lancet() + labs( - title = "Posterior Predictive Densities with Medians\n and 95% Credible Bands", - x = "Days", - y = "Density" + title = "", + x = "Incubation period (days)", + y = "Probability", + colour="Model", + fill="Model" ) + - theme_bw(base_size = 14) + - theme(legend.position = "right") + xlim(c(0, 30)) + + theme_bw(base_size = 28) + + theme(legend.position = "right") + + guides(color=guide_none(), fill=guide_none()) -ggsave('nipah_incubation_period_bsl_model.pdf', plot=bsl_model_plot, width = 11, height = 7) +bsl_main_cdf_plot <- observation_plot / free(bsl_cdf_model_plot) + + plot_layout(heights=c(1,2)) & + theme(plot.margin = margin(0, 2, 0, 2)) + +saveRDS(bsl_main_cdf_plot, "bsl_main_cdf_plot.RDS") +ggsave('bsl_main_cdf_plot.pdf', plot=bsl_main_cdf_plot, width = 8, height = 6) + +bsl_cdf_model_apx_plot <- ggplot(cdf_summary, + aes(x = x, y = mean, color = model, + fill = model)) + + # Posterior densities + geom_ribbon(aes(ymin = low, ymax = high), + alpha = 0.15, colour = NA) + + geom_line(size = 1) + + geom_errorbarh( + data=cdf_quantiles_wide, + aes(xmin = q_low, xmax = q_high, y=y), + height = 0.035 + ) + + geom_point(data=cdf_quantiles_wide, + aes(x = x, y=y), size = 2) + + geom_text( + data = nineth_labels, + aes(x = q_low-2, y = y+0.01, label = label, group=model), + label = expression(90^th~percentile), + color = "black", + size = 5.5, + inherit.aes = FALSE + ) + + geom_text( + data = median_labels, + aes(x = q_low-1, y = y+0.01, label = label, group=model), + label = "Median", + color = "black", + size = 5.5, + inherit.aes = FALSE + ) + + scale_color_lancet() + + scale_fill_lancet() + + labs( + title = "", + x = "Incubation period (days)", + y = "Probability", + colour="Model", + fill="Model" + ) + + facet_wrap(~model, ncol=1, labeller = labeller(model = tools::toTitleCase)) + + xlim(c(0, 40)) + + theme(legend.position = "right") + + theme_minimal(base_size = 20) + + theme(panel.border = element_rect(color = "black", linewidth = 1.25, fill = NA)) + + guides(color=guide_none(), fill=guide_none()) + + scale_color_manual( + values = c( + gamma = "#42B540FF", + weibull = "#00468BFF", + lognormal="#ED0000FF")) + + scale_fill_manual( + values = c( + gamma = "#42B540FF", + weibull = "#00468BFF", + lognormal="#ED0000FF")) + +bsl_apx_plot <- (observation_plot + + xlim(c(0,40))) / free(bsl_cdf_model_apx_plot) + + plot_layout(heights=c(1,3.75)) + plot_annotation(tag_levels="A") & + theme(plot.margin = margin(0, 2, 0, 2), + text=element_text(size=22)) + +saveRDS(bsl_apx_plot, "bsl_apx_cdf_plot.RDS") +ggsave('bsl_apx_cdf_plot.pdf', plot=bsl_apx_plot, width = 7, height = 12) # =============================== # AUTOMATED BSL DIAGNOSTIC REPORT # =============================== - +if (generate_diagnostics){ diagnostic_dir <- "bsl_diagnostics" if (!dir.exists(diagnostic_dir)) dir.create(diagnostic_dir) @@ -443,7 +729,7 @@ ggplot(param_plot_df, aes(x = parameter, y = median, ymin = low_2.5, ymax = high # Save the plot ggsave(file.path(diagnostic_dir, "parameter_forestplot.png"), width = 10, height = 3.5) - +} COMPARITIVE_ANALYSIS <- FALSE From 9eb56f59202acc9a4543d5debce4dcf136786370 Mon Sep 17 00:00:00 2001 From: tristan-myles Date: Mon, 2 Feb 2026 09:32:55 +0000 Subject: [PATCH 07/29] Adds bsl model RDS to shared --- shared/nipah/bsl_model_fits.RDS | Bin 0 -> 374224 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 shared/nipah/bsl_model_fits.RDS diff --git a/shared/nipah/bsl_model_fits.RDS b/shared/nipah/bsl_model_fits.RDS new file mode 100644 index 0000000000000000000000000000000000000000..7e0b073b0ff2eb746448e1409441228341bf57d0 GIT binary patch literal 374224 zcmW){XH*kd_dWcZ@li)W8AJj|8&Oeu2kE1TiuB$D1f)0VErd}7lolxgsgd4$4UnL8 z0wL0SA|(U@BtQruB!rj$dp_Q~*8OnTxo4lf_qjqzfBoP8AJ_KIuf8$Io?L*@zy7)T zy8qWbQ?bX1@?uJ2p~6$&Ge3R&vikp`kfW*FPWI{Ze>pOlB>BB)w^5a#r&t8(oFRzq z1$UD;=XYyELTg(##(#JB-USb5;6)#gkTs!V!l6MqF2UI z$w+ZO>`_3B0gL!8kr(=mp0CCO-Ca>TZ%T}xIn9by>yBM__2}7fSBIXY8e2C7hkkjy z^3*yY`bQDF_xL9+!!*Zhd-+nhB5J{%giNFxHBR%~pUu4~a0WD4UP)ye(~5VOil$mx z-wqS?T^+EKVD@Y}<6SzwI_crO-F*f&n;E zc0jjFxwmn}oi!%%B-%RQ^k#p$4OZ+sq_&3|eb%P87x3rL zd(P}%?9@CuS42`%1Nae@uU}{5+8qmdzcX}0s|e(bLFA}Fzp6S+Jbm}fZ|WGK@f?V{ zQmy2Mwl?2u2YL=O-o>-n;kn^n7bc~hz9U6eS5lrI@jCE1I253b>z_JgrkRw)6sRjZ zZ6{x`2q|Alke%nNA#xs`NJMuY!V^AY&n;-3;0&D}9X)^%SJzBwW1n(BM6Vmx)np*Z zYqL>b`Hl5t#co=Nih-Bu0D0r8^2YU+2#+St&v##z%~*MxY@B{l{3SoP^!&Uf<3??|e_Q6~01lpo znHVH)x{pDLOMtqBvLGLS`)3iQ-h^z*Ter(JY_a5`Rk+I465R4DT6^Y6rUQJ|e$AAR zirNn!wPk{h8f46g!Ngw)iE3;X^y9^g-_sA$UMC&=LSq>^iPBKm(EL2d=p!Y`bns!8 zZd_Cyo*OFquClVQ%jG#7s00;KvZR70j5Q;S>*H2}IgEbE2P_s03&ZZS#_JQCUMd+4 z89^gMWs`dJs7EDXS*+rbw3uzpA77koETKwJ1%H;BYa$jeQdf;LghB9+I?pzKwzWLn zUKp;I`N+2mSQL;?G+|LbfO^$y2V~bp#`|^jut-05UGVMHmrGG{nJshcSLd%AM^-VYt$gfkb^ z@Vzf!f22KX%c{H3(nG9I1w92HEi!V|mvMuxPP4-pTybKBqUgG-V&Xi%!WvjK45j3L zrcf4>0&J-3c3AMu3=DYnqkwwTW>fMnvAfY`hfq^%E|Es;xsR|a_w}$ygqOVmON}PF zI}b|!Xgm2UfiqWH2WX<=Id*e+lTfu2^|g{aFI9*?>E{k+^P5O97=PE0L@x& zXN4{GAg(NRSSd>LY#yvU1?Z^na0XjnI`_`s^e@|u0dtpRyklp;UE z+fMoKRndpB&WHrZ&zlK8Z}3f^QvIVV09;hB`u*BB@E=ju?2EbBg7nY?}IzlUN$9dRBmD@sN8 zI*V$Iv$ubpQA0*{xwir47U)hn$+{Y<$(sN^fHTx zU5c=V%4FjKti{x-v;J3oq(JK^@JJ!PsE3jj%0a@g>Fdg(AztXOP#?|6 zA?JD9QF)<|7T}Nf92a#_KabM~X|_46u_2 z-_vrNJcHiUntkDf%RG8$l<-C03S8wR3hp`tF1&V&^OVo|dy31!+qstLPm@eZCY1QdzHPQal$Ze$bugP1hei`XB$ z5TqI`M&gzu015lauKcinRSL5ohH|1@onomyVGFbx%*tF|HD#zRTMWMpJSAQxDY_BzYGS1-%G$Y87%pId%t zJZ98$a#MlsO@v)8TmV4J^{B#fV~U-k)NBn|dw})X_e=Zezx_Rq>7o<3`>5iHB0}2q zbB_wX2gMc7_@mHq?gThytm&7(-`#HXoQ1Wm+tAuvD&5Em9!b$}(no3tb&uM|e96q~ zJ4vSZ9XwjP6}xWq?eCKjV}K3|)7=6BSGV5e<6D63 zBIv`Jb~~N4EI=ev;3_;6x_2*KP5%87>OcGD$Fx6^uL$KR)S?991M9Lw z7s^i4?p^HHkJ^_To>8CRC(g}II=t}Nf~zk)+T2F+q}3y4mU8ZC`7EP;TUbgQoxt-t z4m>$?HoJEPZnpQ>yB{gIZ|)fJNhS4FlS$Xc3G}BAva`A^gS!eOSjIT{&8!?Wo;_ao z=BZHN>-nR=UGPymOY_&G#qTOy;GkW7jAF^D>zs%%#sz_gKbY)&p&2t0e;V$av$+%B ze#7(Cm+-N9{=YUmXbonDvoa^8x#gd!&)0SbDh1hSB8#zhj+`?&Jq>eP4L#`>_)Y1~ z+wyN)W}(x89?V^nLQ^+y7}hhNOlPy6o6`uJy$JKzA_eGBlRx)HSj^`$6E1c)72C~D z&QxCwKwr!%2Nc}T=$Whcr}O7Yrt=Ava`c21e})!QTm~|4#y}fz^+Z*RhwIY*HFVF% zxfi#EGPfQ>&U;@tLFOm^?J1;+GEVxgmHbq3%o{#6hoHQb>zuu zPJayx0|zBp)(c~tSLyDQWn^(9VIfT|HIVd|PGAZPr%{-c>tk$LvDHu@;`6)`ejAws zwb@iga75CcXUC-UM02QA)lpkUI=q91?Cl^$D@R*qiSWXi3E68|SNOYSGOZ+Kd6>tn z+J%lo-zVw)yhqsz$zIBm`weMfTXTCg?k{CC$@_aiI?*G(XBs%xyelh=iT#*#8qrq-IFp}#nX z|L0USAIpuKv~q3U8ycg*;O#RxWmbaL<%+7}dD@mM{7ujFtarj{nMJ&jrp!+Gh8Y<1 zN!Wo9`-cn4SEpyBJ#Qi$PV+^I<_MwNmn7DWEQl9@Z?D7$qUKvQf~LRRyKtiSH@KHX z8}cX;tLHVHApJ0oeq!%!%jZ$Pg&^9kdrMl^Hj(}wk8%GazYL``+^8{h9M`CvZ4}V~ zYPZ}`ETDYqlN;8}?N{=Wsa0PnwqMqhCyZBfBptm5^6T_iOh@Y|Pa*=*7&Z5JHrKR{ z;|Z>Y_`K(VgEbpj0>>^RB%uwPh#3APmU4}t#9)8fOsz+qnmgy(@quMys!LqKw?*#? zHX{_I9ZJ7qiPJ#1F@VPI$9@T>Qnjj2d_9L(PJ|IVQ5U!BUbZ~$&`3?wi}mx<(+VII@LU)!{S(J#`Lw1E-tXVxU`Pi4%X|uZdnOYyJDM1GL>+yJ4Tk zm~`@%cgC=Y8&dXsJ9oe7TMoWUwnKP%URLyBkOt%JxtXZ;+%>}Su9)fH%BYG7G1*Ky z$8^AX`88oovI|1$AuEp&PVC4LxKgY zDwg(dAJ!%0OTYOozTV-x@#MXxD(Peq3!5zs=Rj*u>;2BZY@jc43o^Q28X{}R>=XDpG0 zvT+sy5!5{IFJ;NwYLg$V8bl^lTcS`MC*3#VEc%jvpQg(8rPpUeqK$epP?YPvthu1{ z(If7>Yj1qQbq4t6m=MIhT*v!@Fdw1AVv#tmmxvIaG`jtObfy0O@twMa8xG2j_=~qD zn&%FB%&+YdM*@0zmL99=I!>%3COsrL??)0`Z z#N{L3p0DE0(%?G87w8Ecz-dBe^)CJ!tEtb@fOtNx60(tO=CqlQr)){c{+r;l{MpVk zHkx2S*R)+l(?}9L(5I`zscy5^7G~aZ3Tb2vf3GUMlKK2-nq5)FZlTChYyvq$n{g`0 zjOgbBc^2O@9g$N2xyz>Olx4QC3e8&G5w+10;GhEbTO4@xYx7gA?|Pog$T(VufR0*B z4Cyy?0+Ba`w0xO&G^l9HEgcl6B^xcvv`p|(k}#%?q7H^g6ZByq|8bFqormdBxhs1G zX*r)`oPYE5>Zm&n8a3f(N;IDLgsP!9K zMg-N9RxGB)yl?ND%|P0&zSLpNe{6E)oJiUXnrmO7*HttYkB)utC)wRjyz>~B6idzi zeCPIUQr)WJ_shiGlJtMF=QHHI({w%RerlmhQ**Wx=6Fu}oyh>HX)mRpXl2@!YNJ_c zm8XT*1O1M&6Id)28ROp9rYBWb8VnoN$Y-hvAxH}!{|FI@sVl@Zv!3fPo(3HWhlT6H z-jibcZy50)EuKgoa0@t@IJII{aayrDj9>5#=@ZKuNZ^a02``H0_}5-1FbK0D{;99} zsi+#AU79;OWRUptcLmzRD z`<+SQju>n6 z0ndI~dYsOe=EIb5{*Y;aaR6D*{5k7R@1tdHnTBSP!luUpH>64zMAprFx=3H zI}yj#5Ou6xV28Yk;3Aq-wH2yo9a;7~-`3+nlCZ<+a$}8@hgSNjb+<`Jq85INH%T;2 zB&jEBPProd*5w^FL0=GAQ)}fT)bL&ST07A^dAg&(5i-8$O>8Ba+m*-1e1IV0$XVUFI?C4`iq=WoI5S!?2}Y*;4LISxf4AJRxSt8|@kZ@I_tS^b~5 z2$LA^k1f1eL(=tsU6}qiL_PS#&9SPZzP6fAv)C{J1RlS0vbM?sSSx z9Qlfs!7qyQR}3y0NB?C%;TT#tcr@_SK=|&>I$aLw?QWz0e!fs>7=N|*RXNMpI<$7@ zqfn5@jj?v+*=4TiRuLVR4_(u0OVVF;oM@xGtRf0_<68@&6=zYc)#5t9hgTqZwjKhf znDU>M-X8G3q`XTB6b#Wa&XK`?>NA&qLf+$E<_*e%r+_nu+HP@kN>ntxvn* zJZp)xR)o1Dq3Ggvs8*SC}@i0}*G@ zFSH4DZvEnp@TtlM+pa%xRs=S(h!odZOui1+iiVcttB2i0bwzjNm42nmjVC8+D=OJy z?Wzl2Q^|i5=;N3S#7$+l!#uxUpNMmkc@g2$f12e_XXdAea<0TG$&5#}6phNwcd%hA z_5=~j+tCF%YDK%`cA|R)v0abSttUYIN!)DbXhEHhp{aC5qnWi^}6-~c`S*Tme861n> z;I~9~>%0n~ej__4{|Qz_Fc8hT_GV)6ZQbm%vzM1hpX4;UV_x$i2*`?QNz3wU@+Xl^NC?KDg z9fJBkE$&)Bw@9Hl9Zhi$-=F-Ej)Q7c21BK_msUMT4u}1>NZ%cD?^*DeMc4s^a%-plz=#WQ$&6 zp*qHYV!~bz>~%3TC?{g)M);z3?u2MC2F@&Xgsg$aA$~(9V2GTJOrY6w61PjGw6c}A zwAWeAAYK>WU=*r4=-iHuZvd-)u^CtOU!;tGeNBieik&A+FG&FWP&Z*4L3DJ4yIv^S zZG$NomI`X@z@*sPG=`H(RCIyyFxW%}SMDlZYq)1)f9aOWubeY}Q<*X58@B@5J7lhu zB}tgqZ9#ez8j0hec!15?e-=Juvx2kdsx_DzTd^1&^E)m)a(`LZ)3W#MclYR{PcMaX z&VR>wvOr%R(3skkDpbp6P`gxA-#TV;M{95Px_nQ%d24@H;XR1vNwZFiS)Fw@){+)J8Xw6>!Fx^*EkvHTYB>rmu^`tR7 zO7_^>Z?s8BBB`ffEiPx5-Qe#h#^`!?4M&w*s2H?=&hGS-Fl8AzET1As5kJ0|#k|p&vfLJ)Kns zB(<(;DDV5<%W&<}ev0HM@G>2@;}pw1|DIzVuwa}}lW@Bf7(SnU@uD}BC{2ekT2!?a z%hG94MltjT6I@8#3%`ypY!=wNd)%})?<#v{8k(}^ya5&yJSP_HzK-HJan3GFNy6Fu z3LAkQn1XsABp>WbkUNy@{te{etXon3lE*p`&JMXovM#(LxEbMoo>z4g@)*Lzo{B5)*|+WZh{eObaE7IRbgBKnZ#)gGwNHJp!* zixVFK^OgqdKO)a-VZ^PCRk@tmvX?}B?BvX{e#)amvT0R96w#4~W}iLyaXJmYQi!5^ zGmd``sa+%xCALm6)3!oC-A}jnVg40z@lB_J^R9L^6)siX9ERNfToEQ{=!+<{keJ)) z;p#X0F{PwRw0pMlw}RuZXZ$P2+P}M`&9|1;;*&^pQ*X733je8@tWo7UVGXSbEgg-B zY2dFV=8Kt{RpmeRhK-5}QMTpHZ;k%1pUS;#y_CBwGtOI6YZW?^GZ;n_f){eFyge+M z>f)!duCZZ5V?fS9Ezfg$mWvp^(J zSGiNGwo$d$|Ho07k*;`VWwf5qA8>JcGIvDZPV=SZ>B*-Y`l(S?c{XsOGT7j(9l!cYSY4|1%0 zVP@w?s8Lr}Dieu&`NE>D?bu`NAo&70-|eqlA(----W>Rr7ZPId)6+PoLtpdnWUs;H z?9l^@#`A&Bk3nrmw1Bl|v&;AkS$UH5jZ>huwRl6=-Mjoa;K#M7<#pd9_Up)Z#9TV* zjbQ`>9q~w9+aqEGe5UTU9qTrYZv573)!Mbzwl(u;?Wcj1=PnQ6ChM!k>O9mzg6)zN z?U%0vs9o}SNn9&-Fk2P516ZWTaSTeDdz!#Qr?`3v0Z z>?GVztKWIHL>k>S0ptIN>g7wNfm+I#Bx^u(i%u#%vjg4QmY+>~%xkjVHB)Tni15Dt1s@#15(krsxA|iDRV7 zycal<-*kyHQ#Jl7ebYrlpDG8me=DOLu>2O_$Y$2KQ`MZmdHW&HU zcYxT4g00ifD+|Qg{G>d8i>x%z)t_54im~b*tudi6Beh(Vdk$@#KLFQ9tk>#&!h{bbeAsFi3~DqNFo2@F21ZH zkPI~FhL6nGK~f%jK|}ugea8(nyQh1)U!L1f z<>p<F*aa&JL_efu})y zIV*~RuceL3!sF57ey6QZ6uu!1#(ZObSKA*nB!Lo#+k?goZKIMX!nLi8M zUH|WO_0--^fLpwVvKVKN)_B?9VGgKj2Gap4Xr+?GIlo4dR+t#&{F*C3 zs|AOx7{O3ZCvMJwTUK96Qmraprvik8Of=6fR6whZl;x$f9P>52OIUhV9CMe!j9sFy z8|uaD;IZ8PWVPer@iPTIbsfEdgZ>F4Y_nNVWd~Sr7`z#U+)M(aw%YP%lGg+>vNJTS zsP5qwh7+Fu4Tca-_@X8QD;`N`zOrVFG@&bEQoxbDEvy&KDTCfiWg65N;>bazS8?6@ zPgZq3x0U2!1__>QS_(#wrMycf4&aPaGgV#?n^s%bH+R`_BFQA|w3${Vz`J$qPy}yN z$X#8L!Y9y39s*PuuF?uN zG5o>y?*N~p1EhPhCC{m74Ni{9;z$m;H`!=IJa##pUXvI7Y2n4zaK0k@L`pt)jSq6NKbs687*?FAe$Zn# zLkm{=AJgWNBD=N54;>~o08o0nABmS&5ANyH*7IOF-FO@7ZP#QlP@_wbbgb}+s^#uS z)2h@q6O!G(<#F}aVasM`DuyCGLxs5sbHgoy%xtyh)7D{pDu`Xk3>u@f7^%tNf>edJ z7)*|oNq~@57vaZWU&Is`ebTdVK-&!+)KlU|hBj<3l&~*S0=0tC@lNe$9YgH5#$