Bayesian time-varying autoregressive models of COVID-19 epidemics.

The COVID-19 pandemic has highlighted the importance of reliable statistical models which, based on the available data, can provide accurate forecasts and impact analysis of alternative policy measures. Here we propose Bayesian time-dependent Poisson autoregressive models that include time-varying coefficients to estimate the effect of policy covariates on disease counts. The model is applied to the observed series of new positive cases in Italy and in the United States. The results suggest that our proposed models are capable of capturing nonlinear growth of disease counts. We also find that policy measures and, in particular, closure policies and the distribution of vaccines, lead to a significant reduction in disease counts in both countries.

INTRODUCTION

Time series of observed counts arise in a variety of contexts in which a number of events are recorded over time, including studies of incidences of a particular disease (Davis, 2003; Zeger, 1988). The reported occurrences are often expressed with daily, weekly, or monthly frequency. Disease modeling has significant practical importance for governments since it provides information that can improve decisionmaking. In addition, studying the correlation between counts and policy variables allows to better gauge the impact of control measures. This paper provides a Bayesian analysis of a time series of counts with the goal of providing reliable short‐term predictions of disease counts and to assess their correlation with explanatory variables. The time series of interest is the incidence of COVID‐19 infectious disease in Italy and Unites States, two paradigmatic countries, and the explanatory variables are the varying Non‐pharmaceutical interventions (NPIs) implemented by most countries to stop the virus propagation.

Since COVID‐19 first emerged in early 2020, it has infected nearly every country and claimed the lives of several thousands of people worldwide. So far, many countries have been locked‐down and stringent social distancing measures have been implemented to contrast the high incidence of the infection over time and its high economic consequences. Previous studies, mostly focused on the first wave of the pandemic (before July 2020), have shown the effectiveness of NPIs in pandemic mitigation (Cowling et al., 2020; Hao et al., 2020). They have also shown that the impact of NPIs may differ among countries with various country characteristics, like health capacity, population density, air temperature, etc. (Esteve et al., 2020; Marcon et al., 2020). Existing studies also indicate that vaccination is the most promising strategy for overcoming the epidemic. However, the unequal distribution and allocation of vaccines among nations and demographic groups may obstruct the path to herd immunity. It is crucial to understand how vaccination combined with NPIs reduces COVID‐19 transmission (Giordano et al., 2021).

The Poisson AutoRegressive (PAR) model is one of the most commonly used models for analyzing disease time series data, as it can well capture time dependence among observations. It can be seen as a special case of the generalized linear model (GLM) methodology (Fahrmeir et al., 2013; Kedem & Fokianos, 2002; Nelder & Wedderburn, 1972) where the link function is a linear function of the covariates, with coefficients that do not change over time. The advantage of this specification is the ease of model implementation and interpretation. However, when the impact of the epidemic changes over time, especially under the effect of policy measures, this type of model may be restrictive, and time‐varying coefficients may be needed.

A flexible way to introduce time‐varying parameters in a GLM model is through spline basis functions, as in the seminal work of Cleveland and Loader (1996). Spline basis functions have also been proposed to specify time‐varying coefficients in AutoRegressive (AR) time series models (Amorim et al., 2008; Cai et al., 1999; Huang & Shen, 2004). For this type of models, Zeger (1988) has suggested a generalized estimating equation methodology, which still is the most commonly used estimation method. Within the Bayesian framework, Chan and Ledolter (1995) proposed a Monte Carlo expectation maximization (MCEM) technique. More recently, Hastie et al. (2001) compared different estimation procedures. Moving to the more complex autoregressive moving average (ARMA) models, the only work that considers time‐varying parameters is that of Andrade et al. (2015).

None of these papers focuses on the time‐varying character of the coefficients. We contribute to fill the gap with a model that extends Fokianos et al. (2008) within a time‐varying coefficient framework. In doing so, we extend the time‐dependent contagion models recently proposed by Agosto et al. (2021), Geir et al. (2022), and Roy and Karmakar (2021) with the inclusion of time‐varying policy covariates. And we also extend the causal contagion models recently proposed by Bo et al. (2021) and Islam et al. (2020) with the inclusion of time‐varying covariate effects, which also allows to compare alternative time lags of policy measures. Note also that our proposed model allows to account for overdispersion of counts, as in the Markov‐switching Poisson GARCH model proposed by Roberts et al. (1998).

More specifically, We consider two classes of time series generalized linear models. The first is a class of AR models, in which the conditional mean of the Poisson process is linked to the past observations and to the potential covariate effects. The second one adds the dependence of the past values on the conditional mean. Since the mean of a Poisson model is the same as the variance, the latter can be seen as an Integer‐Valued Generalized Autoregressive Conditional Heteroskedasticity (INGARCH) model. Originally introduced by Ferland et al. (2006), these models have been thoroughly analyzed in Zhu (2012a, 2012b, 2012c) and Roberts et al. (1998).

We thus compare alternative models along three dimensions: AR versus INGARCH structures; time‐constant versus time‐varying coefficient, and models with or without covariates while testing for different time lags.

The rest of the paper is organized as follows. Section 2 describes the proposed models, Section 3 describes the Bayesian specification and inference. Section 4 describes the data. Section 5 presents the main empirical findings. Section 6 summarizes the highlights of the paper and suggests future research directions.

MODEL SPECIFICATION

The time dependence of COVID‐19 transmission represents a critical factor to understand the epidemic dynamics. The reference Susceptible‐infected‐recovery (SIR) model assumes an exponential growth model, which assumes that the logarithm of the counts follows a linear trend, summarized by R 0, a function of the ratio between the new cases originating on consecutive days. However, the evolution of the epidemics is characterized by time heterogeneity, which leads to a nonlinear growth. This derives, for example, from time‐varying testing and measurement procedures, among others.

