Spatio-temporal Object-Oriented Bayesian Network modelling of the COVID-19 Italian outbreak data.

The spatial epidemic dynamics of COVID-19 outbreak in Italy were modelled by means of an Object-Oriented Bayesian Network in order to explore the dependence relationships, in a static and a dynamic way, among the weekly incidence rate, the intensive care units occupancy rate and that of deaths. Following an autoregressive approach, both spatial and time components have been embedded in the model by means of spatial and time lagged variables. The model could be a valid instrument to support or validate policy makers' decisions strategies.

Introduction

Pandemic due to COVID-19 disease, originated in the Chinese province of Hubei in December 2019, has been rapidly spread worldwide causing the most alarming health emergency of the last two centuries (Li et al., 2020), after the Spanish influenza that occurred in 1918.

Italy was one of the first European countries seriously hit by pandemic, at the beginning of February 2020; the two regions of Lombardy and Veneto were the first and worst affected by COVID-19 outbreak. On 9 March 2020, a national lockdown was imposed by the Italian Government to reduce the exponential growth of daily infections and, consequently, the number of hospitalized patients, mainly those in the intensive care units (Farcomeni et al., 2021), and that of deaths. Social distancing and the other restrictive measures took effect at the beginning of May 2020 when the incidence rate as well as the hospital overflow and the number of deaths strongly decreased.

The second phase of pandemic covered the summer period during which the outbreak seemed under control. The third phase came in October 2020 when Italy faced its second wave, more severe than the first one since widely spread in all regions, in the Southern ones too.

To deal with this new health emergency, the Italian Government introduced new restrictions whose degree of intensity varied among regions and were based on the 21 parameters provided by the expert committee engaged by the national authorities to monitor pandemic’s dynamic. Indeed, since November 2020, Italian regions are classified, weekly, into three (then extended to five) risk profiles. During March 2021, Italy once again faced the pandemic, being in the full third wave.

All COVID-19 data are provided, daily, by the Italian Civil Protection Department and could be freely downloaded from the GitHub repository http://www.protezionecivile.gov.it/attivita-rischi/rischio-sanitario/emergenze/coronavirus. Data are available with a greater detail at regional level; we argue that their reliability is low since they are strongly affected by mis-reporting, notification delays and by the lack of uniformity among regions in collecting and providing information, above all, in the regional planning of the daily number of swabs.

Taking into account all these issues, in this study we focused on the weekly incidence rate, i.e. the number of weekly new infections over regional population per 100 000 inhabitants, on the intensive care units saturation, by considering the weekly median number of patients in the intensive care units over the maximum number of available beds in each region1 and the weekly number of people officially died from Coronavirus over regional population per 100 000 inhabitants. Looking at the corresponding plots shown in Fig. 1, Fig. 2, Fig. 3, the outbreak dynamics with the two waves is fairly evident.

Fig. 1

Weekly new infections over regional population per 100 000 inhabitants.

Fig. 2

Weekly median number of patients in intensive care units over the maximum number of available beds.

Fig. 3

Weekly deaths due to COVID-19 disease over regional population per 100 000 inhabitants.

The aim of the work is twofold since we are interested in the spatial and time dynamics of the pandemic. To this purpose, in the last year, some interesting works have been proposed to deal with the same scope; in Dickson et al. (2020) and Giuliani et al. (2020), an endemic–epidemic time-series mixed-effects generalized linear model for areal disease counts has been applied to predict spatio-temporal diffusion of the phenomenon. In Sartorius et al. (2021), a Bayesian hierarchical space–time SEIR model has been used to assess the spatio-temporal variability of COVID-19 caseloads and deaths at small-area scale in England; they estimated the number of cases and deaths at small-area resolution as well as the impact of mobility restrictions on the COVID-19 contagion and the role of some socio-demographic risk factors on the mortality risk. Also  Bertuzzo et al. (2020) included the mobility among communities, the timing of infection seeding, the mobility restrictions and social distancing to define a spatial model of the COVID-19 spread in Italy. In Bartolucci and Farcomeni (2021), a model based on discrete latent variables, spatially associated and time specific, has been proposed for the analysis of incident cases identifying a common trend and assigning, for each week, each Italian region to one of five risk profiles.

