Least Squares Estimation for the SIRD model

Question

I'm experiencing some difficulties in the estimation of the parameters $\alpha, \beta, \gamma$ for the following discrete-time SIRD (Susceptibles, Infected, Recovered, Dead) model with sampling step of 1 day

$$\tag{1}\begin{cases} S_{t}&=S_{t-1}-\alpha\frac{S_{t-1}I_{t-1}}{N} \\ I_{t}&=I_{t-1}+\alpha\frac{S_{t-1}I_{t-1}}{N}-\beta I_{t-1}-\gamma I_{t-1} \\ R_{t}&=R_{t-1}+\beta I_{t-1} \\ D_{t}&=D_{t-1}+\gamma I_{t-1} \\ \end{cases} \qquad \text{for} \,\, t=1,2,\dots$$

that I've found in this paper. In order to find the unknown $\alpha, \beta, \gamma$, I want to use the least squares regression in his closed-form solution. The parameter $N$ is the size of the population under study, so it is known and hasn't to be estimated.

1 Derivation of the LS estimator

1.1 definitions

1. Let's consider the dataset $D_T\triangleq\{y_0,\dots,y_T\}$ up to the observation horizon $T$, where $y_t\triangleq[S_t, I_t, R_t, D_t]'$ is the vector of the observed values at time $t$ for the variables $S,I,R,D$. Here $'$ denotes the transpose operation, thus $y_t$ is a column vector in $\mathbb{N}^{4\times1}$;

2. Let $\theta\triangleq[\alpha, \beta, \gamma]'$ be the generic vector of parameters. The prediction model $\hat{y}_t(\theta)$ is $(1)$, so $$\tag{2} \hat{y}_t(\theta)\triangleq \begin{bmatrix} S_{t-1}-\alpha\frac{S_{t-1}I_{t-1}}{N} \\ I_{t-1}+\alpha\frac{S_{t-1}I_{t-1}}{N}-\beta I_{t-1}-\gamma I_{t-1} \\ R_{t-1}+\beta I_{t-1} \\ D_{t-1}+\gamma I_{t-1} \end{bmatrix} \qquad \text{for} \,\, t=1,2,\dots $$ with the convention that $\hat{y}_0(\theta)\triangleq 0$;

3. Let $V_T(\theta)\triangleq \frac{1}{2}\sum _{t=0}^T \|y_t-\hat{y}_t(\theta) \|^2$ the quadratic cost up to $T$. Here $\| \cdot \|$ denotes the euclidian norm. The least square estimator $\theta_\text{LS}$ of the 'real' parameter $\bar{\theta}$ is defined as $$\tag{3}\theta_\text{LS}\triangleq \arg\min_{\theta \in \mathbb{R^3}} V_T (\theta)$$ i.e. the minimum for the cost $V_T$.

1.2 analitic solution of $(3)$