From the previous consideration, it appears that the exponential growth assumption needs to be relaxed in favor of a time‐dependent modeling framework. Recently, Agosto et al. (2021) proposed a PAR model for epidemic counts that can take into account both short‐term and long‐term dependence on past contagion shocks. The model is specified by a Poisson conditional distribution for the COVID‐19 counts, along with a logarithmic link function which allows to model the effect of covariates without constraining them to assume positive values, as in a linear link model.

While the model has been quite useful to monitor contagion in the first wave of the pandemic, it has become less performing during the following times since, even if there is a nonlinear growth assumption, it is based on coefficients that do not change over time and, therefore, cannot easily adapt to regime switches.

In this paper, we aim to extend the fixed coefficients PAR model of Agosto et al. (2021) into a time‐varying coefficient framework, more flexible in model specification. Our modeling strategy is motivated to adapt to the rapid changes of the COVID spread. Hence, we propose a nonstationary time‐varying model, which is suitable for our analysis. The constraints of the parameter space in Equations (7) and (18) are motivated by the stationary conditions for the time‐constant case. This is fairly common in the time‐varying AR and INGARCH model literature (Dahlhaus & Rao, 2006; Davis & Mikosch, 2009; Fryzlewicz et al., 2008; Ferreira et al., 2017; Rohan & Ramanathan, 2013).

In detail, we consider two alternative specifications: one based on INGARCH models for counts, equivalent to the PAR models and other one based on the more parsimonious AR models.

AR models

The time series of counts reported at time (day) t is assumed to follow a Poisson distribution, conditional on the information up to : where is the σ‐field that contains all relevant information up to time t. Let = () be a time‐varying r‐dimensional covariate vector and let . A Poisson autoregressive model can then be specified, following Agosto et al. (2021), by the following link function: where , with , and are unknown parameters to be estimated from the data. From an interpretational viewpoint, μ is the intercept, describes the dependence of the expected count of cases on past observations and describes a dependence of the expected counts on r explanatory covariates at time point t.

The above model, as discussed in Agosto et al. (2021), allows to take into account time dependency of COVID‐19 counts, thereby reducing overdispersion. However, it remains the issue of a strong time heterogeneity of the policy measures taken by different countries, which can lead to an incorrect measurement of the effects of the covariates. To solve this issue, we propose to extend the AR model into a time‐varying coefficient version, in which the link function for becomes: where and represent, for , time‐varying splines. We recall that a spline is a continuous function that coincides with a polynomial on every subinterval in which it is defined. The coefficients of the polynomial differ from interval to interval, but the order of the polynomial is the same. We consider, for parsimony, cubic B‐splines.

The introduction of splines allows the effect of covariates to vary along different time intervals, trying to better capture the effect of time‐varying policy measures.

INGARCH models

The PAR model of Agosto et al. (2021) is an INGARCH model which allows to specify the logarithmic link not only as a function of the past observations but also as a function of the past conditional mean of the process.

More formally, the INGARCH model can be specified as follows: where the additional parameter , represents the dependence of the logarithmic link on the conditional mean of the Poisson process.

As done for AR models, we propose to consider a time‐varying extension of INGARCH models, described by the following: where, in addition to the parameters in Equation (3), represents, for , time‐varying parameters specified as functions of cubic B‐splines.

BAYESIAN INFERENCE

To better take model uncertainty into account, we consider a Bayesian specification for both AR (Section 3.1) and INGARCH models (Section 3.2). We specify the priors according to the parameter constraints in the time‐constant framework, putting some transformations on parameters to improve sampling efficiency. At the end, the posterior distribution is specified.

Equations refer to the time‐varying framework. However, for the time‐constant case, they can be easily recovered.

tvBAR models

Let θ = (μ, , ) indicate the vector of parameters of the AR model. Their parameter space is given by

The constraints introduced in Equation (7) on s are motivated by the time constant models where these conditions lead to stationarity and ergodicity. In the time‐varying AR literature, they are not uncommon (Dahlhaus & Rao, 2006; Fryzlewicz et al., 2008).

The prior distributions for the time‐varying Bayesian AutoRegressive (tvBAR) models can be specified, for each B‐spline basis expansion as follows: where are B‐spline basis functions. The subscripts i and k refer to the time‐constant dimension of the parameters (i.e., and ), while s refers to the dimension of the B‐splines functions.

For the unconstrained parameters and (which respectively refer to the intercept and the effect of explanatory covariates) we assume independent Gaussian priors in Equations (11) and (12) with variance equal to 100 to make the prior weakly informative. With reference to the constrained autoregressive parameter, we may specify a uniform prior in order to satisfy the stationarity conditions of time‐constant models.

However, from a computational viewpoint, this implies that the autoregressive parameters () are to be sampled in a bounded space between —1 and 1, reducing the efficiency of the Markov chain Monte Carlo (MCMC) posterior approximation. where . We therefore propose a parameter transformation, mapping the bounded parameter vector Θ to an unbounded parameter vector u to improve sampling efficiency. The proposed mapping function is the following: where . In detail, the final autoregressive parameter () to be sampled in Langevin Monte Carlo (LMC) algorithm is obtained as follows: where and . In this way, we put normal priors on the transformed B‐spline coefficients to be sampled within the LMC algorithm and then transformed back the parameters for the likelihood function (, where is the inverse of function g). This method results in better mixing of the posterior samples.

The log‐likelihood function of the model, combined with the proposed prior distributions, gives rise to the following log posterior distribution () for the tvBAR model:

tvBINGARCH models

Let θ = (μ, , , ), the vector of parameters of the INGARCH model. Their parameter space is given by

