Introduction

COVID-19 (coronavirus disease 2019) is a disease caused by a virus named SARS-CoV-2. It is very contagious and has quickly spread around the world. The virus constantly changes through mutation. Since the beginning of the pandemic, we've seen several prominent variants, including Alpha, Beta, Delta, and Omicron. Among them, Omicron is now the dominant strain of the coronavirus in the world that causes COVID-19. The Omicron variant spreads more easily than earlier variants of the virus that cause COVID-19, including the Delta variant. CDC (Centers for Disease Control and Prevention) expects that anyone with Omicron infection, regardless of vaccination status or whether or not they have symptoms, can spread the virus to others [1]. Learning more about the Omicron variant and preventing its spread helps to fight COVID-19.

In this project, we analyze the confirmed cases of Omicron in South Africa from a time series analysis perspective. In particular, we hope to describe the spread of Omicron given this time interval by using tools in the time analysis aspect. To achieve that, we fit the data using the ARMA model and POMP models and compare their performance.

Exploratory Data Analysis

The data is obtained from the “COVID-19 Data Hub” [2]. We study the Omicron cases in South Africa because it is the place where the Omicron variant of the coronavirus SARS-CoV-2 was first detected and then spreads around the world a little more than two months after [3]. We focus on the data after December \(1^{st}\), 2021 when the Omicron became the dominant version of the confirmed cases in South Africa. The original data for the number of confirmed cases is the cumulative data. In this analysis, we make the difference between the confirmed cases and those in the previous days to achieve the number of daily new confirmed cases. We notice from the time plot that the data skews to the right with the peak appearing in the mid of December. We have 156 data points with a mean of 5160.167 and a standard error of 6593.398. The highest number of daily confirmed cases is 37875 which happened on December \(11^{th}\), 2021. The majority of the data points are between the range of 1278 and 5074. The data is visualized in the following, where the red line represents the 7-days moving average, and the blue line represents the daily new death cases.

We noticed that there are some unusual fluctuations with a period of about 7 days, indicating that there might be some periodic pattern. The autocorrelation plot is the following:

ARMA model

The stationary ARMA(p,q) model with parameters \((\phi_{1:p}, \psi_{1:q}, \mu, \sigma^2)\) is given by \([4]\):

\[\Phi(B)(Y_n - \mu) = \Psi(B)\epsilon_n\]

and functions and parameters are defined as these:

\[\mu = E[Y_n]\]

\[\Phi(x) = 1 - \phi_1 x - \phi_2 x^2 - ... - \phi_p x^p\]

\[\Psi(x) = 1 + \psi_1 x + \psi_2 x^2 + ... + \psi_q x^q\]

\[\epsilon_n \sim N(0, \sigma^2)\]

The optimal \((p,q)\) parameters are decided using AIC criterion and for this purpose, we generate AIC table as the following:

MA0 MA1 MA2 MA3 MA4
AR0 3189.38 3107.23 3070.22 3049.54 3037.27
AR1 3020.99 3001.61 3002.37 3003.95 3005.82
AR2 3008.00 3002.22 3001.48 3005.95 3003.91
AR3 3007.26 3004.18 3002.78 3000.01 3001.73
AR4 3007.66 3005.37 3004.49 3001.71 2991.87

This suggests a model of \(ARMA(3,3)\), and the model fitting summary are the following:

## 
## Call:
## arima(x = data$confirmed.new, order = c(3, 0, 3))
## 
## Coefficients:
##          ar1     ar2     ar3     ma1     ma2      ma3  intercept
##       0.0152  0.0583  0.8211  0.4895  0.3882  -0.5748   4173.565
## s.e.  0.0631  0.0654  0.0492  0.0882  0.0923   0.0843   2969.656
## 
## sigma^2 estimated as 11466978:  log likelihood = -1492,  aic = 3000.01
##          aic loglikelihood
## [1,] 3000.01     -1492.005

We then do a causality and invertibility check. Based on the result, the model is invertible and causal.

## Absolute value of AR polynomial roots are:  1.039101 1.082603 1.082603
## 
## Absolute value of MA polynomial roots are:  1.000002 1.000002 1.739727

We do a small model diagnostic as the following, which seems to suggest that our model performance is not that well. Hence, we continue to conduct the specturm analysis to find some periodic patterns.

Spectrum Analysis

We can find the period by conducting a spectrum analysis. We approximate the spectral density function of the data by using a smoothed periodogram \([4]\). There seem to be a couple of dominant frequencies as is shown in the spectrum plot. The first \(\omega_1\) corresponds to the number of days the data set which does not seem to be interesting. The next \(\omega_2\) of 0.142 corresponds to a period of about 7. This suggests that the periodicity of our data is about one week. There is another local maximum \(\omega_3\) at around 0.283 and this corresponds to the second harmonic of \(\omega_2\) so that we can ignore it. One-week periodicity fits in with our intuition as many institutions post the number of confirmed cases weekly. As a result, the data we have may show the periodicity of a week. The period may also be because there is less staff working on Covid testing for people on weekends so that fewer confirmed cases are noticed.

