1 Introduction
Clinical pathways, also known as disease pathways or integrated care pathways, are care plans outlining standardized processes in the delivery of care for a specific cohort of patients (Campbell et al. 1998). They prescribe a sequence of steps for managing a clinical process with the aim of optimizing patient or populationlevel outcomes such as survival, cost, and wait times. Clinical pathways are widely used across a range of medical domains including cancer (Evans et al. 2013, Schmidt et al. 2018, Ling et al. 2018), mental health (Samokhvalov et al. 2018), and surgical recovery (Van Zelm et al. 2018). With applicability to screening, diagnostic, and treatment processes, they have been shown to be effective in improving patient survival and satisfaction, wait times, inhospital complications, and cost of care (Rotter et al. 2012, Panella et al. 2003, Schmidt et al. 2018). Given the importance of clinical pathways in promoting best practices and standardizing or streamlining care, there is significant interest in developing quantitative metrics to measure the concordance of patienttraversed pathways against the recommended pathways, which we will refer to as reference pathways. Such metrics may assist in monitoring variations in the health system and across the care continuum, identifying bottlenecks or changes in the delivery of care, and ultimately providing datadriven evidence to inform health policy decisions. Understanding the root of any discordance is necessary to fully realize the aforementioned benefits of clinical pathways.
Existing methodologies for measuring similarity between pathways have primarily leveraged process mining, a family of techniques that connect business process management with data mining algorithms (Bose and van der Aalst 2010, Adriansyah et al. 2011, Yang et al. 2017, 2016). General similarity metrics are typically based on edit distance algorithms (van de Klundert et al. 2010, Forestier et al. 2012, Yan et al. 2018), which compare two processes represented as a sequence of distinct activities by counting the minimum number of operations needed to transform one into the other. This type of distance measurement can be weighted, where certain activity deviations are more heavily penalized, or unweighted. Unweighted approaches implicitly assume all deviations are equally disadvantageous, which is limiting since different types of discordance will have differential impact on clinical outcomes. On the other hand, weighted approaches to date have relied on domain knowledge and expert input to determine the weights, which is subjective and sensitive to misspecification. The existing approaches in the literature are either unweighted or use subjectively derived weights. Given these limitations, our focus is on developing a weighted concordance metric with objectively derived or datadriven weights.
There are two sources of information that can be used to derive weights for a concordance metric: reference pathways and data from realworld patientrelated activity. We argue that reference pathways should be the primary information source since they are derived from best available medical evidence and are designed for specific clinical outcomes. Furthermore, concordance measurement in a health system is ideally optimized for a hierarchy of outcomes, with clinical outcomes usually given priority. Thus, it makes sense for reference pathways to be used directly to derive weights for the concordance metric. Realworld patient activity, in the form of individual patient pathways, can serve as an important secondary data source to fine tune the weights. However, it should not be used as a primary source to determine the weights since there can be significant variation between individual patient pathways, and outcomes typically depend on many other contributing and confounding factors besides pathway concordance. Patients with good outcomes may have undertaken activities that are not recommended in the reference pathways or that should not be associated with better outcomes. They may also have missed important activities. Overall, the weights should maximize concordance when a reference pathway is measured against itself while scoring discordant pathways lower. The problem of finding such weights can thus be posed naturally as an inverse optimization problem.
In this paper, our goal is to use inverse optimization to develop a datadriven concordance metric for measuring clinical pathway concordance that has a validated, statistically significant association with patient survival. Such a metric can help pinpoint discordant activity in a complex health system, from geographical regions to individual providers, and therefore support policy makers in crafting informed, datadriven decisions about health system interventions. To support such datadriven decision making, we first need a rigorous approach to measuring concordance.
We model a patient’s journey through the healthcare system as a walk in a directed graph in which each node represents a possible activity that a patient can undertake. Patients incur “costs” by visiting nodes and traversing arcs through the network, which allows us to model both the costs associated with undertaking (or missing) certain activities, as well as sequencing costs for progressing from one activity to the next. A reference pathway can be considered an optimal solution of an appropriately formulated minimum cost network flow problem, a shortest path problem in particular. Given that the “forward” problem is a network flow problem, we can formulate the “inverse” problem tractably using linear programming duality
(Ahuja and Orlin 2001). The goal of the inverse optimization model is to find arc costs such that all reference pathways are optimal for the resulting network. If there does not exist a cost vector that makes all reference pathways simultaneously optimal, then the costs should be such that a measure of aggregate error (e.g., optimality gap) is minimized. Once the arc costs are determined, any input pathway (reference or patienttraversed) can be scored based on the cost of the associated directed walk through the network. In other words, these arc costs form the weights used in the concordance metric.Modern inverse optimization methodologies aim to impute unknown parameters of an optimization problem to make observed decisions minimally suboptimal (Troutt et al. 2006, Keshavarz et al. 2011, Chan et al. 2014, Bertsimas et al. 2015, Aswani et al. 2018, Esfahani et al. 2018, Chan et al. 2019, Babier et al. 2019). Inverse network flow problems are wellstudied in the literature, but primarily focus on a single, noiseless observed decision for which it is possible to determine parameters that make the decision exactly optimal (Burton and Toint 1992, Xu and Zhang 1995, Zhang et al. 1995, Zhang and Ma 1996, Yang et al. 1997, Burton et al. 1997, Zhang and Cai 1998, Ahuja and Orlin 2001, 2002). Inverse network flow problems with multiple noisy observations have received limited attention to date (Faragó et al. 2003, Zhao et al. 2015). We extend these papers by introducing a general inverse optimization framework for standard form linear optimization problems that explicitly accounts for the lower dimensionality of the feasible region due to the equality constraints corresponding to flow balance. This dimensionality issue also requires modifying how goodnessoffit is measured in inverse optimization, which to date has only considered problems with full dimensional feasible regions (Chan et al. 2019). Finally, unlike many previous papers, we do not restrict the sign of the cost vector.
For a practical application of our inverse optimization framework, we study clinical pathway concordance for stage III colon cancer pathways developed by Cancer Care Ontario. Cancer Care Ontario is the Ontario government’s principal advisor on cancer care. It directs and oversees healthcare funds for hospitals and other cancer care providers, enabling them to deliver highquality, timely services and improved access to care. As such, Cancer Care Ontario’s responsibilities include developing guidelines and standards related to the delivery of cancer care, as well as monitoring the performance of the cancer system in Ontario, Canada. Colon cancer is one of the most commonly diagnosed cancers and a leading cause of death from cancer worldwide for both men and women (World Cancer Research Fund 2019). In our numerical experiments, we demonstrate how inverse optimization can be used to learn a cost vector from reference and patient pathways in an appropriately formulated clinical activity network. We develop both a goodnessoffit metric to measure the quality of the inverse optimization solution, as well as a concordance metric to measure the concordance of individual patient pathways. Finally, we provide a rigorous survival analysis to demonstrate the association of our concordance metric with survival, using a cohort of stage III colon cancer patients. Although we focus on colon cancer specifically in this paper, our inverse optimizationbased approach is generalizable to other patient cohorts and even beyond the healthcare context to general process management applications for which datadriven concordance measurement is a critical need.
Our main contributions in this paper are as follows:

We propose the first inverse optimizationbased approach for measuring pathway concordance in any problem context. Methodologically, our approach differs from previous inverse optimization models with the inclusion of a mixture of both “good” and “bad” data points; that is, points that the optimal cost vector should fit well and points that it should not fit well, respectively. Our approach also considers the novel setup where there is a primary (reference pathways) and secondary (patient pathways) dataset, which should be given differential consideration when imputing the cost vector.

We develop two metrics for goodnessoffit assessment. The first metric extends a previously developed metric for insample goodnessoffit measurement between the decision model, the input data, and the imputed parameters. This metric is tailored to our particular inverse optimization model and considers the lower dimensionality of the feasible region in standard form linear optimization problems such as network flow problems. The second metric, used for outofsample measurement, uses the inversely optimized cost vector to calculate a concordance score for a given flow vector (patient pathway). Both metrics are simple to compute and easily interpretable, with intuitive mathematical properties such as taking values between 0 and 1, with 1 representing perfect fit or concordance.

