An Inverse Optimization Approach to Measuring Clinical Pathway Concordance

TL;DR

This paper introduces a novel inverse optimization method to measure how closely patient pathways follow clinical guidelines, demonstrating its effectiveness with real data and linking pathway concordance to patient survival.

Contribution

It develops the first data-driven inverse optimization framework for pathway concordance, integrating multiple data sources and addressing data primacy and alignment.

Findings

Concordance metric significantly correlates with patient survival.

Method effectively distinguishes between aligned and deviated pathways.

Real patient data validates the approach's practical utility.

Abstract

Clinical pathways outline standardized processes in the delivery of care for a specific disease. Patient journeys through the healthcare system, though, can deviate substantially from these pathways. Given the positive benefits of clinical pathways, it is important to measure the concordance of patient pathways so that variations in health system performance or bottlenecks in the delivery of care can be detected, monitored, and acted upon. This paper proposes the first data-driven inverse optimization approach to measuring pathway concordance in any problem context. Our specific application considers clinical pathway concordance for stage III colon cancer. We develop a novel concordance metric and demonstrate using real patient data from Ontario, Canada that it has a statistically significant association with survival. Our methodological approach considers a patient's journey as a walk in a directed graph, where the costs on the arcs are derived by solving an inverse shortest path problem. The inverse optimization model uses two sources of information to find the arc costs: reference pathways developed by a provincial cancer agency (primary) and data from real-world patient-related activity from patients with both positive and negative clinical outcomes (secondary). Thus, our inverse optimization framework extends existing models by including data points of both varying "primacy" and "alignment". Data primacy is addressed through a two-stage approach to imputing the cost vector, while data alignment is addressed by a hybrid objective function that aims to minimize and maximize suboptimality error for different subsets of input data.

Tables

Table 1. Table 1: Comparing concordant and discordant arcs in pathways that have edit distance of one from the reference pathway

Pathway	Operation	Missing concordant arcs	Extra discordant arcs	Pathway representation
A-B-C
	-	0	0
A-A-B-C
	Insertion (duplicated activity)	0	1
X-A-B-C
	Insertion (discordant activity)	1	2
B-C
	Deletion	2	1
X-B-C
	Substitution	2	2
B-A-C
	Transposition	3	3

Table 2. Table 2: Risk of mortality associated with concordance score for categorical and continuous models

		Unadjusted			Adjusted
		HR	95% CI	p-value	HR	95% CI	p-value
Categorical	$\omega$ (Med. vs Low)	0.41	(0.35, 0.48)	<0.0001	0.61	(0.52, 0.72)	<0.0001
	$\omega$ (High vs Low)	0.22	(0.18, 0.26)	<0.0001	0.38	(0.31, 0.47)	<0.0001
Continuous	$\omega$ (unit change)	0.70	(0.68, 0.73)	<0.0001	0.80	(0.77, 0.84)	<0.0001

Table 3. Table 3: Risk of mortality associated with concordance score for categorical and continuous models using a cost vector from the model with patient data ( 𝐈𝐎 pat ​ ( 𝒳 ^ r ) ) superscript 𝐈𝐎 pat superscript ^ 𝒳 𝑟 (\mathbf{IO}^{\textrm{pat}}(\hat{\mathcal{X}}^{r})) where the sets 𝒳 ^ s superscript ^ 𝒳 𝑠 \hat{\mathcal{X}}^{s} and 𝒳 ^ d superscript ^ 𝒳 𝑑 \hat{\mathcal{X}}^{d} are swapped

		Unadjusted			Adjusted
		HR	95% CI	p-value	HR	95% CI	p-value
Categorical	$\omega$ (Med. vs Low)	1.21	(1.01, 1.44)	0.0409	1.03	(0.86, 1.24)	0.7139
	$\omega$ (High vs Low)	1.96	(1.67, 2.32)	<0.0001	1.34	(1.13, 1.59)	0.0009
Continuous	$\omega$ (unit change)	1.35	(1.26, 1.44)	<0.0001	1.15	(1.08, 1.24)	<0.0001

Table 4. Table 4: Risk of mortality associated with ω 𝜔 \omega and three baseline edit distance-based concordance scores for continuous model

	Unadjusted			Adjusted
	HR	95% CI	p-value	HR	95% CI	p-value
$\omega$	0.70	(0.68, 0.73)	<0.0001	0.80	(0.77, 0.84)	<0.0001
LCSD	0.73	(0.70, 0.76)	<0.0001	0.83	(0.80, 0.87)	<0.0001
LD	0.75	(0.72, 0.78)	<0.0001	0.85	(0.82, 0.89)	<0.0001
DLD	0.76	(0.73, 0.79)	<0.0001	0.86	(0.83, 0.90)	<0.0001

Table 5. Table 5: Percentage (%) of detour types originating from different nodes in clinical activity network for the entire patient population

	Start	Chest imaging	Surgery
Concordant	44	37	5
Missing only	21	0	5
Extra only (1-5)	7	52	32
Missing and extra (1-5)	28	0	17
Extra only (6-10)	0	8	12
Missing and extra (6-10)	0	0	12
Extra only (>10)	0	3	5
Missing and extra (>10)	0	0	12
Total	100	100	100

Table 6. Table 6: Attribution of discordance ( 1 − ω 1 𝜔 1-\omega ) to missing and extra activities for the entire patient population

	Start	Chest imaging	Surgery
Missing activities	0.024	0.000	0.047
ED visit	0.028	0.007	0.039
Extra consultation	0.002	0.012	0.037
Extra endoscopy	0.000	0.006	0.003
Extra abdomen/pelvis imaging	0.000	0.020	0.071
Extra chest imaging	0.000	0.017	0.091
Total	0.054	0.062	0.288

Table 7. Table 7: Percentage (%) of detour types originating from different nodes in clinical activity network for screened versus unscreened patients

	Start			Chest imaging			Surgery
	Unscr	Scr	Diff	Unscr	Scr	Diff	Unscr	Scr	Diff
Concordant	41	68	-27	39	29	10	6	7	-1
Missing only	19	18	1	0	0	0	4	3	1
Extra only (1-5)	8	5	3	51	67	-16	35	42	-7
Missing and extra (1-5)	32	9	23	0	0	0	12	11	1
Extra only (6-10)	0	0	0	8	4	4	17	17	0
Missing and extra (6-10)	0	0	0	0	0	0	9	7	2
Extra only (>10)	0	0	0	2	0	2	6	5	1
Missing and extra (>10)	0	0	0	0	0	0	11	8	3
Total	100	100	0	100	100	0	100	100	0

Table 8. Table 8: Attribution of discordance ( 1 − ω 1 𝜔 1-\omega ) to missing and extra activities in screened versus unscreened patients

	Start			Chest imaging			Surgery
	Unscr	Scr	Diff	Unscr	Scr	Diff	Unscr	Scr	Diff
Missing activities	0.025	0.019	0.006	0.000	0.000	0.000	0.035	0.025	0.010
ED visit	0.032	0.009	0.023	0.007	0.004	0.003	0.033	0.030	0.003
Extra consultation	0.003	0.000	0.003	0.012	0.014	-0.002	0.036	0.034	0.002
Extra endoscopy	0.000	0.000	0.000	0.007	0.006	0.001	0.003	0.001	0.002
Extra abdomen/pelvis imaging	0.000	0.001	-0.001	0.020	0.018	0.002	0.068	0.071	-0.003
Extra chest imaging	0.000	0.000	0.000	0.017	0.016	0.001	0.082	0.071	0.011
Total	0.060	0.029	0.031	0.063	0.058	0.005	0.257	0.232	0.025

Table 9. Table 9 : Characteristics of stage III colon cancer cohort included in the survival analysis

Characteristic	Low (n=1414)	Medium (n=1414)	High (n=1414)	Total (n=4242)	p-value
Death, n (%)	508 (35.9)	241 (17.0)	131 (9.3)	880 (20.7)	<0.0001
Concordance score, mean (SD)	0.4 (0.1)	0.6 (0.1)	0.8 (0.1)	0.6 (0.2)	<0.0001
Patient characteristics
Age at diagnosis, mean (SD)	72.4 (13.6)	65.9 (12.8)	63.7 (11.4)	67.3 (13.1)	<0.0001
Female, n (%)	747 (52.8)	680 (48.1)	680 (48.1)	2107 (49.7)	0.0145
Rural residency, n (%)	205 (14.5)	201 (14.2)	196 (13.9)	602 (14.2)	0.8886
Neighbourhood Income Quintile, n (%)
1 (lowest)	239 (16.9)	307 (21.7)	289 (20.4)	835 (19.7)	0.0505
2	283 (20.0)	283 (20.0)	305(21.6)	871 (20.5)
3	274 (19.4)	295 (20.9)	253 (17.9)	822 (19.4)
4	299 (21.1)	283 (20.0)	303 (21.4)	885 (20.9)
5 (highest)	239 (16.9)	307 (21.7)	289 (20.4)	835 (19.7)
Neighbourhood Immigration Tercile, n (%)
1 (lowest)	862 (61.0)	893 (63.2)	824(58.3)	2579 (60.8)	0.0077
2	351 (24.8)	301 (21.3)	333 (23.6)	985 (23.2)
3 (highest)	201 (14.2)	220 (15.6)	257 (18.2)	678 (16.0)
Charlson score = 0, n (%)	1046 (74.0)	1187 (84.0)	1246 (88.1)	3479 (82.0)	<0.0001
Number of ED Visits in the year before cohort entry, mean (SD)	0.7 (1.8)	0.4 (1.0)	0.3 (0.8)	0.5 (1.3)	<0.0001
Number of OHIP Outpatient Visits in the year before cohort entry, mean (SD)	5.1 (6.6)	3.8 (4.8)	3.4 (4.7)	4.1 (5.5)	<0.0001
Number of hospital admissions in the year before cohort entry, mean (SD)	0.2 (0.6)	0.1 (0.3)	0.1 (0.3)	0.1 (0.5)	<0.0001
Cancer-related characteristics
Screening group, n (%)
None	870 (61.5)	739 (52.3)	654 (46.3)	2263 (53.3)	<0.0001
Repeated	34 (2.4)	58 (4.1)	69 (4.9)	161 (3.8)
Diagnostic	169 (12.0)	205 (14.5)	216 (15.3)	590 (13.9)
Sporadic	341 (24.1)	412 (29.1)	475 (33.6)	1228 (28.9)
Stage at diagnosis, n (%)
IIIA (least advanced)	113 (8.0)	163 (11.5)	213 (15.1)	489 (11.5)	<0.0001
IIIB	976 (69.0)	985 (69.7)	946 (66.9)	2907 (68.5)
IIIC (most advanced)	325 (23.0)	266 (18.8)	255 (18.0)	846 (19.9)
Tumor Grade Category, n (%)
Low	1130 (79.9)	1183 (83.7)	1184 (83.7)	3497 (82.4)	0.0359
High	248 (17.5)	196 (13.9)	195 (13.8)	639 (15.1)
Unknown	36 (2.5)	35 (2.5)	35 (2.5)	106 (2.5)
Emergent surgery, n (%)	586 (41.4)	290 (20.5)	104 (7.4)	980 (23.1)	<0.0001
Length of in-hospital stay post surgery, mean (SD)	12.0 (23.3)	6.3 (5.1)	4.7 (2.3)	7.7 (14.2)	<0.0001

