Fair Regression for Health Care Spending

Risk adjustment is a method for correcting payments to health insurers such that they reflect the cost of their enrollees relative to enrollee health. It is implemented by most federally regulated health insurance markets in the United States, including Medicare Advantage and the individual health insurance Marketplaces created by the Affordable Care Act, to prevent losses to insurers who take on sicker enrollees (Pope et al., 2004, McGuire et al., 2013, Kautter et al., 2014). Current risk adjustment formulas use ordinary least squares (OLS) linear regression to predict health plan payments with select demographic information and diagnosis codes from medical claims. These OLS-based formulas are then typically evaluated with overall measures of statistical fit, such as $R^{2}$.

While $R^{2}$ is an important benchmark for evaluating global fit, it lacks information on other dimensions. As a result, risk adjustment has been criticized for not incentivizing efficient payment systems, spending, or population health management (Ash and Ellis, 2012, Layton et al., 2017), and for poorly estimating health costs for some groups by underpredicting their spending relative to average observed spending in the group. Underpredicting spending leads to undercompensation to the insurer, and there is evidence that insurers adjust the prescription drugs, services, and providers they cover (i.e., benefit design) to make health plans less attractive for enrollees in undercompensated groups (Shepard, 2016, Carey, 2017, Geruso et al., 2017). Examples of undercompensated groups include enrollees with specific medical conditions, high-cost enrollees, and partial-year enrollees (van Kleef et al., 2013, Montz et al., 2016, Ericson et al., 2017). Recent research has also shown that health plan insurers have the ability to identify undercompensated groups (Jacobs and Sommers, 2015, Geruso et al., 2017, Rose et al., 2017, Withagen-Koster et al., 2018).

What constitutes a fair or unfair algorithm depends heavily on the context. These fairness concepts and methods have been largely developed in the computer science literature (Chouldechova and Roth, 2018). We will consider risk adjustment formulas unfair if they underpredict spending for a prespecified group of enrollees, which then incentivizes differential treatment for the group via benefit design due to this undercompensation. For example, if average observed spending for individuals with mental health and substance use disorders (MHSUD) is $\$$10,000, but average predicted spending in this group is $\$$8,000, the risk adjustment formula may be unfair for the MHSUD group by ‘substantially’ underpredicting their spending. We define formal metrics for evaluating fairness in risk adjustment formulas using group residual errors in the next section.

Methods for addressing fairness are often separated into three categories based on the point in the learning process at which fairness is addressed: the preprocessing, fitting, or postprocessing phase. If the data are inherently biased, then preprocessing techniques are a possible solution. These methods create fair datasets by transforming or changing the data so that it is no longer biased (Kamiran and Calders, 2009, Zliobaite et al., 2011, Zemel et al., 2013, Calmon et al., 2017, Johndrow and Lum, 2017). It has been shown that current spending patterns among various group may be undesirable due to the plan benefit system, and by using observed spending data, we reinforce these unfair spending patterns. A recent study explored this concept by transferring funds to undercompensated groups in the raw data in order to promote more ideal spending patterns (Bergquist et al., 2019).

One of the most common fitting phase approaches in risk adjustment attempts to fix group undercompensation by adding new variables representative of the groups in the risk adjustment formula (van Kleef et al., 2013). While this is a straightforward idea, it can be problematic if those variables are unavailable, incentivize over- or underutilization of health services, or the risk adjustment formula does not recognize the improvement (Rose and McGuire, 2019). Fitting techniques in fairness include separate formulas for protected classes as well as fairness penalty terms or constraints (Kamishima et al., 2012, Berk et al., 2017a, Zafar et al., 2017a, b, Bechavod and Ligett, 2018, Dwork et al., 2018). We see intersections of these areas in the risk adjustment literature with separate formulas for enrollees with MHSUD (Shrestha et al., 2018, van Kleef et al., 2018) and constrained regression to reduce undercompensation for specific groups (van Kleef et al., 2017). Notably, separate risk adjustment formulas are already used in practice for infants and adults due to known differences in spending patterns. Nonparametric statistical machine learning methods to enhance estimation accuracy in risk adjustment have also been explored for the fitting stage (Rose, 2016, Shrestha et al., 2018, Park and Basu, 2018), but none of these tools are currently deployed in the U.S. health care system.