We provide an indepth case study applying our inverse optimizationbased models and metrics to real patient data from a cohort of stage III colon cancer patients. We show through a rigorous survival analysis that there is a statistically significant association between the concordance score of patienttraversed pathways and survival, validating the meaningfulness of our concordance metric. We further quantify the value of our datadriven approach by demonstrating that the patient datadriven concordance metric is much more strongly associated with survival than a metric developed using only reference pathways and without access to the patient data.
All proofs are included in the Appendix.
2 Motivating a GraphBased Inverse Optimization Approach
In this paper, we leverage network optimization and the corresponding graphbased representation of the problem to measure pathway concordance. To motivate our approach, we briefly highlight the differences between edit distancebased methods (Navarro 2001)
and our proposed approach. In particular, we show that our approach provides more granularity and modeling flexibility in capturing different types of discordance, and provides a framework that facilitates datadriven estimation of weights.
First, consider unweighted edit distance algorithms, which measure distance as the minimum number of distortions or edits required to transform a source pattern into a target pattern (van de Klundert et al. 2010, Forestier et al. 2012, Yan et al. 2018, Riesen 2015). The edit operations are defined on the activities in the pathway and usually include insertions, deletions, substitutions, and transpositions, although some algorithms may count a transposition as two substitutions. In contrast, our graphbased approach measures distance using the difference in the total cost along the arcs traversed by the two pathways. Note that edit distance has been extended to graph structures where the edit operations are applied to nodes and arcs (Gao et al. 2010). In the context of clinical pathways, the entire graph is simply a path representing the clinical events of interest sequentially in time. A key difference is that graph edit distance compares two pathways represented as two distinct graphs, whereas in our approach the two pathways are distinct walks in a single, predefined graph.
Our approach improves on edit distance algorithms by better characterizing pathway deviations and optimally estimating weights for those deviations. We refer to arcs that are in a reference pathway as “concordant” if they connect the desired activities in the order indicated according to clinical guidelines. In contrast, “discordant” arcs are those that do not connect activities contained in the reference pathways, or connect such activities out of sequence. The concordance of an observed pathway is then a function of the discordant arcs present and the concordant arcs absent.
To illustrate the ability of our algorithm to better characterize pathway variation, we present concrete examples shown in Table 1. Consider a network with an artificial start node and end node that bookend all pathways and let be a reference pathway represented by the sequence of arcs , , and . Now consider five example pathways listed in Table 1, all of which have the same edit distance of one from . We can see that despite having the same edit distance, these pathways have different numbers of concordant and discordant arcs. Concordance measurement based on absence of concordant arcs and presence of discordant arcs implies that each pathway can have a different graphbased distance from the reference pathway. Furthermore, when considering pathways with an edit distance of two or more from the reference pathway, differences with graphbased distance measurement are magnified. For example, two adjacent deletions are different from two nonadjacent deletions in terms of the number of concordant and discordant arcs.
Pathway  Edit operation  Missing concordant arcs  Extra discordant arcs  Pathway representation 

