# Effect of a Process Mining based Pre-processing Step in Prediction of the
Critical Health Outcomes

Negin Ashrafi
Department of Industrial and Systems Engineering
University of Southern California
Los Angeles, CA 90089
<EMAIL_ADDRESS>
& Armin Abdollahi
Department of Electrical Engineering
University of Southern California
Los Angeles, CA 90089
<EMAIL_ADDRESS>
Greg Placencia
Department of Industrial and Manufacturing Engineering
California State Polytechnic University-Pomona
Pomona, CA 91768
<EMAIL_ADDRESS>
& Maryam Pishgar
Department of Industrial and Systems Engineering
University of Southern California
Los Angeles, CA 90089
<EMAIL_ADDRESS>

###### Abstract

Predicting critical health outcomes such as patient mortality and hospital
readmission is essential for improving survivability. However, healthcare
datasets have many concurrences that create complexities, leading to poor
predictions. Consequently, pre-processing the data is crucial to improve its
quality. In this study, we use an existing pre-processing algorithm,
concatenation, to improve data quality by decreasing the complexity of
datasets. Sixteen healthcare datasets were extracted from two databases -
MIMIC III and University of Illinois Hospital - converted to the event logs,
they were then fed into the concatenation algorithm. The pre-processed event
logs were then fed to the Split Miner (SM) algorithm to produce a process
model. Process model quality was evaluated before and after concatenation
using the following metrics: fitness, precision, F-Measure, and complexity.
The pre-processed event logs were also used as inputs to the Decay Replay
Mining (DREAM) algorithm to predict critical outcomes. We compared predicted
results before and after applying the concatenation algorithm using Area Under
the Curve (AUC) and Confidence Intervals (CI). Results indicated that the
concatenation algorithm improved the quality of the process models and
predictions of the critical health outcomes.

_Keywords_ Critical health outcomes $\cdot$ concatenation algorithm $\cdot$
deep learning $\cdot$ predictions $\cdot$ process mining

##  1 Introduction

A patient’s health records can be understood as a sequence of observations
including services performed, diagnoses, or lab measurements [1]. In the
process mining field, each of these observations is called an event, and the
sequences of observations are referred to as event logs. Process mining
analyzes and optimizes processes. The process mining approach has been used
widely in healthcare to enhance healthcare processes [2] and to predict
critical health outcomes [1, 3, 4, 5].

Real-life event logs are complex and noisy. Infrequent behaviors and various
concurrences generate inefficient and complex process model through process
discovery algorithms. Healthcare data exacerbates these issues. For instance,
various lab measurements can be taken at the same time causing concurrent
relationships among these different lab measurement events. It is therefore
critical to pre-process raw event logs to improve their quality, to make the
process model more understandable.

A variety of pre-processing algorithms have been used to improve data quality
[6, 7, 8]. The concatenation algorithm is a recent innovation that
significantly improves event logs and process model quality for real-life
benchmark datasets [7]. However, the concatenation algorithm has not been
tested on real-life healthcare event logs.

Early prediction of patient health outcomes is critical to clinical decision
making and effective allocation of medical resources [9]. For example,
Coronary Artery Disease (CAD) causes approximately 610,000 deaths yearly in
the United States alone and 17.8 millions deaths worldwide [10, 11]. Moreover,
approximately, 6 million American adults over the age of 20 had heart failure
(HF) between 2015 to 2018 [12, 13, 14]. This has resulted in the Centers for
Medicare & Medicaid Services penalizing hospitals nearly 2.5 billion USD since
2012 for unplanned 30-day readmission rates [15]. Patient mortality for
Paralytic Ileus (PI) can be as high as 40 in ICU settings [16]. In 2019,
diabetes directly caused 1.5 million deaths. More recently, COVID-19, caused
by the SARS-CoV-2, heavily impacted morbidity and mortality with major
economic consequences [17]. These highlight the need for accurate predictions
models to reduce patient mortality rates and readmission, so as to reduce
financial consequences.

Many frameworks have been used to predict the critical health outcomes [18,
19]. Decay Replay Mining (DREAM) is a recent process mining / deep learning
framework. It has significantly improved predicting healthcare outcomes for
ICU patients with various diseases [3, 5] using raw and complex healthcare
event logs. However, the effectiveness of pre-processing has not yet been
studied using real life healthcare datasets. This work focuses on two issues:
1) Does pre-processing healthcare datasets improve data quality and enhance
the process model? 2) Is the concatenation algorithm an effective pre-
processor that predicts critical health outcomes such as patient mortality and
readmissions?

The structure of this paper is as follows: Section 2 provides the
preliminaries. Section 3 details the methodology used in this study. Section 4
presents the evaluation, covering Datasets, Setup, Results, and Discussion.
Section 5 concludes the findings and directions for future research.

##  2 Preliminaries

The notations used in this section are adopted from [20].

###  2.1 Process Mining

Process mining analyzes and optimizes process sequences called event logs.
Process mining uses the following steps: process discovery, conformance
checking, and process enhancement. The process discovery step is the most
important step in process mining. This step uses process discovery algorithms
to output a process model, mostly called a Petri Net (PN).

###  2.2 Petri Net

A PN is a mathematical model often used to represent a process. It includes a
set of places graphically represented as circles and transitions as
rectangles. Transitions correspond to events. The occurrence of a process is
modeled by firing transitions. Places can be marked or unmarked. Marked places
indicate current process states. A place is marked by putting a dot (token) in
the circle corresponding to that event. For a transition to be fired, all its
input places must contain a token. If a transition fires, the all tokens are
removed from its inputs. Also, each output of the transition receives a token.

###  2.3 Event Logs

An event $a$ $\epsilon$ $A$ is an instantaneous change of a process’ state
where $A$ is the finite set of all possible events. An event $a$ can happen
more than once in a provided process. An event instance $E$ is a vector with
two attributes at minimum: the name of the associated event $a$ and the
corresponding occurrence timestamp, $\tau$. The timestamps $\tau$ of two
events cannot be equal. A trace $g$ $\epsilon$ $G$ is a finite and ordered
sequence of event instances [20].

###  2.4 Decay Replay Mining

Decay Replay Mining (DREAM) [21] is a process mining / deep learning-based
approach used to predict the next event in the sequence of a process. DREAM
uses places from a process model and extends these places with a time decay
function $F(\tau)$. The time decay function $F(\tau)$ uses timestamps $\tau$
information as a parameter when replaying event logs $L$ on a process model
and produces Timed State Samples (TSS).

###  2.5 Timed State Samples

A TSS will be produced by replaying event logs $L$ on a process model [21]. A
TSS, $S(\tau)$ at time $\tau$ results from concatenating time decay function
values $F(\tau)$, token counts $C(\tau)$, and of a marking $M(\tau)$:

$S(\tau)=F(\tau)\oplus C(\tau)\oplus M(\tau)$ (1)

###  2.6 Concatenation Algorithm

The concatenation algorithm [7] is a pre-processing algorithm that removes
some concurrences and self-loops based on a probability algorithm. This
improves the quality of event logs and the resulting process model.

###  2.7 Concurrency Relation

Given any event logs $L$, and any two labels $\ell$, $\ell\prime$ $\epsilon$
$L$, and any two events $e_{i}$, $e_{j}$ $\in$ $E$, a concurrent relation
exists, denoted by $(e_{i}||e_{j})$, given the following conditions:
$|\ell\to\ell\prime|>0$ and $|\ell\prime\to\ell|>0$.

###  2.8 Self-loops

A self-loop exists in the event logs $L$, if $|\ell\to\ell|$ is positive for
some event $e$ $\in$ $E$ with $\lambda(e)$ $=$ $\ell$ [22].

