Network-based time series modeling for COVID-19 incidence in the Republic of Ireland

Network-based time series models have experienced a surge in popularity over the past years due to their ability to model temporal and spatial dependencies, arising from the spread of infectious disease. The generalised network autoregressive (GNAR) model conceptualises time series on the vertices of a network; it has an autoregressive component for temporal dependence and a spatial autoregressive component for dependence between neighbouring vertices in the network. Consequently, the choice of underlying network is essential. This paper assesses the performance of GNAR models on different networks in predicting COVID-19 cases for the 26 counties in the Republic of Ireland, over two distinct pandemic phases (restricted and unrestricted), characterised by inter-county movement restrictions. Ten static networks are constructed, in which vertices represent counties, and edges are built upon neighbourhood relations, such as railway lines. We find that a GNAR model based on the fairly sparse Economic hub network explains the data best for the restricted pandemic phase while the fairly dense 21 -nearest neighbour network performs best for the unrestricted phase. Across phases, GNAR models have higher predictive accuracy than standard ARIMA models which ignore the network structure. For county-specific predictions, in pandemic phases with more lenient or no COVID-19 regulation, the network effect is not quite as pronounced. The results indicate some robustness to the precise network architecture as long as the densities of the networks are similar. An analysis of the residuals justifies the model assumptions for the restricted phase but raises questions regarding their validity for the unrestricted phase. While generally performing better than ARIMA models which ignore network effects, there is scope for further development of the GNAR model to better model complex infectious diseases, including COVID-19.


Introduction
In recent years, statistical models which incorporate networks and thereby acknowledge spatial dependencies when predicting temporal data have experienced a surge in popularity (e.g., Knight et al. 2019Knight et al. , 2016;;Urrutia et al. 2022).Knight et al. (2016) developed a generalised network autoregressive (GNAR) time series model which incorporates a secondary dependence in addition to standard temporal dependence.The secondary dependence is captured in a network.In Knight et al. (2016), the proposed networkbased time series model is leveraged to predict mumps incidence across English counties during the British mumps outbreak in 2005.As graph to be associated with the Mumps network time series, Knight et al. (2016) chose a "county town" for each county and connected all towns which were less than a radius of a fixed number of kilometers away from each other.
Similar to mumps, COVID-19 is a highly infectious disease spread by direct contact between people (Nouvellet 2021).Human movement networks have been extensively relied upon to explain COVID-19 patterns (e.g.Jia 2020;Kraemer 2020;Li et al. 2021;Mo 2021;Nouvellet 2021;Sun et al. 2021;Wu et al. 2020).Therefore, it is a natural conjecture that such movement networks may help predict the spread of COVID-19.To investigate, this paper • Fits GNAR models to predict the weekly COVID-19 incidence for all 26 counties in the Republic of Ireland, exploring different network constructions; • Assesses the prevalence of a network effect in COVID-19 incidence in Ireland and the suitability of GNAR models to predict epidemic outbreaks as complex as COVID-19; • Investigates the influence of changes in inter-county mobility, due to COVID-19 restrictions, on the performance of GNAR models as well as on the model parameters and hyperparameters.
The GNAR model is chosen because multivariate time series are often modelled by vector autoregressive (VAR) models.General VAR models are very flexible but require a large number of parameters to be estimated.The GNAR model which we employ here is a special case of a VAR model.It reduces the number of parameters to be estimated by restricting attention to edges in a network; in the case of a complete graph, the VAR model and the GNAR model coincide.Our overview of network-based time series models, given in Supplementary Material B.2, concludes that many network-based time series models can be conceptualized as a special case of the GNAR model, or are more restrictive with respect to the temporal-spatial dependencies they can model.Moreover, as a VAR-type model, the GNAR model inherits the well-understood VAR model framework, including parameter estimation via least squares, and model selection based on the BIC; for a survey see for example Lütkepohl (2005).These methods yield confidence sets for parameter estimation, which can inform analysis as well as policy development in a quantitative fashion.In contrast, deep learning approaches such as developed in Park et al. (2024) for predicting rental and return patters at bicycle stations do not come with such theoretical guarantees.
In addition to a distance-based network as chosen in Knight et al. (2016), in this paper we explore a collection of network models which are motivated by potential movements of individuals.The movement networks are constructed according to general approaches based on statistical definitions of neighbourhoods as well as approaches specific to the infectious spread of the COVID-19 virus.In abuse of notation, we call these networks COVID-19 networks, although they are only meant to reflect possible transmission routes of the disease.For each network, we select the best performing hyperparameter values, to predict COVID-19 incidence by a GNAR model, using the Bayesian Information Criterion (BIC).By splitting the available Irish data into two phases of the pandemic, restricted and unrestricted, we are able to investigate the potential change in the temporal and spatial dependencies in COVID-19 incidence between the two phases.
Overall, our findings are that while there is a clear network effect, the performance of the optimal GNAR model varies little across different network architectures of similar network density.GNAR models indicate higher predictive accuracy than ARIMA models on a country level, since they account for inter-county dependencies.On an individual county level, the variability of predictive performance is high, resulting in similar performance of ARIMA and GNAR models for some counties, while for others the GNAR model consistently outperforms the ARIMA model.The GNAR model seem better suited to predictive COVID-19 incidence in restricted pandemic phases than in unrestricted pandemic phases; the latter may be related to some of the model assumptions possibly requiring an adaptation as well as an increase in noise during unrestricted pandemic phases and high fluctuation of COVID-19 case numbers.Moreover, the classical VAR model, which is the GNAR model on the complete graph, does not perform as well as the GNAR model with an underlying network that has fewer edges, illustrating the value of using a GNAR model.This paper is organised as follows.Section 2 introduces the data set.The methodology for network construction and for network-based time series modeling is described in Sect.3. Section 4 provides an exploratory data analysis, while the model fit is shown in Sect.5.1.The conclusions for the different pandemic phases are found in Sect.5.2.The results are discussed in Section 6.The Supplementary Material includes a literature review on alternative models for incorporating temporal and spatial dependencies, visualisations of the COVID-19 networks, as well as details on the performance of GNAR models for predicting COVID-19 incidence.