ABC    0  0  
AABC  Insertion (duplicated activity)  0  1  
XABC  Insertion (discordant activity)  1  2  
BC  Deletion  2  1  
XBC  Substitution  2  2  
BAC  Transposition  3  3 
Notes. Missing concordant arcs are dashed. Extra discordant arcs are red.
The preceding discussion assumes that all arcs and operations (in the edit distance setting) are unweighted. If we now consider a weighted distance metric, an issue arises where the cost of edits would need to be determined a priori and exogenously from our reference pathways. In weighted edit distance, weights or costs are associated with each type of operation. The distance between pathways then becomes the minimum cost of edit operations required to transform one into the other. Determining the costs of each operation based on a set of reference and nonreference pathways is complicated by the dependency of the edits on the costs. Expectation maximization can be used, but this approach does not guarantee convergence to a global optimum
(Neuhaus and Bunke 2004, Ristad and Yianilos 1998). In contrast, our inverse optimization approach avoids this dependency and can determine optimal costs for each arc by assuming pathways are optimal solutions to an appropriately formulated optimization problem. The ability to optimally infer arc costs or weights, together with finer characterization of pathway variation using concordant and discordant arcs, are the main strengths of our inverse optimization approach.3 Inverse Optimization for Pathway Concordance Measurement
In this section, we first present a model of the clinical activity network representing all the activities that a cancer patient can undertake during his or her cancer care journey. A patient journey is modeled as a walk through the network. We assume that each arc incurs a cost when a patient takes that arc and that the reference pathways are optimal solutions with respect to an appropriate optimization problem defined on this network. Our goal is to find the costs on each arc and to use it to measure the concordance of a given pathway to the reference pathways.
3.1 The Clinical Activity Network
Let be a set of nodes and let be a set of directed arcs in the network. Nodes represent clinical activities that can be undertaken by a patient. To model the cost of undertaking a specific activity , we split each activity node into two nodes, and , representing the start and end, respectively, of activity , and then add an arc with cost . The remaining arcs connect the end node of some activity with the start node of another activity . Each arc has an associated cost , representing the cost of undertaking activity immediately after activity . The set also includes arcs of the form , which allows the network to model a patient journey that involves consecutive repetitions of activity . Finally, we add two artificial nodes representing the common start and end of each patient journey. Thus, is equal to two times the number of clinical activities plus two.
A reference pathway is represented as the solution of a shortest path problem on this network. We write this shortest path problem as a standard form linear optimization problem representing the corresponding network flow model: . We denote this model the “forward” model and will refer to it as . The decision vector describes the flow on each arc. The matrix is the nodearc incidence matrix, but without the th row, which corresponds to the flow balance constraint of the artificial end node. Thus, it is full rank. The parameter is a vector of zeros except with a at the artificial start node.
A patient journey will be represented as a walk on this network. Note that it may not necessarily be a path since patients may visit an activity or traverse an arc several times. Thus, a patient journey through the network can still be represented by a flow vector , except that the value of may be greater than one, representing the number of times arc is traversed in the corresponding walk. In keeping with the medical terminology, we will refer to a patient journey as a patient “pathway”, which in general should be thought of as a walk through the network. We will reserve the use of the phrase “path” for the mathematical concept of a walk without repeated nodes.
Note that with the current setup, costs are “pathindependent”. That is, the cost of traversing arc is equal to regardless of how the walk arrives at node . To model pathdependent costs, where the cost of completing activity immediately following activity depends on what activities were completed prior to , we could modify the nodes in the network from representing a single activity to representing the sequence of activities from the artificial start node to the current activity. That is, imagine a network organized in “layers”, where the first layer of nodes comprises all activities that can be reached from the start node in one step, the second layer of nodes comprises all sequences of activities that can be reached in two steps, etc. While the costdependent approach provides richer modeling possibilities, the downside is an exponential explosion in the number of nodes in the network. Thus, we proceed with the network as currently defined (i.e., with pathindependent costs) to maintain a parsimonious model, and instead integrate pathdependent behavior into the constraints of the inverse optimization model, which is described next.
3.2 Inverse Optimization Model
Our goal is to identify a cost vector that minimizes the aggregate suboptimality of given reference pathways, assuming they represent shortest paths. We first present a basic inverse optimization model that provides an exact formulation of this problem. Then, we introduce a second formulation that uses patient pathways to refine the optimal cost vector returned by the first model, in order to improve the association between pathway concordance and survival.
3.2.1 A Model Using Only Reference Pathways.
Let denote the set of reference pathways, each assumed to be feasible flow vectors for the forward problem. We follow a standard approach to formulating the inverse optimization problem using duality of linear optimization. Let be the dual vector associated with flow balance constraints. The following formulation is an inverse optimization model that minimizes the sum of squared absolute duality gaps induced by the cost vector and the reference pathways:
(1)  
subject to  
The first constraint represents dual feasibility. The second constraint defines the duality gap for each reference pathway as the difference between the primal and dual objective values. Each duality gap variable is nonnegative since the reference pathways are feasible for the forward problem. The third constraint is a normalization constraint to ensure the cost vector is not trivial (i.e., nonzero). The fourth constraint is needed to ensure that lies in the lower dimensional space defined by the constraints . Since lies in the space defined by the intersection of the equality constraints (), it must be orthogonal to the vectors defining those constraints (). Without the constraint , an optimal solution (i.e., cost vector) to formulation (1) could be orthogonal to the entire forward feasible region, which would render uninformative since every feasible flow vector would be optimal. To our knowledge, the issue of finding an optimal cost vector via inverse optimization for a lower dimensional feasible region (which is characteristic of any inverse network flow problem) has not been studied in the inverse optimization literature previously. Note that the constraint means that the optimal cost vector will be a circulation.
To model additional applicationspecific considerations, such as userspecified rankings between certain activities or preferences for how a patient should begin or end her pathway, we will add a set of constraints on the cost vector, , to formulation (1). Specific examples of such constraints are provided in Section 5.
First, we prove that formulation (1) has an optimal solution for any nontrivial network.
has an optimal solution if and only if there are at least two distinct paths from the start node to the end node.
Given Lemma 3.2.1, we assume going forward that the network is nontrivial (i.e., more than a single path) and always has an optimal solution. Next, we present an important property of an optimal solution to .
Let be an optimal solution to formulation (1). Then, the network with arc costs does not contain any directed negative cost cycles.
This result has an intuitive realworld interpretation. It says that a patient traversing the cancer care network cannot accrue a positive benefit through repeating any sequence of activities. This will certainly be the case if the outcome is measured in terms of monetary cost of the activities. But even viewed through the lens of survival, it makes sense that such negative cost cycles do not exist. Otherwise, there would be a sequence of activities that guarantees improved survival, which is not realistic. We also note that a finite optimal value for (1) is important for the goodnessoffit and concordance metrics, which we define later.
Formulation (1) is nonconvex because of the normalization constraint. However, it can be solved using polyhedral decomposition: solving quadratic optimization problems where the normalization constraint is replaced with for all and with one additional constraint on one being set to either or :
(2)  
subject to  
Furthermore, if in a particular problem context it is clear which arc should be either the most heavily rewarded (set ) or heavily penalized (set ), then only one of these problems needs to be solved. As we demonstrate later, this is the case for our clinical problem.
3.2.2 Refining the Solution using Patient Data.
Let denote a dataset of patient pathways with good clinical outcomes (i.e., survived) and denote a dataset of patient pathways with poor clinical outcomes (i.e., died), all assumed to be feasible flow vectors for the forward problem. If formulation (1) has multiple optimal solutions, we can use these two datasets to determine a cost vector that not only maximizes fit with the reference pathways, but also provides “separation” between and . Once formulation (1) is solved and an optimal duality gap vector, , is generated, we use as input into the following inverse optimization model:
(3)  
subject to  
Since data from represents “good” data, we treat it like input to traditional inverse optimization models (and like the reference pathways), where the goal is to minimize the aggregate suboptimality of these data points. However, since consists of “bad” data points, we do not want to learn a cost vector that fits this data well. Instead, the cost vector should generate data points as far away as possible from those in when solving the forward problem, which is why in formulation (3) the objective maximizes suboptimality with respect to . The weight simply scales the two subobjectives to account for the difference in the number of data points in the two groups. The second constraint forces the cost vector to achieve the optimal duality gap found in model (1). In other words, all feasible solutions of (3) are optimal solutions to (1). The third and fourth constraints define the duality gaps with respect to the pathways in and , respectively. Again, the duality gaps are nonnegative since all patient pathways are feasible for the forward problem. The remaining constraints are as defined in (1). Formulation (3) can be solved using the same polyhedral decomposition technique described previously. The difference is that each subproblem is a linear optimization problem, since the objective function is now linear. A linear objective was chosen for formulation (3) to facilitate the combination of duality gaps that needed to be minimized and maximized in the same objective.
Next, we establish that formulation (3) will always return an optimal cost vector.
has an optimal solution.
In our numerical results, we implement two slightly different versions of formulation (3). In Section 6, we visualize the solution of our inverse optimization model (i.e., the cost vector) and compute goodnessoffit using all of our patient data to form the sets and . These results are thus “insample”. However, the patient concordance scores in that section and in the survival analysis in Section 7 are computed “outofsample”. That is, we use only a subset of the full patient dataset to form and in model (3), and the resulting optimal cost vector is used to calculate a concordance score for the patients that were not used to determine . These outofsample concordance scores then constitute an independent variable in the survival analysis model.
4 GoodnessofFit and Concordance Metrics
In this section, we develop two fitness metrics. The first measures the goodness of fit of the optimal cost vector generated by the inverse optimization model against the input data. The second metric measures the concordance of a given patient pathway against the reference pathways, using the optimal cost vector generated by the inverse optimization model. It should be applied to patient pathways that have yet to be observed and that are not in the “training data” that was used to generate optimal cost vector.
4.1 GoodnessofFit of the Cost Vector
We extend a previously developed goodnessoffit metric for inverse optimization, the coefficient of complementarity (Chan et al. 2019), so that it is appropriate for measuring the fitness of cost vectors produced by model (3). The previous metric required the forward feasible region to be fulldimensional, which is not the case here due to the flow balance equality constraints. Our metric is like
in linear regression in that it measures the fit between the input data and the estimated model, and possesses analogous mathematical properties. First, we present a general form of our goodnessoffit metric that is applicable to any problem, regardless of whether the forward feasible region is fulldimensional, and examine its theoretical properties. Then, we introduce a modified form that is appropriately tailored to the lowerdimensional space in which the forward feasible region resides.
Since our primary goal is to find a cost vector that fits the reference pathways, and since model (3) simply searches among the optimal cost vectors for one that separates patients with good and bad outcomes, we focus on measuring modeldata fit with respect to the reference pathways and not the patient data.
Let be an optimal cost vector for (3) and be a set of cost vectors that are feasible for (1). Let be an optimal solution to and be defined similarly with respect to . We define the goodness of fit of with respect to the reference pathways as
(4) 
Note that and can be used interchangeably since they represent the optimal primal and dual objective values, respectively, which are equal. Therefore, . Similarly, each term in the denominator is equal to a duality gap for a reference pathway with respect to the cost vector .
Next, we show that has attractive mathematical properties that facilitate its use and interpretation. Recall that an optimal for formulation (3) is also optimal for (1). Thus, although Theorem 4.1 focuses on formulation (1), it is equally applicable to cost vectors generated from model (3). Some additional notation is needed before proceeding. Let be a set of arcs from the network, . We define
where is an optimal solution for , defined as with additional constraints for .

Optimality: is maximized by an optimal solution to .

Boundedness: .

Monotonicity: if .
Given the reference pathways and an optimal solution to (1) (or (3)), measures the degree to which the reference pathways are shortest paths with respect to . If all reference pathways are shortest paths with respect to , then .
When the forward feasible region is fulldimensional, an appropriate choice for the baseline cost vectors consist of the vectors normal to the constraints defining the feasible region, i.e., the rows of the coefficient matrix (Chan et al. 2019). Thus, where it is straightforward to work directly in the lowerdimensional space, we should modify the computation of to avoid having to make a subjective choice for . In particular, we can eliminate the equality constraints by variable substitution (details in Appendix 10) and compute the following “fulldimensional” version of in dimension
(5) 
where and defines the th inequality constraint defining the feasible region in the lowerdimensional space (i.e., , ). Importantly, the properties proved in Theorem 4.1 hold for (Chan et al. 2019). In our numerical results, we project the reference pathways and forward feasible region to the lowerdimensional space and calculate when measuring fitness between the optimal cost vector and the reference pathways.
4.2 Concordance of a Patient Pathway
Using the cost vector generated by inverse optimization, our concordance metric measures the cost of a patient pathway against the cost of a shortest path with respect to as
(6) 
where . In other words, is the cost of a longest walk, which can be determined via dynamic programming (see Appendix 11). Note that we use the cost of the shortest path, , instead of the cost of a reference pathway when measuring concordance. Since a reference pathway may not actually be a shortest path under , and since there may be multiple reference pathways with different costs under , we define the concordance metric using to have a consistent and objective baseline. Of course, if the inverse optimization model finds a cost vector that perfectly fits all of the reference pathways (i.e., ), then the cost of all reference pathways will be equal to , and will be measuring concordance directly against the cost of the reference pathways.
Next, we show that has the attractive property of being a unitless metric between 0 and 1, with 1 representing a perfectly concordant pathway (i.e., is a shortest path with respect to ) and 0 representing a perfectly discordant pathway (i.e., is a longest walk with respect to among all walks of at most the same number of steps as ).
Given any feasible solution to the forward problem:

,

if and only if ,