the idea to solve $(3)$ is to use the standard technique by solving with respect to $\theta$ the equation $$\tag{4}\frac{\partial V_T(\theta)}{\partial \theta}=0$$ the solution is a minimum for $V_T$ since $(3)$ is a convex problem under mild assumptions regarding the dataset $D_T$ (invertibility of the next matrix $R_T$ defined below). In order to solve $(4)$, let's start by observing that the prediction model $(2)$ is linear in their parameters. In fact we can write that $$\tag{5}\hat{y}_t(\theta)=\varphi_t \theta + y_{t-1} \qquad \text{for} \,\, t=0, 1, 2,\dots$$ by introducing the regression matrices in $\mathbb{R^{4\times3}}$ $$\tag{6}\varphi_t \triangleq \begin{bmatrix} -\frac{S_{t-1}I_{t-1}}{N} & 0 & 0 \\ \phantom{-}\frac{S_{t-1}I_{t-1}}{N} & -I_{t-1} & -I_{t-1} \\ 0 & \phantom{-}I_{t-1} & 0\\ 0 & 0 & \phantom{-}I_{t-1} \end{bmatrix} \qquad \text{for} \,\, t=1,2,\dots$$ with the conventions that $\varphi_0, y_{-1}=0$. From $(5)$ it follows straightforward that the gradient of the cost $V_T$ is $$\tag{7}\begin{align}\frac{\partial V_T(\theta)}{\partial \theta} &= \sum_{t=0}^T - \frac{\partial \hat{y}_t (\theta)}{\partial \theta}[y_t-\hat{y}_t(\theta)]\\ &=-\sum_{t=0}^T \varphi_t'[y_t-(\varphi_t \theta + y_{t-1})] \\ &=\sum_{t=1}^T \varphi_t'[\varphi_t \theta - \Delta y_t] \\ &=\left(\sum_{t=1}^T \varphi_t '\varphi_t\right)\theta - \sum_{t=1}^T \varphi_t'\Delta y_t \end{align}$$ where $\Delta y_t \triangleq y_t-y_{t-1}$. If we introduce the matrix $R_T\in\mathbb{R}^{3\times3}$ and the vector $\tilde{\theta}_T\in\mathbb{R}^{3}$ $$\tag{8}R_T\triangleq \sum_{t=1}^T \varphi_t '\varphi_t \qquad \tilde{\theta}_T\triangleq \sum_{t=1}^T \varphi_t'\Delta y_t$$ the gradient in $(7)$ gets the following final sintetic expression $$\tag{9}\frac{\partial V_T(\theta)}{\partial \theta} = R_T\theta-\tilde{\theta}_T$$ now, by combining $(4)$ with $(9)$ and by resolving with respect $\theta$, we can finally conclude that the least square estimator that we are searching is $$\tag{10}\boxed{\theta_\text{LS}=R_T^{-1}\tilde{\theta}_T}$$

2 Naive implementation in Python

2.1 dataset

I want to estimate $\bar{\theta}$ for the COVID-19 epidemy in Italy, so I've built the dataset by retrieving from worldometers.info the number of infected $I_t$, recovered $R_t$ and dead $D_t$ individuals day by day. Since $S_t+I_t+R_t+D_t=N$ is costant in time, the number of susceptibles day by day is $S_t=N-(I_t+R_t+D_t)$.

2.2 least squares estimation of the parameters

in order to compute $(10)$, we need:

1. to build $\varphi_t$ and $\Delta y_t$. For the former we can use the definition $(6)$, for the latter we can observe that $$\tag{11} \Delta y_t \triangleq y_t-y_{t-1}=\begin{bmatrix} S_t-S_{t-1} \\ I_t-I_{t-1} \\ R_t-R{t-1} \\ D_t-D_{t-1} \end{bmatrix} \qquad \text{for} \,\, t=1,2,\dots$$
2. to build $R_T$ and $\tilde{\theta}_T$. The idea to compute $(8)$ is to accumulate during the time the products $\varphi_t ' \varphi_t$ and $\varphi_t '\Delta y_t$.

After this 2 simple step the estimation is given by $(10)$.

2.3 simulation

for the simulation we use the prediction model $(1)$ with the least squares parameters that we have just found. For the initial condition of the simulation I consider the situation where in the population there is only one infected individual that spreads the disease to the other people.

$$\begin{cases} S_{0}&=N-1 \\ I_{0}&=1 \\ R_{0}&=0 \\ D_{0}&=N-(S_0+I_0+R_0) \\ \end{cases}$$ here the starting number of dead individuals is obtained by imposing the costraint $S_0+I_0+R_0+D_0=N$.