Table 10. Table 10 : Risk of mortality associated with unnormalized concordance score for categorical and continuous models

		Unadjusted			Adjusted
		HR	95% CI	p-value	HR	95% CI	p-value
Categorical	$ϵ$ (Med. vs Low)	2.07	(1.69, 2.53)	<0.0001	1.60	(1.30, 1.96)	<0.0001
	$ϵ$ (High vs Low)	3.73	(3.09, 4.50)	<0.0001	2.16	(1.75, 2.65)	<0.0001
Continuous	$ϵ$ (unit change)	1.03	(1.02, 1.03)	<0.0001	1.02	(1.01, 1.02)	<0.0001

Equations

$${\mathbf{IO}}^{\text{ref}}({\hat{\mathcal{X}}}^r):\underset{\mathbf{c},\mathbf{p},{\mathbf{ϵ}}^r}{\mathrm{minimize}}$$

$${\mathbf{A}}^{\prime }\mathbf{p}\le \mathbf{c},$$

$${\mathbf{c}}^{\prime }{\hat{\mathbf{x}}}_q^r={\mathbf{b}}^{\prime }\mathbf{p}+{ϵ}_q^r,q=1,\dots ,R,$$

$$\Vert \mathbf{c}{\Vert }_{\infty }=1,$$

$$\mathbf{Ac}=0.$$

$${\mathbf{IO}}_{ij}^{\text{ref}}({\hat{\mathcal{X}}}^r):\underset{\mathbf{c},\mathbf{p},{\mathbf{ϵ}}^r}{\mathrm{minimize}}$$

$${\mathbf{A}}^{\prime }\mathbf{p}\le \mathbf{c},$$

$${\mathbf{c}}^{\prime }{\hat{\mathbf{x}}}_q^r={\mathbf{b}}^{\prime }\mathbf{p}+{ϵ}_q^r,q=1,\dots ,R,$$

$$c_{ij}=-1\vee c_{ij}=1,$$

$$-1\le c_{ij}\le 1,(i,j)\in \mathcal{A},$$

$$\mathbf{Ac}=0.$$

$${\mathbf{IO}}^{\text{pat}}({\hat{\mathcal{X}}}^s,{\hat{\mathcal{X}}}^d,{\mathbf{ϵ}}^{r\ast }):\underset{\mathbf{c},\mathbf{p},{\mathbf{ϵ}}^s,{\mathbf{ϵ}}^d}{\mathrm{minimize}}$$

$${\mathbf{A}}^{\prime }\mathbf{p}\le \mathbf{c},$$

$${\mathbf{c}}^{\prime }{\hat{\mathbf{x}}}_q^r={\mathbf{b}}^{\prime }\mathbf{p}+{ϵ}_q^{r\ast },q=1,\dots ,R,$$

$${\mathbf{c}}^{\prime }{\hat{\mathbf{x}}}_q^s={\mathbf{b}}^{\prime }\mathbf{p}+{ϵ}_q^s,q=1,\dots ,S,$$

$${\mathbf{c}}^{\prime }{\hat{\mathbf{x}}}_q^d={\mathbf{b}}^{\prime }\mathbf{p}+{ϵ}_q^d,q=1,\dots ,D,$$

$$\Vert \mathbf{c}{\Vert }_{\infty }=1,$$

$$\mathbf{Ac}=0.$$

$$\omega (\hat{\mathbf{x}})=1-\frac{{\mathbf{c}}^{\ast }{}^{\prime }\hat{\mathbf{x}}-{\mathbf{c}}^{\ast }{}^{\prime }{\mathbf{x}}^{\ast }}{M(\hat{\mathbf{x}})-{\mathbf{c}}^{\ast }{}^{\prime }{\mathbf{x}}^{\ast }}$$

$$\varphi (({\nu }_1,\dots ,{\nu }_k))=\sum _{i=1}^{k-1}{c}_{{\nu }_i^e{\nu }_{i+1}^s}^{\ast }+\sum _{i=1}^k{c}_{{\nu }_i^s{\nu }_i^e}^{\ast }$$

$$C(\theta )=\frac{\stackrel{\text{extra discordant arcs}}{\overbrace{c_{{\gamma }_i^e{\delta }_1^s}^{\ast }+\varphi ({\Delta }_i)+c_{{\delta }_k^e{\gamma }_{i+l}^s}^{\ast }}}-\stackrel{\text{missing concordant arcs}}{\overbrace{\left(c_{{\gamma }_i^e{\gamma }_{i+1}^s}^{\ast }+\varphi (({\gamma }_{i+1},\dots ,{\gamma }_{i+l-1}))+c_{{\gamma }_{i+l-1}^e{\gamma }_{i+l}^s}^{\ast }\right)}}}{M(\hat{\mathbf{x}})-{\mathbf{c}}^{\ast }{}^{\prime }{\mathbf{x}}^{\ast }}.$$

$$\frac{D}{S}\sum _{q=1}^S{ϵ}_q^s-\sum _{q=1}^D{ϵ}_q^d$$

$$=\frac{D}{S}\sum _{q=1}^S{\mathbf{c}}^{\prime }{\hat{\mathbf{x}}}_q^s-\sum _{q=1}^D{\mathbf{c}}^{\prime }{\hat{\mathbf{x}}}_q^d,$$

$$C(v,k)=\underset{u\in \mathcal{N}}{\max }\{C(u,k-1)+c(u,v)\},k=2,\dots ,N,v\in \mathcal{N},$$

$$C(s,1)=0,$$

$$C(v,1)=\infty ,v\ne s,$$

Full text

\OneAndAHalfSpacedXI\useunder

\ul

\TheoremsNumberedThrough\ECRepeatTheorems\EquationsNumberedThrough\MANUSCRIPTNO

\RUNAUTHOR

Chan et al.

\RUNTITLE

Inverse Optimization for Clinical Pathway Concordance

\TITLE

An Inverse Optimization Approach to Measuring Clinical Pathway Concordance

\ARTICLEAUTHORS

\AUTHOR

Timothy C. Y. Chan \AFFDepartment of Mechanical and Industrial Engineering, University of Toronto, Toronto, Ontario M5S3G8, Canada, \EMAILtcychan@mie.utoronto.ca

\AUTHOR

Maria Eberg \AFFIQVIA, Kirkland, Quebec H9H 5M3, Canada, \EMAILMaria.eberg@mail.McGill.ca

\AUTHOR

Katharina Forster, Claire Holloway, Luciano Ieraci \AFFOntario Health (Cancer Care Ontario), Toronto, Ontario M5G2C1, Canada, \EMAIL{Katharina.Forster, Claire.Holloway, Luciano.Ieraci}@cancercare.on.ca

\AUTHOR

Yusuf Shalaby, Nasrin Yousefi \AFFDepartment of Mechanical and Industrial Engineering, University of Toronto, Toronto, Ontario M5S3G8, Canada, \EMAIL{yusuf.shalaby, nasrin.yousefi}@mail.utoronto.ca

\ABSTRACT

Clinical pathways outline standardized processes in the delivery of care for a specific disease. Patient journeys through the healthcare system, though, can deviate substantially from these pathways. Given the positive benefits of clinical pathways, it is important to measure the concordance of patient pathways so that variations in health system performance or bottlenecks in the delivery of care can be detected, monitored, and acted upon. This paper proposes the first data-driven inverse optimization approach to measuring pathway concordance in any problem context. Our specific application considers clinical pathway concordance for stage III colon cancer. We develop a novel concordance metric and demonstrate using real patient data from Ontario, Canada that it has a statistically significant association with survival. Our methodological approach considers a patient’s journey as a walk in a directed graph, where the costs on the arcs are derived by solving an inverse shortest path problem. The inverse optimization model uses two sources of information to find the arc costs: reference pathways developed by a provincial cancer agency (primary) and data from real-world patient-related activity from patients with both positive and negative clinical outcomes (secondary). Thus, our inverse optimization framework extends existing models by including data points of both varying “primacy” and “alignment”. Data primacy is addressed through a two-stage approach to imputing the cost vector, while data alignment is addressed by a hybrid objective function that aims to minimize and maximize suboptimality error for different subsets of input data.

\KEYWORDS