if and only if .
We end this section with a brief discussion on . The denominator of normalizes the cost difference between the given pathway and an ideal pathway using the cost of the longest walk with length up to the length of . Normalizing in this fashion implies that a longer walk could be, paradoxically, more concordant than a shorter walk, something that would not occur if the normalization term was a constant. However, we choose to normalize using for several reasons. First, this normalization allows the computation of to depend only on the given pathway (and ) and not any other pathways. Normalizing by a fixed amount, , while maintaining the property that for all would require that be longer than all patient pathways that would be considered by , an amount that could not be known a priori. Second, this approach aligns with common normalization approaches for edit distance metrics that are based on string length (Yujian and Bo 2007). Later, in our numerical results, we show that using
results in an approximately normal distribution of
for our dataset.5 Application to Stage III Colon Cancer
In this section, we demonstrate the application of our inverse optimization framework to measure clinical pathway concordance, focusing on stage III colon cancer.
5.1 Problem Background
5.1.1 Disease Pathway Management at Cancer Care Ontario.
The Disease Pathway Management (DPM) program at Cancer Care Ontario takes a systems view of setting priorities for cancer control, planning cancer services and improving the quality of cancer care in Ontario. DPM views patient journeys across the entire cancer continuum as part of an integrated process, as opposed to evaluating individual points of care in isolation. Clinicians are engaged to establish evidencebased best practices for the sequence of activities patients should take through the care network. These best practices are captured in pathway maps (Cancer Care Ontario 2019b), which reflect clinical guidelines for the care that patients should receive to optimize clinical outcomes. Pathway maps cover points of care spanning disease prevention, screening, diagnosis, treatment, recovery and endoflife palliative care.
5.1.2 Stage III Colon Cancer and Pathway Map.
Colon cancer is a malignant tumor of the colon. In 2018, roughly 180,000 people were diagnosed with colorectal cancer in North America and there were 64,000 deaths due to this cancer (World Health Organization 2019); colon cancer constitutes approximately 70% of the colorectal cancer cases (Colorectal Cancer Alliance 2019). Stage III colon cancer is characterized by metastasis of the cancer to the lymph nodes near the colon or by direct invasion into nearby organs. Chemotherapy is recommended as adjuvant therapy for surgically resected stage III disease. Twelve cycles of chemotherapy constitutes a complete treatment course.
The stage III colon cancer pathway map (Cancer Care Ontario 2019a) spans diagnosis and treatment, and contains the following major categories of activities: 1) clinical consultations, 2) endoscopy, 3) diagnostic imaging, 4) surgery, 5) adjuvant chemotherapy. Clinical consultations include consultations with a gastroenterologist (GI), surgeon and medical oncologist (MO). Diagnostic imaging includes imaging of the abdomen, pelvis and chest with any of the following modalities: computed tomography (CT), magnetic resonance imaging (MRI), ultrasound (US), or xray.
All concordant pathways must start with an endoscopy followed by abdominal imaging, chest imaging, surgical resection of the tumor, a MO consultation and finally chemotherapy. The patient will have had a GI or surgical consultation prior to the endoscopy, and a surgical consultation at some point before the surgical resection. Similarly, a MO consultation must occur before chemotherapy, preferably after the surgical resection. Abdominal imaging may take place as either an abdomen CT scan followed by a pelvis CT scan, or as an abdominal ultrasound. Similarly, chest CT scans or chest xrays constitute chest imaging. Chemotherapy is considered concordant if the patient pathway includes the recommended 12 treatment cycles. All of these different combinations are captured in the pathway map. Enumerating all of the possible concordant pathways from the pathway map results in 108 distinct reference pathways for stage III colon cancer.
5.2 Data
5.2.1 Data Sources.
In Ontario, health care encounters are recorded in administrative databases using standardized methods. Using encrypted personspecific identifiers, we linked multiple sources to select the study population, define patient characteristics, and record activity related to screening, imaging, treatments, emergency department (ED) visits, and death. The data sources used in this study are listed in Appendix 12.
5.2.2 Study Population.
We developed a populationbased cohort of individuals newly diagnosed with stage III colon cancer in 2010 in Ontario. Their clinical activities were recorded from six months prior to cancer diagnosis until either death or censoring at four years following diagnosis. Included patients had surgical resection within one year of diagnosis. Patients were excluded based on the following exclusion criteria: no Ontario health insurance number, residence outside of Ontario, nonincidental cancer diagnosis, any prior history of cancer or identification of a new primary cancer during the followup period, cancers diagnosed based on death records, unknown cancer stage, and unknown tumor histology. The final cohort contained 763 patients, 257 of whom died (from any cause) within the fouryear postdiagnosis observation period.
The dataset is an eventlog for each patient in the cohort describing the activities undertaken and their corresponding times relative to the date of diagnosis. However, not all events in the eventlog are relevant to the diagnostic and treatment phases of care, which is our focus. Therefore we refine each patient’s eventlog by filtering out events earlier than 30 days prior to diagnosis and following one year postdiagnosis, since patients are expected to complete diagnostic and treatment activities within one year. If the patient completed chemotherapy (all 12 cycles) within the year, we remove all events after their last chemotherapy cycle. We also remove activities that have no bearing on concordance. For example, minor procedures like fecal occult blood tests and barium enemas were deemed inconsequential or extremely infrequent and were removed. All data refinements were implemented following consultations with clinicians.
5.3 Network Design
Next, we form a directed graph representing the clinical activity network based on activities from the pathway map, observed patient pathways, and discussions with clinical experts from the DPM team.
In clinical practice, GI, surgery and MO consultations always precede endoscopies, surgical resection and chemotherapy treatment, respectively. We therefore do not include nodes for these consultations in the graph and assume that they occur with the corresponding procedure. However, it is possible for patients to have additional consultations, which are costly for the system and deemed discordant. These discordant consultations, like visits to the emergency department, were quite prevalent in our dataset, so we included EXTRA CONSULT and ED VISIT as nodes in the graph. Another reason to include the latter is because reducing unnecessary visits to the ED is a Cancer Care Ontario priority.
Imaging events are distinct from consultations in that they are not implicit with other events, so we keep them as separate standalone activities. Imaging activities that were deemed functionally equivalent from a concordance perspective were merged into a single node. For example, abdomen MRI and ultrasound were combined into the node ABDOMEN MRI/US. We formed a single node for PELVIS MRI/US and CHEST IMAGING (which encompasses chest Xrays and CT scans) for the same reason. CT is preferred over the other imaging modalities for abdomen and pelvis, so we defined separate nodes for abdomen CT and pelvis CT.
Finally, our discussions with experts revealed that the type of treatment received by the patient should have extra bearing on their concordance score. To model the different possible ways each clinical or patient pathway can end, we include an “outcome layer” with four distinct nodes. Each pathway must exit the graph via one of these four nodes. Patients who receive 12 chemotherapy treatments exit through the node CHEMO COMPLETE. Receiving between 1 and 11 chemotherapy treatment cycles leads to an exit through the node CHEMO PARTIAL. The remaining pathways (no chemotherapy) exit either through the nodes MO CONSULT END (if the patient receives an MO consultation) or RESECTION END (if the patient does not). Recall that our patient cohort only includes patients who had surgical resection. The final network design used in our model is shown in Figure 1.
5.4 Inverse Optimization Model Specification
Since our final network design merged together several activities specified in the Cancer Care Ontario pathway map, the set of 108 reference pathways were reduced to two distinct reference pathways. One reference pathway is shown in Figure 1. The other reference pathway replaces ABDOMEN CT and PELVIS CT with ABDOMEN MRI/US. The network, reference pathways, survived patient pathways and died patient pathways constitute the , and , respectively, in formulations (1) and (3). The value of the th component of is 1 if reference pathway includes arc and 0 otherwise. Patient pathways are mapped to the graph similarly, except that each component of the vector represents the number of times that arc is traversed (which may be greater than one) in the patient pathway. Finally, with respect to , we included three sets of contextspecific constraints on the imputed cost vector: 1) activity ranking constraints, 2) subpath constraints, and 3) penalty constraints.
Activity ranking constraints enforce a relative importance between pairs of activities. For example, clinical experts deemed ABDOMEN CT to be a more important activity than PELVIS CT. To account for this relationship in the model, we introduced a constraint that enforced the cost of completing ABDOMEN CT to be at most the cost of completing PELVIS CT (since it is a minimization problem). Specifically, if represents ABDOMEN CT and represents PELVIS CT, then the constraint would be . There were eight such constraints in total between the nine nonoutcome layer nodes.
Subpath rankings generalize activity ranking constraints to a sequence of nodes. For example, the subpath of ENDOSCOPY to ABDOMEN CT to PELVIS CT to CHEST IMAGING was preferred to the subpath of ENDOSCOPY to ABDOMEN MRI/US to CHEST IMAGING. Therefore, the sum of the arc costs along the first subpath (starting from the end node of ENDOSCOPY and ending with the start node of CHEST IMAGING) was constrained to be at most the sum of the arc costs along the second subpath. There were four such constraints in total. The first is the subpath ranking above. The remaining three rank the subpaths exiting the graph through the four outcome layer nodes, starting each subpath at RESECTION. In particular, the nodes in the outcome layer as depicted in Figure 1 are arranged from top to bottom in order of most preferred to least preferred outcome. For example, if represents RESECTION, represents RESECTION END, represents MO CONSULT END, and represents the final end node, then the subpath ranking between RESECTION END and MO CONSULT END would be . Similar subpath rankings were generated for MO CONSULT END versus CHEMO PARTIAL and CHEMO PARTIAL versus CHEMO COMPLETE.
Finally, since CHEMO COMPLETE is the best outcome node through which to exit the graph, we set the cost of the CHEMO COMPLETE to END arc to , anchoring all the arc costs relative to this one. This constraint is also beneficial from a model tractability perspective. Recall that the inverse optimization model is solved via polyhedral decomposition, which amounts to solving several distinct optimization problems, each with a constraint of the form or . Since we have an applicationspecific justification for setting one of the arc costs so that the constraint is automatically satisfied, we only need to solve a single optimization problem.
6 Inverse Optimization Results
We obtained an optimal cost vector, , by first solving formulation (1) followed by formulation (3) with the inputs described in Section 5.4. This cost vector comprises all arc costs and by definition minimizes the total suboptimality error of the two reference pathways, and , and secondarily minimizes (maximizes) the suboptimality error of the patient pathways from survived (died) patients in (). Figure 2 illustrates the values on the graph. Recall that since satisfies , the sum of the costs on the incoming and outgoing arcs for each node are equal. As a consequence, the associated cost of passing through each node in the outcome layer is doubled, since each node in that layer has only one outgoing arc (i.e., the intraactivity arc followed by the arc from the outcome node to END).
Given the optimal cost vector , we found that . In other words, the duality gap associated with the two reference pathways was zero, indicating perfect modeldata fit with the reference pathways. However, the average error from the patient pathways was positive, indicating that they are not shortest paths in general.
Finally, we calculated the concordance score of each patient pathway using (6) following a 10fold crossvalidation approach. The patient data was randomly partitioned into ten subsets (“folds”) that maintained the same class balance as the entire dataset. For each fold, the inverse optimization model was trained using the two reference pathways and 90% (nine of the subsets) of the patient data. The resulting cost vector was used in to measure the concordance of patient pathways in the outofsample 10% (the remaining subset). This process was carried out for all 10 folds and the distribution of outofsample concordance scores for the entire patient cohort is shown in Figure 3. The fact that this distribution is approximately normal stems from the choice of normalization in . As Figure 4 shows, the distributions of the duality gap and pathway length
are both skewed similarly, resulting in an approximately linear relationship between the two quantities. In other words, the lack of fitness of a given patient pathway is associated with its corresponding length, and dividing one by the other creates the shape of the distribution seen in Figure
3.As a crude analysis of association between concordance and survival, we generated smoothed cumulative distributions of the concordance scores, stratified by survival outcome, in Figure 5. Interestingly, the figure depicts firstorder stochastic dominance, with a clear separation between patients who survived versus those who died. That is, for any given concordance score
, the probability of scoring below
is higher for patients who died. Note that stochastic dominance is used in medical decision making to compare the costs or effectiveness of two or more distinct treatments (Leshno and Levy 2004, Sendi et al. 2003).The results shown in Figure 5, while intuitive, do not provide a rigorous measure of the relationship between concordance and survival. Thus, in the next section, we perform a detailed survival analysis to establish the association between pathway concordance and survival.
7 Validating the Concordance Metric against Survival
7.1 Methods
For each patient in our study population, we obtained their characteristics as outlined in Section 5.2
. Patient characteristics were summarized as counts with proportions for categorical data and means with standard deviations for continuous data. We binned the patient concordance scores into terciles and difference in covariate distributions among concordance score terciles was assessed with a chisquared test for categorical variables and a oneway ANOVA for means and standard deviations. A KaplanMeier estimator was used for graphical presentation of (unadjusted) survival over time by the terciles of concordance scores. Difference in survival probabilities was assessed using the logrank test statistic.
A semiparametric Cox proportional hazards model was implemented to evaluate the association between concordance scores and mortality (Cox 1971). Cox regression provides estimates of association between each variable and the hazard of event, i.e., the instantaneous rate of event. A unit increase in a covariate is multiplicative with respect to the hazard. An estimate of the hazard ratio (HR) above indicates a variable that is positively associated with the risk of mortality, i.e., negatively associated with survival.
A range of known potential confounders and predictors of survival was assessed for inclusion in the model. These variables consist of demographic characteristics (age, sex, rural residency, quintiles of neighbourhood median income, terciles of immigrant population), Charlson comorbidity score (Charlson et al. 1987, Quan et al. 2005), cancer diagnostic and treatment characteristics (screening category, stage, tumor grade, emergent surgery status, length of hospital stay following the surgical treatment), and healthcare utilization a year prior to diagnosis (number of outpatient visits, presence of inpatient admissions, presence of ED visits), which is a proxy for overall patient health. Screening category refers to frequency and timing of blood tests (James et al. 2018). Stage specifies the cancer substage, which denotes severity. The emergent surgery indicator refers to patients who had an urgent inhospital admission when the surgical resection was performed.
We employed preliminary variable selection to control for important confounding effects while excluding covariates that would compromise model efficiency. Variable importance in predicting mortality was evaluated using an ensemble tree method for analysis of rightcensored survival data (Ishwaran et al. 2008).
For the selected model, assumption of a linear relationship between a continuous independent variable and the outcome was evaluated by plotting the residuals against the covariate. If nonlinearity was detected either from the plots of Martingale residuals or Schoenfeld residuals, restricted cubic splines were used to model the continuous predictor (Stone and Koo 1985).
Variance of the estimated association between the concordance score and survival can increase in case of overfitting. Bootstrap model validation with 1000 resamples was performed to evaluate the degree of overfitting and potential overadjustment for confounders (Harrell Jr 2015).
Patients were excluded from the final survival model if they did not have all required covariate values, other than tumor grade, recorded. Patients who died within one year of diagnosis were also excluded, as they did not have the opportunity to complete their treatment journey in accordance with the clinical guidelines.
7.2 Results
A total of 665 individuals met the inclusion criteria for survival analysis. The mean age at diagnosis was 67.8 years (standard deviation = 13.2). Of the included patients, 165 (24.8%) died, with a higher proportion of deaths occurring among those with low concordance scores (38.5%, pvalue ).
Table 2 lists the characteristics of the entire cohort and by concordance terciles. Significant differences between the terciles were observed for age, rural residency, comorbidity score, healthcare utilization, cancer stage, and characteristics related to cancer treatment. Sex, neighbourhood income and immigrant population, screening group, and tumor grade were not significantly different across terciles.
Characteristic 




