Finding Regions of Heterogeneity in Decision-Making via Expected Conditional Covariance
Year: 2021

Abstract:
Individuals often make different decisions when faced with the same context, due to personal preferences and background. For instance, judges may vary in their leniency towards certain drug-related offenses, and doctors may vary in their preference for how to start treatment for certain types of patients. With these examples in mind, we present an algorithm for identifying types of contexts (e.g., types of cases or patients) with high inter-decision-maker disagreement. We formalize this as a causal inference problem, seeking a region where the assignment of decision-maker has a large causal effect on the decision. Our algorithm finds such a region by maximizing an empirical objective, and we give a generalization bound for its performance. In a semi-synthetic experiment, we show that our algorithm recovers the correct region of heterogeneity accurately compared to baselines. Finally, we apply our algorithm to real-world healthcare datasets, recovering variation that aligns with existing clinical knowledge.

heading: 1 Introduction
text: Understanding heterogeneity in decision-making is an established problem in medicine (Birkmeyer et al., 2013; Corallo et al., 2014; De Jong et al., 2006), consumer choice (Ortega et al., 2011; Scarpa et al., 2005), and law (Kang et al., 2012; Kleinberg et al., 2018; Arnold et al., 2018). In the context of medicine, this is referred to as the study of practice variation (Atsma et al., 2020; Cabana et al., 1999), where it is often observed that doctors, facing the same clinical context, will make different decisions. Likewise, in a legal context, judges often differ in their leniency in their decisions regarding bail (Kleinberg et al., 2018), juvenile incarceration (Aizer and Doyle Jr, 2013), the use of alternatives to incarceration (Di Tella and Schargrodsky, 2013), and incarceration length (Kling, 2006). In some scenarios this variation may be justified: The best medical treatment may not be obvious. In others, it may be grossly unfair, as in the case of racial bias in bail decisions (Arnold et al., 2018).
∗Equal contribution
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
In this work, we tackle the question of how to find and characterize this variation in the first place. In particular, we present a learning algorithm for identifying a “region of heterogeneity”, defined as a subset of all contexts (e.g., patients, cases) for which the identity of the decision-maker substantially affects the decision. In medicine, a better understanding of treatment variation can inform the development and dissemination of clinical guidelines. In the legal domain, characterizing the cases where judges vary most in their leniency may help with investigating potential issues of fairness.
We formalize characterizing the region of heterogeneity as a causal inference problem: We want to characterize examples where changing the decision-maker would have resulted in a different decision. The challenge is two-fold: First, we only observe a single decision-maker per example, so we cannot directly observe how (for instance) multiple judges would have decided the same case. Second, our data on individual decision-makers is often scarce. For instance, in Section 5, we consider a medical dataset with more than 400 doctors, each of whom has fewer than 9 samples on average.
We will refer to decision-makers as “agents”, and our contributions are as follows: In Section 2, we propose an objective defined in terms of counterfactual decisions across different agents, and show that this objective can be identified from observational data. Moreover, this objective does not require the use of agent-specific statistical models, making it amenable to our sparse setting. In Section 3, we give an iterative algorithm to identify regions of disagreement by maximizing this objective and provide intuition (in the form of a generalization bound) for the factors that drive its performance. In Section 4, we use a semi-synthetic dataset, derived from crowd-sourced recidivism predictions, to demonstrate that our algorithm recovers the correct region of heterogeneity accurately, even when there are many agents. Finally, in Section 5, we apply our algorithm to a real-world healthcare dataset and confirm that it recovers intuitive regions of variation in first-line diabetes treatment. We conclude with a discussion of related work and implications. Our code is available at https://github.com/clinicalml/finding-decision-heterogeneity-regions.
Our algorithm does not determine whether variation is inherently good or bad or how it should be addressed. Rather, more careful study with domain experts would be required to determine if variation can (or should be) reduced and how. In addition, false discovery of variation is possible and could have a negative impact. We expect that validation on independent datasets would be required in real-world applications, using the regions identified by our method as plausible hypotheses to test.

heading: 2 Characterizing Heterogeneity from a Causal Perspective

heading: 2.1 Notation
text: Let the data be drawn from a distribution P(X, A,Y), where X is a random variable representing context (or features), A is a discrete agent, and Y ∈ {0, 1} is the binary decision. The spaces of all X and A are denoted as X andA with realized values as lower case x and a, respectively. Indicator variables 1 [·] are one if the statement inside the brackets is true and zero otherwise. For a subset S ⊆ X, 1 [x ∈ S ] is sometimes written as a function S (x), where S : X → {0, 1}. A subset S may have several disjoint regions. E [Y |X ∈ S ] denotes the average Y across samples in S . For instance, if S = {X : X0 < 10}, then E [Y |X0 < 10] is a scalar average of Y among samples with X0 < 10.

