4. Montecarlo sampling and the Reed-Frost model under intervention#
We saw in the last chapter that the Reed-frost model is a powerful tool for understanding how an outbreak unfolds and the role of uncertainty in outbreak dynamics. However, the Reed-Frost model is flexible. We will see how can add timed interventions to this model. In addition, we will learn about a powerful tool called Monte carlo sampling that can be used to estimate difficult-to-compute probabilities.
Lets look at how we can model public health interventions during the course of a pathogen to alter (hopefully reduce) the number of infectors over time.
First, the general idea of Montecarlo sampling.
Montecarlo Sampling to approximate integrals#
Suppose we have a random variable \(X \sim f_{X}\). In other words, \(X\) has some arbitrary probability density (or mass) function. Further, we may wish to approximate some expected value from this random variable or even the probability density function associated with \(X\), especially if \(f_{X}\) doesnt have a closed form.
To be more concrete, we may want to compute quantities like
or like
where \(I\) is called the indicator function.
If the integral cannot be computed because it does not have a readily available closed form then we must approximate it. There are many different ways to approximate an integral like the two above. We will explore Montecarlo sampling.
The Montecarlo approximation to th above \(E(g(X))\) takes the following steps: (1) sample a set of \(n\) values from the random variable \(X\) according to the probability density \(f\); (2) for the \(i^{\text{th}}\) sampled \(x\) value, called \(x^{i}\), compute \(g(x^{i})\) (3) finally, average over all \(g(x^{i})\) values. This can be written as
Lets look at an example Suppose we aee modeling the number of infcted individuals over time using the Reed-Frost model with parametrs \(N\) and \(p\). We are quite sure about the number of animals in the system, \(N\), but are uncertain about \(p\). Instead of a fixed value for \(p\) of 0.05, we decide \(p \sim \text{Beta}( 50 \cdot 0.05, 50 (1-0.05) )\).
Lets take a quick look at this distribution.
import scipy.stats
import numpy as np
import matplotlib.pyplot as plt
domain = np.linspace(0,1,50)
pdf_vals = scipy.stats.beta( 50*0.05, 50*(1-0.05) ).pdf(domain)
fig,ax = plt.subplots()
ax.plot(domain,pdf_vals)
ax.set_xlabel("p")
ax.set_ylabel(r"$f_{p}$")
Text(0, 0.5, '$f_{p}$')
Though our probability density for \(p\) is peaked around 0.05, there is still considerable uncertainty around that value. Our goal then is to understand statistical properties of the Reed-Frost model given this density around \(p\).
Given \(i0=1\), typically, \(I_{1} \sim \text{Binomial}(S_{0}, 1-(1-p)^{i0}) \) and the expected value is then \(E(I_{1}) = S_{0}(1-(1-p)^{i_{0}})\). But how will our expected value change now that \(p \sim \text{Beta}\) ? How about the distribution of \(I_{1}\)?
Well our model for \(I_{1}\) now depends on \(p\), a probability density, and so then instead of sampling \(I_{1}\) we need to first sample \(p\) and then sample \(I_{1}\).
The fact that \(p\) is now a probability density makes it trickier to see what the density of \(I_{1}\) is because we must now compute a joint probability density over \((p,I_{1})\)
and because now the density for \(I_{1}\) by itself (without \(p\)) must be marginalized
Well this wont be any fun at all! Before, we could enumerate all the \(i\) values and then plug them into the Binomial density. Now we need to handle the additional uncertainity in \(p\) in addition to \(I\).
The Montecarlo approach is to sample \((p^{i},I^{i}_{1})\) for \(i\) from one to some large number and then estimate statistical quantities of interest. Lets try this approach by sampling \(p\) and then sampling \(I_{1} | p\) to collect \(1000\) joint samples of \((p,I_{1})\).
import pandas as pd
samples = { "p":[], "I1":[] }
n=100
i0=10
for sample in range(1000):
p = np.random.beta(50*0.05,50*(1-0.05))
I1_p = np.random.binomial(n-i0, 1 - (1-p)**i0 )
samples["p"].append(p)
samples["I1"].append(I1_p)
samples = pd.DataFrame(samples)
fig,ax = plt.subplots()
ax.hist(samples["I1"],10,density=True)
domain = np.arange(min(samples["I1"]), max(samples["I1"]))
ax.plot( domain, scipy.stats.binom(n-i0, p=1-(1-0.05)**10).pmf(domain), label="p fixed" )
ax.set_xlabel("Number of infected")
ax.set_ylabel("Probability density")
ax.legend()
plt.show()
Above, we plotted the estimated probability density as a histogram in blue and the Binomial probability density if we considered \(p\) fixed at the value 0.05. We note that the density of the number of infected at time unit one \((I_{1})\) has much more uncertainity when \(p\) is also uncertain, compared to when \(p\) is fixed.
Montecarlo sampling the full time series#
We can extend our montecarlo sampling of \(I_{1}\) to \(I_{2}\) and so on by following the same procedure as above.
In this way, we can sample from this process jointly the distrbution \((p,I_{1},I_{2},\cdots, I_{T})\)
import seaborn as sns
samples = { "p":[], "t":[], "I":[], "uncert":[] }
n=1000
p=0.00125
T=25
for sample in range(4000):
p = np.random.beta(50*0.001225,50*(1-0.001225))
S = n-1
R = 0
I = 1
for t in range(T):
I_tp1 = np.random.binomial(S, 1 - (1-p)**I )
S_tp1 = S - I_tp1
R_tp1 = R + I
S = S_tp1
R = R_tp1
I = I_tp1
samples["p"].append(p)
samples["t"].append(t)
samples["I"].append(I_tp1)
samples["uncert"].append(True)
#--without uncertainity
for sample in range(4000):
p = 0.001225
S = n-1
R = 0
I = 1
for t in range(T):
I_tp1 = np.random.binomial(S, 1 - (1-p)**I )
S_tp1 = S - I_tp1
R_tp1 = R + I
S = S_tp1
R = R_tp1
I = I_tp1
samples["p"].append(p)
samples["t"].append(t)
samples["I"].append(I_tp1)
samples["uncert"].append(False)
samples = pd.DataFrame(samples)
samples = samples.rename(columns={"uncert":"Uncertainity"})
sns.lineplot(
x="t",
y="I",
data=samples,
hue="Uncertainity",
estimator="mean",
errorbar=("pi", 95)
)
<Axes: xlabel='t', ylabel='I'>
#--initial conditions
s0 = 100
i0 = 1
r0 = 0
#--probability of transmission
p=0.015
#--create three lists and add the initial number of s,i,r to each list.
s = [s0]
i = [i0]
r = [r0]
def pr_inf(p,i):
return 1 - (1-p)**i
start = 0
end = 32
step = 1
for t in np.arange(start,end,step): #<--Whats this function doing?
#--simulate the number of infectors at time t+delta
num_of_sus_at_time_t = s[-1] ## <- why?
num_of_inf_at_time_t = i[-1] ## <- why?
i__t_plusDelta = np.random.binomial( n=num_of_sus_at_time_t
,p=pr_inf(p,num_of_inf_at_time_t) )
#--compute the number of susceptibles at time t+delta
s__t_plusDelta = num_of_sus_at_time_t-i__t_plusDelta
#--compute the number of removed at time t+delta
r__t_plusDelta = r[-1] + num_of_inf_at_time_t
#--add the results to our lists.
s.append(s__t_plusDelta)
i.append(i__t_plusDelta)
r.append(r__t_plusDelta)
4. Exploring the impact of an intervention using Monte-Carlo#
Mathemtical setup#
Monte-carlo procedures can also be useful for exploring how an intervention, or some modification to the Reed-Frost process, impacts important outcomes like: the total number of infections, the peak number of infections, where the peak number of infections occur, etc etc.
Lets first look at a mathematical setup for an intervention that occurs at time \(t*\) and, when this intervention occurs, the innate probbaility of infections is reduced from \(p\) to \(\rho p\) where \(1-\rho\) describes how effective our intervention is. If \(1-\rho\) is 1 then \(p \to 0\) and our intervention is so effective that after \(t*\) no one can be infected anymore. If, on the other hand, that \(1-\rho\) is 0 then \(p \to p\) and our intervention has no impact.
The Reed-Frost model can be modified as follows
The only change we have now is that \(p\) is a step function. This is one way to model an intervention that is planned for \(t*\) and has a constant effect after. Modified models like this can be useful when discussing plan for when to implement a well-known intervention. In most cases, you would work with decisionmakers to decide what outputs from this model are most important.
A common output that is important is the total number of infections.
An intervention with fixed efficacy \(\rho\) can have a different impact depending on when it is applied to the population.
The goal, then, is to understand the distribution of the total number of infections given a fixed efficacy and differing intervention times \(t*\).
An analytic solution to the number of infections at the end of the outbreak (say at time \(T\)) is time consuming to compute exactly. In fact, it may not even be feasible to compute this exactly. It may take too long. Instead, we can use a Monte-carlo procedure to estimate the total number of infections.
We know that the total number of infections, \(\tau\) , is a deterministic function of the number of infections at every time step, or
The variability in \(\tau\) is not related to how it is computed, but the varability in the individual number of infections per time unit. That is, the total number of infections is a transformation of the folllowing random variables
and we can (a) pick a value \(t*\) (b) sample from \(I_{1},I_{2},I_{3},\cdots,I_{T}\) (c) compute the total number of infections per sample (d) use a histogram to estimate the probability density for \(\tau | t*\). The dataframe we build should look like this
How to encode that mathematics into the computer#
To encode our step function for \(p\) into the Reed-Frost model, we can use an if/else block. As the mode loops through time units \(t=0,1,2,3,4\) and so on, we will check that \(t \ge t*\). When this condition is true we will set \(p \to p \times \rho\) and otherwise set \(p \to p \).
Lets specify the model to have parameters
Covariate |
Value |
|---|---|
\(p\) |
0.01 |
\(\rho\) |
0.50 |
\(S_0\) |
200 |
\(I_0\) |
1 |
\(R_0\) |
0 |
\(t^*\) |
10 |
#--initial conditions
s0 = 200
i0 = 1
r0 = 0
#--probability of transmission
p = 0.01
reduction_by_intervention = 0.50 #<--new line here
time_of_intervention = 10 #<--new line here
p_intervention = p*reduction_by_intervention #<--new line here
#--create three lists and add the initial number of s,i,r to each list.
s = [s0]
i = [i0]
r = [r0]
#--function for the probability of infection
def pr_inf(p,i):
return 1 - (1-p)**i
#--observational period
start = 0
end = 32
step = 1
for t in np.arange(start,end,step): #<--Whats this function doing?
#--is t greater than t*? Is so then reduce p
if t>=time_of_intervention:
ptrans = p_intervention
else:
ptrans = p
#--simulate the number of infectors at time t+delta
num_of_sus_at_time_t = s[-1] ## <- why?
num_of_inf_at_time_t = i[-1] ## <- why?
i__t_plusDelta = np.random.binomial( n=num_of_sus_at_time_t
,p=pr_inf(ptrans,num_of_inf_at_time_t) )
#--compute the number of susceptibles at time t+delta
s__t_plusDelta = num_of_sus_at_time_t-i__t_plusDelta
#--compute the number of removed at time t+delta
r__t_plusDelta = r[-1] + num_of_inf_at_time_t
#--add the results to our lists.
s.append(s__t_plusDelta)
i.append(i__t_plusDelta)
r.append(r__t_plusDelta)
plt.plot(i)
plt.axvline(10,color="black")
<matplotlib.lines.Line2D at 0x7ca5bf212810>
Above, we plotted (for a single instance) the number of infections over time in blue. The black vertical curve is the time of the intervention (t*) that reduced the probability of infections from \(p\) to \(p \rho\). From this single instance, we see that intervening at time unit 10 would be a bit late. However, this only tells us about a single sample from our model. We can learn more about intervention timing by understanding the distribution of (say) total number of infections over time given a specific intervention time.
dataset = { "simulation":[], "time":[], "infections":[], "t*":[] }
#--initial conditions
s0 = 200
i0 = 1
r0 = 0
rho = 0.50
#--probability of transmission
p=0.01
def pr_inf(p,i):
return 1 - (1-p)**i
start = 0
end = 32
step = 1
SIMS = 500 #<- Number of simulations to run.
for sim in range(SIMS): #<-What does range do?
for tstar in [1,3,5,7,9,11]:
#--create three lists and add the initial number of s,i,r to each list.
s = [s0]
i = [i0]
r = [r0]
for t in np.arange(start,end,step): #<--Whats this function doing?
#--simulate the number of infectors at time t+delta
num_of_sus_at_time_t = s[-1] ## <- why?
num_of_inf_at_time_t = i[-1] ## <- why?
#--THIS is the implementation of that step function
if t>= tstar:
pinf = pr_inf(p,num_of_inf_at_time_t) * rho
else:
pinf = pr_inf(p,num_of_inf_at_time_t)
i__t_plusDelta = np.random.binomial( n=num_of_sus_at_time_t ,p=pinf )
#--compute the number of susceptibles at time t+delta
s__t_plusDelta = num_of_sus_at_time_t-i__t_plusDelta
#--compute the number of removed at time t+delta
r__t_plusDelta = r[-1] + num_of_inf_at_time_t
#--add the results to our lists.
s.append(s__t_plusDelta)
i.append(i__t_plusDelta)
r.append(r__t_plusDelta)
#--data collection
dataset["simulation"].append(sim)
dataset["time"].append(t)
dataset["infections"].append(i__t_plusDelta)
dataset["t*"].append(tstar)
dataset = pd.DataFrame(dataset)
dataset #<--default is to print
| simulation | time | infections | t* | |
|---|---|---|---|---|
| 0 | 0 | 0 | 0 | 1 |
| 1 | 0 | 1 | 0 | 1 |
| 2 | 0 | 2 | 0 | 1 |
| 3 | 0 | 3 | 0 | 1 |
| 4 | 0 | 4 | 0 | 1 |
| ... | ... | ... | ... | ... |
| 95995 | 499 | 27 | 0 | 11 |
| 95996 | 499 | 28 | 0 | 11 |
| 95997 | 499 | 29 | 0 | 11 |
| 95998 | 499 | 30 | 0 | 11 |
| 95999 | 499 | 31 | 0 | 11 |
96000 rows × 4 columns
Lets take a quick look at this dataset. We can see that we have individual columns for the number simulation that we run where a simulation was the number of infeceted individuals over time, and the intervention time \(t*\).
Lets first make a summary dataset of the total number of infections per simulation and per t*.
def add_up_infections(x):
return pd.Series({"ttl":np.sum(x["infections"])})
summary_data = dataset.groupby(["t*","simulation"]).apply(add_up_infections).reset_index()
summary_data
| t* | simulation | ttl | |
|---|---|---|---|
| 0 | 1 | 0 | 0 |
| 1 | 1 | 1 | 67 |
| 2 | 1 | 2 | 22 |
| 3 | 1 | 3 | 16 |
| 4 | 1 | 4 | 5 |
| ... | ... | ... | ... |
| 2995 | 11 | 495 | 141 |
| 2996 | 11 | 496 | 156 |
| 2997 | 11 | 497 | 169 |
| 2998 | 11 | 498 | 1 |
| 2999 | 11 | 499 | 162 |
3000 rows × 3 columns
sns.boxplot(x="t*",y="ttl",data=summary_data)
<Axes: xlabel='t*', ylabel='ttl'>
Lets simulate our original Reed-Frost model (without an intervention) 2000 times and store the final epidemic sizes in the list RF_original_epidemic_size.
Next, we’ll simulate our intervention model and store the final epidemic sizes in the list RF_intervention_epidemic_size.
Then we can compare the performance of our intervention.
We will add to our original model, a list called RF_original_epidemic_size and append to this loop the final epidemic size.
Final epidemic size for Original Reed Frost model#
Now we have all the final epidemic sizes stored in RF_original_epidemic_size.
Lets look at the (appeoximate) probability distribution of the final epidemic size.
plt.hist(RF_original_epidemic_size,20,density=True)
plt.ylabel("Proportion")
plt.xlabel("Final epidemic size (RF-model original)")
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 1
----> 1 plt.hist(RF_original_epidemic_size,20,density=True)
2 plt.ylabel("Proportion")
3 plt.xlabel("Final epidemic size (RF-model original)")
4 plt.show()
NameError: name 'RF_original_epidemic_size' is not defined
#--initial conditions
s0 = 100
i0 = 1
r0 = 0
#--probability of transmission
p=0.015
def pr_inf(p,i):
return 1 - (1-p)**i
#--probability of transmission
p = 0.015
reduction_by_intervention = 0.50 #<--new line here
time_of_intervention = 10 #<--new line here
p_intervention = p*reduction_by_intervention #<--new line here
start = 0
end = 32
step = 1
SIMS = 5000 #<- Number of simulations to run.
RF_intervention_epidemic_size = []
for sim in range(SIMS): #<-What does range do?
#--create three lists and add the initial number of s,i,r to each list.
s = [s0]
i = [i0]
r = [r0]
for t in np.arange(start,end,step): #<--Whats this function doing?
#--is t greater than t*? Is so then reduce p
if t>=time_of_intervention:
ptrans = p_intervention
else:
ptrans = p
#--simulate the number of infectors at time t+delta
num_of_sus_at_time_t = s[-1] ## <- why?
num_of_inf_at_time_t = i[-1] ## <- why?
i__t_plusDelta = np.random.binomial( n=num_of_sus_at_time_t
,p=pr_inf(ptrans,num_of_inf_at_time_t) )
#--compute the number of susceptibles at time t+delta
s__t_plusDelta = num_of_sus_at_time_t-i__t_plusDelta
#--compute the number of removed at time t+delta
r__t_plusDelta = r[-1] + num_of_inf_at_time_t
#--add the results to our lists.
s.append(s__t_plusDelta)
i.append(i__t_plusDelta)
r.append(r__t_plusDelta)
#--Model is done running for the observational period.
#--Now we can record the final epidemic size for this run
RF_intervention_epidemic_size.append( r[-1] ) #<-umm, why is it r[-1]?
plt.hist(RF_intervention_epidemic_size,20,density=True)
plt.ylabel("Proportion")
plt.xlabel("Final epidemic size (RF-model intervention)")
plt.show()
Lets compare the two plots. To do this, we will ask to plot a histogram for both the final size of the original model and the final size of the intervention model over the same bins.
bins = np.linspace(0,101,20)
fig,ax = plt.subplots()
ax.hist(RF_original_epidemic_size,bins ,color="red" , density=True, alpha=0.50,label="Orig.")
ax.hist(RF_intervention_epidemic_size,bins ,color="blue", density=True, alpha=0.50,label="Inter.")
plt.ylabel("Proportion")
plt.xlabel("Final epidemic size (RF-model intervention)")
plt.legend()
plt.show()
Homework#
Explore unknown \(N\)
Explore unknown \(R_{0}\)
Explore the interplay of timing and efficacy