pvalue  
Death, n (%)  85 (38.5)  49 (22.1)  31 (14.0)  165 (24.8)  <0.0001  
Concordance score, mean (SD)  0.4 (0.1)  0.5 (0.0)  0.7 (0.1)  0.5 (0.2)  <0.0001  
Patient characteristics  
Age at diagnosis, mean (SD)  74.0 (10.7)  65.9 (13.6)  63.4 (12.8)  67.8 (13.2)  <0.0001  
Female, n (%)  121 (54.8)  119 (53.6)  101 (45.5)  341 (51.3)  0.1044  
Rural residency, n (%)  39 (17.6)  22 (9.9)  25 (11.3)  86 (12.9)  0.0348  
Neighbourhood Income Quintile, n (%)  
1 (lowest)  47 (21.3)  60 (27.0)  37 (16.7)  144 (21.7)  0.4321  
2  40 (18.1)  42 (18.9)  44 (19.8)  126 (18.9)  
3  48 (21.7)  40 (18.0)  46 (20.7)  134 (20.2)  
4  46 (20.8)  42 (18.9)  53 (23.9)  141 (21.2)  
5 (highest)  40 (18.1)  38 (17.1)  42 (18.9)  120 (18.0)  
Neighbourhood Immigration Tercile, n (%)  
1 (lowest)  147 (66.5)  133 (59.9)  140 (63.1)  420 (63.2)  0.2832  
2  49 (22.2)  47 (21.2)  48 (21.6)  144 (21.7)  
3 (highest)  25 (11.3)  42 (18.9)  34 (15.3)  101 (15.2)  
Charlson score = 0, n (%)  193 (87.3)  208 (93.7)  215 (96.8)  616 (92.6)  0.0022  
Number of ED Visits in the year before cohort entry, mean (SD)  0.7 (1.4)  0.4 (1.1)  0.3 (0.7)  0.5 (1.1)  0.0051  
Number of OHIP Outpatient Visits in the year before cohort entry, mean (SD)  27.5 (39.6)  17.7 (18.7)  14.5 (16.0)  19.9 (27.4)  <0.0001  
Cancerrelated characteristics  
Screening group, n (%)  
None  118 (53.4)  120 (54.1)  109 (49.1)  347 (52.2)  0.5905  
Repeated  12 (5.4)  16 (7.2)  12 (5.4)  40 (6.0)  
Diagnostic  31 (14.0)  34 (15.3)  30 (13.5)  95 (14.3)  
Sporadic  60 (27.1)  52 (23.4)  71 (32.0)  183 (27.5)  
Stage at diagnosis, n (%)  
IIIA (least advanced)  17 (7.7)  12 (5.4)  32 (14.4)  61 (9.2)  0.0191  
IIIB  134 (60.6)  136 (61.3)  122 (55.0)  392 (58.9)  
IIIC (most advanced)  70 (31.7)  74 (33.3)  68 (30.6)  212 (31.9)  
Tumor Grade Category, n (%)  
Low  173 (78.3)  172 (77.5)  182 (82.0)  527 (79.2)  0.6545  
High  40 (18.1)  40 (18.0)  30 (13.5)  110 (16.5)  
Unknown  8 (3.6)  10 (4.5)  10 (4.5)  28 (4.2)  
Emergent surgery, n (%)  110 (49.8)  65 (29.3)  19 (8.6)  194 (29.2)  <0.0001  
Length of inhospital stay post surgery, mean (SD)  13.6 (17.0)  6.9 (4.7)  5.4 (2.1)  8.6 (10.8)  <0.0001  
ED visits in the first year post diagnosis, n (%)  192 (86.9)  182 (82.0)  106 (47.7)  480 (72.2)  <0.0001  
Number of ED visits post diagnosis, mean (SD)  2.9 (2.1)  2.3 (1.4)  1.7 (1.6)  2.4 (1.8)  <0.0001  
Chemotherapy in the first year post diagnosis, n (%)  60 (27.1)  138 (62.2)  193 (86.9)  391 (58.8)  <0.0001  
Number of chemotherapy visits, mean (SD)  9.5 (4.0)  10.9 (5.6)  11.3 (5.4)  10.9 (5.3)  0.0674  
Statistically significant differences in (unadjusted) survival over time were observed for concordance score terciles (pvalue , Figure 6).
Table 3 summarizes the association between concordance and survival for both unadjusted and adjusted Cox models with both categorical and continuous measures for the concordance score. The unadjusted Cox model found a significant effect of pathway concordance on survival, with decreased risk of mortality associated with having higher values of the concordance score in both the continuous (HR = 0.72, 95% CI = (0.66, 0.79)) and categorical models (HR = 0.51, 95% CI = (0.36, 0.73) for Medium vs. Low tercile; HR = 0.31, 95% CI = (0.20, 0.47) for High vs. Low tercile). Note that a unit change in the continuous model is scaled to represent a change of 0.1 in .
Along with concordance score, the adjusted Cox regression included age, sex, rural residency, number of outpatient visits a year prior to diagnosis, screening group, cancer stage and grade, emergent surgery status, and length of hospital stay following the surgery. A nonlinear relationship between three covariates (age, number of outpatient visits, and length of hospital stay following the surgical treatment) was observed, so restricted cubic splines with three knots were used to model these covariates. No violation of model assumptions was detected. Results of the adjusted analyses showed statistical significance for the association between the concordance score and mortality. Similar to the unadjusted model, there was a significant effect of concordance on survival in the continuous model, with decreased risk of mortality associated with higher values of the concordance score (HR = 0.80, 95% CI = (0.71, 0.90)). However, only the High vs. Low tercile difference was significant (HR = 0.57, 95% CI = (0.36, 0.92)) for the categorical model.
Unadjusted  Adjusted  