Postprocessing techniques modify the results after fitting by, for example, creating specific classification thresholds for different groups (Bansal et al., 2014, Hardt et al., 2016, Kleinberg et al., 2018, El Mhamdi et al., 2018). These methods separate fit from fairness objectives and allow use of the same prediction function for multiple fairness objectives. Reinsurance, paying insurers for a portion of the costs of high-cost enrollees, can be considered postprocessing for risk adjustment in that it reduces undercompensation for high-risk enrollees (McGuire and van Kleef, 2018).

In this paper, we focus on the fitting phase and expand on concepts from statistics, computer science, and health economics, proposing new estimation methods and measures to improve risk adjustment formulas for undercompensated groups. We develop two new fair regression estimators for continuous outcomes that reduce residual errors for an undercompensated group by building fairness considerations directly into the objective function. We also extend a definition of fairness from the computer science and statistics literature for the risk adjustment setting while additionally considering existing measures.

Our application features the IBM MarketScan Research Databases. This set of databases contains enrollee-level claims, demographic information, and health plan spending for a sample of individuals (and their dependents) insured by private health plans and large employers across the country. In 2014, the IBM MarketScan Research Databases were used by the federal government to develop the risk adjustment formula for the individual health insurance Marketplaces. Thus, this data source is particularly policy relevant. The undercompensated group we focus on for this data application is enrollees with MHSUD. We select this group for two major reasons. First, individuals with MHSUD are known to have substantially undercompensated payments in current risk adjustment formulas (Montz et al., 2016). Second, about 20% of people in the United States have MHSUD, thus it is a priority area for policy change. Although the data are representative of only a subset of the U.S. health insurance market, our methods are appropriate for other markets and different application settings with continuous outcomes. The methods and metrics we present are compared in this data analysis as well as simulation studies.

This section describes our approach to fair regression. It involves a suite of fairness measures for evaluating new and existing regression tools in an effort to improve risk adjustment formulas for undercompensated groups. A typical algorithmic fairness problem has an outcome $Y$ and input vector $\boldsymbol{X}$ that includes a protected group $A\subset\boldsymbol{X}$. The goal is to create an estimator for the function $f(\boldsymbol{X})=Y$ that maps $\boldsymbol{X}$ to $Y$, while aiming to ensure that the function is fair for protected group $A$. Although our main goal is to understand whether estimation methods beyond OLS, including those we newly propose, improve fairness for risk adjustment, we also wish to focus on interpretability for stakeholders, such as government agencies, insurers, providers, and enrollees. Therefore, constrained and penalized regressions were natural choices to enforce fairness in risk adjustment for undercompensated groups.

We selected a random sample of 100,000 enrollees from the IBM MarketScan Research Databases. Age, sex, and diagnosed health conditions, all from the year 2015, were used to predict total annual expenditures in 2016. Diagnosed health conditions took the form of the established Hierarchical Condition Category (HCC) variables created for risk adjustment. HCCs were developed by the Department of Health and Human Services to group a selection of International Classification of Disease and Related Health Problems (ICD) codes into indicators for various health conditions (Pope et al., 2004, Kautter et al., 2014). We considered the 79 HCC variables currently used in Medicare Advantage risk adjustment formulas and retained the 62 HCCs that had at least 30 enrollees with the condition. See Web Appendix C for a list of the 62 HCCs included in the regression formulas. Our sample of enrollees was 52% female and between the ages of 21 and 63, with median age 45. Mean and median annual expenditures per enrollee were $6,651 and $1,511, respectively.

