Probabilistic Graphical Models

The knowledge of probability theory provides the framework for assessing multiple variables and their likelihood; this makes the modeling of complex reasoning more faithful to reality. The reasoning task is to estimate the probability of one or more variables given the observation of other variables. The joint probability distribution over the space of possible states (where each variable has one or multiple states) for the set of random variables \(\mathcal{V}\) need to be constructed. However, computing the joint distribution over a set of dozens or hundreds of random variables is a daunting task, but PGMs can be used to simplify the construction of complex distributions. PGMs as the name implies use a graph-based representation to encode the conditional (in)dependencies between random variables.The graph depicts the causal relationships that exist between variables in the distribution. Consider a distribution modeled as a graph \(\mathcal{G}(\mathcal{V}, \mathcal{E})\) where \(\mathcal{V}\) is the set of variables and \(\mathcal{E}\) is the set of edges; the variable interactions can be expressed in the form \((A \perp B \mid C)\), where \(A, B, C \in \mathcal{V}\) and \((C,A) \in \mathcal{E}\).

The joint probability distribution over a set of variables \(\mathcal{V} = \{V_1, …, V_n\}\) can be expressed as \(P(V_1 = v^k_1, …, V_n = v^k_n)\) where \(v^k_i \in \mathcal{S}_i\) and \(\mathcal{S}_i\) is the set of states for variable \(V_i\). With the graph structure, the joint distribution is decomposed into smaller factors based on the conditional independencies among variables. The joint distribution can be defined as a product of the factors as:

\[\begin{equation} P(V_1, …, V_n) = \prod_{i=1}^{n}P(V_i \mid Pa(\hat{V_i})) \end{equation}\] where \(n\) is the total number of variables in \(\mathcal{G}\), while \(V_i\) and \(Pa(\hat{V_i})\) represent a variable in \(\mathcal{V}\) and the ‘parents’ of the variable respectively. This is also referred to as the conditional probability distribution of the variable set.

Bayesian Networks

To compute the probability distributions from the graph structure, there are two main approaches. One is the Bayesian networks in which case the edges are directed as parent \(\rightarrow\) child, and the other is Markov networks in which case the edges are undirected. Both networks differ in the set of independencies they encode and in the induced factorization from the distribution. Figure 1 shows a three-variable Bayesian network example with the independence assertion that \((A \perp B \mid C)\), where \(A\) is a child of \(C\) and is independent of \(B\). This can be interpreted as the Markov condition (Tsagris 2019) which states that each variable is independent of its non-descendants given its parent(s). In this example, \(B\) is not a descendant of \(A\), and the edge \(C \rightarrow A\) exists if \(C\) is a direct cause of \(A\). Probabilistically, the conditional independence holds if \(P(A,B \mid C)=P(A \mid C)P(B \mid C)\).

Formally, a Bayesian network or BN for short (Pearl 1988), is a directed acyclic graph (DAG) \(\mathcal{G}\) over variables \(V_i \in \mathcal{V}\) and a joint probability distribution \(P\). The directed graph is acyclic because no cycles are allowed. From a BN, any probability of the form \(P(V_i \mid Pa(\hat{V_i}))\) can be determined, i.e the conditional probability on each variable can be computed. Note that if a variable has no parent, only the marginal probability \(P(V_i)\) is computed. In practice, a BN is constructed from a dataset of samples, and once a feasible structure is found, we can run inference on the dataset.

Fig. 1 An example of a Bayesian network structure illustrating conditional independence

Constructing a BN requires learning the conditional independencies in the dataset. There are several structure learning algorithms which search the best graphical structure from the dataset. There are two broad structure learning techniques – score-based structure learning and constraint-based structure learning. For both techniques, a learning algorithm traverses the search space of possible DAGs and selects the one with the optimal score. The most popular ones of these algorithms include the Exhaustive search, HillClimb search, PC algorithm (Tsagris 2019) and MMHC algorithm (Tsamardinos, Brown, and Aliferis 2006). When a DAG has been found, the important task is to estimate the conditional probability distributions (CPDs), defined as the set of parameter values \(\theta_\mathcal{G}\) for the variables \(V_i \in \mathcal{V}\). For each variable \(V_i\), there is a separate multinomial probability distribution over every possible parent’s state \(v^k_i\). Thus, each variable’s CPD is considered a full table.