The Irish COVID-19 data set
By March 2023, the Republic of Ireland (abbreviated in this paper as Ireland) had recorded a total of 1.7 million confirmed COVID-19 cases and 8,719 deaths since the beginning of the pandemic (Health Protection Surveillance Centre 2022a).The Health Protection Surveillance Centre identified four main variants of concern for the COVID-19 virus in Ireland (Health Protection Surveillance Centre 2022b), in addition to the original variant: Alpha from 27.12.2020,Delta from 06. 06.2021, Omicron I from 13.12.2021, and Omicron II from 13.03.2022(Figure 1a in Health Protection Surveillance Centre (2022b)).The open data platform (Government of Ireland 2022;Ordnance Survey Ireland 2024) by the Irish Government provides weekly updated multivariate time series data on confirmed daily cumulative COVID-19 cases for all 26 Irish counties, starting from the beginning of the pandemic in February 2020.A COVID-19 case is attributed to the county in which the patient has their primary residence.1To our knowledge, spatially more detailed COVID-19 data is not available for Ireland.The limited granularity makes it difficult to implement fine-scale spatial models.The GNAR model was originally demonstrated using county-level data, indicating its potential for modeling infectious diseases at lower resolutions.The cumulative case count is given for 100,000 inhabitants.Age profiles vary across counties and COVID-19 infection rates are age dependent.Hence, the cumulative case count is adjusted for age distribution according to the 2016 census of Ireland, to ensure inter-county comparability (Central Statistics Office 2016).In our dataset, the first COVID-19 case was registered in Dublin on 02.03.2020 and the last reported date is 23.01.2023, spanning a total of 152 weeks (Ordnance Survey Ireland 2024).From 20.03.2020 onward, COVID-19 cases were recorded in every Irish county.The daily COVID-19 data is aggregated to a weekly level to avoid modelling artificial weekly effects (Kubiczek and Hadasik 2021;Sartor 2020).Due to delayed reporting during winter 2021/22, the weekly COVID-19 incidences from 12.12.2021 to 27.02.2022are averaged over a window of 4 weeks (Health Protection Surveillance Centre 2021, 2022c;Wei 2006).
The main COVID-19 regulations restricting physical movement and social interaction between Irish counties (Brennan 2022;Loughlin 2022;McQuinn et al. 2021) are used to naturally split the data into five sequential subsets, where the COVID-19 incidence is small at the beginning of the pandemic and shows a clear increasing trend over time.We denote by σ the average standard deviation in COVID-19 incidence across the 26 Irish counties within the considered data subset.