The constraints in Equation (18) are inspired from the stationary condition from time‐constant literature on INGARCH models (Geir et al., 2022).

For each B‐spline basis expansion, prior distributions can be specified as follows: where are B‐spline basis functions. The subscripts i, j, and k refer to the time‐constant dimension of the parameters (i.e., , , and ), while s refers to the dimension of the B‐splines functions.

Priors for μ and can be recovered from Equations (11) and (13). Moreover, we need to additionally estimate the parameter λ0 to start the recursion of the conditional mean (Equation 27).

Similar to what is done in Section 3.1, we propose a transformation of the time‐dependent parameters to ensure the condition about the absolute value of their sum (Equations 24 and 25 ), taking the logarithmic function (15) to map into an unbounded parameter space for sampling. The starting values for time‐dependent parameters (, ), together with their transformed version (, ) can be summarized as follows: where , , and with .

Specifically, , , , , and λ0 are the parameters to be sampled in MCMC iterations. After sampling, the parameters and are transformed back (, ) for the likelihood function. The posterior distribution for the Bayesian INGARCH (BINGARCH) model can then be obtained as follows for the time‐varying case:

To calculate the posterior distribution from the previous expressions, we employ the gradient‐based. The Hamiltonian Monte Carlo (HMC) algorithm is similar to Roy and Karmakar (2021). Tuning of HMC parameters is necessary to obtain a good acceptance rate. HMC has two parameters, the number of leapfrog steps and step length. Since gradient computation is expensive, we consider setting the number of leapfrog steps to one. In this way, HMC becomes LMC (Neal, 2011). According to Roberts et al. (1998), the optimal acceptance rate for LMC is around 0.57. We thus maintain an acceptance rate between 0.45 and 0.7 by adaptively tuning the step‐length parameter. Based on our results, this produces excellent mixing of the posterior samples (see Figure 3 for reference).

FIGURE 3

Trace plots of average likelihood function deviations across the MCMC chain for the BINGARCH model

We want to remark that, due to the increasing complexity of the parameter space for the INGARCH model, , , and λ0 are sampled together.

DATA

The proposed models have been applied to two types of data sources: COVID‐19 impact data (taken from the European Centre for Disease Prevention and Control (ECDC): www.ecdc.europa.eu/en/covid‐19/data) and policy intervention data (taken from Oxford Covid‐19 Government Response Tracker (OxCGRT): www.bsg.ox.ac.uk/research/research‐projects/covid‐19‐government‐response‐tracker).

From ECDC, we collected daily frequency COVID‐19 count data from two representative countries: Italy and the United States. We consider, as target response variable, the change in total positive cases, defined as the difference between the total counts of positives for the current day and those the day before. The starting point of the time series was chosen as the date when the country recorded more than 100 recorded cases.

We also obtained data on the taken policy measures from OxCGRT (Hale et al., 2020), a project from Oxford's Blavatnik School of Government, in which each government response to the pandemic is measured according to an ordinal scale that ranges from no measures to the most severe form of measure. Our primary interventions of interest are the measures defined as C—containment and closure policies as well as those defined as H—health system policies.

We remark that the epidemic and the policy measures taken to tackle the epidemics have been quite heterogeneous across countries, with NPI measures taken at the local authority level, varying across regions and across time. This heterogeneity suggests that these covariates have an intrinsic measurement error.

Tables 1 and 2 show a preliminary analysis of the OxCGRT variables of interest. Table 1 indicates sparse frequency data, mostly in correspondence with health system policies, which have less frequent changes over time, compared to containment and closure policies. Table 2 indicates high collinearity between policy measures, which can be explained by the fact that NPIs have not been implemented independently of each other but, rather, in a multivariate policy mix approach.

TABLE 1

Frequency table of NPI variables obtained from OxCGRT. Variables have been coded on an ordinal scale according to the severity of government response. NPIs with reference to Italy and United States are shown

NPI	Country	0	1	2	3	4	5
C1_School.closing	Italy	0	22	161	246	0	0
	United States	0	0	98	331	0	0
C2_Workplace.closing	Italy	0	33	214	182	0	0
	United States	10	50	265	104	0	0
C3_Cancel.public.events	Italy	0	0	429	0	0	0
	United States	0	53	376	0	0	0
C4_Restrictions.on.gatherings	Italy	28	0	113	17	271	0
	United States	2	1	9	65	352	0
C5_Close.public.transport	Italy	174	233	22	0	0	0
	United States	8	421	0	0	0	0
C6_Stay.at.home.requirements	Italy	88	84	236	21	0	0
	United States	6	175	248	0	0	0
C7_Restrictions.on.internal.movement	Italy	142	0	287	0	0	0
	United States	5	49	375	0	0	0
H2_Testing.policy	Italy	0	0	429	0	0	0
	United States	0	0	5	424	0	0
H3_Contact.tracing	Italy	0	0	429	0	0	0
	United States	0	429	0	0	0	0
H6_Facial.Coverings	Italy	26	0	0	0	403	0
	United States	1	25	2	86	315	0
H7_Vaccination.policy	Italy	293	0	110	26	0	0
	United States	280	75	31	20	16	7

TABLE 2

Most significant partial correlations between NPI covariates obtained from OxCGRT. A p‐value less than 0.001 is flagged with ***. Results with reference to Italy and United States are shown

Country	NPI(1)	NPI(2)	Pcorr
Italy	C1_School.closing	C3_Cancel.public.events	−0.623***
	C3_Cancel.public.events	C4_Restrictions.on.gatherings	−0.955***
	C4_Restrictions.on.gatherings	C5_Close.public.transport	0.464***
	C3_Cancel.public.events	C6_Stay.at.home.requirements	−0.637***
	C6_Stay.at.home.requirements	C7_Restrictions.on.internal.movement	0.58***
	C3_Cancel.public.events	H6_Facial.Coverings	−0.673***