Spassiani et al. (2021) proposed a space–time statistical analysis to deal with the issue of partitioning a spatial area (the Veneto region) into separate components of sub-units sharing a common spread pattern using mathematical morphology, hierarchical clustering, parametric fitting and non-parametric hypothesis testing. In the field of clustering, D’Urso et al. (2021) proposed a spatial robust fuzzy clustering model based on B-splines to cluster the 20 Italian regions taking into account the time and space component. Other clustering approaches have been based on Moran’s I index and the local indicators of spatial association (LISA), see Ramírez-Aldana et al. (2020) and Zhang et al. (2020). In Cuadros et al. (2020), a spatially-explicit mathematical model has been developed to simulate the transmission infection process taking into account the uneven distribution of the healthcare capacity in Ohio (U.S.). In order to predict the ICU occupancy at regional level, Farcomeni et al. (2021) combined two forecasting methods, a generalized linear mixed regression model and a region-specific time-series model for counts.

In our work too, we defined a spatial–temporal model for the Italian COVID-19 outbreak data with a particular focus on the relationships between the infection rate and the ICU occupancy rate using the theoretical formalism of Bayesian Networks (BNs). In particular the Object-Oriented Bayesian Networks (OOBN) have been proposed that embed in a unique framework both components, modelling the (spatial and time) dependence relationships between the domain variables in a multivariate context.

We believe that the added value of using OOBN in this field is really that of modelling the spatio-temporal relationships using a multivariate approach, also providing a pictorial representation of the dependence structure by means of a graph that is one of the main strengths of BNs. Moreover, its inference engine could be used to simulate scenarios, fixing time and or space or moving across them, in order to support public decisions processes.

The paper is organized as follows. In Section 2, the OOBNs are described by a methodological point of view showing all their potentialities; in Section 3, the application results to COVID-19 data are shown while in Section 4, further considerations and new research starting points are addressed.

Methodology

The Object-Oriented Bayesian Networks (Koller and Pfeffer, 1997) are the extension of the Bayesian Networks (Cowell et al., 1999), i.e. probabilistic graphical models able to handle complex dependence structures between random variables. Properly, a Bayesian network is a multivariate statistical model satisfying sets of (conditional) independence statements encoded in a Directed Acyclic Graph (DAG), consisting of a set of nodes (or vertices) and a set of links (or edges) between pairs of the nodes.

Each node corresponds to a random variable of the domain while edges between two nodes represent the probabilistic dependence among corresponding variables. An arrow from to means that is influenced by : is said to be the child of that is, in turn, said to be the parent of . All parents of a node, say , in the graph , are denoted by ; in the DAG of Fig. 4, for instance, , .

Fig. 4

DAG example.

The independence assumptions can be read off directly from the graph looking at the pairs of nodes not connected by an edge: in the DAG of Fig. 4, for example, is conditionally independent from given .

The dependencies are quantified through a set of probability tables (CPTs) associated to each node, that are conditional if the node has parents nodes, marginal otherwise. Furthermore, in a BN, the joint distribution can be factorized according to the DAG structure as: Eq. (1) implies that the joint distribution factorizes according to the DAG structure into the product of local distributions performing inference efficiently in a Bayesian network.

For discrete BNs (Heckerman et al., 1995), the assumption is that and each component follow a multinomial distribution so that: where , the parameter of the local distributions, is the conditional probability of given the th parents’ configuration.

For sake of completeness, we specify that, in the continuous and mixed case, the learning and inference methods work under the assumption of a multivariate Gaussian and conditional Gaussian distribution, respectively.

The parameters of a BN may be continuously updated as a new information becomes available and BN itself can be used to carry out what-if analysis. The BN, in fact, performs probabilistic inference by means of computational efficient algorithms allowing inference rapidly. Lauritzen and Spiegelhalter (1988) and Jensen et al. (1990) developed algorithms by which information (the so-called evidence) is propagated throughout the network thanks to a message passing algorithm defined on a tree structure (the junction tree) derived from the BN itself.

By means of these algorithms built for evidence propagation, one can update the marginal probability distributions of some variables when information on one or more other variable distributions is observed (the so-called belief updating).

