Statistical modeling of COVID-19 deaths with excess zero counts

Introduction

Non-Gaussian time series modeling is a subject that has attracted many researchers recently and a number of studies have emerged in the literature to propose alternative solutions for modeling of the realistic situations like count time series. Two widely known probability distributions have been employed in modeling sequences of observed counts. First, we have the Poisson distribution which is characterized by the property of equal dispersion, or the equality between its mean, and its variance. A random variable would exhibit over-dispersion if its variance is greater than its mean. Secondly, a common choice of probability distribution for unbounded counts with over-dispersion is the negative binomial distribution which is used instead of the binomial distribution to reflect observed over dispersion in the data. Since the traditional ARMA time series models require Gaussian errors which may not be met by count random variables, a class of generalized autoregressive moving average models was developed to handle flexible observation-driven models for non-Gaussian time series models. In these models, the dependent variable can have a conditional exponential family distribution, such as the Poisson or the negative binomial distribution, given past history of the process. Following the original generalized autoregressive moving average (GARMA) models developed by Benjamin, Rigby, and Stasinopoulos (2003), a number of subsequent studies presented special cases and variations of the original modeling, where a link function is introduced and the conditional mean under that link function takes an ARMA structure (Zheng, Hiao, and Chen 2015).

However in some cases, a specific count, usually zero, can occur more often than other series counts and if we overlook the frequent occurrence of zero values, this could result in over-dispersion and biased estimators of the parameter and the standard deviation. Therefore we may obtain misleading statistical inferences. To account for frequent zero values in the data, Lambert (1992) presented zero-inflated models where the distribution is a mixture of some count models, such as the Poisson or the negative binomial, with one that is degenerate at zero. But authors caution that in order to use properly these models, one has to verify whether these zeros are structural or random (Tang et al. 2018).

In order to proceed further with appropriate statistical modeling of zero inflation with time series of counts, a number of research contributions have emerged recently in the literature to address the issue. Some scholars including Zhu and co-authors (2012) introduced the zero-inflated Poisson (ZIP) and the zero-inflated negative binomial (ZINB) models in a context of integer-valued generalized autoregressive conditional heteroskedasticity framework. Others (Yang, Cavanaugh, and Zamba 2015) developed dynamic models to deal with zero-inflated count time series with both ZIP and ZINB predictive distributions. In this paper, we study the empirical performance of these models and assess their forecast accuracy when the data contains frequent zeros with a particular pattern using COVID-19 data from Tunisia.

Materials and methods

The objective of this research is to use alternative count series models to study how they fit COVID-19 death records and to study their prediction accuracy. Several mathematical and statistical models have been introduced recently to study the spread of the coronavirus disease and to forecast its related deaths and infections. For instance, Mohd Yusoff (2020) used a system dynamic methodology to build a COVID-19 confirmed case model. However, as argued in Alahmadi et al. (2020), while these computational models and the use of algorithms and computational resources are useful to predict the spread of infectious diseases and allow disease dynamic models to be parameterized with good accuracy, the lack of identification of model parameters can limit notably the usefulness of these models. In addition, publicly available data on daily infections and deaths from developing countries may not always be reliable. With these challenges in mind, we suggest studying the forecast performance of coronavirus related deaths with competing statistical models based on parametric assumptions.

We consider two classes of models. The first class includes generalized count data models with an ARMA structure to depict observations dependence (ACC). The second class includes zero-inflated models for a time series of counts (ZIM). Although ZIM models seem to be the obvious choice when we have frequent zeros in the data, a number of researchers stated that these models are being routinely applied leading to overuse (Warton 2005) and suggest that a negative binomial distribution for the counts may be sufficient to deal with over-dispersion. The validity of their argument of course is related to whether or not the observed excess zeros come from a different data generating process, and only if so, there will be a gain in using ZIM models.

In heath studies, and specifically in the presence of a deadly infectious disease, one major issue is how to control the pandemic and how to limit hospitalizations and the deaths related to it. The zero death records may suggest an effective health policy or a proper intervention and the opposite is true if the policy or the health protocol is not effective. So with such data, one may need to know if there is dependence between the intervention policy and the observed zero counts and thus whether these zeros are in fact structural.

Results

