Analyzing the effect of duration on the daily new cases of COVID-19 infections and deaths using bivariate Poisson regression: a marginal conditional approach

1.   Introduction

The world has been devastated by the outbreak of COVID-19 caused by SARS-CoV-2 since January 2020. While many countries could control the spread, however, observing the second wave of recurrence many other countries are witnessing the infection at an alarming rate. The developed, developing, and under-developed countries are facing the unforeseen challenges caused by the COVID-19 pandemic, which took a significant toll on various aspects of life on people all across the world. The growing COVID-19 crisis hit developing countries, not only as a health crisis in the short term but will also exhibit a devastating socio-economic impact over the months and years to come. The level of anxiety in Africa, Europe, Bangladesh, Brazil, India, Iran, the USA, South Korea, and many other countries increased due to new cases and fatalities recently. The socio-economic, psychological, and other impacts have already started ^[1,2]. It is difficult to diagnose this disease, and we already observed a delay between the onset of symptoms and an accurate diagnosis ^[3]. Also, there are still many undiagnosed and delayed-diagnosis infections due to a lack of diagnostic kits. Both the undiagnosed and diagnosed infections have a very high ability to transmit the virus to other susceptible people or family members ^[4,5]. The ability to diagnose and identify the infected persons in time and the treatment have a tremendous impact on daily new cases and deaths. The timely and effective isolation of symptomatic and confirmed COVID-19 cases is crucial for controlling this pandemic. The first wave of this pandemic already proved how it could overwhelm the under-resourced hospitals and fragile health systems in many countries.

The COVID-19 data are now available from various online sources. The majority of them are reporting daily new cases and deaths, along with very few other characteristics. For example, through official communications under the International Health Regulations (IHR, 2005), by monitoring the official ministries of health websites and social media accounts, the World Health Organization (WHO) collected the numbers of confirmed COVID-19 cases and deaths from 31 December 2019 to 21 March 2020 ^[6]. Since 22 March 2020, daily global data are compiled through WHO region-specific dashboards. Also, Worldometer ^[7] reports country-specific data along with a few other characteristics. Some sources are also reporting the number of daily tests to detect the infections. These data are mostly used to display daily trends for new cases and deaths–a significant number of research papers have already been published on the COVID-19 pandemic in various journals around the world. The objectives of much of this research are to predict the trends for new cases and deaths and assess the impact of available covariates on new cases and outcomes. Relevant policy makers of countries are using this data to make plans and to reduce the rate of this pandemic. Already, many countries have slowed down the spread of SARS-CoV-2, the virus that causes COVID-19, by taking stringent measures. However, many developing and under-developing countries failed to slow down the infection rate by adopting stringent measures due to local circumstances and culture. Every region in the world needs to make progress against the COVID-19; only then we can be safe. However, globally we are very far from the goal, and the global number of confirmed cases is growing enormously fast.

The daily new cases and deaths, two count outcomes, are naturally correlated and dependent. Analysis of these two dependent count outcomes and assessing the impact of related covariates needs a proper and accurate modeling approach. A marginal-conditional modeling approach for correlated count outcomes will allow us to determine the covariate impact on the responses jointly and better prediction. We employ a generalized linear covariate dependent bivariate Poison model to analyze the dependent daily new cases and death and assessed the impact of some characteristics on the outcomes. The rest of the paper is organized as follows. In Section 2, we present the methodology used in this paper. The results of data analysis are presented in Section 3. Finally, in Section 4 we discuss the conclusions.

2.   Bivariate Poisson regression model

In this section, we discuss the methodology to model the daily new cases of COVID-19 infections and deaths reported by world health organization (WHO). The world health organization regularly updates daily new cases and deaths for each country. It is noteworthy to mention that these daily counts of new cases and deaths are dependent or correlated. Therefore, more new cases can overwhelm the heath care system which in turn can cause more deaths. Moreover, it is vital to study the impact of duration and other available characteristics on these two dependent count outcomes. A bivariate Poisson model can be used to analyze such correlated count outcomes, which will provide in-depth insights and dynamics regarding new cases deaths and covariates.

To analyze these dependent count outcomes along with covariates, we use a bivariate generalized Poisson model ^[8]. This bivariate generalized Poisson regression model uses a marginal-conditional modeling approach based Poisson-Poisson distribution. Let $Y_{i1}$ and $Y_{i2}$ be the count responses of the daily new cases of COVID-19 infections and deaths, respectively, for the $i$th ($i = 1, 2, \ldots, n$) day. It can be shown that the joint distribution of $Y_{i1}$ and $Y_{i2}$ is ^[8]:

The joint mass function of a bivariate generalized Poisson model for count responses $Y_{i1}$ and $Y_{i2}$ can be expressed as,

where, $\ln {\lambda _{i1}} = x'_i{\beta _1}$, and $\ln {\lambda _{i2}} = x'_i{\beta _2}$ are link functions. In (2.1), $x_i' = (1, x_{i1}, \cdots, x_{ip})$ is the vector of the covariates and ${\beta _1}^\prime = ({\beta _{10}}, {\beta _{11}}, \cdots, {\beta _{1p}})$ and $\beta _2^\prime = (\beta _{20}, \beta _{21}, \cdots, \beta _{2p})$ are the regression coefficients corresponding to the new infections and deaths. It can be shown that the marginal means of $Y_{i1}$ and $Y_{i2}$ are $E(Y_{i1}) = \mu _{i1} = \lambda _{i1}$ and $E(Y_{i2}) = \mu _{i2} = \lambda _{i1} \lambda _{i2}$, respectively. It may be noted that the subscript $i$ of $\lambda$ is used to represents different subjects due to varying combinations of covariate values.

The log-likelihood function can be written as:

To estimate the regression parameters, we can take the first and second derivatives of the above log-likelihood function and develop the estimating equations as

and

The second derivatives are:

Consequently, The observed information matrix can be obtained as