2.4 code

  import matplotlib.pyplot as plt
  import numpy as np

  #1 DATASET

  #observed infected 
  oI = np.array([    3,     3,     3,     3,     3,     4,    19,
                    75,   152,   221,   310,   455,   593,   822,
                  1049,  1577,  1835,  2263,  2706,  3296,  3916,
                  5061,  6387,  7985,  8514, 10590, 12839, 14955,
                 17750, 20603, 23073, 26062, 28710, 33190, 37860,
                 42681, 46638, 50418, 54030, 57521, 62013, 66414 ])
  #observed recovered
  oR = np.array([    0,     0,     0,     0,     0,     0,     1,
                     2,     2,     2,     3,     4,    46,    47,
                    51,    84,   150,   161,   277,   415,   524,
                   590,   623,   725,  1005,  1046,  1259,  1440,
                  1967,  2336,  2750,  2942,  4026,  4441,  5130,
                  6073,  7025,  7433,  8327,  9363, 10362, 10951 ]) 
  #observed dead
  oD = np.array([   0,     0,      0,     0,     0,     0,     1,
                    2,     3,      7,    11,    12,     7,    21,
                   29,    41,     52,    79,   107,   148,   197,
                  233,   366,    463,   631,   827,  1016,  1266,
                 1441,  1809,   2158,  2503,  2978,  3405,  4032,
                 4825,  5476,   6077,  6820,  7503,  8215,  9134 ])                         
  #observed susceptibles 
  N = 60*1000000  #population size
  T = oI.size    #observation horizon

  oS = np.zeros((T,))

  for t in range(0, T):
      oS[t] = N-(oI[t]+oR[t]+oD[t])    

  ##############################################################################

  #2 LEAST SQUARES ESTIMATION OF THE PARAMETER

  #initializazion of RT and thetatildeT
  RT = np.zeros((3,3))  
  thetatildeT = np.zeros((3,))  

  #construction of RT and thetatildeT
  for t in range(1, T):
      #definition of phit and Deltayt
      phit = np.array([  [-oS[t-1]*oI[t-1]/N,          0,         0],                        
                         [ oS[t-1]*oI[t-1]/N,   -oI[t-1],  -oI[t-1]], 
                         [                 0,    oI[t-1],         0],
                         [                 0,          0,   oI[t-1]]  ])
 
      Deltayt = np.array([oS[t]-oS[t-1], oI[t]-oI[t-1], 
                          oR[t]-oR[t-1], oD[t]-oD[t-1] ])
 
      #accumulation in RT and thetatildeT
      RT += np.dot(phit.transpose(),phit)
      thetatildeT += np.dot(phit.transpose(), Deltayt)
 
  #least squares estimation
  thetaLS = np.dot(np.linalg.inv(RT), thetatildeT)
 
  ##############################################################################

  #3 PREDICTION

  #prediction model parameters
  alpha = thetaLS[0]
  beta = thetaLS[1]
  gamma = thetaLS[2]

  #initialization of the prediction model variables
  S = np.zeros((T,))
  I = np.zeros((T,))
  R = np.zeros((T,))
  D = np.zeros((T,))

  #initial condition of the prediction
  S[0] = N-1
  I[0] = 1
  R[0] = 0
  D[0] = N-(S[0]+I[0]+R[0])
  
  #simulation
  for t in range(1,T):
      S[t] = S[t-1]-alpha*(S[t-1]*I[t-1]/N)
      I[t] = I[t-1]+alpha*(S[t-1]*I[t-1]/N)-beta*I[t-1]-gamma*I[t-1]
      R[t] = R[t-1]+beta*I[t-1]
      D[t] = D[t-1]+gamma*I[t-1]

  #############################################################################

  #4 PLOTS

  fig, axs = plt.subplots(2, 1, constrained_layout=True)
  axs[0].set_title('Observed Data')
  axs[0].plot(range(0,T), oI)
  axs[0].plot(range(0,T), oR)
  axs[0].plot(range(0,T), oD)
  axs[0].legend("IRD 1",loc="upper left")
  axs[1].set_title('Predicted Data')
  axs[1].plot(range(0,T), I)
  axs[1].plot(range(0,T), R)
  axs[1].plot(range(0,T), D)
  axs[1].legend("IRD 1",loc="upper left")

2.5 results

the prediction model doesn't work well, this is the plot of the prediction errors between the observed data and the predicted data.

I can't understand if somewhere I have made some mistake or if the estimation that I'm using cannot provide good predictions.

Answer 1

1

Your equation number 5 should be

$$\hat{y}_t(\theta)=\varphi_t \theta + \hat{y}_{t-1} \qquad \text{for} \,\, t=0, 1, 2,\dots$$

instead of

$$\hat{y}_t(\theta)=\varphi_t \theta + y_{t-1} \qquad \text{for} \,\, t=0, 1, 2,\dots$$