Methodology
Let G = {V, E} denote a fixed, simple, undirected, unweighted network with vertex set V containing N vertices and edge set E ; an edge between vertices i and j is denoted by i ∼ j .The neighbourhood of a subset of vertices A ⊂ V is defined as the set of neighbours outside of A to the vertices in A, N (A) = i∈A {j ∈ V\A : i ∼ j}.The set of r th -stage neighbours, or the r th -stage neighbourhood, for vertex i ∈ A is defined recursively as N (0) (i) = {i} and

COVID-19 networks: constructions and properties
The key to the GNAR model is the network.The true network underlying the data generating process (in this case, who infected whom in the spread of the disease) is usually unknown.Ideally, expert knowledge can be leveraged to build a network that captures the relationship between vertices, representing the subjects of interest.Networks to model the spread of an infectious disease, such as COVID-19, are frequently modelled off of human mobility patterns, which are considered to have a shaping influence on disease spread (e.g.Colizza et al. 2006;Jia 2020;Li et al. 2021;Mo 2021;Sun et al. 2021).To our knowledge, detailed information on weekly population flow between Irish counties is unavailable.Hence, we construct implicit COVID-19 transmission networks (COVID-19 networks hereafter) based on geographical approaches, in line with Knight et al. (2016).
In the Railway-based network, an edge is established between two counties if there exists a direct train link between the respective county towns (without change of trains) and the county towns are closest to each other on this train connection.The Queen's contiguity network connects each county to the counties it shares a border with (Sawada 2022).The Economic hub network adds an additional edge between each county and its nearest economic hub to the Queen's contiguity network: Dublin, Cork, Limerick, Galway or Waterford (Gardham 2022).To measure the distance to the nearest economic hub we use the Great Circle distance d C (i, j) , the shortest distance between two points on the surface of a sphere (Weisstein 2002).For two points i, j with latitude δ i , δ j and longitude i , j on a sphere of radius r > 0, The K-nearest neighbours network (KNN) connects a vertex with its K nearest neighbours with respect to d C (Bivand 2022;Eppstein et al. 1997).The distance-based neigh- bour network (DNN) constructs an edge between counties if their Great Circle distance d C lies within a certain range [l, r]; this construction is similar to the one used in Knight et al. (2016).For the COVID-19 network, we set l = 0 and consider r a hyperparameter, chosen large enough to ensure that no vertex is isolated.The maximum value for r is determined by the largest distance between any two vertices, for which it returns a fully connected network (Bivand et al. 2013).
In addition to these geographical networks, the Delaunay triangulation constructs geometric triangles between vertices such that no vertex lies within the circumsphere of any constructed triangle (Chen and Xu 2004), thus ensuring that there are no isolated vertices.The Gabriel, Sphere of Influence network and Relative neighbourhood are obtained from the Delaunay triangulation network by omitting certain edges.In a Gabriel network, vertices x and y in Euclidean space are connected if they are Gabriel neighbours; that is, where d(x, y) = n i=1 (x i − y i ) 2 denotes the Euclidean distance.In a Sphere of Influ- ence network (SOI), long distance edges in the Delaunay triangulation network are eliminated and only edges between SOI neighbours are retained, as follows.For x ∈ V and d x the Euclidean distance between x and its nearest neighbour in V , let C x denote the circle centred around x with radius d x .For y ∈ V , the quantities d y and C y are defined analogously.Vertices x and y are SOI neighbours if and only if C x and C y intersect at least twice, preserving the symmetry property of the Delaunay triangulation (Bivand et al. 2013).The Relative neighbourhood network only retains edges between relative neighbours, The Relative neighbourhood network is contained in the Delaunay triangulation, SOI and Gabriel network, and is the sparsest of the four networks (Bivand et al. 2013).Finally, the Complete network represents the homogeneous mixing assumption, where all countries are connected (Bansal et al. 2007).
Figure 1 shows the Economic hub network and the KNN network ( k = 11 and k = 21 ) for Ireland.Figures of the other networks are found in the Supplementary Material A; network summaries are provided in Table 1.
The networks are created based on the literature on spatial modelling; Bivand et al. (2008, pp. 239-251) suggests the Delaunay network and its variants, Queen's contiguity network, distance-based networks, and K-nearest neighbour networks.In De Souza et al. (2021) it was found that infrastructure has an effect on the spread of COVID-19 in Brazil; here we use the railway network as infrastructure network.The Economic hub network is motivated by the idea of including a proxy of commuter flows in the network construction, as commuting to work has been shown in Mitze and Kosfeld (2022) to be related to the spread of COVID-19 in Germany.