inverse optimization, network flow, clinical pathway concordance, colon cancer, survival analysis

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 over a specific period of time (Campbell et al. 1998, De Bleser et al. 2006, Kinsman et al. 2010). They prescribe a sequence of steps for managing a clinical process with the aim of optimizing patient or population-level outcomes such as survival, cost, and wait times. Generally, they are defined by medical experts and multidisciplinary teams based on “evidence from local, national and international clinical practice guidelines” (Cancer Care Ontario 2020a). They support multidisciplinary care and evidence-based clinical practice and facilitate variance management (Vanhaecht et al. 2006). 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, Courtney et al. 2020), 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, in-hospital complications, hospital length of stay, and cost of care (Rotter et al. 2010, 2012b, Panella et al. 2003, Schmidt et al. 2018, Vanounou et al. 2007, van Dam et al. 2013). 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 patient-traversed pathways against the recommended pathways, which we will refer to as reference pathways. Such metrics 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 data-driven evidence to inform health policy decisions.

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, Williams et al. 2014), 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. Existing approaches in the literature are either unweighted, which implicitly assume all deviations are equally disadvantageous, or subjectively weighted, which is sensitive to misspecification. Given these limitations, our focus is on developing a weighted concordance metric with objectively derived or data-driven weights.

There are two sources of information that can be used to derive weights for a concordance metric: reference pathways and data from real-world patient activity. Reference pathways should be the primary information source since they are derived from best available medical evidence and the goal is to measure system performance against recommendations. Furthermore, concordance measurement in a health system is ideally optimized for a hierarchy of outcomes, with clinical outcomes usually given priority. Thus, reference pathways should be used directly to derive weights for the concordance metric and for the metric to be strongly associated with outcomes. Real-world patient activity, in the form of individual patient pathways, can serve as a secondary data source to fine tune the weights. However, patient pathways should not be the primary data source since there can be significant deviations from recommendations, and outcomes typically depend on many other factors besides pathway concordance. 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, we use inverse optimization to develop a data-driven concordance metric for measuring clinical pathway concordance that has a statistically significant association with clinical outcomes. Our modeling approach progresses in two phases: the first phase identifies a cost vector that maximizes fitness against the reference pathways; the second phase fine-tunes the cost vector so it further maximizes fitness against the patient pathways with positive outcomes and minimizes fitness against those with negative outcomes.

We emphasize that our goal is to measure performance of the health system, rather than appropriateness of individual patient decisions. Thus, our metric measures the extent to which patients are concordant with reference pathways. Once we have this measure for different subpopulations (geographic, demographic, intervention-based, etc), policy makers can investigate the roots of discordance. For example, issues like congestion should be investigated after system performance is quantified. It should not be accounted for a priori since that would inaccurately characterize system performance by giving a more optimistic picture of health system performance than otherwise exists. Similarly, patient-specific discordance should be investigated post-hoc to see whether the discordance was inappropriate or not. Overall, our metric can help pinpoint discordant activity in a complex health system, which can support policy makers in crafting informed, data-driven decisions about health system interventions, prioritizing them for implementation, and modeling the impact of these interventions on concordance.

We model a patient’s journey through the healthcare system as a walk in a directed graph in which each node represents an activity that the 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 a specific 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 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.

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. 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). Cancer Care Ontario, a business unit of Ontario Health, 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 high-quality, 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. Indeed, one of the stated initiatives of Cancer Care Ontario’s Disease Pathway Management program is to “measure whether care follows the pathway maps” (Cancer Care Ontario 2020b) – this is the problem that our paper addresses.

1.1 Literature Review

1.1.1 Inverse Optimization.

Modern inverse optimization methodologies aim to impute unknown parameters of an optimization problem to make observed decisions minimally suboptimal in either offline settings (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. 2020) or online settings where data arrives sequentially (Dong et al. 2018, Bärmann et al. 2018, Dong and Zeng 2020). In our clinical pathway concordance problem, all data are available a priori in the form of reference pathways and patient pathways. Thus, we use an “offline” inverse optimization model.

Inverse network flow problems are well-studied in the literature (Zhang and Ma 1996, Yang et al. 1997, Zhang and Cai 1998, Ahuja and Orlin 2001, 2002), in particular inverse shortest path problems (Burton and Toint 1992, Xu and Zhang 1995, Zhang et al. 1995, Burton et al. 1997). However, they primarily focus on a single, noiseless observed decision for which it is possible to determine parameters that make the decision exactly optimal. Inverse network flow problems with multiple noisy observations have received only limited attention to date (Faragó et al. 2003, Zhao et al. 2015). Due to the equality constraints corresponding to flow balance, the feasible region in network flow problems is not full dimensional, which means there may exist non-zero cost vectors that are orthogonal to the entire feasible region, making all feasible solutions optimal. The lower dimensionality of the feasible region has not been discussed in prior literature. Part of the reason might be that previous papers have assumed that the distance to some prior cost vector should be minimized (Burton and Toint 1992, Xu and Zhang 1995, Zhang et al. 1995, Burton et al. 1997, Ahuja and Orlin 2001), and this objective may discourage the inverse optimization model from identifying a cost vector that is orthogonal to the entire feasible region. In inverse optimization problems with multiple noisy observations, a more relevant objective is to minimize some measure of error in the model-data fit, and as a result, we need to directly account for the lower dimensionality of the feasible region in the constraints.

Many previous papers focus on the case of non-negative cost vectors, which either facilitates the solution algorithm or is required for a particular application (Faragó et al. 2003, Burton and Toint 1992, Xu and Zhang 1995, Zhang et al. 1995, Burton et al. 1997). In our case, we do not restrict the sign of the cost vector, since activities can be either rewarding or costly in terms of pathway concordance. Additionally, the idea of a two-stage model, with different input datasets and a subset of data to minimize fitness against, has not been proposed in the inverse optimization literature.

1.1.2 Measuring Pathway Concordance.

Clinical pathway concordance measurement is often limited to a binary approach where patient pathways are audited against a set of predefined criteria or sequences of activities and are classified as either concordant or discordant (Cheng et al. 2012, Bergin et al. 2020, Forster et al. 2020). However, the selection of such criteria or sequences can be subjective and certain types of discordance may not be captured. In addition, the selected criteria may vary in the strength of their relationship to the measured outcome, and a binary application of concordance/discordance would not incorporate this important clinical information. Other studies that explore the effectiveness of the clinical pathways compared to traditional care have further simplified concordance measurement by dividing patients into two groups, before and after clinical pathway implementation, and measure the differences of the two groups in terms of the desired outcomes (Vanounou et al. 2007, Delaney et al. 2003, Panella et al. 2003).

More granular methodologies for measuring similarity between pathways have primarily leveraged edit distance as a metric. Edit distance quantifies the similarity between two strings by counting the minimum number of operations required to transform one string to the other. There are many different variations of edit distance, distinguished by the operations they include (Navarro 2001). For example, the Longest Common Subsequence distance considers only insertions and deletions (Apostolico and Guerra 1987). Levenshtein distance is the most common edit distance metric and allows insertion, deletion, and substitution (Levenshtein 1966). Damerau–Levenshtein distance allows insertion, deletion, substitution, and the transposition of two adjacent activities (Damerau 1964). All these edit distance metrics have been used to measure clinical pathway concordance (van de Klundert et al. 2010, Yan et al. 2018, Williams et al. 2014).

Finally, we note that these studies have focused on a single intervention at a single institution (Rotter et al. 2012a), although they could potentially be used at the population level.

1.2 Contributions

Our main contributions in this paper are as follows:

We propose the first inverse optimization-based approach for measuring pathway concordance in any problem context. Methodologically, our approach considers the novel setup where there are primary (reference pathways) and secondary (patient pathways) datasets, which are given differential consideration when imputing the cost vector. Furthermore, our approach is unique in that our data contains subsets that differ in their intended fitness with the cost vector: one that the cost vector aims to fit well and another that it does not. 2. 2.

We develop a metric to assess patient pathway concordance. This metric uses the inversely optimized cost vector to calculate a normalized concordance score for a given flow vector (patient pathway). The concordance metric is simple to compute and easily interpretable, taking values between 0 and 1, with 1 representing perfect concordance. 3. 3.

We provide an in-depth case study applying our inverse optimization-based models and metrics to real patient data from stage III colon cancer patients. In particular, we:

1. (a)

Undertake a rigorous survival analysis to show that there is a statistically significant association between pathway concordance and survival, validating the meaningfulness of our concordance metric; 2. (b)

Validate that our survival analysis results are robust with respect to an unnormalized version of the metric; 3. (c)

Demonstrate that our results are not highly sensitive to potential unobserved confounders; 4. (d)

Demonstrate the value of our approach by showing that a metric developed using only reference pathways can yield an insignificant association between concordance and survival or even a significant association in the incorrect direction; 5. (e)

Demonstrate that our inverse optimization-based concordance metric outperforms several baseline concordance measures based on edit distance. 4. 4.

We develop a rigorous framework to fully apportion the discordance of any patient pathway to its “detours” from the reference pathways. This approach allows us to characterize the contribution of different types of detours, as well as the individual activities within those detours, to total discordance at the individual and population level. Applied to our data, we identify extra imaging and ED visits as primary sources of discordance, consistent with the literature. We also go deeper to provide new insight into how discordance evolves over the course of a patient’s journey: detours tend to be shorter earlier in the pathway and more complex with many additional discordant activities later in the pathway. Finally, we show how discordance in two different subpopulations, repeatedly screened versus unscreened patients, differs from population-level discordance.

Although we focus on colon cancer specifically in this paper, our inverse optimization-based approach is generalizable to other patient cohorts and even beyond the healthcare context to general process management applications. Proofs are included in the Electronic Companion.

2 Motivating a Graph-Based Inverse Optimization Approach

In this paper, we leverage network optimization and the corresponding graph-based representation of the problem to measure pathway concordance. To motivate our approach, we briefly highlight the differences between edit distance-based 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.

First, consider unweighted edit distance algorithms, which measure distance as the minimum number of operations required to transform a source pattern into a target pattern (van de Klundert et al. 2010, Forestier et al. 2012, Yan et al. 2018, Williams et al. 2014). The 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 graph-based 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 operations are applied to nodes and arcs (Gao et al. 2010). 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, pre-defined 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 $S$ and end node $E$ that bookend all pathways and let $A-B-C$ be a reference pathway represented by the sequence of arcs $(S,A)$, $(A,B)$, $(B,C)$ and $(C,E)$. Now consider five example pathways listed in Table 1, all of which have an edit distance of one from $A-B-C$. 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 graph-based distance from the reference pathway. Furthermore, when considering pathways with an edit distance of two or more from the reference pathway, differences with graph-based 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.