United States	C3_Cancel.public.events	C7_Restrictions.on.internal.movement	0.591***
	C5_Close.public.transport	H2_Testing.policy	0.458***
	C2_Workplace.closing	H6_Facial.Coverings	−0.483***
	C1_School.closing	H7_Vaccination.policy	−0.641***

We have addressed these issues by grouping NPIs into homogeneous categories, in line with Cowling et al. (2020) and Davies et al. (2020) who concluded that combined interventions have been the most effective strategy. Bo et al. (2021) also showed that social distancing in combination with other NPIs has led to a significant decrease in COVID‐19 transmission.

More specifically, following Li et al. (2020), the original policy variables have been transformed from the ordinal to the binary scale, taking a value of “0” in case of no intervention or recommended intervention and a value of “1” in case of required intervention. We also grouped policies into five variables, following Islam et al. (2020): closure policies (X1_closures), restrictions on mass gatherings and events (X2_events), lockdown‐type policies (X3_lockdown), testing, tracing, and masks (X4_masks), and vaccination policy (X5_vaccines). The precise mapping between the OxCGRT variables and ours is shown in Table 3.

TABLE 3

Transformation of OxCGRT NPIs. Variables have been transformed from ordinal to binary scale taking the value of 0 in case of “no intervention” or “recommended intervention” and 1 in case of “required intervention.” Policy variables have then been grouped based on the thematic area of reference

NPI	Description	OxCGRT	Country	Value	Freq
X1_closures	Require closing of schools (or work from home) for some or all levels	c1; c2	Italy	0	247
			Italy	1	192
			United States	0	325
			United States	1	104
X2_events	Require canceling of events and restrictions on gatherings between 11‐100 people or less	c3; c4	Italy	0	141
			Italy	1	298
			United States	0	62
			United States	1	367
X3_lockdown	Require closing of public transport, require not leaving the house with some exceptions and internal movement restrictions in place	c5; c6; c7	Italy	0	172
			Italy	1	267
			United States	0	52
			United States	1	377
X4_masks	Open public testing, comprehensive contact tracing and facial coverings required outside the home	h2; h3; h6	Italy	0	36
			Italy	1	403
			United States	0	28
			United States	1	401
X5_vaccines	Availability of vaccines for two or more groups of people	h7	Italy	0	303
			Italy	1	136
			United States	0	355
			United States	1	74

EMPIRICAL FINDINGS

We first implemented the frequentist AR and INGARCH models using the tsglm function from the R package tscount (Liboschik et al., 2017). As shown in Table 4, the p order of the INGARCH structure has been selected according to the absolute error of the median predictions with respect to the true values (median absolute percentage error (MAPE))). Based on the mean incubation period of the disease, p equal to 1, 7, and 14 have been tested. Results show that p equal to 1 reports the lowest MAPE in both Italy and United States. Then, the time‐dependent model structure has been extended adding exogenous covariates for policy effect.

TABLE 4

Selection of the order of autoregressive dependence on past counts (p) for INGARCH structure based on MAPE

Country	Model	p = 1	p = 7	p = 14
Italy	AR	15%	24%	46%
	INGARCH	15%	24%	46%
United States	AR	11%	15%	29%
	INGARCH	12%	16%	28%

Additionally, we have implemented the same model structures within the Bayesian framework for which we have developed specific R functions. Since we extend the frequentist framework models by introducing B‐splines, we discuss about the selection of the number of knots for the B‐splines (Table 7) before presenting parameter estimates for both time‐constant and time‐varying coefficients.

TABLE 7

Selection of the number of knots for time‐varying Bayesian models based on the MAPE

Country	Model	Knots = 2	Knots = 4	Knots = 6	Knots = 8	Knots = 10
Italy	tvBAR	15%	14%	44%	21%	17%
	tvBINGARCH	100%	15%	100%	100%	18%
United States	tvBAR	39%	14%	10%	18%	12%
	tvBINGARCH	100%	11%	100%	100%	19%

Finally, we compared and discussed the results obtained in terms of predictive performance, using the last week of data as a test set. Point forecasts have been evaluated in terms of MAPE, whereas the whole distribution has been evaluated using scoring rules (Gneiting & Raftery, 2007), implemented in the R package scoringutils (Bosse et al., 2020).

The source code to reproduce the results is available as Supporting Information on the journal's web page.

Frequentist analysis

The results from the frequentist Poisson AR/INGARCH models, grouped by country, are shown in Table 5.

TABLE 5

Parameter estimates and standard errors for frequentist Poisson AR/INGARCH models for COVID‐19 data. When the p‐value against the null hypothesis of a zero coefficient is less than 0.001 the result is flagged with ***

Parameter	Country	Model AR	Model INGARCH
μ	Italy	0.367*** (0.006)	0.428*** (0.007)
	United States	0.526*** (0.003)	0.169*** (0.001)
α	Italy	0.962*** (0.001)	1*** (0.003)
	United States	0.954*** (0)	0.466*** (0.001)
β	Italy		−0.044*** (0.003)
	United States		0.52*** (0.001)

The AR estimates in Table 5 show a strong positive α1 for both countries, which indicates a positive dependence of counts over time. The INGARCH estimates in the same table show that, for Italy, α1 is much greater than β1. This indicates that short‐term dependence on past observed disease counts prevails over long‐term dependence: a result that could be explained by the containment and closure measures introduced very promptly in response to the COVID‐19 spread. Differently, in the United States, the estimates for α1 and β1 have a similar magnitude, revealing the persistence of long memory effects.

The impacts of policy measures are summarized in Table 6.

TABLE 6