Estimating the CPD for each variable naturally follows the use of relative frequencies with respect to the occurrence of the variable states; this approach is called maximum likelihood estimation (MLE). The MLE has the problem of overfitting the data when there are not enough observations in the data. In contrast, the Bayesian parameter estimator starts with some pseudo state counts before observing the data - this reduces the tendency of overfitting or biased estimation.

Data Clustering

Clustering of data is important in knowing the hidden structure of the data, because it is often a complex task to superficially identify similarities between data points in a high-dimensional space and in a large variance data distribution. The goal of clustering is to identify subgroups in the data in a manner that data points that have similar properties or patterns belong to the same cluster, while data points that are not similar belong to different clusters. Clustering algorithms can help find such similarities in the data by using the Euclidean-distance function to determine the assignment of data points to appropriate clusters. Unlike supervised learning, clustering is an unsupervised learning task since no one knows the ground truth to compare and evaluate with the outcomes from clustering (Dabbura 2022).

The most popular and widely used clustering algorithm is the Kmeans algorithm due to its simplicity and robustness in fitting data points. As the name suggests, the Kmeans algorithm partitions the dataset into \(K\) pre-defined clusters where each data point belongs to only one cluster. The clustering algorithm tries to minimize an objective function which is the sum of squared error (SSE) between the data points and the cluster centers as follows:

\[\begin{equation} \underset{\mu_1,\cdots,\mu_K \in Z}{\text{minimize}} \quad \sum_{k=1}^K \sum_{z \in Z_k} \lVert z - \mu_k \rVert^2 \end{equation}\] where \(Z_k\) is the set of data points that are closest to center point \(\mu_k\) more than other data points.

Specifically, the algorithm tries to minimize the intra-cluster distance and maximize the inter-cluster distance, while keeping the SSE as low as possible. The Kmeans algorithm works by initializing \(K\) centroids (i.e cluster centers) randomly and assigning data points to the centroids, this continues to iterate until there is no change in the assignment of data points to centroids, in which case, the algorithm converges. A drawback with this algorithm is that, since the initialization of centroids is random, the algorithm may be stuck in a local minimum and not reach the global minimum, as different initializations usually lead to different clusters. This problem has been solved with the Kmeans++ algorithm (Arthur and Vassilvitskii 2006) which works similarly as the Kmeans algorithm but differs in the centroid initialization process, i.e centroids are initialized smartly, and the quality of clustering is improved.

Time-series Forecasting

The significance of learning from data helps in decision-making, and PGMs can be used to answer many questions of uncertainty from data. However, we need to be able to make predictions about the future from historical data. Machine learning is a powerful tool for such tasks as it helps in forecasting, predicting, and extracting meaningful information from historical data. In several data science fields, time-series forecasting has been considered a pertinent regression task as it pertains to the sequence of observations in time intervals, be it hourly, daily, weekly, monthly, or yearly. The regression task in time-series forecasting is to predict the future values based on current and previous observations. While this may seem to be an interesting, and perhaps simple task, it is by no means trivial, as it requires an understanding of the components of the time-series data such as trends, seasonality, irregularity, and cyclicity (“Using Machine Learning for Time Series Forecasting Project,” n.d.). These components relate to the pattern repetition of cycles (seasonality), increase or decrease in the behavior of the observations (trend) and repetitive changes in the time-series (cyclicity).

Several traditional approaches such as exponential smoothing or ARIMA have been applied to forecasting, but ML approaches have significantly achieved higher accuracy for the prediction task. With ML, the input is formulated in a supervised learning manner, where the input can be designed using time-delay embedding techniques, i.e for an instance \(i\), the feature vector is \(\mathbf{x^{(i)}}=(y^{(i-1)}, …, y^{(i-N)})\) where \(y^{(i-N)}\) is the output of the \(Nth\)-step instance before the \(i\)th instance. The goal of learning is to find the best model that fits the data such that \(\hat{y}^{(i)} = y^{(i)}\), where \(\hat{y}^{(i)}\) and \(y^{(i)}\) are the predicted value and target value respectively.

The loss function for the regression model is defined as the root mean square error (RMSE), which is the standard deviation of the prediction errors. It is interpreted as, the lower the better:

\[\begin{equation} RMSE = \sqrt{\frac{1}{M} \sum_{i=1}^{M}(\hat{y}^{(i)}-y^{(i)})^2} \end{equation}\] where \(M\) is the sample size of the dataset.

Among the popular regression models with strong capability for time-series forecasting are the random forest regression algorithm and linear regression model. The performance comparison for both algorithms on the COVID-19 time-series dataset is presented in the results and discussion section.

