The following code reproduces our analysis of data from Namibia as described in our paper.
full_dat <- moire::namibia_data
epi_dat <- full_dat |>
dplyr::select(sample_id, HealthFacility, HealthDistrict, Region, Country) |>
dplyr::distinct()
all_hfs <- epi_dat |>
dplyr::pull(HealthFacility) |>
unique()
verbose <- FALSE
allow_relatedness <- TRUE
burnin <- 5e3
num_samples <- 1e4
r_alpha <- 1
r_beta <- 1
eps_pos_alpha <- 1
eps_pos_beta <- 1
eps_neg_alpha <- 1
eps_neg_beta <- 1
num_threads <- parallelly::availableCores() - 1
for (hf in all_hfs) {
hf_dat <- full_dat |>
dplyr::filter(HealthFacility == hf) |>
dplyr::select(sample_id, locus, allele) |>
moire::load_long_form_data()
hf_res <- moire::run_mcmc(
hf_dat, hf_dat$is_missing,
allow_relatedness = allow_relatedness,
burnin = burnin, samples_per_chain = num_samples,
pt_chains = 40, pt_num_threads = num_threads, thin = 10,
verbose = verbose, adapt_temp = TRUE, r_alpha = r_alpha, r_beta = r_beta,
eps_pos_alpha = eps_pos_alpha, eps_pos_beta = eps_pos_beta,
eps_neg_alpha = eps_neg_alpha, eps_neg_beta = eps_neg_beta
)
# create an output directory
dir.create("mcmc_output", showWarnings = FALSE)
# format path name
hf <- gsub(" ", "_", hf)
# save the results
saveRDS(hf_res, file.path("mcmc_output", paste0(hf, ".rds")))
}