We defined enrollees with MHSUD, our protected group $A$, using Clinical Classification Software (CCS) categories. This classification system maps each MHSUD-related ICD code to a CCS category, unlike the HCCs, which only map a subset of MHSUD-related ICD codes. Based on CCS categories, 13.8% of the sample had a diagnosis code for MHSUD compared to 2.6% had we used HCCs. We note that we do not capture enrollees with MHSUD who do not have an ICD code for their condition(s). The mean annual expenditures for MHSUD enrollees in our sample were $11,520 versus $5,880 for enrollees without MHSUD (and $3,744 versus $1,274 for median annual expeditures).

We compared each method to determine which estimators were best at reducing undercompensation for enrollees with MHSUD, and at what cost to overall statistical fit. In Table 1, we report the top estimators with respect to fairness for each of the six methods, having selected the hyperparameter value that optimizes the fairness measures (for those that have these parameters). Hyperparameter values were user-specified from the range of plausible values. For example, in the covariance constrained regression, $m$ can range from 0 to 1, and we considered $m\in\{0.2,0.4,0.6,0.8\}$. Comparisons of global fit versus group fairness for the three methods with variation in performance by hyperparameter can be found in Figure 1.

OLS had a cross-validated $R^{2}$ measure of 12.9%, a predictive ratio of 0.837 for individuals with MHSUD, and underestimated average MHSUD spending by -$1,872, with a mean residual difference of -$2,165. The fair covariance measure was 256. Average spending for enrollees without MHSUD was overestimated by $293 with a predictive ratio of 1.050. OLS had the worst performance along all fairness metrics while producing an $R^{2}$ only trivially higher than the competing methods.

We found the best improvement in fairness for MHSUD using the existing average constrained regression and our new covariance constrained regression. These two methods had similar, although not identical performance, and reduced the average undercompensation for enrollees with MHSUD to -$46 (versus -$1,872 in the OLS), a relative improvement of 98%. They also increased the predictive ratio from 0.837 to 0.996. Enrollees without MHSUD were overestimated by only $4 and had a predictive ratio of 1.001. Both methods reduced the fair covariance measure from 256 to 6. Unsurprisingly, these two estimators were also the worst performers on overall fit as measured by $R^{2}$, although it was a loss of only 4%, from 12.9% to 12.4%. This small 0.5 percentage point loss in $R^{2}$ may be tolerable to policymakers.

Recall that the weighted average constrained regression is a compromise estimator between the OLS and average constrained regression. As $\alpha$ approached one in the first panel of Figure 1, the metrics more closely resembled the OLS results. As $\alpha$ approached zero we saw values closer to the average constrained regression results, although $\alpha=0.2$ was not only dominated by the average constrained and covariance constrained regressions, but also the net compensation penalized regression with $\lambda=10000$.

The remaining two methods were regressions with customized penalty terms to punish unfair estimates. Our proposed net compensation penalized regression varied substantially by hyperparameter (see second panel in Figure 1), although was the third best performer overall when $\lambda=$10000. Large $\lambda$ values yielded extremely poor performance on both overall fit and fairness. At $\lambda=20000$, $R^{2}$ dropped by 12% to 11.9%, and when $\lambda$ increased to $30000$, $R^{2}$ dropped to 9%, a relative reduction of 29%. These two $\lambda$ values led to a large overcompensation for enrollees with MHSUD. The covariance was also negative, indicating that the residual value for MHSUD was systematically too high. The mean residual difference penalized regression was less sensitive to hyperparameters compared to the net compensation penalized regression (see third panel in Figure 1). The best performance for mean residual difference penalized regression was at $\lambda=30000$; it improved on the MHSUD predictive ratio for OLS by 7% (from 0.837 to 0.895) with an $R^{2}$ loss of less than 1%. However, the best performing net compensation penalized regression had an 81% improvement over the best performing mean residual difference penalized regression when comparing MHSUD net compensation, as well as large improvements in predictive ratios (0.895 versus 0.980) and fair covariance (164 versus 31).