###  2.9 F-Measure

An F-Measure is an evaluation metric uses to measure the quality of a process
model. The F-Measure is a trade-off value between fitness and precision [23].
The Fitness of a process model represents the behavior of the event logs $L$.
Precision shows the capability of a model to produce the behaviors found only
in the event logs $L$. In this work, the alignment-based method is used to
measure both fitness and precision per [20, 24].

###  2.10 Complexity

Complexity is an evaluation metric that measures the quality of a process
model. Complexity shows how easily a process model can be understood. Several
metrics measure the complexity of a model [22], including Control-Flow
Complexity (CFC) [23], size of the models [25], and structuredness [25]. CFC
shows how many branches are prompted by split gateways in a process model. The
size of a process model calculates the numbers of the nodes and arcs in the
model. A process model is structured, if, for any split gateway in the process
model a corresponding join gateway exits.

###  2.11 Split Miner

Split Miner (SM) is a process discovery algorithm that produce a sound PN from
the event logs $L$. It produces a process model with less complexity and
higher F-Measure value than other process discovery algorithms. The SM
algorithm includes the following steps: First, generate a directly-follow
graph and identify short loops. Second, discover the existing concurrency
relation between events. Third, utilize filtering by introducing 2 threshold
values, to control the filtering process (frequency threshold, $\epsilon$),
the other, $\eta$, to control concurrency relations. Moreover, the algorithm
derives choice and concurrency relations by adding split gateways. In the end,
the join gateways are discovered. The output of these steps is a BPMN, and it
can be converted to a PN through using ProMs’ BPMN Miner package [26]. This
algorithm is publicly available as a java application [26].

##  3 Methodology

This section focuses on a proposed method for predicting several health
outcomes with and without applying a pre-processing step on event logs $L$.
First, feature selection is described. Next, conversion of Electronic Health
Records (EHR) to event logs $L$ is explained followed by the pre-processing
step. Finally, the DREAM algorithm is introduced to predict health outcomes.

The intuition behind this methodology is that applying an appropriate pre-
processing step on complex healthcare event logs $L$ decreases complexity of a
process model. Hence, it predicts critical health outcomes more accurately.
Since healthcare event logs $L$ contain many concurrences and self loops,
applying a concatenation algorithm to concurrent events makes the process
model simpler and generates more accurate predictions. An overview of the
proposed framework is shown in Figure1.

Figure 1: Overview of the Proposed Model

###  3.1 Feature Selection

MIMIC III [27] and University of Illinois Hospital (UIH) databases[18] were
used for predictions. Due to the nature of our approach (process mining / deep
learning approach) every patient must have a medical history available at the
time of prediction. Table1 shows the features extracted from MIMIC III and UIH
databases by disease category.

Table 1: Features Extracted from MIMIC III and UIH Databases

Data Type (Part 1) Type of Disease Admission Insurance Lab Related
Demographics Elixhauser Comorbidities CPT and ICD-9-CM codes Discharge PI x x
x x NA NA x HF x x x x x x x COVID x x Diabetes x x x x x x NA CAD x x x x x
NA NA

Data Type (Part 2) Type of Disease Severity scores Care unit in and out from
ICU Location Encounter-based Report-based Vitals signs ICU-based PI NA x HF x
COVID x x x x x Diabetes x NA CAD NA NA

###  3.2 Conversion of EHR to the Event Logs

Event logs $L$ can be understood as a sequence of events with timestamps
$\tau$ associated to events that have occurred. Events represent activities
performed on patients such as admission, diagnosis, lab measurements, etc.,
also known as patients careflows. Each unique event has a corresponding
timestamp $\tau$ when it occurred. In this paper, the timestamp $\tau$ of
comorbidity and artificial events correspond to the discharge time of the
corresponding hospital admission. Since events in a trace $g$ occur
sequentially, multiple comorbidities or artificial events with the same
timestamps $\tau$ are delayed by multiples of 1 ms to retain their order.
Table2 summarizes the conversion of the EHR to the event logs $L$.

Table 2: Conversion of EHR from Databases to Event Logs

MIMIC III and UIH Tables Events ADMISSIONS, PATIENTS Admission type,
Demographics, patient’s exit type, discharge, insurance type D_LABITEMS,
LABEVENTS, ADMISSIONS Lab events ICUSTAYS, CALLOUT Careunit type LABEVENTS
Abnormal flagged events Admissions, Diagnoses_ICD Elixhauser Comorbidity Score
events Diagnose_ICD, Procedures_ICD, CPTevents 30 Artificial event
abstractions

###  3.3 Pre-processing Step

In Algorithm1, the resultant event logs ($L$) from the previous step were then
used as an input to the concatenation algorithm. Assume that events $e_{i}$,
$e_{j}$, $e_{k}$, and $e_{l}\in L$, if $(e_{i}\rightarrow e_{j})$ and
$(e_{j}\rightarrow e_{i})$; also $(e_{k}\rightarrow e_{l})$ and
$(e_{l}\rightarrow e_{k})$, if there is a concurrent relations between
$e_{i}$, $e_{j}$, $e_{i}||e_{j}$, and between $e_{k}$, $e_{l}$,
$e_{k}||e_{l}$.

Algorithm 1 Event Log Concatenation Algorithm

0: Event logs $L$, threshold $p^{*}$

0: Modified event logs $L^{\prime}$

1: Initialize an empty list of concurrent pairs

2: for all events $e_{i},e_{j}\in L$ do

3: if $(e_{i}\rightarrow e_{j})\land(e_{j}\rightarrow e_{i})$ then

4: Add $(e_{i}||e_{j})$ to the list of concurrent pairs

5: end if

6: end for

7: Initialize an empty list of valid pairs

8: for all pairs $(e_{i}||e_{j})$ in the list of concurrent pairs do

9: if $P(e_{i}\rightarrow e_{j})>p^{*}$ and $P(e_{j}\rightarrow e_{i})>p^{*}$
then

10: $P(e_{i}||e_{j})=P(e_{i}\rightarrow e_{j})+P(e_{j}\rightarrow e_{i})$

11: Add $(e_{i}||e_{j},P(e_{i}||e_{j}))$ to the list of valid pairs

12: end if

13: end for

14: Sort the list of valid pairs in descending order of $P(e_{i}||e_{j})$

15: for all pairs $(e_{i}||e_{j})$ in the sorted list of valid pairs do

16: Remove $e_{i}$ and replace $e_{j}$ with $e_{i}\oplus e_{j}$

17: end for

18: for all events $e_{i}\in L$ not concatenated do

19: if $e_{i}\in e_{i}\oplus e_{j}$ then

20: Remove $e_{i}$ and replace with $e_{i}\oplus e_{j}$

21: end if

22: end for

23: for all events $e_{i}\in L$ do

24: if $e_{i}\rightarrow e_{i}$ then

25: Remove $e_{i}$

26: end if

27: end for

28: Feed the modified event logs $L^{\prime}$ to the SM algorithm for process
modeling

29: Use the resultant process model as input to the DREAM algorithm for
prediction