It is worth noticing that classical BNs are widely used as a valid tool to deal with reasoning under uncertainty but they refer to all situations in which the dependence relationships are evaluated at a given instant, not being able to model temporal sequences or large and complex domains which involve inter-related objects. All these limitations are overcome by the OOBNs that have been already efficiently applied, for example, in forensic genetics (Dawid et al., 2007), in the context of measurement error (Marella and Vicard, 2013) and in Finance (Mortera et al., 2013, De Giuli et al., 2019).

In this application, OOBNs are used to cope with the issue of modelling spatial time series using the BN framework with its inference engine (see Wilkinson et al., 2013 for an interesting application in the same field). Properly, an OOBN allows a hierarchical definition of a BN using building blocks (classes). With respect to the static BN, the OOBN is a network that includes not only ordinary nodes but also objects, the so called instance nodes: they are subnets representing instances of other network classes and they can themselves contain instance nodes.

Objects are connected to other nodes through some of its basic nodes, the so called interface nodes, that can be distinguished into input nodes and output nodes, whose definition will be clarified in the next paragraph.

In particular, when applied to model temporal sequences, each OOBN’s instance node corresponds to a time slice whose input and output nodes play a key role in modelling the temporal relationship between variables.

To get reader familiar with all these definitions, we provided a detailed description of them with reference to our proposed spatio-temporal model in the next section.

The spatio-temporal OOBN modelling Italian COVID-19 data

The OOBNs have been efficiently applied to model COVID-19 time series of the variables listed in Table 1, at a regional level. They span from the last week of February 2020 to the first week of February 2021.

Table 1

Variables description.

Label	Description	Categories
Area	Geographic areas	N-W (North-West); N-E (North-East); C (Centre); S-I (South-Isles)
Deaths	Weekly deaths due to COVID-19 disease over regional population per 100000 inhabitants	0−1; 1−10; ≥10
ICU	Weekly median number of patients in intensive care units over the maximum number of available beds in each region	0−0.1; 0.1−0.3; 0.3−1
newP	Weekly new infections over regional population per 100000 inhabitants	0−20; 20−100; 100−200; ≥200
newP_N	The mean of weekly new infections in the neighbouring areas over regional population per 100000 inhabitants	0−20; 20−100; 100−200; ≥200

Based on an autoregressive approach, the proposed model takes into account three time slices, i.e. three consecutive weeks. Each time slice is defined by a static BN, shown in Fig. 5, whose dependence structure has been built taking into account the logical and temporal relationships among variables. As already mentioned, the interface nodes have to be declared; therefore, the so-called input nodes and output nodes have been identified. The former have been used to model the temporal sequence and they should not be intended as real nodes but as parameters, the latter are instead real nodes that should be bound to the input nodes of the next time slice. To this purpose, the nodes newP_prev, newP_N_prev, ICU_prev and Area, denoted by dotted and dashed outlines, have been declared as input nodes: being not real nodes but placeholder nodes, the first three have been connected by a directed link to the real nodes newP, newP_N and ICU in the previous time slice respectively, and the fourth one to the real node also called Area, put outside the time slices and visible at the higher level of the hierarchy (see Fig. 6).

Fig. 5

The static BN in each time-slice.

Fig. 6

The spatio-temporal OOBN for the Italian COVID-19 data.

In turn, the nodes newP, newP_N and ICU, denoted by solid and shaded outlines, have been declared as output nodes and they play a key role since, through them, the information can be propagated from a time slice to another one. Both input nodes and output nodes are visible outside each class differently from the so-called private node, that are used to identify all nodes not visible outside the class. In our model the only declared private node has been Deaths, visible only within the current time slice.

As already said, at the higher level, the interface nodes are the only ones included and visible in the instance nodes, represented by round-shaped rectangles, corresponding, each one, to a time slice, as shown in Fig. 6.2 The top level network is, then, shown in Fig. 7.

Fig. 7

Top level network.

We argue that, while the temporal component is modelled by means of the time slices and of the probabilities , and , the spatial component is embedded in the model by including the node newP_N, that is the variable newP lagged by means of a spatial weight matrix ,3 whose element is defined as follows: where is the number of neighbours of polygon .