heading: 2.2 Heterogeneity as a Causal Contrast
text: Our conceptual goal is to identify a region S ⊆ X where different agents tend to make different decisions even when faced with the same context. We can formalize this in the language of potential outcomes from the causal inference literature (Pearl, 2009; Hernán and Robins, 2020), which for clarity we will refer to as potential decisions: In particular, we denote Y(a) to be the potential decision of agent a. The fundamental challenge of causal inference is that we do not observe all potential decisions {Y(a) : a ∈ A} for each sample, but only a single decision Y . With this in mind, we will make the following assumptions, standard in the literature on causal effect estimation.
Assumption 1 (Causal Identification Assumptions). (i) Consistency: Y = y, A = a =⇒ Y(a) = y, and (ii) No Unmeasured Confounding (NUC): For all a ∈ A, Y(a) ⊥ A | X.
Consistency links the potential Y(a) to the observed Y , and NUC says that there are no unobserved factors that influence both the assignment of agents and the decision itself. For instance, the quasirandom assignment of cases to judges conditioned on features X satisfies NUC (Kleinberg et al., 2018).
NUC may be violated if key aspects of the case (e.g., misdemeanor vs. felony) are omitted as features. For instance, misdemeanor and felony cases may be seen by different judges and have different decision processes, but this variation is not due to agent preferences. Given these assumptions, we propose a causal measure of agent-specific bias, defined as a contrast between potential decisions. Definition 1 (Conditional Relative Agent Bias). For an agent a ∈ A and a subset S ⊆ X, the conditional relative agent bias is defined as
E[Y(a) − Y(π(x)) | A = a, X ∈ S ] (1) where Y(a) is the potential decision of agent a, and Y(π(x)) B ∑ a′ E[Y(a′) | x]π(a′ | x) denotes the expected potential decision under the agent assignment distribution π(a′ | x) B P(A = a′ | X = x).
Note that under Assumption 1, Y(π(x)) = E[Y | x],2 but here we emphasize its causal interpretation as the expected decision of a random agent. Equation (1) represents the relative difference between the decision of an agent (on their particular distribution of cases in the region) and the potential decision of a random agent. In particular, Equation (1) can be written as follows under Assumption 1
E[Y(a) − Y(π) | A = a, X ∈ S ] = ∫
x∈S E[Y(a) − Y(π) | X = x]p(x | A = a, X ∈ S )dx, (2)
where we shorten Y(π(x)) to Y(π). This is the average difference (over p(x | a), restricted to those x in the set S ) of the conditional expected difference between Y(a) and Y(π). For example, suppose that the agent a is a judge who is particularly lenient on bail decisions for felony arrests (the region S ), and Y = 1 denotes granting bail. Then imagine the following counterfactual: Take the felony cases that are assigned to this judge and reassign each individual case, described by x, to a random judge a′, proportionally to p(a′ | x). We may then observe, on average, that the bail rate would decrease, because most judges are less lenient than judge a, corresponding to a positive value of Equation (1).
Equation (1) has the additional advantage of being easy to estimate: Under Assumption 1, it can be rewritten3 as E[Y−E[Y | X] | A = a, X ∈ S ], the expected residual in predicting (using the conditional expectation E[Y | X]) the decision of an agent a across the context x typically seen by that agent.

heading: 2.3 Formalizing a Causal Objective
text: Our primary goal is to discover a region S where substantial heterogeneity exists across agents. To do so, we define an aggregate objective across a group G of agents, where G(a) ∈ {0, 1} is an indicator function for membership.
Q(S ,G) B ∑
a:G(a)=1
P(A = a | X ∈ S )E[Y(a) − Y(π) | A = a, X ∈ S ], (3)
We now show that this quantity can be identified and estimated from observational data without requiring agent-specific statistical models, before discussing the interpretation of this objective. Theorem 1 (Causal Identification). Under Assumption 1, Q(S ,G) can be identified as
Q(S ,G) = ES [Cov(Y,G | X)] = ES [(Y − E[Y | X])G], (4) where ES [·] B E[· | X ∈ S ] and Cov(Y,G | X) is the conditional covariance.
Theorem 1 and other theoretical results are proven in Appendix A. The result follows from proving that the agent-specific bias (Definition 1) is identifiable using the expected conditional covariance between Y and the binary indicator 1 [A = a]. With this in mind, we optimize the following objective, where the set S is constrained to be at least a certain size β and S is a hypothesis class of functions S .
max S∈S,G Q(S ,G) s.t., P(S ) ≥ β, (5)
Interpretation: Intuitively, this objective measures the disagreement between the agents in the group G(a) = 1 and the overall average E[Y | X] on the region S . Hence, the choice of group is important for interpreting the objective: If G(a) = 1 for all agents, the objective will be zero for any set S , as can be seen from Equation (4), applying the definition of the conditional expectation.
2See Proposition A1 for a short proof, and Proposition A2 for the derivation of Equation (2). 3See Proposition A3 in Appendix A.1.
Accordingly, we seek a region S for which the partially maximized objective L(S ) B maxG Q(S ,G) is large: This partial maximization is obtained by taking G(a) = 1 whenever the conditional relative agent bias of agent a (on the set S ) is non-negative. Thus, Equation (5) can be re-written as
max G Q(S ,G) = ∑ a∈A P(A = a | X ∈ S ) |E[Y(a) − Y(π) | A = a, X ∈ S ]|+ , (6)
where |x|+ B max(x, 0), and this objective becomes an average over agents who have a positive bias. This population objective is also equivalent (up to a constant factor) to the (weighted) sum of the magnitude of each agent’s conditional relative agent bias. See Proposition A4 in Appendix A.1.
Lack of Overlap: We have not made the overlap or positivity assumption that P(A = a | x) > 0 for all x, a. While this assumption is required to identify conditional causal effects E[Y(a) − Y(a′) | X] (Nie and Wager, 2017; Wager and Athey, 2018; Shalit et al., 2017), it is not required for identifying our causal contrast. Our problem only requires each context has a positive probability of being seen by more than one decision maker. For instance, suppose that S contains both misdemeanors and felonies and there are four judges a0, . . . , a3. If judges a0 and a1 make bail decisions exclusively for felonies while judges a2 and a3 make bail decisions exclusively for misdemeanors, our measure captures disagreement between a0 and a1 and between a2 and a3 even though comparisons between a0 and a2 or other pairs are impossible to make. Thus, we have chosen to compare Y(a) to the decisions of viable alternative agents, weighted by their probability p(a′ | x) of being selected for such a case.

