Epidemic modeling - Part 7
Convolution, deconvolution, and daily data
- Motivation for write-up
- Available data
- Convolution
- Finding $h_L[j]$ and $h_I[j]$
- Convolution in practice
- De-convolution
- Deconvolution in practice
- Discussion
#collapse_hide
!pip install plotly==4.9.0
import pandas as pd
import numpy as np
import math
from scipy import signal
import plotly.graph_objects as go
import plotly.express as px
from scipy.stats import expon
from scipy.stats import gamma
from scipy.stats import weibull_min
from numpy.random import default_rng
rng = default_rng()
# Let's build a numerical solution
def seir_model(init, parms, days):
S_0, E_0, I_0, R_0 = init
Epd, Ipd, Rpd = [0], [0], [0]
S, E, I, R = [S_0], [E_0], [I_0], [R_0]
dt=0.1
t = np.linspace(0,days,int(days/dt))
sigma, beta, gam = parms
for _ in t[1:]:
next_S = S[-1] - beta*S[-1]*I[-1]*dt
Epd.append(beta*S[-1]*I[-1]*dt)
next_E = E[-1] + (beta*S[-1]*I[-1] - sigma*E[-1])*dt
Ipd.append(sigma*E[-1]*dt)
next_I = I[-1] + (sigma*E[-1] - gam*I[-1])*dt
Rpd.append(gam*I[-1]*dt)
next_R = R[-1] + (gam*I[-1])*dt
S.append(next_S)
E.append(next_E)
I.append(next_I)
R.append(next_R)
return np.stack([S, E, I, R, Epd, Ipd, Rpd]).T
#collapse_hide
# Need this new function for model below:
def make_df(p,num_E, num_I, num_R):
df = pd.DataFrame(np.full((p,1), 'S').T[0], columns=['State'])
df['Day'] = 0
tochange=df.loc[rng.choice(p, size=num_E+num_I+num_R, replace=False),'State'].index
df.loc[tochange[0:num_E],'State'] = 'E'
df.loc[tochange[num_E:num_I+num_E],'State'] = 'I'
df.loc[tochange[num_E+num_I:num_E+num_I+num_R],'State'] = 'R'
return df
#collapse_hide
def seir_model_stoch(beta, p, num_E, num_I, num_R, days):
# Initialize population dataframe with data given by user
df = make_df(p,num_E, num_I, num_R)
# This variable is used to track daily value of beta if it varies over time
xxbeta=np.array([],dtype=float)
# Initialize the arrays to return
# Below are numbers of S, E, I, R total
S=np.array([],dtype=int)
E=np.array([],dtype=int)
I=np.array([],dtype=int)
R=np.array([],dtype=int)
# Below are the daily additions in S, E, I, R
Spd=np.array([],dtype=int)
Epd=np.array([],dtype=int)
Ipd=np.array([],dtype=int)
Rpd=np.array([],dtype=int)
b=beta
# Stochastic model so use random values to decide on progression
rand = np.random.random(size=(p,days))
# Depending if you want exponential or gamma distribution for T_Latent
EtoI = gamma.rvs(1.8,loc=0.9,scale=(5.2-1.8)/0.9,size=p)
# Depending if you want exponential, gamma, or Weibull distribution for T_Infectious
# Uses distributions found on blog part 3
ItoR = weibull_min.rvs(2.3, loc=2, scale=20.11, size=p)
# Iterate over every day the simulation is run
for j in range(0,days-1):
# Record daily beta values
xxbeta=np.append(beta, b)
# First we get the index of the individuals that will change state today:
# Random number tells you which 'S' have been exposed on this day
StoE_index = df.loc[(df.State == 'S') & (rand[:,j] < b[j]*len(np.where(df.State=='I')[0])/p)].index
# For each row, if a person has been a certain number of days in E, they will go to I
# This follows EtoI variable which is either exponential or gamma distributed according to above
EtoI_index = df.loc[(df.State == 'E') & (j-df.Day >= EtoI)].index
# Similaraly as above
# For each row, if a person has been a certain number of days in I, they will go to R
# This follows EtoI variable which is either exponential or gamma distributed according to above
ItoR_index = df.loc[(df.State == 'I') & (j-df.Day >= ItoR)].index
# Use indexes collected above to populate per day values
Epd = np.append(Epd,len(StoE_index))
Ipd = np.append(Ipd,len(EtoI_index))
Rpd = np.append(Rpd,len(ItoR_index))
# Now we use the indexes collected above randomly to change the actual population dataframe to the new states
df.loc[ItoR_index, 'State'] = 'R'
df.loc[EtoI_index, 'State'] = 'I'
df.loc[StoE_index, 'State'] = 'E'
df.loc[ItoR_index, 'Day'] = j
df.loc[EtoI_index, 'Day'] = j
df.loc[StoE_index, 'Day'] = j
df.loc[ItoR_index, 'DayR'] = j
df.loc[EtoI_index, 'DayI'] = j
df.loc[StoE_index, 'DayE'] = j
# Append the S, E, I, and R arrays
S=np.append(S,len(np.where(df.State=='S')[0]))
E=np.append(E,len(np.where(df.State=='E')[0]))
I=np.append(I,len(np.where(df.State=='I')[0]))
R=np.append(R,len(np.where(df.State=='R')[0]))
# Code below for control measures to reduce beta values
# if ((I[-1] > 1000) & (Ipd[-1] > 399)):
# b = beta2
# elif ((I[-1] > 1000) & (Ipd[-1] < 400)):
# b = beta3
Epd[0]+=num_E
Ipd[0]+=num_I
Rpd[0]+=num_R
return S,E,I,R, Epd, Ipd, Rpd, xxbeta, df
This is the 7th part of a multi-part series blog post on modeling in epidemiology.
The goal of this 7th installment is to expand on the notions seen in part 5.
Most notably, the idea is to expand on the notion of convolution and deconvolution and to see how it can be useful to describe an epidemic.
If you look at the the COVID-19 trackers around the web, or even mine, you can get a sense of what data is available for study.
Generally speaking we have the following:
- Total number of positive individuals
- Daily number of newly diagnosed individuals
- Daily number of recovered individuals
- Daily number of deaths
Research has also shown some vague numbers on both $T_{Latent}$ and $T_{Infectious}$.
About the data:
- We don't have data for exposure - how can we get it using deconvolution ??
- The data for some regions can be very sparse
- $h_L[j]$ is simply the probability of an individual having a latent period of j days
- $h_I[j]$ is similarly the probability of an individual having an infectious period of j days
#collapse_show
days = np.arange(100)
cdf = pd.DataFrame({
'T_Latent': gamma.cdf(days, 1.8,loc=0.9,scale=(5.2-1.8)/0.9),
'T_Infectious': weibull_min.cdf(days, 2.3,loc=2,scale=20.11)
})
h_L = cdf.diff().T_Latent
h_I = cdf.diff().T_Infectious
h_L[0] = 0
h_I[0] = 0
First run the SEIR model to obtain the actual $I_{pd}$:
#collapse_hide
# Define parameters for stochastic model
days = 300
p = 10000
num_E = 1
num_I = 0
num_R = 0
beta_stoch = 0.5*np.ones(days)
# Run 2 stochastic simulations
results_stoch0 = seir_model_stoch(beta_stoch, p, num_E, num_I, num_R, days)
Now obtain $h_L[j]\circledast E_{pd}[j]$:
#collapse_hide
Ipd=signal.fftconvolve(h_L, results_stoch0[4], mode='full')
#collapse_hide
fig = go.Figure(data=[
go.Scatter(name='Ipd_Actual', x=np.arange(len(results_stoch0[5])), y=results_stoch0[5]),
go.Scatter(name='Ipd_convolved', x=np.arange(len(Ipd)), y=Ipd)
])
fig.update_layout(
xaxis_title = 'Day',
yaxis_title = 'Count',
title={
'text':r'$\text{Actual } I_{pd} \text{ vs. } h_L[j]\circledast E_{pd}[j]$',
'x':0.5,
'xanchor':'center'
}
)
fig.show()
We can use the actual $R_{pd}$ from the SEIR model above.
Now obtain $h_I[j]\circledast I_{pd}[j]$:
#collapse_hide
Rpd=signal.fftconvolve(h_I, results_stoch0[5], mode='full')
#collapse_hide
fig = go.Figure(data=[
go.Scatter(name='Rpd_Actual', x=np.arange(len(results_stoch0[6])), y=results_stoch0[6]),
go.Scatter(name='Rpd_convolved', x=np.arange(len(Rpd)), y=Rpd)
])
fig.update_layout(
xaxis_title = 'Day',
yaxis_title = 'Count',
title={
'text':r'$\text{Actual } R_{pd} \text{ vs. } h_I[j]\circledast I_{pd}[j]$',
'x':0.5,
'xanchor':'center'
}
)
fig.show()
Problems with de-convolution
Illl-posed
Unfortunately de-convolution is not as straightforward as the convolution.
A quick recap of the math, in the frequency domain we have (Fourier or z transform) $$F \{ f ∗ g \} = F \{ f \} ~ F \{ g \}$$ And so if we have: $$y=\ f ∗ g$$ Then: $$F \{ y \} = F \{ f \} ~ F \{ g \}$$ And: $$g = F^{-1} \left(\frac{F\{y\}}{F\{f\}} \right) $$
But this last equation is often ill-posed.
Noise in data
We can see above that daily new data varies a lot from day to day due to the stochasticity of the model. This results in high-frequency decomposition in the frequency domain.
Initial conditions and boundary-value problems
Boundary values refer to the values on day 0 and the last day we have data for. Due to the way data is reported on JHU, we someimes have crayz values here, most commonly a value of 0 for the latest day if the data isn0t reported on time for example. These gaps can create issues for deconvolution.
In the next blog post I go around this by simplz averaging the last 3 days instead of using only the last day.
Iterative deconvoluton
Multiple iterative deconvolution algorithms exist: Lucy_Richardson, Gold, or Van Cittert among others.
LR has already been used to get daily incidence of the Spanish flu based on death records in Philadelphia at the time for example 1.
An adaptation of those above is used here with the basics below (see code for details):
Let's use the equation for $I_{pd}$ as an example: $$I_{pd}[j] = \sum_{n_L=0}^{M_L-1}h_L[n_L]~E_{pd}[j-n_L]= h_L[j]\circledast E_{pd}[j]$$
The idea is that after an initial guess for $E_{pd}$ we can iteraively find a better one.
With the $n^{th}$ guess for $E_{pd}$ written as $E_{pd, n}$ we have: $$I_{pd}[j]= h_L[j]\circledast E_{pd, n}[j]$$ $$\leftrightarrow 0 = I_{pd}[j] - h_L[j]\circledast E_{pd, n}[j]$$ $$\leftrightarrow E_{pd,n+1} = E_{pd,n} + I_{pd}[j] - h_L[j]\circledast E_{pd, n}[j]$$
Where $O[j] = I_{pd}[j] - h_L[j]\circledast E_{pd, n}[j]$ is the error term.
Hence we want to minimze O[j].
Initial guess
We know our $h_L$ and $h_I$ are in part simply delay functions, so an initial first guess is to simply use the same signal delayed by the time difference caused by the impulse response.
To go around those high-frequency notes in the data, we can siply use a lowpass filter implemented as a butterworth filter below:
The basic effect is that it smoothes the daily values as we have when convolving from the previous state daily values.
#collapse_show
def lowpass(x, fc=0.05):
fs = 1 # Sampling frequency
t = np.arange(len(x)) #select number of days done in SEIR model
signala = x
#fc = 0.05 # Cut-off frequency of the filter
w = fc / (fs / 2) # Normalize the frequency
b, a = signal.butter(5, w, 'low')
return signal.filtfilt(b, a, signala)
#collapse_show
# Let's define an iteration function:
def iter_deconv(alpha, impulse_response, input_signal, delay, comparator):
conv=signal.fftconvolve(impulse_response, input_signal, mode='full')
correction=np.roll(comparator-conv[:len(comparator)], delay)
input_signal=np.floor(lowpass(alpha*correction+input_signal))
input_signal[input_signal<0]=0
return input_signal
# Define a function to return MSE between two signals as a measure of goodness of fit
def msecalc(A, B):
return ((A - B)**2).mean(axis=0)
#collapse_show
#regularization parameter
alpha=2
# Setup up the resultant Ipd we want to compare our guess with
Ipd=np.floor(lowpass(results_stoch0[5]))
Ipd[Ipd<0]=0
# Find delay caused by h_L
delay=Ipd.argmax()-signal.fftconvolve(Ipd, h_L, mode='full').argmax()
# We want initial guess to simply be the result of the convolution delayed
initial_guess=np.roll(Ipd,delay)
Enext = initial_guess
# AN array to record MSE between result we want and our iterated guess
mse=np.array([])
mse=np.append(mse, 100000)
mse=np.append(mse, msecalc(Ipd, signal.fftconvolve(h_L, Enext, mode='full')[:len(Ipd)]))
itercount=0
while mse[-1] < mse[-2]:
itercount=itercount+1
Enext=iter_deconv(alpha, h_L, Enext, delay, Ipd)
mse=np.append(mse, msecalc(Ipd, signal.fftconvolve(h_L, Enext, mode='full')[:len(Ipd)]))
print("Iteration #" + str(itercount) +": MSE= "+str(mse[itercount]))
print("Iteration #" + str(itercount+1) +": MSE= "+str(mse[-1])+" so we use the result of the previous iteration.")
#collapse_hide
fig = go.Figure(data=[
go.Scatter(name='E_pd', x=np.arange(150), y=results_stoch0[4]),
go.Scatter(name='Epd=lowpass(E_pd)', x=np.arange(150), y=lowpass(results_stoch0[4])),
go.Scatter(name='Epd=deconv(I_pd)', x=np.arange(150), y=Enext),
go.Scatter(name='I_pd', x=np.arange(150), y=results_stoch0[5]),
go.Scatter(name='Ipd=lowpass(I_pd)', x=np.arange(150), y=lowpass(results_stoch0[5])),
go.Scatter(name='Ipd=conv(E_pd)', x=np.arange(150), y=signal.fftconvolve(h_L, lowpass(results_stoch0[4]), mode='full')),
go.Scatter(name='Ipd=conv(deconv(Ipd))', x=np.arange(150), y=signal.fftconvolve(h_L, Enext, mode='full')[:len(Ipd)])
])
fig.update_layout(
xaxis_title = 'Day',
yaxis_title = 'Count',
title={
'text':r'$\text{Actual } E_{pd} \text{ vs. deconvolution of } I_{pd}$',
'x':0.5,
'xanchor':'center'
}
)
fig.show()
The iterative deconvolution seems to work nicely and we get close results.
#collapse_hide
#regularization parameter
alpha=2
# Setup up the resultant Ipd we want to compare our guess with
Rpd=np.floor(lowpass(results_stoch0[6]))
Rpd[Rpd<0]=0
# Find delay caused by h_I
delay=Rpd.argmax()-signal.fftconvolve(Rpd, h_I, mode='full').argmax()
# We want initial guess to simply be the result of the convolution delayed
initial_guess=np.roll(Rpd,delay)
Inext = initial_guess
# AN array to record MSE between result we want and our iterated guess
mse=np.array([])
mse=np.append(mse, 100000)
mse=np.append(mse, msecalc(Rpd, signal.fftconvolve(h_I, Inext, mode='full')[:len(Rpd)]))
itercount=0
while mse[-1] < mse[-2]:
itercount=itercount+1
Inext=iter_deconv(alpha, h_I, Inext, delay, Rpd)
mse=np.append(mse, msecalc(Rpd, signal.fftconvolve(h_I, Inext, mode='full')[:len(Rpd)]))
print("Iteration #" + str(itercount) +": MSE= "+str(mse[itercount]))
print("Iteration #" + str(itercount+1) +": MSE= "+str(mse[-1])+" so we use the result of the previous iteration.")
#collapse_hide
fig = go.Figure(data=[
go.Scatter(name='Ipd', x=np.arange(150), y=results_stoch0[5]),
go.Scatter(name='Ipd=lowpass(Ipd)', x=np.arange(150), y=lowpass(results_stoch0[5])),
go.Scatter(name='Ipd=deconv(Rpd)', x=np.arange(150), y=Inext),
go.Scatter(name='Rpd', x=np.arange(150), y=results_stoch0[6]),
go.Scatter(name='Rpd=conv(Inext)', x=np.arange(150), y=signal.fftconvolve(h_I, Inext, mode='full')),
go.Scatter(name='Rpd=lowpass(Rpd)', x=np.arange(150), y=lowpass(results_stoch0[6])),
go.Scatter(name='Rpd=conv(deconv(Rpd))', x=np.arange(150), y=signal.fftconvolve(h_I, Inext, mode='full')[:len(Rpd)])
])
fig.update_layout(
xaxis_title = 'Day',
yaxis_title = 'Count',
title={
'text':r'$\text{Actual } I_{pd} \text{ vs. deconvolution of } R_{pd}$',
'x':0.5,
'xanchor':'center'
}
)
fig.show()
Again our algorithm works nicely.
There are a lot of papers and research done in the field of deconvolution in general but I have found only one using it in epidemic models (link above).
The theory is skipped here but it is required if we want to use more advanced techniques. I may write another blog post with more details.
After doing this you might be left wondering, well what's the point ?
A few things:
- If we have daily data of infectious $I_{pd}$ we can then use the deconvolution technique to find an approximation for $E_{pd}$
- If we have a close estimate for $E_{pd}$ we can in turn estimate the value of $\beta$ or $R$
- We can verify the robusteness of $h_L$ and $h_I$ against real world data and how these vary geographically in a pandemic
These points will be studied a bit more in the next blog post.