Since each element of th row is divided by the number of its neighbours (the row normalization of ), the effect of any individual neighbour decreases as the number of neighbours increases. We argue that also the node Area, accounting for the geographical stratification, plays a key role in explaining similarities and differences among territories.

Obviously, that based on boundaries is one of the two main approaches used to define spatial connectivity; the other one is based on suitable functions of the distance between two spatial units.

To this purpose, it is worth noticing a recent proposal suggesting the use of alternative time-dependent spatial linkage structures like the number of flight connections and the relationships in international trade (Krisztin et al., 2020).

We specify that all CPTs have been estimated directly from data by means of the EM algorithm implemented in the Hugin software under the assumption that the data structure and the tables of the time slices are identical, therefore each transition probability is identical for all time slices (which is often the case in real-world applications).

As we can see in the next section, the inference process too, used to simulate scenarios of interest, is particularly efficient nonetheless the complexity of the structure.

Results

In this section, we used the estimated OOBN for inference, thus simulating some scenarios of interest by inserting and propagating the available information (called evidence) on some nodes, through the network, to update the marginal probabilities of all the remaining nodes. We argue that an OOBN has a structure comparable to that of a BN, therefore it is possible to use the same inference algorithms for belief updating proposed for classical BNs, already briefly cited in Section 2.

In Table 2, the marginal probabilities before inserting evidences have been reported too.

Table 2

Marginal probability distribution for each node of the OOBN and each time slice.

Week	Area				ICU			Deaths			newP				newP_N
	N-W	N-E	C	S-I	0-0.1	0.1–0.3	0.3–1	0–1	1–10	>= 10.0	0–20	20–100	100–200	>= 200	0–20	20–100	100–200	>= 200
i-3	0.200	0.200	0.200	0.400	0.590	0.249	0.161				0.528	0.196	0.126	0.150	0.508	0.192	0.144	0.156
i-2	0.200	0.200	0.200	0.400	0.560	0.322	0.119	0.472	0.432	0.096	0.512	0.201	0.134	0.152	0.435	0.255	0.169	0.142
i-1	0.200	0.200	0.200	0.400	0.563	0.312	0.125	0.501	0.422	0.077	0.498	0.204	0.143	0.155	0.434	0.244	0.180	0.143
i	0.200	0.200	0.200	0.400	0.561	0.301	0.139	0.509	0.412	0.079	0.486	0.206	0.150	0.158	0.439	0.230	0.183	0.148

The first scenario simulates a typical emergency situation in which, for three consecutive weeks, the number of weekly new infections per 100 000 inhabitants has been greater than or equal to 200.

As we can see from Table 3, the probability of exceeding the critical threshold of ICU bed occupancy rate, fixed by the Italian authorities to 30%, considerably increased, week by week, with the risk, if any mitigation measures will be imposed, to have the 53.2% of probability to observe a heavy overflow of intensive care units at the week . The corresponding updated marginal probabilities of the node newP_N, give evidence of a strong positive spatial dependence in the transmission process.

Table 3

First scenario: marginal probabilities.

Week	Area				ICU			Deaths			newP				newP_N
	N-W	N-E	C	S-I	0-0.1	0.1–0.3	0.3–1	0–1	1–10	>= 10.0	0–20	20–100	100–200	>= 200	0–20	20–100	100–200	>= 200
i-3	0.200	0.200	0.200	0.400	0.590	0.249	0.161				0.000	0.000	0.000	1.000	0.508	0.192	0.144	0.156
i-2	0.195	0.239	0.178	0.393	0.091	0.651	0.258	0.012	0.669	0.319	0.000	0.000	0.000	1.000	0.158	0.161	0.306	0.375
i-1	0.190	0.239	0.178	0.393	0.040	0.538	0.422	0.004	0.671	0.325	0.000	0.000	0.000	1.000	0.069	0.076	0.294	0.561
i	0.190	0.239	0.178	0.393	0.020	0.448	0.532	0.007	0.636	0.356	0.002	0.007	0.138	0.853	0.037	0.046	0.291	0.626

To highlight the differences among North, Centre and South of Italy, the same scenario has been replicated for each geographic area, as reported in Table 4, Table 5, Table 6, Table 7. As we can observe, the South and Isles had a lower increase of the probability of heavy ICU overflow than the other areas, characterized by a sharp growth of the number of patients admitted to intensive care units.

