diff --git a/nipah_workflow.R b/nipah_workflow.R index ea65967e..0ec2d1dd 100644 --- a/nipah_workflow.R +++ b/nipah_workflow.R @@ -42,10 +42,23 @@ orderly_run("db_cleaning",list(pathogen="NIPAH", debug_mode=TRUE)) orderly_run("nipah_latex_tables", list(pathogen="NIPAH")) # *---------------------------- Plots and analysis ----------------------------* +orderly_run("nipah_serology", list(pathogen="NIPAH")) + +# orderly_run("nipah_map", list(pathogen="NIPAH")) + orderly_run("nipah_transmission", list(pathogen="NIPAH")) orderly_run("nipah_severity", list(pathogen="NIPAH")) +orderly_run("nipah_bsl_data_synthesis", list(pathogen="NIPAH")) + +# I assume the issue below is caused by the BSL library and other packages will +# explicitly reference MASS when a function is needed +# MASS::select masks dplyr::select +# MASS::area masks patchwork::select +select <- dplyr::select +area <- patchwork::area + orderly_run("nipah_delays", list(pathogen="NIPAH")) orderly_run("nipah_summary", list(pathogen="NIPAH")) diff --git a/shared/bsl_data_synthesis.R b/shared/bsl_data_synthesis.R index 53987dab..e105e2b2 100644 --- a/shared/bsl_data_synthesis.R +++ b/shared/bsl_data_synthesis.R @@ -188,7 +188,7 @@ bsl_summarise_posteriors <- function(posterior_samples_list, L = 50) { } -# Posterior predictive draw samples +# Posterior predictive draw samples - confirm not same as below bsl_make_density_summary <- function(dist_name, post, x_seq = seq(0, 20, length.out = 300), n_draws = 200, L = 20) { @@ -231,9 +231,11 @@ bsl_make_density_summary <- function(dist_name, post, # -------------------------- # Posterior predictive density summaries (for plotting) # -------------------------- -bsl_make_density_summary <- function(dist_name, post, - x_seq = seq(0, 20, length.out = 400), - n_draws = 200, L = 20) { +bsl_make_posterior_summary <- function(dist_name, post, + x_seq = seq(0, 20, length.out = 400), + n_draws = 200, L = 20, + posterior_cdf=FALSE){ + post <- as.matrix(post) draws <- sample(1:nrow(post), min(n_draws, nrow(post))) dens_mat <- matrix(NA, nrow = length(draws), ncol = length(x_seq)) @@ -251,17 +253,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/shared/cleaned_outbreak_data.RDS b/shared/cleaned_outbreak_data.RDS new file mode 100644 index 00000000..a92124bc Binary files /dev/null and b/shared/cleaned_outbreak_data.RDS differ diff --git a/shared/nipah/bsl_model_fits.RDS b/shared/nipah/bsl_model_fits.RDS new file mode 100644 index 00000000..7e0b073b Binary files /dev/null and b/shared/nipah/bsl_model_fits.RDS differ diff --git a/shared/nipah_functions.R b/shared/nipah_functions.R index c840230d..1805ed8c 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 @@ -113,7 +113,8 @@ curation <- function(articles, outbreaks, models, parameters, plotting) { # function to produce forest plot for given dataframe forest_plot <- function(df, label, color_column, lims, text_size = 11, show_label = FALSE, custom_colours = NA, - segment_show.legend=NA, sort=FALSE, qa_alpha=1) { + segment_show.legend=NA, sort=FALSE, qa_alpha=1, + point_size=3) { stopifnot(length(unique(df$parameter_unit[!is.na(df$parameter_unit)])) == 1)#values must have same units if (sort){ @@ -132,8 +133,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]])) @@ -142,22 +143,22 @@ forest_plot <- function(df, label, color_column, lims, text_size = 11, y = urefs, yend = urefs, color = .data[[color_column]],), linewidth=3, alpha = df$segment_alpha, show.legend = segment_show.legend) + geom_errorbar(aes(xmin=parameter_uncertainty_lower_value, xmax=parameter_uncertainty_upper_value, - y = urefs), + y = urefs, linetype="Uncertainty"), width = 0.25, lwd=0.5, color = "black", alpha=df$plot_alpha) + geom_errorbar(data= df[!df$uncertainty_present,], aes(xmin=parameter_2_lower_bound, xmax=parameter_2_upper_bound, - y = urefs), - width = 0.25, lwd=0.5, color = "black", linetype="dashed", + y = urefs, linetype="Variability"), + width = 0.25, lwd=0.5, color = "black", lineend = "square", alpha=df[!(df$uncertainty_present),]$plot_alpha) + geom_errorbar(data= df[df$uncertainty_present,], aes(xmin=parameter_2_lower_bound, xmax=parameter_2_upper_bound, - y = urefs), - width = 0.25, lwd=0.5, color = "black", linetype="dashed", + y = urefs, linetype="Variability"), + width = 0.25, lwd=0.5, color = "black", lineend = "square", position = position_nudge(y=-0.25), alpha=df[df$uncertainty_present,]$plot_alpha) + geom_point(aes(x = parameter_value, y = urefs, shape = parameter_value_type, fill = .data[[color_column]]), - alpha=df$plot_alpha, size = 3, stroke = 1, color = "black") + alpha=df$plot_alpha, size = point_size, stroke = 1, color = "black") if (all(df$parameter_class=="Reproduction number")) { gg <- gg + @@ -172,6 +173,9 @@ forest_plot <- function(df, label, color_column, lims, text_size = 11, Other = 23, `Central - unspecified`=25), breaks = c("Mean", "Median", "Unspecified", "Other", "Central - unspecified")) + + scale_linetype_manual(name = "Variation Type", + values = c("Uncertainty" = "solid","Variability" = "dashed"), + breaks = c("Uncertainty", "Variability")) + scale_x_continuous(limits = lims, expand = c(0, 0)) + scale_y_discrete(labels = setNames(df$refs, df$urefs)) + labs(x = label, y = NULL) + @@ -187,6 +191,9 @@ forest_plot <- function(df, label, color_column, lims, text_size = 11, Other = 23, `Central - unspecified`=25), breaks = c("Mean", "Median", "Unspecified", "Other", "Central - unspecified")) + + scale_linetype_manual(name = "Variation Type", + values = c("Uncertainty" = "solid","Variability" = "dashed"), + breaks = c("Uncertainty", "Variability")) + scale_x_continuous(limits = lims, expand = c(0, 0)) + scale_y_discrete(labels = setNames(df$refs, df$urefs)) + labs(x = label, y = NULL) + @@ -196,9 +203,13 @@ forest_plot <- function(df, label, color_column, lims, text_size = 11, } if (cats == 1) { - gg <- gg + guides(fill = "none", color="none", shape = guide_legend(title = NULL,order = 1)) + gg <- gg + guides(fill = "none", color="none", + shape = guide_legend(title = NULL,order = 1), + linetype=guide_legend(title = NULL,order = 2)) } else { - gg <- gg + guides(fill = "none", color = guide_legend(title = NULL,order = 1), shape = guide_legend(title = NULL,order = 2))} + gg <- gg + guides(fill = "none", color = guide_legend(title = NULL,order = 1), + shape = guide_legend(title = NULL,order = 2), + linetype=guide_legend(title = NULL, order = 3))} if(show_label) gg <- gg + geom_text_repel(aes(x = coalesce(parameter_value), y = urefs, label = population_country_ISO), nudge_y = 0.5, segment.color = "grey50" ) @@ -332,7 +343,7 @@ metamean_wrap <- function(dataframe, estmeansd_method, digits = digits, digits.sd = digits, digits.weight = digits, col.diamond.lines = "black",col.diamond.common = colour, col.diamond.random = colour, weight.study = "same", col.square.lines = "black", col.square = colour, col.study = "black", col.inside = "black", - at = seq(lims[1],lims[2],by=2), xlim = lims, xlab = label, fontsize = 10, colgap.forest.left = paste0( colgap_shift,"cm")) + at = seq(lims[1],lims[2],by=2), xlim = lims, xlab = label, fontsize = 13, colgap.forest.left = paste0( colgap_shift,"cm")) dev.off() } else { mtan <- metamean(data = dataframe, @@ -357,7 +368,7 @@ metamean_wrap <- function(dataframe, estmeansd_method, digits = digits, digits.sd = digits, digits.weight = digits, col.diamond.lines = "black",col.diamond.common = colour, col.diamond.random = colour, weight.study = "same", col.square.lines = "black", col.square = colour, col.study = "black", col.inside = "black", - at = seq(lims[1],lims[2],by=2), xlim = lims, xlab = label, fontsize = 10) + at = seq(lims[1],lims[2],by=2), xlim = lims, xlab = label, fontsize = 13) dev.off() } @@ -497,7 +508,7 @@ metagen_wrap <- function(dataframe, estmeansd_method, digits = digits, digits.sd = digits, digits.weight = digits, col.diamond.lines = "black",col.diamond.common = colour, col.diamond.random = colour, weight.study = "same", col.square.lines = "black", col.square = colour, col.study = "black", col.inside = "black", - at = seq(lims[1],lims[2],by=2), xlim = lims, xlab = label, fontsize = 10) + at = seq(lims[1],lims[2],by=2), xlim = lims, xlab = label, fontsize = 11.5) dev.off() } else { mtan <- metagen(data = dataframe, @@ -525,7 +536,7 @@ metagen_wrap <- function(dataframe, estmeansd_method, digits = digits, digits.sd = digits, digits.weight = digits, col.diamond.lines = "black",col.diamond.common = colour, col.diamond.random = colour, weight.study = "same", col.square.lines = "black", col.square = colour, col.study = "black", col.inside = "black", - at = seq(lims[1],lims[2],by=2), xlim = lims, xlab = label, fontsize = 10) + at = seq(lims[1],lims[2],by=2), xlim = lims, xlab = label, fontsize = 11.5) dev.off() } @@ -562,15 +573,24 @@ metaprop_wrap <- function(dataframe, subgroup, method.tau = "ML") png(file = "temp.png", width = width, height = height, res = resolution) + par(mar = c(2, 2, 2, 1)) forest(mtan, layout = "RevMan5", overall = plot_pooled, pooled.events = TRUE, print.subgroup.name = FALSE, sort.subgroup = sort_by_subg, study.results = plot_study, digits = digits, - col.diamond.lines = "black",col.diamond.common = colour, col.diamond.random = colour, - col.subgroup = "black", col.inside = "black", - weight.study = "same", #col.square.lines = "green", col.square = "blue", #not working - at = at, xlim = xlim, xlab="Case Fatality Ratio", fontsize=11) + col.diamond.lines = "black",col.diamond.common = colour, + col.diamond.random = colour, + col.square = colour, col.square.lines = "black", + col.study = "black", col.subgroup = "black", + col.inside = "black", weight.study = "same", + at = at, xlim = xlim, xlab="Case Fatality Ratio", + fs.predict.labels = 11.5, + fs.hetstat=11, + fs.test.subgroup = 11, + fs.axis = 11, + fontsize = 14, + plotwidth = "72.5mm") dev.off() } else { mtan <- metaprop(data = dataframe, @@ -586,10 +606,17 @@ metaprop_wrap <- function(dataframe, subgroup, overall = plot_pooled, pooled.events = TRUE, study.results = plot_study, digits = digits, - col.diamond.lines = "black",col.diamond.common = colour, col.diamond.random = colour, - col.subgroup = "black", col.inside = "black", - weight.study = "same", #col.square.lines = "green", col.square = "blue", #not working - at = at, xlim = xlim, xlab="Case Fatality Ratio", fontsize=11) + col.diamond.lines = "black",col.diamond.common = colour, + col.diamond.random = colour, + col.square = colour, col.square.lines = "black", + col.subgroup = "black", col.inside = "black", weight.study = "same", + at = at, xlim = xlim, xlab="Case Fatality Ratio", + fs.predict.labels = 11.5, + fs.hetstat=11, + fs.test.subgroup = 11, + fs.axis = 11, + fontsize = 14, + plotwidth = "72.5mm") dev.off() } diff --git a/src/db_cleaning/nipah/nipah_cleaning.R b/src/db_cleaning/nipah/nipah_cleaning.R index f0927552..9e69dea9 100644 --- a/src/db_cleaning/nipah/nipah_cleaning.R +++ b/src/db_cleaning/nipah/nipah_cleaning.R @@ -279,6 +279,12 @@ param_cleaning <- function(df){ paste0("Proportion asymptomatic is reported in the paper.", "Denominator is contacts of Nipah patients who gave blood specimen") + # Sero uncert issue: Reported as 8.35 in the paper but this is clearly a typo; + # other estimates in the table have the same estimate but a 95% CI upper bound + # of 38.35 + sero_2892_sero_uncert_filter <- df$access_param_id=="008_003" + df[sero_2892_sero_uncert_filter, "parameter_uncertainty_upper_value"] <- 38.35 + # CovID: 3057 delay_3057_incp_filter <- df$access_param_id=="151_001" df[delay_3057_incp_filter, "parameter_2_value_type"] <- "Range (paired)" @@ -702,32 +708,69 @@ param_cleaning <- function(df){ df[df$access_param_id=="093_003", "parameter_unit"] <- "Percentage (%)" # Split 2931 incubation row - hd_2931_row_filter <- df$access_param_id=="138_016" - df[hd_2931_row_filter, "distribution_type"] <- NA - df[hd_2931_row_filter, "distribution_par1_type"] <- NA - df[hd_2931_row_filter, "distribution_par1_value"] <- NA - df[hd_2931_row_filter, "distribution_par1_uncertainty"] <- NA - df[hd_2931_row_filter, "distribution_par2_type"] <- NA - df[hd_2931_row_filter, "distribution_par2_value"] <- NA - df[hd_2931_row_filter, "distribution_par2_uncertainty"] <- NA - - new_2931_row <- df[hd_2931_row_filter, ] - new_2931_row$parameter_data_id <- generate_new_id(df, "parameter_data_id", 10) + hd_2931_incp_row_filter <- df$access_param_id=="138_016" + + new_2931_incp_row <- df[hd_2931_incp_row_filter, ] + new_2931_incp_row$parameter_data_id <- generate_new_id( + df, "parameter_data_id", 10) + # No corresponding redcap entry so make an ID + new_2931_incp_row$access_param_id <- "138_3141" + new_2931_incp_row$parameter_value <- 9.7 + new_2931_incp_row$parameter_value_type <- "Mean" + + new_2931_incp_row$parameter_statistical_approach <- "Estimated model parameter" + new_2931_incp_row$parameter_paired <- "No" + new_2931_incp_row$parameter_2_unit <- NA + new_2931_incp_row$method_2_from_supplement <- NA + new_2931_incp_row$parameter_2_statistical_approach <- NA + + new_2931_incp_row$parameter_2_value_type <- NA + new_2931_incp_row$parameter_2_lower_bound <- NA + new_2931_incp_row$parameter_2_upper_bound <- NA + + # Remove dist param values + df[hd_2931_incp_row_filter, "distribution_type"] <- NA + df[hd_2931_incp_row_filter, "distribution_par1_type"] <- NA + df[hd_2931_incp_row_filter, "distribution_par1_value"] <- NA + df[hd_2931_incp_row_filter, "distribution_par1_uncertainty"] <- NA + df[hd_2931_incp_row_filter, "distribution_par2_type"] <- NA + df[hd_2931_incp_row_filter, "distribution_par2_value"] <- NA + df[hd_2931_incp_row_filter, "distribution_par2_uncertainty"] <- NA + + df <- rbind(df, new_2931_incp_row) + + # Split 2931 serial interval row + hd_2931_si_row_filter <- df$access_param_id=="138_015" + + new_2931_si_row <- df[hd_2931_si_row_filter, ] + new_2931_si_row$parameter_data_id <- generate_new_id( + df, "parameter_data_id", 10) + # No corresponding redcap entry so make an ID - new_2931_row$access_param_id <- "138_3141" - new_2931_row$parameter_value <- 9.7 - new_2931_row$parameter_value_type <- "Mean" - - new_2931_row$parameter_statistical_approach <- "Estimated model parameter" - new_2931_row$parameter_paired <- "No" - new_2931_row$parameter_2_unit <- NA - new_2931_row$method_2_from_supplement <- NA - new_2931_row$parameter_2_statistical_approach <- NA - - new_2931_row$parameter_2_value_type <- NA - new_2931_row$parameter_2_lower_bound <- NA - new_2931_row$parameter_2_upper_bound <- NA - df <- rbind(df, new_2931_row) + new_2931_si_row$access_param_id <- "138_2718" + new_2931_si_row$parameter_value <- 13 + new_2931_si_row$parameter_value_type <- "Median" + + new_2931_si_row$parameter_statistical_approach <- "Estimated model parameter" + new_2931_si_row$parameter_paired <- "No" + new_2931_si_row$parameter_2_unit <- NA + new_2931_si_row$method_2_from_supplement <- NA + new_2931_si_row$parameter_2_statistical_approach <- NA + + new_2931_si_row$parameter_2_value_type <- NA + new_2931_si_row$parameter_2_lower_bound <- NA + new_2931_si_row$parameter_2_upper_bound <- NA + + # Remove dist param values + df[hd_2931_si_row_filter, "distribution_type"] <- NA + df[hd_2931_si_row_filter, "distribution_par1_type"] <- NA + df[hd_2931_si_row_filter, "distribution_par1_value"] <- NA + df[hd_2931_si_row_filter, "distribution_par1_uncertainty"] <- NA + df[hd_2931_si_row_filter, "distribution_par2_type"] <- NA + df[hd_2931_si_row_filter, "distribution_par2_value"] <- NA + df[hd_2931_si_row_filter, "distribution_par2_uncertainty"] <- NA + + df <- rbind(df, new_2931_si_row) # Labels for IQR and Range are different for variability so copying from # uncertainty results in different labels diff --git a/src/db_compilation/redcap_compilation.R b/src/db_compilation/redcap_compilation.R index 907d0cbe..6d22408a 100644 --- a/src/db_compilation/redcap_compilation.R +++ b/src/db_compilation/redcap_compilation.R @@ -200,7 +200,7 @@ orderly_dependency( # Manually fixed files and "cleaning" script - these need to be in the # src/db_compilation folder -orderly_resource(fixing_files) +orderly_resource(setNames(fixing_files, fixing_files)) ## Outputs orderly_artefact( diff --git a/src/nipah_bsl_data_synthesis/bsl_data_synthesis.R b/src/nipah_bsl_data_synthesis/bsl_data_synthesis.R deleted file mode 100644 index 3bd2ea49..00000000 --- a/src/nipah_bsl_data_synthesis/bsl_data_synthesis.R +++ /dev/null @@ -1,405 +0,0 @@ -library(BSL) -library(parallel) -library(doParallel) -library(MASS) -library(foreach) -library(ggplot2) -library(tidyverse) -library(ggsci) -library(coda) -library(bridgesampling) - -# -------------------------- -# 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 -# -------------------------- - - -# -------------------------- -# Dataset simulator -# -------------------------- -bsl_simulate_dataset <- function(dist, loc, phi, n) { - x <- switch(dist, - lognormal = rlnorm(n, meanlog = loc, sdlog = phi), - gamma = { mean_d <- exp(loc); shape <- phi; rate <- shape / mean_d; rgamma(n, shape = shape, rate = rate) }, - weibull = { scale <- exp(loc); shape <- phi; rweibull(n, shape = shape, scale = scale) }, - stop("Unknown dist")) - x[!is.finite(x)] <- 1e-6 - x -} - -# -------------------------- -# Simulator returning matrix -# -------------------------- -bsl_make_simulator_matrix <- function(dist, datasets, summaries_fun, summary_names) { - function(theta, M = 1) { - mu0 <- theta[1]; tau <- exp(theta[2]); phi <- exp(theta[3]) - - sims <- matrix(NA, nrow = M, ncol = length(summary_names) * length(datasets)) - colnames(sims) <- rep(summary_names, length(datasets)) - - for (m in 1:M) { - out_list <- lapply(datasets, function(d) { - n <- d$n - loc_d <- rnorm(1, mu0, tau) - x <- bsl_simulate_dataset(dist, loc_d, phi, n) - s <- summaries_fun(x) - s_vec <- numeric(length(summary_names)) - names(s_vec) <- summary_names - s_vec[names(s)] <- s - s_vec - }) - sims[m, ] <- unlist(out_list) - } - - sims[!is.finite(sims)] <- 1e-6 - sims - } -} - -# -------------------------- -# Prior function -# -------------------------- -bsl_fnPrior_hier <- function(theta) { - mu0 <- theta[1]; log_tau <- theta[2]; log_phi <- theta[3] - dnorm(mu0, 0, 5, log = TRUE) + - dnorm(log_tau, log(0.5), 1.5, log = TRUE) + - dnorm(log_phi, 0, 1.5, log = TRUE) -} - - -# -------------------------- -# Create MODEL objects -# -------------------------- -bsl_create_model <- function(dist, theta0) { - newModel( - fnSim = function(theta, M = 1) bsl_make_simulator_matrix(dist, datasets, summaries_fun, summary_names_used)(theta, M), - fnSum = summaries_fun, - fnLogPrior = bsl_fnPrior_hier, - theta0 = theta0, - test = FALSE - ) -} - - -# -------------------------- -# Bridge estimator for model comparison -# -------------------------- -bsl_bridge_estimate <- function(loglik) { - m <- max(loglik) - m + log(mean(exp(loglik - m))) -} - - - -# Convert fits_multi into coda::mcmc.list for each model -bsl_to_mcmc_list_multi <- function(fits_multi) { - lapply(fits_multi, function(chain_list) { - chains <- lapply(chain_list, function(f) { - as.mcmc(as.matrix(f@theta)) - }) - mcmc.list(chains) - }) -} - -# -------------------------- -# posterior samples for further analysis -# -------------------------- -# Helper function to compute mean and 90th percentile for each posterior draw -# corrected get_summary_stats for hierarchical theta = (mu0, log_tau, log_phi) -bsl_get_summary_stats <- function(dist_name, post, L = 50) { - # post: data.frame or matrix of posterior draws (columns in order theta1, theta2, theta3) - post <- as.matrix(post) - n_draws <- nrow(post) - - # storage: predictive mean and q90 per posterior draw - pred_mean <- numeric(n_draws) - pred_median <- numeric(n_draws) - pred_q90 <- numeric(n_draws) - - for (i in seq_len(n_draws)) { - theta <- post[i, ] - mu0 <- as.numeric(theta[1]) - tau <- exp(as.numeric(theta[2])) - phi <- exp(as.numeric(theta[3])) # phi is distribution-specific parameter in your simulator - - # draw L study-level locations (loc_d) to integrate over the hierarchical random effect - locs <- rnorm(L, mean = mu0, sd = tau) - - if (dist_name == "lognormal") { - # simulator: rlnorm(n, meanlog = loc_d, sdlog = phi) - # mean of lognormal with (meanlog = m, sdlog = s): exp(m + s^2/2) - means <- exp(locs + phi^2 / 2) - median <- qlnorm(0.5, meanlog = locs, sdlog = phi) - q90s <- qlnorm(0.9, meanlog = locs, sdlog = phi) - - } else if (dist_name == "gamma") { - # simulator used: mean_d <- exp(loc); shape <- phi; rate <- shape/mean_d - # hence mean = mean_d, where mean_d = exp(loc) - # 90th percentile: qgamma(0.9, shape=phi, rate = phi/mean_d) - mean_d <- exp(locs) - means <- mean_d - # shape = phi, rate = phi / mean_d - median <- qgamma(0.5, shape = phi, rate = phi / mean_d) - q90s <- qgamma(0.9, shape = phi, rate = phi / mean_d) - - } else if (dist_name == "weibull") { - # simulator used: scale <- exp(loc); shape <- phi - # mean of Weibull: scale * gamma(1 + 1/shape) - scales <- exp(locs) - means <- scales * gamma(1 + 1/phi) - median <- qweibull(0.5, shape = phi, scale = scales) - q90s <- qweibull(0.9, shape = phi, scale = scales) - - } else { - stop("Unknown distribution name in get_summary_stats") - } - - # integrate over locs by averaging - pred_mean[i] <- mean(means, na.rm = TRUE) - pred_q90[i] <- mean(q90s, na.rm = TRUE) - pred_median[i] <- mean(median, na.rm = TRUE) - } - - # return data.frame of predictive summaries per posterior draw - data.frame(mean = pred_mean, q90 = pred_q90, median = pred_median) -} - -# wrapper to compute posterior point estimate and 95% CrI -bsl_summarise_posteriors <- function(posterior_samples_list, L = 50) { - out <- lapply(names(posterior_samples_list), function(model_name) { - post <- posterior_samples_list[[model_name]] - stats <- bsl_get_summary_stats(model_name, post, L = L) - - data.frame( - model = model_name, - mean_mean = mean(stats$mean), - mean_low = quantile(stats$mean, 0.025), - mean_high = quantile(stats$mean, 0.975), - median_mean = mean(stats$median), - median_low = quantile(stats$median, 0.025), - median_high = quantile(stats$median, 0.975), - q90_mean = mean(stats$q90), - q90_low = quantile(stats$q90, 0.025), - q90_high = quantile(stats$q90, 0.975) - ) - }) - do.call(rbind, out) -} - - -# Posterior predictive draw samples -bsl_make_density_summary <- function(dist_name, post, - x_seq = seq(0, 20, length.out = 300), - n_draws = 200, L = 20) { - post <- as.matrix(post) - draws <- sample(1:nrow(post), min(n_draws, nrow(post))) - dens_mat <- matrix(NA, nrow = length(draws), ncol = length(x_seq)) - - for (i in seq_along(draws)) { - theta <- post[draws[i], ] - mu0 <- theta[1]; tau <- exp(theta[2]); phi <- exp(theta[3]) - - # integrate over L loc_d draws - locs <- rnorm(L, mu0, tau) - dens_l <- sapply(locs, function(loc_d) { - if (dist_name == "lognormal") { - dlnorm(x_seq, meanlog = loc_d, sdlog = phi) - } else if (dist_name == "gamma") { - mean_d <- exp(loc_d) - shape <- phi - rate <- shape / mean_d - dgamma(x_seq, shape = shape, rate = rate) - } else if (dist_name == "weibull") { - scale <- exp(loc_d) - dweibull(x_seq, shape = phi, scale = scale) - } - }) - dens_mat[i, ] <- rowMeans(dens_l) # average over L loc_d draws - } - - data.frame( - x = x_seq, - mean = apply(dens_mat, 2, mean, 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 - ) -} - - -# -------------------------- -# Posterior predictive density summaries (for plotting) -# -------------------------- -bsl_make_density_summary <- function(dist_name, post, - x_seq = seq(0, 20, length.out = 400), - n_draws = 200, L = 20) { - post <- as.matrix(post) - draws <- sample(1:nrow(post), min(n_draws, nrow(post))) - dens_mat <- matrix(NA, nrow = length(draws), ncol = length(x_seq)) - for (i in seq_along(draws)) { - theta <- post[draws[i], ] - mu0 <- theta[1]; tau <- exp(theta[2]); phi <- exp(theta[3]) - locs <- rnorm(L, mu0, tau) - dens_l <- sapply(locs, function(loc_d) { - if (dist_name == "lognormal") { - dlnorm(x_seq, meanlog = loc_d, sdlog = phi) - } else if (dist_name == "gamma") { - mean_d <- exp(loc_d); shape <- phi; rate <- shape / mean_d - dgamma(x_seq, shape = shape, rate = rate) - } else if (dist_name == "weibull") { - scale <- exp(loc_d); dweibull(x_seq, shape = phi, scale = scale) - } - }) - dens_mat[i, ] <- rowMeans(dens_l) - } - data.frame( - x = x_seq, - mean = apply(dens_mat, 2, mean, 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 - ) -} - - -# =============================== -# AUTOMATED BSL DIAGNOSTIC REPORT -# =============================== - -bsl_to_mcmc <- function(fit) { - if (!is.null(fit@theta)) as.mcmc(as.matrix(fit@theta)) - else stop("fit@theta not found in object.") -} - -bsl_combine_chains <- function(fit_list) { - mcmc.list(lapply(fit_list, bsl_to_mcmc)) -} - -bsl_run_diagnostics <- function(fits_multi, diagnostic_dir = "bsl_diagnostics") { - results <- list() - - for (model_name in names(fits_multi)) { - message("\n🔍 Checking model: ", model_name) - model_dir <- file.path(diagnostic_dir, model_name) - if (!dir.exists(model_dir)) dir.create(model_dir) - - fit_list <- fits_multi[[model_name]] - mcmc_obj <- bsl_combine_chains(fit_list) - - # === TRACEPLOT === - png(file.path(model_dir, paste0("traceplot_", model_name, ".png")), - width = 1800, height = 1000, res = 150) - traceplot(mcmc_obj, main = paste("Traceplots for", model_name)) - dev.off() - - # === AUTOCORRELATION PLOTS === - png(file.path(model_dir, paste0("acf_", model_name, ".png")), - width = 1600, height = 800, res = 150) - par(mfrow = c(1, min(3, ncol(as.matrix(mcmc_obj[[1]])))), mar = c(4,4,2,1)) - for (p in seq_len(ncol(as.matrix(mcmc_obj[[1]])))) { - acf(as.numeric(as.matrix(mcmc_obj[[1]])[, p]), - main = paste("ACF -", model_name, "param", p)) - } - dev.off() - - # === DIAGNOSTICS === - gel <- gelman.diag(mcmc_obj) - ess <- effectiveSize(mcmc_obj) - gew <- geweke.diag(mcmc_obj) - - df <- data.frame( - parameter = names(ess), - Rhat = gel$psrf[, 1], - ESS = as.numeric(ess), - Geweke_Z = sapply(gew, function(x) x$z) - ) - df$model <- model_name - write.csv(df, file.path(model_dir, paste0("diagnostics_", model_name, ".csv")), row.names = FALSE) - - results[[model_name]] <- df - } - - # Combine results across models - all_results <- bind_rows(results) - write.csv(all_results, file.path(diagnostic_dir, "all_diagnostics_summary.csv"), row.names = FALSE) - - # === 4. Plots across models === - - # Effective sample size - p1 <- ggplot(all_results, aes(x = parameter, y = ESS, fill = model)) + - geom_bar(stat = "identity", position = position_dodge()) + - geom_hline(yintercept = 200, linetype = "dashed", color = "red") + - theme_bw() + - labs(title = "Effective Sample Size (ESS)", y = "ESS per parameter") - - # R-hat - p2 <- ggplot(all_results, aes(x = parameter, y = Rhat, color = model, group = model)) + - geom_point(size = 3, position = position_jitter(width = 0.1, height = 0)) + - geom_hline(yintercept = 1.05, linetype = "dashed", color = "red") + - theme_bw() + - labs(title = "Gelman–Rubin R̂ Diagnostic", y = "R̂ value") - - ggsave(file.path(diagnostic_dir, "effective_sample_size.png"), p1, width = 9, height = 6) - ggsave(file.path(diagnostic_dir, "gelman_rhat.png"), p2, width = 9, height = 6) - - message("\n✅ Diagnostics complete! Files saved in: ", normalizePath(diagnostic_dir)) - return(invisible(all_results)) -} - - -bsl_make_trace_df <- function(fits_multi) { - trace_dfs <- lapply(names(fits_multi), function(mdl) { - chain_list <- fits_multi[[mdl]] - bind_rows(lapply(seq_along(chain_list), function(ch) { - df <- as.data.frame(chain_list[[ch]]@theta) - colnames(df) <- paste0("theta", seq_len(ncol(df))) - df$iter <- seq_len(nrow(df)) - df$chain <- paste0("chain_", ch) - df$model <- mdl - df - })) - }) - bind_rows(trace_dfs) -} - - -bsl_summarise_parameters <- function(posterior_samples_list) { - require(dplyr) - out_list <- lapply(names(posterior_samples_list), function(model_name) { - post <- as.matrix(posterior_samples_list[[model_name]]) - if (ncol(post) < 3) stop("Posterior matrix must have 3 columns: mu0, log_tau, log_phi") - - mu0_draws <- post[, 1] - tau_draws <- exp(post[, 2]) # on natural scale - phi_draws <- exp(post[, 3]) # on natural scale - - df <- data.frame( - model = model_name, - parameter = c("mu0", "tau", "phi"), - mean = c(mean(mu0_draws, na.rm = TRUE), - mean(tau_draws, na.rm = TRUE), - mean(phi_draws, na.rm = TRUE)), - median = c(median(mu0_draws, na.rm = TRUE), - median(tau_draws, na.rm = TRUE), - median(phi_draws, na.rm = TRUE)), - low_2.5 = c(quantile(mu0_draws, 0.025, na.rm = TRUE), - quantile(tau_draws, 0.025, na.rm = TRUE), - quantile(phi_draws, 0.025, na.rm = TRUE)), - high_97.5 = c(quantile(mu0_draws, 0.975, na.rm = TRUE), - quantile(tau_draws, 0.975, na.rm = TRUE), - quantile(phi_draws, 0.975, na.rm = TRUE)) - ) - - # Add readable formatted CI string - df <- df %>% - mutate( - ci_95 = paste0(round(low_2.5, 3), " — ", round(high_97.5, 3)), - mean = round(mean, 3), - median = round(median, 3) - ) %>% - select(model, parameter, mean, median, ci_95, low_2.5, high_97.5) - df - }) - do.call(rbind, out_list) -} 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 diff --git a/src/nipah_delays/nipah_delays.R b/src/nipah_delays/nipah_delays.R index 26378465..c1aa0d61 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_main_cdf_plot.RDS")) + orderly_shared_resource("nipah_functions.R" = "nipah_functions.R") source("nipah_functions.R") @@ -65,12 +69,6 @@ qa_scores <- articles |> dplyr::select(covidence_id,qa_score) parameters <- dfs$parameters |> left_join(qa_scores) -parameters <- parameters |> - filter(parameter_class=="Human delay") - -parameters |> - write_csv("delay_temp_check.csv") - # Symptom Onset/Fever -> Recovery/Death parameters[parameters$access_param_id=="121_003", "parameter_type"] <- "Human delay - symptom onset>recovery/death" @@ -88,8 +86,6 @@ parameters <- parameters |> mutate(parameter_type = str_replace(parameter_type, "Human delay - ", ""), parameter_type = str_to_sentence(parameter_type)) -# Favour uncertainty over variability: - # *--------------------------------- Summary ----------------------------------* num_delays <- NROW(parameters) @@ -112,8 +108,6 @@ parameters |> select(parameter_type, other_delay_start, other_delay_end) |> print() - - # Eight rows with variability only varb_only_rows <- parameters |> filter(is.na(parameter_value) & @@ -184,12 +178,12 @@ d6 <- d3 |> # Convert to factors d3 <- d3 |> mutate(parameter_type=factor(parameter_type, - levels=c("Admission to care>death", + levels=c("Time in care (length of stay)", "Admission to care>discharge/recovery", - "Time in care (length of stay)"), - labels=c("Death", + "Admission to care>death"), + labels=c("Time in care", "Discharge/recovery", - "Time in care"))) + "Death"))) d4 <- d4 |> mutate(parameter_type=factor(parameter_type, @@ -239,10 +233,16 @@ p1_incb_plots <- list("all"=list(), "qa"=list()) p2_oa_plots <- list("all"=list(), "qa"=list()) p3_ao_plots <- list("all"=list(), "qa"=list()) p4_oo_plots <- list("all"=list(), "qa"=list()) +p5_si_plots <- list("all"=list(), "qa"=list()) 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_main_cdf_plot <- readRDS("bsl_main_cdf_plot.RDS") +bsl_main_cdf_plot <- bsl_main_cdf_plot & + theme(text = element_text(size = 28)) + for (i in seq_along(qa_thresh_vec)){ label <- labels[i] qa_threshold <- qa_thresh_vec[i] @@ -319,6 +319,26 @@ for (i in seq_along(qa_thresh_vec)){ colour_col_label, ".pdf")), plot = p3_ao_plots[[plot_type]][[colour_col]], width = 11, height = 15) + }else if(plot_type=="all"){ + p3_ao_plots[[plot_type]][[colour_col]] <- forest_plot( + d3 |> filter(qa_score>qa_threshold), 'Hospitalisation-to-outcome (days)', + colour_col, c(0,45), text_size = text_size, sort=TRUE, + custom_colours = custom_colours, qa_alpha=qa_alpha) + + ggforce::facet_col(facets = vars(parameter_type), + scales = "free_y", + space = "free") + + ggsave(file.path(plot_type, + paste0("figure_5", label, "_admis_outcome_facet_", + colour_col_label, ".pdf")), + plot = p7_oo_reduced_plots[[plot_type]][[colour_col]], + width = 15, height = 15) + + p7_oo_reduced_plots[[plot_type]][[colour_col]] <- + p7_oo_reduced_plots[[plot_type]][[colour_col]] + + guides(shape = guide_legend(title = "Parameter type", order=1), + color = guide_legend(title = "Outcome"), + linetype = guide_legend(title = "Variation type")) } # Potential typo in paper so no upper bound for the range is recorded, set to @@ -386,6 +406,12 @@ for (i in seq_along(qa_thresh_vec)){ plot = p4_oo_plots[[plot_type]][[colour_col]] , width = 15, height = 15) + # Serial interval + p5_si_plots[[plot_type]][[colour_col]] <- forest_plot( + d5, 'Serial interval (days)', colour_col, c(0,22), + text_size = text_size, sort=TRUE, custom_colours = custom_colours, + qa_alpha=qa_alpha) + if (colour_col== "parameter_type"){ # Do we want consistent colours across the SI and main plot? # If so, remove the filter @@ -417,6 +443,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")) } @@ -439,7 +466,8 @@ for (i in seq_along(qa_thresh_vec)){ d7 |> filter(qa_score>qa_threshold), 'Symptom onset-to-outcome (days)', colour_col, c(0,85), text_size = text_size, sort=TRUE, - custom_colours = custom_colours, qa_alpha=qa_alpha) + custom_colours = custom_colours, qa_alpha=qa_alpha) + + facet_wrap(~parameter_type, ncol=1, scales="free_x") ggsave(file.path(plot_type, paste0("figure_5", label, "_onset_outcome_reduced_", colour_col_label, ".pdf")), @@ -448,8 +476,36 @@ 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")) + }else if(plot_type=="all"){ + # Note:: removing two estimates + d4_filtered <- d4 |> + filter(!(parameter_type %in% c("Recovery/death"))) + + # update x-lim to 85 if including the above + p7_oo_reduced_plots[[plot_type]][[colour_col]] <- forest_plot( + d4_filtered, 'Symptom onset-to-outcome (days)', + colour_col, c(0,85), text_size = text_size, sort=TRUE, + custom_colours = custom_colours, qa_alpha=qa_alpha) + + ggforce::facet_col(facets = vars(parameter_type), + scales = "free_y", + space = "free") + + theme(strip.text.y = element_text(angle=0)) + + ggsave(file.path(plot_type, + paste0("figure_5", label, "_onset_outcome_facet_", + colour_col_label, ".pdf")), + plot = p7_oo_reduced_plots[[plot_type]][[colour_col]], + width = 15, height = 15) + + p7_oo_reduced_plots[[plot_type]][[colour_col]] <- + p7_oo_reduced_plots[[plot_type]][[colour_col]] + + 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 +514,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,10 +523,10 @@ 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")) } - } common_left_legend <- theme( legend.position = "right", @@ -494,8 +551,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, @@ -503,38 +561,67 @@ for (i in seq_along(qa_thresh_vec)){ ggsave(paste0("figure_5", label,"_delays.png"), plot = delays_plot, width = 25, height = 13) }else{ - # If stacking plots in a single row with three cols, need to adjust the - # legend to have the same size plot area? - p3_ao_plots[[plot_type]][["parameter_type"]] <- - p3_ao_plots[[plot_type]][["parameter_type"]] + - guides(shape = guide_none(), - color = guide_legend(title = "Outcome")) + common_left_legend - # p7_oo_reduced_plots[[plot_type]][["parameter_type"]] <- - # p7_oo_reduced_plots[[plot_type]][["parameter_type"]] + common_left_legend + p1_incb <- p1_incb_plots[[plot_type]][["population_country"]] + + guides(shape = guide_legend(title = "Parameter type", order=1), + color = guide_legend(title = "Outcome"), + linetype = guide_legend(title = "Variation type")) + + theme(legend.position = c(0.8,0.375)) + p3_ao <- p3_ao_plots[[plot_type]][["population_country"]] + + guides(shape = guide_none(), + linetype = guide_none(), + color = guide_none()) - # 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_annotation(tag_levels = 'A') + p5_si <- p5_si_plots[[plot_type]][["population_country"]] + + guides(shape = guide_none(), + linetype = guide_none(), + color = guide_none()) - 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') + p7_oo <- p7_oo_reduced_plots[[plot_type]][["population_country"]] + + guides(shape = guide_none(), + linetype = guide_none(), + color = guide_none()) +# +# left_col <- p1_incb / p5_si / free(bsl_model_plot) + +# plot_layout(heights = c(21, 2.75, 20)) + +# plot_annotation(tag_levels = "A") +# +# right_col <- p7_oo / p3_ao + +# plot_layout(heights = c(31, 8)) +# delays_plot <- (left_col | right_col) + +# plot_layout(widths = c(1, 1)) + + + + # p1_incb / p5_si / free(bsl_model_plot) | p7_oo / p3_ao + design <- " + A#D + B#D + B#D + B#E + C#E" + + delays_plot <- + p1_incb + free(bsl_main_cdf_plot) + p5_si + p7_oo + p3_ao + + plot_layout( + design = design, + widths = c(1, 0.05, 1), + heights = c(21, 8.25, 5, 5.25, 2.75) # makes D = 21+3+7 taller + ) + plot_annotation( + # tag_levels="A" + tag_levels = list(c("A", "B", "", "C", "D", "E")) + ) & theme(plot.tag.position = c(0, 1), + plot.tag = element_text(size = 30)) ggsave(paste0("figure_5", label,"_delays.pdf"), plot = delays_plot, - width = 28, height = 15) + width = 26, height = 20) ggsave(paste0("figure_5", label,"_delays.png"), plot = delays_plot, - width = 28, height = 15) + width = 26, height = 20) } } - # ============================================================================== +# *--------------------------------- Not used ---------------------------------* # Incubation facet: incubation_pc_facet <- forest_plot(d1, "Incubation period (days)", @@ -545,16 +632,5 @@ incubation_pc_facet <- forest_plot(d1, facet_wrap("population_country", scales="free") -# Excluding: QA score too low anyway -d4 |> - filter(covidence_id%in% c(2979, 4057)) |> - select(parameter_type, qa_score) - -# Only a single serial interval so plot not saved, but the extracted row contains -# both variability and uncertainty (derived from the gamma dist - not directly -# stated in the paper -p5_aout <- forest_plot(d5, - 'Serial interval (days)', - 'parameter_type', - c(0,20), - text_size = text_size) +ggsave(paste0("figure_5_incubation_facet.pdf"), plot = incubation_pc_facet, + width = 15, height = 9) diff --git a/src/nipah_latex_tables/nipah_latex_tables.R b/src/nipah_latex_tables/nipah_latex_tables.R index e43f953e..8d32dea1 100644 --- a/src/nipah_latex_tables/nipah_latex_tables.R +++ b/src/nipah_latex_tables/nipah_latex_tables.R @@ -920,6 +920,10 @@ sero_params <- parameters |> population_sample_type, population_group, parameter_notes, refs, central, id, access_param_id) +sero_params <- sero_params |> + mutate(in_CSF = case_when(str_detect(parameter_notes,'CSF')~TRUE, + TRUE~FALSE)) + sero_params[, c("parameter_value", "unc_type")] <- sapply( c("parameter_value", "unc_type"), function(col) gsub(" \\\\%", "", sero_params[[col]])) @@ -938,6 +942,9 @@ sero_params <- sero_params |> TRUE~ str_replace_all(parameter_value, "(\\d+(?:\\.\\d+)?)", function(x) sprintf("%.1f", as.numeric(x)))), + parameter_value=ifelse(in_CSF==TRUE, + paste0(parameter_value, "$^+$"), + parameter_value), # CI to 1 decimal unc_type=str_replace_all(unc_type, "(\\d+(?:\\.\\d+)?)", diff --git a/src/nipah_map/nipah_map.R b/src/nipah_map/nipah_map.R index ad3795cf..89abd601 100644 --- a/src/nipah_map/nipah_map.R +++ b/src/nipah_map/nipah_map.R @@ -226,6 +226,8 @@ subcolumns_outbreak <- outbreaks %>% dplyr::select(outbreak_country,outbreak_loc total_cases = total_cases/num_loc, deaths = deaths / num_loc) # for location where we split report the AVERAGE number of cases +saveRDS(subcolumns_outbreak, "cleaned_outbreak_data.RDS") + location_agg <- subcolumns_outbreak %>% group_by(outbreak_location) %>% summarise(tot_cases = sum(total_cases), tot_deaths = sum(deaths)) %>% @@ -371,7 +373,7 @@ gg_bangladesh <- ggplot() + geom_sf(data = om, lwd = 0.001, col = "lightgrey", fill = 'lightblue',alpha=0.3) + geom_sf(data = l0_in, lwd = 0.5, col = "black", fill = NA) + coord_sf(xlim = c(87.9, 93), ylim = c(21, 26.7), expand = FALSE) + - geom_point(data = manual_coords, aes(x = long, y = lat, color = Surveillance.period), size = 2, shape = 18) + + geom_point(data = manual_coords, aes(x = long, y = lat, color = Surveillance.period), size = 3.5, shape = 18) + geom_text(data = manual_coords, aes(x = long, y = lat, label = name), hjust = 0, vjust = 0, nudge_y = 0.05, size = 2.5) + scale_fill_gradient(low='palegreen', high="darkblue", trans = scales::trans_new("log2p", @@ -381,7 +383,8 @@ gg_bangladesh <- ggplot() + name = 'Total cases',na.value="grey95",limits=c(0,240)) + ggsci::scale_color_lancet() + xlab('') + ylab('') + labs(title='Bangladesh IEDCR Surveillance Data',color = "Surveillance Period") + - theme_bw() + theme_bw() + + guides(color = guide_legend(override.aes = list(size = 5))) # North India + Bangladesh from outbreak data gg_northern_india_bangladesh <- ggplot() + @@ -397,7 +400,8 @@ gg_northern_india_bangladesh <- ggplot() + breaks = c(0, 1, 4, 8, 16, 64, 240), name = 'Total cases',na.value="grey95",limits=c(0,240)) + labs(title='Northern India & Bangladesh\nreported outbreaks' ) + - theme_bw() + theme_bw() + guides(fill = guide_none(), + color=guide_none()) # Kerala gg_kerala <- ggplot() + @@ -413,7 +417,8 @@ gg_kerala <- ggplot() + breaks = c(0, 1, 4, 8, 16, 64, 240), name = 'Total cases',na.value="grey95",limits=c(0,240)) + labs(title='Kerala\nreported outbreaks' ) + - theme_bw() + theme_bw() + guides(fill = guide_none(), + color=guide_none()) @@ -432,7 +437,8 @@ gg_malaysia_singapore <- ggplot() + breaks = c(0, 1, 4, 8, 16, 64, 240), name = 'Total cases',na.value="grey95",limits=c(0,240)) + labs(title='Malaysian peninsular\nreported outbreaks' ) + - theme_bw() + theme_bw() + guides(fill = guide_none(), + color=guide_none()) # Philippines @@ -449,18 +455,72 @@ gg_philippines <- ggplot() + breaks = c(0, 1, 4, 8, 16, 64, 240), name = 'Total cases',na.value="grey95",limits=c(0,240)) + labs(title='Philippines\nreported outbreaks' ) + - theme_bw() - + theme_bw() + guides(fill = guide_none(), + color=guide_none()) layout_design <- - "AAAEE - AAAEE - AAAEE - AAAEE - AAAEE - #BCD# - #BCD#" + "AAABBBFFF + AAABBBFFF + AAABBBFFF + AAABBBFFF + AAABBBFFF + CCDDEEFFF + CCDDEEFFF" + +sero_studies <- parameters %>% + filter(parameter_class == 'Seroprevalence') %>% + mutate(population_group = factor( + population_group, + levels = c(sort(setdiff(unique(population_group), + c("Other", "Unspecified"))), + "Other", "Unspecified"))) |> + filter(qa_score>0.5) %>% + mutate(parameter_unit = 'Percentage (%)', #replace_na(parameter_unit, 'Percentage (%)'), + population_group = replace_na(population_group, 'Other'), + sero_CSF = case_when(str_detect(parameter_notes,'CSF') ~ 'CSF', + TRUE ~ 'Other')) |> + filter(parameter_type=="Seroprevalence - IgG") |> + mutate(parameter_value = coalesce(parameter_value,central)) + +sero_forest_pop_group <- forest_plot( + sero_studies, sort=TRUE, 'Serology (%)','population_country', c(-4,104), + text_size = 13) + + ggforce::facet_col(facets = vars(population_group), + scales = "free_y", + space = "free") + + guides(shape = guide_legend(title="Parameter type"), + linetype = guide_none(), + color=guide_legend(title="Country")) + theme( + legend.position = c(0.78, 0.325) + ) + +map_theme <- theme( + plot.margin = margin(2,2,2,2), + legend.position = "none" +) + +text_size = 13 +gg_bangladesh <- gg_bangladesh+ + theme(legend.title = element_text(size = 14), + legend.text = element_text(size = 12)) + +map_plot <- gg_northern_india_bangladesh + gg_bangladesh + gg_kerala + + gg_malaysia_singapore + gg_philippines + sero_forest_pop_group + + plot_layout(design = layout_design) + + plot_annotation(tag_levels = 'A') + + plot_layout(byrow = FALSE) -map_plot <- gg_northern_india_bangladesh + gg_kerala + gg_malaysia_singapore + gg_philippines + gg_bangladesh + - plot_layout(guides = 'collect', design = layout_design) + plot_annotation(tag_levels = 'A') ggsave("nipha_outbreaks_map.png", plot = map_plot, width = 17, height = 13) + +ggsave("sero_forest_pop_group.png", + plot = sero_forest_pop_group, + width = 6, height = 13) + + +ggsave("northern_india_bangladesh.png", + plot = gg_northern_india_bangladesh, + width = 10.4, height = 10) +ggsave("bangladesh.png", plot = gg_bangladesh, width = 10, height = 10) +ggsave("kerala.png", plot = gg_kerala, width = 6.4, height = 8) +ggsave("malaysia_singapore.png", plot = gg_malaysia_singapore, width = 8.4, height = 10) +ggsave("philippines.png", plot = gg_philippines, width = 9.9, height = 10) diff --git a/src/nipah_serology/nipah_serology.R b/src/nipah_serology/nipah_serology.R index 4fbfbc50..147b7c04 100644 --- a/src/nipah_serology/nipah_serology.R +++ b/src/nipah_serology/nipah_serology.R @@ -1,27 +1,15 @@ -#################### -## SEROPREVALENCE ## -#################### - -library(orderly2) -library(tidyverse) -library(stringr) -library(metafor) -library(meta) -library(estmeansd) -library(mixdist) +# *============================== Nipah Serology ==============================* +library(dplyr) library(ggplot2) library(ggsci) -library(sf) -library(ragg) -library(ggspatial) -library(ggforce) -library(png) library(grid) +library(orderly2) library(patchwork) -library(gridExtra) +library(readr) +library(stringr) +library(tidyr) -# *============================================================================* -# *------------------------------ Orderly config ------------------------------* +# *--------------------------------- Orderly ----------------------------------* orderly_strict_mode() orderly_parameters(pathogen = NULL) @@ -32,17 +20,10 @@ orderly_dependency( c("articles.csv", "models.csv", "params.csv", "outbreaks.csv")) orderly_shared_resource("nipah_functions.R"="nipah_functions.R") -orderly_shared_resource("nipah_bangladesh_district_data.csv"="nipah_bangladesh_district_data.csv") -orderly_shared_resource("NIPAH_Bangladesh_IEDCR.csv"="NIPAH_Bangladesh_IEDCR.csv") source("nipah_functions.R") -orderly_artefact("nipah-specific figures",c("figure_2.png","figure_2.pdf")) - -################### -## DATA CURATION ## -################### - +# *------------------------------ Data curation -------------------------------* articles <- read_csv("articles.csv") outbreaks <- read_csv("outbreaks.csv") models <- read_csv("models.csv") @@ -51,342 +32,91 @@ parameters <- read_csv("params.csv") dfs <- curation(articles,outbreaks,models,parameters, plotting = TRUE) articles <- dfs$articles - articles <- epireview::assign_qa_score(articles = articles)$articles -qa_scores <- articles %>% dplyr::select(covidence_id,qa_score) - -outbreaks <- dfs$outbreaks -models <- dfs$models -parameters <- dfs$parameters %>% left_join(qa_scores) %>% - mutate(article_label = make.unique(refs)) %>% - mutate(article_label = factor(article_label,levels=rev(unique(article_label)))) %>% - mutate(in_CSF = case_when(str_detect(parameter_notes,'CSF')~TRUE, - TRUE~FALSE), - population_location = coalesce(population_location, population_country)) - -nipah_bangladesh <-read_csv('nipah_bangladesh_district_data.csv') - -# *============================= LOCATIONS ===============================================* -parameter_locs <- str_to_title(unique(str_trim(na.omit(unlist(lapply(parameters$population_location, FUN=function(x) unlist(str_split(x,';')))))))) -outbreak_locs <- str_to_title(unique(str_trim(na.omit(unlist(lapply(outbreaks$outbreak_location, FUN=function(x) unlist(str_split(x,';')))))))) -locations <- unique(c(parameter_locs, outbreak_locs)) - -#prepare shapefiles for maps -l0_in <- read_sf('../../shared/World_Bank_Official_Boundaries_adm0/WB_GAD_ADM0.shp') %>% #this is the shapefile with country boundaries - rename(COUNTRY = NAM_0) #store country names in column COUNTRY -# - -# -l1_in <- read_sf('../../shared/World_Bank_Official_Boundaries_adm1/WB_GAD_ADM1.shp') %>% #this is the shapefile with level 1 regions - rename(COUNTRY = NAM_0) %>% - mutate(COUNTRY = case_when( #country names must be consistent between shapefiles - COUNTRY == "Cabo Verde" ~ "Cape Verde", - COUNTRY == "Democratic Republic of Congo" ~ "Democratic Republic of the Congo", - COUNTRY == "Guinea Bissau" ~ "Guinea-Bissau", - TRUE ~ COUNTRY)) %>% - rename(REG_CODE = ADM1CD_c) #store region codes, e.g. SL01, in column REG_CODE - -l2_in <- read_sf('../../shared/World_Bank_Official_Boundaries_adm2/WB_GAD_ADM2.shp') %>% #this is the shapefile with level 1 regions - rename(COUNTRY = NAM_0) - -om <- read_sf('../../shared/World_Bank_Official_Boundaries_Ocean_Mask/WB_GAD_ocean_mask.shp') - -mapped_locations <- locations[(locations %in% l2$NAM_2)|(locations %in% l2$NAM_1)|(locations %in% l2$COUNTRY)] -unmapped_locations <- locations[!(locations %in% l2$NAM_2) & !(locations %in% l2$NAM_1) & !(locations %in% l2$COUNTRY)] - - -location_mapping <- tribble( - ~location, ~iso3, ~district, ~division_or_state, ~notes, - # ---- Bangladesh ---- - "Haripur Upazila", "BGD", "Thakurgaon", "Rangpur", "Subdistrict of Thakurgaon", - "Haripur Upazila (Subdistrict) Of Thakurgaon District", "BGD", "Thakurgaon", "Rangpur", "Full location specification", - "Sadar Upazila", "BGD", NA, NA, "Requires specific district context", - "Northwest Bangladesh", "BGD", NA, "Rangpur; Rajshahi", "Regional designation", - "Manikgonj", "BGD", "Manikganj", "Dhaka", NA, - "Rangpur District", "BGD", "Rangpur", "Rangpur", NA, - "Rajshahi District", "BGD", "Rajshahi", "Rajshahi", NA, - "Lalmohirhat", "BGD", "Lalmonirhat", "Rangpur", "Common misspelling of Lalmonirhat", - "Comilla", "BGD", "Cumilla", "Chittagong", "Alternate spelling", - "Joypurhat", "BGD", "Joypurhat", "Rajshahi", NA, - "Bogra", "BGD", "Bogra", "Rajshahi", NA, - "Jessore", "BGD", "Jashore", "Khulna", "Official spelling update", - "Goalando", "BGD", "Rajbari", "Dhaka", "Subdistrict of Rajbari", - "7 Other Northwestern Districts", "BGD", NA, "Rangpur; Rajshahi", "Regional grouping", - "Unspecified 7 Districts", "BGD", NA, NA, "Northwest region unspecified", - "Barishal Division", "BGD", NA, "Barishal", "Division-level entry", - "Rajshahi Division", "BGD", NA, "Rajshahi", "Division-level entry", - "Chattogram Division", "BGD", NA, "Chittagong", "Division-level entry", - "Rangpur Division", "BGD", NA, "Rangpur", "Division-level entry", - "Dhaka Division", "BGD", NA, "Dhaka", "Division-level entry", - "Khulna Division", "BGD", NA, "Khulna", "Division-level entry", - "Mymensingh Division", "BGD", NA, "Mymensingh", "Division-level entry", - - # ---- India ---- - "Siliguri", "IND", "Darjeeling", "West Bengal", NA, - "Nearby Districts Of Kozhikode", "IND", "Kozhikode", "Kerala", "Primary reference district", - "Kerela", "IND", NA, "Kerala", "Common misspelling", - "West Bangal", "IND", NA, "West Bengal", "Common misspelling", - - # ---- Malaysia ---- - "Seremban Hospital", "MYS", "Seremban", "Negeri Sembilan", NA, - "Kuala Lumpur Hospital", "MYS", "Kuala Lumpur", "Kuala Lumpur", NA, - "Ipoh Hospital", "MYS", "Kinta", "Perak", "District of Kinta", - "University Of Malaya Medical Center", "MYS", "Kuala Lumpur", "Kuala Lumpur", NA, - "University Malaya Medical Centre", "MYS", "Kuala Lumpur", "Kuala Lumpur", "Alternate spelling", - "Kelang Hospita", "MYS", "Klang", "Selangor", "Misspelling of Klang", - "Tioman Island", "MYS", "Rompin", "Pahang", NA, - "Kampung Sungai Nipah", "MYS", "Seremban", "Negeri Sembilan", NA, - "Malaysia", "MYS", NA, NA, "Country-level entry", - "Peninsular Malaysia", "MYS", NA, NA, "Regional designation", - "Negri Sembilan State", "MYS", NA, "Negeri Sembilan", "Common alternate spelling", - "University Hospital", "MYS", "Kuala Lumpur", "Kuala Lumpur", "Generic reference", - "Bukit Pelandok", "MYS", "Port Dickson", "Negeri Sembilan", NA, - "Negeri Sembalin", "MYS", NA, "Negeri Sembilan", "Misspelling", - "Fatimah Hospital Ipoh", "MYS", "Kinta", "Perak", NA, - "Sibu Hospital Sarawak", "MYS", "Sibu", "Sarawak", NA, - - # ---- Other Countries ---- - "Villages Across The South Of Cameroon", "CMR", "Océan", "South", "Approximate district", - "Senator Ninoy Aquino", "PHL", "Sultan Kudarat", "Sultan Kudarat", NA, - "Wat Luang", "THA", NA, NA, "Needs precise district info", - - # ---- Empty Entry ---- - "", NA, NA, NA, "Blank entry" -) - -# Remove empty row if needed -location_mapping <- location_mapping %>% filter(location != "") - -location_mapping <- location_mapping %>% - mutate(across(where(is.list), ~map_chr(., ~paste(., collapse = "; ")))) %>% - separate_rows(district, sep = ";") %>% - mutate(district=str_trim(district)) - -location_mapping[!((location_mapping$district %in% l2_in$NAM_2)|(location_mapping$district %in% l2_in$NAM_1)| - location_mapping$division_or_state %in% l2_in$NAM_1),] - -# *=============== DEDUPLICATE OUTBREAKS =============================================================* -outbreaks <- outbreaks %>% mutate(outbreak_location = coalesce(outbreak_location, outbreak_country), - outbreak_end_year = coalesce(outbreak_end_year,outbreak_start_year), - outbreak_duriation_years = outbreak_end_year - outbreak_start_year, - type_cases_sex_disagg = replace_na(type_cases_sex_disagg,'Unspecified'), - outbreak_notes = replace_na(outbreak_notes,''), - outbreak_location_type = str_to_title(replace_na(outbreak_location_type,'Unspecified'))) -outbreaks$is_duplicate <- FALSE -outbreaks[outbreaks$outbreak_duriation_years>5,]$is_duplicate <- TRUE -outbreaks[outbreaks$outbreak_country == 'Malaysia,Singapore',]$is_duplicate <- TRUE -outbreaks[outbreaks$outbreak_country == 'Malaysia' & outbreaks$outbreak_location_type != 'State',]$is_duplicate <- TRUE #we take the individual state outbreaks rather than the national aggregate one. -outbreaks[outbreaks$outbreak_country == 'Singapore' & (is.na(outbreaks$deaths) | outbreaks$type_cases_sex_disagg=='Confirmed'),]$is_duplicate <- TRUE #we take the report with confirmed cases, also not the sero study -outbreaks[str_detect(outbreaks$outbreak_notes,'cluster'),]$is_duplicate <- TRUE #we take the individual outbreaks rather than the cluster - -exclude_ids <- c( - # India outbreaks to exclude - "1150|India|Siliguri|2001|1", - "292|India|Siliguri; West Bangal|2001|1", - - # Bangladesh outbreaks to exclude - "1150|Bangladesh|Meherpur|2001|4", - "851|Bangladesh|Meherpur|2001|4", - "1150|Bangladesh|Naogaon|2003|1", - "851|Bangladesh|Naogaon|2003|1", - "192|Bangladesh|Rajbari|2004|1", - "851|Bangladesh|Rajbari; 7 Other Northwestern Districts|2004|1", - "851|Bangladesh|Faridpur|2004|2", - "1150|Bangladesh|Faridpur|2004|4", - "186|Bangladesh|Rajbari; Faridpur|2004|NA", - "186|Bangladesh|Unspecified 7 Districts|2004|NA", - "1150|Bangladesh|Tangail|2005|1", - "960|Bangladesh|Haripur Upazila (Subdistrict) Of Thakurgaon District|2007|1", - "976|Bangladesh|Sadar Upazila|2007|3", - "1150|India|Nadia|2007|4", - "1110|Bangladesh|Manikgonj|2008|2", - "1110|Bangladesh|Rajbari|2008|2", - "947|Bangladesh|Faridpur|2010|1", - "956|Bangladesh|Lalmonirhat; Dinajpur; Rajbari; Rangpur|2010|12", - "2709|Bangladesh|Chattogram Division|2011|NA", - "2709|Bangladesh|Mymensingh Division|2013|NA", - - # Other outbreaks to exclude - "2886|India|Kerala|2018|5", - "2979|India|Kozhikode|2018|5", - "4358|India|Kerela|2018|5" -) - -outbreaks <- outbreaks %>% filter( !is_duplicate ) %>% - mutate( unique_id = paste(covidence_id, outbreak_country, outbreak_location, - outbreak_start_year, str_remove(outbreak_start_month, "^0+"), sep = "|"), - EXCLUDE = ifelse(unique_id %in% exclude_ids, 1, 0)) %>% - filter(!EXCLUDE) - -subcolumns_outbreak <- outbreaks %>% dplyr::select(outbreak_country,outbreak_location,outbreak_source,cases_confirmed,cases_suspected,cases_asymptomatic,cases_unspecified,outbreak_probable,deaths,outbreak_start_month,outbreak_start_year,refs) %>% - mutate( cases_confirmed_raw = cases_confirmed) %>% - mutate( cases_confirmed = coalesce(cases_confirmed,cases_suspected)) %>% # only use cases suspected if we don't have confirmed cases - mutate( cases_confirmed = case_when(cases_confirmed% - mutate_if(is.numeric, list(~replace_na(., 0))) %>% - mutate( total_cases = cases_confirmed + cases_asymptomatic + cases_unspecified + outbreak_probable) %>% - arrange(outbreak_start_year,outbreak_start_month) %>% - mutate(num_loc = str_count(outbreak_location, ";") + 1 ) %>% - separate_rows(outbreak_location, sep = ";") %>% - mutate(outbreak_location=str_trim(outbreak_location), - total_cases = total_cases/num_loc, - deaths = deaths / num_loc) # for location where we split report the AVERAGE number of cases +qa_scores <- articles |> dplyr::select(covidence_id,qa_score) +parameters <- dfs$parameters |> + left_join(qa_scores) - -# *============================================================================* -#extract data to be plotted -sero_studies <- parameters %>% filter(parameter_class == 'Seroprevalence') %>% - filter(qa_score>0.5) %>% - mutate(parameter_unit = 'Percentage (%)', #replace_na(parameter_unit, 'Percentage (%)'), - population_group = replace_na(population_group, 'Other'), - sero_CSF = case_when(str_detect(parameter_notes,'CSF') ~ 'CSF', - TRUE ~ 'Other')) - -# adjusted CFR -# overlap of outbreak location with sero data -# time to seroconversion, anti-body duration... https://www.neurology-asia.org/articles/20001_023.pdf -# also IgG still positive 10yrs post illness https://www.researchgate.net/profile/Victor-Chong-3/publication/287575780_Ten_year_clinical_and_serological_outcomes_of_Nipah_virus_infection/links/56f11b8c08ae519284fbd500/Ten-year-clinical-and-serological-outcomes-of-Nipah-virus-infection.pdf - -cfr_bangladesh <- read_csv('NIPAH_Bangladesh_IEDCR.csv') - -cfr <- parameters %>% filter(parameter_class == 'Severity') %>% - filter(qa_score>0.5) %>% - mutate(parameter_unit = 'Percentage (%)', #replace_na(parameter_unit, 'Percentage (%)'), - population_group = replace_na(population_group, 'Other')) - -cfr <- cfr %>% mutate(parameter_context_location_type=replace_na(parameter_context_location_type,'Unspecified'), - parameter_notes=replace_na(parameter_notes,'')) - -# DONT DEDUPLICATE AS THIS HAS CONTEXT INFORMATION -# cfr$is_duplicate <- FALSE -# #cfr[cfr$outbreak_duriation_years>5,]$is_duplicate <- TRUE -# cfr[cfr$population_country == 'Malaysia,Singapore',]$is_duplicate <- TRUE -# cfr[cfr$population_country == 'Malaysia' & cfr$parameter_context_location_type != 'State',]$is_duplicate <- TRUE #we take the individual state outbreaks rather than the national aggregate one. -# cfr[str_detect(cfr$parameter_notes,'cluster'),]$is_duplicate <- TRUE #we take the individual outbreaks rather than the cluster -# -# cfr <- cfr %>% filter( !is_duplicate ) %>% -# mutate( unique_id = paste(covidence_id, population_country, population_location, -# population_study_start_year, str_remove(population_study_start_month, "^0+"), sep = "|"), -# EXCLUDE = ifelse(unique_id %in% exclude_ids, 1, 0)) %>% -# filter(!EXCLUDE) - -cfr_from_outbreaks <- subcolumns_outbreak %>% - mutate(outbreak_source = replace_na(outbreak_source, 'Unknown')) %>% - mutate(cfr_ifr_denominator = round(total_cases), - cfr_ifr_numerator = round(deaths), - CFR = cfr_ifr_numerator / cfr_ifr_denominator, - refs = paste(outbreak_country, outbreak_location, sep = " |> "), - parameter_value = CFR, - parameter_unit = 'Percentage') %>% arrange(desc(CFR)) %>% - mutate(outbreak_source = case_when(str_detect(outbreak_source, 'Domestic animal' ) ~ 'Domestic animal', - str_detect(outbreak_source, 'Wild animal' ) ~ 'Wild animal', - TRUE ~ outbreak_source)) - -cfr_outbreak_ma <- metaprop_wrap(cfr_from_outbreaks, subgroup = 'outbreak_country', - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = 'red', - #plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = 'red', - width = 5000, height = 4750, resolution = 500) - -cfr_cfr_studies_group <- metaprop_wrap(cfr, subgroup = 'population_group', - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = 'red', - width = 4000, height = 4750, resolution = 500) - -cfr_cfr_studies_country <- metaprop_wrap(cfr, subgroup = 'population_country', - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = 'red', - width = 4000, height = 4750, resolution = 500) - -cfr_from_bangladesh_surveillance <- cfr_bangladesh %>% - filter(Cases!=0) %>% - mutate(cfr_ifr_denominator = Cases, - cfr_ifr_numerator = Death, - CFR = cfr_ifr_numerator / cfr_ifr_denominator, - refs = Year, - parameter_value = CFR, - parameter_unit = 'Percentage') %>% arrange(desc(CFR)) - -cfr_from_bangladesh_surveillance_ma <- metaprop_wrap(cfr_from_bangladesh_surveillance, subgroup = NA, - plot_pooled = TRUE, sort_by_subg = FALSE, plot_study = TRUE, digits = 3, colour = 'red', - width = 5500, height = 3750, resolution = 500) - -cfr_forest_country <- forest_plot(cfr |> mutate(parameter_value = coalesce(parameter_value,central)) |> - arrange(population_country,desc(parameter_value)), - 'Serology (%)','population_country',c(0,100),text_size = 13) - -cfr_forest_pop_group <- forest_plot(cfr |> mutate(parameter_value = coalesce(parameter_value,central)) |> - arrange(population_group,desc(parameter_value)), - 'Serology (%)','population_group',c(0,100),text_size = 13) - -# Serological studies -sero_study_type <- metaprop_wrap(sero_studies, subgroup = 'population_study_type', - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = 'red', - width = 4000, height = 4750, resolution = 500, - at = seq(0,0.3,by=0.05), xlim = c(0,0.35)) - -sero_group <- metaprop_wrap(sero_studies, subgroup = 'population_group', - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = 'red', - width = 4000, height = 4750, resolution = 500, - at = seq(0,0.3,by=0.05), xlim = c(0,0.35)) - -sero_sample_type <- metaprop_wrap(sero_studies, subgroup = 'population_sample_type', - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = 'red', - width = 4000, height = 4750, resolution = 500, - at = seq(0,0.3,by=0.05), xlim = c(0,0.35)) - -sero_country <- metaprop_wrap(sero_studies, subgroup = 'population_country', - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = 'red', - width = 4000, height = 4750, resolution = 500, - at = seq(0,0.3,by=0.05), xlim = c(0,0.35)) - -sero_CSF <- metaprop_wrap(sero_studies, subgroup = 'sero_CSF', - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = 'red', - width = 4000, height = 4750, resolution = 500, - at = seq(0,0.3,by=0.05), xlim = c(0,0.35)) - -sero_sero_type <- metaprop_wrap(sero_studies, subgroup = 'parameter_type', - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = 'red', - width = 4000, height = 4750, resolution = 500, - at = seq(0,0.3,by=0.05), xlim = c(0,0.35)) - -sero_forest_type <- forest_plot(sero_studies |> mutate(parameter_value = coalesce(parameter_value,central)) |> - arrange(parameter_type,desc(parameter_value)), - 'Serology (%)','parameter_type',c(0,100),text_size = 13) - -sero_forest_pop_type <- forest_plot(sero_studies |> mutate(parameter_value = coalesce(parameter_value,central)) |> - arrange(population_sample_type,desc(parameter_value)), - 'Serology (%)','population_sample_type',c(0,100),text_size = 13) - -sero_forest_pop_group <- forest_plot(sero_studies |> mutate(parameter_value = coalesce(parameter_value,central)) |> - arrange(population_group,desc(parameter_value)), - 'Serology (%)','population_group',c(0,100),text_size = 13) - -sero_forest_sero_CSF <- forest_plot(sero_studies |> mutate(parameter_value = coalesce(parameter_value,central)) |> - arrange(sero_CSF,desc(parameter_value)), - 'Serology (%)','sero_CSF',c(0,100),text_size = 13) - -sero_forest_country <- forest_plot(sero_studies |> mutate(parameter_value = coalesce(parameter_value,central)) |> - arrange(population_country,desc(parameter_value)), - 'Serology (%)','population_country',c(0,100),text_size = 13) - -#why is seroprevalence so low in people under investigation?? - -d0 <- parameters %>% filter(parameter_class == 'Seroprevalence') %>% - filter( (population_sample_type == 'Population-Based' | population_sample_type == 'Community-Based') & - population_group == 'General population') - -forest_plot(d0 |> mutate(parameter_value = coalesce(parameter_value,central)) |> - arrange(population_country,desc(parameter_value)), - 'Serology (%)','population_country',c(0,25),text_size = 13) - - -d1 <- parameters %>% filter(parameter_class == 'Seroprevalence' & - parameter_type != 'IgG' & - parameter_type != 'IgM' & - (population_sample_type == 'Population-Based' | - population_sample_type == 'Community-Based') & - population_group == 'General population') -d2 <- parameters %>% filter(parameter_type == 'IgG' & - (population_sample_type == 'Population-Based' | - population_sample_type == 'Community-Based') & - population_group == 'General population') +# *----------------------------- Data preparation -----------------------------* +parameters <- parameters |> + mutate(article_label = make.unique(refs)) |> + mutate(article_label = factor(article_label,levels=rev(unique(article_label)))) |> + mutate(in_CSF = case_when(str_detect(parameter_notes,'CSF')~TRUE, + TRUE~FALSE)) + +sero_studies <- parameters |> + filter(parameter_class == 'Seroprevalence') + +sero_no_unit_row_count <- sero_studies |> + filter(parameter_unit=="Unspecified" | is.na(parameter_unit)) |> + NROW() + +cat("Number of extracted serology rows with no unit: ", sero_no_unit_row_count) + +# Assumption: parameter unit is always a percentage +sero_studies <- sero_studies |> + mutate(parameter_value = coalesce(parameter_value,central), + parameter_unit="Percentage (%)", + population_group = factor( + population_group, + levels = c(sort(setdiff(unique(population_group), + c("Other", "Unspecified"))), + "Other", "Unspecified"))) + +# CSF studies +sero_studies |> + filter(in_CSF) |> + select(article_id, covidence_id, parameter_type, parameter_value, + qa_score, refs) |> + print() +# *---------------------------------- Plots -----------------------------------* +p1 <- forest_plot(sero_studies, 'Serology (%)', 'parameter_type', c(-4,104), + qa_alpha = 0.3, text_size = 28, sort=TRUE, point_size=6) + + ggforce::facet_col(facets = vars(population_country), + scales = "free_y", + space = "free") + + guides(fill = guide_none(), + linetype = guide_none(), + color = guide_legend(title = "Assay", order =2), + shape = guide_legend(title = "Parameter type", order =1)) + + theme(legend.position = c(0.84, 0.94)) +ggsave(paste0("sero_apx_col_assay_pop_country.pdf"), + plot = p1, width = 19, height = 35) + +p2 <- forest_plot(sero_studies, 'Serology (%)', 'parameter_type', c(-4,104), + qa_alpha = 0.3, text_size = 28, sort=TRUE, point_size=6) + + ggforce::facet_col(facets = vars(population_group), + scales = "free_y", + space = "free") + + guides(fill = guide_none(), + linetype = guide_none(), + color = guide_legend(title = "Assay", order =2), + shape = guide_legend(title = "Parameter type", order =1)) + + theme(legend.position = c(0.84, 0.4)) +ggsave(paste0("sero_apx_col_assay_pop_group.pdf"), + plot = p2, width = 19, height = 35) + +p3 <- forest_plot(sero_studies, 'Serology (%)', 'population_group', c(-4,104), + qa_alpha = 0.3, text_size = 28, sort=TRUE, point_size=6) + + ggforce::facet_col(facets = vars(population_country), + scales = "free_y", + space = "free") + + guides(fill = guide_none(), + linetype = guide_none(), + color = guide_legend(title = "Population group", order =2), + shape = guide_legend(title = "Parameter type", order =1)) + + theme(legend.position = c(0.86, 0.93)) +ggsave(paste0("sero_apx_col_country_pop_country.pdf"), + plot = p3, width = 19, height = 35) + +# SI figures +# Split figures to allow for larger panels and figures to split over two pages +# Figure 1 by assay type +sero_apx_1 <- (p1+ theme(legend.position = c(0.76, 0.94)) | + p2 + guides(color=guide_none(), shape=guide_none())) + + plot_annotation(tag_levels = "A") + +ggsave(paste0("sero_apx_col_assay_patchwork.pdf"), + plot = sero_apx_1, width = 32, height = 38) diff --git a/src/nipah_severity/nipah_severity.R b/src/nipah_severity/nipah_severity.R index 2cf323b7..026a22df 100644 --- a/src/nipah_severity/nipah_severity.R +++ b/src/nipah_severity/nipah_severity.R @@ -9,6 +9,7 @@ library(orderly2) library(patchwork) library(readr) library(stringr) +library(tidyr) # *--------------------------------- Orderly ----------------------------------* orderly_parameters(pathogen = NULL) @@ -17,11 +18,14 @@ 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_Bangladesh_IEDCR.csv"="NIPAH_Bangladesh_IEDCR.csv") +orderly_shared_resource("cleaned_outbreak_data.RDS"="cleaned_outbreak_data.RDS") + source("nipah_functions.R") orderly_artefact("Nipah severity figures", - c("figure_severity.png", "figure_severity.pdf")) - + c(file.path("figures", "figure_severity.png"), + file.path("figures","figure_severity.pdf"))) # *------------------------------ Data curation -------------------------------* articles <- read_csv("articles.csv") outbreaks <- read_csv("outbreaks.csv") @@ -37,13 +41,64 @@ qa_scores <- articles |> dplyr::select(covidence_id,qa_score) parameters <- dfs$parameters |> left_join(qa_scores) +subcolumns_outbreak <- readRDS("cleaned_outbreak_data.RDS") # *----------------------------- Data preparation -----------------------------* +# From extracted outbreaks +iso_lookup <- c( + "India" = "IND", + "Bangladesh" = "BGD", + "Malaysia" = "MYS", + "Philippines" = "PHL", + "Singapore" = "SGP" +) + +cfr_from_outbreaks <- subcolumns_outbreak |> + mutate(outbreak_source = replace_na(outbreak_source, 'Unknown')) |> + mutate(cfr_ifr_denominator = round(total_cases), + cfr_ifr_numerator = round(deaths), + CFR = cfr_ifr_numerator / cfr_ifr_denominator, + article_refs = refs, + refs = paste(outbreak_country, outbreak_location, sep = " |> "), + parameter_value = CFR, + parameter_unit = 'Percentage') |> arrange(desc(CFR)) |> + mutate(outbreak_source = case_when(str_detect(outbreak_source, 'Domestic animal' ) ~ 'Domestic animal', + str_detect(outbreak_source, 'Wild animal' ) ~ 'Wild animal', + TRUE ~ outbreak_source), + refs_full = refs, + refs = paste0( + iso_lookup[trimws(str_extract(refs, "^[^|]+"))], " |> ", + # first 3 letters of each word in location + str_to_upper(str_replace_all(str_extract(refs, "(?<=\\|>).*$"), + "\\b([A-Za-z]{1,3})[A-Za-z]*\\b", + "\\1"))) + ) + +# IEDCR +cfr_bangladesh <- read_csv('NIPAH_Bangladesh_IEDCR.csv') + +cfr_from_bangladesh_surveillance <- cfr_bangladesh |> + mutate(year_cat = case_when( + Year %in% 2000:2004 ~ "2000-2004", + Year %in% 2005:2009 ~ "2005-2009", + Year %in% 2010:2014 ~ "2010-2014", + Year %in% 2015:2019 ~ "2015-2019", + Year %in% 2020:2029 ~ "2020-Present", + TRUE ~ "Unspecified")) |> + filter(Cases!=0) |> + mutate(cfr_ifr_denominator = Cases, + cfr_ifr_numerator = Death, + CFR = cfr_ifr_numerator / cfr_ifr_denominator, + refs = Year, + parameter_value = CFR, + parameter_unit = 'Percentage') |> arrange(desc(CFR)) + +# Extracted CFRs parameters d1 <- parameters |> filter(parameter_type == 'Severity - case fatality rate (CFR)') # Unspecified parameter type for CFR with only numerator and denom -d1 <- d1 |> - filter(parameter_unit != "Unspecified") +# d1 <- d1 |> +# filter(parameter_unit != "Unspecified") # Create subgroups here # Lassa examples for study midyear, population_group, cfr_denom_cat @@ -54,8 +109,11 @@ d1 <- d1 |> population_group == "Healthcare workers" ~ "Healthcare workers", population_group == "Abattoir workers" ~ "Animal workers", population_group == "Animal workers" ~ "Animal workers", + is.na(population_group) ~ "Other", TRUE ~ "Mixed Groups")) |> - mutate(population_study_start_year = as.numeric(population_study_start_year), + mutate(parameter_unit = 'Percentage (%)', + parameter_value = coalesce(parameter_value, central), #using central where no % was reported + population_study_start_year = as.numeric(population_study_start_year), population_study_end_year = as.numeric(population_study_end_year), study_midyear = ifelse(!is.na(population_study_start_year) & !is.na(population_study_end_year), round((population_study_start_year + population_study_end_year) / 2), @@ -69,44 +127,443 @@ d1 <- d1 |> mutate(cfr_denom_cat = case_when( cfr_ifr_denominator %in% 1:29 ~ "Reported Cases < 30", cfr_ifr_denominator %in% 30:99 ~ "Reported Cases = 30-99", - cfr_ifr_denominator %in% 100:399 ~ "Reported Cases = 100-399", + cfr_ifr_denominator %in% 100:329 ~ "Reported Cases = 100-329", TRUE ~ "Unspecified")) |> # Deduplicating CFRs mutate(duplicate_cfr = case_when( - covidence_id %in% c() ~ "Known", - covidence_id %in% c() ~ "Assumed", - TRUE ~ "False")) #only identified for estimates passed to meta-analysis (i.e. denominator not NA) + access_param_id %in% c("037_006", "039_001", "081_002", + "172_004", "172_002") ~ "Known", + access_param_id %in% c( + "002_001", "025_006", "030_003", "037_005", "038_001", "040_008", + "044_002", "044_007", "052_002", "093_001", "052_001", "054_001", + "103_001", "109_001", "109_002", "109_003", "109_006", "109_004", + "109_005", "109_007", "109_008", "109_011", "109_012", "109_013", + "109_014", "109_015", "109_016", "109_017", "109_018", "109_019", + "109_020", "109_021", "113_001", "129_002", "121_001", "129_003", + "138_004", "151_004", "171_004", "173_001", "172_001", "179_001", + "190_001") ~ "Assumed", + TRUE ~ "False")) |> + mutate(population_group = factor(population_group, + levels = c(sort(setdiff(unique(population_group), + c("Other", "Unspecified"))), + "Other", "Unspecified"))) + +# proportion of symptomatic cases +d2 <- parameters |> + filter(parameter_type == "Severity - proportion of symptomatic cases") # *------------------------------ Meta-analysis -------------------------------* -da <- d1 |> filter(duplicate_cfr == "False") - -m1 <- metaprop_wrap(dataframe = da, subgroup = "population_country", - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = "red", - width = 4000, height = 4000, resolution = 500) -m2 <- metaprop_wrap(dataframe = d1, subgroup = "study_midyear_cat", - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = "red", - width = 4000, height = 4000, resolution = 500) -m3 <- metaprop_wrap(dataframe = da |> arrange(cfr_ifr_denominator), subgroup = "cfr_denom_cat", - plot_pooled = TRUE, sort_by_subg = FALSE, plot_study = FALSE, digits = 3, colour = "red", - width = 4000, height = 4000, resolution = 500) -m4 <- metaprop_wrap(dataframe = da, subgroup = "population_group", - plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = FALSE, digits = 3, colour = "red", - width = 4000, height = 4000, resolution = 500) - -p1 <- m1$plot -p2 <- m2$plot -p3 <- m3$plot -p4 <- m4$plot - -patchwork <- (p1 + p2 + p3 + p4) + plot_layout(ncol = 2, widths = c(1,1)) -patchwork <- patchwork + plot_annotation(tag_levels = 'A') -ggsave("figure_severity.png", plot = patchwork, width = 18, height = 12) -ggsave("figure_severity.pdf", plot = patchwork, width = 18, height = 12) +# Plot file structure - many plots created in this task so better to create a +# folder structure +filepath_vec <- c(file.path("figures"), + file.path("figures", "extracted_parameters", "all"), + file.path("figures","extracted_parameters", "qa_filtered")) + +for (filepath in filepath_vec){ + if (!dir.exists(filepath)) { + dir.create(filepath, recursive = TRUE) + } +} + + +# Plot colour +imperial_khaki <- "#EFE58B" +text_size <- 15 +lanonc_colours <- ggsci::pal_lancet("lanonc")(9) +meta_digits <- 2 + +# *----------------------- CFR from extracted outbreaks -----------------------* +# Overall +cfr_outbreak_country <- metaprop_wrap( + cfr_from_outbreaks, subgroup = 'outbreak_country', plot_pooled = TRUE, + sort_by_subg = TRUE, plot_study = FALSE, digits = meta_digits, + colour = imperial_khaki, + width = 9500, height = 6000, resolution = 1000) + +ggsave(file.path("figures", "figure_3_meta_country_extracted_outbreak.pdf"), + cfr_outbreak_country$plot, width = 12, height = 6) +ggsave(file.path("figures", "figure_3_meta_country_extracted_outbreak.png"), + cfr_outbreak_country$plot, width = 12, height = 6) + +# With study breakdown +# Update refs for nicer study printing +# max_width <- max(nchar(paste(cfr_from_outbreaks$outbreak_location, +# cfr_from_outbreaks$outbreak_start_year))) +# cfr_from_outbreaks <- cfr_from_outbreaks |> +# mutate(refs = paste0(str_pad(paste(outbreak_location, outbreak_start_year), +# width = 40, side = "right"), " | ", +# article_refs)) + +cfr_outbreak_country_study <- metaprop_wrap( + cfr_from_outbreaks, + subgroup = 'outbreak_country', plot_pooled = TRUE, + sort_by_subg = TRUE, plot_study = TRUE, digits = meta_digits, + colour = imperial_khaki, + width = 11000, height = 17000, resolution = 1000) + +# Update plot to allow for new format +cfr_outbreak_country_study$result$article_refs <- cfr_from_outbreaks$article_refs +cfr_outbreak_country_study$result$studlab <- paste0( + cfr_from_outbreaks$outbreak_start_year, ", ", cfr_from_outbreaks$outbreak_location) + +colour <- imperial_khaki +png(file = "temp.png", width = 13000, height = 17000, res = 1000) +par(mar = c(2, 2, 2, 1)) +cfr_outbreak_country_study_plot <- forest(cfr_outbreak_country_study$result, layout = "Revman5", + leftcols = c("studlab", "article_refs", "event", "n", "effect.ci"), + leftlabs = c("Outbreak", "Study", "Events", "Total"), + just.addcols = "left", + colgap.forest.left = "3mm", + overall = TRUE, pooled.events = TRUE, + print.subgroup.name = FALSE, sort.subgroup = TRUE, + study.results = TRUE, + digits = 2, + col.diamond.lines = "black",col.diamond.common = colour, + col.diamond.random = colour, + col.square = colour, col.square.lines = "black", + col.study = "black", col.subgroup = "black", + col.inside = "black", weight.study = "same", + at = seq(0,1,by=0.2), xlim = c(0,1), xlab="Case Fatality Ratio", + fs.predict.labels = 11.5, + fs.hetstat=11, + fs.test.subgroup = 11, + fs.axis = 11, + fontsize = 14, + plotwidth = "72.5mm") +dev.off() + +pg <- png::readPNG("temp.png", native = TRUE) +file.remove("temp.png") +cfr_outbreak_country_study_plot <- wrap_elements( + plot = rasterGrob(pg, interpolate = TRUE)) + +ggsave(file.path("figures", "SI_CFR_meta_country_extracted_outbreak.pdf"), + cfr_outbreak_country_study_plot, width = 5, height = 6.5) +ggsave(file.path("figures", "SI_CFR_meta_country_extracted_outbreak.png"), + cfr_outbreak_country_study_plot,width = 5, height = 6.5) + +# *--------------------- CFR from Bangladesh surveillance ---------------------* +# Overall +cfr_from_bangladesh_surveillance_yc <- metaprop_wrap( + cfr_from_bangladesh_surveillance, subgroup = "year_cat", plot_pooled = TRUE, + sort_by_subg = TRUE, plot_study = FALSE, digits = meta_digits, + colour = imperial_khaki, + width = 9500, height = 7000, resolution = 1000) + +ggsave(file.path("figures", "figure_3_meta_year_cat_IEDCR.pdf"), + cfr_from_bangladesh_surveillance_yc$plot, width = 10, height = 8) +ggsave(file.path("figures", "figure_3_meta_year_cat_IEDCR.png"), + cfr_from_bangladesh_surveillance_yc$plot, width = 10, height = 8) + +# With study breakdown +cfr_from_bangladesh_surveillance_yc_study <- metaprop_wrap( + cfr_from_bangladesh_surveillance, subgroup = 'year_cat', + plot_pooled = TRUE, sort_by_subg = TRUE, plot_study = TRUE, + digits = meta_digits, colour = imperial_khaki, + width = 10400, height = 11000, resolution = 1000) + +ggsave(file.path("figures", "SI_CFR_meta_year_cat_IEDCR.pdf"), + cfr_from_bangladesh_surveillance_yc_study$plot, + width = 7.5, height = 8.25) +ggsave(file.path("figures", "SI_CFR_meta_year_cat_IEDCR.png"), + cfr_from_bangladesh_surveillance_yc_study$plot, + width = 7.5, height = 8.25) + +# *---------------------- CFR from extracted parameters -----------------------* +# Extracted CFRs +qa_thresh_vec <- c("all"=-1, "qa"=0.5) +qa_alpha_vec <- c(0.3, 1) +list_label_vec <- c("all", "qa_filtered") +plot_list <- list("all"=list("meta"=list(), "forest"=list()), + "qa_filtered"=list("meta"=list(), "forest"=list())) +labels <- c("SI_allqa", "") +cfr_duplicates <- list("no_dups"="False", + "no_known_dups"=c("False", "Assumed"), + "all"=c("False", "Assumed", "Known")) + +for (i in seq_along(qa_thresh_vec)){ + list_label <- list_label_vec[i] + qa_threshold <- qa_thresh_vec[i] + qa_alpha <- qa_alpha_vec[i] + plot_type <- file.path("figures", "extracted_parameters", list_label) + + d1_filtered <- d1 |> filter(qa_score>qa_threshold, + duplicate_cfr=="False") + d2_filtered <- d2 |> filter(qa_score>qa_threshold) + + # Colours: + all_pop_groups <- d1_filtered |> + distinct(population_group) |> + arrange(population_group == "Other", population_group) |> + pull() + + custom_colour_pop_groups <- lanonc_colours[seq_along(all_pop_groups)] + names(custom_colour_pop_groups) <- all_pop_groups + + all_countries <- d1_filtered |> + distinct(population_country) |> + arrange(population_country) |> + pull() + + custom_colour_countries <- lanonc_colours[seq_along(all_countries)] + names(custom_colour_countries) <- all_countries + + plot_list[[list_label]][["meta"]][["m1"]] <- metaprop_wrap( + dataframe = d1_filtered, subgroup = "population_country", plot_pooled = TRUE, + sort_by_subg = TRUE, plot_study = TRUE, digits = meta_digits, + colour = imperial_khaki, width = 10500, height = 9000, resolution = 1000) + + ggsave(file.path(plot_type, + paste0("figure_3_meta_population_country.pdf")), + plot_list[[list_label]][["meta"]][["m1"]]$plot, + width = 5.8, height = 5) + ggsave(file.path(plot_type, + paste0("figure_3_meta_population_country.png")), + plot_list[[list_label]][["meta"]][["m1"]]$plot, + width = 5.8, height = 5) + +plot_list[[list_label]][["meta"]][["m2"]] <- metaprop_wrap( + dataframe = d1_filtered, subgroup = "study_midyear_cat", + plot_pooled = TRUE, sort_by_subg = TRUE, + plot_study = TRUE, digits = meta_digits, + colour = imperial_khaki, + width = 10500, height = 9000, resolution = 1000) + + ggsave(file.path(plot_type, + paste0("figure_3_meta_study_midyear_cat.pdf")), + plot_list[[list_label]][["meta"]][["m2"]]$plot, + width = 8.5, height = 7) + ggsave(file.path(plot_type, + paste0("figure_3_meta_study_midyear_cat.png")), + plot_list[[list_label]][["meta"]][["m2"]]$plot, + width = 8.5, height = 7) + + # Update to account for dedup: + d1_filtered <- d1_filtered |> + mutate(cfr_denom_cat = ifelse(cfr_denom_cat=="Reported Cases = 100-329", + "Reported Cases = 100-250", cfr_denom_cat)) + + plot_list[[list_label]][["meta"]][["m3"]] <- metaprop_wrap( + dataframe = d1_filtered, subgroup = "cfr_denom_cat", plot_pooled = TRUE, + sort_by_subg = FALSE, plot_study = TRUE, digits = meta_digits, + colour = imperial_khaki, width = 10500, height = 8000, resolution = 1000) + + ggsave(file.path(plot_type, + paste0("figure_3_meta_cfr_denom_cat.pdf")), + plot_list[[list_label]][["meta"]][["m3"]]$plot, + width = 8, height = 6) + ggsave(file.path(plot_type, + paste0("figure_3_meta_cfr_denom_cat.png")), + plot_list[[list_label]][["meta"]][["m3"]]$plot, + width = 8, height = 6) + + plot_list[[list_label]][["meta"]][["m4"]] <- metaprop_wrap( + dataframe = d1_filtered, subgroup = "population_group", plot_pooled = TRUE, + sort_by_subg = TRUE, plot_study = TRUE, digits = meta_digits, + colour = imperial_khaki, width = 10000, height = 9000, resolution = 1000) + + ggsave(file.path(plot_type, + paste0("figure_3_meta_population_group.pdf")), + plot_list[[list_label]][["meta"]][["m4"]]$plot, + width = 7.3, height = 6) + ggsave(file.path(plot_type, + paste0("figure_3_meta_population_group.png")), + plot_list[[list_label]][["meta"]][["m4"]]$plot, + width = 7.3, height = 6) + + # Forest plot + plot_list[[list_label]][["forest"]][["p_cfr_1"]] <- forest_plot( + d1_filtered, "Case-Fatality Ratio (%)", "population_country", + c(-10,110), custom_colours = custom_colour_countries, + text_size=text_size, qa_alpha=qa_alpha, sort=TRUE) + + guides(shape = guide_legend(title = "Parameter type", order=1), + fill = guide_none(), + linetype = guide_none(), + color = guide_legend(title = "Population country", order=3)) + ggsave(file.path(plot_type, + paste0("figure_3_forest_cfr_population_country.pdf")), + plot = plot_list[[list_label]][["forest"]][["p_cfr_1"]], + width = 10, height = 12) + + plot_list[[list_label]][["forest"]][["p_cfr_2"]] <- forest_plot( + d1_filtered, "Case-Fatality Ratio (%)","population_group", + c(-10,110), custom_colours = custom_colour_pop_groups, + text_size=text_size, qa_alpha=qa_alpha, sort=TRUE) + + guides(shape = guide_legend(title = "Parameter type", order=1), + fill = guide_none(), + linetype = guide_none(), + color = guide_legend(title = "Population group", order=2)) + + ggsave(file.path(plot_type, + paste0("figure_3_forest_cfr_population_group.pdf")), + plot = plot_list[[list_label]][["forest"]][["p_cfr_2"]], + width = 10, height = 12) + + plot_list[[list_label]][["forest"]][["p_prop_1"]] <- forest_plot( + d2_filtered, "Percentage of Symptomatic Cases (%)", "population_country", + c(-10, 110), custom_colours = custom_colour_countries, + text_size=text_size, qa_alpha=qa_alpha, sort=TRUE) + + guides(color = guide_legend(title = "Population country", order=2), + linetype = guide_none(), + shape = guide_legend(title = "Parameter type", order=1)) + + ggsave(file.path(plot_type, + paste0("figure_3_forest_prop_country.pdf")), + plot = plot_list[[list_label]][["forest"]][["p_prop_1"]], + width = 8, height = 5) + + plot_list[[list_label]][["forest"]][["p_prop_2"]] <- forest_plot( + d2_filtered, "Percentage of Symptomatic Cases (%)", "population_group", + c(-10, 110), custom_colours = custom_colour_pop_groups, + text_size=text_size, qa_alpha=qa_alpha, sort=TRUE) + + guides(color = guide_legend(title = "Population group", order=2), + linetype = guide_none(), + shape = guide_legend(title = "Parameter type", order=1)) + + ggsave(file.path(plot_type, + paste0("figure_3_forest_prop_population_group.pdf")), + plot = plot_list[[list_label]][["forest"]][["p_prop_2"]], + width = 8, height = 5) +} + +# Additional CFR from extracted parameters results +# Can't do year since the estimates may relate to a range and outbreak location +# is not clean (possible that there's a single estimate for multiple locations) +# No duplicates +d1_dup_false <- d1 |> filter(duplicate_cfr=="False") + +meta_dup_false <- metaprop_wrap( + dataframe = d1_dup_false, subgroup = "population_country", plot_pooled = TRUE, + sort_by_subg = TRUE, plot_study = TRUE, digits = meta_digits, + colour = imperial_khaki, width = 10500, height = 9000, resolution = 1000) + +ggsave(file.path("figures", "extracted_parameters", + "SI_CFR_meta_country_param_dup_eq_false.pdf"), + meta_dup_false$plot, + width = 5.8, height = 5) + +# Assumed duplicates included +# Error when trying to fit this? +d1_dup_assumed <- d1 |> + filter(duplicate_cfr!="Known") + +meta_dup_assumed <- metaprop_wrap( + dataframe = d1_dup_assumed, subgroup = "population_country", plot_pooled = TRUE, + sort_by_subg = TRUE, plot_study = TRUE, digits = meta_digits, + colour = imperial_khaki, width = 10500, height = 18000, resolution = 1000) + +ggsave(file.path("figures", "extracted_parameters", + "SI_CFR_meta_country_param_dup_neq_known.pdf"), + meta_dup_assumed$plot, + width = 5.8, height = 10) + +# All +meta_all <- metaprop_wrap( + dataframe = d1, subgroup = "population_country", plot_pooled = TRUE, + sort_by_subg = TRUE, plot_study = TRUE, digits = meta_digits, + colour = imperial_khaki, width = 10500, height = 20000, resolution = 1000) + +ggsave(file.path("figures", "extracted_parameters", + "SI_CFR_meta_country_param_dup_all.pdf"), + meta_all$plot, + width = 5.8, height = 10.8) + +# Combine figures +p1 <- plot_list[["all"]][["meta"]][["m1"]]$plot + + theme(plot.margin = margin(-50, -250, -250, -50)) +p2 <- cfr_outbreak_country$plot + + theme(plot.margin = margin(-250, -250, -250, -250)) +p3 <- cfr_from_bangladesh_surveillance_yc$plot + + theme(plot.margin = margin(-50, -250, -50, -250)) +p4 <- plot_list[["all"]][["forest"]][["p_cfr_1"]] + + theme(legend.position=c(0.275, 0.775), legend.direction = "vertical") +p5 <- plot_list[["all"]][["forest"]][["p_prop_1"]] + + guides(color = guide_none(), + linetype = guide_none(), + shape = guide_none()) + +left_col <- p4 / p5 + plot_layout(heights = c(26, 3)) +right_col <- p1 / p2 / p3 + + plot_layout(heights = c(2.4, 2.4, 2.45)) + +patchwork <- (left_col | right_col) + + plot_layout(widths = c(2.5, 4)) + + plot_annotation(tag_levels = "A") & + theme(plot.tag.position = c(0, 1), plot.margin = margin(5.5, 0, 0, 5.5), + plot.tag = element_text(size = 14)) + +ggsave(file.path("figures", "figure_severity_forest.png"), + plot = left_col, width = 8, height = 10) +ggsave(file.path("figures","figure_severity_forest.pdf"), + plot = left_col, width = 8, height = 10) + +ggsave(file.path("figures", "figure_severity.png"), + plot = patchwork, width = 11.5, height = 10, dpi=300) +ggsave(file.path("figures", "figure_severity.pdf"), + plot = patchwork, width = 11.5, height = 10, dpi=300) # Additional plots for SI -#figure_S6-S10: meta-analysis with all estimates plotted -#figure_S11: meta-analysis with only known duplicates excluded -db <- d1 |> filter(duplicate_cfr %in% c("False","Assumed")) -#figure_S12: meta-analysis without de-duplication -dc <- d1 +# Since loop only uses deduplicated CFR, recreate the plot (inefficient) +all_pop_groups <- d1 |> + distinct(population_group) |> + arrange(population_group == "Other", population_group) |> + pull() + +custom_colour_pop_groups <- lanonc_colours[seq_along(all_pop_groups)] +names(custom_colour_pop_groups) <- all_pop_groups + +all_countries <- d1 |> + distinct(population_country) |> + arrange(population_country) |> + pull() + +custom_colour_countries <- lanonc_colours[seq_along(all_countries)] +names(custom_colour_countries) <- all_countries + +p1 <- forest_plot( + d1, "Case-Fatality Ratio (%)","population_country", + c(-10,110), custom_colours = custom_colour_countries, + text_size=text_size, qa_alpha=0.3, sort=TRUE) + + guides(shape = guide_legend(title = "Parameter type", order=1), + fill = guide_none(), + linetype = guide_none(), + color = guide_legend(title = "Country", order=2)) + +p2 <- forest_plot( + d1, "Case-Fatality Ratio (%)","population_group", + c(-10,110), custom_colours = custom_colour_pop_groups, + text_size=text_size, qa_alpha=0.3, sort=TRUE) + + guides(shape = guide_legend(title = "Parameter type", order=1), + fill = guide_none(), + linetype = guide_none(), + color = guide_legend(title = "Population group", order=2)) + +patchwork_si <- p1 + + theme(legend.position=c(0.175, 0.7), legend.direction = "vertical") + + p2 + + theme(legend.position=c(0.25, 0.7), legend.direction = "vertical") + + guides(shape=guide_none()) + + plot_layout(ncol = 2) + plot_annotation(tag_levels = 'A') + +ggsave(file.path("figures", "figure_SI_severity.png"), + plot = patchwork_si, width = 15, height = 18) +ggsave(file.path("figures", "figure_SI_severity.pdf"), + plot = patchwork_si, width = 15, height = 18) + +p3 <- forest_plot( + d1, "Case-Fatality Ratio (%)","population_group", + c(-10,110), custom_colours = custom_colour_pop_groups, + text_size=text_size,qa_alpha=0.3, sort=TRUE) + + guides(shape = guide_legend(title = "Parameter type", order=1), + fill = guide_none(), + linetype = guide_none(), + color = guide_legend(title = "Population group", order=2)) + + ggforce::facet_col(facets = vars(population_country), + scales = "free_y", + space = "free") + + theme(legend.position=c(0.2, 0.875)) +ggsave(file.path("figures", "figure_SI_severity_facet.png"), + plot = p3, width = 10, height = 18) +ggsave(file.path("figures", "figure_SI_severity_facet.pdf"), + plot = p3, width = 10, height = 18) diff --git a/src/nipah_transmission/nipah_transmission.R b/src/nipah_transmission/nipah_transmission.R index 3b098b0d..dd1d0b57 100644 --- a/src/nipah_transmission/nipah_transmission.R +++ b/src/nipah_transmission/nipah_transmission.R @@ -35,12 +35,22 @@ 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") d4 <- parameters |> filter(parameter_class == "Attack rate") -d5 <- parameters |> filter(parameter_type == "Severity - proportion of symptomatic cases") -d6 <- parameters |> filter(parameter_class == "Reproduction number") +d5 <- parameters |> filter(parameter_class == "Reproduction number") # Add prop symptomatic # arrange data and format for plotting @@ -52,20 +62,26 @@ variables_to_mutate <- c("parameter_value", d1 <- d1 |> mutate(across(all_of(variables_to_mutate), ~ . * 10^4)) #multiply by 10^4 d2 <- d2 |> mutate(across(all_of(variables_to_mutate), ~ . * 10^4)) #multiply by 10^4 -d4 <- d4 |> mutate(across(all_of(variables_to_mutate), - ~ ifelse(parameter_unit == "No units", . * 100, .)) - ) |> +d4 <- d4 |> + mutate(across(all_of(variables_to_mutate), + ~ ifelse(parameter_unit == "No units", . * 100, .))) |> 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 d2 <- d2 |> mutate(across( c(parameter_upper_bound, parameter_lower_bound), - ~ ifelse(covidence_id == 2760, NA, .)) - ) + ~ ifelse(covidence_id == 2760, NA, .))) d2 <- d2 |> arrange(genome_site,-central) @@ -77,21 +93,12 @@ d3 <- d3 |> arrange(-central) d4 <- d4 |> mutate(arate=c("Primary","Primary")) |> - arrange(arate,-central) + arrange(arate,-central) -# Update central d5 <- d5 |> - mutate(parameter_value = ifelse(is.na(parameter_value), - cfr_ifr_numerator/cfr_ifr_denominator, - parameter_value) - ) -d5 <- d5 |> - arrange(-central) - -d6 <- d6 |> arrange(parameter_type, -central) -d6 <- d6 |> +d5 <- d5 |> mutate(parameter_type = factor(parameter_type, levels = unique(parameter_type), labels = c("Basic (R0)"))) @@ -103,7 +110,7 @@ text_size <- 12 # Get custom colours so that genome has different colours lanonc_colours <- ggsci::pal_lancet("lanonc")(9) -all_pop_groups <- bind_rows(d3, d4, d5, d6) |> +all_pop_groups <- bind_rows(d3, d4, d5, d5) |> distinct(population_group) |> # arrange alphabetically but put other last arrange(population_group == "Other", population_group) |> @@ -133,78 +140,85 @@ 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(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_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, 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)), + linetype = 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()) + # 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)", + "Overdispersion (max nr. of secondary cases)", "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(), + linetype = guide_none(), shape = guide_none()) p4 <- forest_plot(d4 |> filter(qa_score>qa_threshold), - "Primary Attack Rate (%)", + "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(), + linetype = guide_none(), shape = guide_none()) p5 <- forest_plot(d5 |> filter(qa_score>qa_threshold), - "Proportion of Symptomatic Cases (%)", - "population_group", - c(-5, 110), custom_colours = custom_colour_pop_groups, - text_size=text_size, qa_alpha=qa_alpha) + - guides(color = guide_none(), - shape = guide_none()) - - p6 <- forest_plot(d6 |> filter(qa_score>qa_threshold), "Basic Reproduction Number", "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, limits=names(custom_colour_pop_groups)) + guides(fill = guide_none(), - color = guide_legend(title = "Population type", order = 1, - override.aes = list(fill = custom_colour_pop_groups)), - shape=guide_none()) - + linetype = guide_none(), + color = guide_legend( + title = "Population type", order = 2, + override.aes = list(fill = custom_colour_pop_groups)), + shape=guide_legend(title = "Parameter type", order=1)) + + p5 <- p5 + theme(legend.position = c(0.835, 0.3)) + p1 <- p1 + theme(legend.position = c(0.875, 0.85)) # Save plots - patchwork <- (p6 + p3 + p4 + p5 + p1 + p2) + - plot_layout(ncol = 2, widths = c(1,1), guides = "collect") + patchwork <- (p5 + p3 + p4 + p1) + + plot_layout(ncol = 2, widths = c(1,1)) + + patchwork <- patchwork + + plot_annotation(tag_levels = "A") + + plot_layout(byrow = FALSE) - patchwork <- patchwork + plot_annotation(tag_levels = "A") ggsave(paste0("figure_",label,"trans.png"), - plot = patchwork, width = 14, height = 10) + plot = patchwork, width = 17, height = 10) ggsave(paste0("figure_",label,"trans.pdf"), - plot = patchwork, width = 14, height = 10) + plot = patchwork, width = 17, height = 10) } # *============================================================================*