We also examined the HCC variable coefficients for the best performing estimators, the average constrained and covariance constrained regressions, in comparison to OLS. Risk adjustment coefficients communicate incentives to insurers and providers related to prevention and care. For example, coefficients that do not reflect costs can impact an insurer’s incentives in creating their plan offerings. Coefficients for the average constrained and covariance constrained regressions were nearly identical when rounded to the nearest whole dollar, thus we display OLS versus covariance constrained regression in Figure 2. We considered the largest five increases and largest five decreases from OLS to covariance constrained regression, and observed sizable increases in the estimated coefficients associated with MHSUD. The largest relative increase was 180% for “Schizophrenia.” Relative decreases were much smaller.

A set of simulation scenarios was developed to explore how these regression methods perform in other settings. We generated a population of 100,000 observations with two continuous outcomes $Y_{1}$ and $Y_{2}$ that were each a function of covariates in $\boldsymbol{X}=(X_{1},X_{2},...,X_{9})$ and two distinct yet partially overlapping protected classes ($A_{1}$ and $A_{2}$) that depended on variables in $\boldsymbol{X}$. Scenario 1 considered a complex functional form for $Y_{1}$ and regression estimators that were misspecified, including omitted $\boldsymbol{X}$ variables. Scenario 2 examined a less complex functional form in $Y_{2}$ and regression estimators that were misspecified, including additional noise variables but no omitted $\boldsymbol{X}$ variables. A third scenario is discussed in Web Appendix B of the Supplementary Material, along with complete details for the simulated population and first two scenarios. For each scenario, we drew 500 samples of $N=$1,000 and $N=$10,000 observations from the simulated population of 100,000 observations. As in the data analysis, hyperparameter values were user-specified from the range of plausible values.

Selected results are presented in Figure 3, which includes OLS and those methods that improved fairness measures for protected class $A_{1}$ with a relative $R^{2}$ loss $\leq 10$%. Notably, average constrained and covariance constrained regression, the tied top estimators in our data analysis, do not appear. This was common across settings; average constrained and covariance constrained regression often struggled with functional form misspecification. However, net compensation penalized regression, which performed well in our data analysis, also performed well in the simulations with respect to achieving metric balance between global fit decreases and group fit increases. Additional results are available in Web Appendix B of the Supplementary Material.

We proposed new fair regression methods aiming to improve risk adjustment for undercompensated groups and asserted that a broader set of metrics is needed. As expected, there was no single method that performed the best across all the measures. One of our newly proposed techniques, net compensation penalized regression, had strong performance with respect to fairness and global fit in both the data analysis and simulations. Selecting the ‘best’ method relies on subjective decisions regarding how to balance group fairness versus overall fit tradeoffs. Improvements in fairness resulted in subsequent decreases in $R^{2}$. However, for many estimators, particularly in our data analysis, improvements in fairness were larger than the subsequent decreases in overall fit. This suggests that if we allow for a slight drop in overall fit, we could greatly increase compensation for MHSUD. Policymakers need to consider whether they are willing to sacrifice small reductions in global fit for large improvements in fairness.

We used a sample of enrollees in our demonstration. At scale in a policy implementation, data from millions of enrollees would be used to estimate health spending. Solutions to group undercompensation must be scalable, and current software may or may not yet be capable of handling the sample sizes required. We tested the CVXR optimization package on larger samples and found that it was able to find solutions on a sample of 1,000,000 observations over the span of 3 days (versus 7 hours for the 100,000 enrollee sample). While the optimization results were not within the ideal optimal threshold, they still converged and the results were similar to those presented in this paper, which is promising. Future work includes additional studies regarding scalability. In our analyses, we also considered a selection of user-specified hyperparameter values. A more thorough approach, with possibly improved results, would explore the hyperparameter space in an automated way to select values that optimize over joint fairness and fit objectives. As a general guideline, we found that $\lambda=N/10$ yielded reasonable metric balance for our newly proposed net compensation penalized regression.