Probabilistic Graph Structure Learning on CDC Dataset

The scope of this work is limited to COVID-19 cases in Chicago, the largest city in Cook County, Illinois. Two separate datasets are used for the time-series forecasting and BN with data clustering. For the forecasting, a publicly available COVID-19 dataset in Google Health COVID-19 Open Data Repository (“Covid-19 Open Data - Google Health,” n.d.) is used (the dataset consists of the record of the number of daily confirmed cases and death cases from January 25, 2020, to May 13, 2022), while for the BN and data clustering, the CDC COVID-19 Case Surveillance Data is used (“Covid-19 Case Surveillance Restricted Access Detailed Data” 2021). Note that the CDC does not take responsibility for the scientific validity or accuracy of methodology, results, or statistical analysis (“Covid-19 Case Surveillance Restricted Access Detailed Data” 2021).

The CDC dataset consists of features that are sub-classified into symptoms, demographic information, and severity, for patients who are infected by COVID-19. There are a total of 24 features in the dataset and only 7,517 samples with these features are used for the Chicago case study. Table 2 shows the alphabetic coding of each variable, and their cardinality in the dataset (where cardinality refers to the number of variable states). The final six variables without codes in Table 2 are removed from the processed dataset because their values are either ‘missing’ or ‘unknown’ for most of the cases. Also, since there is a strong correlation between whether a patient is hospitalized and whether the patient is taken to the intensive care unit (ICU), to avoid redundant information, the variable ‘hospitalized’ is removed (this helps the BN to find a good structure). This is because a patient who is in the ICU is also hospitalized, and for many instances in the dataset, patients who are hospitalized are taken to the ICU. Hence, a total of 17 variables are used in constructing the BN. All of the variables except age, sex, and race, have two binary states – yes(1) and no(0). These binary values (alongside the variable codes) will be used to interpret the CPD of the variables in the results and discussion section.

To learn the BN from the dataset requires the selection of a search strategy, i.e an algorithm that finds the best structure that encodes the conditional independencies in the data. The search space of DAGs is always super-exponential with respect to the number of variables in the dataset. The exhaustive search algorithm always attempts to identify the ideal structure, but only if the dataset consists of no more than five variables (i.e nodes). Since this is not realistic in practice, as we need to handle data of high-dimensional variable space, a heuristic search algorithm helps to find a good structure that implements a greedy local search and terminates when a local maximum is reached. Figure 2 shows the constructed BN for the CDC COVID-19 dataset using the heuristic Hill-Climb search algorithm. The algorithm returns a DAG with 17 nodes and 47 edges.

Fig. 2 Directed acyclic graph using Hill-Climb search algorithm

The method used for learning the parameters, i.e estimating the CPDs for each variable in the dataset, is the Bayesian parameter estimation method. The MLE which uses the relative frequencies has the problem of overfitting the data and is considered not to benefit the objective of optimal CPDs. The following lists show the BN’s output of the causal relationship between each variable and its parents, as well as the arrow pointing to the local independencies associated with that variable, if any. The value in front of each variable represents the number of states. For example, A:2 means the variable \(A\) has two possible states.

P(A:2)

P(G:2 | A:2)

P(M:2 | A:2, E:2, G:2, K:2) \(\rightarrow\) (M ⟂ Q, N, B, H, L | A, G, K, E)

P(B:2 | A:2, G:2, L:2)

P(N:2 | B:2, H:2, K:2, L:2, Q:2) \(\rightarrow\) (N ⟂ A, R, G, D, E, M | K, Q, B, H, L)

P(Q:2 | B:2, G:2, K:2) \(\rightarrow\) (Q ⟂ A, L | K, G, B)

P(H:2 | B:2, E:2, K:2, L:2) \(\rightarrow\) (H ⟂ A, R, G, Q, M | K, B, E, L)

P(E:2 | B:2, K:2, Q:2) \(\rightarrow\) (E ⟂ A, G, L | K, B, Q)

P(K:2 | B:2, G:2, L:2) \(\rightarrow\) (K ⟂ A | G, B, L)

P(C:9 | I:2, K:2, Q:2) \(\rightarrow\) (C ⟂ R, G, B, E, L, M, A, D, F, N, H | I, K, Q)

P(O:6 | C:9) \(\rightarrow\) (O ⟂ R, G, Q, B, E, L, M, A, D, K, F, N, I, H | C)