heading: 3 Identifying Regions with Heterogeneity
text: In Section 3.1, we introduce an iterative optimization algorithm for a finite sample version of Objective (5) that alternatingly optimizes S and G. In Sections 3.2 and 3.3, we discuss practical heuristics for choosing the region size parameter β on training data and validating if the resulting region generalizes to held-out data. Finally, we build intuition for the factors that influence performance of this algorithm via a generalization bound in Section 3.4 under simplifying assumptions.

heading: 3.1 Iterative Optimization Algorithm
text: We let Q̂(S ,G) be the empirical analog of Q(S ,G) (Equation 4), which we can write as follows,
Q̂(S ,G) B 1∑ a, j 1 [ xa j ∈ S ] ∑ a, j (ya j − f (xa j)) ·G(a) · 1 [ xa j ∈ S ] . (7)
where f (x) denotes a model of the conditional expectation f (x) ≈ E[Y | X = x]. For simplicity of notation, we assume that there are R samples (indexed by j) for each of a finite set of N agents (indexed by a), giving N · R samples in total. Our algorithm (Algorithm 1) takes as input the data {(xa j, ya j)} and a minimum region size β, and outputs a model h(x) and a threshold value b that describe a region of heterogeneity S = {x ∈ X; h(x) ≥ b}. Starting with S = X (the entire space), the algorithm identifies the grouping that maximizes Q̂(S ,G), then uses that grouping to identify the region maximizing the same quantity, repeating this process until convergence. The algorithm uses a classifier f (x) to estimate E[Y | X = x] and a regression model h(x) to estimate the conditional covariance of the decision Y and the grouping G at X = x. Note that we can use any supervised learning algorithms for f and h, allowing us to learn interpretable regions as part of the algorithm if h(x) is interpretable (e.g., decision trees). If sufficient data is available, samples can be split into three parts for estimating f (x) in line 2, computing G(a) in lines 5-8, and training h(x) and estimating the (1 − β)-th quantile in line 10. We do not perform this sample splitting because our sample sizes are already small. Under-fitting f (x) by further restricting the sample size could lead to false discovery if f (x) does not capture the variation explained by X.
Optimizing over G given S . Given a region S , our first result identifies the grouping G : A → {0, 1} that maximizes Q̂(S ,G) and shows that it can be expressed in terms of Q̂(S , 1 [A = a]). Proposition 1. Given S ⊆ X, Q̂(S ,G) is maximized over the space of functions G : A → {0, 1} at GS , where GS (a) = 1 [ Q̂(S , 1 [A = a]) ≥ 0 ] .
Intuitively, this proposition states that to maximize the empirical expected conditional covariance of the decision and grouping on a region, we must group agents by whether their residuals ya j − f (xa j) are (on average) positive or negative on S .
Algorithm 1 Identifying regions with variation
1: Input: Data {{xa j, ya j}Rj=1} N a=1 , minimum region size β. 2: Fit a model f (x) to E(Y | X = x). 3: Initialize S = X. 4: repeat 5: for a = 1, . . . ,N do 6: Compute Q̂(S , 1 [A = a]), where Q̂(S , 1 [A = a]) B 1∑ j 1[xa j∈S ] ∑ j(ya j − f (xa j))1 [ xa j ∈ S ] , 7: Set G(a) = 1 if Q̂(S , 1 [A = a]) ≥ 0 and 0 otherwise. 8: end for 9: Compute ba j = (ya j − f (xa j))G(a), a = 1, . . . ,N, j = 1, . . . ,R.
10: Fit a model h(x) to predict ba j from xa j, and let b be the (1 − β)-th quantile of h(xa j). 11: S ′ ← S . 12: S ← {xa j; h(xa j) ≥ b}. 13: until S = S ′ or iteration limit reached. 14: Output: Model h and threshold b, defining a region S = {x ∈ X; h(x) ≥ b}.
Optimizing over S given G. To optimize Q̂(S ,G) for a fixed grouping G over the hypothesis class S, we train a model h(x) to predict (ya j − f (xa j))G(a) given xa j, where h ∈ H . Using h as an estimate in Eq. 7, we find a set S to maximize the quantity 1∑ x 1[x∈S ] ∑
x∈S h(x), the empirical conditional expectation of h(x) over S . This quantity is maximized (subject to our β constraint) by including the largest β-fraction of the h(xa j) in S . Hence, we pick b as the (1 − β)-th quantile of h(xa j) and choose our region as Ŝ G = {x ∈ X; h(x) ≥ b}.