We focused on one group that risk adjustment is known to disadvantage, but it is important to extend such strategies to multiple groups. Improvements for one group could result in subsequent undercompensation for other groups, and balancing fairness across an increasing number of groups is an as yet unsolved problem in risk adjustment. Our simulations examined two protected classes, and we found that improving fairness for one group did not generally help or harm the second group. Earlier research developing methods for the preprocessing phase found that reducing undercompensation for enrollees with MHSUD improved fairness measures for other groups, including enrollees with multiple chronic conditions but without MHSUD. Among the groups included in their comparisons, only enrollees with heart disease had slight reductions in fairness (Bergquist et al., 2019). But even the act of defining the groups poses a problem, as this can be subjective, potentially favoring larger groups with well-funded advocacy organizations. Undercompensation could be undetected in many other lesser-known groups. However, we can only measure undercompensation for groups that are identified by available data, and socioeconomic information, such as poverty and housing, are not available at the individual level for risk adjustment (Ellis et al., 2018).

Broadly, data-driven decisions have come under scrutiny for perpetuating human biases and disparities, which certainly exists in risk adjustment. Arguments for a more comprehensive view of research results is increasing among scientific researchers today (O’Neil, 2017, Gibney, 2018). Recent work argues that evaluating methods from a purely statistical standpoint can lead to negative consequences, and that policy aims should be better incorporated into our research (Corbett-Davies and Goel, 2018). Our article follows in this spirit, and we presented additional estimators and comparisons across multiple measures for the numerous (sometimes competing) goals of risk adjustment. While we worked within the specific context of risk adjustment, the fairness methods and measures discussed here have implications for other settings with continuous outcomes, which have been understudied relative to binary outcomes.

Web Appendix A presents the Group Residual Difference, a new extension of a fairness measure for comparing individual residual errors that is not currently practical to implement at scale in risk adjustment. Web Appendix B, referenced in Sections 2 and 4, describes the simulation study details and results. Web Appendix C is a list of HCCs included in the data analysis regression formulas from Section 3. The IBM MarketScan Research Databases used in Section 3 are not available for public dissemination as they contain protected patient information and we were granted access via a restricted data use agreement. Instead, we generated a simulated version of the analysis data that preserves important relationships while protecting the original content, as described in Web Appendix D. We provide this simulated analysis data with code as well as code to reproduce the simulation study from Section 4 with this submission and online: https://github.com/zinka88/Fair-Regression.

This work was supported by the National Institutes of Health through an NIH Director’s New Innovator Award DP2-MD012722. The authors thank the Health Policy Data Science Lab, Thomas G. McGuire, Berk Ustun, José Zubizarreta, and the anonymous reviewers for helpful comments.

Group Residual Difference. In the fairness literature, one definition for continuous outcomes is that persons with similar $Y$ should have similar predicted outcomes $\hat{Y}$ regardless of their protected class (Berk et al., 2017a). This relies on a user-defined distance function $d$, such as $|Y_{i}-Y_{j}|$, to ensure that people who are ‘close’ have similar outcomes. We extend this definition for risk adjustment by comparing residuals rather than predicted outcomes for the two groups: $[{1}/(n_{g}n_{c})\sum_{i\in g,j\in g^{c}}d(Y_{i},Y_{j})\{Y_{i}-\hat{Y}_{i}-(Y_{j}-\hat{Y}_{j})\}]^{2}.$ We refer to this new measure as the group residual difference. However, this measure is not practical to implement at scale in risk adjustment, which often involves millions of enrollees. The group residual difference requires comparing the residual of every enrollee in the undercompensated group to every other enrollee in the complement group. This scaling issue was also noted in the earlier work our metric extends upon (Berk et al., 2017a). Our group residual difference metric can be useful in settings where $N$ is smaller.