Parameter estimates and standard errors of NPI covariates for the frequentist Poisson ARX/INGARCHX models applied to the COVID‐19 data. When the p‐value against the null hypothesis of a zero coefficient is less than 0.001 the result is flagged with ***

Parameter	Country	Model ARX	Model INGARCHX
X1_closure	Italy	−0.03*** (0.001)	−0.028*** (0.001)
	United States	−0.037*** (0.001)	−0.013*** (0)
X2_events	Italy	0.189*** (0.005)	0.224*** (0.005)
	United States	−0.043*** (0.002)	−0.058*** (0.001)
X3_lockdown	Italy	0.026*** (0.003)	0.05*** (0.003)
	United States	0.055*** (0.002)	0.059*** (0.001)
X4_masks	Italy	0.068*** (0.004)	0.102*** (0.004)
	United States	−0.04*** (0.002)	−0.124*** (0.001)
X5_vaccines	Italy	−0.035*** (0.001)	−0.042*** (0.001)
	United States	−0.048*** (0.001)	−0.02*** (0)

Table 6 shows that, for Italy, closure policies (X1_closures) and vaccines (X5_vaccines) have significant negative coefficients, meaning that these policy measures have significantly reduced infection counts. In contrast, restrictions on gatherings and events (X2_events), lockdown measures (X3_lockdown), and testing, tracing, and masks (X4_masks) reveal a significant positive impact on counts. This could be explained by the fact that these policy measures are typically imposed when the number of cases is high, leading to a “reverse effect” that leads to stronger policies when counts are higher.

The results for the United States also indicate negative significant effects for closure policies (X1_closures) and vaccines (X5_vaccines), but also for restrictions on gatherings and events (X2_events) and for testing, tracing, and masks policies (X4_masks). Only lockdown‐type policies (X3_lockdown) have a positive coefficient, which can be attributed to the “reverse effect” found for Italy. These results are in line with Hsiang et al. (2020) who have found effectiveness of most anticontagion policies in slowing down the exponential growth of the pandemic in the United States.

Bayesian analysis

To implement the proposed Bayesian time‐varying models, a preliminary selection of the number of spline knots is necessary. A large number of knots can cause overfitting, whereas a small number can result in underfitting. We select the number of knots on the basis of the results from the MAPE of the tvBAR and time‐varying Bayesian INGARCH (tvBINGARCH) models, applied to the training data, as shown in Table 7.

From Table 7, the best number of knots, which minimizes the MAPE for the tvBAR model is 4 for Italy and 6 for the United States; for the tvBINGARCH model, it is 4 for both countries. Figures 1 and 2 show the fitted cubic B‐splines of the daily new observed cases of COVID‐19 for Italy and for the United States, applying tvBAR models with, respectively, four and six knots.

FIGURE 1

Cubic spline fit of daily new observed cases of COVID‐19 in linear scale for Italy

FIGURE 2

Cubic spline fit of daily new observed cases of COVID‐19 in linear scale for the United States

Figure 2 shows that the peaks and the evolution of COVID‐19 data are well captured by a spline with six knots for the United States, which identifies three main peaks. Figure 1 shows that a simpler four knots spline is a good approximation for Italy, which identifies two main peaks, smoothing the two observed peaks in the second wave of the pandemic.

To derive the posterior distributions, the LMC algorithm has been run for 10,000 iterations, with the first 5000 discarded as burn‐in. All computations have been validated with convergence diagnostic assessments. For illustration, we report here the trace plots of the average likelihood function deviations along the MCMC chain for the coefficient functions (θ), calculated with the following formula: where represents the tth postburn samples with T = 5000 post burn iterations. For explanatory purposes, the trace plots for θ = with regard to the BINGARCH model are reported in Figure 3. They exhibit efficient mixing of the posterior samples.

The posterior distributions of the time constant parameters are given in Tables 8 and 9 for both countries.

TABLE 8

Posterior mean and 5%, 25%, 75%, and 95% posterior quantiles using time‐constant BAR/BINGARCH models for Italy

Parameter	Model	Posterior mean	5%	25%	75%	95%
μ	BAR	0.3837	0.3773	0.3807	0.3863	0.3903
	BINGARCH	0.4356	0.4259	0.4324	0.4399	0.4422
α	BAR	0.9608	0.9601	0.9605	0.9611	0.9614
	BINGARCH	0.4776	0.4773	0.4774	0.4778	0.4781
β	BINGARCH	0.5873	0.5824	0.5849	0.5895	0.5921

TABLE 9

Posterior mean and 5%, 25%, 75%, and 95% posterior quantiles using time‐constant BAR/BINGARCH models for the United States

Parameter	Model	Posterior mean	5%	25%	75%	95%
μ	BAR	0.5326	0.531	0.532	0.5329	0.5342
	BINGARCH	0.18	0.1784	0.18	0.1803	0.1806
α	BAR	0.954	0.9539	0.954	0.9541	0.9542
	BINGARCH	0.4666	0.4664	0.4665	0.4667	0.4668
β	BINGARCH	0.5183	0.5182	0.5183	0.5184	0.5185

From Tables 8 and 9, note that, for the United States, the frequentist findings in Section 5.1 are confirmed by the Bayesian estimates: both short‐term and long‐term dependence of counts are present. The same result is found for Italy, which is different from what occurred with the frequentist estimates, and consistent with the similarity between the two countries.

We now examine the posterior distributions of the time‐varying coefficients. Figures 4 and 5 depict the time evolution of the estimated μ and α1 from the Bayesian AutoRegressive (BAR) models for Italy and the United States. Figure 4 suggests that, for Italy, the estimated autoregressive coefficient of order 1 has a high impact on disease counts at the beginning of the period of analysis: the first wave of the pandemic. The impact starts to decrease when the number of cases decrease and increases in correspondence to the second wave of the pandemic. An almost opposite behavior is shown by the mean trend, which depends on the stock of cases, rather than on their flows. In April 2021, decreasing autoregressive coefficient estimates emerge, together with a low mean trend, resulting in a decreasing magnitude of the pandemic.