heading: 3.2 Tuning the Region Size Parameter
text: For real datasets, we need to choose β without knowledge of the true value. Given that our objective can be calculated on held-out data using the functions S ,G, a seemingly obvious approach would be to compute Q(S ,G) on a validation set and select the β that leads to the highest Q(S ,G). However, for a fixed data distribution, smaller values of β tend to produce higher values of Q(S ,G), and there is a trade-off between finding a smaller region of higher variation and a larger region that may include areas of lower (but still meaningful) variation. This motivated our original constraint P(S ) ≥ β. To select β, we propose a heuristic inspired by permutation-based hypothesis testing (Wasserman, 2004). We compare the training objective to a reference distribution of values (for the same β) that might be seen even if all agents followed the same policy. For each candidate β, we (i) run our algorithm and compute the objective on training data qobs B Q̂(Ŝ , Ĝ). (ii) For T iterations, we randomly shuffle the agents and re-run the algorithm to get a new objective value. This gives us a distribution P̂null over Q(S ,G) from a distribution where P(X,Y) and P(A) are unchanged but Y, X ⊥ A. (iii) Finally, we compute a p-value pβ = P̂null(Q > qobs) and choose the β with the smallest p-value. In Section 4.2, we find that this heuristic empirically recovers the true β value in semi-synthetic settings.

heading: 3.3 Validation of the Region
text: We may wish to validate the learned region Ŝ independently of the grouping Ĝ. In particular, finding Ĝ is not our main goal, and we observe in our semi-synthetic experiments that our algorithm can find the true region S even when the grouping Ĝ is fairly poor (due to few samples per agent), as shown in Appendix B.2. We can optimize over G in Q(S ,G) to obtain an objective that depends only on S and can be used to compare regions. By Proposition 1, we obtain an empirical analog of Equation (6) as
L̂(S ) B max G
Q̂(S ,G) = 1∑ a, j 1 [ xa j ∈ S ] ∑ a ∣∣∣∣∣∣∣∣ ∑ j (ya j − f (xa j))1 [ xa j ∈ S ]∣∣∣∣∣∣∣∣ + , (8)
where |x|+ is equal to the positive part [x]+ = max(x, 0) as before. We then use this objective L̂(S ) to answer the following question: Does our chosen region Ŝ yield a significantly higher objective value on test data than a randomly selected region of the same size? An example of this analysis is given in Table 1 for the real-data experiment in Section 5.

heading: 3.4 Generalization Error
text: We give a generalization bound for Algorithm 1 to build intuition for the factors that influence performance. To derive this bound, we consider a simplified setting, where there exists a set S ′,G′ such that the following set of assumptions hold. Assumption 2 (Group-based variation). For all x ∈ S ′, E[Y | X = x, A = a] = E[Y | X = x,G′(a)] and for all x < S ′, E[Y | X = x, A = a] = E[Y | X = x] Assumption 3 (Non-zero relative biases). There exists a constant α > 0 such that for all x ∈ S ′, E[Y | X = x,G′(A) = 1] − E[Y | X = x] > α, and E[Y | X = x,G′(A) = 0] − E[Y | X = x] < −α, Assumption 4 (All agents see samples in S ′). There exists a constant ω > 0, such that for every a ∈ A, P(X ∈ S ′ | A = a) > ωP(X ∈ S ′).
Note that under these assumptions, S ′,G′ maximize the objective Q(S ,G) (see Appendix A.4), so we will refer to them as S ∗,G∗ for the remainder of this section. Assumption 2 says that there are two groups of agents, who follow two distinct decision policies within a region S ∗ but follow an identical decision policy outside of S ∗. Assumption 3 says that one group has a positive bias across all of S ∗, relative to the average over both groups, and the other group has a negative bias. To simplify the analysis, we also make Assumption 4 that every agent has some non-zero chance of observing some contexts X in the region, but note that we do not require that p(x | a) > 0 for all x ∈ S ′. Under these assumptions, we demonstrate that the first iteration of Algorithm 1 will find, with high probability, a region Ŝ whose value Q(Ŝ ,G∗) (for the same grouping G∗ defined above) is close to that of the optimal S ∗. Note that we do not claim that the iterative algorithm finds the globally optimal solution. For simplicity, we assume that f (x) perfectly recovers E[Y | X]. This can be relaxed at the cost of additional terms in the bound that go to zero as the overall sample size increases. Under Assumptions 2, 3, and 4, assume that P(S ∗) = β. For the informal version presented here, we assume that P(S ∗) = P(Ŝ ) = β, where Ŝ is returned by our algorithm, and that exactly a β-fraction of our samples fall into S ∗ and Ŝ (in the Appendix, we give a version without these simplifications). Theorem 2. Under the assumptions above, if S ∗ ∈ S and R > 2 ln 2
α2β2ω2 , the first iteration of Algorithm 1
returns Ŝ such that, with probability at least 1 − δ, Q(S ∗,G∗) − Q(Ŝ ,G∗) ≤ , where
= √ 2 ln(3/δ) βN · R + 2 β η + √3η(1 − η)δ · N  + 1β 2R(S,N · R) + 4 √2 ln(12/δ)N · R  ,
where R(S,N · R) is the Rademacher complexity of S, and η B exp ( −Rα2β2ω2/2 ) .
The term η plays an important role: It bounds the expected misclassification error P(Ĝ(a) , G∗(a)). For sufficiently large R, we have that η < 1/2 with high probability, i.e., we have a better-than-random chance of identifying the group for an individual agent. Moreover, η decreases as we increase the number of samples R for each agent, the separation α between the two groups on S ∗, the region size β, and the constant ω. For sufficiently small η, our algorithm discovers a region whose value (in terms of Q(S ,G∗)) is close to that of the optimal region. The generalization bound improves as the number of agents N increases, the number of samples R for each agent increases, or the complexity of the hypothesis class decreases. The latter is measured here by the Rademacher complexity R(S,N · R) of our hypothesis class S, which can be bounded by standard arguments. In conclusion, under some additional assumptions, Algorithm 1 identifies an approximately optimal solution with high probability after one iteration. We show via semi-synthetic experiments in Appendix B.3 that convergence is generally fast in practice.