P(P:2 | C:9, O:6) \(\rightarrow\) (P ⟂ R, G, Q, B, E, L, M, A, D, K, F, N, I, H | C, O)

P(D:2 | G:2, H:2, M:2, R:2) \(\rightarrow\) (D ⟂ A, K, Q, N, B, E, L | R, G, H, M)

P(I:2 | D:2, G:2, N:2, R:2) \(\rightarrow\) (I ⟂ Q, B, E, L, M, A, K, H | R, G, N, D)

P(R:2 | E:2, M:2, Q:2) \(\rightarrow\) (R ⟂ A, G, K, N, B, H, L | E, Q, M)

P(L:2 | G:2) \(\rightarrow\) (L ⟂ A | G)

P(F:2 | I:2, K:2, N:2) \(\rightarrow\) (F ⟂ R, G, Q, B, E, C, L, M, A, D, P, H, O | I, K, N)

The causal relationships between variables can guide us in knowing how to run inference using the query (child) variable and evidence variables (parents). For example, from the lists, we can observe that abdominal pain (A), diarrhea (G), and abnormal chest X-ray (B) have no local independencies, and that abdominal pain (A) is not influenced by any variable, but it causes diarrhea (G), nausea vomit (M) and abnormal chest X-ray (B).

Kmeans++ Clustering on CDC Dataset

Further to the probabilistic graph structure learning, there is the need to extract meaningful knowledge such as symptoms pattern from the dataset. Running the Kmeans++ algorithm requires that the number of \(K\) classes be specified. However, the major challenge is knowing what \(K\) value will give the best clustering result, as choosing the wrong \(K\) value can cause the data points to be misclassified. In designing the Kmeans++ algorithm, the elbow method (Kodinariya and Makwana 2013) is used in the experiments to serve as a guide in the selection of the possible \(K\) value for optimal clustering.

Considering the context of the clustering, which is to find symptoms patterns in the dataset, the \(K\) value suggested by the elbow method is used for the initial clustering, and then an optimal \(\hat{K}\) value is computed as:

\[\begin{equation} \hat{K} = \sum_{i=1}^{K} n_i(P)n_i(O)n_i(C) \end{equation}\]

where \(n_i(P)\), \(n_i(O)\), and \(n_i(C)\) are the cardinality of the demographic variables - sex, race, and age respectively, in class \(i\).

The clustering algorithm outputs the data points and their associated class labels. As part of the symptoms pattern learning process, the labeled data is used to train a deep neural network model in a supervised learning manner. The deep neural network architecture with an input layer, consists of four hidden layers and an output layer with ReLU activation and SoftMax activation respectively. The configurations of the hyperparameters for training the deep neural network model are listed in Table 3. The dimension m by n of the weight matrix \(\textbf{W}_i\) correspond to number of neurons by number of input features.

Regression Models for COVID-19 Time-series Dataset

The first step in designing a regression model for time-series forecasting is to know the properties of the data, e.g whether it is stationary or not. This understanding provides proper guidance on feature preprocessing steps and the parameter settings for a regression model. The COVID-19 time-series dataset, which captures the trend of daily confirmed cases and death cases is shown in Figure 3. There is a trend in the daily death cases, as more death cases are recorded during the spring and winter season than in the summer season, however, there is no such trend in the daily confirmed cases.

Fig. 3 Left-Right: Trend of daily confirmed cases and daily death cases

With the trend in the time-series data in Figure 3, the random forest, and linear regression models are designed using the time-delay embedding technique. Each input instance in the dataset is a 7-dimensional feature vector which corresponds to the number of cases recorded in the past seven days. This means that the next day depends on the cases seven days prior. Here, the random forest and linear regression models are trained separately using two dataset configurations, i.e date+confirmed cases; and date+deceased cases.

Conditional Probability Distributions

The output of learning the parameters of the BN is the CPDs for each variable. The Markov condition: “a variable is independent of non-descendants given parents” is used to derive the CPDs. Figure 4 – Figure 13 show the CPDs of each variable, where the row corresponds to the target variable states {yes(1), no(0)} and the columns correspond to the different parents’ state {yes(1), no(0)} configurations. For example, [row1, column1] in Figure 4 is interpreted P(G=0 | A=0) = P(diarrhea=’No’ | abdominal pain=’No’) = 0.90. Variables with more than two states or whose state are not yes(1) and no(0) are clearly defined. The variable ‘age group’ has nine states and it is defined as {(0-9yrs), (10-19yrs), (20-29yrs), (30-39yrs), (40-49yrs), (50-59yrs), (60-69yrs), (70-79yrs), (80+yrs)} = {0, 1, 2, 3, 4, 5, 6, 7, 8}, the variable ‘sex’ has two states and it is defined as {‘Female’, ‘Male’} = {0, 1}, and the variable ‘race and ethnicity’ has six states which is defined as:

