Formative Assessment
2026-03-11
Ebola
A growing number of Ebola cases have been reported in a small town. The following table shows the number of reported cases each day for the last 11 days.
Day 1 2 3 4 5 6 7 8 9 10 11
Cases 1 2 2 3 3 6 4 6 13 13 20
You have been given additional information that will help you to model the disease
The population of the town is 100,000
An inf#ected person will infect others at a rate of 3 per day
The latent period of the disease is 10 days
The average infectious duration is assumed to be 2 days, since after 2 days most infected people will be too ill to have contact with others
Some individuals will die from infection while others will recover and become immune. In either case we can assume they are “removed” from the system
The mayor of the town is considering two possible ways to keep the disease under control
Force anyone suspected of having latent Ebola into quarantine. Each day, 20% of infected people who are in the latent period of their infection are discovered and moved into quarantine. The quarantine period is 14 days.
Part 1. Introduce the problem
Give a brief description of the disease, who it affects, and where it is a problem. Talk about the factors relevant to how it spreads and how it can be controlled. Use online sources like the World Health Organisation, the NHS and the various centres for disease control.
# Creating data frame
day <- 1:11
cases <- c(1,2,2,3,3,6,4,6,13,13,20)
data <- data.frame(day, cases)A plot of the case numbers over time.
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.5.2ggplot(data, aes(x = day, y = cases)) +
geom_point(size = 4) +
geom_line(linewidth = 1.2, color="red") +
labs(
title = "Daily Ebola Cases Over Time",
x = "Day",
y = "Number of Cases"
) +
scale_x_continuous(limits = c(1,12), breaks = 1:12) +
theme_minimal()
A plot of the line of best fit on the same axes as the data and the rate of exponential growth (the slope)
data$log_cases <- log(data$cases)
model <- lm(log_cases ~ day, data = data)
summary(model)##
## Call:
## lm(formula = log_cases ~ day, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.40048 -0.15546 -0.01162 0.18206 0.27549
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.10674 0.15410 -0.693 0.506
## day 0.27050 0.02272 11.905 8.24e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2383 on 9 degrees of freedom
## Multiple R-squared: 0.9403, Adjusted R-squared: 0.9337
## F-statistic: 141.7 on 1 and 9 DF, p-value: 8.238e-07r <- coef(model)[2]
data$fit_cases <- exp(predict(model))
library(ggplot2)
ggplot(data, aes(x = day, y = cases)) +
geom_point(size = 4) +
geom_line(linewidth = 1.2) +
geom_line(aes(y = fit_cases), color = "red", linewidth = 1.3) +
labs(
title = "Ebola Outbreak: Exponential Growth Fit",
x = "Day",
y = "Number of Cases"
) +
scale_x_continuous(limits = c(1,12), breaks = 1:12) +
annotate("text",
x = 3,
y = max(data$cases),
label = paste("Growth rate r =", round(r,3)),
size = 5,
color = "red") +
theme_minimal()
# Load the required package
library(deSolve)## Warning: package 'deSolve' was built under R version 4.5.2# Define the SEIR model
ODE_system <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
# Force of infection
lambda <- beta * I / N
# Differential equations
dS <- -lambda * S
dE <- lambda * S - sigma * E
dI <- sigma * E - gamma * I
dR <- gamma * I
list(c(dS, dE, dI, dR))
})
}# Define initial conditions
state <- c(
S = 99999,
E = 0,
I = 1,
R = 0
)# Define parameters
parameters <- c(
beta = 3, # transmission rate per day
sigma = 1/10, # latent -> infectious
gamma = 1/2, # infectious -> removed
N = 100000
)# Define time steps
times <- seq(0, 100, by = 1)# Solve with deSolve
out <- ode(
y = state,
times = times,
func = ODE_system,
parms = parameters
)
out <- as.data.frame(out)
head(out)## time S E I R
## 1 0 99999.00 0.000000 1.0000000 0.0000000
## 2 1 99996.52 2.350496 0.7136612 0.4128424
## 3 2 99994.46 4.091527 0.6934755 0.7570302
## 4 3 99992.22 5.832185 0.8154804 1.1295925
## 5 4 99989.47 7.906637 1.0389306 1.5892236
## 6 5 99985.89 10.562292 1.3609281 2.1848420# Plot the epidemic curves
plot(out$time, out$S, type="l", col="blue", ylim=c(0,100000),
xlab="Time (days)", ylab="Population",
main="SEIR Model of Ebola Outbreak")
lines(out$time, out$E, col="orange")
lines(out$time, out$I, col="red")
lines(out$time, out$R, col="green")
legend("right",
legend=c("Susceptible","Exposed","Infectious","Removed"),
col=c("blue","orange","red","green"),
lty=1)
library(ggplot2)
library(tidyr)## Warning: package 'tidyr' was built under R version 4.5.2out_long <- pivot_longer(out,
cols = c(S,E,I,R),
names_to = "Compartment",
values_to = "Count")
ggplot(out_long, aes(x=time, y=Count, colour=Compartment)) +
geom_line(size=1.2) +
labs(title="SEIR Model Simulation of Ebola Outbreak",
x="Time (days)",
y="Population") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Policy A: social restrictions
The mayor of the town is considering two possible ways to keep the disease under control a) Impose social restrictions on the whole town. This will reduce the transmission rate from 3 to 2
# Load the required package
library(deSolve)# Define the SEIR model
ODE_system_restriction <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
# Force of infection
lambda <- beta * I / N
# Differential equations
dS <- -lambda * S
dE <- lambda * S - sigma * E
dI <- sigma * E - gamma * I
dR <- gamma * I
list(c(dS, dE, dI, dR))
})
}
# Define initial conditions
state <- c(
S = 99999,
E = 0,
I = 1,
R = 0
)
# Define parameters
parameters_restriction <- c(
beta = 2,
sigma = 1/10,
gamma = 1/2,
N = 100000
)
# Define time steps
times <- seq(0, 100, by = 1)
# Solve with deSolve
out <- ode(
y = state,
times = times,
func = ODE_system_restriction,
parms = parameters_restriction
)
out <- as.data.frame(out)
head(out)## time S E I R
## 1 0 99999.00 0.000000 1.0000000 0.0000000
## 2 1 99997.37 1.541590 0.6773451 0.4063216
## 3 2 99996.14 2.563665 0.5776240 0.7137862
## 4 3 99994.99 3.417353 0.5888597 1.0021109
## 5 4 99993.75 4.276110 0.6622722 1.3128577
## 6 5 99992.31 5.235707 0.7781906 1.6713959library(ggplot2)
library(tidyr)
out_long <- pivot_longer(out,
cols = c(S,E,I,R),
names_to = "Compartment",
values_to = "Count")
ggplot(out_long, aes(x=time, y=Count, colour=Compartment)) +
geom_line(size=1.2) +
labs(title="SEIR Model Simulation of Ebola Outbreak",
x="Time (days)",
y="Population") +
theme_minimal()
Policy B: quarantine of latent cases
- Force anyone suspected of having latent Ebola into quarantine. Each day, 20% of infected people who are in the latent period of their infection are discovered and moved into quarantine. The quarantine period is 14 days.
library(deSolve)
# Define the SEIR model
ODE_system_quarantine <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
lambda <- beta * I / N
dS <- -lambda * S
dE <- lambda * S - sigma * E - q * E
dI <- sigma * E - gamma * I
dQ <- q * E - omega * Q
dR <- gamma * I + omega * Q
list(c(dS, dE, dI, dQ, dR))
})
}
# Define initial state
state_quarantine <- c(
S = 99999,
E = 0,
I = 1,
Q = 0,
R = 0
)
# Define parameters
parameters_quarantine <- c(
beta = 3,
sigma = 1/10,
gamma = 1/2,
q = 0.2,
omega = 1/14,
N = 100000
)
# Time step
times <- seq(0, 100, by = 1)
# 5. Solve the model
out_quarantine <- ode(
y = state_quarantine,
times = times,
func = ODE_system_quarantine,
parms = parameters_quarantine
)
out_quarantine <- as.data.frame(out_quarantine)
library(ggplot2)
library(tidyr)
out_long <- pivot_longer(
out_quarantine,
cols = c(S, E, I, Q, R),
names_to = "Compartment",
values_to = "Population"
)
ggplot(out_long, aes(x = time, y = Population, color = Compartment)) +
geom_line(linewidth = 1.3) +
labs(
title = "SEIQR Model with Quarantine Policy",
x = "Time (days)",
y = "Population"
) +
theme_minimal() +
theme(legend.position = "right")
ODE_system_quarantine <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
lambda <- beta * I / N
dS <- -lambda * S
dE <- lambda * S - sigma * E - q * E
dI <- sigma * E - gamma * I
dQ <- q * E - omega * Q
dR <- gamma * I + omega * Q
list(c(dS, dE, dI, dQ, dR))
})
}
times <- seq(0, 100, by = 1)
parameters_base <- c(
beta = 3,
sigma = 1/10,
gamma = 1/2,
N = 100000
)
parameters_social <- c(
beta = 2,
sigma = 1/10,
gamma = 1/2,
N = 100000
)
parameters_quarantine <- c(
beta = 3,
sigma = 1/10,
gamma = 1/2,
q = 0.2, # MUST exist
omega = 1/14,
N = 100000
)
library(deSolve)
out_quarantine <- ode(
y = state_quarantine,
times = times,
func = ODE_system_quarantine,
parms = parameters_quarantine
)
out_quarantine <- as.data.frame(out_quarantine)
parameters_quarantine## beta sigma gamma q omega N
## 3.000000e+00 1.000000e-01 5.000000e-01 2.000000e-01 7.142857e-02 1.000000e+05out_base <- ode(y = state, times = times, func = ODE_system, parms = parameters_base)
out_social <- ode(y = state, times = times, func = ODE_system, parms = parameters_social)
out_quarantine <- ode(y = state_quarantine, times = times, func = ODE_system_quarantine, parms = parameters_quarantine)
out_base <- as.data.frame(out_base)
out_social <- as.data.frame(out_social)
out_quarantine <- as.data.frame(out_quarantine)
compare_data <- data.frame(
time = out_base$time,
Baseline = out_base$I,
Social_Restrictions = out_social$I,
Quarantine = out_quarantine$I
)
library(tidyr)
compare_long <- pivot_longer(
compare_data,
cols = c(Baseline, Social_Restrictions, Quarantine),
names_to = "Scenario",
values_to = "Infectious"
)
library(ggplot2)
ggplot(compare_long, aes(x = time, y = Infectious, color = Scenario)) +
geom_line(linewidth = 1.4) +
labs(
title = "Comparison of Ebola Control Policies",
x = "Time (days)",
y = "Number of Infectious Individuals"
) +
scale_color_manual(values = c(
"Baseline" = "red",
"Social_Restrictions" = "blue",
"Quarantine" = "purple"
)) +
theme_minimal() +
theme(
legend.position = "right",
plot.title = element_text(hjust = 0.5)
)