We then calculate $P(e_{i}\rightarrow e_{j})$ and $P(e_{j}\rightarrow e_{i})$;
and $P(e_{k}\rightarrow e_{l})$ and $P(e_{l}\rightarrow e_{k})$. Moreover,
threshold value $p^{*}$ is introduced, if $P(e_{i}\rightarrow e_{j})>p^{*}$
and $P(e_{j}\rightarrow e_{i})>p^{*}$, the combination is chosen. Similarly,
if $P(e_{k}\rightarrow e_{l})>p^{*}$ and $P(e_{l}\rightarrow e_{k})>p^{*}$,
the combination is chosen for the next step. Additionally, the probabilities
of the selected combinations were added and sorted in the descending order,
$P(e_{i}||e_{j})=P(e_{i}\rightarrow e_{j})+P(e_{j}\rightarrow e_{i})$ also
$P(e_{k}||e_{l})=P(e_{k}\rightarrow e_{l})+P(e_{l}\rightarrow e_{k})$.
Furthermore, if $P(e_{i}||{e_{j}})=P(e_{k}||e_{l})$ a re-position step was
applied. After the order of combinations for concatenation is finalized, the
concatenation algorithm selects the ordered combinations one by one for the
concatenation. Assume that $(e_{i},$ $e_{j})$ is chosen for concatenation,
event $e_{i}$ is removed and replaced with the concatenated event,
$e_{i}\oplus e_{j}$. This step is repeated until all selected combinations are
concatenated in the event logs $L$. After concatenating all combinations in
the event logs $L$ based on descending order, some events will remain that are
not concatenated. Each non concatenated event will be examined and checked for
event, $e_{i}\in$ $e_{i}\oplus e_{j}$. In that case, $e_{i}$ will be removed
and replaced with $e_{i}\oplus e_{j}$. If $e_{i}$ $\notin$ $e_{i}\oplus
e_{j}$, while the final trace remains unchanged. Finally, if $e_{i}\rightarrow
e_{i}$ exist, it is removed. The modified event logs $L$ was then fed to the
SM algorithm to produce a process model. The resultant process model was input
to the DREAM algorithm for prediction as detailed in the next subsection.

###  3.4 DREAM algorithm

The algorithm has three different phases [21]. The first phase is to create a
PN model from an event log $L$ and assign a decay function to each PN location
$f_{p}(\tau)$. The discovery log is then replayed, and feature arrays
comprising decay function response values, token movement counts, and resource
use are extracted. Lastly, we train a Neural Network (NN) to predict the next
event, using these feature arrays as training data. To simulate data values
that decrease over time, one can use decay functions. Financial domains,
physical systems, and population trend modeling are among the common
applications for these functions. We assign a linear decay function
$f_{p}(\tau)$ to every place in the PN produced by SM on an event log $L$.

Algorithm 2 PN Discovery and Decay Function Initialization

0: Event log $L$ with traces $\\{g_{1},g_{2},\ldots,g_{n}\\}$

0: Decay function parameters $\alpha,\beta$

0: Discovered PN with initialized decay functions

1: Step 1: PN Model Discovery

2: Use SM to discover a PN model from event log $L$

3: for each place $p$ in the PN model do

4: Associate place $p$ with a linear decay function
$f_{p}(\tau)=\max(\beta-\alpha\cdot\Delta_{p},0)$

5: end for

6: Step 2: Decay Function Enhancement

7: for each place $p$ in the PN model do

8: Initialize the decay function $f_{p}(\tau)$ to $0$ by setting
$\Delta_{p}=\infty$

9: end for

10: Step 3: Calculate Decay Rates

11: for each place $p$ in the PN model do

12: Calculate the maximum trace duration $\Delta_{\text{max}}(L)=\max_{1\leq
i\leq|L|}(\nu_{d_{ts}}(L_{i,\gamma(L_{i})})-\nu_{d_{ts}}(L_{i,1}))$

13: Define $v_{p}(g)$ as the number of tokens entering place $p$ during replay
of trace $g$

14: if $\max_{g\in L}v_{p}(g)\leq 1$ then

15: Set $\alpha_{p}=\frac{\beta}{\Delta_{\text{max}}(L)}$

16: else

17: Set $\alpha_{p}=\frac{\beta}{\text{mean}_{g\in L}\delta_{p}(g)}$

18: end if

19: end for

20: Step 4: Update Decay Functions with Calculated Rates

21: for each place $p$ in the PN model do

22: Update the decay function
$f_{p}(\tau)=\max(\beta-\alpha_{p}\cdot(\tau-\tau_{p}),0)$

23: end for

24: Return the discovered PN model with initialized decay functions

By using $\Delta_{p}=\tau-\tau_{p}$, we can express the time difference
between the current time, $\tau$, and the most recent time a token entered
place p, $\tau_{p}$. With the help of the two decay function parameters,
$\alpha$ and $\beta$, we can adjust the significance level. Based on the
reactivation durations of a place p, $\alpha$ should ideally be set so that
the slope of $f_{p}(\tau)$ covers the entire range from $\beta$ to 0. In other
words, the slope should not be excessively steep for a small $\Delta_{p}$ such
that $f_{p}(\tau)=0$ or too flat for a large $\Delta_{p}$ such that
$f_{p}(\tau)=B$. This suggests varying $\alpha$ values with each place in the
PN model. In the next step of Algorithm2 we try to estimate these $\alpha_{p}$
values. Using the event log $L$ and the corresponding PN found by the SM on
$L$, we estimate $\alpha_{p}$. There are a limited number of event instances
in each trace $g$ in $L$. The $j^{\text{th}}$ event instance of an event log
L’s $i^{\text{th}}$ trace is referred to by $L_{i,j}$. $\lambda(L_{i})$
expresses the number of event instances of the $i^{\text{th}}$ trace of the
event log $L$. The value of the attribute d for the event instance E is
$\nu_{d}(E)$. When an attribute d is absent from an event instance E, the
result is $\nu_{d}(E)=0$, an empty set. We use $d_{ts}$ to represent the
attribute timestamp. $\Delta_{max}(L)$, the Maximum observed trace duration in
an event log $L$, is defined in the algorithm and based on the value of
$\max_{g\in L}v_{p}(g)$ the value of $\alpha_{p}$ is determined. If the value
is less than one then $\alpha_{p}$ will be set to a value that ensures that
before the final event instance of a corresponding trace $g$ in $L$ occurs,
the response of $f_{p}(\tau)$ will never equal 0. By doing this, we ensure
that, until the end of a trace, we will have information about the occurrence
of a particular event in the decay function’s response. However, when the
value of $\max_{g\in L}v_{p}(g)$ is larger than one we calculate the average
time it takes for a place p to reactivate. With this knowledge, we adjust the
decay rate so that, for the average time between reactivations, $f_{p}(\tau)$
yields a level of recent token movement importance.

Algorithm 3 Event Log Replay and Feature Extraction

0: Event log $L$, PN with places $P$, decay functions $f_{p}(\tau)$ for each
place $p\in P$

0: Set of timed state samples $S$

1: Initialize counting vector $C(\tau)$ to zero

2: for each trace $g\in L$ do

3: Reset $F(\tau)$, $C(\tau)$, and marking $M(\tau)$ to initial states

4: for each event $e_{i}$ in trace $g$ do

5: Calculate decay function responses $F(\tau)=\\{f_{p}(\tau)\ |\ p\in P\\}$

6: Update $C(\tau)$ for token movements

7: Update PN marking $M(\tau)$ based on transition fired by $e_{i}$

8: Construct timed state sample $S(\tau)=F(\tau)\oplus C(\tau)\oplus M(\tau)$

9: Add $S(\tau)$ to set of timed state samples $S$

10: end for

11: end for

12: Return set of timed state samples $S$

