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)s
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", "hanam_data.RDS")) # hcwpre
get_model_info_hanam(hanam)

````

# install new custom_ab functions

```r

load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS"))) #all_models_hcw_pre_cross

# Cross sectional data
setup_run_serosolver(
    all_models_hanam[[1]],
    chains = 4,
    iterations = 10000,
    pt = TRUE,
    filename = "dis_dep",
    cross_sectional = TRUE,
    continue_run = FALSE)

# Non crossectional data
model_no <- 2
setup_run_serosolver(
    all_models_hanam[[model_no]],
    chains = 4,
    iterations = 10000,
    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_", hanam$study_name_short, ".RDS"))) #all_models_hcw_pre_cross

loadSavePost <- function(model_no, cross_sec, fileoutput) {
    all_post_test_cross_hanam <- load_mcmc_post(all_models_hanam[[model_no]],
        fit_place = "local",
        folder_name_1 = "dis_dep",
        cross_sectional = cross_sec,
        burnin = 0,
        thin = 1)

    all_post_test_cross_trim_hanam <- vector(mode = "list", length = 1)
    samp_no_list <- vector(mode = "list", length = 1)
    for (p in 1:1) { 
        L <- nrow(all_post_test_cross_hanam$theta_list_chain[[1]])
        B <- round(L / 2)
        chain_list <- all_post_test_cross_hanam$theta_list_chain
        pars_plot <- setdiff(all_models_hanam[[model_no]]$other$pars_plot, "tau_vac" )
        all_post_test_cross_trim_hanam[[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_hanam <- vector(mode = "list", length = 1)
    samp_no_list <- vector(mode = "list", length = 1)
    for (p in 1:1) { 
        L <- nrow(all_post_test_cross_hanam$inf_list_chains[[1]])
        B <- round(L / 2)
        inf_hist <- all_post_test_cross_hanam$inf_list_chains
        all_post_test_cross_trim_inf_hist_hanam[[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_hanam, post_trim = all_post_test_cross_trim_hanam, inf_hist_trim = all_post_test_cross_trim_inf_hist_hanam)
    save(post_cross, file = here::here("outputs", "hanam", "fits", "dis_dep", fileoutput))
}

loadSavePost(1, TRUE, "cross_sec/post_cross.RData")
loadSavePost(2, FALSE, "full/LN_post_full.RData")
loadSavePost(3, FALSE, "full/SP_post_full.RData")

Comparison figures and statistics in sum_figs folder

load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS"))) #all_models_hcw_pre_cross

loadPlotFigs <- function(model_no, folder_name, fileoutput) {
    load(file = here::here("outputs", "hanam", "fits", "dis_dep", fileoutput))

    plot_traces(post_cross$post_raw, all_models_hanam[[model_no]], file = folder_name)
    plot_histogram(post_cross$post_raw, all_models_hanam[[model_no]], file = folder_name)
    plot_titre_post(post_cross$post_raw, all_models_hanam[[model_no]], file = folder_name)
}

loadPlotFigs(1, "dis_dep/cross_sec", "cross_sec/post_cross.RData")
loadPlotFigs(2, "dis_dep/full", "full/LN_post_full.RData")
loadPlotFigs(3, "dis_dep/full", "full/SP_post_full.RData")

# Best fit of the two models according to WAIC 
save_plot_waic_cross_hanam(all_models_hanam[c(2, 3)], file = "dis_dep", FALSE)

dchodge/pilot.hcw documentation built on Feb. 2, 2022, 11:29 p.m.