In dchodge/pilot.hcw: What the Package Does (One Line, Title Case)

Make classes out of each data set

library(latex2exp)
library(patchwork)
library(here)
require(plyr)
require(data.table)
require(tidyverse)
library(devtools)
library(loo)
# Required for this analysis
require(reshape2)
require(foreach)
require(doParallel)
require(bayesplot)
require(coda)
require(ggplot2)
require(viridis)
require(ggpubr)
library(mvoutlier)

require(foreach)
require(doParallel)
require(Rcpp)

require(extrafont)
require(showtext)
font_add("Avenir", "/Users/davidhodgson/Documents/personal/fonts/avenir_ff/AvenirLTStd-Book.otf")  
showtext_auto()
theme_set(theme_light())
library(wesanderson)

Load packages and data

Currently using the serosolver

#install_github("https://github.com/dchodge/serosolver", ref = "custom_ab")
#install("/Users/davidhodgson/Documents/research/flu/serosolver_vac")
#install(here::here("rcppfunchcw"))

devtools::load_all()
require(serosolver)
require(rcppfunchcw)

load(here::here("data", "hcwpre_data.RDS")) # hcwpre

get_model_info_hcw_pre_cross(hcwpre)
get_model_info_hcw_pre_full(hcwpre)


````

# install new custom_ab functions
```r

get_new_custom_ab <- function() {
install(here::here("rcppfunchcw"))
require(rcppfunchcw)
devtools::load_all()
get_model_info_hcw_pre_cross(hcwpre)
get_model_info_hcw_pre_full(hcwpre)

load(file = here::here("data", paste0("modelinfo_cross_", hcwpre$study_name_short, ".RDS"))) #all_models_hcw_pre_cross
load(file = here::here("data", paste0("modelinfo_full_", hcwpre$study_name_short, ".RDS"))) #all_models_hcw_pre_full
}

Fitting is done in a separate R script

load(file = here::here("data", paste0("modelinfo_cross_", hcwpre$study_name_short, ".RDS"))) #all_models_hcw_pre_cross
1:9 %>% map(~setup_run_serosolver(
    all_models_hcw_pre_cross[[.x]],
    chains = 4,
    iterations = 1000000,
    pt = TRUE,
    filename = "dis_dep",
    cross_sectional = TRUE,
    continue_run = FALSE)
)


load(file = here::here("data", paste0("modelinfo_full_", hcwpre$study_name_short, ".RDS"))) #all_models_hcw_pre_full
1 %>% map(~setup_run_serosolver(
    all_models_hcw_pre_full[[.x]],
    chains = 4,
    iterations = 5000,
    pt = TRUE,
    filename = "dis_dep",
    cross_sectional = FALSE,
    continue_run = FALSE)
)

Load the mcmc and plot the traces and schematics of titre loss

load(file = here::here("data", paste0("modelinfo_cross_", hcwpre$study_name_short, ".RDS"))) #all_models_hcw_pre_cross
all_post_test_cross <- all_models_hcw_pre_cross %>%
    map(~load_mcmc_post(.x,
        fit_place = "local",
        folder_name_1 = "dis_dep",
        cross_sectional = TRUE,
        burnin = 0,
        thin = 1))

all_post_test_cross_trim <- vector(mode = "list", length = 9)
samp_no_list <- vector(mode = "list", length = 9)
for (p in 1:9) { 
    L <- nrow(all_post_test_cross[[p]]$theta_list_chain[[1]])
    B <- round(L / 2)
    chain_list <- all_post_test_cross[[p]]$theta_list_chain
    pars_plot <- setdiff(all_models_hcw_pre_cross[[p]]$other$pars_plot, "tau_vac" )
    all_post_test_cross_trim[[p]] <- chain_list %>% map(~(.x[B:L, pars_plot] %>% as.data.frame))
    samp_no_list[[p]] <- chain_list %>% map(~.x[B:L, c("sampno")])
}

all_post_test_cross_trim_inf_hist <- vector(mode = "list", length = 9)
samp_no_list <- vector(mode = "list", length = 9)
for (p in 1:9) { 
    L <- nrow(all_post_test_cross[[p]]$inf_list_chains[[1]])
    B <- round(L / 2)
    inf_hist <- all_post_test_cross[[p]]$inf_list_chains
    all_post_test_cross_trim_inf_hist[[p]] <- inf_hist %>%
        map(~(.x %>% filter(sampno == 1200001) %>% as.data.frame %>%
        rename(yr = i, id = j)))
}

post_cross <- list(post_raw = all_post_test_cross, post_trim = all_post_test_cross_trim, inf_hist_trim = all_post_test_cross_trim_inf_hist)
save(post_cross, file = here::here("outputs", "hcw_pre", "fits", "dis_dep", "cross_sec", "trimmed_post_cross.Rdata"))

Comparison figures and statistics in sum_figs folder

purrr::map2(all_post_test_cross, all_models_hcw_pre_cross, ~plot_traces(.x, .y, file = "dis_dep/cross_sec"))
purrr::map2(all_post_test_cross, all_models_hcw_pre_cross, ~plot_histogram(.x, .y, file = "dis_dep/cross_sec"))
purrr::map2(all_post_test_cross, all_models_hcw_pre_cross, ~plot_titre_post(.x, .y, file = "dis_dep/cross_sec"))

# Plots box plots of the 8 different models
plot_post_compare_cross(all_post_test_cross_trim, all_models_hcw_pre_cross, file = "dis_dep/cross_sec")

# Plots the mpsrf for the models to check convergences
mpsrf_df <- purrr::map2(all_post_test_cross_trim[2:9], all_models_hcw_pre_cross[2:9], ~calc_mpsrf(.x, .y)) %>% bind_rows
save_plot_mpsrf_cross(mpsrf_df, all_models_hcw_pre_cross[[1]], file = "dis_dep/cross_sec")

# Plots the waics for the models and compare
save_plot_waic_cross_hcwpre(all_models_hcw_pre_cross, file = "dis_dep")

Get the priors for the full version

# Best fit model is mt according to waic
post_best <- all_post_test_cross_trim[[3]] %>% bind_rows
post_best %>% apply(2, median) %>% as.data.frame %>% t
post_best %>% apply(2, quantile, c(0.01, 0.99)) %>% as.data.frame

Load all the data from the full simulations

load(file = here::here("data", paste0("modelinfo_full_", hcwpre$study_name_short, ".RDS"))) #all_models_hcw_pre_full
all_post_test_full <- all_models_hcw_pre_full %>%
    map(~load_mcmc_post(.x,
        fit_place = "local",
        folder_name_1 = "dis_dep",
        cross_sectional = FALSE,
        burnin = 0,
        thin = 1))

num <- 4

# Trim the the data, use last 50% and throw out outliers 
all_post_test_full_trim <- vector(mode = "list", length = num)
samp_no_list <- vector(mode = "list", length = num)
for (p in 1:num) { 
    L <- nrow(all_post_test_full[[p]]$theta_list_chain[[1]])
    B <- round(L / 2)
    chain_list <- all_post_test_full[[p]]$theta_list_chain
    pars_plot <- setdiff(all_models_hcw_pre_full[[p]]$other$pars_plot, c("tau_vac"))
    all_post_test_full_trim[[p]] <- chain_list %>% map(~(.x[B:L, pars_plot] %>% as.data.frame))
    samp_no_list[[p]] <- chain_list %>% map(~.x[B:L, c("sampno")])
}

all_post_test_full_trim_inf_hist <- vector(mode = "list", length = num)
samp_no_list <- vector(mode = "list", length = num)
for (p in 1:num) { 
    L <- nrow(all_post_test_full[[p]]$inf_list_chains[[1]])
    B <- round(L / 2)
    inf_hist <- all_post_test_full[[p]]$inf_list_chains
    all_post_test_full_trim_inf_hist[[p]] <- inf_hist %>%
        map(~(.x %>% filter(sampno == 1300001) %>% as.data.frame %>%
        rename(yr = i, id = j)))
}

post_full <- list(post_raw = all_post_test_full, post_trim = all_post_test_full_trim, inf_hist_trim = all_post_test_full_trim_inf_hist)
save(post_full, file = here::here("outputs", "hcw_pre", "fits", "dis_dep", "full", "trimmed_post_full.Rdata"))

Plot all the figures

purrr::map2(all_post_test_full, all_models_hcw_pre_full, ~plot_traces(.x, .y, file = "dis_dep/full"))
purrr::map2(all_post_test_full, all_models_hcw_pre_full, ~plot_titre_post(.x, .y, file = "dis_dep/full"))
purrr::map2(all_post_test_full, all_models_hcw_pre_full, ~plot_histogram(.x, .y, file = "dis_dep/full"))

# Plots box plots of the 8 different models
plot_post_compare_full(all_post_test_full_trim, all_models_hcw_pre_full, file = "dis_dep/full")

# Plots the mpsrf for the models to check convergences
mpsrf_df <- purrr::map2(all_post_test_full_trim, all_models_hcw_pre_full, ~calc_mpsrf(.x, .y)) %>% bind_rows
save_plot_mpsrf_full(mpsrf_df, all_models_hcw_pre_full[[1]], file = "dis_dep/full")

# Plots the waics for the models and compare
save_plot_waic_full(all_models_hcw_pre_full, file = "dis_dep", cross_sec = FALSE)