POMP - SIR model

In this section we tried to analyse the omicron data using the basic POMP SIR model. We treat the omicron pandemic as a partially observed Markov process, with only three basic states:

  1. Susceptible - All individuals at the start

  2. Infectious - Individuals that has been infected by Omicron

  3. Recovered - Individuals that has been recovered

Here is the state diagram of the model:

Here is the parameters defined for our model

  1. \(N = S(t) + I(t) + R(t)\) is the total number of the population. This variable is given as a constant and it is not estimated

  2. \(\mu_{SI} = \dfrac{\beta I(t)}{N}\) is the infection rate. We are able to change \(\beta\) which is defined as the force of infection. The variable \(B\) is estimated.

  3. \(\mu_{IR} = \gamma\) is the recovery rate. This variable is estimated.

  4. \(\rho\) is the rate in which cases are reported. This variable is estimated.

  5. \(\eta\) is the another parameter we are going to estimate. It is defined as the fraction of population that are susceptible. This variable is estimated.

  6. We also have one additional variable \(H\) which is an accumulator variable. It is the number of individuals who have moved from State \(I\) to the state \(R\) over the course of that week. Accumulator variables are not estimated.

We model the data by a negative binomial variable \([4]\): \(reports \sim NegBin(\rho H(t), k)\). The variable k is also estimated.

We simulated our model with different parameters to find a suitable initialization for our pomp model. After several simulations, the model is initialized with values of \(N = 5*10^7\) (fixed), \(\beta = 9\) (will be estimated), \(\mu_{IR} = 0.04\) (will be estimated), \(\rho = 0.5\) (will be estimated), \(k=10\) (will be estimated), \(\eta = 0.04\) (will be estimated). It is worth to mention that reaching these values also involved particle filtering runs. Here it the simulations for initial values:

After running the POMP model (mif2 function), we get the following results. The log likelihood soon converges.

It seems that all parameters converge to specific intervals while widths of intervals change. Here is the simulation of the model with estimated group of parameters \(\beta=16.49857\), \(\mu_{IR}=0.02254904\), \(\rho=0.4073849\), \(k=3.991062\), \(\eta=0.04078927\), \(N=50000000\),the best log likelihood is -1997, with a standard error of 6.98.

Then, we conduct a global search over the parameter space \(\beta \sim [0.1, 50], \mu_{IR}\sim[0.01, 1], \rho\sim[0.01, 0.99], k\sim[0.1, 10], \eta \sim [0.01,0.99]\). The likelihood for the local and global result is the following, with a simulation of the best model estimated with parameters \(\beta=0.932\), \(\mu_{IR}=0.0216\), \(\rho=0.964\), \(k=0.402\), \(\eta=0.575\), \(N=50000000\),the best log likelihood is -1677, with a standard error of 0.243.

Since the simulation result is not that good, we are also going to utilize more complex POMP model.

POMP - SEAPIRD model

Since the SIR model fails to give a fairly good performance, we implemented a more complicated framework comparing to the basic SIR model, the SEAPIRD (Susceptible-Exposed-Asymptomatic-Presymptomatic-Infectious-Recovered-Deceased) model \([6]\). Comparing to the basic SIR and SEIR model, it adds three mode states: asymptomatic, presymptomatic and deceased. Since some omicron patients do not show any clear symptoms while remaining infectious to other people, it may be proper for us to include that state in our model, and the presymptomatic state acts like a bridge state from exposed to infectious (with symptoms). Also, even though omicron variant is not as fatal as the delta variant, death cases are still reported in major countries, and hence we also include the deceased states to denote death. Graphically, our model can be represented in this figure:

Though the basic logic is the same with the SIR model, we still introduced several new parameters here. \(\beta\) is the transmission rate from Susceptible people to exposed people, with the number of change modeled by \(\Delta N_{SE} \sim Binomial(S, 1-e^{-\beta \frac{I+A+P}{N}\Delta t})\). \(\alpha\) denotes the proportion of asymptomatic patients, therefore we set the asymptomatic change follows \(\Delta N_{EA} \sim \alpha \times Binomial(E, 1-e^{-\mu_{EI}\Delta t})\). Consequently, the presymptomatic change follows \(\Delta N_{EP} \sim (1-\alpha) \times Binomial(E, 1-e^{-\mu_{EI}\Delta t})\). The rest parts are similar with basic SIR model \([5]\):

