From e2ce73a3c4add5097ed5e2bfc460884e2f5bbe48 Mon Sep 17 00:00:00 2001 From: Nick Croucher Date: Mon, 2 Aug 2021 06:33:20 +0100 Subject: [PATCH] Update serotype analysis --- .../s-pneumoniae-serotype-invasiveness.Rmd | 145 +++++++++++++----- 1 file changed, 105 insertions(+), 40 deletions(-) diff --git a/examples/s-pneumoniae-serotype-invasiveness.Rmd b/examples/s-pneumoniae-serotype-invasiveness.Rmd index 03e326f..09b4068 100644 --- a/examples/s-pneumoniae-serotype-invasiveness.Rmd +++ b/examples/s-pneumoniae-serotype-invasiveness.Rmd @@ -374,9 +374,20 @@ fit_list <- list(null_poisson_fit, adjusted_type_specific_poisson_fit, adjusted_type_specific_negbin_fit) -rhat_plots <- lapply(fit_list, plot, plotfun = "rhat", binwidth = 0.00005) +plot_rhat <- function(fit) { + rhat_plot <- plot(fit, plotfun = "rhat", binwidth = 0.00005) + + theme(axis.text.x = element_text(angle = 90)) + return(rhat_plot) +} + +rhat_plots <- lapply(fit_list, plot_rhat) -cowplot::plot_grid(plotlist = rhat_plots, ncol = 4, nrow = 2, labels = "AUTO") +combined_rhat_plots <- cowplot::plot_grid(plotlist = rhat_plots, ncol = 4, nrow = 2, labels = "AUTO") + +ggsave(combined_rhat_plots, + file = "infant_combined_rhat_plots.png", + width = 8, + height = 6) ``` @@ -386,13 +397,19 @@ Convergence can also be checked using trace plots of the log likelihood values: plot_lp <- function(fit) { pl <- rstan::traceplot(fit, pars = "lp__") - no_leg_pl <- pl + theme(legend.position = "none") + no_leg_pl <- pl + theme(legend.position = "none", + axis.text.x = element_text(angle = 90)) return(no_leg_pl) } trace_plots <- lapply(fit_list, plot_lp) -cowplot::plot_grid(plotlist = trace_plots, ncol = 4, nrow = 2, labels = "AUTO") +combined_trace_plots <- cowplot::plot_grid(plotlist = trace_plots, ncol = 4, nrow = 2, labels = "AUTO") + +ggsave(combined_trace_plots, + file = "infant_combined_trace_plots.png", + width = 8, + height = 6) ``` @@ -408,13 +425,28 @@ model_outputs <- lapply(fit_list, prediction_plots <- lapply(model_outputs, progressionEstimation::plot_case_carrier_predictions, - legend = FALSE) + legend = FALSE, + n_label = 0) -no_leg_plot <- cowplot::plot_grid(plotlist = prediction_plots, ncol = 1, nrow = 8, labels = "AUTO") +no_leg_plot <- cowplot::plot_grid(plotlist = prediction_plots, + ncol = 1, + nrow = 8, + labels = "AUTO", + vjust = -0.01) + + theme(plot.margin = unit(c(0.5,0,0,0), "cm")) -prediction_plot_legend <- progressionEstimation::plot_case_carrier_predictions(model_outputs[[1]], just_legend = TRUE) +prediction_plot_legend <- progressionEstimation::plot_case_carrier_predictions(model_outputs[[1]], just_legend = TRUE) -cowplot::plot_grid(no_leg_plot, prediction_plot_legend, ncol = 1, rel_heights = c(0.9,0.1)) +pred_v_obs_plot <- + cowplot::plot_grid(no_leg_plot, + prediction_plot_legend, + ncol = 1, + rel_heights = c(1,0.1)) + +ggsave(pred_v_obs_plot, + file = "infant_serotypes_pred_v_obs_plot.png", + height = 20, + width = 10) ``` @@ -558,13 +590,19 @@ full_adjusted_type_specific_poisson_fit <- rhat_plot <- plot(full_adjusted_type_specific_poisson_fit, plotfun = "rhat", binwidth = 0.00005) -trace_plot <- rstan::traceplot(full_adjusted_type_specific_poisson_fit, pars = "lp__") +trace_plot <- rstan::traceplot(full_adjusted_type_specific_poisson_fit, pars = "lp__") + + theme(axis.text.x = element_text(angle = 90)) -cowplot::plot_grid(plotlist = list(rhat_plot,trace_plot), +best_fitting_model_validation_plot <- cowplot::plot_grid(plotlist = list(rhat_plot,trace_plot), nrow = 1, ncol = 2, labels = "AUTO") +ggsave(best_fitting_model_validation_plot, + file = "best_fitting_infant_model_validation_plot.png", + width = 8, + height = 6) + ``` ``` {r Plot invasiveness across all serotypes} @@ -676,12 +714,20 @@ pv_full_adjusted_type_specific_poisson_fit <- rhat_plot <- plot(pv_full_adjusted_type_specific_poisson_fit, plotfun = "rhat", binwidth = 0.00005) -trace_plot <- rstan::traceplot(pv_full_adjusted_type_specific_poisson_fit, pars = "lp__") +trace_plot <- rstan::traceplot(pv_full_adjusted_type_specific_poisson_fit, pars = "lp__") + + theme(axis.text.x = element_text(angle = 90)) -cowplot::plot_grid(plotlist = list(rhat_plot,trace_plot), +best_fitting_pv_model_validation_plot <- + cowplot::plot_grid(plotlist = list(rhat_plot,trace_plot), nrow = 1, ncol = 2, labels = "AUTO") + +ggsave(best_fitting_pv_model_validation_plot, + file = "best_fitting_pv_infant_model_validation_plot.png", + width = 8, + height = 6) + ``` ``` {r Analyse the pv output} @@ -888,7 +934,8 @@ combined_invasiveness_df <- dplyr::left_join(sample_sizes, by = c("Serotype" = "type")) # Plot comparison of invasiveness estimates -ggplot(combined_invasiveness_df, +(fixed_effects_comparison <- + ggplot(combined_invasiveness_df, aes(x = Invasiveness_adjusted_type_specific_poisson, xmin = Invasiveness_lower_adjusted_type_specific_poisson, xmax = Invasiveness_upper_adjusted_type_specific_poisson, @@ -915,39 +962,57 @@ ggplot(combined_invasiveness_df, dplyr::filter((Invasiveness_adjusted_type_specific_poisson < 0.0005 & Invasiveness_FE > 0) | Serotype == "46"), aes(label = Serotype) ) +) ``` ``` {r Plot comparison with mixed effects} # Plot comparison of invasiveness estimates -ggplot(combined_invasiveness_df, - aes(x = Invasiveness_adjusted_type_specific_poisson, - xmin = Invasiveness_lower_adjusted_type_specific_poisson, - xmax = Invasiveness_upper_adjusted_type_specific_poisson, - y = Invasiveness_REML, - ymin = Invasiveness_lower_REML, - ymax = Invasiveness_upper_REML, - colour = classification)) + - geom_errorbar(alpha = 0.25) + - geom_errorbarh(alpha = 0.25) + - geom_point() + - scale_x_continuous(trans = "log10") + - theme_bw() + - scale_color_manual(values = vaccine_colours) + - theme(legend.position = "bottom") + - ylab("Invasiveness estimated from log odds ratio") + - xlab("Invasivess estimate from adjusted type-specific Poisson (cases per carrier per year)") + - facet_wrap(~Sample_count) + - stat_cor(aes(x = Invasiveness_adjusted_type_specific_poisson, - y = Invasiveness_REML), - inherit.aes = FALSE, - method = "kendall", - cor.coef.name = "tau") + - ggrepel::geom_label_repel(data = combined_invasiveness_df %>% - dplyr::filter((Invasiveness_adjusted_type_specific_poisson < 0.0005 & Invasiveness_FE > 0) | Serotype == "46"), - aes(label = Serotype) - ) +(random_effects_comparison <- + ggplot(combined_invasiveness_df, + aes(x = Invasiveness_adjusted_type_specific_poisson, + xmin = Invasiveness_lower_adjusted_type_specific_poisson, + xmax = Invasiveness_upper_adjusted_type_specific_poisson, + y = Invasiveness_REML, + ymin = Invasiveness_lower_REML, + ymax = Invasiveness_upper_REML, + colour = classification)) + + geom_errorbar(alpha = 0.25) + + geom_errorbarh(alpha = 0.25) + + geom_point() + + scale_x_continuous(trans = "log10") + + theme_bw() + + scale_color_manual(values = vaccine_colours) + + theme(legend.position = "bottom") + + ylab("Invasiveness estimated from log odds ratio") + + xlab("Invasivess estimate from adjusted type-specific Poisson (cases per carrier per year)") + + facet_wrap(~Sample_count) + + stat_cor(aes(x = Invasiveness_adjusted_type_specific_poisson, + y = Invasiveness_REML), + inherit.aes = FALSE, + method = "kendall", + cor.coef.name = "tau") + + ggrepel::geom_label_repel(data = combined_invasiveness_df %>% + dplyr::filter((Invasiveness_adjusted_type_specific_poisson < 0.0005 & Invasiveness_FE > 0) | Serotype == "46"), + aes(label = Serotype) + ) +) + +combined_odds_ratio_plot <- + cowplot::plot_grid( + plotlist = list(fixed_effects_comparison + theme(legend.position = "none"), + random_effects_comparison), + ncol = 1, + nrow = 2, + rel_heights = c(1,1.2), + labels = "AUTO" + ) + +ggsave(combined_odds_ratio_plot, + file = "combined_odds_ratio_plot.png", + width = 8, + height = 12) ```