and the approximate variance-covariance matrix for ${\hat{{\beta }}}' = \left({\hat{{{\beta }'}}_1, \hat{{{\beta }'}}_2 } \right)$ is $\widehat{Var}(\hat \beta) = {I_o}^{ - 1}$. Using the Newton-Raphson method, we can obtain the estimated parameters.

Now, using the Poison-multinomial relationship, we can predict the probabilities for bivariate Poison outcomes ^[9]. For notational convenience, the subscript $i$ is omitted in what follows. To estimate the joint probability, we do it in two steps. First, we need to calculate the marginal probability for each count ($m$) of outcome $Y_1$ $(m = 0, \cdots, k_1)$. Similarly, we need to estimate the conditional probability for each count ($s$) of $Y_2$ $(s = 0, \cdots, k_2)$ for a given value ($m$) of $Y_1$. The estimate of the marginal probability from the Poisson distribution can be obtained as:

Here, $\hat P(Y_1 = m\mid x)$ is the estimated probability for specific value $m$ for $Y_1$ and for given covariate value $x$. In other words, we are estimating the probabilities of each value of $Y_1$ by considering all the combinations of different covariate patterns.

For $Y_1 = m$, let $y_{m1}+\cdots+y_{ml}+\dots+y_{mn_m} = n_m$, where $y_{ml} = 1$ if $Y_1 = m$, $y_{ml} = 0$ otherwise, $m = 0, 1, \cdots, k_1$, $l = 1, \cdots, n_m$, and $\sum\limits_{m = 0}^{k_1}n_m = n$. The estimate of $P_m$ is

For the conditional probabilities of $Y_2 = s$ for any given value of $Y_1 = m$ the corresponding estimates of multinomial probabilities are

For $Y_2 = s$, let $y_{s1\mid m}+\cdots+y_{sh\mid m}+\dots+y_{sn_m\mid m} = n_{sm}$, where $y_{sh\mid m} = 1$ if $Y_1 = m$, $Y_2 = s$, $y_{sh\mid m} = 0$ otherwise, $m = 0, 1, \cdots, k_1$, $h = 1, \cdots, n_{sm}$, and $\sum\limits_{s = 0}^{k_2}n_{sm} = n_m$. For more comprehensive illustrations, readers are referred to ^[9]. Consequently, the joint probability of $Y_1 = m$ and $Y_2 = s$ can be estimated as follows:

We can easily predict the marginal, conditional, and joint probabilities using the fitted marginal and conditional models for new cases with different scenarios of outcomes and covariates. For marginal probabilities, we use the following equation

and for conditional probabilities we use

Then the joint probability can be estimated using the relation between Poisson and Multinomial described in Eqs (2.4) and (2.5).

3.   Analysis of World Health Organization COVID-19 Data

As of July 7, 2020, the World Health Organization (WHO) reported, globally, there have been 11,500,302 confirmed cases of COVID-19, including 535,759 deaths. We downloaded this map data from the WHO Coronavirus Disease (COVID-19) Dashboard. This data set includes daily new cases and deaths from 215 countries around the world by WHO region. The two-count outcomes, daily new cases and deaths, are defined as $Y_{i1}$ and $Y_{i2}$, respectively. The seven WHO regions EURO (European), AFRO (African), AMRO (American), EMRO (Eastern Mediterranean), SEAR (South East Asian), WPRO (Western Pacific), and Other are used as categorical covariates by considering EURO as the reference category. The duration (in days) from the first reporting date along with population density per $Km^2$ are also used as covariates in the bivariate Poisson regression model. To model the non-linear relationship with outcomes and assess the effect of duration more accurately, we also used the squared duration as a covariate in the bivariate Poisson model. The duration ranges from 0 to 176 days, and the density ranges from 0 to 26337 per $Km^2$. The ranges of these variables are vast, with long gaps between the values. For modeling purposes and to avoid convergence problems, we used deciles for each of these covariates. Besides, we recorded the number of daily new cases above 50 as fifty, and the daily number of deaths above 5 as five. Again, the reason is to avoid convergence problems due to the wide variation in the worldwide data. We fitted two models; in model 1, we only included the following covariates as main effect terms,

and in model 2, we included first-order interaction terms between covariates along with the main effects as follows:

To analyze the data for this paper, we used the bpglm R package to fit all the models ^[10].

The analysis results of the proposed bivariate Poisson model with the main effects only (Model 1) are presented in Table 1. Our results in Table 1 indicates that all the WHO regional indicators are showing a statistically significant relationship with daily new cases. Except for the EMRO region, all other areas are showing fewer daily new cases, on average, than the reference EURO region. The EMRO region is showing more daily new cases, on average, than the EURO region. This means the daily new cases of the COVID-19 infections in each region except ERMO are lower as compared to the EURO after the first case was diagnosed. This also suggests that compared to EURO, in the other regions the COVID-19 outbreak spread more slowly but in the ERMO region, infection spread more rapidly. Table 1 indicates that both linear and quadratic terms in duration are statistically significant in the model. The estimates tell us that, with the other variables fixed, that the number of new cases rises as the duration increases up to about 6, then decreases. Our results also indicates that population density has a negative significant effect on the number of cases. This is an indication that under proper safety measures it is possible to control the new cases even in highly populated areas. For the conditional model, with daily death counts as an outcome, except for three WHO region indicators, the other regions had fewer daily deaths, on average, than the EURO region. We are considering the adjusted p-values which adjust for over (or under) dispersion. The other covariates show a significant relationship with the daily number of deaths. However, we believe this may vary for country-specific data.

In Table 2, we present the results using the proposed bivariate Poisson model including the interaction terms as covariates with the main effects (Model 2). In this model, along with the main effects, we included interaction terms between duration and density and the interaction terms between duration and all the WHO region indicators. All the WHO region indicators (main effect terms) showed a significant ($p < 0.01$) reduction of daily new cases compared to the reference indicator EURO, except for Other regions. However, we need to be cautious with the interpretation of these terms because of the inclusion of the interaction terms involving region, which are all significant.

To compare the performance of the models, we have calculated the log-likelihood, Akaike information criterion (AIC), Bayesian information criterion (BIC) and Deviance using both Model 1 and 2. The results in Table 3 indicate that Model 2 is performing better. Figure 1 displays the predicted conditional probabilities for the given daily number of new cases. The number of new cases varies from 0 to 50 (50 represents fifty or more) and is shown on the x-axis. Each line in the figure represents the number of deaths, which ranges from 0 to 5 (5 represents five or more). The top line shows the probability of no deaths for the number of new cases. The probability shows a steady downward trend and goes close to zero as the number of new cases increases, as one would expect. The solid line for predicted probability for five or more deaths conditional on the number of daily new new cases rises sharply for fifty or more new cases, and the risk is around 0.29. This sharp increase is due to collapsing of 50 or more new cases as fifty to avoid estimation problems. We note that for some days with no new cases, deaths are also reported. Figure 2 presents the trajectories of joint probabilities for new cases and deaths. Overall the risk for both the new cases and deaths is gradually decreasing.

3.1.   Country specific analysis

We observed that the effect of duration on both the outcomes of daily new cases and daily deaths are positive for both the models (Tables 1 and 2). We hypothesized that these relationships might not remain the same for country-specific data. The reason is that the developed countries with advanced healthcare systems and with proper planning to handle the COVID-19 pandemic may observe different results. For example, the impact of duration on new cases and deaths may have negative relationships. Also, a longer duration may have a positive impact on daily new cases but a negative effect on daily deaths. However, they may have a positive impact on both the outcomes as we observed in model 1 and model 2.

To assess the country-specific impact of duration on daily new cases and deaths, we fitted separate models for selected countries. We used duration as the only covariate because the WHO region and density for each country are fixed values. Table 3 presents the bivariate Poisson regression results for the selected countries. For China, increased duration decreases new cases and deaths significantly. Similar patterns are found for Myanmar, but it was not significant for the marginal model. In the cases of Malaysia, Sri Lanka, Republic of Korea, Spain, Australia, and Morocco, duration reduces mortality ($p < 0.01$), but it is opposite for new cases. Duration reduces daily new cases ($p < 0.01$), for the African country Tanzania. The majority of countries observed increased in daily new cases and deaths.

From the country-specific analysis, we plotted the predicted risk of the number of deaths for fifty or more new cases for the selected countries (Figure 3 to Figure 6). The x-axis represents the number of deaths, and the predicted risks are on the y-axis. In Figure 3, the conditional probability (conditional on 50+ new cases) for five or more deaths showed an upward trend and much higher (0.22) for five or more fatalities for Italy, followed by Germany, Spain, and the UK, respectively. The risk for five or more fatalities (Figure 4) is highest for Australia, followed by France, the Netherlands, and New Zealand, respectively. The trend for the Netherlands is mostly flat and for New Zealand it is close to zero. Figure 5 presents the same trajectory for Canada, Japan, the Republic of Korea, and the USA. The path for Canada started ascending from 3 fatalities for the given new cases. The risk of deaths for five or more deaths is highest, around 0.12, followed by Japan, the Republic of Korea, and the USA, respectively. For China, the risk for five or more fatalities is highest, around 0.27 (Figure 6), followed by Malaysia, Morocco, and South Africa, respectively.

4.   Conclusions

The entire world is trying hard to reach the same goal: the new cases of COVID-19 need to go to zero. We can achieve this goal only if we can end the pandemic everywhere. To reduce COVID-19 infection, countries around the world are implementing various restrictions. These measures included travel restrictions and quarantine of both suspected individuals and subjects who have had close contact with suspected cases. Reducing new infections toward zero will reduce the spread of this disease and the number of deaths. In this paper, we have analyzed the COVID-19 infected new cases and deaths (two correlated count outcomes) and assessed the impact of duration, region, and density. The results showed that increasing duration significantly increased the number of daily new cases and deaths, which is observed for the world data. The impact of population density showed a negative effect on both outcomes. We believe this could only happen as some densely populated countries may have taken rigorous measures to fight this pandemic. We further investigated the impact of density, duration along with regional impact variable by introducing the interaction effect in the model. The impact of density showed a positive relationship with the number of deaths using the interaction model.

Regarding regional impact compared to the EURO, the rest of the regions showed a reduction of new cases except for Others. However, the AMRO and EMRO regions showed an increase in deaths compared to the reference EURO region. Regarding regional impact compared to the EURO, the rest of the regions showed a reduction of new cases except for Others. However, AMRO and EMRO regions showed an increase in deaths compared to the reference EURO region. These findings are based on the world's daily reported new cases and deaths.

The country-specific analyses reveal variations in conclusions regarding the relationship between correlated outcomes and duration. With the introduction of various stringent measures to fight this pandemic by many countries, we expect for an extended period, both new cases and death should go down. For some countries with a longer duration, both the new cases and deaths are reducing, e.g., China and Myanmar. On the other hand, a longer duration increases the daily new cases for the Republic of Korea, Malaysia, Sri Lanka, Morocco, and Iran, but reduces the number of deaths. I took out reference to Tanzania because P-values are large. Many other countries showed an increase in both new cases and deaths. These findings may be an indication of the measures to fight this pandemic are more advantageous for some countries but not for all. We believe that analysis with the availability of more detailed data along with related risk factors may provide an in-depth understanding regarding the dynamics.

Acknowledgments

This research was supported by grants from the Natural Sciences and Engineering Research Council of Canada (NSERC). The authors are grateful to the referees for their helpful comments on the original draft of our paper, which have served to greatly improve the presentation of our results.

Conflict of interest

The authors declare no conflict of interests.