Linking dynamic patterns of COVID-19 spreads in Italy with regional characteristics: a two level longitudinal modelling approach

1.   Introduction

Much of the effort on statistical modeling of coronavirus (COVID-19) spread has contributed to spreading patterns and forecasting of COVID-19 development; these patterns have been found to vary among locations ^[1,2]. As the survival time of coronaviruses on surfaces has found to depend on temperature, some researchers have explored the association of daily or cumulative confirmed cases with environmental factors ^[3,4,5]. To forecast COVID-19 related new cases and deaths, various authors have employed deep learning, artificial intelligence (AI) and time series approaches. For example, recurrent neural network (RNN) based deep learning techniques is proposed by ^[6] for forecasting COVID-19 related new cases. For predicting the dynamical behavior of COVID-19 cases various artificial intelligence (AI) based modeling techniques such as Bayesian regression neural network, cubist regression, k-nearest neighbors, quantile random forest, and support vector regression were employed in ^[7]. Autoregressive integrated moving average (ARIMA) based time series forecasting techniques were also incorporated in the literature to envision COVID-19 related new cases as well as deaths ^{[8,9,10,11,12]}. To predict the short-term spread of COVID-19, Singhal et al. ^[13] have proposed Gaussian mixture model-based techniques. Batista ^[14] has used susceptible-infected-recovered (SIR) model for estimating the final size of the COVID-19 pandemic spread. To accommodate time varying transmission and removal rate, and temporal trend of COVID-19 disease, Hong and Li ^[15] proposed time-dependent Poisson model. Bertozzi et al. ^[16] incorporated three models such as exponential growth, self-existing branching process and SIR models to predict the behavior of COVID-19 transmission in various stages of the diseases. These models in ^{[13,14,15,16]} are suitable for analyzing time series data without accommodating any covariate. But incorporating covariate information for predicting the COVID-19 disease spread may reveal important information. Furthermore, some researchers have studied the link between daily fatality rates with regional factors such as health resources, but found no significant factors ^[17]. As the spreading patterns of COVID-19 development vary a lot among locations, fitting regression models of daily confirmed cases or fatality rates directly with regional factors might not always reveal important relationships. In this study, we investigate the link between regional spreading patterns of COVID-19 development in Italy and regional factors in two steps. First, we extend the method of Zhang et al.^[1] to handle longitudinal count in order to characterize regional spreading patterns of COVID-19 daily confirmed cases; the varying growth and declining patterns as well as turning points among regions in Italy have been well captured by regional regression parameters of this model. We then associate these regional regression parameters with regional factors. The effects of regional factors on regional spreading patterns of COVID-19 daily confirmed cases have been effectively evaluated.

The coronavirus disease 2019 (COVID-2019) has so far spread to over 200 countries causing more than 20 million confirmed cases globally by August 15, 2020. In addition to countermeasures such as lockdowns and social distancing, it is also important to explore and reveal the relationship between virus transmission and potential environmental covariates such as temperature, to help us better understand the behavior of the virus, and further implement more efficient public interventions to contain the spread. Italy, as one of the countries first affected by the virus with a total number of more than 240,000 confirmed cases by the end of June, has experienced the transmission periods from exponential growth, turning point to flattened tail, which provides us an ideally mature data set for relevant regression analysis. Furthermore, potential covariates such as demographic and economic status vary from region to region in the country, and specially, as Italy has a geographically narrow and long shape from north to south, plus the mountainous terrain in midland, temperature conditions also vary at regional level. The variability of covariates makes regression analysis approach more feasible. In addition, the fact that Italy has carried out a unified countermeasure policy across the country helps minimize various effects brought by human intervention at regional level.

Our proposed hierarchical linear model deals with the analysis of the daily confirmed cases from February 24 to June 30, 2020 from various regions in Italy. To accommodate the heterogeneity among various regions, we consider regional level random effects. To depict infection-Time relationship, we adopt a power law with exponential cutoff (PLEC) function ^[1,18], which has exponential increase and decay terms. By incorporating the PLEC function into the hierarchical linear model, we can achieve the two major goals in this study: First, we construct a new hierarchical model to fit and predict the COVID-19 transmission trend in Italy, and compare it with the PLEC model ^[1,18]. Second, we conduct regression analysis on the relationship between COVID-19 transmission in Italy and spatially-varied covariates, including population density, GDP per capita and temperature. Our main interest in this study is the relationship between transmission and temperature.

Construction of the manuscript is as follows. After the proposed methodology in Section 2, we discuss our analysis results in Section 3. Some further discussions and conclusions of the proposed research are in Sections 4 and 5, respectively.

