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

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)

Plot sensible things

load(file = here::here("data", paste0("modelinfo_cross_", hcwpre$study_name_short, ".RDS"))) #all_models_hcw_pre_cross
load(file = here::here("outputs", "hcw_pre", "fits", "dis_dep", "cross_sec", "trimmed_post_cross.Rdata"))
all_post_test_cross <- post_cross$post_raw

chain_short <- as.data.frame(all_post_test_cross[[6]]$theta_chain) %>% filter(chain_no == 1)
infection_histories <- all_post_test_cross[[6]]$inf_chain %>% filter(chain_no == 1)
titre_dat <- all_models_hcw_pre_cross[[6]]$create_post$titre_data %>% filter(sample_time == 0)
vaccination_histories <- all_models_hcw_pre_cross[[6]]$create_post$vaccination_histories

antigenic_map <- all_models_hcw_pre_cross[[6]]$create_post$antigenic_map
strain_isolation_times <-  all_models_hcw_pre_cross[[6]]$create_post$strain_isolation_times
par_tab <- all_models_hcw_pre_cross[[6]]$create_post$par_tab

vaccination_histories <- all_models_hcw_pre_cross[[6]]$create_post$vac_history
custom_ab_kin_func <- all_models_hcw_pre_cross[[6]]$custom_ab_kin_func
custom_ab_kin_func_1 <- all_models_hcw_pre_cross[[1]]$custom_ab_kin_func

custom_antigenic_maps_func <- all_models_hcw_pre_cross[[6]]$custom_antigenic_maps_func

expand_titredat <- FALSE
nsamp <- 100
titre_before_infection <- FALSE
measurement_indices_by_time <- NULL
individuals <- unique(titre_dat$individual)
mu_indices <- NULL

plot_infection_histories_long(
                        chain = chain_short,
                        infection_histories = infection_histories,
                        vaccination_histories = vaccination_histories,
                        titre_dat = titre_dat,
                        individuals = 34,
                        antigenic_map = antigenic_map,
                        strain_isolation_times = strain_isolation_times,
                        par_tab = par_tab,
                        custom_ab_kin_func = custom_ab_kin_func,
                        custom_antigenic_maps_func = custom_antigenic_maps_func)

load(file = here::here("data", paste0("modelinfo_full_", hcwpre$study_name_short, ".RDS"))) #all_models_hcw_pre_cross
load(file = here::here("outputs", "hcw_pre", "fits", "dis_dep", "full", "trimmed_post_full.Rdata"))

all_post_test_full_raw <- post_full$post_raw
all_post_test_full <- post_full$post_trim
post_out_trim <- all_post_test_full[[1]][[2]]

mu_short_vac_v <- post_out_trim$mu_short_vac
wane_vac_v <- post_out_trim$wane_vac
mu_v <- post_out_trim$mu

sero_data <- all_models_hcw_pre_full[[1]]$create_post$titre_data 
sero_data_base <- sero_data %>% 
        group_by(virus, sample_time) %>%
        summarise(titre_mean = mean(titre, rm.na = TRUE)) %>%
        pivot_wider(names_from = sample_time, values_from = titre_mean) %>%
        mutate(day30 = (`30` - `0`), day180 = (`180` - `0`)) %>%
        select(virus, base = `0`, day30, day180) %>%
        pivot_longer(c(day30, day180), names_to = "time", values_to = "abs") %>%
        mutate(time = case_when(time == "day30"~"30", time == "day180"~"180"))

dataplot <- sero_data_base %>% filter(virus == 2016 * 12) %>% mutate(time = as.numeric(time) / 30)

func_linear <- function(x, boost, wane, mu, mu_vac) boost - wane * x + (mu * mu_vac)
ggplot() +
    xlim(0, 12) +
    lapply(500:550,
        function(i) geom_function(
            fun = func_linear, 
            args = list(boost = mu_short_vac_v[i], wane = wane_vac_v[i], mu = mu_v[i], mu_vac = mu_vac_v[i]), 
            alpha= 0.7, size = 0.2, color = "red")
    ) + 
    geom_point(data = dataplot, aes(x = time, y = abs), size = 3) +
    labs(x = "Time (months)", y = "Titre increase")













all_post_test_full <- post_full$post_raw
chain_full <- as.data.frame(all_post_test_full_raw[[pos]]$theta_chain) %>% filter(chain_no == 2)
infection_histories <- all_post_test_full_raw[[pos]]$inf_chain %>% filter(chain_no == 2)
titre_dat <- all_models_hcw_pre_full[[pos]]$create_post$titre_data 
vaccination_histories <- all_models_hcw_pre_full[[pos]]$create_post$vaccination_histories

antigenic_map <- all_models_hcw_pre_full[[pos]]$create_post$antigenic_map
strain_isolation_times <-  all_models_hcw_pre_full[[pos]]$create_post$strain_isolation_times
par_tab <- all_models_hcw_pre_full[[pos]]$create_post$par_tab