FIGURE 4

Time‐varying posterior coefficient distribution (tvBAR model) for Italy

FIGURE 5

Time‐varying posterior coefficient distribution (tvBAR model) for the United States

Figure 5 suggests that, for the United States, no evident turning points in the trend of estimated functions appear; the autoregressive coefficient is almost always close to 1 with a high mean trend. In the last period of analysis, the autoregressive component decreases and the mean trend becomes the prevailing factor.

Figures 6 and 7 report the time evolution of the μ, α1, and β1 estimated parameters from the BINGARCH models.

FIGURE 6

Time‐varying posterior coefficient distribution (tvBINGARCH model) for Italy

FIGURE 7

Time‐varying posterior coefficient distribution (tvBINGARCH model) for the United States

For Italy, Figure 6 confirms the results from the BAR model, also showing an increasing contribution of the conditional mean in the initial period of analysis, reaching the peak in the summer, where the contagion was decreasing and quite stable.

For the United States, Figure 7 reports a higher contribution of α1 over β1 during the first wave of the pandemic and the opposite behavior in the second wave. A growing mean trend is still present in April 2021.

The further step of the analysis involves adding NPIs covariates within the proposed modeling framework, comparing the results from the BAR and BINGARCH models. Tables 10 and 11 show the results, when time constant parameters are considered.

TABLE 10

Posterior mean and 5%, 25%, 75%, and 95% posterior quantiles of NPI covariates using time‐constant BARX/BINGARCHX models for Italy

Parameter	Model	Posterior mean	5%	25%	75%	95%
X1_closures	BARX	−0.0243	−0.0268	−0.0259	−0.0229	−0.0218
	BINGARCHX	−0.0223	−0.0237	−0.0229	−0.0217	−0.0207
X2_events	BARX	0.1916	0.189	0.1903	0.193	0.1942
	BINGARCHX	0.2327	0.2294	0.2312	0.2343	0.2352
X3_lockdown	BARX	0.0307	0.0271	0.0296	0.032	0.0329
	BINGARCHX	0.0478	0.0457	0.0469	0.0487	0.0503
X4_masks	BARX	0.0731	0.0703	0.0722	0.0743	0.0752
	BINGARCHX	0.1088	0.1064	0.1076	0.1101	0.1124
X5_vaccines	BARX	−0.0301	−0.0319	−0.0308	−0.0293	−0.0282
	BINGARCHX	−0.0378	−0.0399	−0.0388	−0.0369	−0.0354

TABLE 11

Posterior median and 5%, 25%, 75%, and 95% posterior quantiles of NPI covariates using time‐constant BARX/BINGARCHX models for the United States

Parameter	Model	Posterior mean	5%	25%	75%	95%
X1_closures	BARX	−0.0369	−0.0378	−0.0374	−0.0365	−0.036
	BINGARCHX	−0.0129	−0.0132	−0.0131	−0.0128	−0.0126
X2_events	BARX	−0.0432	−0.0469	−0.0443	−0.0417	−0.0405
	BINGARCHX	−0.056	−0.0565	−0.0561	−0.0559	−0.0557
X3_lockdown	BARX	0.04	0.0366	0.0384	0.0424	0.0432
	BINGARCHX	0.0517	0.0513	0.0515	0.0517	0.0525
X4_masks	BARX	−0.0421	−0.045	−0.044	−0.0405	−0.0386
	BINGARCHX	−0.1234	−0.1239	−0.1237	−0.1231	−0.1219
X5_vaccines	BARX	−0.0476	−0.0492	−0.0482	−0.047	−0.0463
	BINGARCHX	−0.0202	−0.0207	−0.0203	−0.02	−0.0197

The frequentist results in Table 6 are confirmed by the Bayesian posterior estimates in Tables 10 and 11. Table 10 suggests that for Italy the closures of schools and workplaces (X1_closures) and vaccination policy (X5_vaccines) are both effective in reducing COVID‐19 counts; whereas other variables have a positive sign. Table 11 confirms that, in the United States, most policies (with the exception of lockdown measures) reduce COVID‐19 counts.

When time‐varying parameters are considered, the median point estimates in Tables 10 and 11 are replaced by a range of estimates, corresponding to different time points. Figures 8 and 9 show the results for Italy and the United States, comparing the median estimates of the Bayesian AutoRegressive eXogenous (BARX) models with the estimates of the time‐varying Bayesian AutoRegressive eXogenous (tvBARX) models. In addition, we consider different lag length windows of 1, 7, and 14 days to identify a possible time lag between policy interventions and their epidemiological effect.

FIGURE 8

Time‐constant (BARX) versus time‐varying (tvBARX) posterior coefficient estimates for NPI covariates (Italy). The range of variation has been computed from the min‐max range of the posterior median over the time interval of interest

FIGURE 9

Time‐constant (BARX) versus time‐varying (tvBARX) posterior coefficient estimates for NPI covariates (United States). The range of variation has been computed from the min‐max range of the posterior median over the time interval of interest

From Figure 8, we can gauge that a time lag of about 1–2 weeks seems to increase the impact of covariates representing containment and closure policies (X1_closure), confirming that a large number of days are needed between the activation of this type of policy and the realization of its epidemiological effect. Figure 9 suggests a more immediate time impact of policies for the United States, as they lead to more negative estimates without time lag or with a delay of at most 1 day.