Table 4

First scenario with reference to North-West of Italy: marginal probabilities.

Week	Area				ICU			Deaths			newP				newP_N
	N-W	N-E	C	S-I	0-0.1	0.1–0.3	0.3–1	0–1	1–10	>= 10.0	0–20	20–100	100–200	>= 200	0–20	20–100	100–200	>= 200
i-3					0.590	0.249	0.161				0.000	0.000	0.000	1.000	0.508	0.192	0.144	0.156
i-2	1.000	0.000	0.000	0.000	0.295	0.401	0.304	0.000	0.753	0.247	0.000	0.000	0.000	1.000	0.175	0.175	0.175	0.475
i-1	1.000	0.000	0.000	0.000	0.147	0.319	0.533	0.000	0.685	0.315	0.000	0.000	0.000	1.000	0.088	0.088	0.088	0.737
i	1.000	0.000	0.000	0.000	0.069	0.206	0.725	0.008	0.598	0.394	0.000	0.000	0.172	0.828	0.044	0.059	0.107	0.790

Table 5

First scenario with reference to North-East of Italy: marginal probabilities.

Week	Area				ICU			Deaths			newP				newP_N
	N-W	N-E	C	S-I	0-0.1	0.1–0.3	0.3–1	0–1	1–10	>= 10.0	0–20	20–100	100–200	>= 200	0–20	20–100	100–200	>= 200
i-3					0.590	0.249	0.161				0.000	0.000	0.000	1.000	0.508	0.192	0.144	0.156
i-2	0.000	1.000	0.000	0.000	0.000	0.754	0.246	0.000	0.651	0.349	0.000	0.000	0.000	1.000	0.175	0.190	0.277	0.358
i-1	0.000	1.000	0.000	0.000	0.000	0.492	0.508	0.000	0.680	0.320	0.000	0.000	0.000	1.000	0.091	0.120	0.294	0.494
i	0.000	1.000	0.000	0.000	0.012	0.333	0.655	0.008	0.593	0.399	0.006	0.006	0.058	0.929	0.057	0.084	0.285	0.574

Table 6

First scenario with reference to Centre of Italy: marginal probabilities.

Week	Area				ICU			Deaths			newP				newP_N
	N-W	N-E	C	S-I	0-0.1	0.1–0.3	0.3–1	0–1	1–10	>= 10.0	0–20	20–100	100–200	>= 200	0–20	20–100	100–200	>= 200
i-3					0.590	0.249	0.161				0.000	0.000	0.000	1.000	0.508	0.192	0.144	0.156
i-2	0.000	0.000	1.000	0.000	0.197	0.369	0.435	0.066	0.650	0.285	0.000	0.000	0.000	1.000	0.175	0.175	0.297	0.353
i-1	0.000	0.000	1.000	0.000	0.066	0.321	0.614	0.022	0.592	0.386	0.000	0.000	0.000	1.000	0.088	0.088	0.342	0.483
i	0.000	0.000	1.000	0.000	0.023	0.201	0.775	0.019	0.563	0.419	0.000	0.012	0.185	0.802	0.049	0.045	0.406	0.499

Table 7

First scenario with reference to South Italy and the Islands: marginal probabilities.

Week	Area				ICU			Deaths			newP				newP_N
	N-W	N-E	C	S-I	0-0.1	0.1–0.3	0.3–1	0–1	1–10	>= 10.0	0–20	20–100	100–200	>= 200	0–20	20–100	100–200	>= 200
i-3					0.590	0.249	0.161				0.000	0.000	0.000	1.000	0.508	0.192	0.144	0.156
i-2	0.000	0.000	0.000	1.000	0.000	0.837	0.163	0.000	0.648	0.352	0.000	0.000	0.000	1.000	0.131	0.131	0.392	0.347
i-1	0.000	0.000	0.000	1.000	0.000	0.770	0.230	0.000	0.695	0.305	0.000	0.000	0.000	1.000	0.038	0.038	0.372	0.552
i	0.000	0.000	0.000	1.000	0.001	0.747	0.252	0.001	0.715	0.284	0.000	0.008	0.150	0.842	0.016	0.016	0.332	0.636

