mac-theobio / McMasterPandemic

SEIR+ model
GNU General Public License v3.0
20 stars 5 forks source link

Potential issue with calibration with do_hazard=TRUE #49

Open matthewcso opened 3 years ago

matthewcso commented 3 years ago

As Ben and I previously discussed during our meeting, there are some potential problems with using do_hazard and calibration...

# -----------------------------------------------------------------------------------------------
# lots of parameter-setting code, imports, etc. is be placed here - can be found in ontario/ont_cal.R
# ------------------------------------------------------------------------------------------
no_hazard <- calibrate(data=ont_all_sub
                       , base_params=params
                       , opt_pars = opt_pars
                       , time_args = list(break_dates = bd)
                       , sim_args = list(step_args=list(do_hazard=FALSE)))

do_hazard <- calibrate(data=ont_all_sub
                       , base_params=params
                       , opt_pars = opt_pars
                       , time_args = list(break_dates = bd)
                       , sim_args = list(step_args=list(do_hazard=TRUE)))

end_date <- as.Date(no_hazard$mle2@data$end_date) + 30
prediction_no_hazard <- predict(no_hazard,end_date=end_date) %>% filter(! (var %in% c('cumRep')))

end_date <- as.Date(do_hazard$mle2@data$end_date) + 30
prediction_do_hazard <- predict(do_hazard,end_date=end_date) %>% filter(! (var %in% c('cumRep')))

replace_na <- function(x) {x[is.na(x)]<- 0; x }
stopifnot(all(prediction_do_hazard$date == prediction_no_hazard$date))
stopifnot(all(prediction_do_hazard$var == prediction_no_hazard$var))
stopifnot(all(prediction_do_hazard$vtype == prediction_no_hazard$vtype))

merged <- dplyr::inner_join(prediction_do_hazard, prediction_no_hazard, by=c("date", "var", "vtype"), suffix=c("_do_hazard", "_no_hazard")) %>%
  mutate(value_do_hazard = replace_na(value_do_hazard), value_no_hazard = replace_na(value_no_hazard)) %>%
  mutate(abs_diff = abs(value_do_hazard-value_no_hazard))

ggplot_df <- merged %>% filter(var=='incidence')
ggplot(ggplot_df) +
  theme(text=element_text(size=20))+
  geom_line(aes(x=date, y=value_do_hazard, color='do_hazard=TRUE'), size=1) +
  geom_line(aes(x=date, y=value_no_hazard, color='do_hazard=FALSE'), size=1) +
  labs(x='date', y='incidence', title='do_hazard prediction discrepancies', col='')
ggsave('discrepancies.png', width=16, height=10)

discrepancies

matthewcso commented 3 years ago

Sorry, I realized that I forgot to comment on this issue. Changing the code to use params_timevar results in a completely different-looking figure, which also differs between do_hazard=TRUE and do_hazard=FALSE.

Code to replicate this figure can be found in refactor/recalibration/comparison.R.

discrepancies_break