2.   Materials and methods

2.1.   Data

In this section, we briefly discuss the dataset used in the manuscript first before introducing the model. We collected daily confirmed cases data in Italy from Github ^[19] website from February 24, 2020 to June 30, 2020, from each of the 21 regions in Italy (including autonomous areas). There are situations when daily new cases are negative in certain regions due to data correction by the Italian authority, and we treat these observations as 0. Coordinates of each region were also gathered for further use. Italian demographic and economic data were collected from Wikipedia ^[20], including each region's population density and Gross Domestic Product (GDP) per capita. We also collected each region's historical daily average temperature from January 1 to June 30, 2020 on Wunderground website ^[21], and then calculated an average temperature in Fahrenheit during this period for each region. To better understand the geography, we present the map of Italy with various regions presented in Figure 1.

We consider regional factors such as population density, GDP per capita, and average temperature on the transmission of COVID-19 virus in Italy. To accommodate the different longitudinal patterns among the regions, we introduce three region-specific random effects, corresponding the three parameters in the PLEC function. Figure 2 demonstrates how fixed and random effects affect the transmission of the virus in different regions. Curves in different regions generally follow patterns that can be described by some PLEC curves, experiencing exponential growth and decay process, but the starting and turning point, as well as the growth rate of each curve differ across the regions in Italy, due to the influence from both fixed and random effects. By adding legends on average temperature and population density, some fixed effects from these two covariates can be visualized: Regions with higher average temperature appear to have curves at lower position, while curves of higher population density regions are at higher positions in the graph.

2.2.   Model formulation and assumptions

2.2.1.   Level-one model with PLEC function

In this section we present the PLEC function which will be used in the level-one model. It is noteworthy to point out that the PLEC function has a satisfactory performance on depicting COVID-19 transmission curves in different countries ^[18,22]. Following Ma ^[18], Wei and Zhang ^[22], the PLEC function has the following form

where $T$ is the time(in days) since the transmission begins, $I$ is the number of daily new infections, $c$, $w$, and $d$ are three positive parameters to be estimated from the data. Each of the parameters has its biological interpretation: $cT^w$ term is a power function with $w$ dominating the transmission growth with a multiplier $c$, and the exponential term $\exp(-dT)$ characterizes the declining trend of the transmission curve. Following ^[1], we can get the maximum value of infection by taking the derivative of $df(I)/dT$ and setting it to zero as,

which occurs at the turning point of

In order to incorporate the PLEC function with spatially-varied covariates in Italy, we assume that the daily confirmed case, which is a count response, follows a Poisson distribution as,

where

is the PLEC function in an alternative form, and $e^{c_j}$ corresponds to $c$ in Eq (2.1). In (2.4), $\mu_{ij}$, the mean value of daily new cases, on the ith level-one unit ($T$), representing 128 longitudinal daily observations from February 24 to June 30, nested within jth level-two (regional level) unit, where $j$ = 1, 2, ..., 21, representing 21 regions in Italy. The maximum value of infection and its turning point in jth region can be derived as

If we further apply log-transformation, Eq (2.4) can be expressed as

In (2.6), $T$, and its log form $\log(T)$ are level-one (day) predictors, which take identical day serial values across regions, with $-d_j$ and $w_j$ as coefficients, and $c_j$ can be treated as the intercept term.

2.2.2.   Level-two model

In level-two model, we consider the regression coefficients $c_j$, $w_j$, and $-d_j$ from level-one models as dependent variables incorporated with corresponding level-two covariates.

In our case, $\mathbf{x}_j$ is the vector of spatially-varied covariates at regional level, incorporated with level-one regression coefficients, and it may contain multiple covariates of interest. In addition, the $\mathbf{x}_j$ expression can be different for the three equations in (2.7). $\boldsymbol{\beta}_0$, $\boldsymbol{\beta}_1$ and $\boldsymbol{\beta}_2$ are coefficients of level-two covariates and $\alpha_0$, $\alpha_1$, $\alpha_2$ are corresponding intercepts. The residual terms $e_{0j}$, $e_{1j}$, and $e_{2j}$ are the random effects of the jth level-two unit (region) on the level-one parameters $c_j$, $w_j$, and $-d_j$ respectively. Following the derivation of Sullivan et al.^[23], We assume that the error terms $\mathbf{e}_j = (e_{0j}, e_{1j}, e_{2j})^\intercal$ follow a tri-variate normal with mean (0, 0, 0) and an unknown covariance matrix $G$.

2.2.3.   Combined model