\[ \begin{aligned} \Delta N_{AR} &\sim Binomial(A, 1-e^{-\mu_{AR}\Delta t}) \\ \Delta N_{PI} &\sim Binomial(P, 1-e^{-\mu_{PI}\Delta t}) \\ \Delta N_{IR} &\sim Binomial(I, 1-e^{-\mu_{IR}\Delta t}) \\ \Delta N_{ID} &\sim Binomial(I, 1-e^{-\mu_{ID}\Delta t}) \\ \end{aligned} \]

Also, since local government will possibly take procedures to keep virus from spreading, hence, we add three parameters \(c_1, c_2, c_3\) into the transition from susceptible from exposed, each measures 50 days of the Markov process. \(\Delta N_{SE} \sim Binomial(S, 1-e^{-\beta \times c_i \frac{I+A+P}{N}\Delta t})\).

Another different thing from the SIR model is that we model the weekly reported cases as normal distribution rather than negative binomial. Since our reported cases is fairly large and the mode is around 1925. Thus, we parametrize the reported cases as \(Y_{cases} \sim Normal(\rho H, \tau \rho H (1-\rho))\).

The final step is parameter configuration and initialization. we set our \(\Delta_t\) as 1 day, and \(S=500000, I=169\) according to our data. We tried to first simulate using a parameter setting: \(\beta=1.83,\mu_{IR}=0.02,\mu_{ID}=0.0002,\mu_{EI}=0.86,\alpha=0.67,\mu_AR=0.85,\mu_{PI}=0.81,c_1=1.06,c_2=4.79,c_3=1.45,\rho=0.98,N=500000,\tau=154735\). Another note is that in this part, we used the 7-day average smoothed version of the daily new cases.

Simulation result is poor using this parameter setting. But the shape somehow looks similar, and therefore we conduct a local search and a global search. For the local search, we can see that the log likelihood soon converges, and some of the parameter like \(\mu_{EI}, \rho\) also converge quickly. The other parameters seem not converge to a point, but they remain stable in a range. Also, to show how the resulting parameter performs, we put the parameter with the best log likelihood result.

We can see that both of them fits well, and generally, according to the pairwise scatter plot, generally the global search gives a better log likelihood result, but the results looks pretty similar. For the local search, the best parameter setting is \(\beta=1.15,\mu_{IR}=0.697,\mu_{ID}=0.000325,\mu_{EI}=0.0372,\alpha=0.710,\mu_AR=0.180,\mu_{PI}=0.232,c_1=4.56,c_2=1.59,c_3=0.830,\rho=0.984,N=500000,\tau=262000\), with a log likelihood of -1448.337 and standard error 0.006812. For the global search, the best parameter setting is \(\beta=2.56,\mu_{IR}=0.0385,\mu_{ID}=0.00000864,\mu_{EI}=0.266,\alpha=0.0285,\mu_AR=3.49,\mu_{PI}=0.178,c_1=0.565,c_2=10.6,c_3=4.43,\rho=0.993,N=500000,\tau=624000\), with a log likelihood of -1446.122 and standard error 0.006424.

The log likelihood convergence plot is the following. We can see that the log likelihood does not converge to a point due to the small number of iterations, but it becomes steady near the end of \(Nmif\) iterations, which is an improvement comparing to the SIR model.

Conclusion

The purpose of this analysis is to analyze the confirmed cases of Omicron in South Africa from a time series analysis perspective. We tried modelling this time series using ARMA, SIR and SEAPIRD models and compared their performances. Our final conclusion is that SIR model is not that competitive comparing to the basic simple ARMA model, judging from the likelihood perspective. However, SEAPIRD model did a great job but it somehow reaches the bottleneck since both local and global search gives approximately the same result. We may investigate other models further to see if it gives better performance.

References

[1] Centers for Disease Control and Prevention. Omicron Variant: What You Need to Know, Apr 16, 2022. URL https://www.cdc.gov/coronavirus/2019-ncov/variants/omicron-variant.html?s cid=11734:omnicron%20variant:sem.ga:p:RG:GM:gen:PTN:FY22.

[2] Ardia D. Guidotti, E. COVID-19 Data Hub, Apr 16, 2022. URL https://covid19datahub.io/articles/r.html

[3] Smriti Mallapaty. Where did Omicron come from? Three key theories, Apr 16, 2022. URL https://www.nature.com/articles/d41586-022-00215-2.11

[4] Analysis of Time Series lecture notes 2022. URL https://ionides.github.io/531w22/

[5] Final Project 13 in 2021: An Investigation into COVID-19 in California. URL https://ionides.github.io/531w21/final_project/project13/blinded.html

[6] Mahmoud Harmouch, Disease Spread Modeling With The “SEIRD” Model. URL https://blog.devgenius.io/covid-19-modeling-using-the-sierd-model-and-visualization-using-plotly-and-ipywidgets-e6d5fbfc07aa