Generalised network autoregressive models
Network-based time series models incorporate non-temporal dependencies in the form of networks in addition to temporal dependencies commonly modeled in time series models (Knight et al. 2019(Knight et al. , 2016;;Zhu et al. 2017).In contrast to standard time series methodology and spatial models (Box et al. 2015;Hamilton 2020;Wei 2006), network-based time series models are not limited to geographic relationships but can incorporate any generic network.As COVID-19 is an infectious disease with spatial spreading behaviour, warranting the construction of networks based on spatial information, we use terms relating to spatial dependence in our exposition.Other types of dependence could easily be incorporated in the model through networks which reflect the hypothesised dependence.
The global α generalised network autoregressive models GNAR (p-s 1 , ..., s p ) models the observation X i,t for a vertex i at time t as the weighted linear combination of an autore- gressive component of order p and a network neighbourhood autoregressive component of a certain order, also called neighbourhood stage; for i = 1, . . ., p , the entry s i gives the largest neighbourhood stage considered for vertex i when regressing on up to p past values.In our analysis, X i,t denotes the 1-lag difference in weekly COVID-19 incidence over time t for county i.The effect of neighbouring vertices depends on some weight ω i,q .A GNAR model GNAR(p-s 1 , . . ., s p ) has the following form, where ε i,t ∼ N (0, σ2 i ) are uncorrelated. 2As weights ω i,q we choose the normalised inverse shortest path length (SPL) weight, where d i,q denotes the shortest path length (SPL) (Knight et al. 2016); in connected networks, 1 ≤ d i,q < ∞ for i = q .For i ∈ V and q ∈ N (r) (i) , we thus take, If the network is dynamic over time instead of static, the weights can be constructed to depend on time (Knight et al. 2016).
The general GNAR model relies on vertex specific coefficients α i,j , instead of vertex unspecific autoregressive coefficients, α j , indicating vertex specific temporal depend- ence.This modification is comparable to including individual random effects in regression models.
(1) To fit a GNAR model, we must choose two hyperparameters, the lag p, or α-order, and the vector of neighbourhood stages, s = (s 1 , ..., s p ) , also called β-order.They can be determined either through expert knowledge, e.g. on the spread of infections, or through a criterion-based search (Knight et al. 2016).The model coefficients are computed via Estimated Generalised Least Squares (EGLS) estimation (Knight et al. 2016;Leeming 2019;Lütkepohl 1991). 3

GNAR model selection and predictive accuracy
For our analysis of the Irish COVID-19 data, model selection, i.e. the choice of α -and β -order, is performed by minimizing the Bayesian Information Criterion (BIC) (Knight et al. 2016).The BIC avoids overfitting by penalizing the observed likelihood L by the dimensionality of the required parameters (Schwarz 1978).For a sample X of size n and a parameter The GNAR package assumes Gaussian errors (R Package Documentation 2022); under this assumption, the BIC is consistent.This assumption could be weakened; it can be shown that the BIC is consistent for the GNAR model (1) if the error term is i.i.d. with bounded fourth moments (Leeming 2019;Lütkepohl 2005;Lv and Liu 2014).
The predictive accuracy of a GNAR model is measured by the mean absolute scaled error (MASE), due to its insensitivity towards outliers, its scale invariance and its robustness (Hyndman and Koehler 2006).MASE is defined for each county i as the ratio of absolute forecasting error εi,t = |X i,t − Xi,t | divided by the mean absolute error between true and a naive 1-lag random walk forecast for the entire observed time period [1, T] (Hyndman and Koehler 2006;Urrutia et al. 2022);

The weekly incidence differences
The GNAR model requires stationary data.Stationary data has no trend over time and is homogeneous, i.e. has time-independent variance (Knight et al. 2016).The weekly COVID-19 count is clearly not stationary, as it shows an increasing trend in Fig. 2. To remove any linear trend, we perform 1-lag differencing on the weekly COVID-19 incidence for the 26 Irish counties, resulting in the incidence difference, (1-lag) COVID-19 ID, between two subsequent weeks (Montgomery et al. 2015).We assess the stationarity by applying a Box-Cox transformation to each data subset.Figure 12 in the Supplementary Material C, indicate that no further transformation is required.

Constructed networks
For the COVID-19 KNN network, neighbourhood sizes sequencing from k = 1 to the fully connected network, k = 25 , with step size 2, are considered.The minimal distance There is considerable variability in network characteristics, Table 1, in particular regarding the network density.The KNN and DNN networks for the abovementioned hyperparameters have much larger average degree than the other networks; the sparsest network is the Relative neighbourhood network.Consequently, the SPL is shortest in the denser DNN and KNN networks.The Railway-based network has the longest average SPL due to its vertex chains and the low number of shortcuts between counties.For the Economic hub network, the introduction of shortcuts to the economic hubs leads to a decrease in average SPL, i.e. the disease spreads quicker, compared to the Queen's contiguity network.The Gabriel network is sparser than the SOI network, with slightly longer average SPL.Deleting long edges in the Delaunay triangulation network to obtain the SOI network decreases the average degree and the average local clustering coefficient, but increases the average SPL.The Queen's network, the Economic hub network, the Delaunay network, the Gabriel network, and the SOI network show small world behaviour, i.e. high clustering with short SPL.To assess small world behaviour, the average SPL and average local clustering for a network is compared to a Bernoulli Random graph G(n, m) with identical size n and number of edges m as the network under investigation.The Railway network has much larger average SPL than G(n, m), while the dense KNN and DNN networks have almost the same average SPL and local clustering coefficient as the G(n, m) network.
Although there are differences in the detailed summary statistics, the networks can be clustered according to density and average local clustering coefficient; we use the k-means algorithm, running 10 randomisations to ensure robustness over a range of k = 1, . . ., 10 .The corresponding elbow plot implies two clusters of COVID-19 networks.As evident in Fig. 3, one cluster has high density and high average local clustering coefficient, while the second cluster has low density and low to medium average local clustering coefficient.

Spatial effects
Intuitively, if spatial correlation is present in a network, the closer in SPL two vertices are, the more highly correlated their COVID-19 incidences.Moran's I quantifies spatial correlation by estimating the average weighted correlation across space (Cliff and Ord 1981;Moran 1950;Zhou and Lin 2008).Let t ∈ T and x (t) i denote the COVID-19 ID for county i at time t, where W 0 = N i,j=1 w ij for normalisation.For non-neighbours, the weights are zero, i.e. ∀r : j � ∈ N (r) (i) : w ij = 0 .Here we use weights w ij = e −d ij where d ij denotes the SPL between vertex i and j (Coscia 2021).The spatial dependency between counties varies strongly over time for every network, see Fig. 4 and Figure 11 in Supplementary Material A.1.Peaks in Moran's I coincide with peaks in the 1-lag COVID-19 ID at the beginning of the pandemic as well as during the winters 2020/21 and 2021/22.The introduction of restrictive regulations, e.g.lockdowns, shows a decreasing trend in Moran's I while the ease of restrictions from summer 2021 onward has lead to an increasing trend in Moran's I.This indicates a network effect in the data, which is associated with the intercounty mobility, and becomes particularly evident after the official end of pandemic restriction in March 2022.To further assess the presence of non-linear spatial correlation, we also apply Moran's I to the ranks of the COVID-19 ID at each time point t over the duration of data observation.The rank-based Moran's I follows the same trajectory, with less extreme peaks, as evident in Fig. 5.

GNAR model fitting
To assess the benefit of accounting for network effects, the GNAR model is compared to a standard county-specific ARIMA model which is allowed to include seasonal effects. 5 The GNAR model allows for more flexible spatial-temporal dependencies than other network-based time series models as detailed in Supplementary Material B.2, including the ARIMA models, with a network specific selection of α -and β-order.The models are selected by choosing the model with the lowest BIC.
On average, the ARIMA model achieves a BIC = 1846.88on the entire data set, BIC = 534.58for the restricted data set and BIC = 670.27for the unrestricted data set.Optimal 6 GNAR models for each COVID-19 network achieve much lower BIC.For the restricted phase, the best phase-specific GNAR model yields BIC = 58.91 on the Eco- nomic hub network, and for the unrestricted phase BIC = 190.07 on the KNN ( k = 21 ) network, see Table 2.When fitted on the entire data set, the GNAR-5-11110 model with the KNN network ( k = 21) 7 achieves the lowest BIC = 193.95.8All BICs are considera- bly smaller than those obtained from the ARIMA fit, with the minimal BIC for the entire data set much larger than the BIC for the restricted phase, and also larger than the BIC for the unrestricted phase, thus justifying the use of GNAR models, as well as the split of the data.
The nature of the virus suggests that the transmission of COVID-19 between Irish counties may depend strongly on the population flow between counties (Lotfi et al. 2020).Protective COVID-19 restrictions taken by the Irish Government restricted and at times forbade inter-county travel in Level 3-5 lockdowns (Department of the Table 2 Overview over the best performing model and network for restricted and unrestricted pandemic phases average residual ε and average (av.)MASE indicated for the predicted 5 weeks at the end of the observed time period across all counties, 11.04.2021-09.05.2021 for the restricted dataset and 25.12.2022-23.01.2023  5 Fitted for each county individually, the ARIMA models might differ in orders and parameters.
6 Optimal describes the best performing combination of α -and β-order as well as global-α setting which obtain the minimal BIC value.The a-priori range of α-order spans {1, ..., 7} .The maximum lag to consider follows from Schwert's rule (Ng and Perron 1995), applied to the minimum number of weeks across the individual five datasets ( w = 18 ).The possible choices for the β-order are listed in Supplementary Material B.4.The maximum neighbourhood stage that can be included in the GNAR model is determined by the smallest maximum SPL across most networks, which is 5.For the complete network, only 1st stage neighbourhood can be modelled, while for the Economic hub network the maximum neighbourhood stage is 4. 7 From hereon referred to as the KNN network.
Taoiseach 2020; McQuinn et al. 2021).As supported by the positive and negative trends in Moran's I, the spatial dependence of COVID-19 incidence across counties is likely to have decreased during lockdowns and increased during periods in which inter-county travel was allowed (Wang 2022).This provides additional subject specific motivation to train a pandemic phase specific GNAR model.

Pandemic phases
Table 2 summarises the optimal GNAR models and COVID-19 network for the restricted and unrestricted data set.For both phases, the best performing GNAR model select an autoregressive component of order 7.The average MASE are smaller for the restricted than the unrestricted pandemic phases, implying that GNAR models are more suited to predicting periods with strict regulations than periods with fewer or no restrictions.The variance for residuals and MASE is smaller for the GNAR model than the ARIMA model.The optimal network for the unrestricted pandemic phase is much denser than the optimal network for the restricted phase.As evident from Tables 3, 4 in Supplementary MaterialD.1, the BIC values for the optimal GNAR model lie within the range [58.91, 68.36] for the restricted data set and within the range [190.07, 192.56] for the unrestricted data set.Figure 13 in Supplementary Material D.1 illustrates that denser networks perform better for the unrestricted data set while sparse networks achieve lower BIC for the restricted data set, on average.A decrease in inter-county dependence due to COVID-19 restrictions should result in decreasing values for the β-coefficients in the GNAR model.This hypothesis can only be partially verified, see Fig. 6.The absolute value of β-coefficients increases from the restricted to the unrestricted phase, implying increased spatial dependence after COVID-19 restrictions have been eased or lifted.We note that the GNAR model picks up a decrease in temporal dependence in COVID-19 ID.As a disease spreads more freely due to lenient or no restrictions, it has been observed in other data studies that case numbers can grow more erratic and become less dependent on historic data (Firth 2020;Kraemer 2020).This effect, in addition to peaks and high volatility in COVID-19 ID observed during pandemic phases with less strict regulations, might contributed to the negative α-coefficient values for the unrestricted data set.In gen- eral, an increase in noise during unrestricted pandemic phases might contribute to the decrease in predictive performance of the GNAR model.
Identical observations can be made when considering how the coefficients develop between the restricted and unrestricted phase for the GNAR model that is optimal for the entire data set, namely, GNAR(5-1,1,1,1,0) with the KNN ( k = 21 ) network; the β-coefficients increase in absolute value for the unrestricted phase compared to the restricted phase, see Fig. 14 in Supplementary Material D.2.
The predictive accuracy for both datasets is comparable and varies from county to county, see Figs. 7 and 8 for 9 example counties; MASE values for the remaining counties follow similar patterns.
For the restricted phase, GNAR models achieve lower MASE than the ARIMA models except for counties Cavan, Galway, Leitrim, Louth, Sligo, Tipperary, Roscommon, for which the ARIMA model performs equally well.For the unrestricted phase, the ARIMA model predicts particularly poorly for counties Carlow, Kilkenney, Louth and Waterford, Fig. 7 MASE values for the restricted pandemic phase, for 9 selected counties and for the restricted phase for counties Cavan, Clare, Limerick and Wexford.We note that Cavan, Leitrim and Sligo have a border with Northern Ireland, which could introduce some confounding factors.
For the restricted phase, the predicted COVID-19 ID values differ more strongly between the GNAR model and the ARIMA model, see Figure 16 in Supplementary Material D.3.For the unrestricted phase, the GNAR and ARIMA models estimate roughly the same trajectory while the former achieves smaller residuals for most counties.

Assessing the model assumptions
The above models assume that the observations follow a Gaussian i.i.d.error structure.To assess this assumption, we test whether the residuals εi,t follow a normal distribu- tion with a county-specific Kolmogorov-Smirnov test, aggregated over time.In contrast, we obtain a majority of significant p-values across counties for the unrestricted phase (the counts of p-values are #p ≤ 0.025 = 21 , #p > 0.025 = 5), 9 raising doubts about the Gaussianity assumption in the unrestricted phase.We obtain primarily insignificant p-values across counties for the restricted phase ( #p ≤ 0.025 = 8 , #p > 0.025 = 18)10 , so that the Gaussian assumption is not rejected.
Table 5 in Supplementary Material D.4 details the average MASE, average residual and p-value for each county, resulting for the two optimal GNAR models for the restricted and unrestricted data set.The Gaussian nature of residuals indicate suitability of the GNAR model to model restricted pandemic phases and ensure consistency in coefficient estimates.For the unrestricted phase, the Gaussianity in the model assumptions could not be statistically verified.These conclusions are supported by the county-specific QQplots in Supplementary Material D.4.The hyperparameters α -and β-order were set data-adaptively to minimize the BIC criterion, which assumes Gaussian error.Under the assumption of Gaussian error, the BIC would be minimal for the correct higher order network dependence.As we examined a large range of β-order choices, this deviance from Gaussianity leads us to propose investigating alternative error structures as future work (Fig. 9).
The GNAR model further assumes that the errors are uncorrelated.To assess this assumption, the residuals are investigated according to their temporal as well as spatial autocorrelation using the Ljung-Box test and Moran's I based permutation test (Ljung and Box 1978;Moran 1950).The former concludes significant temporal correlation for short-term lags in the GNAR residuals for each county.Thus, there is evidence that the GNAR model insufficiently accounts for temporal dependence in COVID-19 incidence in subsequent weeks.The residuals show remaining spatial autocorrelation.The Moran's I based permutation test counts N m = 8 Moran's I values outside the corresponding 95% credibility interval (expected 0.05 • 45 ≈ 2 ) for the restricted phase and N m = 16 for the unrestricted phase (expected 0.05 • • • 107 ≈ 5 ).The reduction in spatial correlation for the restricted phases and the Economic hub network is greater ( N m = 10 on COVID-19 cases to N m = 8 for residuals) than for the unrestricted phases and the KNN network ( N m = 16 for COVID-19 cases and residuals).We conclude that there is evidence that the GNAR model may not sufficiently incorporate the spatial relationship in COVID-19 case numbers across counties.These possible violations of the model assumptions have to be taken into account when interpreting the model fit.

Discussion
In general, a network model could can be a powerful tool to inform the spread of infectious diseases, see for example (Britton 2019;Overton 2020).In this paper, we modelled the COVID-19 incidence across the 26 counties in the Republic of Ireland by fitting GNAR models, leveraging different networks to represent spatial dependence between the counties.While we do not assume that the disease only spreads along the network, we consider the edges to represent the main trajectory of the infection.The analysis shows that there is a clear network effect, but networks of similar density perform similarly in predictive accuracy.GNAR models perform better on data collected during pandemic phases with inter-county movement restrictions than data gathered during less restricted phases.Sparse networks perform better for the restricted data set, while denser networks achieve lower BIC for the unrestricted data set.
There are some caveats relating to the model.First, the time series for the restricted phase and for the unrestricted phase are actually concatenated time series; Fig. 2 shows the original time series.The concatenation was carried out because of data availability, but it is possible that it obfuscates some potentially interesting phase-specific signals.
As seen in Fig. 2, even after differencing the time series do not display clear stationarity.Further, COVID-19 is subject to "seasonal" effects, e.g.systematic reporting delays due to weekends and winter waves (Kubiczek and Hadasik 2021;Nichols 2021;Sartor 2020).The GNAR model does not have a seasonal analogue which can incorporate seasonality in data, like SARIMA for ARIMA models (Shumway et al. 2000).Future work could introduce a seasonal component to the GNAR model, improving its applicability to infectious disease modelling.There may be other spatio-temporal patterns such as non-linear effects which the GNAR model currently does not include.Moreover, the COVID-19 pandemic had a strong influence on mobility patterns (COVID-19 Community Mobility Report 2022; Manzira et al. 2022), in particular due to restrictions of movement and an increased apprehension towards larger crowds.Considering only static networks may introduce a bias to the model (Bansal et al. 2007;Mo 2021;Perra 2012).Future work could therefore explore how GNAR models can include dynamic networks to incorporate a temporal component of spatial dependency.Alternative weighting schemes for GNAR models could be investigated to account for differences in edge relevance across time and network.
Regarding the theory of GNAR models, alternative error distributions, such as a Poisson distributed error term as in Armillotta and Fokianos (2021), could be explored given the indication of non-Gaussian residuals for the unrestricted pandemic phase.The stability of parameter estimation in GNAR models also warrants further investigation.The network constructions themselves could also be refined.Simulations have shown that GNAR models are sensitive to network misspecifications.Omitting edges may result in bias in the GNAR coefficients.While this paper has carried some robustness analysis regarding network choice, future analysis could focus on more content-based approaches to constructing networks, e.g.building a network based on the intensity of inter-county trade, computed according to the gravity equation theory (Chaney 2018).Many researchers have successfully modelled the initial spread of COVID-19 from Wuhan across China based on detailed mobility patterns, e.g.Jia (2020); Kraemer (2020).Finally, our statistical analysis did not include information about the dominant strain.
With COVID-19 being an evolving disease, different strains may display different transmission patterns.If more detailed data become available then this question would also be of interest for further investigation.However, this study has illustrated that it may be of benefit to use a GNAR model for the spread of an infectious disease, in particular during movement restrictions, so that the spread is mainly local.It has also detailed a range of possible network choices, and it has provided a set of tests to assess the performance of the model and its fit.

Fig. 1
Fig. 1 Map of Ireland and COVID-19 networks; economic hub towns marked in blue for the COVID-19 DNN network measures 90.3 km, between Kerry and Cork, and the maximal value 338.5 km, between County Cork and Donegal.The KNN and DNN network parameters are chosen to minimise the BIC of the associated GNAR model.For the restricted pandemic phase, the KNN network has k = 11 and the DNN network d = 325 .For the unrestricted pandemic phase, it is k = 21 and d = 325.

Fig. 2
Fig.2Weekly (1-lag) COVID-19 incidence difference (ID) and weekly COVID-19 incidence (dark grey dashed lines) for 100,000 inhabitants, from the start of the pandemic in March 2020 to mid January 2023, red for restricted phases, green for unrestricted; predominant COVID-19 virus variants shown by color scale at the bottom, in order: Original, Alpha, Delta, Omicron I and Omicron II

Fig. 4
Fig. 4 Moran's I across time, with weights based on SPL; main COVID-19 regulations by the Irish Government indicated by vertical lines; in order: initial lockdown, county-specific restrictions, Level-5 lockdown, allowance of inter-county travel, official end of all restrictions; 95% credibility interval in grey dashed

Fig. 6
Fig. 6 Development of GNAR model coefficients for the restricted and unrestricted pandemic phase; restricted phase with Economic hub network, unrestricted phase with KNN ( k = 21 ) network

Fig. 8
Fig. 8 MASE values for the unrestricted pandemic phase, for selected

Fig. 9
Fig.9QQ-plot for the residuals from the best performing GNAR model and network (Economic hub network for restricted phase, KNN (k= 21) network for unrestricted phase) for restricted and unrestricted pandemic phase for county Dublin (left) andDonegal (right)