heading: 4 Semi-Synthetic Experiment: Recidivism Prediction
text: For conceptual motivation in the introduction, we discussed the legal system: As a potential application of our method, one could determine types of cases for which the idiosyncratic preferences of judges have a significant impact on their decisions. Lacking data on judge decisions with sufficient context, we turn to the more controlled setting of human predictions of recidivism.
Dataset: We use publicly available data from Lin et al. (2020),4 who ask participants on Amazon’s Mechanical Turk platform to make recidivism predictions based on information present in the
4Available at https://github.com/stanford-policylab/recidivism-predictions
“Correctional Offender Management Profiling for Alternative Sanctions” (COMPAS) dataset for Broward County, FL (Dressel and Farid, 2018). Participants (or “agents”) are shown 5 risk factors: age, gender, number of prior convictions, number of juvenile felony charges, and number of juvenile misdemeanor charges. The charge in question is also given, as well as whether the charge is a misdemeanor or felony. The dataset contains 4550 cases evaluated by 87 participants.
Semi-Synthetic Policy Generation: To benchmark our method, we generate semi-synthetic data where we have access to a “ground truth” region of heterogeneity. We retain the features presented to the original participants and construct two policies, which we refer to as the “base” and “alternative” policies: For the base policy, we learn a logistic regression model on the binary decisions across the whole dataset. For the alternative policy, an extra positive term is added to the logistic regression for samples within the region. We construct two scenarios with different regions of variation: (1) all drug possession charges, and (2) all misdemeanor charges where the individual is 35 years old or younger. These make up 22% and 21% of the data, respectively. Then, we generate synthetic agents (randomly assigned to cases) and assign half of the agents to the base policy and half to the alternative. Synthetic decisions are then sampled from the logistic regressions. For each scenario, the two groups of agents follow the same stochastic policy outside of the region, and one group systematically prefers Y = 1 within the region. More details can be found in Appendix B.1.

heading: 4.1 Performance versus Baselines
text: Baselines: We compare how well our approach identifies the true region of heterogeneity with several baselines. To our knowledge, the problem of finding regions of heterogeneity with a large number of agents has not been studied before. Many causal inference methods for treatment effect estimation are designed for a single, binary treatment. However, naively estimating the treatment effect between each pair of providers would scale O ( |A|2 ) . Therefore, we develop new baselines. Some (including
the causal forest and U-learner adaptations described in Appendix B.2) are based on causal inference methods augmented to identify a region of heterogeneity and grouping of agents where possible.
Direct models: This baseline measures how much adding the agent to the feature set improves prediction of decisions. We fit logistic regressions with and without the agent feature to estimate E[Y | A, X] and E[Y | X]. For each sample (x, y, a), we compute |y − E[Y | X = x]| − |y − E[Y | A = a, X = x]| to quantify how much the model with agents outperforms the model without agents. Then, we fit a “region model” to predict this quantity from X. This region model is either a ridge regression, decision tree, or random forest model. Finally, we compute the top β quantile of predictions from the region model in the training and validation sets and use this cut-off to select points in the test set.
TARNet: A treatment-agnostic representation network (Shalit et al., 2017) models the outcomes of all treatments for each sample by learning a shared representation of the sample features and then
having a separate prediction head for each treatment. We implement the shared representation model using a 2-layer neural network with ReLU and dropout. Each prediction head is a linear layer with a sigmoid function. TARNet predicts E[Y | X, A] for every agent for each sample. VarA[E[Y | X, A]] measures the variation across all agents if they had seen context X and is analogous to our objective without grouping. We predict this quantity with the region models as in the direct model baseline.
Results: We evaluate how well the algorithms identify the samples within the region of heterogeneity when we vary the number of agents among 2, 5, 10, 20, 40, and 87, where 87 is the number of agents in the original dataset. Figure 1 shows the best overall region models for each set-up, with the other models deferred to Appendix B.2. The metric in Figure 1 is the region AUC, defined as how well the model classifies whether samples belong in the region of heterogeneity when compared to the true region. Algorithm 1 consistently performs well for both semi-synthetic set-ups, especially when the number of agents increases to a realistic level (and the number of samples per agent decreases). The direct baseline deteriorates very rapidly as the number of agents increases in the drug possession set-up, while the TARNet baseline deteriorates rapidly in the misdemeanor under age 35 set-up. Refer to Appendix B.2 for additional baseline details, region models, and evaluation metrics. We also show that our method is robust in a set-up with more than 2 agent groups in Appendix B.4.