The preceding discussion assumes unweighted distance measurement. If we now consider a weighted distance metric, an issue arises where the cost of operations would need to be determined a priori and exogenously from the 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 operations required to transform one into the other. Determining the costs of each operation based on a set of reference and non-reference pathways is complicated by the dependency of the operations 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

We first present a model of the clinical activity network representing 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, accumulating costs along each arc traversed. Then, we present an inverse optimization approach to determining arc costs such that reference pathways are optimal solutions to an appropriate optimization problem defined on this network.

3.1 The Clinical Activity Network

Let $\mathcal{N}$ be a set of $m$ nodes and let $\mathcal{A}\subseteq\mathcal{N}\times\mathcal{N}$ be a set of $n$ 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 $i$, we split each activity node into two nodes, $i_{s}$ and $i_{e}$, representing the start and end, respectively, of activity $i$, and then add an arc $(i_{s},i_{e})$ with cost $c_{i_{s}i_{e}}$. The remaining arcs connect the end node of some activity $i$ with the start node of another activity $j$. Each arc $(i_{e},j_{s})$ has an associated cost $c_{i_{e}j_{s}}$, representing the cost of undertaking activity $j$ immediately after activity $i$. The arc costs are unrestricted in sign. The set $\mathcal{A}$ also includes arcs of the form $(i_{e},i_{s})$, which allows the network to model a patient journey that involves consecutive repetitions of activity $i$. Finally, we add two artificial nodes representing the common start and end of each patient journey. Thus, $m$ 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: $\underset{\mathbf{x}}{\text{min}}\{\mathbf{c}^{\prime}\mathbf{x}:\mathbf{A}\mathbf{x}=\mathbf{b},\mathbf{x}\geq\mathbf{0}\}$. We denote this model the “forward” model and will refer to it as $\mathbf{FO}(\mathbf{c})$. The decision vector $\mathbf{x}\in\mathbb{R}^{n}$ describes the flow on each arc. The matrix $\mathbf{A}\in\mathbb{R}^{(m-1)\times n}$ is the node-arc incidence matrix, but without the $m$th row, which corresponds to the flow balance constraint of the artificial end node. Thus, it is full rank. The parameter $\mathbf{b}\in\mathbb{R}^{m-1}$ is a vector of zeros except with a $+1$ at the artificial start node. The objective function quantifies the minimum cost path from the artificial start to the artificial end 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 traverse an arc several times. Each journey through the network can be represented by a flow vector $\mathbf{x}\in\mathbb{R}^{n}$, where the value of $x_{ij}$ is the number of times arc $(i,j)$ is traversed in the walk. In keeping with the medical terminology, we will refer to a patient journey as a “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 “path-independent”. That is, the cost of traversing arc $(i,j)$ is equal to $c_{ij}$ regardless of how the walk arrives at node $i$. To model path-dependent costs, where the cost of completing activity $j$ immediately following activity $i$ depends on what activities were completed prior to $i$, 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 path-dependent 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 path-independent costs) to maintain a parsimonious model, and instead integrate path-dependent behavior into the constraints of the inverse optimization model.

3.2 Inverse Optimization Model

Our goal is to identify a cost vector $\mathbf{c}\in\mathbb{R}^{n}$ 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.

3.2.1 A Model Using Only Reference Pathways.

Let $\hat{\mathcal{X}}^{r}=\{\hat{\mathbf{x}}^{r}_{1},\ldots,\hat{\mathbf{x}}^{r}_{R}\}$ denote the set of $R$ 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 $\mathbf{p}\in\mathbb{R}^{m-1}$ 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 $\mathbf{c}$ and the $R$ reference pathways:

[TABLE]

The first constraint represents dual feasibility. The second constraint defines the duality gap $\epsilon^{r}_{q}$ for each reference pathway $\hat{\mathbf{x}}^{r}_{q}$ 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., not zero). The fourth constraint is needed to ensure that $\mathbf{c}$ lies in the lower dimensional space defined by the constraints $\mathbf{A}\mathbf{x}=\mathbf{b}$. Since $\mathbf{c}$ lies in the space defined by the intersection of the equality constraints ($\mathbf{a}_{i}^{\prime}\mathbf{x}=b_{i},\;i=1,\ldots,m-1$), it must be orthogonal to the vectors defining those constraints ($\mathbf{a}_{i}$). Without the constraint $\mathbf{A}\mathbf{c}=\mathbf{0}$, an optimal solution (i.e., cost vector) to formulation (1) could be orthogonal to the entire forward feasible region, which would render $\mathbf{c}$ uninformative since every feasible flow vector would be optimal. Note that the constraint $\mathbf{A}\mathbf{c}=\mathbf{0}$ means that the optimal cost vector will be a circulation.

To model additional application-specific considerations, such as user-specified 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, $\mathbf{c}\in\mathcal{C}$, 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 non-trivial network.

Lemma 3.1

$\mathbf{IO}^{\emph{ref}}(\hat{\mathcal{X}}^{r})$* has an optimal solution if and only if there are at least two distinct paths from the start node to the end node.*

Going forward, we assume that the network is non-trivial. Next, we present an important property of an optimal solution to $\mathbf{IO}^{\textrm{ref}}(\hat{\mathcal{X}}^{r})$.

Proposition 3.2

Let $(\mathbf{c}^{*},\mathbf{p}^{*},\boldsymbol{\epsilon}^{r*})$ be an optimal solution to formulation (1). Then, the network with arc costs $\mathbf{c}^{*}$ does not contain any directed negative cost cycles.

This result has an intuitive real-world 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 concordance metric, which we define later.

Formulation (1) is non-convex because of the normalization constraint. However, it can be solved using polyhedral decomposition: solving $2n$ quadratic optimization problems where the normalization constraint is replaced with $c_{ij}\in[-1,1]$ for all $(i,j)\in\mathcal{A}$ and with one additional constraint on one $c_{ij}$ being set to either $+1$ or $-1$:

[TABLE]

Furthermore, if in a particular problem context it is clear which arc should be either the most heavily rewarded (set $c_{ij}=-1$) or heavily penalized (set $c_{ij}=+1$), then only one of these $2n$ 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 $\hat{\mathcal{X}}^{s}=\{\hat{\mathbf{x}}^{s}_{1},\ldots,\hat{\mathbf{x}}^{s}_{S}\}$ denote a dataset of $S$ patient pathways with positive clinical outcomes (i.e., survived) and $\hat{\mathcal{X}}^{d}=\{\hat{\mathbf{x}}^{d}_{1},\ldots,\hat{\mathbf{x}}^{d}_{D}\}$ denote a dataset of $D$ patient pathways with negative 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 this patient data to determine a cost vector that not only maximizes fit with the reference pathways, but also provides “separation” between $\hat{\mathcal{X}}^{s}$ and $\hat{\mathcal{X}}^{d}$. Once formulation (1) is solved and an optimal duality gap vector, $\boldsymbol{\epsilon}^{r*}$, is generated, we use $\boldsymbol{\epsilon}^{r*}$ as input into the following inverse optimization model:

[TABLE]

This model treats $\hat{\mathcal{X}}^{s}$ 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. In contrast, the optimal cost vector from (3) purposefully aims to not fit the data from $\hat{\mathcal{X}}^{d}$. Instead, the cost vector should generate data points as far away as possible from those in $\hat{\mathcal{X}}^{d}$ when solving the forward problem, which is why in formulation (3) the objective maximizes suboptimality with respect to $\hat{\mathcal{X}}^{d}$. The weight $D/S$ simply scales the two sub-objectives 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 $\hat{\mathcal{X}}^{s}$ and $\hat{\mathcal{X}}^{d}$, 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.

Proposition 3.3

$\mathbf{IO}^{\emph{pat}}(\hat{\mathcal{X}}^{s},\hat{\mathcal{X}}^{d},\boldsymbol{\epsilon}^{r*})$* has an optimal solution.*

Looking ahead, our numerical results will use two distinct patient datasets. One dataset comprises both $\hat{\mathcal{X}}^{s}$ and $\hat{\mathcal{X}}^{d}$ and is input to (3) to find an optimal cost vector (see Sections 5.4 and 5.5). Using this cost vector, we calculate concordance scores for patients in the second dataset. These out-of-sample concordant scores form the basis of the subsequent numerical experiments, including the survival analysis where they constitute an independent variable.

4 Concordance Metric

In this section, we develop a metric that measures the concordance of a given patient pathway against the reference pathways. 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 the optimal cost vector.

Using the cost vector $\mathbf{c}^{*}$ generated by inverse optimization, our concordance metric $\omega(\hat{\mathbf{x}})$ measures the cost of a patient pathway $\hat{\mathbf{x}}$ against the cost of a shortest path with respect to $\mathbf{c}^{*}$ as

[TABLE]

where $M(\hat{\mathbf{x}})=\underset{\mathbf{x}}{\text{max}}\{\mathbf{c}^{*}{{}^{\prime}}\mathbf{x}:\mathbf{A}\mathbf{x}=\mathbf{b},\mathbf{x}\geq\mathbf{0},\|\mathbf{x}\|_{1}\leq\|\hat{\mathbf{x}}\|_{1},\mathbf{x}\textrm{ is a walk}\}$ is the cost of a longest walk with up to $\|\hat{\mathbf{x}}\|_{1}$ steps, which can be determined via dynamic programming (see 11). Normalizing in this fashion implies that a walk with more steps 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 $M(\hat{\mathbf{x}})$ for several reasons. First, this normalization allows the computation of $\omega(\hat{\mathbf{x}})$ to depend only on the given pathway $\hat{\mathbf{x}}$ (and $\mathbf{c}^{*}$) and not any other pathway. Normalizing by a fixed amount, $M$, while maintaining the property that $\omega(\hat{\mathbf{x}})\geq 0$ for all $\hat{\mathbf{x}}$ would require that $M$ be longer than all patient pathways that would be considered by $\omega(\hat{\mathbf{x}})$, 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 $M(\hat{\mathbf{x}})$ helps to “un-skew” the distribution of duality gaps over the patient pathways, resulting in an approximately normal distribution of $\omega(\hat{\mathbf{x}})$ for our dataset.