We apply these two methods to a published COVID-19 death data for the North African country of Tunisia available from (https://www.statista.com/statistics/1185550/number-of-coronavirus-related-deaths-tunisia/). The interesting part about these data is the specific pattern of the zero counts. Most death counts were zero prior to July 31st, 2020, before the government decided to intervene and to change its policy by relaxing the strict travel restrictions in order to boost the economy. Then the post-intervention sample data shows a significant increase in daily deaths related to the pandemic.

Tunisia had an early success in dealing with the COVID-19 pandemic. In fact the Tunisian government had an effective policy to control and to limit the spread of the disease when the pandemic first hit the country in the early spring of 2020. Right after the first cases were recorded in Tunisia in March 2020, the government responded swiftly with a comprehensive set of measures aiming to slow the spread of the virus. The government also suspended all travel, mandated working from home for non-essential workers, imposed mandatory confinement and nightly curfews, and shut down schools and businesses and banned public gatherings. In addition, the military and the police forces were asked to enforce these strict measures to prevent and control the spread of the virus. However the economic situation in the country was sluggish and researchers suggested that the pandemic is likely to exacerbate Tunisia’s development challenges so more people will probably fall below the poverty line (Kokas et al. 2020) and the economic situation in the country could go from bad to worse if the government does not loosen the restrictions. This might explain why the Tunisian government in late July 2020 decided to lift travel restrictions and lockdown. However, after opening its borders in order to boost the tourism sector and to fix its stalled economy, Tunisia saw a significant spike in the number of its COVID-19 infection cases and more alarming was the number of deaths resulting from the progression of the pandemic. The latest statistics are displayed in Figure 1 and show the updated cumulative numbers of COVID-19 deaths in Tunisia, from March 15, 2020 to September 10, 2021.

Figure 1:

Number of COVID-19 related deaths in Tunisia (March 15, 2020–September 10, 2021. Source: Statista).

The official statistics on deaths related to the pandemic from March 15, 2020 to February 11, 2021, which represents the period covered by this study, indicate that the Tunisian data shows 148 zero values out of a total of 334 sample daily observations. This represents 44.3% of the total sample. Table 1 lists some descriptive statistics for the death counts before and after the end of July 2020, a date which marks government intervention to lift travel restrictions.

Figure 2 displays the empirical probability distribution of deaths and clearly shows a significant proportion of zero death count (44.3%) occurring before the travel restriction relaxing policy which was implemented by the Tunisian government in late July 2020.

Figure 2:

The probability distribution of COVID-19 deaths in Tunisia (March 15, 2020 to February 11, 2021).

Before presenting the statistical findings of the count time series models and their performance in forecasting COVID-19 deaths, we test the data for possible unit roots using the Augmented Dickey Fuller test (Said and Dickey 1984) and the Phillips–Perron test (Phillips and Perron 1988). The results are displayed in Table 2 below. The choice of lag lengths in the test regressions is based on AIC and BIC information criteria.

Both ADF and PP tests indicate that the death count series is stationary at the 10% level of significance.

Given that death counts sharply increased after July, 31st 2020 following a government policy to relax travel restrictions, we include an intervention covariate (X _t ) in the conditional mean model to assess its significance for the observed surge of deaths and we compare it with similar count models that do not include intervention as an additional covariate.

Discussion

We use R codes to perform the statistical analysis for this study. First, we start with a model with intervention (relaxation of travel restrictions). Based on the government policy to opt for travel restrictions removal in late July 2020, and referring to Liboschik et al. (2016) who studied the covariate effects in time series count data, we include in our models a covariate to account for the intervention policy.

We estimate alternative models with different lag lengths and based on AIC and BIC information criteria we select the model that yields the best fit and which is given by:

where λ ̂ t is the estimated conditional mean at time t, Y _t is the death count at time t, and the intervention variable X _t is modeled as one level shift that takes 1 after July 31st, 2020 and the value 0 prior to that. In Table 3, we report the estimation results of model Eq. (7) with Poisson and negative binomial conditional distributions. The models are estimated with the maximum likelihood method (Fokianos 2012) and the results show statistical significance of the coefficients in the models with travel intervention. The coefficients beta_2 and beta_3 correspond to the regression on the second lagged and the third lagged observations, respectively. The coefficient eta_1 refers to the single intervention policy variable (X). We can conclude that all these coefficients are statistically significant at the 5% level because their confidence intervals do not contain zero.

The quasi information criterion (QIC) has been proposed by Pan (2001) as an alternative to Akaike’s information criterion (AIC) which is properly adjusted for regression analysis based on the generalized estimating equations. As noted by the author, QIC computes the quasi information criteria of a generalised linear model for time series of counts. In case of models with the Poisson distribution the QIC has approximately the same value as the AIC. However, in case of models with another distribution it can be a more adequate alternative to the AIC. It should also be reminded that the lower value of an information criterion (AIC, BIC, or QIC) indicates a better-fit model.

Next, we run the same models without intervention and the results are displayed in Table 3. It is shown that the coefficients in the models with intervention are more significant than those without intervention. In particular, for the negative binomial distribution the reported confidence intervals include the value zero and thus suggest that the coefficients are not statistically significant at the 5% level.

Also, for the nested models, we test the significance of the intervention policy and we obtain the following result:

Chisq-statistic: 312.8429 on 1 degree(s) of freedom, p-value < 0.0001. This proves that the relaxation of travel restrictions had a significant effect to increase COVID-19 deaths in Tunisia.

Furthermore, and in order to evaluate the adequacy of these ACC models where frequent zeros were not dealt with in a specific way, we analyse the statistical performance of the models. First, we refer to the Probability Integral Transform (PIT) which relates to the result that data values that are modeled as being random variables from any given continuous distribution can be converted to random variables having a standard uniform distribution on the interval [0, 1]. Thus, if the model is adequate (the predictive distribution is correct) then PIT should be close to the uniform distribution as discussed in Czado, Gneiting, and Held (2009) and also in Gneiting, Balabdaoui, and Raftery (2007). Figure 3 shows the residuals’ correlogram obtained from the ACC model with the relaxation of travel restrictions.

Figure 3:

Autocorrelation function of the residuals from the autoregressive conditional count model (ACC) with intervention.

Figures 4 and 5 display the probability integral transform of the Poisson distribution and the negative binomial distribution, respectively. It is shown that the histograms of the probability integral transform (PIT) have a non-uniform shape for both distributions, suggesting over dispersion.

Figure 4:

Probability integral transform of the autoregressive conditional count model with intervention, for the Poisson distribution.

Figure 5:

Probability integral transform of the autoregressive conditional count model with intervention, for the negative binomial distribution.

The statistical results show that both predictive distributions are not well calibrated and they do not fit closely to the uniform distribution, likely because of the excess zero counts.

Another additional criterion to assess their fit is the marginal calibration, which is the difference between the average predictive cumulative density function and the empirical distribution function of the observations. If the predictions from a given model are appropriate then marginal calibration should be close to zero as explained in Christou and Fokianos (2015). It can be noticed from Figures 6 and 7 that there is a significant difference between predictive and observed cumulative density functions for both distributions. This result in fact, proves our previous finding and again points to a poor fit of ACC models for the COVID-19 death counts.

Figure 6:

Marginal calibration for the Poisson probability. The plot displays a significant difference between predictive and empirical distributions specifically for counts close to zero, suggesting a poor fit of the ACC model.

Figure 7:

Marginal calibration for the negative binomial probability. The plot displays a significant difference between predictive and empirical distributions in particular for death counts close to zero, which suggests a poor fit of the ACC model.

In addition to the criteria discussed above regarding the fit of these models, one can assess both calibration and sharpness by computing a single score from alternative scoring rules as described by Gneiting and Katzfuss (2014). In this paper, we choose only four scoring rules; from each we compute a mean score in order to evaluate the probabilistic forecast of the Poisson and the negative binomial distributions. The notion of sharpness refers to the concentration of the predictive distribution and can be measured by the prediction interval width.

The statistical analysis in this study is based on the tscount package in R (Liboschik et al. 2017). We estimate the models and evaluate the model adequacy for each of the predictive distributions by computing mean difference scores. The model with the lowest score is preferable. These scores are summarized in Table 4 and show in fact inconclusive results regarding which distribution has an overall better fit for COVID-19 deaths.

In the next step, we estimate the zero-inflated Poisson mixed auto regression model (AR-ZIP) and the zero-inflated negative binomial mixed auto regression model (AR-ZINB) with intervention for the COVID-19 death data. Specifically, we use a first order autoregression to represent the dependence between the responses. The results are reported in Table 5 and show a good fit for both count specifications given by each distribution. The AR(1) coefficient and the intervention policy are both significant in the ZIP model.

The statistical results displayed in Table 5 show the estimation of models with and without the travel ban lift policy. With respect to AR-ZIP models with and without intervention, we did a simple computation of BIC from the reported log likelihood functions and the number of parameters in each model to confirm that the model with intervention has a better fit. Also, we find that the relaxation of travel restrictions in Tunisia has a significant impact and the forecasts for COVID-19 deaths based on models with intervention are more accurate than those without intervention, as indicted by the root mean square error values which are reported in Table 6.

In order to compare the empirical performance of the zero-inflation model with intervention and the autoregressive conditional probability model with intervention, we compute root mean square forecast errors for 40 days prediction of deaths for each model. The results in Table 6 show that the AR(1)-ZIP model with intervention has the best forecast performance.

To make a final assessment regarding the superiority of zero-inflation models, we additionally compute the Vuong test (Vuong 1989) for the two non-nested models (Model 1 and Model 3) and we report the test result in Table 7.

The test proves that the two competing models are distinguishable and that Model 1 outperforms Model 3 in fitting COVID-19 death counts. This presents an additional proof that zero-inflation models are better than autoregressive conditional count models in fitting and forecasting deaths related to infectious diseases in the presence of excess zero counts.

Conclusions

In this paper, we conduct an empirical analysis in order to compare the statistical performance of two types of models for count time series. Both models use non-Gaussian conditional distributions to account for possible over-dispersion in the data. The models however differ in how they treat excess zero counts. We apply these models to evaluate their fit and forecast performances using Tunisian data on COVID-19 daily deaths from March 2020 to February 2021. It is found that the zero-inflated models have more significant coefficients and present lower root mean square forecast errors. The statistical results displayed in Table 6 cast an overall outperformance of the zero inflation models. In particular, the root mean square forecast error for the AR(1) zero-inflation model with travel ban policy is nearly 16 percent less than the autoregressive conditional Poisson model with intervention, and it is 34 percent smaller than the ACP model without travel intervention. Therefore, in order to model the spread of an infectious disease and to forecast its related deaths accurately we recommend the application of zero-inflation models when excess zero counts are present in the data. Our findings also prove the statistical significance of an intervention policy implemented by the government to lift travel restrictions and which resulted in a change in the course of the pandemic death count trajectory.