As a result, the slope will not be overly flat or steep. In Algorithm3
$F(\tau)$ is a vector of decay functions for all places. $C(\tau)$ shows the
number of tokens that have entered each place from time 0 to $\tau$. When a
token enters place p and time $\tau$ the corresponding value in $C(\tau)$
vector is updated. Here $M(\tau)$ is a vector representing the marking of a
PN. So for each trace $g$ and events in traces, the value of the
aforementioned vectors are updated, and then the timed state sample $S(\tau)$
is constructed by concatenating the vectors. Consequently, a timed state
sample, $S(\tau)$ provides a time-dependent description of a PN process state.
It includes data regarding time-based token movements, such as token counts
per location (loop information), the current PN state using the marking, and
the last time a token entered a place in relation to the current time. In the
last step of the algorithm, the $S(\tau)$ is input to the different NN with
different architecture for predicting the next event.

###  3.5 Prediction

DREAM was used to predict the critical health outcomes for both pre-processed
and non-processed event logs $L$. DREAM replays the event logs $L$ on a
process model and produces time information related to the variables which are
called TSS. The TSS, demographic information, and severity scores were fed
into the dense NN for prediction. The NN architecture differs by disease but
remaine the same for both pre-processed and unprocessed event logs $L$.

####  3.5.1 CAD disease

TSS are fed into a unique branch of one hidden layer with $200$ neurons as
shown in Figure2. The Dropout (DO) rate after the first layer was regularized
at DO = 0.2. Demographic information was fed into a single hidden layer with
$20$ neurons and the same dropout rate. These were concatenated into a second
layer with $90$ neurons and the same dropout rate. The critical outcome which
is predicted in this case is the mortality of ICU patients with CAD.

Figure 2: Architecture of the NN to predict CAD related outcome

####  3.5.2 HF disease

TSS, demographics information, and severity scores were fed separately to
three branches each with three hidden layers (Figure3). A batch normalization
layer was added after the first hidden layer for each branch with a dropout
rate of 0.2. The output layer included a softmax activation function to
predict unplanned 30-day readmission of ICU HF patients. The critical outcome
predicted in this case was unplanned 30-day readmission of ICU patients with
HF.

Figure 3: Architecture of the NN to predict HF related outcome

####  3.5.3 Diabetes

Each input layer is fed to an individual hidden layer branch before a hidden
layer concatenates the branches and feeds it to two further layers. All hidden
layers use a Rectified Linear Unit activation function (ReLU). The branched
first hidden layers also have a batch normalization layer. Dropout layers with
a rate of 0.4 for regularization are used. The output layer consists of a
softmax activation function to output the patient’s probability of in-hospital
mortality for diabetes patients. This architecture is shown in Figure4.

Figure 4: Architecture of the NN to predict Diabetes related outcome

####  3.5.4 COVID- 19

TSS, demographics information and comorbidities were fed separately to two
branches where the first branch contained three hidden layers with 90, 50 and
20 neurons respectively (Figure5). The first and second hidden layers each had
a dropout layer with a rate of 0.2. The second branch also contained one
hidden layer with 5 neurons. The two branches were concatenated to a branch
with three hidden layers, containing 90, 50, and 20 neurons respectively. The
dropout layer after the second concatenated hidden layer had a rate of 0.3.
The output layer used a softmax activation function to predict mortality of
COVID- 19 patients every 6-hour during the first 72 hours after hospital
admission.

Figure 5: Architecture of the NN to predict COVID-19 related outcome

####  3.5.5 PI

In Figure6, TSS were fed into two hidden layers. The first had 76 neurons and
a subsequent dropout rate of DO = 0.5. The second hidden layer contained 20
neurons. The demographic information was fed into a single hidden layer of 5
neurons. Both layers were concatenated into two further layers with 96 and 10
neurons respectively. Between these two layers was a layer with a dropout rate
of DO = 0.5. The outcome predicted in this case was the mortality of ICU
patients diagnosed with PI 24 hours after being admitted to the ICU.

Figure 6: Architecture of the NN to predict PI related outcome

##  4 Evaluation

This section discusses experimental evaluation of the model to predict the
critical health outcomes. Subsection one describes the datasets. Subsection
two describes the modeling set up. Subsection three highlights the results.

###  4.1 Datasets

Data was obtained from MIMIC III (Medical Information Mart for Intensive Care)
[27], a large database containing information relating to patients admitted to
Beth Israel Deaconess Medical Center (BIDMC). This data includes vital signs,
medications, laboratory measurements, observations and notes charted by care
providers, fluid balance, procedure codes, diagnostic codes, imaging reports,
hospital length of stay, survival data, and more [27]. Data was also obtained
from UIH, a tertiary, teaching hospital in Chicago. This study was approved by
The University of Illinois at Chicago (UIC) Institutional Review Board [18]. A
summary of the datasets for each disease category is shown in Table3.

Table 3: Summary of datasets by disease category

Disease Total Numbers of Patients Train Test Validation HF 3,411 2,422 555 434
COVID-19 508 303 104 101 PI 1,153 681 336 136 CAD 2,176 1,202 674 300 Diabetes
2,435 1,552 609 274

###  4.2 Setup

RapidProM was used to run experiments [24] extending RapidMiner with process
mining analysis capabilities. A concatenation algorithm’s threshold value of
$p^{*}$ was used to select the combinations of concurrent events for further
concatenation steps.

The optimal value of $p^{*}$ was based on the highest F-Measure value. The
hyper-parameter optimization was conducted on 16 datasets using steps of 0.1.
Of these, 13, found an optimal value of $p^{*}=0.7$.

SM takes two threshold values: $\eta$ and $\epsilon$. These values were
optimized based on the highest F-Measure values using steps of 0.1. Thresholds
selected were $\eta=0.4$ and $\epsilon=0.1$.

Train and validation sets are required to discover a process model and train
the NN. The validation set is used to select the best model. A test set is
then used to evaluate model performance. A summary of the set up for training
the NN for each disease category is shown in Table4.

Table 4: Summary of setup for training NN

Disease Activation Function Epoch Batch Size Optimizer HF ReLU 100 12 Adam
COVID-19 ReLU 350 10 Adam PI ReLU 350 50 RMSprop CAD ReLU 300 56 RMSprop
Diabetes ReLU 200 256 Adam

We measured the quality of the process models with and without pre-processing
using common metrics of fitness, precision, and F-Measure as proxies for
accuracy, size, CFC, and structuredness, which indicate complexity.

The metric used to evaluate the model prediction is the AUC score which is
equal to the probability that a classifier ranks a randomly chosen positive
instance higher than a randomly chosen negative one. AUC is a better
classification estimate than other common classification performance metrics
[28]. Higher AUC scores indicate that a model is better at distinguishing
between discharged patients and those who die. Furthermore, the 95% confidence
interval (CI) for the AUC score is calculated using DeLong’s method [29].

###  4.3 Results

In our first experiments without pre-processing, we used the optimal threshold
values found for the SM algorithm, plus raw event logs to discover process
models. The model was fed to the DREAM algorithm to produce a TSS. The TSS,
severity score and demographics were then fed to a NN for prediction.

For our second experiments, with pre-processing we applied the concatenation
algorithm to raw event logs $L$ using optimal thresholds values for the
algorithm to produce modified event logs $L$. The modified event logs $L$ were
then fed to the SM algorithm with the aforementioned optimal threshold values
to discover a process model, which was then fed to the DREAM algorithm to
produce the TSS. The TSS, severity score and demographics were then fed to a
NN for prediction.

Results of the evaluation are summarized in Table5. The F-Measure is
calculated for different types of disease before and after the concatenation
algorithm. In Table6, AUC and CI are calculated and in Table7 complexity is
measured by evaluating the size of the models, sructuredness and CFC. We ran
the Wilcoxon Test to determine whether improvements observed for evaluation
metrics applied after concatenation were statistically significant or not
($p=0.05$). Of 6 pre-processed datasets, 9 showed statistically significant
differences in AUC scores. All 16 showed statistically significant differences
in F-Measures and complexity metrics.