Note that we use the cost of the shortest path, $\mathbf{c}^{*}{{}^{\prime}}\mathbf{x}^{*}$, instead of the cost of a reference pathway when measuring concordance. Since a reference pathway may not actually be a shortest path under $\mathbf{c}^{*}$, and since there may be multiple reference pathways with different costs under $\mathbf{c}^{*}$, we define the concordance metric using $\mathbf{c}^{*}{{}^{\prime}}\mathbf{x}^{*}$ to ensure a consistent baseline. Note that if the inverse optimization model finds a cost vector that perfectly fits all of the reference pathways, then the cost of all reference pathways will be equal to $\mathbf{c}^{*}{{}^{\prime}}\mathbf{x}^{*}$, and $\omega(\hat{\mathbf{x}})$ will be measuring concordance directly against the reference pathways.

It follows from the definition of $\omega(\hat{\mathbf{x}})$ that it is unitless metric between 0 and 1, with 1 representing a perfectly concordant pathway (proof omitted).

Proposition 4.1

Given any feasible solution $\hat{\mathbf{x}}$ to the forward problem:

$\omega(\hat{\mathbf{x}})\in[0,1]$, 2. 2.

$\omega(\hat{\mathbf{x}})=1$* if and only if $\mathbf{c}^{*}{{}^{\prime}}\hat{\mathbf{x}}=\mathbf{c}^{*}{{}^{\prime}}\mathbf{x}^{*}$,* 3. 3.

$\omega(\hat{\mathbf{x}})=0$* if and only if $\hat{\mathbf{x}}\in\underset{\mathbf{x}}{\arg\max}\{\mathbf{c}^{*}{{}^{\prime}}\mathbf{x}:\mathbf{A}\mathbf{x}=\mathbf{b},\mathbf{x}\geq 0,\|\mathbf{x}\|_{1}\leq\|\hat{\mathbf{x}}\|_{1},\mathbf{x}\textrm{ is a walk}\}$.*

Finally, we note that the unnormalized cost difference between the patient pathway ($\hat{\mathbf{x}}$) and a shortest path ($\mathbf{x}^{*}$) using the optimal cost vector ($\mathbf{c}^{*}$), which is exactly the duality gap $\epsilon$ associated with the patient pathway, is another possible measure of pathway concordance. However, we believe the unnormalized metric is less interpretable, since it requires the decision maker to put it in context of the underlying network and arc costs. By definition, the normalized metric does this automatically, which facilitates the comparison of concordance across different patient cohorts or different diseases that may be modeled using different networks. We believe the interpretability of the normalized metric makes it an attractive tool for non-experts. Nevertheless, we will demonstrate robustness of our survival analysis results with respect to the unnormalized metric.

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 evidence-based 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. A pathway map is essentially a flowchart that provides a high level view of the care pathways for a specific disease and shows all the different ways a patient can navigate the system. A single pathway map may encompass many clinical pathways. Pathway maps cover points of care spanning disease prevention, screening, diagnosis, treatment, recovery and end-of-life 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.

From a systems perspective, stage III colon cancer patients are considered a uniform cohort and the optimal management of stage III colon cancer follows a single pathway map (Cancer Care Ontario 2019a). This pathway map spans diagnosis and treatment (see 16), 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 x-ray. The set of concordant pathways for stage III colon cancer are derived from its pathway map.

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 x-rays constitute chest imaging. Chemotherapy is considered concordant if the patient pathway includes at least six 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 person-specific 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 12.

5.2.2 Study Population.

We gathered two datasets for our numerical experiments, which were used as “training” and out-of-sample “testing” sets. The first one is a population-based 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, non-incidental cancer diagnosis, any prior history of cancer or identification of a new primary cancer during the follow-up 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 four-year post-diagnosis observation period. The second dataset is similar to the first one, except that it is for newly diagnosed patients from 2012-2016. After applying the same inclusion criteria, the final 2012-2016 cohort included 4826 patients, 1464 of whom died (from any cause) within their follow-up period.

Each dataset consists of event logs for each patient in the cohort that describe the activities undertaken and their corresponding times relative to the date of diagnosis. However, not all events in the event log are relevant to the diagnostic and treatment phases of care, which is our focus. Therefore we refine each patient’s event log by filtering out events earlier than 30 days prior to diagnosis and following one year post-diagnosis, since patients are expected to complete diagnostic and treatment activities within one year. If the patient completed chemotherapy within the year, we remove all events after their last chemotherapy cycle. We also remove activities that have no bearing on concordance, such as minor procedures like barium enemas. All data refinements were implemented in consultation with our clinician team.

5.3 Network Design

Next, we form 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 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 kept 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 x-rays 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 pathway can end, we include an “outcome layer” with four distinct activities; each pathway must exit the graph via one of these four. Patients who receive six chemotherapy treatments exit through the node CHEMO COMPLETE. Receiving between one and five 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 and died patient pathways from the 2010 dataset constitute the $(\mathbf{A},\mathbf{b}),\hat{\mathcal{X}}^{r}:=\{\hat{\mathbf{x}}_{1}^{r},\hat{\mathbf{x}}_{2}^{r}\}$, $\hat{\mathcal{X}}^{s}$ and $\hat{\mathcal{X}}^{d}$, respectively, in formulations (1) and (3). The value of the $(i,j)$-th component of $\hat{\mathbf{x}}_{q}^{r}$ is 1 if reference pathway $q$ includes arc $(i,j)$ 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 in the pathway. Finally, with respect to $\mathcal{C}$, we included three sets of application-specific constraints: 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. If $i$ represents ABDOMEN CT and $j$ represents PELVIS CT, then the constraint would be $c_{i_{s}i_{e}}\leq c_{j_{s}j_{e}}$. There were eight such constraints in total between the nine non-outcome 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. 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 $i$ represents RESECTION, $j$ represents RESECTION END, $k$ represents MO CONSULT END, and $E$ represents the final end node, then the subpath ranking between RESECTION END and MO CONSULT END would be $c_{i_{e}j_{s}}+c_{j_{s}j_{e}}+c_{j_{e}E}\geq c_{i_{e}k_{s}}+c_{k_{s}k_{e}}+c_{k_{e}E}$. 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 activity through which to exit the graph, we set the cost of the CHEMO COMPLETE to END arc to $-1$, anchoring all the arc costs relative to this one. This constraint is also beneficial from a model tractability perspective because instead of using polyhedral decomposition to solve (1), only a single optimization problem is solved (see Section 3.2.1).

5.5 Inverse Optimization Results

We obtained an optimal cost vector, $\mathbf{c}^{*}$, 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, $\hat{\mathbf{x}}^{r}_{1}$ and $\hat{\mathbf{x}}^{r}_{2}$, and secondarily minimizes (maximizes) the suboptimality error of the patient pathways from survived (died) patients in $\hat{\mathcal{X}}^{s}$ ($\hat{\mathcal{X}}^{d}$). Figure 2 illustrates the $\mathbf{c}^{*}$ values on the graph. Recall that since $\mathbf{c}^{*}$ satisfies $\mathbf{A}\mathbf{c}^{*}=\mathbf{0}$, 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 intra-activity arc $(i_{s},i_{e})$ followed by the arc from the outcome node $i_{e}$ to END).

We calculated the (in-sample) concordance score of each patient pathway from the 2010 dataset using (4). The distribution of concordance scores is shown in Figure 3. The fact that this distribution is approximately normal stems from the choice of normalization in $\omega(\hat{\mathbf{x}})$. As Figure 9 shows, the distributions of the duality gap $\epsilon$ and pathway length $\|\hat{\mathbf{x}}\|_{1}$ 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.

Figure 4 shows smoothed cumulative distributions of the concordance scores, stratified by survival outcome, depicting first-order stochastic dominance. That is, for any given concordance score $\omega_{0}\in[0,1]$, the probability of scoring below $\omega_{0}$ is higher for patients who died. Formal validation of this observation requires a rigorous analysis of the relationship between concordance and survival, which we carry out in the next section.

6 Validating the Concordance Metric against Survival

In this section, we investigate the association between concordance and survival using the 2012-2016 dataset. However, we use the concordance metric developed using the 2010 dataset. As such, this analysis provides out-of-sample validation.

6.1 Methods

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 chi-squared test for categorical variables and a one-way ANOVA for means and standard deviations. A Kaplan-Meier 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 log-rank test statistic.

A semi-parametric 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 $1$ 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, number of in-patient admissions, number 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 sub-stage, which denotes severity. The emergent surgery indicator refers to patients who had an urgent in-hospital 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 right-censored survival data (Ishwaran et al. 2008).

For the survival model, assumption of a linear relationship between continuous independent variables and the outcome was evaluated by plotting the residuals against each covariate. If nonlinearity was detected either from the plots of Martingale residuals or Schoenfeld residuals, restricted cubic splines were used to model the corresponding 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 5000 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.

6.2 Results

A total of 4242 individuals met the inclusion criteria for survival analysis. The mean age at diagnosis was 67.3 years (standard deviation = 13.1). Of the included patients, 880 (20.7%) died. A higher proportion of deaths occurred among those with low concordance scores (35.9%, p-value $<0.0001$).

Table 9 lists the characteristics of the patients in the 2012-2016 dataset and by concordance terciles. Significant differences between the terciles were observed for age, sex, immigrant population, comorbidity score, healthcare utilization, screening group, cancer stage, tumor grade, and length of hospital stay following the surgery. Rural residency and neighbourhood income were not significantly different across terciles at the 0.05 level.

Statistically significant differences in (unadjusted) survival over time were observed for concordance score terciles (p-value $<0.0001$, Figure 5).

Table 2 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 concordance on survival in both the continuous (HR = 0.70, 95% CI = (0.68, 0.73)) and categorical models (HR = 0.41, 95% CI = (0.35, 0.48) for Medium vs. Low tercile; HR = 0.22, 95% CI = (0.18, 0.26) for High vs. Low tercile). Note that a unit change in the continuous model is scaled to represent a change of 0.1 in $\omega(\hat{\mathbf{x}})$.

