# -----------------------------
# SchistoDynAgg Weighted posterior distributions from ABC estimation
# Chris Hoover
# -----------------------------
library(tidyverse)
library(abc)
devtools::load_all()
# Get options passed from BASH -----------------
# Number of weighted draws from initial posterior to construct final posterior
opts <- commandArgs(TRUE)
if(length(opts) > 0){
post_samp_size <- as.numeric(opts[1])
} else {
post_samp_size <- 1000
}
# Load data and additional setup ----------------
abc_sims <- readRDS(here::here("Data/Derived/abc_fit_rejection_4cases_hrg_fixed_1e+06iterations.rds"))
yO <- readRDS(here::here("Data/Derived/ABC_yO_data.rds"))
# First remove shehia-years with failed convergence. Generally because of no infection in community
abc_successes <- sapply(1:length(abc_sims), function(s){
ifelse(length(abc_sims[[s]]) == 4, NA, s)
})
abc_sims2 <- abc_sims[na.omit(abc_successes)]
# Get posterior summaries ----------------------
abc_posteriors <- bind_rows(lapply(abc_sims2, function(abc_run){
shehia <- as.character(abc_run[[1]]) ; cat(shehia)
year <- as.integer(abc_run[[2]])
pop <- as.character(abc_run[[3]])
# Process case 1 runs
case1_post_dists <- abc_run[[4]]$dist[abc_run[[4]]$region] # Distance metrics of inputs reaching tolerance
case1_weights <- max(case1_post_dists)/(case1_post_dists+1/post_samp_size)+1 # Weights normalized based on largest distance, truncated to avoid Inf
case1_init_post <- abc_run[[4]]$unadj.values # Initial posterior as values reaching tolerance with rejection method
case1_post_samps <- sample( # Weighted sample
x = 1:nrow(case1_init_post),
size = post_samp_size,
prob = case1_weights,
replace = TRUE
)
case1_wgtd_post <- case1_init_post[case1_post_samps,] # Final weighted posterior
# Process case 2 runs
case2_post_dists <- abc_run[[5]]$dist[abc_run[[5]]$region] # Distance metrics of inputs reaching tolerance
case2_weights <- max(case2_post_dists)/(case2_post_dists+1/post_samp_size)+1 # Weights normalized based on largest distance
case2_init_post <- abc_run[[5]]$unadj.values # Initial posterior as values reaching tolerance with rejection method
case2_post_samps <- sample( # Weighted sample
x = 1:nrow(case2_init_post),
size = post_samp_size,
prob = case2_weights,
replace = TRUE
)
case2_wgtd_post <- case2_init_post[case2_post_samps,] # Final weighted posterior
# Process case 3 runs
case3_post_dists <- abc_run[[6]]$dist[abc_run[[6]]$region] # Distance metrics of inputs reaching tolerance
case3_weights <- max(case3_post_dists)/(case3_post_dists+1/post_samp_size)+1 # Weights normalized based on largest distance
case3_init_post <- abc_run[[6]]$unadj.values # Initial posterior as values reaching tolerance with rejection method
case3_post_samps <- sample( # Weighted sample
x = 1:nrow(case3_init_post),
size = post_samp_size,
prob = case3_weights,
replace = TRUE
)
case3_wgtd_post <- case3_init_post[case3_post_samps,] # Final weighted posterior
# Process case 4 runs
case4_post_dists <- abc_run[[7]]$dist[abc_run[[7]]$region] # Distance metrics of inputs reaching tolerance
case4_weights <- max(case4_post_dists)/(case4_post_dists+1/post_samp_size)+1 # Weights normalized based on largest distance
case4_init_post <- abc_run[[7]]$unadj.values # Initial posterior as values reaching tolerance with rejection method
case4_post_samps <- sample( # Weighted sample
x = 1:nrow(case4_init_post),
size = post_samp_size,
prob = case4_weights,
replace = TRUE
)
case4_wgtd_post <- case4_init_post[case4_post_samps,] # Final weighted posterior
return(list("Shehia" = list(shehia),
"Year" = list(year),
"Pop" = list(pop),
"Case1_Pstr" = list(case1_wgtd_post),
"Case2_Post" = list(case2_wgtd_post),
"Case3_Post" = list(case3_wgtd_post),
"Case4_Post" = list(case4_wgtd_post)))
}))
saveRDS(abc_posteriors,
file = here::here("Data/Derived/ABC_Weighted_Posteriors.rds"))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.