Table 5: Results in terms of F-Measure on event logs, and the Wilcoxon test
results on improved cases to statistically identify significance level

Evaluation Metrics Type of Disease F-Measure before concat after concat
p-value HF 0.821 0.830 0.048 PI 0.842 0.861 0.039 COVID 6hr 0.844 0.871 0.037
COVID 12hr 0.844 0.871 0.037 COVID 18hr 0.844 0.871 0.037 COVID 24hr 0.844
0.871 0.037 COVID 30hr 0.844 0.871 0.037 COVID 36hr 0.844 0.871 0.037 COVID
42hr 0.844 0.871 0.037 COVID 48hr 0.844 0.871 0.037 COVID 54hr 0.844 0.871
0.037 COVID 60hr 0.844 0.871 0.037 COVID 66hr 0.844 0.871 0.037 COVID 72hr
0.844 0.871 0.037 Diabetes 0.791 0.792 0.055 CAD 0.832 0.854 0.045

###  4.4 Discussion

To reduce concurrences and self-loops of the complex healthcare data, we
applied a concatenation algorithm as a pre-processing step to improve data
quality to improve the trustworthiness of prediction models for clinicians. We
observed significant statistical improvements in AUC values, F-Measure and
complexity metrics.

This was especially evident for the CAD and COVID-19 datasets which were more
complex compared to other datasets. In those cases the predicted results
improved by significant statistical amount of 2% for the AUC metric from our
experiments.

One of the main advantages of a process mining approach is that it discovers a
detailed process model which allows clinicians to comprehensively visualize
the entire process patients (careflow) go through in a medical system. This is
precisely why process mining has been preferred over other traditional machine
learning methods. Process mining also accurately models temporal time-related
information related to variables and utilizes them effectively as an input to
the neural networks when other traditional techniques cannot. Lastly, process
mining leverages the extensive and detailed medical history of patients from
prior hospital admissions.

A limitation of the proposed approach is that it requires patients’ medical
history that small clinics and hospitals might not have. It also excludes
patients with no prior hospital admissions like new patients and patients who
are not admitted before the current hospital visit. Moreover, hospitals tend
not to share patient data across outside hospital networks. Thus, the proposed
approach best matches large hospital networks. In addition, the training,
validation, and test sets all came from the MIMIC-III and UIH datasets.
Therefore, using an independent dataset from a different health center or
hospital would better test the performance of the model for future research.

Table 6: Results of the proposed model in terms of AUC and CI on event logs,
and the Wilcoxon test results on the improved cases to statistically identify
the level of significance

Evaluation Metrics AUC CI Type of Disease before concat after concat p-value
before concat after concat HF 0.930 0.947 0.045 [0.898- 0.960] [0.910- 0.971]
PI 0.810 0.823 0.047 [0.768- 0.840] [0.788- 0.850] COVID 6hr 0.776 0.779 0.055
[0.678- 0.876] [0.701- 0.790] COVID 12hr 0.782 0.79 0.043 [0.685- 0.880]
[0.708- 0.850] COVID 18hr 0.806 0.802 NA [0.719- 0.901] [0.788- 0.850] COVID
24hr 0.799 0.791 NA [0.698- 0.890] [0.718- 0.870] COVID 30hr 0.814 0.841 0.034
[0.718- 0.910] [0.786- 0.855] COVID 36hr 0.814 0.799 NA [0.718- 0.900] [0.778-
0.890] COVID 42hr 0.817 0.817 NA [0.701- 0.870] [0.789- 0.860] COVID 48hr
0.806 0.82 0.042 [0.738- 0.890] [0.784- 0.880] COVID 54hr 0.853 0.869 0.039
[0.758- 0.910] [0.788- 0.890] COVID 60hr 0.843 0.842 NA [0.768- 0.870] [0.748-
0.890] COVID 66hr 0.875 0.881 0.046 [0.778- 0.940] [0.798- 0.950] COVID 72hr
0.900 0.915 0.045 [0.868- 1.00] [0.870- 0.990] Diabetes 0.873 0.861 NA [0.851-
0.940] [0.788- 0.890] CAD 0.871 0.890 0.044 [0.831- 0.913] [0.870- 0.923]

Table 7: Results of in terms of complexity metrics on event logs, and the
Wilcoxon test results on the improved cases to statistically identify the
level of significance

Evaluation Metrics Size Struct. CFC Type of Disease before concat after concat
p-value before concat after concat p-value after concat before concat p-value
HF 70 46 0.032 142 121 0.034 71 66 0.032 PI 151 89 0.037 301 190 0.036 88 67
0.035 COVID 6hr 1001 401 0.031 2325 998 0.033 142 71 0.031 COVID 12hr 1001 401
0.031 2325 998 0.033 142 71 0.031 COVID 18hr 1001 401 0.031 2325 998 0.033 142
71 0.031 COVID 24hr 1001 401 0.031 2325 998 0.033 142 71 0.031 COVID 30hr 1001
401 0.031 2325 998 0.033 142 71 0.031 COVID 36hr 1001 401 0.031 2325 998 0.033
142 71 0.031 COVID 42hr 1001 401 0.031 2325 998 0.033 142 71 0.031 COVID 48hr
1001 401 0.031 2325 998 0.033 142 71 0.031 COVID 54hr 1001 401 0.031 2325 998
0.033 142 71 0.031 COVID 60hr 1001 401 0.031 2325 998 0.033 142 71 0.031 COVID
66hr 1001 401 0.031 2325 998 0.033 142 71 0.031 COVID 72hr 1001 401 0.031 2325
998 0.033 142 71 0.031 Diabetes 68 35 0.034 121 101 0.043 70 69 0.035 CAD 98
79 0.032 132 97 0.031 73 66 0.038

##  5 Conclusion

Healthcare data is complex and have many concurrences and loops. This makes
predicting healthcare outcomes difficult, because data quality is poor. By
using concatenation algorithm as a pre-processing step, data complexity
decreased, improving the performance of process discovery algorithms by
recommending a probability-based concatenation of events, with concurrent
relations. Thus, it can predict critical healthcare outcomes more accurately.

Predicting critical health outcomes allows medical stakeholders to provide
optimal treatment for patients to increase life expectancy. For instance, if a
patient is predicted to die or readmit, investing many resources on them will
not be optimal. Moreover, in cases when we can predict a patient’s outcome
sooner, medical stakeholders can expedite treatment to save a patient’s life.
Furthermore, knowing outcomes with high AUC would allow medical teams improve
treatment decisions and allocate medical resources more effectively.

We used 16 healthcare datasets and fed them to the SM algorithm to discover
the process model. Moreover, the pre-processing step was also applied to the
same event logs, and the resultant event logs were used by the SM algorithm to
discover a process model. We then measured the F-Measure and complexity of the
models for both cases. We fed both pre-processed and raw data to the DREAM
algorithm for prediction, and compared the results. In cases where
improvements were observed in evaluation metrics, the Wilcoxon Test was run to
measure the level of significance. Results indicated that the pre-processing
step increased F-Measure values and decreased the complexity of the models
statistically. Prediction of critical health outcomes improved as well. These
results demonstrate that using a pre-processing step improves data quality,
hence improves the process model and prediction modeling results. Using an
effective process mining / deep learning algorithm like the DREAM algorithm
also, ensures medical history and time information related to variables are
not ignored which can predict outcomes more accurately.

In future work, we plan to use other pre-processing approaches to see their
effects on data quality and predicting critical health outcomes. We plan to
find techniques that optimally adjust parameter values automatically, rather
than using hyper-parameter optimization currently in favor. We would also
validate the proposed approach using healthcare datasets. Finally, MIMIC-III
and UIH datasets contains useful information like clinical notes and images
that can be converted into event logs and fed to models as inputs.