HR  95% CI  pvalue  HR  95% CI  pvalue  
Categorical  (Med. vs Low)  0.51  (0.36, 0.73)  0.0002  0.70  (0.48, 1.03)  0.0684 
(High vs Low)  0.31  (0.20, 0.47)  <0.0001  0.57  (0.36, 0.92)  0.0216  
Continuous  (unit change)  0.72  (0.66, 0.79)  <0.0001  0.80  (0.71, 0.90)  0.0003 
Figure 7
summarizes the influence of each covariate in the adjusted continuous Cox model on survival. A hazard ratio below (above) 1 indicates a beneficial (harmful) effect on survival. If the confidence interval contains 1, the effect is not statistically significant at the 0.05 level. Concordance score and female sex have a statistically significant beneficial effect on survival, while cancer stage, cancer grade and emergent surgery have significant harmful effects. The corresponding figure for the adjusted categorical Cox model is provided in Appendix
13.Model validation with 1000 bootstrap resamples showed acceptable coefficient shrinkage (0.84) and model optimism (Cstatistic reduction by 0.024), indicating minor evidence of overfitting. Based on the average of the bootstrap samples, we estimated a bias of 0.003 for the concordance score coefficient (log(HR) = 0.2207 or HR = 0.802 in our model versus log(HR) = 0.2237 or HR = 0.799 for the average of the bootstrap samples). The empirical confidence intervals of the concordance score estimates had similar coverage (0.70, 0.92) as the confidence intervals from the original model. The results of resampling imply model consistency in estimating the effect of the inverse optimizationbased concordance score.
7.3 The Value of Patient Data
Finally, we repeat the above survival analysis using a cost vector (and associated concordance metric) that is derived from the model with only reference pathways as input () in order to quantify the value of our datadriven approach that uses patient pathway data to supplement the reference pathways. We focus on the Cox models (both categorical and continuous) and use the same methodological approach and dataset.
Table 4 summarizes the association between concordance and survival based on the referenceonly cost vector. Similar to the datadriven approach, the unadjusted Cox model found a significant effect of pathway concordance on survival in both the continuous (HR = 0.75, 95% CI = (0.66, 0.86)) and categorical models (HR = 0.66, 95% CI = (0.46, 0.94) for Medium vs. Low tercile; HR = 0.52, 95% CI = (0.36, 0.77) for High vs. Low tercile). However, after adjusting for other relevant covariates, the association weakened significantly. In particular, there is no longer any significant association between concordance score and survival for both the categorical or continuous models. For illustrative purposes, we include the forest plot for only the continuous model in Appendix 13.
Unadjusted  Adjusted  