2

Also you compute the derivative $\varphi$ based on the matrix that contains values of $S, I ,R, D$ that are the observed values, but $\varphi$ should relate to the modeled values.

---

3

I am not sure whether you can continue your attempt to solve this analytically after correcting those mistakes. It looks a bit like how people solve equations with the finite element method and potentially you could solve it in that way as well (but it will be an approximation in terms of polynomial functions, and is not exact).

Another way to solve it is to put the equations as a function, and have some solver optimize it (you can have the solver estimate the gradient). You can read about that here: Fitting SIR model with 2019-nCoV data doesn't conververge

In addition: You can recast the equations into a single differential equation. For the SIR model this is demonstrated here:

Tiberiu Harko, Francisco S. N. Lobo, M. K. Mak Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates arXiv:1403.2160 [q-bio.PE]

The SIRD model is almost analogous. It is nearly the same model, with only the R split up in two. So you can use this differential equation to make initial estimates of parameters.

---

Also

Fitting this kind of data to some model might be a bad idea. The SIR type of models are some sort of logistic growth type models where the growth starts approximately exponentially but eventually the growth rate decreases. It is due to such terms like $dI/dt = I * (factor)$ where the factor is decreasing as $I$ (and $R$ and $D$) grow (in the case of logistic growth the factor is $1-I$, for the SIRD model it a bit more sophisticated but not much different).

However, in the case of the corona epidemic, you get a decrease in the growth rate for a multitude of reasons.

• Weather changes ($R_0$ is not a constant)

• Spatial distribution (This virus spreads in from place to place, and should not be considered with compartmental models that assume homogeneous mixing; a person in Milan is much more likely to infect their family, neighbors, co-workers than a random person in the rest of Lombardy)

• Stochastic Time effects. The article that you refer to tries to bring autocorrelation into the mixture, but you also have some stochastic behavior, people are not gonna get sick exactly at the same time. Some people will get sick earlier than others and this will be according to some function that increases in time and that will make an increase of cases or deaths that might appear as an exponential growth that relates to a transmission model, but it might be not.

• Sampling bias. We can also see rapid increase in the sampling due to biassed sampling. Definitions of the disease are changing (this gave a rapid bump in the curve for Chinese), tests might be limited (several countries are limiting their testing which might give false ideas of reduction in the growth of the cases), positive reinforcement (once people have discovered the disease suddenly many other cases might become assigned to the same cause, and this may potentially occur inaccurate because a single cause of death is not always possible to assign)

  The last to points Sampling Bias and Stochastic time behavior might have occured in the outbreak of SARS (2003) in Amoy Gardens where hundreds of people got sick in a very short time frame. Instead of fitting a model to it, one could also assume that all these hundreds of cases were infected by a single person (and that might be a more likely scenario). Possibly such a situation may have occurred in Italy as well, an initial heavy seeding by unnoticed cases that is now spreading with some time effect and causes the initial exponential decrease (currently the growth looks more like a quadratic curve).

• Last but not least, people respond to the virus which may cause it's spread to increase/decrease. Currently heavy measures have been taken and this restricts to a large extend the ability of the virus to spread. You can not model this with a model that has parameters that are constant in time (well, you can, but the outcome will be meaningless)

The logistic SIR type models will interpret all those reasons for a reduction in the growth rate as a reproduction rate very close to 1 or a low population parameter (you fixed it at the size of the population under study, but this is arbitrary, and also not everybody is gonna be susceptible, possibly many people might have some sort of immunity and get little infected, e.g. a Hoskin's effect or some other effect might make only/mostly the elderly population susceptible).

This makes the seemingly mechanistic model, meaningless regarding the parameters. The outcome will be unrealistic.

Answer 2

Regarding Attempt 0.1: I ran your Python code, and it produces parameters that are about half of what the paper produced, and as you said, the simulation produces vastly smaller cases (about 20,000 times less than observed values). However, even when I used the parameters from the paper to replace the values from thetaLS (at the start of the code section“#3 Prediction”) , I still got vastly smaller number of cases. Do you know what that is caused by? I reread your simulation code repeatedly, and it looks right to me.