Along with concordance score, the adjusted Cox regression included age, sex, rural residency, number of outpatient visits a year prior to diagnosis, comorbidity score, screening group, cancer stage and grade, emergent surgery status, and length of hospital stay following the surgery. A nonlinear relationship between two covariates (age 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, with decreased risk of mortality associated with higher values of the concordance score in the continuous model (HR = 0.80, 95% CI = (0.77, 0.84)) and categorical models (HR = 0.61, 95% CI = (0.52, 0.72) for Medium vs. Low tercile; HR = 0.38, 95% CI = (0.31, 0.47) for High vs. Low tercile).

Figure 6 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 has a statistically significant beneficial effect on survival, while age, comorbidity score, cancer stage, emergent surgery, and length of hospital stay following the surgery have significant harmful effects. The corresponding figure for the adjusted categorical Cox model is shown in 13.

Model validation with 5000 bootstrap resamples showed acceptable coefficient shrinkage (0.95) and model optimism (C-statistic reduction by 0.005), indicating no evidence of overfitting. Based on the average of the bootstrap samples, we estimated a bias of $-$0.0011 for the concordance score coefficient (log(HR) = $-$0.2213 or HR = 0.8015 in our model versus log(HR) = $-$0.2224 or HR = 0.8006 for the average of the bootstrap samples). The empirical confidence intervals of the concordance score estimates had similar coverage (0.77, 0.84) as the confidence intervals from the original model. The results of resampling imply model consistency in estimating the effect of the inverse optimization-based concordance score.

6.3 The Value of Patient Data

The purpose of model (3) is to identify an optimal cost vector (with respect to model (1)) that results in a significant association between concordance and survival. Without using the patient data, the association may be insignificant, or significant but with an incorrect directionality of effect. To illustrate this last point, we repeat the survival analysis using a different optimal solution for model (1), one that is generated by considering the extreme case where the two datasets in model (3) are reversed (recall that all solutions to (3) are optimal for (1)).

Table 3 summarizes the association between concordance and survival. Similar to before, the unadjusted Cox models all find a significant association between pathway concordance and survival. However, the effect is going in the opposite direction now, with a higher concordance score associated with an increased mortality risk (e.g., continuous model HR = 1.35, 95% CI = (1.26, 1.44)). After adjusting for other relevant covariates, the association is either insignificant, or remains significant but again in the opposite direction. The effects of the other covariates remain unchanged. For illustrative purposes, we include the forest plot of the continuous model in 13.

6.4 Robustness Checks

In 14 and 15, we include two robustness checks. The first repeats the survival analysis with unnormalized concordance scores. The second tests for sensitivity to potential unobserved confounders. Overall, we find that the results presented above are robust.

7 A Comparison with Edit Distance

In this section, we compare our inverse optimization-based concordance metric, $\omega$, against three baselines for measuring pathway concordance from the literature (van de Klundert et al. 2010, Forestier et al. 2012, Yan et al. 2018, Williams et al. 2014), all based on edit distance.

As baselines, we consider Longest Common Subsequence distance (only insertion and deletion operations), Levenshtein distance (insertions, deletions, and substitutions) and Damerau-Levenshtein distance (insertions, deletions, substitutions, and adjacent transpositions). We refer to these three distance metrics as LCSD, LD, and DLD, respectively. As is most common in the literature, we assume all edits have unit cost.

Next, we repeated the survival analysis from Section 6.1 with the same covariates for each metric. We use the normalized scores, based on the number of activities in both the patient and reference pathways, for all baselines so their values reside in $[0,1]$, similar to $\omega$. Table 4 summarizes the results of survival analysis for both adjusted and unadjusted continuous models. While there is a significant association between concordance and survival for each approach, the hazard ratio is lower for $\omega$ and the difference to all other baselines is statistically significant (p-value $<0.0001$), based on a bootstrap calculation with 1000 samples.

Given this result, we hypothesized that as the sample size is decreased, the baselines would lose their statistically significant association with survival before $\omega$. To study this issue, we report the median HR and p-values based on 1000 samples of the original dataset for sample sizes ranging from 300 to 1700 patients (Figure 7). For each sample size, $\omega$ has the lowest HR and lowest p-value, indicating a stronger association with survival. In addition, at 300 patients, the association for the LD and DLD metrics becomes insignificant at the 0.05 level. If we consider a stricter measure for statistical significance such as 0.01, then sample sizes above 600 are needed for LD and DLD, and 400 for LCSD, whereas 300 patients remain sufficient for $\omega$.

Overall, we believe our concordance metric has modest advantages over existing approaches based on the survival analysis. Next, we demonstrate an additional advantage of $\omega$. Namely, we show that our approach provides a mathematically rigorous link between a pathway’s concordance score and its decomposition into discordant subpathways.

8 Points of Discordance

In this section, we demonstrate how our framework allows for the rigorous decomposition of a patient pathway into its constitute concordant and discordant parts, which in turn allows for a clean attribution of a pathway’s concordance score to interpretable subpathways. When aggregated over a population, we can identify key activities or sequences of activities that are the root of major discordances in the system. We can also compare subpopulations and investigate differences in their points of discordance. Our analyses are conducted using the 2012-2016 dataset.

8.1 Pathway Decomposition

Let $\mathbf{x}^{*}$ be a reference pathway, assumed to be a shortest path in the clinical activity network, defined by $K$ ordered concordant activities $\gamma_{1},\ldots,\gamma_{K}$. We represent $\mathbf{x}^{*}$ as $(\gamma_{0},\gamma_{1},\gamma_{2},\ldots,\gamma_{K},\gamma_{K+1})$, where $\gamma_{0}$ and $\gamma_{K+1}$ are the artificial START and END nodes of the network. Let $\hat{\mathbf{x}}$ be an arbitrary patient pathway written as $(\gamma_{0},\Delta_{0},\Gamma_{1},\Delta_{1},\Gamma_{2},\Delta_{2},\ldots,\Gamma_{K},\Delta_{K},\gamma_{K+1})$, where $\Gamma_{i}$ represents potential concordant activities and $\Delta_{i}$ represents potential discordant activities. Each $\Gamma_{i},i=1,\ldots,K$, is either $\emptyset$ or $\gamma_{i}$. It is $\emptyset$ if concordant activity $\gamma_{i}$ is missing or if a downstream concordant activity $\gamma_{j},j>i$, is completed before activity $\gamma_{i}$ is visited for the first time. Our modeling framework adheres to the principle of optimality, so anytime a concordant activity $\gamma_{j}$ is reached in a pathway, then the reference pathway from $\gamma_{j}$ onwards is the lowest cost subpathway to END, even if concordant activities $\gamma_{i},i<j$, are missing. In other words, it costs more to go back and complete missed activities than it does to simply continue from $\gamma_{j}$ onwards. Each $\Delta_{i}$ is either $\emptyset$ or an ordered set of discordant activities, which may include discordant activities like an emergency department visit or repeated imaging. As a convention, we assume that $\Delta_{i}=\emptyset$ if $\Gamma_{i}=\emptyset$, to avoid having to make a subjective choice of which set, $\Delta_{i}$ or $\Delta_{i-1}$, is assigned potential discordant activities if $\gamma_{i}$ is missing.

A detour is formed any time the patient pathway deviates from the reference pathway. The detour ends when it rejoins the reference pathway at a later concordant activity (or the END node). The simplest example of a detour would be when there is a set of discordant activities $\Delta_{i}\neq\emptyset$ between concordant activities $\gamma_{i}$ and $\gamma_{i+1}.$ The detour comprises the extra arcs between $\gamma_{i}$ and $\gamma_{i+1}$. Other examples of detours include missing a concordant activity $\gamma_{i}$ ($\Gamma_{i}=\Delta_{i}=\Delta_{i-1}=\emptyset$) or combinations of missing and extra steps (e.g., $\Gamma_{i}=\Delta_{i}=\emptyset$ but $\Delta_{i-1}\neq\emptyset$).

Next, we define the cost of a detour. Given a sequence of activities $(\nu_{1},\ldots,\nu_{k})$, we define

[TABLE]

to be the cost along all arcs formed by this sequence using the optimal cost vector, where $\nu_{i}^{s}$ and $\nu_{i}^{e}$ are the start and end node for activity $\nu_{i}$ (recall construction of clinical activity network in Section 3.1). Now, consider a detour $\theta$ that starts at concordant activity $\gamma_{i}$, ends at concordant activity $\gamma_{i+l}$ (skipping activities $\gamma_{i+1},\ldots,\gamma_{i+l-1}$), and includes discordant activities $\Delta_{i}=(\delta_{1},\ldots,\delta_{k})$ in between. We define the cost of this detour, $C(\theta)$, based on the cycle that is induced by the extra discordant activities and missing concordant activities between $\gamma_{i}$ and $\gamma_{i+l}$. That is,

[TABLE]

If detour $\theta$ has only missing concordant activities but no extra discordant activities, then the cost of extra discordant arcs is simplified to $c^{*}_{\gamma_{i}^{e}\gamma_{i+l}^{s}}$. On the other hand, if detour $\theta$ has no missing concordant activities and only extra discordant activities following activity $\gamma_{i}$ ($l=1$), then the cost of missing concordant arcs is simplified to $c^{*}_{\gamma_{i}^{e}\gamma_{i+1}^{s}}$.

With this definition of detour cost, we can fully partition the discordance of any patient pathway $\hat{\mathbf{x}}$, which is captured by its duality gap $\epsilon$, into the cost of its detours.

Theorem 8.1

Let $\Theta$ be the set of all detours in patient pathway $\hat{\mathbf{x}}$. Then $\omega(\hat{\mathbf{x}})=1-\sum_{\theta\in\Theta}C(\theta).$

Theorem 8.1 states that the discordance associated with a patient pathway is entirely attributable to its detours. In a similar fashion, detour costs can be further broken down by the individual missing and/or extra activities in the detour.

In the next two subsections, we illustrate how this concept of detour costs can be used to identify the major sources of discordance in the population and between different subpopulations. We note that some of the subsequent analysis, particularly descriptive statistics around prevalence of discordant activities and detours, can be accomplished without inverse optimization. However, these results are firmly grounded in our network-based approach to modeling the concordance problem. Modeling the concordance problem as an inverse shortest path problem allows us to leverage the graph structure to very precisely define detours. And as shown in Theorem 1, discordance is exactly equivalent to detour costs. All of this leverages the principle of optimality, which makes detours interpretable and gives the concordance metric an additive property over the detours. Interpretability is a critical application of this work. The purpose of identifying points of discordance is to address them subsequently. To the extent that minimizing discordance necessitates change management in the clinical environment, stakeholder buy-in requires that the data prompting change is accurate and interpretable.

8.2 Population-level Discordance Analysis

Figure 8 provides a high-level view of how patients flow through the clinical activity network. Roughly half of the patients start with the first concordant activity, endoscopy. The majority of the discordant individuals go directly to abdomen or chest imaging. About half of the patients complete chemotherapy and overall only 3% of patients follow a shortest path.

Figure 12 breaks down the proportion of the population discordance by the starting node of detours. We observe that the majority of discordance occurs following surgery. This means that detours after surgery are more prevalent and/or more costly (in terms of discordance score) than after other activities. The next two nodes where the greatest discordance originates is after chest imaging and at the START node (i.e., before completing a single concordant activity). We focus on detours originating from these three nodes in the remainder of the analysis, since they comprise about 94% of the total population discordance.

Table 5 shows that detours that originate closer to the start of the pathway tend to be shorter, while those that originate towards the end of the pathway, around the surgery event, tend to be longer and more complex. From the START node, 56% of the patient pathways begin with a detour. These detours are relatively short, with at most five extra discordant activities, and many that only have missing activities with no extra activities. The chest imaging node, which is visited by 84% of the patient population, generates a mix of concordant transitions, short detours (i.e., fewer than five extra discordant activities), and a small fraction of long detours (i.e., six or more extra discordant activities). The most complex detours arise from the surgery node. After surgery, 41% of patients take a detour with over five extra discordant activities. This behavior is nonexistent at the start of the pathway (0%) and infrequent after chest imaging (11%).

Finally, we decompose the cost of a detour into costs associated with individual missing and extra activities based on the proportion of the detour cost in the arcs associated with individual activities. Table 6 clearly shows that the detour costs originating from the surgery node are dominant. Following surgery, discordance is primarily due to extra imaging activities, then missing chemotherapy, ED visits, and extra clinical consultations. Following a concordant chest imaging step, the largest detour costs are due to extra imaging, followed by extra consultations, ED visits, and extra endoscopies. Finally, from the START node, the two major contributors to discordance are ED visits and missing concordant activities (primarily endoscopies).

A key insight from these results is that the distribution of detour types and discordant activities that comprise these detours (both missing and extra activities) varies as patients progress through the pathway map. While extra imaging as a source of discordance has been demonstrated previously (Simos et al. 2015), our results highlight that they primarily contribute to higher discordance scores after surgery has been completed. ED visits, which are generally undesirable from a health system perspective (Doan et al. 2019, Yun et al. 2017, Sun et al. 2013) account for a much higher proportion of total discordance at the diagnostic phase of patient pathways. As patients progress on their cancer journey, they tend to accumulate extra discordant activities, especially after surgery. This result suggests that extra follow-ups may be needed with patients following surgery, to make sure they are following the appropriate treatment course to complete their cancer care journey. Since the intended use of concordance measurement is quality improvement, monitoring variation in patient pathways is critical to identifying and eventually reducing detours.

Absence of endoscopy should be investigated as a potential quality gap. Policy makers can then tailor quality improvement initiatives based on the possible reasons for missing endoscopy. For example, if the absence of endoscopy is due to emergency presentation, then initiatives to increase screening participation should be considered. If the reason is difficulty in accessing endoscopy services, then process or capacity improvement initiatives should be undertaken. If the reason is an education gap (e.g., a CT scan was requested to evaluate symptoms rather than an endoscopy), the initiative should include improved knowledge translation and exchange for primary care providers.

ED visits are another major point of discordance. To address this quality gap, the first step is to determine which ED visits are avoidable. ED visits that did not lead to hospital admission could potentially be managed in another setting like outpatient clinics. Identifying the presenting symptoms leading to an ED visit may yield targeted quality improvement opportunities. For example, a telephone triage service or nurse-led clinics could be tapped to handle non-emergent issues arising between scheduled follow-up visits. Moreover, additional points of discordance (e.g., extra imaging) may occur as a direct result of ED attendance. Identifying these subsequent discordant activities is critical in assessing the total impact of interventions to reduce ED visits. For example, suppose it was determined that a large number of ED visits could be avoided with the establishment of a telephone triage service. Our concordance metric can be used to model the impact of that service on the measure of concordance as well as the net change in cost from reduced ED visits, including avoidance of subsequent additional imaging that is a direct result of those ED visits.

8.3 Comparing Screened and Unscreened Patients

This final subsection studies how two subpopulations of patients – those who get repeated screening versus those who do not get any screening – differ in their discordance. Comparing these two particular subpopulations is meaningful because organized screening programs are associated with a decreased risk of developing advanced stages of cancer (James et al. 2018). Patients with repeated screening, defined as those with “two or more FOBTs [fecal occult blood tests] where the latest FOBT was performed within nine months of diagnosis and a second was performed 12–24 months prior to the latest FOBT” were significantly less likely to present with stage IV colorectal cancer later. In both groups, we only include patients with stage IIIB colon cancer (controlling for substage with the largest of the three substages), no comorbidities, and that are between 60 and 79 years of age. After applying these inclusion criteria, there were 76 patients in the screened group and 593 patients in the unscreened group.

Table 7 summarizes pathways of the screened and unscreened groups following the start, chest imaging, and surgery nodes. From the start node, 68% of patients in the screened group are concordant compared to only 41% in the unscreened group; unscreened patients have significantly more detours with both missing concordant activities and extra discordant activities. Unscreened patients also tend to have longer detours following chest imaging and surgery. Overall, the pathways of patients in the unscreened group are less likely to complete endoscopy and chemotherapy, and are more likely to have extra activities.

Table 8 breaks down discordance scores $(1-\omega)$ for the two subpopulations. For each node, the average value of $1-\omega$ is greater in the unscreened group than the screened group. At the start node, the difference (0.060 vs. 0.029) is largely due to two factors: unscreened patients are more likely to miss concordant activities (19% of the net difference in $1-\omega$) and unscreened patients are incurring more ED visits (74% of the net difference in $1-\omega$). The latter appears to be a key difference between these two groups. ED visits followed by missing endoscopy are the two major contributors to discordance at the start node in the general population (Table 6) as well as the unscreened group. However, in the screened group, the discordance due to both are significantly diminished, especially ED visits, and their order is swapped. From chest imaging, the higher discordance of the unscreened group (0.063 vs. 0.058) is largely due to extra imaging (60% of the net difference in $1-\omega$). Finally, at the surgery node, the higher discordance of the unscreened group (0.257 vs. 0.232) is largely due to extra imaging activities (32% of the net difference in $1-\omega$) and missing chemotherapy (40%) with additional discordance primarily attributed to ED visits and extra consultations.

James et al. (2018) hypothesized that the reason for improved patient outcomes among patients who are part of an organized screening program is that they are more likely to have access to other program features and their associated benefits. These patients have better access to care and are more likely to seek it. Our results support this hypothesis and provide more granular insight into where the potential benefit arises at different parts of the pathway.

9 Conclusions

This paper proposes the first data-driven inverse optimization approach to measuring pathway concordance in any problem context. Our specific development focuses on clinical pathway concordance for cancer patients. 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 work focuses on population-level health system monitoring and reporting. Given its statistically significant association with clinical outcomes, our metric provides meaningful quantitative measures of system efficiency and variation. Using our metric to identify the points of discordance in a health system allows officials to evaluate complex practice patterns at the population level and describe clinical outcomes in relation to variation in care. 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.

As a concrete application, we perform an in-depth case study using stage III colon cancer patient data, demonstrating how our concordance metric identifies population-level discordance, as well as points of discordance between patients who had repeated screening versus those who did not. Our framework rigorously attributes discordance to different sources throughout the pathway, ultimately providing targets for quality improvement initiatives and interventions to improve system performance. The development of our population-level metric is enabled by access to population-level administrative health data and pathway maps that are designed for the population, rather than individuals. Future work will focus on the application of the metric for identifying and prioritizing quality gaps that can be addressed with specific quality improvement initiatives, and on modeling the impact of different interventions on population concordance.

\ECSwitch\ECHead

Electronic Companion

10 Proofs

Proof 10.1

*Proof of Lemma 3.1. Since the objective function has a lower bound of zero, it suffices to show that $\mathbf{IO}^{\textrm{ref}}(\hat{\mathcal{X}}^{r})$ is feasible if and only if there are at least two distinct paths from the start node (denoted $s$) to the end node (denoted $t$).

($\Leftarrow$) Let $P_{1}$ and $P_{2}$ be two distinct paths from $s$ to $t$. Let $P_{1}\setminus P_{2}=\{(i,j)\in\mathcal{A}\,|\,(i,j)\in P_{1}\text{ and }(i,j)\notin P_{2}\}$ and $P_{2}\setminus P_{1}=\{(i,j)\in\mathcal{A}\,|\,(i,j)\in P_{2}\text{ and }(i,j)\notin P_{1}\}$. If we set $c_{ij}=-1$ for $(i,j)\in P_{1}\setminus P_{2}$, $c_{ij}=1$ for $(i,j)\in P_{2}\setminus P_{1}$, and $c_{ij}=0$ for all the other arcs in the network, then $\mathbf{c}$ satisfies $\|\mathbf{c}\|_{\infty}=1$ and $\mathbf{A}\mathbf{c}=\mathbf{0}$. Note that $P_{1}$ is the shortest path for the designed cost vector $\mathbf{c}$ which has finite total cost. From strong duality, its dual is also optimal meaning that there exists a $\mathbf{p}$ that satisfies constraint $\mathbf{A}^{\prime}\mathbf{p}\leq\mathbf{c}$ for $\mathbf{c}$. Finally, we build a feasible solution for $\mathbf{IO}^{\textrm{ref}}(\hat{\mathcal{X}}^{r})$ by setting $\epsilon^{r}_{q}=\mathbf{c}^{\prime}\hat{\mathbf{x}}^{r}_{q}-\mathbf{b}^{\prime}\mathbf{p}$ for $q=1,\ldots,R$.

($\Rightarrow$) Suppose to the contrary that there is only one path from $s$ to $t$, which means the entire network itself is the path from $s$ to $t$. Since there is a feasible $\mathbf{c}$, there is some arc $(i,j)$ with cost $|\mathbf{c}_{ij}|=1$. The constraint $\mathbf{A}\mathbf{c}=\mathbf{0}$ implies that the incoming arc to $i$ and outgoing arc from $j$ also has cost equal to 1 in absolute value. Continuing this logic, the same will be true of the outgoing arc from $s$, but this violates $\mathbf{A}\mathbf{c}=\mathbf{0}$, meaning $\mathbf{c}$ cannot be feasible. \Halmos*

Proof 10.2

Proof of Proposition 3.2. The dual of $\mathbf{FO}(\mathbf{c}^{*})$ is feasible since $\mathbf{c}^{*}$ and $\mathbf{p}^{*}$ satisfy the dual feasibility conditions. Hence, $\mathbf{FO}(\mathbf{c}^{*})$ is not unbounded, implying no negative cost cycles. \Halmos

Proof 10.3

Proof of Proposition 3.3. Let $(\mathbf{c}^{*},\mathbf{p}^{*},\boldsymbol{\epsilon}^{r*})$ be an optimal solution to $\mathbf{IO}^{\textrm{ref}}(\hat{\mathcal{X}}^{r})$. We can construct a feasible solution to $\mathbf{IO}^{\textrm{pat}}(\hat{\mathcal{X}}^{s},\hat{\mathcal{X}}^{d},\boldsymbol{\epsilon}^{r*})$ by setting $\epsilon^{s}_{q}=\mathbf{c}^{*}{{}^{\prime}}\hat{\mathbf{x}}^{s}_{q}-\mathbf{b}^{\prime}\mathbf{p}^{*}$ for $q=1,\ldots,S$ and $\epsilon^{d}_{q}=\mathbf{c}^{*}{{}^{\prime}}\hat{\mathbf{x}}^{d}_{q}-\mathbf{b}^{\prime}\mathbf{p}^{*}$ for $q=1,\ldots,D$. Hence, what remains is to show that the objective is bounded. Substituting the third and fourth constraints into the objective, we have

[TABLE]

which is clearly bounded below since $\hat{\mathbf{x}}^{s}_{q}$ and $\hat{\mathbf{x}}^{d}_{q}$ are data and $\|\mathbf{c}\|_{\infty}=1$. In particular, if $M=\max_{q,q^{\prime}}\{\|\hat{\mathbf{x}}^{s}_{q}\|,\|\hat{\mathbf{x}}^{d}_{q^{\prime}}\|\}$, then by Holder’s inequality the absolute value of the objective is bounded by $2DM$. \Halmos

Proof 10.4

Proof of Theorem 8.1.

Taking the patient pathway $\hat{\mathbf{x}}$ and subtracting the reference pathway $\mathbf{x}^{*}$ results in a set of disjoint cycles, where the forward arcs in each cycle correspond to extra discordant activities (in the patient pathway but not in the reference pathway) and the backward arcs correspond to missing concordant activities (in the reference pathway but not in the patient pathway). If there are no repeated discordant activities within a detour, then there is a 1-1 correspondence between detours and cycles. Otherwise, a detour will generate one cycle with forward and backward arcs, as well as one or more cycles with only forward arcs corresponding to the repeated discordant activities. By construction (see equation (5)), the numerator in the cost of a detour $\theta$ is equal to the cost around the cycle(s) that it generates. Thus, the cost difference $\mathbf{c}^{*}{{}^{\prime}}\hat{\mathbf{x}}-\mathbf{c}^{*}{{}^{\prime}}\mathbf{x}^{*}$ equals the sum of the costs around all cycles. Since each cycle belongs to only one detour $\theta\in\Theta$, $\omega(\hat{\mathbf{x}}):=1-(\mathbf{c}^{*}{{}^{\prime}}\hat{\mathbf{x}}-\mathbf{c}^{*}{{}^{\prime}}\mathbf{x}^{*})/(M(\hat{\mathbf{x}})-\mathbf{c}^{*}{{}^{\prime}}\mathbf{x}^{*})=1-\sum_{\theta\in\Theta}C(\theta)$. \Halmos

11 Calculation of the Cost of Longest Walk ($M(\hat{\mathbf{x}})$)

We use the Bellman-Ford equations to calculate the longest walk for walks with up to a given number of steps. Let the node set be $\mathcal{N}$, the start node be $s$, and the end node be $t$. Let $C(v,k)$ be the cost of the longest walk from $s$ to $v\in\mathcal{N}$ that takes $k$ steps, computed via backward induction over the equations

[TABLE]

where $c(u,v)=c_{uv}$ if $(u,v)\in\mathcal{A}$, and $c(u,v)=\infty$ otherwise. We solve this system of equations once up to a value of $N$ that equals to the cost of the longest patient pathway ($\max_{q,q^{\prime}}\{\|\hat{\mathbf{x}}^{s}_{q}\|,\|\hat{\mathbf{x}}^{d}_{q^{\prime}}\|\}$). Then, the costs of the longest walks with the number of steps for each integer up to $N$ are stored and can simply be looked up when computing $\omega(\hat{\mathbf{x}})$ 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 Long-Term 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 Tables and Figures

Table 9 lists the characteristics of the patients in the 2012-2016 dataset and by concordance terciles.

Figure 9 shows the distributions of the duality gap and pathway length for 2010 dataset. Figure 10 summarizes the influence of each covariate in the adjusted categorical Cox model from Section 6.2 on survival. Figure 11 depicts a similar forest plot, but for the adjusted continuous Cox model with the concordance score based on the reference pathway-only model from Section 6.3. Figure 12 shows the distribution of discordance due to detours starting from different nodes in the pathway.

14 Unnormalized Concordance Scores

To test the robustness of our survival analysis results, we repeated the analysis using the duality gap $\epsilon(\hat{\mathbf{x}}):=\mathbf{c}^{*}{{}^{\prime}}\hat{\mathbf{x}}-\mathbf{c}^{*}{{}^{\prime}}\mathbf{x}^{*}$ as an unnormalized concordance metric, instead of $\omega(\hat{\mathbf{x}})$. As shown in Table 10, we observe a significant association between concordance and survival. Note that smaller values of $\epsilon$ imply higher concordance, which is why an HR value greater than 1 indicates the appropriate association between concordance and survival. Furthermore, the same covariates that were significant in the original analysis remain significant now.

15 Sensitivity Analysis

As a final sensitivity analysis, we investigate how the association between concordance and survival can be affected by a potential unobserved confounder. For example, patient behavioral characteristics are a well-known confounder (Mackian et al. 2004, Norman and Connor 2015). Patients with appropriate health-seeking behavior may be more likely to follow guidelines (Sewitch et al. 2007) and thus have higher concordance. But they are also more likely to engage in healthy behaviours that may be associated with higher survival due to reductions in cancer and non-cancer-related risk factors such as physical activity, smoking and diet.

Our approach follows an established Monte Carlo simulation method (Lash and Fink 2003, Albert et al. 2012, Corrao et al. 2014). We represent the unobserved confounder as a binary variable, $U$, where a value of 1 indicates a trait that affects both concordance and survival (e.g., health-seeking behavior). Each patient is placed in one of six “concordance-mortality” categories: a combination of the concordance ($\omega$) terciles from Section 6.2 (low, medium, and high) and the clinical outcome (survived or died). The prevalence of $U=1$ differs by category.

We assume the prevalence of $U=1$ in patients who survived is $b$-fold higher than those who died. We also assume the prevalence of $U=1$ in high concordant (resp. medium concordant) patients is $RR$-fold higher than medium concordant (resp. low concordant) patients, where $RR$ denotes the Risk Ratio. A Risk Ratio of 1 is equivalent to no relationship between $U$ and concordance, representing a scenario when the estimated coefficient for the concordance score is unchanged.

We consider multiple scenarios, distinguished by the base prevalence assigned to the medium-survived category. The prevalence in the other categories is determined by $b$ and $RR$. In each scenario, we use the prevalence of $U=1$ in each category as the parameter of a Bernoulli random variable to generate a value of $U$ for each patient in that category. After generating $U$ for all the patients, we repeat the survival analysis from Section 6.2 using the same models as before but with the additional covariate $U$. We report the HR for the concordance metric and the corresponding confidence intervals based on 100 simulations for each scenario.

Figure 13 presents the results of the sensitivity analysis for four scenarios: base prevalence values of 0.1, 0.2, 0.3, and 0.4. We observe that for any base prevalence and $b$, HR increases towards 1 as $RR$ increases. Overall, only for high values of all key parameters – $b$, $RR$, base prevalance – does the significant association between concordance and survival vanish. For example, assuming a base prevalence of 10%, we require a 5-10 times higher prevalence in the survived group over the died group and a 10 times higher prevalence in adjacent concordance terciles for the association to become insignificant. Given that the distribution must be very highly skewed towards the high concordance-low mortality risk individuals for our previous results to be affected, the sensitivity analysis suggests that our association is not highly sensitive to unobserved confounders.

16 Colon Cancer Pathway Maps

The diagnosis and treatment pathway maps for colon cancer is shown below. Note that colon and rectum cancers (colorectal) have the same diagnosis pathway map, but the treatment pathway maps are different.

*See pages - of pathway1.pdf *See pages - of pathway2.pdf