heading: 4.2 Tuning the Region Size
text: We validate the proposed approach of tuning β (discussed in Section 3.2), by applying the methodology to our semi-synthetic setting here. We sample 30 semi-synthetic datasets and consider candidate values of β in [0.02, 0.42] in increments of 0.04. For each proposed value of β we use T = 40 random permutations of the agents. Figure 2 presents results for the misdemeanor under age 35 set-up. As the candidate value of β increases (up to the true value of β), the p-value decreases, and the distribution of selected β values are centered on the true value of β.

heading: 5 Real-data Experiment: First-Line Diabetes Treatment
text: We apply our algorithm on a real-world dataset consisting of first-line (initial) treatment for type 2 diabetes and examine how the treatment variation we discover aligns with clinical knowledge. We present an additional real-world experiment (on Parkinson’s disease) in Appendix D.
Data and Setup: We use an observational, de-identified, dataset provided by a large health insurer in the United States, spanning from 2010 to 2020. The task is to classify first-line treatment decisions between metformin (Y = 0)–the typical recommendation from the American Diabetes Association– and other common first-line treatments such as sitagliptin, glipizide, glimepiride, or glyburide (Y = 1) (American Diabetes Association, 2010; Hripcsak et al., 2016). As relevant clinical features, we include the patients’ most recent eGFR (mL/min/1.73m2) and creatinine (mg/dL) measurements,
incidence of heart failure, and treatment date. Because treatment date does not define a type of patient, we omit it from the region model. However, including it in the outcome model is essential because of increasing use of metformin over time. The agent indicator A is the group practice of the doctor responsible for the patient’s treatment. 3,980 patients and 447 group practices are included in our cohort. After requiring at least 4 patients per agent, 3,576 patients and 176 group practices are included. This filter ensures each group practice has at least 1 sample in the training and validation sets and at least 2 samples in the test set. In this experiment, we choose β = 0.25 as input to our algorithm. See Appendix C for additional cohort definition and set-up details.
Interpretation of Results: To interpret the region, we use decision trees as our region model h(x). The tree is visualized in Appendix C. The decision tree identifies the region of heterogeneity as the union of (i) eGFR below 71.5 and (ii) eGFR above 98.5 and creatinine above 0.815. These regions align with clinical intuition. In the first region, low eGFR values indicate impaired renal function (Group et al., 2009), which is a contraindication for metformin since it is traditionally associated with increased risk of lactic acidosis (Tahrani et al., 2007). Still, treatment decisions can vary here because guidelines for managing patients with eGFR below 45 are lacking (Group et al., 2015). Note that this region provides an example of how our algorithm works when overlap is not satisfied. Although 33 of 176 group practices do not see patients with these features, we can still conclude that this is a region of heterogeneity among the 143 agents with cases. In the second region, there are no obvious contraindications for metformin. Thus, understanding why some doctors on average only prescribe metformin 78% of the time to patients in this region may help us identify whether this is a region in which we can standardize practice.
Assessing Significance: In Table 1, we perform a sanity check, assessing whether our algorithm identifies a region S whose variation in held-out data is higher than that of a randomly selected region, using the partially optimized objective L(S ) = maxG Q(S ,G) laid out in Section 3.3 to compare regions directly. Table 1 shows that L(Ŝ ) is similar on the training and test data. Furthermore, we compare L(Ŝ ) on the test set to the distribution of L(S rand), where S rand are random subsets of the test data of the same size as Ŝ . We compute the latter distribution using 100 random subsets and observe that the test value of L(Ŝ ) is more than two standard deviations above the mean of the latter. This gives us confidence that the discovered region S generalizes to a region of heterogeneity beyond the training set. We direct the reader to Appendix C for additional analyses, such as evaluating the stability of the region over multiple folds of splitting the data.