## Appendix: Notations

$a$

Event

$A$

Finite set of all events

$\alpha$

Decay rate

$\alpha_{p}$

Decay rate for a specfic place p

$\beta$

Constant parameter of a decay function

$C(\tau)$

Token counting vector from time 0 to $\tau$ each element represents the number
of tokensthat entered a specific place

$d$

Attribute

$d_{ts}$

Timestamp attribute

$D$

Finite set of all possible attributes

$\Delta_{max}(L)$

Maximum observed trace duration in an event log $L$

$\Delta_{p}$

Difference of $\tau$ and $\tau_{p}$

$\delta_{p}(g)$

function of average time between a token is consumed in place p until a new
token is produced in p based on an input trace g

$E$

Event instance vector

$e$

An instance of E

$\eta$

SM threshold which controls concurrency relations

$\epsilon$

SM threshold which controls filtering process

$f_{p}(\tau)$

Decay function of place $p$

$F(\tau)$

Decay function response vector

$g$

Case or trace

$G$

Finite set of all possible traces

$L$

Event log which is a set of traces

$L_{i,j}$

$j^{\textit{th}}$ event instance in the $i^{\textit{th}}$ trace of an event
log $L$

$\ell$

label of each event

$\lambda$

function that retrieves the label of each event

$\lambda_{p}(g)$

Number of tokens a place $p$ produces when replaying a trace $g$

$M(\tau)$

Vector representing the marking of a PN

${mean()}$

Arithmetic mean function

$N$

Set of all possible event instances

$\nu_{d}(E)$

Function returning value of attribute $d$ of event instance $E$

$p$

Place

$P$

Set of all places

$PN$

Mathematical definition of a labeled PN

$p\star$

Threshold value for concatenation algorithm

$S(\tau)$

Timed state sample at time $\tau$

$S$

Set of timed state samples

$\tau$

Time

$\tau_{p}$

Most recent time that a token entered place $p$

$v_{p}(g)$

number of tokens a place p produces when replaying a trace g

$e_{i}\oplus e_{j}$

Concatenated event for $e_{i}$ and $e_{j}$

$p(e_{i}\rightarrow e_{j})$

The probability that ej happens instantly after ei

$p(e_{i}||e_{j})$

probability of the frequency that event ej which occurs immediately after
event ei

## Acknowledgment

The authors would like to thank all contributors of MIMIC-III for providing
the public EHR dataset.

## References

  * [1] Pishgar, M., Theis, J., Del Rios, M. et al. Prediction of unplanned 30-day readmission for ICU patients with heart failure. BMC Med Inform Decis Mak 22, 117 (2022). https://doi.org/10.1186/s12911-022-01857-y
  * [2] Ghasemi M, Amyot D. Process mining in healthcare: a systematised literature review. Int J Electron Healthc. 2016;9:60.
  * [3] M. Pishgar, M. Razo, J. Theis and H. Darabi, "Process Mining Model to Predict Mortality in Paralytic Ileus Patients," 2021 International Conference on Cyber-Physical Social Intelligence (ICCSI), 2021, pp. 1-6, doi: 10.1109/ICCSI53130.2021.9736217.
  * [4] J. Theis, W. Galanter, A. Boyd and H. Darabi, "Improving the In-Hospital Mortality Prediction of Diabetes ICU Patients Using a Process Mining/Deep Learning Architecture," in IEEE Journal of Biomedical and Health Informatics.
  * [5] Pishgar, M., Harford, S., Theis, J. et al. A process mining- deep learning approach to predict survival in a cohort of hospitalized COVID-19 patients. BMC Med Inform Decis Mak 22, 194 (2022). https://doi.org/10.1186/s12911-022-01934-2
  * [6] M. Sani, S. van Zelst, and A. Van Der, “Improving the performance of process discovery algorithms by instance selection,” Comput. Sci. Inf. Syst., vol. 17, no. 3, pp. 927–958, 2020.
  * [7] M. Pishgar, M. Razo and H. Darabi, "Improving Process Discovery Algorithms Using Event Concatenation," in IEEE Access, vol. 10, pp. 69072-69090, 2022, doi: 10.1109/ACCESS.2022.3185235.
  * [8] S. Baharlouei, S. Patel, and M. Razaviyayn, “f-FERM: A scalable framework for robust fair empirical risk minimization,” in Proc. Int. Conf. Learn. Representations (ICLR), 2024.
  * [9] N. Ashrafi, Y. Liu, X. Xu, Y. Wang, Z. Zhao, and M. Pishgar, “Deep learning model utilization for mortality prediction in mechanically ventilated ICU patients,” medRxiv, Mar. 2024. [Online]. Available: https://doi.org/10.1101/2024.03.20.24304653.
  * [10] JC. Brown, TE. Gerhardt, E. Kwon . Risk Factors For Coronary Artery Disease. 2021 Jun 5. In: StatPearls [Internet]. Treasure Island (FL): StatPearls Publishing; 2022 Jan–. PMID: 321192
  * [11] N. Ashrafi, A. Abdollahi, G. Placencia, and M. Pishgar, “Process mining/deep learning model to predict mortality in coronary artery disease patients,” medRxiv, Jun. 2024. [Online]. Available: https://www.medrxiv.org/content/10.1101/2024.06.26.24309553v1.
  * [12] Mozaffarian DBE, Go A, Arnett D, Blaha M, Cushman M, et al. Heart disease and stroke statistics-2016 update: a report from the American Heart Association. Circulation. 2016;133:e38–360.
  * [13] J. Zhang, H. Li, N. Ashrafi, Z. Yu, G. Placencia, and M. Pishgar, “Prediction of in-hospital mortality for ICU patients with heart failure,” medRxiv, Jun. 2024. [Online]. Available: https://doi.org/10.1101/2024.06.25.24309448.
  * [14] S. Mohan, C. Thirumalai and G. Srivastava, "Effective Heart Disease Prediction Using Hybrid Machine Learning Techniques," in IEEE Access, vol. 7, pp. 81542-81554, 2019, doi: 10.1109/ACCESS.2019.2923707.
  * [15] Desautels T, Das R, Calvert J, Trivedi M, Summers C, Wales DJ, et al. Prediction of early unplanned intensive care unit readmission in a UK tertiary care hospital: a cross-sectional machine learning approach. BMJ Open. 2017;7(9): e017199.
  * [16] L.Wei, X.Yongtao, H. Jianhu, L. Lina, T. Yiqing, Y. Weihui, C.Yi, and C. Wei, ”Causes and prognosis of chronic intestinal pseudo-obstruction in 48 subjects: a 10-year retrospective case series,” Medicine 97, 2018.
  * [17] Taubenberger JK, Morens DM. 1918 Influenza: the mother of all pandemics. Rev Biomed. 2006;17(1):69–79.
  * [18] Galanter, W., Rodríguez-Fernández, J.M., Chow, K. et al. Predicting clinical outcomes among hospitalized COVID-19 patients using both local and published models. BMC Med Inform Decis Mak 21, 224 (2021). https://doi.org/10.1186/s12911-021-01576-w.
  * [19] M. Chen, Y. Hao, K. Hwang, L. Wang, and L. Wang, "Disease prediction by machine learning over big data from healthcare communities," _IEEE Access_ , vol. 5, pp. 8869–8879, 2017.
  * [20] W. van der Aalst, “Data science in action,” in Process Mining: Data Science in Action, W. van der Aalst, Ed., Berlin, Germany: Springer, 2016, pp. 3–23.
  * [21] J.Theis, and H.Darabi,"Decay Replay Mining to Predict Next Process Events," IEEE Access, 7, 119787–119803, 2019.
  * [22] J. Mendling, Metrics for Process Models: Empirical Foundations of Verifi- cation, Error Prediction, and Guidelines for Correctness. Cham, Switzer- land: Springer, 2008.
  * [23] J. Cardoso, “Business process control-flow complexity: Metric, evalua- tion, and validation,” Int. J. Web Services Res., vol. 5, no. 2, pp. 49–76, Apr. 2008.
  * [24] W.M.vanderAalst,A.Bolt,andS.J.vanZelst,“RapidProM:Mineyour processes and not just your data,” 2018, arXiv:1703.03740.
  * [25] J. Mendling, H. A. Reijers, and W. M. P. van der Aalst, “Seven process modeling guidelines (7PMG),” Inf. Softw. Technol., vol. 52, pp. 127–136, Feb. 2010, doi: 10.1016/j.infsof.2009.08.004.
  * [26] R.Confort.BPMNMiner: Automated Discovery of BPMN Process Models with Hierarchical Structure. Accessed: Aug. 1, 2022. [Online]. Available: https://github.com/raffaeleconforti/ResearchCode
  * [27] Johnson AE, Pollard TJ, Shen L, Lehman LWH, Feng M, Ghassemi M, Mark RG. MIMIC-III, a freely accessible critical care database. Sci Data. 2016;3(1):1–9.
  * [28] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, “Dropout: A simple way to prevent neural networks from overfitting,” J. Mach. Learn. Res., vol. 15, no. 1, pp. 1929–1958, Jan. 2014.
  * [29] S. M. Lundberg, and L.Su-In. “A unified approach to interpreting model predictions.” Advances in Neural Information Processing Systems, 2017.