Both Figures 8 and 9 show a large variation in policy effects within the considered period. In particular, the effects of lockdowns (X3_lockdowns) and testing, tracing, and masks (X4_masks) vary much more than any of the other policy measures. This can be attributed to the presence of reverse causality effects, as policies have been implemented in response to the rise in pandemic counts.

Overall, it can be concluded that Bayesian estimates are in line with frequentist ones since the model structure is the same in the time constant case, but, additionally, incorporate greater uncertainty by means of the posterior distribution and the time‐varying coefficient estimates.

Model comparison

All models have been assessed in terms of the absolute error of the median forecasts, calculated as the absolute value of the difference between the true values and median posterior forecasts obtained from MCMC iterations. A high MAPE means that the median forecasts tend to be far away from the true values.

The accuracy of model prediction intervals has also been evaluated in terms of overprediction, underprediction, and bias, calculated from the scoring rules applied to the predictive Monte Carlo samples. They are defined with respect to the upper bound (overprediction) or lower bound (underprediction) of the prediction interval. The bias is a measure between —1 and 1 that expresses the tendency to underpredict (—1) or overpredict (1) the true values. All posterior predictions have been evaluated both on the train and on the test set (one week ahead). Tables 12 and 13 give the results of model comparison for Italy, whereas Tables 14 and 15 show the results for the United States.

TABLE 12

Model comparison for Italy with the Bayesian AR models

Model	Set	MAPE	Under	Over	Bias
BAR	Train	15%	672	672	0.13
	Test	36%	0	3576	1
BARX	Train	15%	646	647	0.053
	Test	28%	0	2934	1
tvBAR	Train	15%	362	369	−0.019
	Test	10%	10	950	0.357
tvBARX	Train	12%	128	155	0.002
	Test	15%	11	128	0.057

TABLE 13

Model comparison for Italy with the Bayesian INGARCH models. Note that the BINGARCH model specification has not been reported due to low‐performance metrics

Model	Set	MAPE	Under	Over	Bias
BINGARCHX	Train	15%	446	463	0.073
	Test	28%	0	1671	0.686
tvBINGARCH	Train	15%	572	579	0.004
	Test	9%	27	1718	0.554
tvBINGARCHX	Train	17%	360	399	−0.118
	Test	21%	440	479	−0.074

TABLE 14

Model comparison for the United States with the Bayesian AR models

Model	Set	MAPE	Under	Over	Bias
BAR	Train	11%	5352	5354	0.046
	Test	50%	1520	18098	0.714
BARX	Train	11%	5088	5088	0.047
	Test	31%	3003	12170	0.714
tvBAR	Train	10%	4426	4472	−0.034
	Test	32%	3455	6700	0.443
tvBARX	Train	11%	3942	3977	−0.037
	Test	29%	4410	7816	0.436

TABLE 15

Model comparison summary for the United States with the Bayesian INGARCH models

Model	Set	MAPE	Under	Over	Bias
BINGARCH	Train	12%	5110	5116	−0.05
	Test	33%	2684	12016	0.714
BINGARCHX	Train	12%	4998	5022	−0.013
	Test	26%	3706	8168	0.671
tvBINGARCH	Train	11%	2304	2219	−0.033
	Test	60%	131	4976	0.414
tvBINGARCHX	Train	10%	247	219	−0.017
	Test	35%	1353	9421	0.717

The results show that most of the proposed models fit the data quite well. For Italy, the tvBARX model seems to be the best, with a small MAPE (equal to 12%) on the training set and a bias almost equal to 0. The model also performs quite well on test data, with a 15% MAPE and a bias equal to 0.057. Similar conclusions are obtained for the United States, although the bias on the test set is higher (0.436), consistent with the higher number of cases with respect to Italy.

To further demonstrate the accuracy of our proposed models, Figures 10 and 11 plot the posterior median of the fitted and predicted values, together with the number of observed cases of COVID‐19 for Italy and the United States. The interval range (50–90%) has been added to both figures.

FIGURE 10

Predictions versus observed values over train and test data, with daily frequency (tvBARX model for Italy)

FIGURE 11

Predictions versus observed values over train and test data, with daily frequency (tvBARX model for the United States)

From Figures 10 and 11 both models appear to perform quite well.

We finally remark that our proposed TvBINGARCHX model was employed within the ECDC European Forecast Hub initiative and compared with many other models that attempted to predict counts of positive cases from 1 to 4 weeks ahead. The latter model was chosen since allows to account for both short‐term and long‐term dependence together with policy effects. The countries of interest were France, Germany, and Italy. Table 16 summarizes the evaluation metrics for the submission period that goes from March 2021 to November 2021, showing the best performance for Germany, with a bias equal to 0, slightly higher for Italy (0.6), and France (0.7). Real‐time forecasts of cases are really challenging, given data issues in COVID‐19 databases and intrinsic model error together with changing policies, population behavior, and testing procedures.

TABLE 16

Evaluation metrics of weekly forecasts made with the tvBINGARCHX model for Germany, France, and Italy for the European Forecast Hub project summarized throughout all the submission period (from 8 March 2021 tol 29 November 2021)

Location	Under	Over	Bias
DE	2	0	0
FR	0	1858	0.7
IT	0	3124	0.6

DISCUSSION

Motivated by the need to deliver reliable estimates of the COVID‐19 contagion dynamics, in this paper we have provided novel Bayesian Poisson autoregressive models, with time‐varying coefficients. While time‐varying coefficients allow to better take into account of nonlinear growth patterns, the Bayesian framework allows to incorporate the uncertainty of model outputs. The addition of policy intervention variables allows to better understand the timing, sign, and magnitude of the policy measures undertaken by the governments to face the pandemic.