vaccination_histories <- all_models_hcw_pre_full[[pos]]$create_post$vac_history
custom_ab_kin_func <- all_models_hcw_pre_full[[pos]]$custom_ab_kin_func
custom_ab_kin_func_1 <- all_models_hcw_pre_full[[pos]]$custom_ab_kin_func

custom_antigenic_maps_func <- all_models_hcw_pre_full[[pos]]$custom_antigenic_maps_func

expand_titredat <- FALSE
nsamp <- 100
titre_before_infection <- FALSE
measurement_indices_by_time <- NULL
individuals <- unique(titre_dat$individual)
mu_indices <- NULL

y <- get_titre_predictions(
    chain = chain_full,
    infection_histories = infection_histories,
    vaccination_histories = vaccination_histories,
    titre_dat = titre_dat,
    individuals = 1:49,
    antigenic_map = antigenic_map,
    strain_isolation_times = strain_isolation_times,
    par_tab = par_tab,
    custom_ab_kin_func = custom_ab_kin_func,
    custom_antigenic_maps_func = custom_antigenic_maps_func,
    expand_titredat = TRUE,
    nsamp = 100
)

plot_infection_histories_long(
    chain = chain_full,
    infection_histories = infection_histories,
    vaccination_histories = vaccination_histories,
    titre_dat = titre_dat,
    individuals = 5:10,
    antigenic_map = antigenic_map,
    strain_isolation_times = strain_isolation_times,
    par_tab = par_tab,
    custom_ab_kin_func = custom_ab_kin_func,
    custom_antigenic_maps_func = custom_antigenic_maps_func,
    size_plot = 5
)

chain_full$mu_vac %>% mean
chain_full$mu_short_vac %>% mean

Working!!!

sum_data <- titre_dat %>% group_by(sample_time, titre) %>% summarise(n = n()) 
weighted_sum_data <- sum_data %>% summarise(w_mean = weighted.mean(titre, n, na.rm = TRUE))

#predictions
titre_pred <- y$predicted_observations
titre_pred <- y$predictions

sum_data <- titre_dat %>% group_by(sample_time, titre) %>% summarise(n = n()) 
weighted_sum_data <- sum_data %>% summarise(w_mean = weighted.mean(titre, n, na.rm = TRUE))

sum_pred <- titre_pred %>% group_by(samples, median) %>% summarise(n = n()) 
weighted_sum_pred <- sum_pred %>% summarise(w_mean = weighted.mean(median, n, na.rm = TRUE))

Post plots

all_post_test_full_trim <- post_full$post_trim
postplots_long <- all_post_test_full_trim[[5]] %>% bind_rows
dec <- 1:12 %>% map(~
    (postplots_long[["mu_short_vac"]] * (1.0 - ((1 / 12 + postplots_long[["wane_vac"]]) * .x )))
) 

df_dec <- dec %>% setNames(as.character(c(1:12))) %>%
    as.data.frame %>% pivot_longer(c(1:12), values_to = "titre", names_to = "time") %>%
    mutate(time = str_extract(time, "[0-9,.]+")) %>%
    mutate(time = factor(time, levels = as.character(c(1:12))))
df_dec <- df_dec[df_dec$titre > 0, ]

sero_data <- hcwpre$sero_data
sum_data <- sero_data %>% group_by(sample_time, titre_value) %>% summarise(n = n()) 
weighted_sum_data <- sum_data %>% summarise(w_mean = weighted.mean(titre_value, n, na.rm = TRUE)) %>%
    mutate(sample_time = sample_time / 30)

ggplot() +
        geom_point(data = sum_data, aes(x = sample_time, y = titre_value, size = n), fill = "black", alpha = 0.5) +
        geom_line(data = weighted_sum_data, aes(x = sample_time, y = w_mean), color = "red", size = 1) +
        geom_point(data = weighted_sum_data, aes(x = sample_time, y = w_mean), shape = 21, fill = "red", size = 5) + 
        labs(x = "Time post vaccination", y = "Titre value", size = "Number of samples")

df_dec %>%
    ggplot() + 
        geom_boxplot(aes(x = time, y = titre)) +
        labs(x = "Months post vaccination") +
        geom_point(data = weighted_sum_data, aes(x = sample_time, y = w_mean), shape = 8, color = "red", size = 5) +
        geom_hline(yintercept = postplots_long[["mu"]] * postplots_long[["mu_vac"]] %>% mean)

antigenic_map_melted <- c(melt_antigenic_coords(antigenic_map[, c("x_coord", "y_coord")]))

chain$sigma1 %>% mean
s <- 0.01
lapply(antigenic_map_melted, function(x) max(1 - x * s, 0)) %>% unlist