It is worth noting that all areas showed a clear positive spatial dependence in terms of incidence rate, the North-East in particular.

Next, we simulated a typical scenario of lockdown with a constant decrease of new infections across three consecutive weeks (from the upper class to the two lower ones).

As it was foreseeable, looking at the probabilities shown in Table 8, the ICU occupancy rate decreased but slowly, as the number of deaths. If we look at the node NewP_N, the presence of a positive spatial effect is evident also in this case. It is worth noting that the positive effect of a lockdown on ICU beds saturation becomes visible the week after, as expected.

Table 8

Second Scenario: marginal probabilities.

Week	Area				ICU			Deaths			newP				newP_N
	N-W	N-E	C	S-I	0-0.1	0.1–0.3	0.3–1	0–1	1–10	>= 10.0	0–20	20–100	100–200	>= 200	0–20	20–100	100–200	>= 200
i-3	0.200	0.200	0.200	0.400	0.590	0.249	0.161				0.000	0.000	0.000	1.000	0.508	0.192	0.144	0.156
i-2	0.177	0.133	0.143	0.547	0.208	0.402	0.390	0.200	0.514	0.286	0.000	0.000	1.000	0.000	0.173	0.207	0.371	0.249
i-1	0.177	0.133	0.143	0.547	0.276	0.467	0.257	0.067	0.767	0.167	0.000	1.000	0.000	0.000	0.163	0.312	0.380	0.145
i	0.177	0.133	0.143	0.547	0.281	0.499	0.220	0.291	0.604	0.105	0.127	0.696	0.173	0.004	0.199	0.372	0.360	0.069

In Table 9, another hypothetical lockdown scenario has been simulated, with a drastic reduction of contagion and, then of the ICU occupancy. Indeed, at the week , the probability to be in the last class of ICU is less than 4% while the probability to record a weekly incidence of 0–20 new infections over 100 000 inhabitants is of 92% (the latter directly depends on the infection rate of the previous week ): a further strong evidence of the positive effect of the infection control policy.

Table 9

Fourth scenario.

Week	Area				ICU			Deaths			newP				newP_N
	N-W	N-E	C	S-I	0-0.1	0.1–0.3	0.3–1	0–1	1–10	>= 10.0	0–20	20–100	100–200	>= 200	0–20	20–100	100–200	>= 200
i-3	0.200	0.200	0.200	0.400	0.590	0.249	0.161				0.000	0.000	0.000	1.000	0.508	0.192	0.144	0.156
i-2	0.000	0.207	0.296	0.497	0.265	0.422	0.313	0.368	0.399	0.233	0.000	1.000	0.000	0.000	0.304	0.232	0.232	0.232
i-1	0.000	0.207	0.296	0.497	0.320	0.591	0.088	0.332	0.609	0.059	1.000	0.000	0.000	0.000	0.569	0.230	0.101	0.101
i	0.000	0.207	0.296	0.497	0.546	0.421	0.033	0.542	0.427	0.031	0.922	0.078	0.000	0.000	0.712	0.187	0.051	0.051

The simulated scenarios have been used to show the main effects of the lockdowns’ strategy and to give evidence of a strong positive spatial dependence among territories. The model could be improved embedding other variables of interest and other temporal relationships. It is worth noting that the model’s estimates could be updated when new data becomes available, enhancing their reliabilities.

Conclusions

The aim of the work has been that of modelling the interactions among the main factors related to the COVID-19 outbreak in Italy taking into account both temporal and spatial information.

The methodological framework of the Object-Oriented Bayesian networks has been used to define a dynamic model in a multivariate context. By simulating scenarios, all the OOBNs potentialities have been explored proving that they could be a valid instrument to support decision making processes and strategies. Moreover, other model specifications could be defined in the future. Indeed, we believe that a more interesting application could regard the study of the epidemic dynamic taking into account the parameters chosen by the Scientific Committee to define the regional risk profiles. This could be useful to evaluate the efficacy of the strategy based on local lockdowns and different levels of social restrictions.

All this could further increase the informational gain and, therefore, the knowledge useful for decision-making and government actions to deal with the COVID-19 pandemic emergency even more consciously and then for strategies capable to mitigate and neutralize the negative effects of the COVID-19 outbreak in Italy.