HR  95% CI  pvalue  HR  95% CI  pvalue  
Categorical  (Med. vs Low)  0.66  (0.46, 0.94)  0.0215  1.07  (0.72, 1.60)  0.7239 
(High vs Low)  0.52  (0.36, 0.77)  0.0008  1.01  (0.65, 1.58)  0.9518  
Continuous  (unit change)  0.75  (0.66, 0.86)  <0.0001  0.94  (0.81, 1.10)  0.4477 
Overall, this analysis shows that a datadriven approach, using patient data to refine the cost vector so alignment with patient outcomes is improved, produces a pathway concordance metric that is much more relevant and statistically meaningful for populationbased health outcome monitoring.
8 Conclusions
This paper proposes the first datadriven inverse optimization approach to measuring pathway concordance in any problem context. Our specific development focuses on clinical pathway concordance for cancer patients. We model the care network as a directed graph and use both recommended pathways and actual patient pathways to develop and refine a novel pathway concordance metric. We confirm the value of our datadriven approach by showing that the concordance scores are strongly associated with survival for patients with stage III colon cancer, even after adjusting for other patient covariates.
In contrast to previous work in clinical pathway concordance, which focus on single events or courses of treatment by patients in a single institution, our concordance score can support populationlevel health system monitoring and reporting. Given its statistically significant association with clinical outcomes, our metric provides meaningful quantitative measures of system efficiency and variation. Improving concordance measurement in a datadriven and generalizable manner allows officials to evaluate complex practice patterns at the population level and describe clinical (e.g., survival) and system (e.g., cost) outcomes in relation to variation in care. These rich data can then be explored for providing decision support to clinicians and health system administrators. Measuring pathway concordance may help identify gaps and bottlenecks in the health system, facilitate modeling of system changes and resulting costs, and allow for interjurisdictional comparison while accommodating system differences. These activities, in turn, can help support planning and priority setting. Concordance metrics optimized for clinical outcomes are also useful in promoting and implementing best practice, helping providers focus on appropriate interventions and reduce inappropriate resource use. These topics are the focus of future work, where concordance measurement will be explored to identify opportunities for improvement along the entire care continuum and among different patient populations.
Our methodological contribution includes a twostage inverse optimization approach to handling primary and secondary data sources, as well as a model that tries to separate good and bad (secondary) data points as much as possible. These methodological advances, motivated by our realworld problem context, suggest avenues for future research in inverse optimization. To keep our models tractable, temporal aspects of concordance measurement were accounted for in how we processed the event logs to form patient pathways. Future work may also consider modifying the clinical activity network so the state space captured by the nodes also include temporal characteristics, such as receipt of treatment before or after a certain time frame.
9 Proofs
Proof.
Proof of Lemma 3.2.1. Since the objective function has a lower bound of zero, it suffices to show that is feasible if and only if there are at least two distinct paths from the start node (denoted ) to the end node (denoted ).
() Let and be two distinct paths from to . Let and . If we set for , for , and for all the other arcs in the network, then satisfies and . Note that is the shortest path for the designed cost vector which has finite total cost. From strong duality, its dual is also optimal meaning that there exists a that satisfies constraint for . Finally, we build a feasible solution for by setting for .
() Suppose to the contrary that there is only one path from to , which means the entire network itself is the path from to . Since there is a feasible , there is some arc with cost . The constraint implies that the incoming arc to and outgoing arc from also has cost equal to 1 in absolute value. Continuing this logic, the same will be true of the outgoing arc from , but this violates , meaning cannot be feasible.
∎
Proof.
Proof of Proposition 3.2.1. The dual of is feasible since and satisfy the dual feasibility conditions. Hence, is not unbounded, implying no negative cost cycles. ∎
Proof.
Proof of Proposition 3.2.2. Let be an optimal solution to . We can construct a feasible solution to by setting for and for . Hence, what remains is to show that the objective is bounded. Substituting the third and fourth constraints into the objective, we have
which is clearly bounded below since and are data and . In particular, if , then by Holder’s inequality the absolute value of the objective is bounded by . ∎
Proof.
Proof of Theorem 4.1.

Given and , the denominator term in is fixed. An optimal solution to , by definition of its objective function, minimizes the numerator of , thus maximizing .

It is clear that since is the ratio of sums of squares, which are nonnegative. To show , note that for all because is a feasible solution to .

An optimal solution to is feasible for , since the latter problem is a relaxation of the former. Invoking the first statement in this theorem, .
∎
Proof.
Proof of Theorem 4.2. Note that , where the first inequality is by definition of and second inequality is by optimality of . It is straightforward to then show that , with equality at the lower or upper bound when the first or second inequalities above on , respectively, are tight. ∎
10 LowerDimensional Transformation for Calculation
Because and has rank , we can partition to two matrices and in which is a square matrix comprising the linear independent columns of while comprises the remaining columns of (Bertsekas 1999). We denote the corresponding decomposition of and by and , respectively. Substituting in the shortest path problem transforms the feasible region to a lower dimensional space where is the variable, and the feasible region is represented by inequality constraints only: where and . The transformed problem is a fulldimensional problem in dimension .
11 Calculation of Longest Walk Length
We use the BellmanFord equations to calculate the longest walk for walks up to a given length . Let the node set be , the start node be , and the end node be . Let be the cost of the longest walk from to that takes steps, computed via backward induction over the equations
where if , and otherwise. We solve this system of equations once up to a value of that equals the length of the longest patient pathway (). Then, the costs of the longest walks with length for each integer up to are stored and can simply be looked up when computing for each patient pathway.
12 Data Sources

Cancer Care Ontario

Ontario Cancer Registry (OCR)

Collaborative Stage Dataset (CSI)

Colonoscopy Interim Reporting Tool (CIRT)

Lab Reporting Tool (LRT)

Screening HUB

Activity Level Reporting (ALR)


Ministry of Health and LongTerm Care (MoHLTC)

Registered Persons Database (RPDB)

Ontario Health Insurance Plan (OHIP) Claims


Canadian Institute for Health Information (CIHI)

National Ambulatory Care Reporting System (NACRS)

Discharge Abstract Database (DAD)


Statistics Canada

Postal Code Conversion file (PCCF+6B)

13 Supplementary Figures for Survival Analysis
References
 Adriansyah et al. (2011) Adriansyah A, van Dongen BF, van der Aalst WM (2011) Conformance checking using costbased fitness analysis. Enterprise Distributed Object Computing Conference (EDOC), 2011 15th IEEE International, 55–64 (IEEE).
 Ahuja and Orlin (2001) Ahuja RK, Orlin JB (2001) Inverse optimization. Operations Research 49(5):771–783.
 Ahuja and Orlin (2002) Ahuja RK, Orlin JB (2002) Combinatorial algorithms for inverse network flow problems. Networks: An International Journal 40(4):181–187.
 Aswani et al. (2018) Aswani A, Shen ZJ, Siddiq A (2018) Inverse optimization with noisy data. Operations Research 66(3):870–892.
 Babier et al. (2019) Babier A, Chan TC, Lee T, Mahmood R, Terekhov D (2019) A unified framework for model fitting and evaluation in inverse linear optimization. arXiv preprint arXiv:1804.04576v2 .
 Bertsekas (1999) Bertsekas D (1999) Nonlinear Programming (Athena Scientific), 2 edition, ISBN 1886529000.
 Bertsimas et al. (2015) Bertsimas D, Gupta V, Paschalidis IC (2015) Datadriven estimation in equilibrium using inverse optimization. Mathematical Programming 153(2):595–633.
 Bose and van der Aalst (2010) Bose RJC, van der Aalst W (2010) Trace alignment in process mining: opportunities for process diagnostics. International Conference on Business Process Management, 227–242 (Springer).
 Burton et al. (1997) Burton D, Pulleyblank W, Toint PL (1997) The inverse shortest paths problem with upper bounds on shortest paths costs. Network Optimization, 156–171 (Springer).
 Burton and Toint (1992) Burton D, Toint PL (1992) On an instance of the inverse shortest paths problem. Mathematical Programming 53(13):45–61.
 Campbell et al. (1998) Campbell H, Hotchkiss R, Bradshaw N, Porteous M (1998) Integrated care pathways. Bmj 316(7125):133–137.
 Cancer Care Ontario (2019a) Cancer Care Ontario (2019a) Colorectal cancer pathway map. URL https://www.cancercareontario.ca/en/pathwaymaps/colorectalcancer.
 Cancer Care Ontario (2019b) Cancer Care Ontario (2019b) Pathway maps. URL https://www.cancercareontario.ca/en/pathwaymaps.
 Chan et al. (2014) Chan TC, Craig T, Lee T, Sharpe MB (2014) Generalized inverse multiobjective optimization with application to cancer therapy. Operations Research 62(3):680–695.
 Chan et al. (2019) Chan TC, Lee T, Terekhov D (2019) Inverse optimization: Closedform solutions, geometry, and goodness of fit. Management Science 65(3):1115–1135.

