# SPARKLE Modelling Workshop
# Malaysia, 2026
##################################
## MODELLING INTERVENTIONS IN R ##
##################################Modelling interventions in R
SPARKLE workshop, Malaysia, 2026
Presenter: Michael Lydeamore, Department of Econometrics and Business Statistics, Monash University
Contributors
- With thanks to Lisa White, Nuffield Department of Medicine, Oxford University
- With thanks to Michael Meehan, Australian Institute of Tropical Health & Medicine, James Cook University
- With thanks to Wirichada Pan-Ngum, Faculty of Tropical Medicine, Mahidol University
Summary
One of the most useful aspects of a mathematical model for infectious disease transmission is its capability to explore future intervention scenarios: what could be.
We will explore some R code for Malaysia’s first major COVID-19 outbreak in 2020. This will involve the following parts:
- A whiteboard discussion of how to include interventions in a model
- Coding a model to capture the first two phases of the Movement Control Order (MCO)
- Estimating the effectiveness of the 1st wave interventions
- Exploring counterfactuals
- Practice building models to explicitly incorporate interventions of your choice in breakout groups
Modelling Interventions
As always, we should set up a new script with a good header:
Open your course RStudio Project, create a new R script named session9_interventions.R, save it with your course files and copy the header into it. Install any missing packages once with install.packages(c("dplyr", "deSolve", "ggplot2", "readr")).
Now, let’s load some packages.
library(dplyr)
library(deSolve)
library(ggplot2)We’ll be making use of Malaysia’s national 2020 outbreak data from the Ministry of Health Malaysia’s official open-data repository. Download malaysia_covid_2020.csv and save it in the same directory as your R script. Here is a reminder of the outbreak and the timing of the first two MCO phases:
Let’s load back in the data, and have a look at a plot:
first_wave <- readr::read_csv(
"malaysia_covid_2020.csv",
col_types = readr::cols(
Date = readr::col_date(),
Cases = readr::col_double(),
Cumulative_cases = readr::col_double()
)
)
ggplot(first_wave, aes(x = Date, y = Cases)) +
geom_col(fill = "blue") +
ylab("Daily cases") +
xlab("") +
ggtitle("Malaysia's 2020 COVID-19 Outbreak, Jan-Jun") +
theme_bw()
Coding Interventions
We’re going to model two intervention stages: the nationwide MCO beginning 18 March and MCO Phase 2 beginning 1 April. These dates follow the Department of Statistics Malaysia chronology. This retains the two-stage structure of the original exercise while using Malaysia-specific policy dates.
fitting_start_date <- as.Date("2020-02-28") # start of fitting window
mco_start_date <- as.Date("2020-03-18") # Start of nationwide MCO
mco_phase2_start_date <- as.Date("2020-04-01") # Start of MCO Phase 2
fitting_end_date <- as.Date("2020-05-13") # end of fitting windowWe should filter down our data to the time period we want to fit to.
outbreak <- first_wave |>
filter(Date >= fitting_start_date, Date <= fitting_end_date)One of the most important things when you’re doing an optimisation is a decent initial guess at your parameter values. Luckily, in a previous session, we got a good estimate of R0, and we know a few things about how effective the various interventions might be.
parameters <- c(
R0 = 5.235132, # from the Malaysia NM fitting session
latent_period = 5, # Rounded from Trauer, J.M., Lydeamore, M.J., Dalton, G.W. et al. Understanding how Victoria, Australia gained control of its second COVID-19 wave. Nat Commun 12, 6266 (2021). https://doi.org/10.1038/s41467-021-26558-4
infectious_period = 6, # Rounded from Trauer et al
mco_phase1_efficacy = 0.70, # guess, assuming Reff = (1-efficacy) * R0
mco_phase2_efficacy = 0.80 # guess
)Since we’re using the same model from Sessions 4 and 5, we also need the same initial conditions, state variables and time window.
# Initial conditions
Total_population <- 3.27e7 # Approximate population of Malaysia in 2020
Initial_exposed <- 0 # Simplifying assumption
Initial_infectious <- 8.997775 # from the Malaysia NM fitting session
Initial_recovered <- 0 # simplifying assumption
Initial_susceptible <- Total_population - Initial_exposed - Initial_infectious - Initial_recovered
# State variables
state <- c(
Susceptible = Initial_susceptible,
Exposed = Initial_exposed,
Infectious = Initial_infectious,
Recovered = Initial_recovered
)
# Time window
times <- seq(fitting_start_date, fitting_end_date, by = 1)We’re going to use a slightly different implementation of the SEIR model from the previous session. See if you can spot the differences.
covid_intervention <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
# Calculate the total population size
Total_population <- Susceptible + Exposed + Infectious + Recovered
# Calculate intervention efficacy
if (t < as.numeric(mco_start_date - fitting_start_date)) {
intervention_efficacy <- 0
} else if (t >= as.numeric(mco_start_date - fitting_start_date) &&
t < as.numeric(mco_phase2_start_date - fitting_start_date)) {
intervention_efficacy <- mco_phase1_efficacy
} else {
intervention_efficacy <- mco_phase2_efficacy
}
# Calculate the effective reproduction number in the presence of interventions
Reff <- (1 - intervention_efficacy) * R0 # assuming everyone starts susceptible
# Calculate the average force of infection imposed on each susceptible individual
force_of_infection <- Reff * Infectious / (Total_population * infectious_period)
# Calculate the net (instantaneous) change in each state variable
Susceptible_change <- -force_of_infection * Susceptible
Exposed_change <- force_of_infection * Susceptible - Exposed / latent_period
Infectious_change <- Exposed / latent_period - Infectious / infectious_period
Recovered_change <- Infectious / infectious_period
# Return net changes as list
return(list(
c(
Susceptible_change,
Exposed_change,
Infectious_change,
Recovered_change
)
))
})
}Now, we wrap this function in a call to ode to solve our model, and give us the outputs we are interested in.
# Wrapper function to solve model and tidy up output
solve_intervention_model <- function(initial_state = state,
model_times = times,
model_function = covid_intervention,
model_parameters = parameters) {
out <- ode(
y = initial_state,
times = as.numeric(model_times - model_times[1]),
func = model_function,
parms = model_parameters
)
# Calculate the prevalence, incidence and cumulative incidence (for comparison with data)
out <- tibble::as_tibble(out) |>
mutate(
Incidence = Exposed / model_parameters["latent_period"],
Cumulative_incidence = cumsum(Incidence) + Incidence[1],
Population = Susceptible + Exposed + Infectious + Recovered,
Prevalence = (Exposed + Infectious)/Population,
Date = model_times
)
return(tibble::as_tibble(out))
}We can solve this for our initial guess parameters, and see how it looks
# Run the model for out initial set of parameters
out_init <- solve_intervention_model()
# Make sure to add observational uncertainty
out_init <- out_init |>
mutate(
lower50 = qpois(p = 0.25, lambda = Incidence), # 50% confidence interval (i.e., 25 - 75th centiles)
upper50 = qpois(p = 0.75, lambda = Incidence),
lower95 = qpois(p = 0.025, lambda = Incidence), # 95% confidence interval (i.e., 2.5 - 97.5th centiles)
upper95 = qpois(p = 0.975, lambda = Incidence)
)
# Plot initial estimate
ggplot(outbreak) +
geom_col(aes(x = Date, y = Cases), width = 1, fill = "dodgerblue2", colour = "blue") +
geom_ribbon(
data = out_init[-1, ], aes(x = Date, ymin = lower50, ymax = upper50),
fill = "firebrick2", colour = "firebrick2", alpha = 0.8
) +
geom_ribbon(
data = out_init[-1, ], aes(x = Date, ymin = lower95, ymax = upper95),
fill = "firebrick2", colour = "firebrick2", alpha = 0.5
) +
ylab("Daily cases") +
xlab("") +
ggtitle("Malaysia's 2020 COVID-19 Outbreak") +
theme_bw()
TipNeed to catch up?
Download the code for this session so far. Keep malaysia_covid_2020.csv in the same directory before running it.
Model calibration
So we have written code that can solve our model for a given set of parameters. Now, just like we did last time, we can use maximum likelihood estimation to determine the optimal parameters for our model.
NoteParameter scaling
Optimisation algorithms rely on parameters being able to take on any value (i.e. \((-\infty, \infty)\))
Our effectiveness parameters are defined as “reductions in transmission”, which means they are bounded between 0 and 1. To put them in the form the optimiser is expecting, we’ll use the logit transform.
logit <- function(x) {
return(log(x / (1 - x)))
}
# And its inverse
inverse_logit <- function(x) {
return(exp(x) / (1 + exp(x)))
}And now we’re ready to write our negative log likelihood function.
negative_log_likelihood <- function(transformed_parameters, # assumes logit form of parameters -- per Wednesday arvo session
data = outbreak$Cases[-1],
initial_state = state,
model_times = times,
model_function = covid_intervention, # note this is our ode solver
model_parameters = parameters) {
# Back-transform model parameters (inverse logit)
model_parameters[c("mco_phase1_efficacy", "mco_phase2_efficacy")] <-
inverse_logit(transformed_parameters[c("mco_phase1_efficacy", "mco_phase2_efficacy")])
# Solve model with updated parameters
out <- solve_intervention_model(
initial_state,
model_times,
model_function,
model_parameters
)
return(-sum(dpois(
x = data,
lambda = out$Incidence[-1],
log = TRUE
))) # Poisson observation model
}The complicated part of this is transforming back and forth between the logit-transformed parameters, and the “actual” parameters that we care about. We have to make sure to transform before we start.
init_parameters <- c(logit(parameters[c(
"mco_phase1_efficacy",
"mco_phase2_efficacy"
)]))
# Calculate initial negative log-likelihood
negative_log_likelihood(init_parameters)[1] 3327.117negative_log_likelihood_initial <- -sum(dpois(
x = outbreak$Cases[-1],
lambda = out_init$Incidence[-1], # This is from the previous section!
log = TRUE
))
negative_log_likelihood_initial[1] 3327.117Now we’re in a place to use optim again, just like we’ve done previously.
optim_NM <- optim(
par = init_parameters,
fn = negative_log_likelihood,
method = "Nelder-Mead",
hessian = TRUE
)
# Check for convergence
optim_NM$convergence # 0 - converged; nonzero - failed[1] 0# Back-transform parameters
c(inverse_logit(optim_NM$par[c(
"mco_phase1_efficacy",
"mco_phase2_efficacy"
)]))mco_phase1_efficacy mco_phase2_efficacy
0.7280899 0.8692998 # Inspect the optimal negative log-likelihood
optim_NM$value[1] 633.1059This looks pretty good, and so we can be reasonably happy we’ve reached an optimal solution. Now that we have a solution, let’s plot it.
optimal_parameters <- parameters
optimal_parameters[c(
"mco_phase1_efficacy",
"mco_phase2_efficacy"
)] <- inverse_logit(optim_NM$par[c(
"mco_phase1_efficacy",
"mco_phase2_efficacy"
)])
# Solve the model given the optimal parameters and initial conditions
optimal_solution <- solve_intervention_model(
initial_state = state,
model_times = times,
model_function = covid_intervention,
model_parameters = optimal_parameters
)
# Calculate the observational confidence intervals
optimal_solution <- optimal_solution |>
mutate(
lower50 = qpois(p = 0.25, lambda = Incidence), # 50% confidence interval (i.e., 25 - 75th centiles)
upper50 = qpois(p = 0.75, lambda = Incidence),
lower95 = qpois(p = 0.025, lambda = Incidence), # 95% confidence interval (i.e., 2.5 - 97.5th centiles)
upper95 = qpois(p = 0.975, lambda = Incidence)
)
ggplot(outbreak) +
geom_col(
aes(x = Date, y = Cases),
width = 1.0,
fill = "dodgerblue2",
colour = "blue"
) +
geom_ribbon(
data = optimal_solution[-1, ],
aes(x = Date, ymin = lower50, ymax = upper50),
fill = "firebrick2", colour = "firebrick2", alpha = 0.8
) +
geom_ribbon(
data = optimal_solution[-1, ],
aes(x = Date, ymin = lower95, ymax = upper95),
fill = "firebrick2", colour = "firebrick2", alpha = 0.5
) +
ylab("Daily cases") +
xlab("") +
ggtitle("Intervention model fit (optimized)") +
geom_vline(
aes(xintercept = mco_start_date),
linetype = "dashed"
) +
geom_vline(
aes(xintercept = mco_phase2_start_date),
linetype = "dashed"
) +
annotate(
geom = "text",
label = c("MCO begins", "MCO Phase 2"),
x = c(mco_start_date, mco_phase2_start_date),
y = c(175, 175),
angle = 90,
vjust = 1.5
) +
theme_bw()
Great! The fit looks reasonable, at least visually. To report on parameter uncertainty, we use the same tricks as last time, with the square root of the diagonals of the hessian matrix.
covar_matrix <- solve(optim_NM$hessian)
std_errors <- sqrt(diag(covar_matrix))
estimates_transformed <- tibble::tibble(
estimate = optim_NM$par,
std_error = std_errors
) |>
dplyr::mutate(
lower95CI = estimate - 1.96 * std_error,
upper95CI = estimate + 1.96 * std_error
)
estimates_transformed# A tibble: 2 × 4
estimate std_error lower95CI upper95CI
<dbl> <dbl> <dbl> <dbl>
1 0.985 0.0178 0.950 1.02
2 1.89 0.0154 1.86 1.92# Back-transform estimates and confidence intervals
estimates <- inverse_logit(estimates_transformed[, -2])
estimates estimate lower95CI upper95CI
1 0.7280899 0.7211254 0.7349444
2 0.8692998 0.8658420 0.8726816
TipNeed to catch up?
Download the code for this session so far. Keep malaysia_covid_2020.csv in the same directory before running it.
Exploring counterfactuals
Now we have a model for COVID-19 transmission in Malaysia, and estimates of how effective the first two MCO stages were, we can explore what could have happened, or the counterfactuals.
What could have happened if no intervention was applied? What could have happened if the MCO had started earlier? What other interesting counterfactual scenarios could we explore? Which of the scenarios you suggested previously might be of interest to the general public and why? Would the scenarios of interest to decision makers be different? If so, how and why?