As described in Section 4 of the main text, our simulation study population of 100,000 observations considered covariates $\boldsymbol{X}=(X_{1},X_{2},...,X_{9})$, two protected class indicator variables ($A_{1}$ and $A_{2}$), and two continuous outcome variables ($Y_{1}$ and $Y_{2}$). $X_{1}$ was generated from a Normal distribution with mean 70 and standard deviation 15. Both $X_{2}$ and $X_{3}$ had Poisson distributions, with $\lambda$ values of 10 and 35, respectively. The last six covariates ($\boldsymbol{X}_{4:9}$) were drawn from Bernoulli distributions with probabilities 0.5, 0.1, 0.05, 0.8, 0.03, and 0.2. $A_{1}$ and $A_{2}$ were also drawn from Bernoulli distributions, but depended on other generated variables in the population:

	$\displaystyle\centering A_{1}\@add@centering$	$\displaystyle\sim$	$\displaystyle\text{Bernoulli}(X_{4}\times X_{9}/2+.01)$
	$\displaystyle A_{2}$	$\displaystyle\sim$	$\displaystyle\text{Bernoulli}(X_{4}^{2}/3+.05).$

They had prevalence rates of $6\%$ and $22\%$, respectively, with 2.1% overlap. Both outcomes, $Y_{1}$ and $Y_{2}$, depended on variables in $\boldsymbol{X}$ as well as $A_{1}$ and $A_{2}$:

	$\displaystyle\centering Y_{1}\@add@centering$	$\displaystyle=$	$\displaystyle(X_{1}\times X_{2}\times X_{4})+(A_{1}\times X_{2}\times X_{7})+(X_{3}\times X_{5}\times X_{6})+2^{(X_{8}\times X_{9})}$
			$\displaystyle+\phantom{1}(A_{1}\times X_{1}\times X_{5})+(A_{2}\times X_{3}\times X_{5})$
	$\displaystyle Y_{2}$	$\displaystyle=$	$\displaystyle X_{1}+X_{2}+(X_{3}\times X_{4}\times X_{5})+(A_{1}\times X_{3})+(A_{1}\times A_{2}\times X_{1}).$

We estimated regressions in three scenarios representing differing types of functional form misspecification:

	$\displaystyle\text{Scenario 1: }Y_{1}$	$\displaystyle=$	$\displaystyle\beta_{1}X_{1}+\beta_{2}X_{2}+\beta_{3}X_{3}+\beta_{4}X_{5}+\beta_{5}X_{6}+\beta_{6}X_{7}+\beta_{7}X_{8}+\beta_{8}X_{9}$
	$\displaystyle\text{Scenario 2: }Y_{2}$	$\displaystyle=$	$\displaystyle\gamma_{1}X_{1}+\gamma_{2}X_{2}+\gamma_{3}X_{3}+\gamma_{4}X_{4}+\gamma_{5}X_{5}+\gamma_{6}X_{6}+\gamma_{7}X_{7}+\gamma_{8}X_{8}+\gamma_{9}X_{9}$
	$\displaystyle\text{Scenario 3: }Y_{2}$	$\displaystyle=$	$\displaystyle\zeta_{1}X_{1}+\zeta_{2}X_{4}+\zeta_{3}X_{6}+\zeta_{4}X_{7}+\zeta_{5}X_{8}+\zeta_{6}X_{9}.$

Complete results for 500 draws from the population with $N=1,000$ and $N=10,000$ are given in Web Tables 1 and 2. Simulation data and complete analytic code to reproduce the simulation analyses are available online: https://github.com/zinka88/Fair-Regression.

