Skip to content

Commit

Permalink
Merge pull request #273 from HopkinsIDD/fix-snapshot
Browse files Browse the repository at this point in the history 
Fixing aggregation and some aesthetics in automatic postprocessing
  • Loading branch information
@saraloo
saraloo authored Aug 5, 2024
2 parents a27a733 + 7350d30 commit 1666804
Show file tree
Hide file tree
Showing 2 changed files with 172 additions and 57 deletions.
3 changes: 1 addition & 2 deletions flepimop/R_packages/inference/R/inference_slot_runner_funcs.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,8 +51,7 @@ aggregate_and_calc_loc_likelihoods <- function(
## Filter to this location
dplyr::filter(
modeled_outcome,
!!rlang::sym(obs_subpop) == location,
date %in% unique(obs$date[obs$subpop == location])
!!rlang::sym(obs_subpop) == location
) %>%
## Reformat into form the algorithm is looking for
inference::getStats(
Expand Down
226 changes: 171 additions & 55 deletions postprocessing/postprocess_snapshot.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,25 +170,71 @@ if("hosp" %in% model_outputs){
# pdf(fname, width = 20, height = 18)
# pdf(fname)
fit_stats <- names(config$inference$statistics)
subpop_names <- unique(outputs_global$hosp %>% .[,subpop])

for(i in 1:length(fit_stats)){
statistics <- purrr::flatten(config$inference$statistics[i])
cols_sim <- c("date", statistics$sim_var, "subpop","slot")
cols_data <- c("date", "subpop", statistics$data_var)
hosp_outputs_global_tmp <- copy(outputs_global$hosp)[,..cols_sim]

# aggregate based on what is in the config
df_sim <- lapply(subpop_names, function(y) {
lapply(unique(outputs_global$hosp$slot), function(x)
purrr::flatten_df(inference::getStats(
hosp_outputs_global_tmp %>% .[subpop == y & slot == x] ,
"date",
"sim_var",
stat_list = config$inference$statistics[i],
start_date = config$start_date_groundtruth,
end_date = config$end_date_groundtruth
)) %>% dplyr::mutate(subpop = y, slot = x)) %>% dplyr::bind_rows()
}) %>% dplyr::bind_rows()

df_data <- lapply(subpop_names, function(x) {
purrr::flatten_df(
inference::getStats(
gt_data %>% .[subpop == x,..cols_data],
"date",
"data_var",
stat_list = config$inference$statistics[i],
start_date = config$start_date_groundtruth,
end_date = config$end_date_groundtruth
)) %>% dplyr::mutate(subpop = x) %>%
mutate(data_var = as.numeric(data_var)) }) %>% dplyr::bind_rows()


## summarize slots
print(outputs_global$hosp %>%
.[, ..cols_sim] %>%
# print(outputs_global$hosp %>%
# .[, ..cols_sim] %>%
# .[, date := lubridate::as_date(date)] %>%
# .[, as.list(quantile(get(statistics$sim_var), c(.05, .25, .5, .75, .95), na.rm = TRUE, names = FALSE)), by = c("date", "subpop")] %>%
# ggplot() +
# geom_ribbon(aes(x = date, ymin = V1, ymax = V5), alpha = 0.1) +
# geom_ribbon(aes(x = date, ymin = V2, ymax = V4), alpha = 0.1) +
# geom_line(aes(x = date, y = V3)) +
# geom_point(data = gt_data %>%
# .[, ..cols_data],
# aes(lubridate::as_date(date), get(statistics$data_var)), color = 'firebrick', alpha = 0.1) +
# facet_wrap(~subpop, scales = 'free', ncol = gg_cols) +
# labs(x = 'date', y = fit_stats[i], title = statistics$sim_var) +
# theme_classic()
# )
print(
df_sim %>%
setDT() %>%
.[, date := lubridate::as_date(date)] %>%
.[, as.list(quantile(get(statistics$sim_var), c(.05, .25, .5, .75, .95), na.rm = TRUE, names = FALSE)), by = c("date", "subpop")] %>%
.[, as.list(quantile(sim_var, c(.05, .25, .5, .75, .95), na.rm = TRUE, names = FALSE)), by = c("date", "subpop")] %>%
setnames(., paste0("V", 1:5), paste0("q", c(.05,.25,.5,.75,.95))) %>%
ggplot() +
geom_ribbon(aes(x = date, ymin = V1, ymax = V5), alpha = 0.1) +
geom_ribbon(aes(x = date, ymin = V2, ymax = V4), alpha = 0.1) +
geom_line(aes(x = date, y = V3)) +
geom_point(data = gt_data %>%
.[, ..cols_data],
aes(lubridate::as_date(date), get(statistics$data_var)), color = 'firebrick', alpha = 0.1) +
facet_wrap(~subpop, scales = 'free', ncol = gg_cols) +
labs(x = 'date', y = fit_stats[i], title = statistics$sim_var) +
geom_ribbon(aes(x = date, ymin = q0.05, ymax = q0.95), alpha = 0.1) +
geom_ribbon(aes(x = date, ymin = q0.25, ymax = q0.75), alpha = 0.1) +
geom_line(aes(x = date, y = q0.5)) +
# if inference, plot gt along side
geom_point(data = df_data,
aes(lubridate::as_date(date), data_var), color = 'firebrick', alpha = 0.2, size=1) +
facet_wrap(~subpop, scales = 'free') +
labs(x = 'date', y = fit_stats[i]) +
theme_classic()
)

Expand All @@ -206,23 +252,44 @@ if("hosp" %in% model_outputs){
# )

## plot cumulatives
print(outputs_global$hosp %>%
.[, ..cols_sim] %>%
.[, date := lubridate::as_date(date)] %>%
.[, csum := cumsum(get(statistics$sim_var)), by = .(subpop, slot)] %>%
.[, as.list(quantile(csum, c(.05, .25, .5, .75, .95), na.rm = TRUE, names = FALSE)), by = c("date", "subpop")] %>%
ggplot() +
geom_ribbon(aes(x = date, ymin = V1, ymax = V5), alpha = 0.1) +
geom_ribbon(aes(x = date, ymin = V2, ymax = V4), alpha = 0.1) +
geom_line(aes(x = date, y = V3)) +
geom_point(data = gt_data %>%
.[, ..cols_data] %>%
.[, csum := cumsum(replace_na(get(statistics$data_var), 0)) , by = .(subpop)]
,
aes(lubridate::as_date(date), csum), color = 'firebrick', alpha = 0.1) +
facet_wrap(~subpop, scales = 'free', ncol = gg_cols) +
labs(x = 'date', y = fit_stats[i], title = paste0("cumulative ", statistics$sim_var)) +
theme_classic()
print(
df_sim %>%
setDT() %>%
.[, date := lubridate::as_date(date)] %>%
.[, .(date, subpop, sim_var, slot)] %>%
data.table::melt(., id.vars = c("date", "slot", "subpop")) %>%
# dplyr::arrange(subpop, slot, date) %>%
.[, csum := cumsum(value), by = .(slot, subpop, variable)] %>%
.[, as.list(quantile(csum, c(.05, .25, .5, .75, .95), na.rm = TRUE, names = FALSE)), by = c("date", config$subpop_setup$subpop)] %>%
setnames(., paste0("V", 1:5), paste0("q", c(.05,.25,.5,.75,.95))) %>%
ggplot() +
geom_ribbon(aes(x = date, ymin = q0.05, ymax = q0.95), alpha = 0.1) +
geom_ribbon(aes(x = date, ymin = q0.25, ymax = q0.75), alpha = 0.1) +
geom_line(aes(x = date, y = q0.5)) +
geom_point(data = df_data %>% setDT() %>%
.[, csum := cumsum(data_var) , by = .(subpop)],
aes(lubridate::as_date(date), csum), color = 'firebrick', alpha = 0.2, size=1) +
facet_wrap(~subpop, scales = 'free') +
# facet_wrap(~get(subpop), scales = 'free') +
labs(x = 'date', y = paste0("cumulative ", fit_stats[i])) +
theme_classic()
# outputs_global$hosp %>%
# .[, ..cols_sim] %>%
# .[, date := lubridate::as_date(date)] %>%
# .[, csum := cumsum(get(statistics$sim_var)), by = .(subpop, slot)] %>%
# .[, as.list(quantile(csum, c(.05, .25, .5, .75, .95), na.rm = TRUE, names = FALSE)), by = c("date", "subpop")] %>%
# ggplot() +
# geom_ribbon(aes(x = date, ymin = V1, ymax = V5), alpha = 0.1) +
# geom_ribbon(aes(x = date, ymin = V2, ymax = V4), alpha = 0.1) +
# geom_line(aes(x = date, y = V3)) +
# geom_point(data = gt_data %>%
# .[, ..cols_data] %>%
# .[, csum := cumsum(replace_na(get(statistics$data_var), 0)) , by = .(subpop)]
# ,
# aes(lubridate::as_date(date), csum), color = 'firebrick', alpha = 0.1) +
# facet_wrap(~subpop, scales = 'free', ncol = gg_cols) +
# labs(x = 'date', y = fit_stats[i], title = paste0("cumulative ", statistics$sim_var)) +
# theme_classic()
)

}
Expand All @@ -231,47 +298,96 @@ if("hosp" %in% model_outputs){

## hosp by highest and lowest llik

if("llik" %in% model_outputs){
llik_rank <- copy(outputs_global$llik) %>%
.[, .SD[order(ll)], subpop]
high_low_llik <- rbindlist(list(data.table(llik_rank, key = "subpop") %>%
.[, head(.SD,5), by = subpop] %>%
.[, llik_bin := "top"],
data.table(llik_rank, key = "subpop") %>%
.[, tail(.SD,5), by = subpop]%>%
.[, llik_bin := "bottom"])
)
}

fname <- paste0("pplot/hosp_by_llik_mod_outputs_", opt$run_id,".pdf")
# pdf_dims <- data.frame(width = gg_cols*2, length = num_nodes/gg_cols * 2)
# pdf(fname, width = pdf_dims$width, height = pdf_dims$length)
pdf(fname, width = 20, height = 20)
pdf(fname, width = 20, height = 10)

for(i in 1:length(fit_stats)){
statistics <- purrr::flatten(config$inference$statistics[i])
cols_sim <- c("date", statistics$sim_var, "subpop","slot")
cols_data <- c("date", "subpop", statistics$data_var)
hosp_outputs_global_tmp <- copy(outputs_global$hosp)[,..cols_sim]

if("llik" %in% model_outputs){
llik_rank <- copy(outputs_global$llik) %>%
.[, .SD[order(ll)], subpop]
high_low_llik <- rbindlist(list(data.table(llik_rank, key = "subpop") %>%
.[, head(.SD,5), by = subpop] %>%
.[, llik_bin := "top"],
data.table(llik_rank, key = "subpop") %>%
.[, tail(.SD,5), by = subpop]%>%
.[, llik_bin := "bottom"])
)
# high_low_hosp_llik <- copy(outputs_global$hosp) %>%
# .[high_low_llik, on = c("slot", "subpop"), allow.cartesian = TRUE]

high_low_hosp_llik <- copy(outputs_global$hosp) %>%
.[high_low_llik, on = c("slot", "subpop"), allow.cartesian = TRUE]
# aggregate simulation output and data by time based on what is in the config
df_sim <- lapply(subpop_names, function(y) {
lapply(unique(outputs_global$hosp$slot), function(x)
purrr::flatten_df(inference::getStats(
hosp_outputs_global_tmp %>% .[subpop == y & slot == x] ,
"date",
"sim_var",
stat_list = config$inference$statistics[i],
start_date = config$start_date_groundtruth,
end_date = config$end_date_groundtruth
)) %>% dplyr::mutate(subpop = y, slot = x)) %>% dplyr::bind_rows()
}) %>% dplyr::bind_rows() %>% setDT()

df_data <- lapply(subpop_names, function(x) {
purrr::flatten_df(
inference::getStats(
gt_data %>% .[subpop == x,..cols_data],
"date",
"data_var",
stat_list = config$inference$statistics[i],
start_date = config$start_date_groundtruth,
end_date = config$end_date_groundtruth
)) %>% dplyr::mutate(subpop = x) %>%
dplyr::mutate(data_var = as.numeric(data_var)) %>%
dplyr::mutate(date = lubridate::as_date(date)) }) %>%
dplyr::bind_rows() %>% setDT()

# add likelihood ranking to simulation output
high_low_hosp_llik <- df_sim %>%
.[high_low_llik, on = c("slot", "subpop"), allow.cartesian=TRUE] %>% # right join by "on" variables
.[subpop != "Total"]

hosp_llik_plots <- lapply(unique(high_low_hosp_llik %>% .[, subpop]),
function(e){
high_low_hosp_llik %>%
.[, date := lubridate::as_date(date)] %>%
.[subpop == e] %>%
ggplot() +
geom_line(aes(lubridate::as_date(date), get(statistics$data_var),
group = slot, color = ll))+#, linetype = llik_bin)) +
# scale_linetype_manual(values = c(1, 2), name = "likelihood\nbin") +
.[, date := lubridate::as_date(date)] %>%
ggplot() +
geom_line(aes(x = date, y = sim_var, group = slot, color = ll)) +
scale_linetype_manual(values = c(1, 2), name = "likelihood\nbin") +
scale_color_viridis_c(option = "D", name = "log\nlikelihood") +
geom_point(data = gt_data %>%
.[, ..cols_data] %>%
.[subpop == e] ,
aes(lubridate::as_date(date), get(statistics$data_var)), color = 'firebrick', alpha = 0.1) +
facet_wrap(~subpop, scales = 'free', ncol = gg_cols) +
labs(x = 'date', y = fit_stats[i]) + #, title = paste0("top 5, bottom 5 lliks, ", statistics$sim_var)) +
theme_classic() +
guides(linetype = 'none')
geom_point(data = df_data %>% .[subpop == e],
aes(lubridate::as_date(date), data_var), color = 'firebrick', alpha = 0.2, size=1) +
facet_wrap(~subpop, scales = 'free') +
labs(x = 'date', y = fit_stats[i]) +
theme_classic() +
theme(legend.key.size = unit(0.2, "cm"))
# high_low_hosp_llik %>%
# .[, date := lubridate::as_date(date)] %>%
# .[subpop == e] %>%
# ggplot() +
# geom_line(aes(lubridate::as_date(date), get(statistics$data_var),
# group = slot, color = ll))+#, linetype = llik_bin)) +
# # scale_linetype_manual(values = c(1, 2), name = "likelihood\nbin") +
# scale_color_viridis_c(option = "D", name = "log\nlikelihood") +
# geom_point(data = gt_data %>%
# .[, ..cols_data] %>%
# .[subpop == e] ,
# aes(lubridate::as_date(date), get(statistics$data_var)), color = 'firebrick', alpha = 0.1) +
# facet_wrap(~subpop, scales = 'free', ncol = gg_cols) +
# labs(x = 'date', y = fit_stats[i]) + #, title = paste0("top 5, bottom 5 lliks, ", statistics$sim_var)) +
# theme_classic() +
# guides(linetype = 'none')
}
)

Expand All @@ -290,8 +406,8 @@ if("hnpi" %in% model_outputs){

gg_cols <- 4
num_nodes <- length(unique(outputs_global$hnpi %>% .[,subpop]))
pdf_dims <- data.frame(width = gg_cols*3, length = num_nodes/gg_cols * 2)
#pdf_dims <- data.frame(width = 20, length = 5)
# pdf_dims <- data.frame(width = gg_cols*3, length = num_nodes/gg_cols * 2)
pdf_dims <- data.frame(width = 20, length = 10)

fname <- paste0("pplot/hnpi_mod_outputs_", opt$run_id,".pdf")
pdf(fname, width = pdf_dims$width, height = pdf_dims$length)
Expand Down

0 comments on commit 1666804

Please sign in to comment.