Skip to content

Commit

Permalink
Update serotype analysis
Browse files Browse the repository at this point in the history
  • Loading branch information
@nickjcroucher
nickjcroucher committed Aug 2, 2021
1 parent 92cf94f commit e2ce73a
Showing 1 changed file with 105 additions and 40 deletions.
145 changes: 105 additions & 40 deletions examples/s-pneumoniae-serotype-invasiveness.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
```

Expand All @@ -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)
```

Expand All @@ -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)
```

Expand Down Expand Up @@ -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}
Expand Down Expand Up @@ -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}
Expand Down Expand Up @@ -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,
Expand All @@ -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)
```

Expand Down

0 comments on commit e2ce73a

Please sign in to comment.