In danielward27/sbizombies: Simulation-based inference on zombie outbreak models

knitr::opts_chunk$set(echo = TRUE)

Imports

library(tidyverse)
theme_set(theme_bw())

Read in pseudo-observed data

pod_list = list()
for (i in 1:100){
  pod_file = paste0("pods/pod_theta_", i, ".txt")
  data = read.table(pod_file, col.names = c("beta", "gamma"))
  pod_list = append(pod_list, list(data))
}

true_parameters = do.call("rbind", pod_list)
true_parameters$pod_idx = 1:100
true_parameters = true_parameters %>%
  pivot_longer(cols = c(beta, gamma), names_to = "parameter", values_to = "true_parameters")

Read in predicted means

algorithms = c("rej", "mcmc", "sa")

data_list = list()
for (algorithm in algorithms){
  data_file = paste0("posterior_means/", algorithm, "_posterior_means.txt")
  alg_df = read.table(data_file, col.names = c("beta", "gamma"))
  alg_df$algorithm = algorithm
  alg_df$pod_idx = 1:100
  data_list = append(data_list, list(alg_df))
}

posterior_means = do.call("rbind", data_list)
posterior_means = posterior_means %>%
  pivot_longer(cols = c(beta, gamma), names_to = "parameter", values_to = "inferred_parameter_value")
posterior_means

Combine both data frames

df = left_join(posterior_means, true_parameters, by = c("pod_idx", "parameter"))
df$algorithm[df$algorithm == "rej"] = "ABC-Rej"
df$algorithm[df$algorithm == "mcmc"] = "ABC-MCMC"
df$algorithm[df$algorithm == "sa"] = "ABC-SA"

For plooting change factors

df$algorithm = factor(df$algorithm, levels = c("ABC-Rej", "ABC-MCMC", "ABC-SA"))

MCMC acceptance rate

acceptance_rates = c()
for (i in 1:100){
  data_file = paste0("mcmc_results/mcmc_posterior_", i, ".txt")
  df_ = read.table(data_file, col.names = c("beta", "gamma"))
  acceptance_rates = c(acceptance_rates, length(unique(df_$beta))/length(df_$beta))
}

print(mean(acceptance_rates, trim = 0.05))

Plot predicted vs observed

p1 = ggplot(df, aes(true_parameters, inferred_parameter_value)) +
  geom_point() +
  facet_grid(parameter~algorithm, scales = "fixed") +
  geom_abline(slope=1, colour = "brown") +
  labs(x = "True parameter value", y = "Inferred parameter value")
p1

ggsave("plots/true_vs_predicted.png", p1, width = 6, height = 3.5, dpi = 600)

Plot RMSE

df$absolute_error = abs(df$inferred_parameter_value - df$true_parameters)

p2 = ggplot(df, aes(y = absolute_error, x = algorithm)) +
  geom_boxplot(fill = "lightgrey") +
  facet_wrap(~parameter) +
  stat_summary(fun=mean, geom="point", shape=4, size=3, color="black", fill="red") +
  labs(x = "Algorithm", y = "Absolute error")
p2

ggsave("plots/absolute_error_boxplot.png", p2, width = 6, height = 3.5, dpi = 600)

Logprob metric

logprobs_df = read.csv("metrics/logprobs.csv")

logprobs_df = logprobs_df %>% rename("ABC-Rej" = rejection, "ABC-MCMC" = mcmc, "ABC-SA" = semi_automatic)

logprobs_df = logprobs_df %>%
  pivot_longer(everything(), names_to = "algorithm", values_to = "negative_log_prob")

logprobs_df$algorithm = factor(logprobs_df$algorithm, levels = c("ABC-Rej", "ABC-MCMC", "ABC-SA"))

p3 = ggplot(logprobs_df, aes(y = negative_log_prob, x = algorithm)) +
  geom_boxplot(fill = "lightgrey") +
  stat_summary(fun=mean, geom="point", shape=4, size=3, color="black", fill="red") +
  labs(x = "Algorithm", y = "Negative log-probability")

p3

ggsave("plots/neg_log_boxplot.png", p3, width = 6, height = 3.5, dpi = 600)

Tables for report

df %>%
  group_by(parameter, algorithm) %>%
  summarise("mean RMSE" = mean(rmse, na.rm = TRUE)) %>%
  pivot_wider(names_from = parameter, values_from = `mean RMSE`)

logprobs_df %>%
  group_by(algorithm) %>%
  summarise("mean negative log-probability" = mean(negative_log_prob, na.rm = TRUE))

Average length of the true summary statistics

n_s = c()
for (i in 1:100){
  pod_file = paste0("pods/pod_s_", i, ".txt")
  sum_stats = read.table(pod_file)
  n_s = c(n_s, nrow(sum_stats))
}
print(min(n_s))
print(max(n_s))
print(mean(n_s))