library(odin2)
library(dust2)
library(monty)
14 Projections and counterfactuals
14.1 Introduction
We have seen how to fit odin
models with monty
to produce estimates of fitted parameters and the trajectories of the underlying model. It is common to use such outputs for either of the following:
- Projections - running the model forward from the end of fitted time series to estimate what would happen under given choices of parameters.
- Counterfactuals - estimating what would have happened had parameters been different from those estimated from or assumed in the fits. This could be run from the initial time used in fitting, or from a point in the data time series.
Generally, both of these require that you configure your likelihood object created with dust_likelihood_monty
to save additional information about the filtering process, via the save_state
or save_snapshots
inputs. As with saving fitted trajectories via the save_trajectories
input, the outputs can only be extracted after running your sampler if you specified to save them when setting up the likelihood object. It is not possible to extract such information retrospectively so it is important that you set things up correctly.
In the following we’ll illustrate how to correctly setup the likelihood object and then how to use the sampled output from fitting to run both projections and counterfactuals. We begin with some basic setup for an example of an SIS model incorporating school opening and closing, before moving onto the specific implementation for projections and then counterfactuals.
<- odin({
sis update(S) <- S - n_SI + n_IS
update(I) <- I + n_SI - n_IS
update(incidence) <- incidence + n_SI
initial(S) <- N - I0
initial(I) <- I0
initial(incidence, zero_every = 1) <- 0
<- interpolate(schools_time, schools_open, "constant")
schools <- parameter()
schools_time <- parameter()
schools_open dim(schools_time, schools_open) <- parameter(rank = 1)
<- ((1 - schools) * (1 - schools_modifier) + schools) * beta0
beta
<- 1 - exp(-beta * I / N * dt)
p_SI <- 1 - exp(-gamma * dt)
p_IS <- Binomial(S, p_SI)
n_SI <- Binomial(I, p_IS)
n_IS
<- parameter(1000)
N <- parameter(10)
I0 <- parameter()
beta0 <- parameter()
gamma <- parameter()
schools_modifier
<- data()
cases ~ Poisson(incidence)
cases })
We will fit this model to some case data.
<- read.csv("data/schools.csv")
data plot(data, pch = 19, col = "red")
We’ll assume we know the times of schools opening and closing
<- c(0, 50, 60, 120, 130, 170, 180)
schools_time <- c(1, 0, 1, 0, 1, 0, 1) schools_open
and we will fit the parameters beta0
, gamma
and schools_modifier
<- monty_packer(c("beta0", "gamma", "schools_modifier"),
packer fixed = list(schools_time = schools_time,
schools_open = schools_open))
<- monty_dsl({
prior ~ Exponential(mean = 0.3)
beta0 ~ Exponential(mean = 0.1)
gamma ~ Uniform(0, 1)
schools_modifier })
We will create our filter object, starting at time = 0
with dt = 0.25
, and setup our sampler
<- 0
time_start <- 0.25
dt <- dust_filter_create(sis, time_start = time_start, dt = dt,
filter data = data, n_particles = 200)
<- diag(c(2e-4, 2e-4, 4e-4))
vcv <- monty_sampler_random_walk(vcv) sampler
14.2 Projections
We want to run projections from the state of the model at the end of the time series. This is not saved by default so in order to do this we need to specify save_state = TRUE
when setting up our dust likelihood object. After obtaining our samples we will thin them, resulting in a sample size of 200 in this case.
<- dust_likelihood_monty(filter, packer, save_trajectories = TRUE,
likelihood save_state = TRUE)
<- likelihood + prior
posterior
<- monty_sample(posterior, sampler, 500, initial = c(0.3, 0.1, 0.5),
samples n_chains = 4)
<- monty_samples_thin(samples, burnin = 100, thinning_factor = 8)
samples <- dust_unpack_state(filter, samples$observations$trajectories) traj_fit
The saved state
can be found in samples$observations
. First we will flatten out the chains dimension of this and the fitted parameters, and unpack each vector of fitted parameters to produce a list (of length equal to the sample size) of list of parameters required by our odin model.
<- array(samples$observations$state, c(3, 200))
state <- array(samples$pars, c(3, 200))
pars <- lapply(seq_len(200), function(i) packer$unpack(pars[, i])) pars
We create a dust system using our compiled odin model and these parameters, to start at time = 150
(the final timepoint in the data time series). We set n_particles = 1
to run once per parameter set.
<- dust_system_create(sis, pars, time = 150, n_particles = 1,
sys n_groups = length(pars), dt = dt)
We then set the state of the system with the end state from fitting and simulate
dust_system_set_state(sys, state)
<- seq(150, 200)
t <- dust_system_simulate(sys, t)
res <- dust_unpack_state(sys, res) traj_projection
We can then plot our projections
<- array(traj_fit$incidence, c(150, 200))
incidence_fit matplot(data$time, incidence_fit, type = "l", col = "#00000044", lty = 1,
xlab = "Time", ylab = "Incidence", xlim = c(0, 200))
matlines(t, t(traj_projection$incidence), col = "blue")
points(data, pch = 19, col = "red")
14.3 Counterfactuals
Counterfactuals are used to examine what would have happened if certain parameters (fitted or fixed) had been different. Broadly there are two types of counterfactuals
- Counterfactuals run from the beginning of the model - when a change of parameters would result in the counterfactual diverging from the fitted trajectories from the very beginning
- Counterfactuals run from a timepoint in the data timeseries - when a change of parameters would result in the counterfactual diverging from the fitted trajectories only from that given timepoint onwards, and everything before that should be the same as the fitted trajectories.
14.3.1 Running counterfactuals from the beginning
Suppose we want to examine what would happen if the infection rate was reduced by 20%. This would affect the model from the beginning - we can just run a simulation from the start time. There is nothing that we are required to save while fitting in order to do this (although we will save trajectories from fitting to compare those to the counterfactual trajectories).
<- dust_likelihood_monty(filter, packer, save_trajectories = TRUE)
likelihood
<- likelihood + prior
posterior
<- monty_sample(posterior, sampler, 500, initial = c(0.3, 0.1, 0.5),
samples n_chains = 4)
<- monty_samples_thin(samples, burnin = 100, thinning_factor = 8)
samples <- dust_unpack_state(filter, samples$observations$trajectories) traj_fit
We will unpack the fitted parameters into a list (of length equal to the sample size) of lists of parameters
<- array(samples$pars, c(3, 200))
pars <- lapply(seq_len(200), function(i) packer$unpack(pars[, i])) pars
We will then use a function to multiply beta0
by 0.8 in each parameter list.
<- function(p) {
f $beta0 <- 0.8 * p$beta0
p
p
}<- lapply(pars, f) pars
We create a dust system setting the time to the start time used in the fitting. We set the state to the initial state as defined in the odin code, and simulate.
<- dust_system_create(sis, pars, n_particles = 1, n_groups = length(pars),
sys time = time_start, dt = dt)
dust_system_set_state_initial(sys)
<- seq(0, 150)
t <- dust_system_simulate(sys, t)
res <- dust_unpack_state(sys, res) traj_counterfactual
We can then plot the results to compare the fitted trajectories to the counterfactual ones.
<- array(traj_fit$incidence, c(150, 200))
incidence_fit matplot(data$time, incidence_fit, type = "l", col = "#00000044", lty = 1,
xlab = "Time", ylab = "Incidence")
matlines(t, t(traj_counterfactual$incidence), col = "blue")
points(data, pch = 19, col = "red")
Note that the trajectories from fitting are produced by filtering to the data, whereas the counterfactual trajectories are produced from unfiltered simulation. This is why we see much more variance in the counterfactual trajectories.
14.3.2 Running counterfactuals from a timepoint in the data time series
Suppose now we want to examine what would have happened if schools had not opened at time = 60
(as they had in our fitting) and had instead stayed closed until time = 130
. This counterfactual would diverge from the fitted trajectories at time = 60
. To implement this counterfactual we need to save a “snapshot” of the state at time = 60
.
<- dust_likelihood_monty(filter, packer, save_trajectories = TRUE,
likelihood save_snapshots = 60)
<- likelihood + prior
posterior
<- monty_sample(posterior, sampler, 500, initial = c(0.3, 0.1, 0.5),
samples n_chains = 4)
<- monty_samples_thin(samples, burnin = 100, thinning_factor = 8)
samples <- dust_unpack_state(filter, samples$observations$trajectories) traj_fit
We will unpack the fitted parameters into a list (of length equal to the sample size) of lists of parameters
<- array(samples$pars, c(3, 200))
pars <- lapply(seq_len(200), function(i) packer$unpack(pars[, i])) pars
and use a function to change schools_time
and schools_open
within each list of parameters
<- function(p) {
f $schools_time <- c(0, 50, 130, 170, 180)
p$schools_open <- c(1, 0, 1, 0, 1)
p
p
}<- lapply(pars, f) pars
We create a dust system setting the time to the start time of our counterfactual, set the state of the system using our snapshot (reshaping to remove the chains dimension first) and simulate.
<- dust_system_create(sis, pars, n_particles = 1, n_groups = length(pars),
sys time = 60, dt = dt)
<- array(samples$observations$snapshots, c(3, 200))
snapshot dust_system_set_state(sys, snapshot)
<- seq(60, 150)
t <- dust_system_simulate(sys, t)
res <- dust_unpack_state(sys, res) traj_counterfactual
We can then plot the results to compare the fitted trajectories to the counterfactual ones.
<- array(traj_fit$incidence, c(150, 200))
incidence_fit matplot(data$time, incidence_fit, type = "l", col = "#00000044", lty = 1,
xlab = "Time", ylab = "Incidence")
matlines(t, t(traj_counterfactual$incidence), col = "blue")
points(data, pch = 19, col = "red")