key: cord-1051970-x5z64uh1 authors: Augusto, Adriano; Deitz, Timothy; Faux, Noel; Manski-Nankervis, Jo-Anne; Capurro, Daniel title: Process Mining-Driven Analysis of COVID-19’s Impact on Vaccination Patterns date: 2022-05-04 journal: J Biomed Inform DOI: 10.1016/j.jbi.2022.104081 sha: b541f3f57ae80c5e82cacaef1a75927c2d56f29a doc_id: 1051970 cord_uid: x5z64uh1 Process mining is a discipline sitting between data mining and process science, whose goal is to provide theoretical methods and software tools to analyse process execution data, known as event logs. Although process mining was originally conceived to facilitate business process management activities, research studies have shown the benefit of leveraging process mining in healthcare contexts. However, applying process mining tools to analyse healthcare process execution data is not straightforward. In this paper, we show a methodology to: i) prepare general practice healthcare process data for conducting a process mining analysis; ii) select and apply suitable process mining solutions for successfully executing the analysis; and iii) extract valuable insights from the obtained results, alongside leads for traditional data mining analysis. By doing so, we identified two major challenges when using process mining solutions for analysing healthcare process data, and highlighted benefits and limitations of the state-of-the-art process mining techniques when dealing with highly variable processes and large data-sets. While we provide solutions to the identified challenges, the overarching goal of this study was to detect differences between the patients‘ health services utilization pattern observed in 2020–during the COVID-19 pandemic and mandatory lock-downs –and the one observed in the prior four years, 2016 to 2019. By using a combination of process mining techniques and traditional data mining, we were able to demonstrate that vaccinations in Victoria did not drop drastically–as other interactions did. On the contrary, we observed a surge of influenza and pneumococcus vaccinations in 2020, as opposed to other research findings of similar studies conducted in different geographical areas. • To overcome the identified challenges, we have proposed: i) a novel method to identify trace boundaries in incomplete process traces having unknown start and end events; and ii) a novel method to fix timestamp equivalence issues in event log data by leveraging domain expertise. • Our process mining analysis on the largest real-life dataset of general practice healthcare processes, containing more than 31 million entries over a time of 45 months, demonstrates the benefits and the limitations of applying process mining in general practice day-to-day healthcare processes. • From a clinical perspective, the results of our analysis show a surge in influenza and pneumococcus vaccinations of adults and elderly people, in 2020, during the COVID-19 pandemic. These findings are in contrast with other research studies conducted in different geographical areas. The discipline of Process Mining [1] was born with the goal to design automated data analysis techniques that could support the phases of the business process management lifecycle [2] , especially those phases where the data analysis plays a central role, e.g., process discovery and process monitoring. Over the past two decades, research in 5 the area of process mining has generated a number of methodologies and software tools (henceforth, process mining techniques) [3, 4, 5, 6] . Process mining techniques usually require process execution data, which is known as event logs. It is possible to distinguish two major families of process mining techniques [2] : operational techniques and tactical techniques. The former family encompasses techniques whose goal is to gen-10 erate insights in real-time during the process execution, e.g., estimating the probability of a negative event to happen; or the likelihood of a specific process outcome. The latter family encompasses techniques whose goal is to help analysts to discover, analyse, and periodically monitor the process execution in order to understand how the process is performed, what are its weaknesses, and how the process can be improved. Two ing the COVID-19 pandemic in Victoria (Australia). Our research goal was to identify 25 differences between the patients' health services utilization pattern observed in 2020during the COVID-19 pandemic and mandatory lock-downs -and the one observed in the prior four years, 2016 to 2019. Given that health services are provided via the enacting of healthcare processes, process mining techniques are ideal for achieving our research goal, in particular, process discovery and process variant analysis techniques. 30 To this end, we analysed process execution data extracted from more than 100 general practice (GP) clinics in Victoria. This data included more than 30 million events capturing the GP healthcare processes of more than one million patients in Victoria, over a time-span of approximately five years. While, in this study, we do not build on top of our findings to lead process improvement initiatives, given that no clinics or hospitals 35 were directly involved, it is worth noting that the outcome of our process analysis can be exploited to improve healthcare services during emergency periods, for instance, enhancing healthcare resources allocation (including personnel, medications, vaccines, and services) and ensuring short processing times. The main contributions of this study are the following. 40 -The showcasing of a methodology to: i) prepare general practice healthcare process data for conducting a process mining analysis; ii) select and apply suitable process mining solutions for successfully executing the analysis; and iii) extract valuable insights from the results of the process mining analysis alongside leads for traditional data mining analysis. 45 -A critical assessment of the capabilities of state-of-the-art process mining techniques as well as their limitations when dealing with large data-sets recording highly variable processes -which is typical of the healthcare processes. This 55 allowed us to draw directions for future research in the area of applied process mining in healthcare. -From a medical perspective, the results of our analysis show that vaccinations in Victoria did not drop as drastically as other clinical interactions did. On the contrary, we observed a surge of influenza and pneumococcus vaccinations in 60 2020. Also, these findings differ from other research findings of similar studies conducted in different geographical areas in the equivalent seasonal periods [13, 14, 15] , providing a different perspective. The remainder of the paper is structured as follows. In Section 2 we discuss related work and background. In Section 3, we describe the data, the analysis we ran, 65 the challenges we faced, and the solution we adopted. In Section 4, we review the findings of the data analysis, providing a medical interpretation and considering their consequences. Lastly, Section 5 summarises our results and draws the conclusion. 70 range of process mining techniques such as automated process discovery, conformance checking, and process variant analysis. Automated process discovery techniques [20, 21, 22, 23, 3] allow one to discover patients' clinical pathways from the recordings of their healthcare process activities captured by the hospitals and clinics information systems [24, 16] . Conformance 85 checking techniques [25, 5] allow one to automatically compare the observed healthcare process behaviour (in the form of process execution data) against a prescribed process behaviour to identify differences between actual and normative healthcare behaviour. The latter is usually provided in the form of an imperative process model or as a set of declarative process rules [26] , which rather than capturing the full process be-90 haviour may describe clinical guidelines. Process variant analysis techniques [27, 28, 4, 29] allow one to automatically compare two or more sets of healthcare process executions exhibiting different outcomes (or performance) to identify relevant differences between the executions that may have had an impact on the outcome or performance of the healthcare process. These type of techniques are applied to answer questions 95 such as: what were the differences between the healthcare treatments provided by two different hospitals to patients having the same diagnosis? One of the earliest application of process mining in healthcare dates back to 2008, Mans et al. [24] used Heuristics Miner [20] to extract insights from healthcare process data, both from clinical and administrative perspective, including process han-100 dovers analysis by leveraging the process mining analytics platform ProM. 2 Poelmans et al. [30] used a combination of process mining and data mining techniques to detect and analyse differences in the healthcare pathways of patients treated for breast cancer and how they would respond to different therapies. Lakshmanan et al. [31] proposed an approach for discovering patients healthcare pathways and correlate them to techniques to analyse the quality and the costs of the healthcare services provided to 110 patients at one South Australian hospital. Roviani et al. [26] reported a case study on how to leverage declarative process mining techniques to identify divergences between clinical guidelines and the observed execution of clinical processes, at the urology department of the Isala hospital in the Netherlands. Leonardi et al. [8] proposed a method to abstract low-level process execution data (in the form of simple actions), turning it 115 into high-level data that can be used for process mining applications. They validated their method by discovering process models from healthcare services, showing that their method improved the graphical representation of the healthcare processes, and facilitated the clustering of similar process executions. Alvarez et al. [9] applied process mining techniques to discover process models capturing how healthcare professionals 120 operate within emergency rooms, analysing them to identify opportunities for process improvement. Chen et al. [10] proposed a framework to extract high-level descriptions of medical treatment processes from electronic medical records by applying clustering techniques on doctor order set sequences. Their framework allows to enrich the extracted process descriptions with additional information regarding the process per-125 formance (e.g., cost, length), providing support for improvement. Yang et al. [11] designed a process mining approach to automatically and in real-time detect process deviations from recommended clinical guidelines. They validated their approach on a set of pediatric trauma resuscitation procedures, demonstrating the effectiveness of their solution. All these studies on process mining in healthcare represent only a fraction of the existing ones, but reporting on all of them would require a separate study and it would be outside the scope of this one. Hence, we refer the interested reader to the latest literature reviews [19, 18 ]. Since the early months of the COVID-19 pandemic, the main drivers behind lockdowns and stay-at-home measures were the need to reduce face-to-face interactions to 140 prevent the virus from spreading uncontrollably, the subsequent increase in morbidity, mortality, and overwhelming of healthcare service providers. COVID-19 apparently increasing the risk of this diseases and attributed the change to reduced access to health services [34] . These findings were confirmed in the USA [35] . Similar effects were described for patients with acute myocardial infarction [36] , and cancer [37, 38] , among other conditions. This phenomenon was also observed for preventative care services such as cancer screening [39, 40, 41] , and maternal and child 155 health services [42] . In particular, there were growing concerns that a significant reduction in immunizations would result in an increase in vaccine-preventable conditions [43, 13] . The goal of this study was to analyse changes in health services utilization patterns during the 2020 COVID-19 pandemic and associated lock-downs in Victoria (Aus-we were able to solve some of these challenges, by proposing approaches that can be reused in different contexts, other challenges remain open or partially addressed and should be considered in future research work in the area of process mining. Before discussing our analysis, we provide some formal definition for the concepts we refer to throughout section. While we contextualised these definitions within our study, we remark that these are well-known definitions and concepts in the area of process mining [1] . Definition 1. Event -An event e captures the execution of an activity within a process 175 instance. An event can be represented as a tuple (x 1 , x 2 , . . . , x n ), where each element x i captures an attribute of the event, and at least three attributes are present: the process instance ID (c -event ID); the label of the activity the event refers to (aevent activity); and the timestamp (t -event timestamp). Additional attributes usually capture the process resource who executed the activity, customer information, etc. In 180 the following, given an event e, we will refer to its three required attributes with the notation e |c , e |a , e |t . Definition 2. Event Log -An event log L is a sequence of events e 1 , e 2 , . . . , e n , such that all the events are ordered by their timestamp. Formally, ∀e i ∈ L | i ∈ 185 Definition 3. Trace 3 -Given an event log L, a trace of the event log τ ∈ L is a sequence of events, τ = e 1 , e 2 , . . . , e n , such that all the events belong to the event log, all the events are ordered by their timestamp, and all the events have the same event We note that, according to Definition 3, we can also consider an event log as a 190 multiset of traces. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 Definition 4. Directly-follows Relation -Given an event log L = e 1 , e 2 , . . . , e n , we say that a directly-follows relation holds between any two events e i , e j ∈ L if and only if e i and e j belong to the same trace and j = i + 1, in other words, the two events follow each other in (at least) one trace. We indicate such a relation with the notation 195 e i → e j . Formally, given e i , e j ∈ L, e i → e j ⇐⇒ ∃τ ∈ L | e i , e j ∈ τ ∧ j = i + 1. We extend the concept of directly-follows relation to the event activities, i.e., if e i → e j then we say that also e i|a → e j |a holds. Definition 5. Directly-Follows Graph (DFG) -Given an event log L, its Directly- In other words, each node of the DFG represents a unique activity recorded in the event log, and each edge of the DFG represents a directly-follows relation between two activities -represented by the source node and target node of the edge. activities, and decisions involving actors and data objects triggered by a specific start event and leading to a specific end event (i.e., process outcome) that delivers value to a customer. In this study we used the Patron dataset [44] . This dataset stores de-identified patient data from the Patron primary care data repository (extracted from consenting general practices), that has been created and is operated by the Department of General Practice at The University of Melbourne [44, 45] . This dataset is aggregated from more than 100 General Practice (GP) clinics in Victoria (Australia) and includes both 215 administrative and clinical data, including all interactions between patients and their GPs, for more than one million patients. Access to the data was approved by the Mel- Medications; Investigations (Pathology and Imaging). While the first three tables contain information regarding the patient and their clinical history; the last three tables contain information regarding the patient healthcare processes, respectively: information on patient visits to and interactions with their GP doctor(s); information on patient drugs prescriptions; and information on patient pathology and imaging tests and re-225 sults. Looking at the latter three tables through the lens of Definition 1, an event ID corresponds to a patient ID, which identifies a unique patient accessing GP services across all the tables. An event activity corresponds to a medical activity the patient underwent. From the three tables capturing the patient healthcare process, it is possible to extract 230 seven medical activities, which are reported in Table 1 . These activities capture all the recorded interactions of a patient with their GP, including the drugs they have been prescribed, their pathology and imaging tests and results, and their vaccinations. For simplicity, we will refer to each of these seven activities by using a letter A to G (following the mapping in Table 1 ). Lastly, the event timestamp corresponds to the time a 235 medical activity was completed. We note that the Patron dataset does not record information regarding activities' lifecycle, e.g., the start and the completion of the activities, and that the timestamp granularity is at day-level (i.e., the smallest difference between timestamps is at day-level). Such a timestamp granularity is frequent in the healthcare contexts, and (at least in our case) it is related to how the system records events into the 240 database. Consequently, it was virtually impossible to infer a better timestamp granularity (e.g., hours and minutes) or the duration of a single medical activity (e.g., how long a GP visit would last). In light of this, in the Patron dataset, a trace captures a unique patient accessing GP services over the time, i.e., a process instance of the GP day-to-day healthcare process. The de-identified data was stored in a secure virtual machine. While this was a strict requirement for analysing the data, such a secure environment posed some challenges during the data analysis stage (discussed later in this section), mostly related to the fact that it did not allow for internet access. To conduct our analysis, we adhered to the methodology proposed by van Eck et al. [46] , adapting it to our context. The PM 2 methodology [46] has six stages: planning; data extraction; data processing; data mining and analysis; evaluation; and process improvement and support. We thoroughly executed all the stages with the exception of the last stage. Given that this study did not involve GP clinics and healthcare practitioners, 255 we did not have the means to implement a redesigned process, besides, it would have been outside the scope of this study. Following the PM 2 methodology, we started from the planning, which includes three steps: i) selecting the process to analyse; ii) determining the process analysis goal; 260 iii) and assembling a team. Indirectly, the Patron dataset drove the process selection. Given that it captures ambulatory patients' interactions with their GPs, we selected for our analysis the GP day-to-day healthcare process. Our analysis objective was to identify differences between the GP healthcare services provided in 2020-during the COVID-19 pandemic and mandatory lock-downs -and those observed in the prior four 265 years, 2016 to 2019. The authors of this paper composed the research team, bringing expertise in process mining, data mining, and (medical) general practice. Once the scope of our analysis was set, we moved to data extraction and processing. The Patron dataset, as mentioned above, already included all the data we required to analyse the selected process. The extraction of this data was performed outside this 270 study, and it is not our contribution. However, healthcare process data rarely comes 11 in the form of ready-to-use event logs [16, 33] , which is the required data format for conducting a process mining analysis [1, 46] , and the Patron dataset was no exception. During this stage, we focused on transforming the available data into an event log that could allow us to achieve our analysis goal. This required us to identify what entries of 275 the relational database were suitable to be turned into events. As mentioned above, we extract all the entries from three tables out of six, which captured the medical activities shown in Table 1 . Each entry of a table included the patient ID and the timestamp, hence, the conceptual mapping from table entries to events was straightforward. We note that this mapping was facilitated by the existing of a very extensive data dictionary 280 describing the Patron dataset, which often is not available. The data extracted captured a time-span of (almost) five years, from January 2016 to November 2020, but we reduced this time-span to keep only the data collected be- The data was extracted via an ad-hoc R-script and saved in the form of CSV event logs. These CSV logs were then converted in the standard XES format via Apromore 295 (academic version) 4 , which can be used without internet access and does not have limits on the amount of data to be processed, as opposed to Disco 5 or Celonis 6 . Alternatively, we could have converted the CSV event logs into XES format via ProM 7 4 https://apromore.org/academic-alliance/ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 Once we obtained the event logs from the Patron dataset, we proceeded to the data mining and analysis stage. By looking at the data through the lens of Definition 1, 2, and 3, we could summarise its characteristics as shown in Tables 2 and 3 . rarely observed. In all the logs (GP16-GP20), 98% of the distinct traces are observed between 1 and 5 times (see Table 3 ), and less than 0.03% of distinct traces are observed more than a thousand times. The traces include 31.8 million events, which -to the best of our knowledge -dwarf any of the real-life public logs used in automated process discovery research [3] . The trace length varies widely, with minimum, average, and 310 maximum length of 1, 12, and 2317 events (respectively). Given that our goal was to compare the patients behaviour in the months between March and November 2020 against the patients behaviour in the same timeframe of the past four years, we divided the log into five sublogs (namely, GP20, GP19, GP18, GP17, GP16), each of them capturing the 9-month timeframe in one of the five years 315 under analysis. Such an approach is common for performing process behavioural comparison -known in the area of process mining as process variant analysis [4] . Looking at Tables 2 and 3, we notice that dividing the GP16-20 log into five sublogs does not affect much the variety of the process behaviour. Although the absolute number of events and traces reduces, each of the five (sub)logs maintains remarkable characteris-320 tics; i.e., 5.1 million (GP20 log) to 6.9 million (GP19) events, and 401 thousand (GP20 13 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 log) to 520 thousand (GP17 log) traces (on average, 41% distinct). As a comparison, the largest real-life event log used in the series of business process intelligence challenges had 1.6 million events. 8 By analysing the characteristics of these five logs, we can immediately draw some initial observations. healthcare process execution was observed approximately two times over 9 months in all the five years) but also, and perhaps most importantly, by the distribution of distinct traces over a set frequency buckets. Table 3 shows that in each year, the percentage of 335 distinct traces that were observed between 1 and 5 times is constant at 98.2%. Observation 1 was expected, given that a strict lockdown was enforced in Victoria from 16-March to 21-June and from 04-July to 09-Nov, that possibly deterred patients from accessing healthcare for what they considered minor issues. Observation 2 had a surprising nature, in fact, intuitively, we would expect that the combination of 340 lockdown and pandemic would foster standardization in healthcare processes (i.e., less variability). We explored the distribution of the activities over time, their frequencies, and how they varied over the five years. This information is shown in Figure 1 . 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 Table 3 : Distribution of distinct traces over a set of frequency buckets recorded in the event log -in this case, the total number of activities is equivalent to the total number of events recorded in the event log. Figure in 2020 with respect to another year. From the data captured in these plots, we can observe the following. Observation 3. In 2020, the relative frequency of activity B (GP records a measurement) dropped to 9.0% from an average of 12.1%. Although this seems a small variation, we note that in 2019, 2018, 2017, and 2016, the relative frequency of activity B 360 was remarkably stable at 12.2%, 12.0%, 12.1%, and 12.2% (respectively). Observation 4. In 2020, the relative frequency of activity D (GP prescribes a refill) increased to 9.2% from an average of 5.0%. Also in this case, we note that in 2019, 2018, 2017, and 2016, the relative frequency of activity D was somewhat stable at 5.8%, 5.1%, 4.6%, and 4.4% (respectively). 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 However, taking into account the changes of absolute frequency for activity D 385 (see Figure 1j ), we can observe that drug prescriptions have increased steadily in the past four years with an average increase of 14.2%. Given that also a similar trend can be observed for activity C (capturing a first-time drug prescription), we cannot conclude that the increase observed in activity D derived exclusively from the COVID-19 pandemic context. Lastly, Observation 5 is probably the most interesting one, also because it is in contrast with research findings of similar studies conducted in different geographical areas during the equivalent seasons [13, 14] . The data clearly shows that vaccinations were not substantially impacted in 2020, with a decrease in absolute frequency that is lower than the one of other activities (see Figure 1j ). Leaving aside medication-related Table 2 -total events), the drop of vaccination activities is well below the average drop of the other medical activities. In addition, there is a noticeable shift in the vaccination timeline for the 2020 year, bringing the vaccinations forward of one month. Observation 5 set a direction for additional analysis, which led us to additional findings that we will discuss in depth in Section 4. Recalling the event log definition (Definition 2), given that the order of the events in 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 an event log is imposed by the order of their timestamps, incorrect or imprecise times-415 tamps can have a significant (negative) impact on the identification of directly-follows relations and, consequently, on the output of process mining tools that rely on directlyfollows relations. This is a well-known problem in the field of process mining [51, 52] , especially in healthcare [16] , where activities are documented manually. We recall that also in our case the event timestamps had a day-granularity. To give an idea of the issue, let us consider a patient visiting a GP doctor (activity A), the doctor measures the blood pressure of the patient (activity B), and then prescribes a medication for the first time (activity C). The activities order is A, B, C . However, they will be recorded in the information systems having all the same timestamp (i.e., the day of the visit), and not necessarily in the order they have been executed. For example, the fact that the patient has visited the doctor may be recorded at the end of a consultation, and the doctor may log activities B and C after they really occurred (inputting them manually on a computer software). As a result, the actual recording may read as follow C, A, B . The more the activities to be recorded, the more are the users involved in their (manual) logging, the greater is the amount of errors. Past research studies in process mining have addressed the problem of cleaning (or repairing) imprecise timestamps and timestamps errors [53, 54, 55, 56] , however, three of the proposed methods require as input a reference process model [53, 54, 55] , while the method of Conforti et al. [56] requires to have at least a subset of the events recorded in the event log that are not affected by imprecise timestamps. In our case, 435 we could not rely on any of these existing methods, missing their requirements. While recent work [12] called for improving the quality of the data captured by healthcare information systems, with the goal to fix the problem at its root, we would like to highlight the opportunity (and the need) for additional research addressing the problem of automated repairing and the cleaning of event log data errors -especially 440 timestamps. To continue our analysis and ensure the most reliable outcome, we devised an effective solution to deal with the imprecise timestamps. We imposed a standard order among the activities (matching the alphabetical order of their labels, see Table 1 ), and we reordered the events in the event log based on two attribute values: the event times- 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 tamp and the event activity. The latter attribute used as a tie-breaker on timestamps equality. For example, let us consider the following sequence of events e 1 , e 2 , e 3 , e 4 , e 5 , and let us assume that the five events (e 1 to e 5 ) have all the same timestamp and that the corresponding sequence of activities is D, A, G, F, E . In such a case, we would re-450 order the events as e 2 , e 1 , e 5 , e 4 , e 3 , yielding the sequence of activities A, D, E, F, G . Note that the event IDs do not play a role in the ordering. Events having the same ID will be ordered correctly, while events having different IDs would not be affected by the reordering. Our solution is based on the idea that, in most of the scenarios (and especially in 455 healthcare), certain activities have logical order constraints, e.g., a GP doctor cannot take and record a patient blood pressure (activity B) if the patient is not attending a visit (activity A). Yet, our solution has limitations, given that not all the activities have a logical order constraint, e.g., a patient may be administered a vaccine (activity G) either before or after she is prescribed a medication (activity C or D). In fact, there 460 are only three strict logical order constraints in our case, and they are: A before B, C before D, and E before F . The order we imposed satisfies the three constraints, but also enforces others. We note that, while enforcing additional constraints may distort the factual reality, it homogenise the data allowing for a correct and fair comparison. To describe the effects of our solution, let us consider two traces A, B, G and 465 A, G, B , and let us assume that all the events within each trace have the same timestamp. Comparing the two traces as they are would tell us that they are different, but according to the data they are not (i.e., the timestamps are equal, so any order is valid in principle). Enforcing a standard order over the activities as a tie-breaker on timestamp equality ultimately leads to data standardization and a correct interpretation. Our approach for fixing imprecise timestamps due to high-level granularity can be generalized to virtually any other context when the objective of the process analysis is the comparison of process variants, so it should not be considered as an ad-hoc approach for our specific scenario. However, we acknowledge that to define the logical order on the activities, the input of domain experts may be required. In our case, we 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 Dr Manski-Nankervis. Lastly, we note that the time complexity of our approach is linear on the number of events contained in the event log, making it not only effective but also efficient. Once we solved the problem of imprecise timestamps, we focused on the process behaviour, analysing how the process activities follow one another and what their execution leads to. However, we note that our healthcare process instances do not perfectly fit the traditional definition of process [2] (see Definition 6), because they miss both a specific start event and a specific end event, making these process instances unbounded. In our context, a patient may consult their GP doctor to discuss several health issues at once, each of them may lead to different outcomes and some of them may never reach an outcome (e.g., a chronic disease, which requires to be indefinitely monitored), forcing the customer to indefinite follow-ups. At the same time, while following health issues up, new health issues may arise. As one can see, the GP day-to-day healthcare 490 process is conceptually unbounded. In particular, when we look at the activities of a patient within a specific timeframe, the first activity we observe is not necessarily the one that started their GP day-to-day healthcare process, and the only way to determine that with 100% accuracy would be to have a timeframe at least equal to the patient age -which is an unrealistic requirement for most of the patients. Existing process mining techniques for automated process discovery and variant analysis (e.g., [22, 21, 27, 29] ) are not very effective when dealing with unbounded process instances, because by design they would implicitly (and erroneously, in our context) assume the first event of a trace in the input event log to be the start of the process instance, and the last event of a trace to be the end of the process instance. We can, 500 however, identify the most appropriate start and end events given a process instance. This can be achieved by narrowing down the scope of an unbounded process instance, for example, by focusing on a single GP visit or a single health issue/procedure. To do that we devised an algorithm that leverages domain experts knowledge, once again, the co-authors Dr Capurro and Dr Manski-Nankervis. 505 We started from the assumption that a process instance should begin with a visit 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 to the GP doctor (i.e., activity A), effectively making activity A the only possible start event of a trace. Any subsequent activity different than activity A (i.e., activities B to G) is assumed to be a follow-up of the initial visit to the GP doctor. However, when a second activity A is observed for the same process instance, we have to distinguish 510 two cases: i) the new activity A is a follow-up of the past activities; ii) the new activity A is not related to the past activities (i.e., this would trigger a new process instance). We distinguished the two cases on a time basis. Precisely, if the new activity A is more than six months away from the first observed activity A and more than one month away from the last observed activity of the current process instance, we are in case ii); If the event ID (e |c ) is not yet in the map Π and the event activity (e |a ) is in the set α, we create a new empty trace (τ ), we append e to τ , we add the event ID and the trace to the map Π, we save the timestamp of e in T 0 and in T n (lines 5 to 10). If e |c is already mapped in Π and e |a is not an allowed start activity (line 12), we retrieve the trace linked to the event ID (Π(e |c )) and we append e to that trace (line 13). Then, we update the timestamp information by overwriting the last observed event timestamp in the map T n (line 14). If e |c is already mapped in Π and e |a is an allowed start activity (line 12), we 535 distinguish the two possible cases mentioned above. Case ii), if e |t is less than or equal to ∆ 0 or less than or equal to ∆ n , then we append e to the already existing trace Π(e |c ) 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 (as just described above -see lines 16 to 18). Otherwise, Case i), we create a new event ID (that is not present in the event log), 9 we link the new event ID to the existing trace in Π that is mapped to e |c , we create a new empty trace (τ ), we append e to τ , we add We note that the information shown in Table 2 is the one obtained after the execution of Algorithm 1. The column filtered events reports the number of events that were removed by applying Algorithm 1, i.e., events that are not preceded by an activity A. On average, we removed 3.6% of events from the data, which is a negligible amount. At this stage, we can finally turn our attention to the process behaviour analysis, by leveraging process mining techniques [46] . Since we are interested in identifying process behavioural differences over five different timeframes (each captured in an event 555 log), the appropriate process mining techniques are in the class of automated process discovery [3] and process variant analysis [4] . Automated process discovery techniques receive in input an event log and automatically produce a process model, which is a graphical representation of the process behaviour, such as a workflow chart, a Petri net, or a BPMN model 10 . By looking at different process models, it is possible to 560 detect behavioural differences. On the other hand, process variant analysis techniques receive two event logs and automatically produce an artifact that highlights the process behavioural differences. Differences captured by variant analysis techniques are either at control-flow level (i.e., process behavioural differences in terms of executed 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 Add (e |c , e |t ) to T0; 10 Add (e |c , e |t ) to Tn; First, we attempted to discover a process model from each of the five event logs, by running each of the three automated process discovery tools. Fodina was not able to output a sound process model from any of the five event logs within the timeout. Inductive Miner discovered exactly the same process model from all the five logs, which 585 is the one capture in Figure 2a . Split Miner discovered the same process model from four logs (GP16-GP19, Figure 2b ), but a different one from the GP20 log ( Figure 2c) . The process models shown in these figures are modelled using the BPMN 2.0 standard notation. 12 For the readers who are unfamiliar with this notation, we remind that: the round graphical elements with thin and thick border represent respectively the start 590 and the end events of the process; the squared boxes represent the activities that can be executed during a process instance; a squared box marked with the symbol represent 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 an activity that can be consecutively executed several times, before allowing the process execution to continue; the diamonds with an X represent decision points (when having multiple outgoing arcs) or passive join points (when having multiple incoming 595 arcs); lastly each arc in a BPMN 2.0 model captures the allowed process execution flow. The process models we discovered are not structurally complex (e.g., spaghettimodels), but they clearly show that a large variety of behaviour is allowed -with few constraints and many cyclical patterns. This finding highlights that the GP day-to-day Miner from the GP20 log (Figure 2c ), the process model highlights that activity A must 610 be the first to be executed, but it can be repeated after executing other activities. The model correctly captures other behavioural constraints, for example, activity C must precede activity D, and activity G can follow only after activity D. However, the model allows all the activities to be repeated any number of times. The difference between the model discovered by Split Miner from the GP20 log and any of those discovered 615 from the GP16-GP19 logs is minimal, and mostly involving activity G. The latter is not anymore following activity D. From this automated process discovery exercise, we can see that the models dis- 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 We ran each of the three tools four times, providing in input four pairs of logs: (GP16, GP20); (GP17, GP20); (GP18, GP20); (GP19, GP20). By doing so, each tool 625 would identify the differences in the observed process behaviour between each pair of years, perfectly addressing our research and analysis goal. Also in this case, one of the three tools [28] was unable to provide an output within the timeout, or it was unable to identify any statistically significant difference. 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 an intuitive idea of the magnitude of the difference (the darker the shade the greater the difference), but not an exact estimate. An identical decoding reasoning applies for interpreting the color-coding of Figures 3c and 3d , with the only difference being the comparison metric, which is not the waiting time but the frequency of activities (i.e., 640 color-coded nodes) and the frequency of the directly-follow relations (i.e., color-coded arcs). The process comparator was able to identify statistically significant differences (p < 0.05) both in terms of time elapsed between two consecutive activities and the frequency of observing two consecutive activities. We remind that the process com- Interestingly, the pair (A, B) was observed much less in 2020 than in the previous two years. While these differences support Obsersations 3, 4, and 5 (see Section 3.5), the 655 process comparator detected many other differences that we could not justify or link to our initial observations. Once again, this can be related to the amount and variability 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 of the behaviour recorded in the event logs. The output of the process comparator requires some interpretation from the analyst, which sometimes is not straightforward. Instead, the tool of Cecconi et al. [29] 660 provides natural language descriptions of the identified differences. The top-10 differences identified by this tool [29] when comparing the logs GP20 and GP19 are reported in Table 4 . This output, similar to the previous one, refines Observation 3 and Observation 4 (see Section 3.5), highlighting a decrease in the number of observations of activity B in the healthcare processes of 2020 (see Table 4 , rows 5, 7 and 10), and 665 an increase in the number of observations of activity D in the healthcare processes of 2020 (see Table 4 , rows 1-4 and row 9). Similar results were obtained when comparing the data from the 2020 against the data from the 2018, 2017, and 2016. of space, clarity, and simplicity, we report the DFG of only two event logs (GP20 and GP19) and in their matrix form, where each matrix row (and column) represents a node of the DFG -i.e., an activity; and each cell of the matrix captures the frequency of the edge between the two nodes -i.e., how many times we observe in the event log a directly-follows relation between two activities. Tables 5 and 6 report the DFGs in matrix form of the event logs GP20 and GP19 (respectively). For clarity, in each DFG matrix cell, we have also reported in brackets the relative frequency of the directly-follows relation (i.e., the ratio of the number of times we observe a directly-follows relation over the total number of directly-follows 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 relations observed in an event log). We note that any directly-follows relation can 680 be observed in the two DFGs. Although some of them are rare (e.g., C → C, with a frequency in the order of hundreds), the vast majority can be observed with a frequency in the order of thousands -although they represent a tiny percentage of the total number of the observed directly-follows relations. According to this process data, the GP day-to-day healthcare process allows for executing medical activities in any 685 order. Consequently, the data generated from such kind of process is very complex and challenging to analyse with state-of-the-art process mining techniques -as we have seen with Fodina [21] and Taymouri's tool [28] , both unable to output sound results. An alternative would be to remove some behaviour by applying filtering techniques [59, 60, 61] , however, we recall that all the automated process discovery algo-690 rithms that we used already apply a filter [21, 23, 22] , as well as two of the three variant analysis approaches [27, 29] . While in a business context some infrequent behaviour may be a violation of com- 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 pliance or internal business rules, in our context, all behaviour is actually allowed. As such, we are not interested in removing behaviour, but rather in narrowing our focus, 695 and consider only a portion of behaviour that can be fruitfully analysed. With that in mind, we considered only the most frequent behaviour. Table 7 reports the top-20 most frequent traces that we could observe in each of the five event logs. Scanning carefully through Table 7 , we notice that in 2020, traces containing the activity G were more frequent than other years (for clarity, we reported these traces in 700 Table 8 ). In 2020, not only the traces containing the activity G were more frequent, but they accounted for the 25% percent of the most frequent behaviour (5 traces out of 20). This finding is remarkable, and when paired it with Observation 5 (discussed in Section 3.5) clearly hints to a variation in the behaviour involving vaccinations during the 2020. A similar reasoning can also be done for traces containing the activity D (see A, D , across the five logs). Further investigation of the most frequent process behaviour may reveal several additional differences, but within the scope of this study, we decided to investigate the specific behavioural difference related to activity G and its traces among the top-20. Our process mining analysis highlights two limitations of the state-of-the-art pro-710 cess mining techniques that we used in this study: 1. Process mining techniques for automated process discovery and process variant analysis suffer of scalability and/or quality issues when they deal with too much and too variable behaviour. 2. Automated process discovery techniques try to capture as much behaviour as 715 possible from the event log, filtering infrequent behaviour only when it is strictly necessary to either simplify the process model or increase its accuracy. However, depending on the context, one may be interested in capturing very little process behaviour from the event log -requiring special filtering techniques. While limitation 1 is ground for future research directions and studies. Limitation 720 2, at the moment, can be addressed manually, by applying ad-hoc filters of the process behaviour recorded in the event logs (as we did). The best ad-hoc filters must be identified by domain experts, often on a trial-and-error basis, and applied either via 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 ProM plugins or commercial tools such as Celonis, Disco, or Apromore (which we used). Future process mining techniques should allow the user to automatically design 725 such filters, without relying on domain experts knowledge. For instance, by automatically analysing the outputs of a set of process mining techniques (e.g., both automated process discovery and variant analysis) -as we did manually. Lastly, we analysed each of the process traces by looking into the distribution of 740 the executed activities over time, we reported this information in Figures 5a to 5f . At this point, it is evident that the differences in behaviour were not only in terms of how the activities were executed (i.e., their order and frequencies) but also when. We summarised these findings in the following observation. 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 Observation 6. In 2020, we can observe a clear (left-)shift (i.e., towards March) and 745 early peak in the distribution of the activities executed within the most frequent behaviour of the vaccination process, as well as a different trend when compared to the past four years, which holds for all the activities involved in the vaccination process. Furthermore, the most frequent vaccination process in 2020 was more complex than the previous four years, allowing for more behavioural variants with frequently requiring 750 additional activities (activity A, and B). The observed change can be explained by the recommendation that Australians receive their influenza vaccinations before the normal season (April-May). This recommendation was broadly advertised to minimize a possible double hit to the healthcare system: an epidemic of SARS-CoV-19 in addition to the usual Fall/Winter influenza 755 season. 13 In the next section, we will discuss more in depth this observation from a medical perspective. 13 https://www.abc.net.au/news/health/2020-04-01/australians-urged-to-get-flu-vaccination-coronaviruscovid- 19/12107264 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 Figure 4 , the order of the activities is captured in brackets. Activities that are optional (i.e., may not be executed, according to the process behaviour) are explicitly mentioned. Our process mining analysis, including the extraction of process data, its plotting (e.g., Figure 1 ), and the application of state-of-the-art process discovery and process 760 variant analysis tools, provided us with a clear picture of the GP day-to-day healthcare process as well as a lead to follow. It allowed us to identify relevant differences in the behaviour of the patients in 2020, in particular, when considering the vaccination process. Accordingly, we refined our original research question into two sub-questions. To answer these research questions, process mining techniques can provide little help. Process mining and, more in general, process science and process thinking can 770 be a lighthouse in a ocean of data. However, it is difficult to dig deeper by only relying on process mining techniques, given that at the current stage they do not take into account rich perspectives surrounding the process behaviour. In fact, to the best of our knowledge, there are no reliable and effective process mining techniques -in the area of automated process discovery [3] and process variant analysis [4] -that give a global 775 picture of the process, taking into account all the additional data recorded in the event attributes available in the event log. For instance, to answer our research questions, To continue with our analysis, we extracted all the data regarding vaccination events (activity G) from the original dataset (GP16-20 event log), and we analysed the different vaccine types and the immunity they provide, Table 9 shows a mapping be-790 tween the vaccines and the labels we will use to simplify the presentation of the data. Figure 6a shows the absolute number of vaccines we observed in 2020, grouped by the provided immunisation (see Table 9 ). Figure 6b to 6e report the change in the absolute 795 number of vaccines observed in 2020, when compared to the past four years. We 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 following observation. Observation 7. In 2020, there was a surge of influenza (V8) and pneumococcus (V13) vaccinations (see Figure 6 , vaccine V8 and V13), predominant in adults and elderly people, and in contrast with a decrease of these vaccinations for young people (see In the next section, we explore more in depth these implications from a medical point of view. In this section, we described how we have analysed the data and the observations we could draw from it. To analyse the data, we followed a well-known methodology, PM 2 [46] . We note that the latter was designed to be applied in a business-context. However, we argue that this does not pose a threat to its applicability in healthcare, 825 in fact, we were able to adhere to its stages from start to end, with the exception of omitting the execution of the process improvement stage, since it was out of the scope of this study. 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 As broadly discussed in this section, the data we had access to, i.e., the Patron dataset, was far from ready for a process mining analysis. Not only the data was het-830 erogeneous both within a given year and across different years, but its structural quality did not meet the minimum requirements for applying and running state-of-the-art process mining tools, for instance, the timestamps had a day-granularity rather than hour and minute. To overcome the data quality issues, we proposed and applied two novel solutions. While these solutions allowed us to proceed further with our analysis, we 835 recognise that different solution may have been designed which may have, in theory, led to different outcomes. Hence, the observations and the justifications we provide in this study should be interpreted in light of the procedures we applied from the data cleaning and filtering to the process mining analysis. To execute the process mining analysis, we relied on a subset of the existing state-840 of-the-art process mining techniques for automated process discovery and variant analysis, which we selected according to the findings of the most recent literature reviews. While, in theory, applying other techniques may have yielded different or better results, we recall that the process mining techniques we used were the latest and the most reliable. 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 The observations reported in this study cannot be generalised to Australia, nor the 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 state of Victoria. However, we note that the analysed data captured the behaviour of approximately 400 thousand patients (per year), which account for almost 6% of the entire population of the state of Victoria -a remarkable percentage, especially when we consider that not the whole population regularly visit GP clinics. While the ob-850 servations reported in this study are derived from the data and, hence, objective, their analysis and our discussion represent our interpretation. We note that when providing an explanation for a specific observation we considered findings of other similar studies and the experience of two domain experts who co-authored this study (Dr Capurro and Dr Manski-Nankervis). In theory, alternative interpretations for some of our obser-855 vations may be possible but, to the best of our knowledge, the one we provided in this study are the most reasonable and realistic. The study presented here represents the first use of process mining techniques to analyze the impact of the COVID-19 pandemic in health services utilization patterns 860 in primary care. Using a combination of process mining techniques we were able to highlight several relevant changes in health services utilization patterns associated with the disruptions seen in 2020. In addition to these, we were able to highlight some limitations of the process mining tools available today, in particular, when applying them to analyze healthcare process data. Overall, we observed a widespread reduction of GP activities during the period included in our study, when compared to the same period in the four preceding years, concordant with what has been reported in other countries [62] . It is expected that such a reduction of GP activities led to a reduction of specialists visits -given that the Australian is a referral-based system. The consequences of such additional potential 870 reduction of healthcare activities remain to be seen. From the process perspective, in such a situation, we would have expected a reduction in the number of distinct healthcare process execution, instead, the degree of variety of process behaviour remained almost unchanged during this period -with 98% of distinct healthcare process executions observed only between 1 and 5 times. 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 One activity that showed a different behavior were drug prescriptions. We observed an increase in drug refill prescriptions, with peaks in March, April, July and September. These peaks are associated with periods immediately before lock-downs and might represent overstocking of chronic medications. This observation is in line with what has been observed in Australian national drug prescription databases [63] . From the process mining perspective we faced several challenges related to the 895 problem of analysing a vast amount of process execution data. To the best of our knowledge, the event log analysed in this study represents the largest real-life event log used for automated process discovery and process variant analysis, especially, in the healthcare context. We showed that traditional process mining tools present some limitations when attempting to analyze processes with high behavioural variability. The first challenge consisted of imprecise timestamps, since the time granularity was limited to day-level, a recurrent problem in the healthcare context that yet has to be solved. In our case, we relied on clinical knowledge to address this issue by defining a sequence of clinically meaningful activities as a tie-breaker for activities that had identical timestamps. The second challenge involved the identification of start and end events for a pro- 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 cess that is, by nature, unbounded. Once again, we relied on domain expertise to overcome this problem and we presented a generalisation of our solution, suitable for various contexts, in the algorithm described in Section 3. The third challenge was the amount of data itself, and its high behavioral variabil-910 ity, which disarmed state-of-the-art process mining techniques for automated process discovery and process variant analysis. Although the scope of this study was not to devise novel variants of these techniques to deal with such type of data, we highlighted possible directions for future research addressing the improvement of these techniques. Lastly, this study reminds us that state-of-the-art process mining techniques (such 915 as those we used [21, 23, 22, 27, 28, 29] ) do not yet automatically analyse the event log information that is not related to the process behaviour and control-flow (e.g., patient age, medications, etc). This is a straightforward consequence of the existing process mining algorithms design, which do not take into account the additional information even when available. While first steps have been made towards the next generation of 920 process mining algorithms [64] , robust solutions have yet to be proposed. The existing limitation requires process analysts to integrate the process mining analysis with a data mining analysis. While this problem could be solved by further analysing the data from a different perspective, we call for future analysis methodologies and tools that automatically integrate both process and data perspectives. 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 influenza and pneumococcus vaccinations -in contrast with research findings from different geographical areas [13, 15, 14] . The size of the data-set under analysis -counting 31-million events -and the variability of the observed process behavior were unique, and the challenges we faced and overcame during the process mining analysis clearly highlighted the need for improving existing process mining techniques, drawing directions for future work. In particular, future process discovery techniques should integrate in the discovered pro-940 cess models also data surrounding the process behaviour and its control flow. In the healthcare context, this data is the information capturing a patient profile (e.g., age, gender, etc) and their medical procedure (e.g., type of vaccination or prescribed medication). Furthermore, existing process mining techniques are not tailored to deal with large amount of data that captures highly variable behaviour. Future research should 945 consider the design of methods that can automatically filter process execution data to detect and extract the most relevant/interesting process behaviour (not necessarily the most frequent) by analysing the outputs of a set of process mining techniques (e.g., a combination of automated process discovery and process variant analysis). Lastly, as process mining applicability in the healthcare context gains momentum, novel process 950 mining techniques should be tailored for such a context and leverage domain expertise to increase their effectiveness . 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 Fundamentals of business 955 process management Automated discovery of process models from event logs: Review and benchmark Business process variant anal-960 ysis: Survey and classification, Knowledge-Based Systems Conformance checking: a state-ofthe-art literature review Survey and cross-965 benchmark comparison of remaining time prediction methods in business process monitoring Process mining in healthcare: A literature review Leveraging semantic labels for multi-level abstraction in medical process mining and trace comparison Discovering role interaction models in the emergency room using 975 process mining A data-driven framework of typical treatment process extraction and evaluation An approach to automatic process deviation detection in a time-critical clinical process Recommendations for enhancing the usability and understandability of process mining in healthcare Effects of the covid-19 pandemic on routine pediatric vaccine ordering and administration-united states, 2020, MMWR. Morbidity and mortality weekly report 69 The impact of the 990 covid-19 pandemic on immunization campaigns and programs: a systematic review Impact of 995 covid-19-related disruptions to measles, meningococcal a, and yellow fever vaccination in 10 countries Process mining in health-1000 care: Data challenges when answering frequently posed questions, in: Process Support and Knowledge Representation in Health Care Process-aware information systems: bridging people and software through process technology Systematic mapping of process mining studies in healthcare Process mining in healthcare: a systematic review Flexible heuristics miner (FHM), in: Computational Intelligence and Data Mining (CIDM), 2011 IEEE Symposium on Fodina: a robust and flexible heuristic process discovery technique Split miner: automated discovery of accurate and simple business process models from event logs Discovering block-structured process models from event logs containing infrequent behaviour Application of process mining in healthcare-a case study in a dutch hospital Conformance checking Declarative process mining in healthcare Process variant comparison: using event logs to detect differences in behavior and business rules Business process variant analysis based on mutual fingerprints of event logs Detection of statistically significant differences between process variants through declarative rules Investigating clinical care pathways correlated with outcomes Measuring patient flow variations: A cross-organisational process mining approach Process mining for 1050 clinical processes: a comparative analysis of four australian hospitals Covid-19 and stroke-a global world stroke organization perspective Decrease in stroke diagnoses 1055 during the covid-19 pandemic: Where did all our stroke patients go? Covid-19 pandemic and the reduction in st-elevation myocardial infarction admissions Impact of the covid-19 pandemic on cancer care: A global collaborative study Access to cancer surgery in a universal health care system during the covid-19 pandemic Impact of the covid-19 pandemic on lung cancer screening program and subsequent lung cancer Impact of covid-19 pandemic on colorectal cancer screening program Review of the impact of covid-19 on medical services and pro-1075 cedures in australia utilising mbs data: Skin, breast and colorectal cancers, and telehealth services Early estimates of the indirect effects of the covid-19 pandemic on maternal and child mortality in low-income and middle-1080 income countries: a modelling study Special feature: immunization and covid-19 Data for Decisions and the Patron Program Gathering data for decisions: best 1090 practice use of primary care electronic records for research Pm 2 : a process mining project methodology Process Mining -Discovery, Conformance and Enhancement of Business Processes Optimization framework for dfg-based automated process discovery approaches, Software and Systems Modeling Multi-perspective comparison of business process variants based on event logs Event log imperfection patterns for process mining: Towards a systematic approach to cleaning event logs Wanna improve process mining results?, in: 1110 2013 IEEE (CIDM) Improving documentation by repairing event logs Cleaning structured event logs: A graph 1115 repair approach Cleaning timestamps with temporal constraints Automatic repair of same-timestamp errors in business process event logs Metaheuristic optimization for automated business process discovery Measuring fitness and precision of automatically discovered process models: 1125 A principled and scalable approach Filtering out infrequent behavior from business process event logs Discovering more precise process models from event logs by filtering out chaotic activities Improving process discovery results by filtering outliers using conditional behavioural probabilities Reduced in-person and increased telehealth outpatient visits during the covid-19 pandemic Increased dispensing of prescription medications in australia early in the covid-19 pandemic Cocomot: Conformance checking of multi-perspective processes via smt