We can derive our combined model by substituting (2.7) into (2.6):

Now both level-one and level-two covariates (time and spatially-varied covariates), interaction terms, and error terms are included in one model. It is then clear that the hierarchical formulation is equivalent to the random-intercept-and-random-slope model.

2.3.   Computational method

2.3.1.   Matrix form of model formulation

We derive our computational method in matrix form. First, we rewrite Eq (2.6) into the following expression,

where $\log(\boldsymbol{\mu}_j) = (\log(\mu_{1j}), \log(\mu_{2j}), \ldots, \log(\mu_{128, j}))^\intercal$, is the log form of mean response in 128 days in jth region, and $\boldsymbol{\gamma}_j = (c_j, w_j, -d_j)^\intercal$, is the parameter vector at level-one. $D_j$ is a known design matrix with following form identical across regions,

The first column in $D_j$ is a designed column for intercept term $c_j$, and third and second columns take day serial values and their log forms. We further rewrite level-two equation in the following matrix form

where $\boldsymbol{\beta}$ = $(\alpha_0, \boldsymbol{\beta}^{\intercal}_0, \alpha_1, \boldsymbol{\beta}^{\intercal}_1, \alpha_2, \boldsymbol{\beta}^{\intercal}_2)^\intercal$ is a vector of level-two coefficients, including both intercept and slope terms, which need to be estimated in the regression, and $\mathbf{e}_j = (e_{1j}, e_{2j}, e_{3j})^{\intercal}$ is a residual vector for random effects associated with three covariates from level-one model. As the known designed matrix of $\boldsymbol{\beta}$, $W_j$ takes the following form as

where $\boldsymbol{\beta}$, $\mathbf{x}_j$ and $\mathbf{0}$ have the same length, depending on the number of level-two covariates we select in $\mathbf{x}_j$. For a combined model, it has the following matrix form,

where

2.3.2.   Fixed effects estimation

We define $A_j = D_j W_j$, and it has a matrix expression as following:

Then we can use least squares method to estimate $\boldsymbol{\beta}$ by applying following formula:

$A$, $D$ are design matrices across regions, and $\hat{V}$ can be achieved with $G$'s maximum likelihood estimates. Variances of $\hat{\boldsymbol{\beta}}$ can be estimated by:

As our data is balanced with identical day counts in each region, G can be estimated by using closed-form maximum likelihood formulae. Random effects can be predicted by using best linear unbiased prediction (BLUP) method.

3.   Results

3.1.   Model selection

In this section we present our analysis results. We first include all potential spatially-varied covariates of interest into the level-two model, including regional average temperature, GDP per capita and population density. By inserting these level-two covariates into (2.7) and (2.8), we have the expressions of the first model as follows:

We use R function "glmer" in "lme4" package ^[24] to estimate the coefficients and predict random effects by specifying a nested structure at regional level on all three level-one covariates, which are the intercept term, day and log(day). Before applying "glmer" function, all level-two covariates were scaled so that the algorithm can converge more easily. We present the results of Model 1 in Table 1.

Based on the p-values, the three spatially-varied covariates or their interaction terms have significant effects at 0.1 significance level. Our results indicate that population density and GDP Per capita related terms have similar fixed effects on estimated mean value of daily cases, which is a sign that these two regional covariates could be correlated, thus we might need to drop one of them in further analysis. Moreover, we need to take the geographic factor into consideration when studying virus transmission, and certain geographic transmission trend could produce confounding variables. When COVID-19 broke out in Italy, it was believed that the virus hit the north first, and then it started to spread to the south. Therefore, if level-two covariates have a corresponding geographic trend, they could be confounded.

Figure 3 helps visualize the confounding effects caused by the geographic pattern. The size of the circles, with gradient color, represent number of cumulative cases, and the circles' growth trend can be observed by adding animation effects. Figure 3 implies that transmission associated with both average temperature and GDP per capita also follows a certain geographic pattern from north to south. To incorporate this potential confounding geographic effect, we introduce a new covariate "Distance" into the level-two model. The coordinates are used to calculate the geographic distance between each region and region "Lombardia", where people believe the virus first started to spread. The level-two covariate matrix in this "full model" becomes:

Results from the full model are displayed in Table 2. After bringing "Distance" covariate into the model, all GDP per capita related terms become insignificant at any conventional level of significance, so the GDP per capita covariate could be confounded with distance, and thus should be dropped.

After dropping the level-two covariate of GDP per capita, and using backward stepwise method to drop AverageTemp*$\log(T)$ term that is not statistically significant (P-value $= 0.257$), we have our final level-two model:

By substituting (3.4) into (2.8), we have the combined final model:

The final model is also equivalent to a random-intercept-and-random-slope model.

3.2.   Data analysis results

The values of estimated fixed effects using the above-mentioned final model are given in Table 3. It turns out that all remaining terms are significantly associated with COVID-19 transmission at 0.1 significance level, including the one of our main interests, temperature. Comparing the full model, final model and results in Table 1 from Model 1, we find that bringing the geographic covariate "Distance" makes all GDP per capita related terms insignificant, and dropping GDP per capita related terms does not change much on coefficients of population density related terms. As of our main interest, the sign and magnitude of coefficients related to average temperature effect have changed, compared with results in Table 1.

We also present predicted random effects results from "lme4" package ^[24] using BLUP method in Table 4. Random effects have a range of -3.98 to 4.22 on intercept, a range from -0.075 to 0.063 on "day" covariate, and a range from -2 to 1.63 on "logday" covariate. We can visualize the ranges of random effects on each of the level-one covariates in Figure 4, and further detect possible outliers of Region Molise, P.A. Trento and Umbria, as indicated in the figure. To check whether our predicted random effects satisfy the assumptions of the multivariate normal distribution, we draw normal probability plots and present results in Figure 5. Our graphs in Figure 5 show that the predicted random effects satisfy the multivariate normality assumption, and possible outliers from normality plots are consistent with the results in box plots.

3.3.   Model performance

To check the performance of the proposed model we use first days of observations as the training data to construct the model to predict the daily new cases in the rest of the days, and then check the deviation. With the variance-covariance matrix of both fixed effects and random effects given in "lme4" ^[24] results of our final model, We can further calculate a prediction interval for our proposed model by using R package "merTools" ^[25]. Figure 6 demonstrates the results using the first 60 days of data to construct the model and then to predict the rest of the days in each region.

Once we sum the fitted and predicted values from every region, we can compare our proposed model with models only using the PLEC function at national level. Based on the result from Figure 7, the proposed approach has less deviation in expected prediction than only using the PLEC function.

4.   Discussion

Results in Table 3 show that all three estimated effects associated with population density are significant at 0.01 significance level, the effect of "AverageTemperature" is significant at 0.01 significance level, and that associated with the interaction term "day*AverageTemperature" is significant at 0.1 significance level. So we have decided to keep "day*AverageTemperature" term in the model as it is an interaction term of our main interest and its p-value is close to 0.05. In order to link results from combined model to PLEC parameters, we group coefficients in Table 3 by their corresponding level-one parameters. Thus, we can directly interpret these covariates' fixed effects on parameters $c$, $w$ and $-d$ respectively.

