Confusing behaviour when adding or removing time-varying parameters that are not calibrated

[stevencarlislewalker]

Here is how to do it properly, but it really is confusing.

We use examples from the user guide.

library(McMasterPandemic)
library(lubridate)
library(ggplot2)
library(dplyr)
sir = flexmodel(
  params = c(beta = 0.1, gamma = 0.01, N = 100),
  state = c(S = 99, I = 1, R = 0),
  start_date = "2020-03-11",
  end_date = "2020-12-01"
)
sir = (sir
 %>% add_rate("S", "I", ~ (I) * (beta) * (1/N))
 %>% add_rate("I", "R", ~ (gamma))
)
random_timevar = data.frame(
  Date = ymd(20200515),
  Symbol = 'beta',
  Value = 0.01,
  Type = 'abs'
)
random_timevar

sir_with_timevar = (sir
  %>% update_piece_wise(random_timevar)
  %>% update_error_dist(
    S ~ poisson(),
    I ~ poisson(),
    R ~ poisson()
  )
)

timevar_sims = (sir_with_timevar
  %>% simulation_history(include_initial_date = FALSE, obs_error = TRUE)
  %>% tidyr::pivot_longer(-Date, names_to = "var")
  %>% rename(date = Date)
  %>% mutate(value = round(value))
  %>% filter(var %in% c("S", "I", "R"))
)
(ggplot(timevar_sims)
  + geom_line(aes(date, value, colour = var))
)

calibrate_timevar = (random_timevar
  %>% mutate(Value = NA)
)

sir_to_cal_tv = (sir_with_timevar
  %>% update_observed(timevar_sims)
  %>% update_piece_wise(calibrate_timevar)
  %>% add_opt_params(log_beta ~ log_flat(-2))
  %>% add_opt_tv_params(tv_type = "abs"
    , log_beta ~ log_flat(-4)
  )
)

sir_to_cal_tv$do_sim_constraint = TRUE

sir_cal_tv = calibrate_flexmodel(sir_to_cal_tv)
c(
  before = pars_time_series(sir_with_timevar)$Value,
  during = pars_time_series(sir_to_cal_tv)   $Value,
  after  = pars_time_series(sir_cal_tv)      $Value
)

(sir_cal_tv
  %>% fitted
  %>% mutate(var = factor(var, topological_sort(sir_cal_tv)))
  %>% ggplot
  + facet_wrap(~var)
  + geom_line(aes(date, value))
  + geom_line(aes(date, value_fitted), colour = 'red')
)

Now that we are set up, we can add this record.

new_record = data.frame(Date = "2020-05-16", Symbol = "beta", Value = 0.1, Type = "abs")

But we need to add it in a very specific and non-intuitive way.

## so that simulation history works
tv_opt = rbind(pars_time_opt(sir_cal_tv), new_record)
sir_cal_tv2 = update_piece_wise(sir_cal_tv, tv_opt) 

## so that simulate_ensemble works
tv_spec = rbind(pars_time_spec(sir_cal_tv$model_to_calibrate), new_record)
sir_cal_tv2$model_to_calibrate = update_piece_wise(sir_cal_tv2$model_to_calibrate, tv_spec) 

Now we can do the simulations and view the different caused by the new record.

set.seed(1L)
ens = simulate_ensemble(sir_cal_tv)
set.seed(1L)
ens2 = simulate_ensemble(sir_cal_tv2)
sh = simulation_history(sir_cal_tv2)
all.equal(ens, ens2)

## should show the effect of the new record in model 2 only
bind_rows(ens, ens2, .id = "model") |> 
  filter(var == "I") |> ggplot() + geom_line(aes(Date, value, colour = model))

## should be similar, but to noise
plot(filter(ens2, var == "I")$value, sh$I)
abline(a = 0, b = 1)