heading: 6 Related Work
text: Beyond previously mentioned connections to causal effect estimation, we highlight a few areas of research that share similar goals to our own.
Agent-specific Models of Decisions: Prior works have modeled agent-specific decision-making processes by estimating a separate model for each agent. Abaluck et al. (2016) model heterogeneity in physician tendency to run diagnostic tests. Chan Jr et al. (2019) estimate radiologist skill based on diagnosis and miss rates. In our setting, unlike in diagnosis, there is no “correct” decision that can be incorporated into the model. Ribers and Ullrich (2020) assume there is provider-specific noise in determining patient type, which affects the pay-off functions for deciding whether to prescribe antibiotics. When only a few decisions are observed per agent, these agent-specific models cannot be estimated reliably. Currie and MacLeod (2017) also incorporate physician beliefs and skill into a pay-off function. They estimate an aggregate logistic choice model (for C-sections) across all physicians and then learn how individual physicians deviate from this model. They do not learn the regions where this deviation occurs, as they focus on how heterogeneity is associated with
downstream outcomes. Norris (2020) looks for disagreement between agents but relies on some cases being seen by multiple agents. We assume each case is seen by only one decision-maker.
Conditional Independence Testing: While our objective maximizes a causal notion of dependence, one could instead ask if Y is conditionally independent of A given X. Many metrics exist for testing conditional independence, such as the Hilbert-Schmidt independence criterion (HSIC) (Fukumizu et al., 2007; Zhang et al., 2012), conditional mutual information (Runge, 2018), conditional correlation (Ramsey, 2014), and expected conditional covariance (Shah and Peters, 2020). We give the last a causal interpretation under some assumptions and seek a region that maximizes it, in lieu of testing.
Hierarchical /Mixture Models: Bayesian methods are often used to estimate group-level effects, such as a per-physician effect on patient outcomes (Tuerk et al., 2008), where group identifiers are included as a categorical feature in a multi-level generalized linear model (Gelman and Hill, 2006). Alternatively, one could assume a conditional mixture model (Bishop, 2006), where agents belong to latent clusters that each have their own policy. However, both of these methods require parametric assumptions on the distribution of P[Y | x, a], and even so, the optimal mixture model is not necessarily identifiable when both the clusters and policies are unknown (Dasgupta and Schulman, 2007). By contrast, our method does not require making particular parametric assumptions about E[Y | x, a] and seeks to learn the region of heterogeneity directly. Feature Evaluation: Checking for heterogeneity can also be framed as feature evaluation, where we would like to evaluate whether adding the agent identifier will increase predictive power. Feature evaluation methods typically maximize dependence between selected features and labels, utilizing measures similar to those in conditional independence testing, such as the HSIC (Song et al., 2007). Other methods use the correlation of the new feature with the loss gradient of the current predictor as a measure of utility (Koepke and Bilenko, 2012). In contrast, we focus not on checking marginally for the predictive power of agent identifiers, but rather identifying a region.
Crowdsourcing: Our work is conceptually related to identifying which samples are difficult to label via crowdsourcing annotations (Karger et al., 2014; Whitehill et al., 2009). Most crowdsourcing models are generative with latent variables for the correct sample labels, sample difficulty, agent expertise, and agent bias. They then optimize for the likelihood of the observed labels. The set of difficult samples is analogous to our region of heterogeneity. The main difference with our problem is that we do not require any notion of the “correct” label.

heading: 7 Discussion
text: In this work, we take a causal perspective on finding regions where agents (i.e., decision-makers) have heterogeneous preferences, formalizing this heterogeneity in terms of counterfactual (or “potential”) decisions across agents. We propose a causal measure of agent-specific biases and give an objective that aggregates this bias over agents. This objective can be identified from observational data and written in terms of an expected conditional covariance. Importantly, for our applications, this does not require building agent-specific models or assuming overlap across all agents.
We give an iterative algorithm to find a region that maximizes this objective. Then, we demonstrate in semi-synthetic experiments that our algorithm accurately recovers the region of heterogeneity and scales well with the number of agents. In contrast, performance of baslines deteriorates when the number of agents increases. Although our experiments have low-dimensional spaces, we hypothesize our algorithm would scale well to high-dimensional feature spaces since the average policy and region models can handle high-dimensional input spaces. Finally, on a real-world medical dataset, we show that our algorithm can yield insights that align with clinical knowledge.
Our work is motivated by understanding variation in human decision-making. In the judicial domain, our method may help unearth types of cases where decisions are highly dependent on the judge assigned to the case. In the medical domain, our approach may identify types of patients where new guidelines may be required to help doctors make decisions more consistent with standard of care. Domain expertise is required to determine the implications of the regions discovered by our method. Beyond these applications, we see our approach as a useful data science tool for understanding heterogeneity in decisions that appears to be driven by individual-level preferences.
Acknowledgements: We would like to thank Aaron Smith-McLallen, James Denyer, Luogang Wei, Johnathon (Kyle) Armstrong, Neil Dixit, Aditya Sai, and the rest of the data science group at Independence Blue Cross, whose expertise, data, and support enabled the diabetes experiment. We would also like to thank Rebecca Boiarsky for converting the diabetes data to the OMOP common data model format, Monica Agrawal for her helpful comments, and other members of the lab for insightful discussions. We are also grateful to Charles Venuto and Monica Javidnia for their advice on Parkinson’s. This work was supported in part by Independence Blue Cross, Office of Naval Research Award No. N00014-17-1-2791, an Abdul Latif Jameel fellowship, and a NSF CAREER award.