Charlson et al. (1987)
Charlson ME, Pompei P, Ales KL, MacKenzie CR (1987) A new method of classifying prognostic comorbidity in longitudinal studies: development and validation.
Journal of chronic diseases 40(5):373–383.  Colorectal Cancer Alliance (2019) Colorectal Cancer Alliance (2019) Determine your risk and practice prevention. URL https://www.ccalliance.org/colorectalcancerinformation/statisticsriskfactors.
 Cox (1971) Cox D (1971) Regression models and lifetables (with discussion). Journal of the Royal Statistical Society 34(2):187–220.
 Esfahani et al. (2018) Esfahani PM, ShafieezadehAbadeh S, Hanasusanto GA, Kuhn D (2018) Datadriven inverse optimization with imperfect information. Mathematical Programming 167(1):191–234.
 Evans et al. (2013) Evans WK, Ung YC, Assouad N, Chyjek A, Sawka C (2013) Improving the quality of lung cancer care in ontario: the lung cancer disease pathway initiative. Journal of Thoracic Oncology 8(7):876–882.
 Faragó et al. (2003) Faragó A, Szentesi Á, Szviatovszki B (2003) Inverse optimization in highspeed networks. Discrete Applied Mathematics 129(1):83–98.
 Forestier et al. (2012) Forestier G, Lalys F, Riffaud L, Trelhu B, Jannin P (2012) Classification of surgical processes using dynamic time warping. Journal of biomedical informatics 45(2):255–264.
 Gao et al. (2010) Gao X, Xiao B, Tao D, Li X (2010) A survey of graph edit distance. Pattern Analysis and applications 13(1):113–129.
 Harrell Jr (2015) Harrell Jr FE (2015) Regression modeling strategies: with applications to linear models, logistic and ordinal regression, and survival analysis (Springer).
 Ishwaran et al. (2008) Ishwaran H, Kogalur UB, Blackstone EH, Lauer MS, et al. (2008) Random survival forests. The annals of applied statistics 2(3):841–860.
 James et al. (2018) James PD, Rabeneck L, Yun L, Paszat L, Baxter NN, Govindarajan A, Antonova L, Tinmouth JM (2018) Repeated faecal occult blood testing is associated with decreased advanced colorectal cancer risk: A populationbased study. Journal of medical screening 25(3):141–148.
 Keshavarz et al. (2011) Keshavarz A, Wang Y, Boyd S (2011) Imputing a convex objective function. Intelligent Control (ISIC), 2011 IEEE International Symposium on, 613–619 (IEEE).
 Leshno and Levy (2004) Leshno M, Levy H (2004) Stochastic dominance and medical decision making. Health Care Management Science 7(3):207–215.
 Ling et al. (2018) Ling DC, Karukonda P, Smith RP, Heron DE, Beriwal S (2018) Declining brachytherapy utilization for highrisk prostate cancer  can clinical pathways reverse the trend? Brachytherapy 17(6):895–898.
 Navarro (2001) Navarro G (2001) A guided tour to approximate string matching. ACM computing surveys (CSUR) 33(1):31–88.

Neuhaus and Bunke (2004)
Neuhaus M, Bunke H (2004) A probabilistic approach to learning costs for graph
edit distance.
Proceedings of the 17th International Conference on Pattern Recognition, 2004. ICPR 2004.
, volume 3, 389–393 (IEEE).  Panella et al. (2003) Panella M, Marchisio S, Di Stanislao F (2003) Reducing clinical variations with clinical pathways: do pathways work? International Journal for Quality in Health Care 15(6):509–521.
 Quan et al. (2005) Quan H, Sundararajan V, Halfon P, Fong A, Burnand B, Luthi JC, Saunders LD, Beck CA, Feasby TE, Ghali WA (2005) Coding algorithms for defining comorbidities in icd9cm and icd10 administrative data. Medical care 43(11):1130–1139.

Riesen (2015)
Riesen K (2015) Structural pattern recognition with graph edit distance.
Advances in computer vision and pattern recognition, Cham
.  Ristad and Yianilos (1998) Ristad ES, Yianilos PN (1998) Learning stringedit distance. IEEE Transactions on Pattern Analysis and Machine Intelligence 20(5):522–532.
 Rotter et al. (2012) Rotter T, Kinsman L, James E, Machotta A, Willis J, Snow P, Kugler J (2012) The effects of clinical pathways on professional practice, patient outcomes, length of stay, and hospital costs: Cochrane systematic review and metaanalysis. Evaluation & the health professions 35(1):3–27.
 Samokhvalov et al. (2018) Samokhvalov AV, Probst C, Awan S, George TP, Le Foll B, Voore P, Rehm J (2018) Outcomes of an integrated care pathway for concurrent major depressive and alcohol use disorders: a multisite and prospective study. BMC Psychiatry 18(189).
 Schmidt et al. (2018) Schmidt I, Thor J, Davidson T, Nilsson F, Carlsson C (2018) The national program on standardized cancer care pathways in sweden: Observations and findings halfway through. Health Policy 122(9):945–948.
 Sendi et al. (2003) Sendi P, Al MJ, Gafni A, Birch S (2003) Optimizing a portfolio of health care programs in the presence of uncertainty and constrained resources. Social Science & Medicine 57(11):2207–2215.
 Stone and Koo (1985) Stone CJ, Koo CY (1985) Additive splines in statistics. Proceedings of the American Statistical Association. Original pagination is p 45:48.
 Troutt et al. (2006) Troutt MD, Pang WK, Hou SH (2006) Behavioral estimation of mathematical programming objective function coefficients. Management science 52(3):422–434.
 van de Klundert et al. (2010) van de Klundert J, Gorissen P, Zeemering S (2010) Measuring clinical pathway adherence. Journal of biomedical informatics 43(6):861–872.
 Van Zelm et al. (2018) Van Zelm R, Janssen I, Vanhaecht K, de Buck van Overstraeten A, Panella M, Sermeus W, Coeckelberghs E (2018) Development of a model care pathway for adults undergoing colorectal cancer surgery: Evidencebased key interventions and indicators. Journal of evaluation in clinical practice 24(1):232–239.
 World Cancer Research Fund (2019) World Cancer Research Fund (2019) Colorectal cancer statistics: Colorectal cancer is the third most common cancer worldwide. URL https://www.wcrf.org/dietandcancer/cancertrends/colorectalcancerstatistics.
 World Health Organization (2019) World Health Organization (2019) Colorectal cancer. URL https://gco.iarc.fr/today/data/factsheets/cancers/10_8_9Colorectumfactsheet.pdf?fbclid=IwAR3xpp3qDLrfrCAKXjR1XAM6HUzZYpUPqk8kvaY6llk5eyAs1CrZ4k0Me4.
 Xu and Zhang (1995) Xu S, Zhang J (1995) An inverse problem of the weighted shortest path problem. Japan journal of industrial and applied mathematics 12(1):47.
 Yan et al. (2018) Yan H, Van Gorp P, Kaymak U, Lu X, Ji L, Chiau CC, Korsten HH, Duan H (2018) Aligning event logs to tasktime matrix clinical pathways in bpmn for variance analysis. IEEE Journal of Biomedical and Health Informatics 22(2):311–317.
 Yang et al. (1997) Yang C, Zhang J, Ma Z (1997) Inverse maximum flow and minimum cut problems. Optimization 40(2):147–170.
 Yang et al. (2017) Yang S, Dong X, Sun L, Zhou Y, Farneth RA, Xiong H, Burd RS, Marsic I (2017) A datadriven process recommender framework. Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2111–2120 (ACM).
 Yang et al. (2016) Yang S, Zhou M, Webman R, Yang J, Sarcevic A, Marsic I, Burd RS (2016) Durationaware alignment of process traces. Industrial Conference on Data Mining, 379–393 (Springer).
 Yujian and Bo (2007) Yujian L, Bo L (2007) A normalized levenshtein distance metric. IEEE transactions on pattern analysis and machine intelligence 29(6):1091–1095.
 Zhang and Cai (1998) Zhang J, Cai MC (1998) Inverse problem of minimum cuts. Mathematical Methods of Operations Research 47(1):51–58.

Zhang and Ma (1996)
Zhang J, Ma Z (1996) A network flow method for solving some inverse combinatorial optimization problems.
Optimization 37(1):59–72.  Zhang et al. (1995) Zhang J, Ma Z, Yang C (1995) A column generation method for inverse shortest path problems. Zeitschrift für Operations Research 41(3):347–358.
 Zhao et al. (2015) Zhao Q, Stettner A, Reznik E, Segrè D, Paschalidis IC (2015) Learning cellular objectives from fluxes by inverse optimization. Decision and Control (CDC), 2015 IEEE 54th Annual Conference on, 1271–1276 (IEEE).
Comments
There are no comments yet.