# Non-Uniform Gaussian Blur of Hexagonal Bins in Cartesian Coordinates

Reinier Vleugels
Department of Computer Science,
IBIVU Centre for Integrative Bioinformatics,
Vrije Universiteit Amsterdam
and
Magnus Palmblad
Center for Proteomics and Metabolomics,
Leiden University Medical Center  The first author is an MSc student at the
Department of Computer Science, IBIVU Centre for Integrative Bioinformatics,
Vrije Universiteit Amsterdam. This work was conducted as part of his thesis.

###### Abstract

In a recent application of the Bokeh Python library for visualizing physico-
chemical properties of chemical entities text-mined from the scientific
literature, we found ourselves facing the task of smoothing hexagonally binned
data in Cartesian coordinates. To the best of our knowledge, no documentation
for how to do this exist in the public domain. This short paper shows how to
accomplish this in general and for Bokeh in particular. We illustrate the
method with a real-world example and discuss some potential advantages of
using hexagonal bins in these and similar applications.

Keywords: Binning, blurring, Bokeh, histograms, tesselation

##  1 Introduction

Hexagonal binning is a popular alternative for creating two-dimensional
histograms that captures most shapes better than the more common rectilinear
binning methods (Carr et al., 1987). Hexagons are more similar to circles than
squares which means data is aggregated more tightly around the bin center. The
minimum distance from a bin center to the nearest bin border is ~7.5% longer
in hexagonal bins than in equiareal square bins, which likely improves
accuracy of bin selection when interacting with histograms using a computer
mouse or touchscreen. The major disadvantage is the added complexity of
computation and visualization. It is also impossible to exactly subdivide or
aggregate hexagonal bins into smaller or larger hexagons, which is easy to do
with rectangular bins. However, in some applications, these drawbacks are not
critical.

One of the added complexities when working with hexagonal bins is the
coordinate system and transformations between the hexagonal and Cartesian
coordinate systems. There are several hexagonal grid systems, such as the
axial or trapezoidal coordinates used by the Bokeh Python library (Bokeh
Development Team, 2019), the offset coordinates used by the TikZ Shapes LaTeX
library (also to generate Figures 1-5 in this paper) and cube coordinates with
three axes. Other software, including the R tess package (Höhna, 2013), use a
one-dimensional offset enumeration to address individual hexagons. In this
paper we will focus on offset and axial coordinates, though any hexagonal
coordinates or numbering system can be transformed into Cartesian coordinates.

It is sometimes advantageous to smooth images or histograms to reduce noise
and emphasize features of interest. Blurring allows data visualization at a
higher resolution without requiring a large numbers of counts in each bin.
Blurring can also be used to represent uncertainly in measurements or
predictions. When the hexagonal tiles represent spatial data such as a map or
image, blurring can be done in the hexagonal coordinate system, with the
distance between two tiles being the minimum number of steps between them.
Each hexagon thus have 6 neighboring tiles with distance 1, 12 tiles of
distance 2, 18 of distance 3, etc - in general $6n$ tiles in the $n$-th shell
around a given tile (Figure 1). Blurs based on this distance metric are
approximately radial.

(0,0)(0,1)(0,2)(1,0)(1,1)(2,1)(2,0)(0,-1)(0,-2)(1,-1)(1,-2)(2,-1)(-1,-1)(-1,-2)(-2,-1)(-1,0)(-1,1)(-2,1)(-2,0)$x$$y$
Figure 1: Hexagonal tiles with offset hexagonal coordinates $(x,y)$ in two
shells around (0,0). The $x$ axis is here drawn in the middle of the tiles
with $y=0$.
(0,0)(0,-1)(0,-2)(1,-1)(1,-2)(2,-2)(2,-1)(0,1)(0,2)(1,0)(1,1)(2,0)(-1,1)(-1,2)(-2,2)(-1,0)(-1,-1)(-2,0)(-2,1)$q$$r$
Figure 2: Hexagonal tiles with axial (trapezoidal) coordinates $(q,r)$ in the
same two shells around (0,0). This is the hexagonal coordinate system used by
Bokeh.

However, in one of our use cases, based on previous work (Palmblad, 2019), the
two dimensions do not represent spatial locations as in an map or image, but
two different physical variables, with different underlying uncertainties. To
smooth the histogram and visualize these uncertainties, we have to use the
original physical dimensions or coordinates. As we still wish to use hexagonal
bins, for the reasons mentioned above, we have to evaluate the smoothing
function in the Cartesian coordinates of the bins, e.g. the bin centers. Note
that we do not wish to smooth the image. The hexagonal bins should remain
distinct to allow the user to interact with the binned data. Here we show how
this can be done, and demonstrate the result using Bokeh. For a more
comprehensive and interactive overview of the many different hexagonal grid
systems, see Amit Patel’s excellent blog post on the subject (Patel, 2013).

##  2 Methods

A hexagon can be seen as made up by six equilateral triangles. From elementary
trigonometry, we know the ratio of the height to the side in an equilateral
triangle is $\sin{60^{\circ}}$ or $\sqrt{3}/2$. If we define the length of the
sides of the regular hexagons as 1, the center-to-center distance between two
adjacent hexagons, that is two hexagons sharing one edge, is therefore
$\sqrt{3}$. Between hexagonal bins in offset rows or columns, the distance is
3/2 in one dimension and $\sqrt{3}/2$ in the other (Figure 2). The
transformations from offset hexagonal ($H$) coordinates (Figure 1) to
Cartesian ($C$) ones are thus:

$\quad\begin{bmatrix}x\\\
y\end{bmatrix}_{C}=\begin{cases}\begin{bmatrix}3/2&0\\\
0&\sqrt{3}\end{bmatrix}\begin{bmatrix}x\\\ y\end{bmatrix}_{H}&\text{if
}x_{H}\equiv 0{\@displayfalse\pmod{2}}\\\ \begin{bmatrix}3/2&0\\\
0&\sqrt{3}\end{bmatrix}\begin{bmatrix}x\\\
y\end{bmatrix}_{H}+\begin{bmatrix}0\\\ 1/2\end{bmatrix}_{C}&\text{if
}x_{H}\equiv 1{\@displayfalse\pmod{2}}\end{cases}$ (1)