The application of the model to two paradigmatic cases as Italy and the United States reveals that some combinations of policy interventions together with the first phase of the vaccination campaign are associated on average with a reduction in the incidence of COVID‐19. In detail, findings suggest beneficial effects of the closure of schools and workplaces for both countries. We also found that earlier implementation of restrictions results in a greater reduction in the incidence of COVID‐19 for Italy. The same conclusion seems not to be shared by the United States.

Accuracy evaluation of our model prediction intervals shows good in‐sample and out‐of‐sample performance. With respect to the latter, since our tvBINGARCHX model has also been employed within the ECDC European Forecast Hub initiative, low error with respect to 1–4 weeks ahead forecasts confirm the previous results.

Literature comparison

Agosto et al. (2021) extended the log‐linear Poisson autoregression (Geir et al., 2022) to epidemiologic contagion. The frequentist framework of our paper (Section 5.1) refers to this model, being the starting point of our analysis and the model against which comparisons have been made with the Bayesian framework. Unlike Agosto et al. (2021), policy covariates have been introduced.

The Bayesian framework instead introduces the time‐varying models of Roy and Karmakar (2021) and evaluates their performances against the extended structure with time‐varying policy covariates. Time constant coefficient models have also been added to have a direct comparison with Agosto et al. (2021) in the Bayesian framework. So, the model structure has been recovered and extended from these two contributions in terms of policy covariates.

On the other hand, in terms of policy evaluation, Bo et al. (2021) and Islam et al. (2020) have been recalled as a starting point for policy definition and evaluation. The effectiveness of NPIs was evaluated by Bo et al. (2021) in terms of time‐varying reproduction number (Rt) while Islam et al. (2020) evaluated policies in terms of counts, more similar to our approach. We confirm results of Bo et al. (2021) and Islam et al. (2020) on the effectiveness of distancing policies but we extend their work by introducing and quantifying the effectiveness of vaccines. Similar Islam et al. (2020), we introduce lag effects and a combination of policies, resulting in a better evaluation of policy effects.

Strengths and limitations of this study

The strength of our study relies on both methodological extension and epidemic monitoring. The structure of our proposed methods shows both extensions of previous works in terms of time‐dependence structure (i.e., time‐varying coefficients) and exogenous components of the model (i.e., policy covariates). Our model formulation is more flexible since it allows to consider both time dependence structure in terms of autocorrelation with previous counts but also with exogenous covariates. Alternative time lags allow to capture the timing of the effect of public interventions on the trajectory of the pandemic. Furthermore, the time‐varying coefficient framework allows both to model rapid changes of COVID‐19 spread and address overdispersion issues.

At the same time, our study is also subject to several limitations. First, data quality could be an issue for forecasting models, especially in terms of the measurement of policy covariates. We relied on OxCGRT variables to represent government interventions, but it is really challenging to collect information on the exact date, nature, and extent of policies undertaken by the different countries. Measurement with error could result in higher model bias. The discussion could be extended in the same way with regard to pandemic counts since day‐to‐day reports are usually made up of errors. In the end, real‐time forecasts of cases are really challenging, given data issues in COVID‐19 databases and intrinsic model error together with changing policies, population behavior, and testing procedures.

Implications for policy‐makers

Quantifying the impact of NPIs is a critical step in ensuring a strong public health response to COVID‐19 in the first phase of the vaccination campaign.

Risk factors for cluster formation are expected to be similar across countries, including closed and poorly ventilated locations together with crowded places where people have limited possibility for physical distancing. High‐risk environments include venues, events, and activities with certain environmental conditions. As a result, the findings suggest beneficial effects of the closure of schools and workplaces for both countries, the latter representing a possible option for reducing transmission but with a large‐scale economic and social impact.

Mobility plays an important role in human‐to‐human transmission. However, given the difference in terms of population size, limiting people's movement and gathering by requiring public transport closures in a combination of stay‐at‐home orders, and internal movement restrictions allows to decrease municipal (and regional) mobility fluxes and, consequently, the infection rate in Italy. On the other hand, requiring canceling gatherings and events at the stricter level (between 11 and 110 people or less) is effective in large countries such as the United States.

We also found that earlier implementation of restrictions results in a greater reduction in the incidence of COVID‐19 in Italy. The impact on the trajectory of infection changes for the United States, being more immediate. However, more factors need to be accounted for to obtain more accurate conclusions for both countries.

While NPIs reduce transmission rates, vaccination primarily reduces symptomatic disease, hospitalization, and death from infection, as well as infection rates to a lesser amount. Our findings show that vaccination is effective on epidemic evolution but in combination with some NPIs, advocating for the need to keep NPIs in place during the first phase of the vaccination campaign.

However, our results regarding several health system interventions need to be cautiously interpreted. The noneffectiveness of some NPIs can be attributed to several factors. For example, if the outcome of the study is related to the number of cases, as case detection improves through more effective testing and tracing, the number of cases reported will increase without representing a real rise of cases. Finally, for both countries, we have found a reverse effect of some policies, consistent with the fact that policy measures are triggered by the increase in COVID‐19 counts.

Further studies

Further extensions of the work in the paper should consider the inclusion of finer policy measurements, as new databases become available; and, possibly, the inclusion of spatial dependence in the model, especially for large countries such as the United States, as longer time series become available. For all countries, including Italy and the United States, it would also be important to include policy and contagion data at the subnational and regional level, if available, to better capture spatial variation.

CONFLICT OF INTEREST

The authors have declared no conflict of interest.

OPEN RESEARCH BADGES

This article has earned an Open Data badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available in the Supporting Information section.

This article has earned an open data badge “Reproducible Research” for making publicly available the code necessary to reproduce the reported results. The results reported in this article could fully be reproduced.

Supporting Information for this article is available upon request to the corresponding author.

Click here for additional data file.