Our results indicates that average temperature is negatively associated with parameter $c$. Holding other variables constant, the increase of every one scaled unit average temperature will lead to a 0.7244 decrease on parameter $c$, which downsizes the scale of the gamma curve. Also, due to the fact that average temperature is positively associated with parameter $-d$, the same increase on average temperature will cause $-d$ to increase 0.0069. As a result, with time parameter $T$ as a multiplier, $d*T$ will counter off the shrinkage effect from parameter $c$. Thus, the effect from change of average temperature on mean value of daily case is also related to time predictor, and we can derive its pattern by following approach. Based on (2.4), we set $\mu_0 = e^{c_0-d_0 T} T^{w_0}$ and $\mu_1 = e^{c_1-d_1 T} T^{w_1}$, where $\mu_1$ is the new function by only increasing $x$ scaled units average temperature ($x > 0$) while other variables are all held constant, thus, $c_1 - c_0 = -0.7244x$, $(-d_1)-(-d_0) = d_0-d_1 = 0.00699x$, and $w_0 = w_1$. Further, we take $\mu_1/\mu_0 = e^{(c_1-c_0)+T(d_0-d_1)} = e^{-0.7244x+0.00699xT}$. It is clear that $\mu_1 > \mu_0$ when $-0.7244+0.00699T > 0$, and it is free from $x$, which means at day 104, the effect from parameter $-d$ will exceed that from parameter $c$, and before this, it is always true that $\mu_1 < \mu_0$. This pattern is applicable to all regions. Similarly, we can also derive the effect of temperature change on the maximum value of infection and its turning point based on the result in (2.5). We set $\mu_{\max, 0} = e^{c_0-w_0}(w_0/d_0)^{w_0}$, $T_{\max, 0} = w_0/d_0$, and $\mu_{\max, 1} = e^{c_1-w_1}(w_1/d_1)^{w_1}$, $T_{\max, 1} = w_1/d_1$. Thus, $T_{\max, 1}/T_{\max, 0} = d_0/d_1 = d_0/(d_0-0.00699x) > 1$, which means higher average temperature leads to a delayed turning point, and $\mu_{\max, 1}/\mu_{\max, 0} = e^{c_1-c_0}(T_{\max, 1}/T_{\max, 0}) = e^{-0.7244x}(d_0/(d_0-0.00699x))^{w_0}$, suggesting that whether or not the ratio is larger than 1 depends on the value of parameter $-d$ and $w$ in each region. However, from the previous derived property, we know that as long as $T_{\max, 0}$ occurs before the 104th day, $\mu_{\max, 1} < \mu_{\max, 0}$, which means in every region in Italy, increase of average temperature will reduce maximum infections and postpone the turning point. Similarly, we can derive the properties of effects from population density change. From Table 3, we know that population density is positively associated with parameter $c$ and $-d$, and negatively associated with parameter $w$. Therefore, the final effect from population density on mean response depends on a comprehensive results from all three parameters. We assume an increase of $x$ scaled units ($x > 0$) on population density while all other variables are held constant. Thus $c_1-c_0 = 2.2358x$, $w_1-w_0 = -0.5811x$, and $d_0-d_1 = 0.0175x$. We then derive $\mu_1/\mu_0 = e^{(c_1-c_0)+(d_0-d_1)T}T^{w_1-w_0} = e^{2.2358x+0.0175xT}T^{-0.5811x}$. If we set $d(\mu_1/\mu_0)/dT = 0$ and solve for $T$, it can be shown that $\mu_1/\mu_0$ takes its minimal value when $T = 33.206$. Thus, $(\mu_1/\mu_0)_{\min} = e^{2.8169x}33.2057^{-0.5811x}$. We can further derive $(d(\mu_1/\mu_0)/dx)_{\min} = 0.7815e^{2.8169x}33.2057^{-0.5811x} > 0$ (strictly increasing) at $x$'s domain, and $(\mu_1/\mu_0)_{\min} = 1$ when $x = 0$, so it's proved that $\mu_1 > \mu_0$ when $x > 0$, which suggests the number of mean infections will always become larger when population density is increased.

Figure 8 demonstrates the pattern illustrated above in Region Lombardia. The pattern is consistent across regions in the country, as shown in Figures 9, 10 and 11.

One of the goals of Covid-19 modelling is to predict the turning point of the transmission curve before it occurs. In fact, PLEC model starts to predict reasonable turning points after detecting a "slowing down" trend due to its mathematical property, meanwhile, our proposed model may make this detection earlier in some regions by borrowing strength from other regions, which is an advantage of mixed-effects model.

5.   Conclusions

In this study, we have examined the association between regional spreading patterns of COVID-19 development in Italy and regional factors in two steps. In the first step, we have symbolized the regional spreading patterns of the daily confirmed cases of COVID-19 in Italy by a special patterned Poisson regression model for repeated daily counts. Our proposed regression model can capture the varying growth and declining patterns as well as turning points among regions in Italy by regional regression parameters. We then incorporate these regional regression parameters with regional factors as a second step to effectively evaluate the affect of regional factors on regional spreading patterns of COVID-19 daily confirmed cases.

Our two level longitudinal model performs better on prediction than models only using the PLEC function. This is because of the fact that two level model can incorporate regional covariates into the modelling approach, increasing the goodness of fit and predictive power. It also gives some mathematical properties linking to the level-one PLEC parameters, inferring that regional covariates average temperature and population density are associated with COVID-19 infection. This methodology could be useful in future epidemiological study when transmission pattern is potentially associated with social or environmental covariates. Neither our modelling approach nor PLEC function alone works well if we use data before the turning point to construct the model, due to the fact that lockdowns policy plays a crucial role affecting the curve. To incorporate this effect, we can use segmented models suggested by Zhang et al. ^[1] and Wei and Zhang ^[22]. Our analysis results indicate that geographical patterns perform an important role as it could confound other variables in epidemiological studies. More complicated models specifying geographic correlation could be explored on similar topics. A final point of caution is that this modelling procedure is not a designed experiment, which means that it can only explore associations between potential covariates and response, but can not reveal causal relationship, and there might be more confounding or lurking variables to be found.

Acknowledgements

This research was supported in part by grants from the Natural Sciences and Engineering Research Council of Canada. The authors thank the editor and three anonymous referees for their helpful comments that greatly improved the presentation of these results.

Conflict of interest

The authors declare there is no conflicts of interest.