HCC	Description
1	HIV/AIDS
2	Septicemia, Sepsis, Systemic Inflammatory Response Syndrome/Shock
6	Opportunistic Infections
8	Metastatic Cancer and Acute Leukemia
9	Lung and Other Severe Cancers
10	Lymphoma and Other Cancers
11	Colorectal, Bladder, and Other Cancers
12	Breast, Prostate, and Other Cancers and Tumors
17	Diabetes with Acute Complications
18	Diabetes with Chronic Complications
19	Diabetes without Complication
21	Protein-Calorie Malnutrition
22	Morbid Obesity
23	Other Significant Endocrine and Metabolic Disorders
27	End-Stage Liver Disease
28	Cirrhosis of Liver
29	Chronic Hepatitis
33	Intestinal Obstruction/Perforation
34	Chronic Pancreatitis
35	Inflammatory Bowel Disease
39	Bone/Joint/Muscle Infections/Necrosis
40	Rheumatoid Arthritis and Inflammatory Connective Tissue Disease
46	Severe Hematological Disorders
47	Disorders of Immunity
48	Coagulation Defects and Other Specified Hematological Disorders
54	Drug/Alcohol Psychosis
55	Drug/Alcohol Dependence
57	Schizophrenia
58	Major Depressive, Bipolar, and Paranoid Disorders
72	Spinal Cord Disorders/Injuries
75	Myasthenia Gravis/Myoneural Disorders, Inflammatory and Toxic Neuropathy
77	Multiple Sclerosis
78	Parkinson’s and Huntington’s Diseases
79	Seizure Disorders and Convulsions
80	Coma, Brain Compression/Anoxic Damage
84	Cardio-Respiratory Failure and Shock
85	Congestive Heart Failure
86	Acute Myocardial Infarction
87	Unstable Angina and Other Acute Ischemic Heart Disease
88	Angina Pectoris
96	Specified Heart Arrhythmias
99	Cerebral Hemorrhage
100	Ischemic or Unspecified Stroke
103	Hemiplegia/Hemiparesis
107	Vascular Disease with Complications
108	Vascular Disease
111	Chronic Obstructive Pulmonary Disease
112	Fibrosis of Lung and Other Chronic Lung Disorders
114	Aspiration and Specified Bacterial Pneumonias
122	Proliferative Diabetic Retinopathy and Vitreous Hemorrhage
134	Dialysis Status
135	Acute Renal Failure
136	Chronic Kidney Disease, Stage 5
137	Chronic Kidney Disease, Severe (Stage 4)
161	Chronic Ulcer of Skin, Except Pressure
167	Major Head Injury
169	Vertebral Fractures without Spinal Cord Injury
170	Hip Fracture/Dislocation
173	Traumatic Amputations and Complications
176	Complications of Specified Implanted Device or Graft
186	Major Organ Transplant or Replacement Status
188	Artificial Openings for Feeding or Elimination

The IBM MarketScan Research Databases analyzed in Section 3 of the manuscript cannot be distributed online due to their proprietary nature. They also contain protected patient information. Thus, we created a simulated data set with similar properties using key features and relationships from the original data for reproducibility analyses of our code. The simulated analysis data described below and accompanying code to complete the analyses are available online: https://github.com/zinka88/Fair-Regression.

First, we simulated demographic variables, female and age, by sampling from a Bernoulli distribution, $\text{female}\sim\text{Bernoulli}(0.52)$, and truncated Normal distribution with lower bound $a$ and upper bound $b$: $\text{age}\sim\text{Normal}(44,12,a=21,b=63).$ Next, we generated the 62 binary health variables $\boldsymbol{H}=(H_{1},...,H_{T})$, each drawn from a Bernoulli distribution and dependent on the demographic variables female and age with coefficients determined by the relationships in the original data. To create the indicator for MHSUD, $A$, we generated 15 binary MHSUD CCS variables $\boldsymbol{C}=(C_{1},...,C_{15})$ dependent on age, female, and the top six HCCs correlated with MHSUD in the original data. We defined $A=1$ for all observations with at least one MHSUD CCS; 15.7% of observations in the simulated analysis data had MHSUD compared to 13.8% in the original data.

The results from the simulated analysis data are shown in Web Table 3. As demonstrated in our data analysis presented in Section 3 of the main text, we likewise find that the constrained and penalized estimation methods improve fairness measures without a significant decrease in $R^{2}$. The relative rankings of all the methods are similar, although we highlight that net compensation penalized regression performs even more similarly to average constrained and covariance constrained regression methods here.