{(American Indian/Alaska Native, Non-Hispanic), (Asian Non-Hispanic), (Black Non-Hispanic), (Hispanic/Latino), (Native Hawaiian/Other Pacific Islander), (White Non-Hispanic)} = {0, 1, 2, 3, 4, 5}

Fig. 4 CPD table for diarrhea (G)

The CPD table in Figure 4 shows that a patient is less likely to have diarrhea (G) if there is no diagnosis of abdominal pain (A) in the patient and vice-versa.

Fig. 5 CPD table for abnormal chest X-ray (B)

Whether a patient suffers from abnormal chest X-ray (B) or not depends on abdominal pain (A), diarrhea (G) and use of mechanical ventilation (L) as shown in Figure 5. There is a high probability (up to 0.94) that a patient has abnormal chest pain if the patient is assisted with the mechanical ventilation even with or without the presence of diarrhea and abdominal pain.

Fig. 6 CPD table for age group (C)

As shown in Figure 6, we can observe the CPD of the age group (C). For example, patients who are diagnosed of COVID-19 and fall within the age ranges (0-9 years), (10-19 years), (20-29 years) and (30-39 years) are less affected by headache (I), use of ICU (K) and shortness of breadth (Q), while the older population have a high probability of suffering from COVID-19 given that they are diagnosed of headache and/or are taken to the ICU for treatment.

Fig. 7 CPD table for cough (E)

A patient is most likely to have cough (E) if diagnosed of abnormal chest X-ray (B), shortness of breadth (Q) and already taken to the ICU (K). It is possible that if the patient is treated of these symptoms, the cough symptom will either be reduced or prevented. For example, as shown in Figure 7, the probability of having cough is as high as 0.86 when these parent symptoms are present in a patient, and in their absence the cough probability can be as low as 0.34.

Fig. 8 CPD table for death (F)

As shown in Figure 8, whether a patient will die of COVID-19 or not depends on headache (I), ICU (K) and pneumonia (N). The diagnosis of headache or use of ICU or both in a patient can cause a death probability of up to 0.61, 0.98 and 0.75 respectively. This means that a patient taken to the ICU is likely to die with a probability of 0.98, while a COVID-19 patient that has headache has a death probability of 0.61. The extreme case is when a patient has pneumonia, the death probability is as high as 0.96 and 0.99. This reflects the severity of these symptoms in causing COVID-19 death and can help the health offices such as clinics, medical laboratories and hospitals to administer emergency treatment to a patient with headache and/or pneumonia during the COVID-19 pandemic. Individuals themselves can use these symptoms as an indicator to seek immediate medical help.

Fig. 9 CPD table for intensive care unit (K)

The ICU (K) variable depends on abdominal chest X-ray (B), diarrhea (G) and mechanical ventilation (L). The CPD table in Figure 9 reveals that there is a high probability (up to 0.9) for the use of an ICU if a patient can’t breathe properly and needs the mechanical ventilation to breathe. The probability increases to 0.99 if the patient is diagnosed of abdominal chest X-ray and diarrhea.

Fig. 10 CPD table for mechanical ventilation (L)

Figure 10 supports the CPD table in Figure 9. As seen, the use of the mechanical ventilation (L) does not strongly depend on the diagnosis of diarrhea (G) in a patient.

Fig. 11 CPD table for race and ethnicity (O)

The outcome of the BN structure shown by the CPD table in Figure 11 does not accurately reveal a strong relationship between a patient’s race and ethnicity (O) and the age group (C). This is expected because there is no such rule or process to detect a person’s race from their age. We expect that these variables should be independent of each other, but the BN has learned the parameters from the ancestors of the parent variables.

Fig. 12 CPD table for shortness of breadth (Q)

Whether a patient has shortness of breadth (Q) or not depends on if the patient has abnormal chest X-ray (B), diarrhea (G) and is taken to the ICU (K). From Figure 12, the presence of abnormal chest X-ray greatly causes shortness of breadth in a patient. When it is present, the probability is always high compared to when it is low. In combination with the remaining two parent variables, the probability that a patient has shortness of breadth is as high as 0.94.