With an angle $60^{\circ}$ between the axes, the transformation from the
hexagonal coordinates $q$ and $r$ in Figure 2 to Cartesian coordinates is
independent of coordinate parity:

$\quad\begin{bmatrix}x\\\ y\end{bmatrix}_{C}=\begin{bmatrix}3/2&0\\\
-\sqrt{3}/2&-\sqrt{3}\end{bmatrix}\begin{bmatrix}q\\\ r\end{bmatrix}_{H}\quad$
(2)

This transformation is illustrated in Figure 3 below. The Cartesian
coordinates here refer to the centers of the hexagonal bins:

(0,0)(0,-1)(1,-1)(2,-1)(0,0)(0,$\sqrt{3}$)($\frac{3}{2}$,$\frac{\sqrt{3}}{2}$)(3,0)$\sqrt{3}$1
Figure 3: Four hexagonal tiles in axial hexagonal coordinates as in Figure 2
(left) and in corresponding Cartesian coordinates in units of hexagon side
length (right). The Cartesian coordinates corresponding to $(q,r)$ in axial
coordinates are $(3q/2,\sqrt{3}(q/2-r))$. The Cartesian coordinates
corresponding to $(x,y)$ in offset hexagonal coordinates (Figure 1) are
$(3x/2,\sqrt{3}y)$ for even columns ($x\mod{}2=0$) and
$(3x/2,\sqrt{3}y+\sqrt{3}/2)$ for odd ($x\mod{}2=1$).

As the two Cartesian dimensions in our use case represent different physical
dimensions with different uncertainties of measurement or estimation, the
smoothing should not have rotational symmetry like a radial blur. If the
errors are independent, the Gaussian kernel or signal $S$ that should be
transferred from $(0,0)$ to $(x,y)$ in Cartesian coordinates is:

$\vspace{12pt}S(x,y)=\frac{1}{{\sigma_{x}\sigma_{y}\sqrt{2\pi}}}e^{{{-{x}^{2}{y}^{2}}\mathord{\left/{\vphantom{{-\left({x}\right)^{2}}{2\sigma_{x}^{2}}{2\sigma_{y}^{2}}}}\right.\kern-1.2pt}{2\sigma_{x}^{2}}{\sigma_{y}^{2}}}}\vspace{12pt}$
(3)

where $\sigma_{x}$ and $\sigma_{y}$ are the standard deviations in the
Cartesian coordinates corresponding to the two physical dimensions. We only
need to evaluate $S(x,y)$ once for any pair of $\sigma_{x}$ and $\sigma_{y}$
and in practise only for a limited number hexagonal bins (Figure 4):

10039.339.32.42.422.322.30.20.27.27.27.27.20.20.20.20.251.951.951.951.92.62.62.62.65.55.55.55.50.30.30.30.30.10.10.10.10.10.10.10.1
Figure 4: Gaussian blur of a single bin (red circle) in the case
$\sigma_{x}=2$ and $\sigma_{y}=1$ in units of hexagon side lengths. The
numbers represent the transferred signal relative to the blurred value of the
original bin (100%). All other bins are below 0.1% in this example.

We implemented a generic function combining Equations 2 and 3 for creating a
blurred histogram data matrix in Python as part of SCOPE (Search and Chemical
Ontology Plotting Environment). The SCOPE project includes tools for executing
literature searches, collecting text-mined chemical entities of biological
interest and creating interactive visualizations. The histograms were
constructed using pandas and visualized by Bokeh version 1.4.0. Bokeh is
particularly popular for creating interactive visualizations for web browsers
and scales well to large datasets, including the millions of named entities
that can be retrieved from a single literature search. The Bokeh interface in
SCOPE includes sliders to adjust the blurring and scaling, and displays the
most frequent chemical classes in any selected bin.

##  3 Verifications

To visually verify that the Gaussian blurring does what it was designed to do,
we can look at the effect of varying $\sigma_{x}$ and $\sigma_{y}$ in the
kernel (Eq. 3):

39.339.32.42.422.322.30.20.27.27.27.27.251.951.951.951.92.62.62.62.65.55.55.55.50.30.30.30.30.10.10.10.10.10.10.10.10.20.20.20.2100

79.279.239.339.368.768.722.322.351.951.951.951.984.984.984.984.940.140.140.140.148.448.448.448.422.822.822.822.815.715.715.715.716.816.816.816.822.322.322.322.3100

Figure 5: Examples of a single blurred hexagonal bin using the numbers in
Figure 4 (top) and with $\sigma_{x}=4$ and $\sigma_{y}=2$ (bottom),
restricting the blur to the same cells. The color mapping in both cases is
linear from white at 0% to black at 100% of the blurred value of the central
bin.

We tested our Gaussian blurring implementation for Bokeh using SCOPE and
searching Europe PMC for scientific papers mentioning particular techniques
from analytical chemistry in their section-tagged materials and methods
sections (Palmblad, 2019). Hexagonal bins more faithfully traced the shapes of
log P/mass distributions of chemical classes and the resulting histograms
exhibited different distributions as expected from what is known about the
applicability of the analytical techniques:

Figure 6: Distributions of mass (Da) and polarity (log P) of small molecules
retrieved from papers in Europe PMC mentioning hydrophilic interaction
chromatography (HILIC) in their methods section. The distributions are
visualized by SCOPE before (left) and after (right) applying a non-uniform
Gaussian blur of size 2.5. The saturation was set to 2.5, term frequency-
inverse document frequency normalization applied and the Viridis colormap
selected.

The sliders for adjusting the blurring and saturation improve the user
experience and help create clear visualizations. The blurring convey the
uncertainty in the log P predictions. We also found interactive bin selection
to be easier in hexagonal bins than in equiareal rectangular bins. Other
features and applications of SCOPE are beyond the scope of the this paper (no
pun intended), but will be described elsewhere.

##  4 Conclusion

In this short paper we have shown how to calculate and apply non-uniform
Gaussian blurs in Cartesian coordinates in hexagonal bins, and illustrated
this with a real-world example. For simplicity, we evaluated the Gaussian
kernel in the center of each hexagonal bin. It may be argued that the average
kernel value in the bin more accurate represents the underlying distribution.
However, this would require numerical integration over all these hexagonal
bins for what would likely be a very small visual or practical benefit. All
source code is available on GitHub (ReinV/SCOPE) under the Apache 2.0 license.

## References

  * Bokeh Development Team (2019) Bokeh Development Team (2019).  Bokeh: Python library for interactive visualization.  Bokeh Development Team.  https://bokeh.org/.
  * Carr et al. (1987) Carr, D. B., R. J. Littlefield, W. L. Nicholson, and J. S. Littlefield (1987).  Scatterplot matrix techniques for large n.  Journal of the American Statistical Association 82(398), 424–436.
  * Höhna (2013) Höhna, S. (2013).  Fast simulation of reconstructed phylogenies under global time-dependent birth–death processes.  Bioinformatics 29(11), 1367–1374.
  * Palmblad (2019) Palmblad, M. (2019).  Visual and semantic enrichment of analytical chemistry literature searches by combining text mining and computational chemistry.  Analytical Chemistry 91(7), 4312–4316.
  * Patel (2013) Patel, A. (2013).  Hexagonal grids blog.  https://www.redblobgames.com/grids/hexagons/.