Review summaries:
Summary:  This paper considers settings where decisions for individuals are each made by different decision makers (agents). The paper defines the problem of finding a region in feature space in which the agents administer highly varying decisions. For example, perhaps judges mete out different sentences for misdemeanors by young men. The paper defines new causal criteria to characterize such regions in the feature space, where the choice of agent has an effect on the decision. The paper introduces an iterative algorithm to maximize this criteria to identify regions of heterogeneity, and provides theoretical results about the optimality of the proposed algorithm. The paper empirically demonstrates the algorithm on a semi-synthetic dataset and real clinical data about doctor recommendations for diabetic patients.  Contributions:  + A new criteria for defining regions where there is heterogeneity in decisions administered across agents. + Results for estimating the criteria from observational data and an iterative algorithm for optimizing the criteria. + Theoretical results about the algorithm and empirical validation of the algorithm.

In this paper, the authors develop an algorithm for identifying cases (e.g., defendants or patients) on which decision makers (e.g., judges and doctors) "disagree" on the decision.   ## Identification The authors consider a setting in which cases are summarized by features $X$. Each case is associated with a decision maker $A$ and a binary decision $Y$. In the pretrial release setting, each case is a defendant, the features $X$ are observable characteristics of the defendant (such as demographics or charge information), the decision maker $A$ is the assigned judge and the binary decision $Y$ is to either release or detain the defendant. The authors define the potential decision Y(a), which describes the decision that would be made by decision maker $a$.   The authors propose a summary measure of agent disagreement using average contrasts of potential decisions across agents given sets of the observable features. In particular, the authors focus on the "conditional relative agent bias" as  $$ \mathbb{E}[Y(a) - Y(\pi(x) | A = a, X \in S], $$ where $Y(\pi(x))$ denotes the expected decision under the assignment mechanism of decision makers to cases. The authors key identification assumption is that decision makers are as-if randomly assigned to cases conditional on the features (i.e., $A \perp Y(a) \mid X$). This measures the extent to which decision maker a's average decisions over $X \in S$ differs from the average decision that is made under the assignment mechanism. In other words, this measure of disagreement depends on the assignment mechanism of agents.  The conditional relative agent bias measures the extent to which a single agent a's decision disagrees relative to the average decision over $X \in S$. The authors then propose aggregating this measure across agents. For some grouping of agents G, the authors define the estimand $$ Q(S, G) = \sum_{a : G(a) = 1} P(A = a \mid X \in S) \mathbb{E}[Y(a) - Y(\pi(x) | A = a, X \in S]. $$ The authors main identification result (Theorem 1) is to show that if agents are as-if randomly assigned to cases, then this estimand can be identified as the conditional covariance between the observed decision and the agent assignment $$ Q(S, G) = \mathbb{E}[Cov(Y, G \mid X ) \mid X \in S]. $$  ## Estimating Regions of Heterogeneity. Using this identification result, the authors goal is to search over regions of the feature space $S$ and groupings over agents $G$ to find the combination that maximizes this proposed measure of disagreement. To do so, the authors propose to solve  $$ \max_{S, G} Q(S,G) \mbox{ s.t. } P(X \in S) \geq \beta, $$ which requires that the region have probability at least $\beta$ (a tuning parameter to be specified). The authors develop an iterative algorithm to optimize this objective function, and provide some conditions under which they are able to derive a generalization bound on the solution return after the first iteration (Assumptions 2-4 and Theorem 2).  ## Empirical Applications The authors apply their proposed algorithm in two empirical applications.   First, the authors analyze the behavior of their proposed algorithm in a semi-synthetic experiment that is calibrated to data from Lin et al. (2020) on recidivism risk predictions made by MTurk workers. The authors simulate data from two types of agents that follow different decision rules for detention vs. release. The authors compare their algorithm against direct models that fit logistic regressions for the decision on agent identifiers and observable features (direct models) and TARNet. The authors report that their algorithm does better at recovering the "region AUC."  Second, the authors analyze an observational dataset on medical treatment decisions for diabetes. The data contain information on 3,956 patients and 458 group practices, and the authors search for the region of disagreement across these group practices. 

The authors are concerned with identifying regions in a covariate space. A region is characterised by having high variation in (potential) treatment outcome, as compared to the expected (potential) outcome. A method such as this could be used to identify variation of practice in the legal or medical domain (as indicated by the authors).   This is a relevant topic which the authors rightfully relate to the causal literature. Specifically, the authors chose the potential outcomes model, where the treatment is now a clinician or judge; the outcome is binary (amounting to a decision, e.g. bail or no bail); and the condition is an individual. Contrasting potential outcomes, the authors are interested in learning subsets of X, rather than Y(a).

This paper proposes a new method to identify "regions of disagreement" in a setting with human decision-makers assigned to individuals (eg. judges assigned to court cases). The "regions" are framed as subgroups of the data (eg. males over 30 years old) which can be interpretable (based on the chosen model). The paper formalizes this task as that of identifying the subgroup where the identity of the specific decision-maker has a causal effect on the decision, proposes an optimization algorithm to find this subgroup, develops generalization bounds (under some simplifying assumptions), and demonstrates compelling performance on semi-synthetic and real-world datasets.