Fig. 13 CPD table for sore throat (R)

The CPD table in Figure 13 shows the effect of cough (E), nausea vomit (M) and shortness of breadth (Q) on sore throat (R). There is a low probability that a patient has sore throat given the parent symptoms. Note that this is as a result of the samples in the dataset, which does not capture the actual distribution over a large population.

Efficiency of Clustering

The learning of the symptoms pattern in the CDC COVID-19 dataset is efficient, as we can observe from Figure 14 where the top-3 principal components of the data points are projected in 3D for six classes for ease of visualization. Recall, the data points correspond to patients’ symptoms vector, and we can see that based on clustering, the data points are separated into distinguishable classes without overlapping. This means that there is a hidden structure in the data which can be discovered through data clustering.

Fig. 14 Illustration of the clustering of patients’ samples

To show that the clustering result is efficient, a deep neural network model is used to train the labeled clustered data which was obtained from clustering. For this task, the CDC dataset is split into training (80%), validation (15%), and testing (5%) sets. The model’s performance on the testing set is 99.47%. This means that when given a patient’s symptoms, the deep neural network model can predict the symptoms class with high accuracy. Figure 15 shows the accuracy (left) and loss (right) plots of the deep neural network model over the training epochs.

The symptoms structure in the top-5 classes from the clustering are presented in Table 4, which represents the most observed combination of patients’ symptoms in each class. For example, the patients in class 0 are Black Non-Hispanic that are 80+ years and are mostly diagnosed of abnormal chest X-ray, pneumonia and shortness of breadth, while the patients in class 3 are White Non-Hispanic that are within 70 years to 79 years and are mostly diagnosed of abnormal chest X-ray, cough, fever, pneumonia and shortness of breadth.

Fig. 15 Deep neural network model performance on clustered data

Performance Comparison of Time-series Forecasting Models

For the COVID-19 time-series dataset, the random forest and linear regression models are compared with the traditional ARIMA model in predicting the number of future confirmed cases and death (deceased) cases. The ARIMA model uses a p-value for auto regression, d-value for differencing, and q-value for moving average, and from observation in the experiment, it is best configured as (p,d,q) = (10,1,5) and (p,d,q) = (10,1,3) for the date+confirmed cases and date+deceased cases datasets respectively. The datasets are split into 0.9/0.1 for training/testing respectively.

To evaluate performance of the three models, the logarithm function of the RMSE is used (the lower the value, the better). As shown in Figure 16, the \(log(RMSE)\) value for the confirmed cases is high for all the models, while for the deceased cases, the \(log(RMSE)\) value is low for all the models because of the repeating trend in the dataset which is obvious in Figure 3 (right). On the date+confirmed cases dataset, the linear regression and random forest models significantly outperform the ARIMA model, which justifies the capability of ML over statistical or mathematical models in fitting complex time-series data.

Fig 16. Comparing model performance on time-series data

The random forest model and linear regression model can be smartly used to forecast the future impact of COVID-19 in any municipality and can also help to predict the stopping time of the virus if sufficient and accurate historical data is available.

Research Impact to SoReMo Mission

With the proposed methodology of ML with BN structure and parameter learning (i.e CPD) from the COVID-19 dataset, there is so much information to extract that can help the society, government and health agencies to be more proactive to the pandemic. First, as shown in this paper, the society can now be aware of the relationship between symptoms that cause the COVID-19 disease and their likelihood. This means that when a person experiences certain symptoms, they can predict the disease or other symptoms that can cause it and the probability. Second, regarding the death variable in the dataset, when a patient suffers from pneumonia there is need to immediately take the patient to an ICU with mechanical ventilation to save the patient. Third, policy makers can provide certain demographic health-related guidelines to patients based on the clustering result in this paper which shows that patients within the same age group and race have similar symptoms responsible for COVID-19. This can also help in the administering of vaccines to diverse groups within the population, especially in a city like Chicago. Lastly, all the results from the BN learning, data clustering and time-series forecasting can help medical professionals to effectively study other diseases that cause respiratory disorder or death.

In the context of SoReMo mission, this paper has attempted to improve on the health and safety awareness of the people of Chicago by employing a detailed modeling and algorithm design for COVID-19. This is a continuing work, as data of other diseases related to COVID-19 will be analyzed in future works and explained in a similar fashion.