Policy Learning with Observational Data: The Case of Hepatitis C Treatment for HIV/HCV Co-Infected Patients

Policy Learning with Observational Data : The Case of Hepatitis C Treatment for HIV/HCV Co-Infected Patients

\thispagestyle{empty}

Introduction

\setstretch{1.3} \pagenumbering{arabic} When choosing from a menu containing a finite number of alternatives, economic agents do not necessarily have all the information they would like. Entities that have more, or at least different, information may establish rules that aim to guide decisions, often in the hope of increasing welfare or reducing costs. For instance, in clinical settings, physicians make decisions about the treatment they will provide to a given patient in the context of uncertainty regarding the expected benefits or side effects.\footnote{For simplicity, I refer to physicians as the “decision-makers" throughout the paper, although treatment choice is often the product of a joint decision between the patient and the physician.} Third-party health organizations (referred to as the “policymakers" hereafter) produce clinical practice guidelines (CPGs) to guide and inform physicians’ treatment decisions. For example, the U.S. Centers for Disease Control and Prevention (CDC) produces periodic recommendations concerning the prevention and control of seasonal influenza with vaccines grohskopf_prevention_2024, whereas the American College of Gastroenterology (ACG) produced in 2020 a clinical guideline on the management of irritable bowel syndrome lacy_acg_2021. In Canada, the Canadian Association for the Study of the Liver (CASL) updated its guidelines on the management of chronic Hepatitis C (HepC) in 2018 by advocating “against the use of any interferon-based treatment regimens" due to the availability of direct-acting antivirals for nearly all patients shah_management_2018. CPGs were initially defined in 1990 by the Institute of Medicine as “systematically developed statements to assist practitioner and patient decisions about appropriate health care for specific clinical circumstances" field_clinical_1990. Nowadays, CPGs are mostly designed by multidisciplinary development groups through consultations with experts and patients, analysis of other CPGs, and extensive literature reviews kredo_guide_2016, de_leo_approaches_2023. Even though several attempts to standardize the development of CPGs have been performed recently, there is no standard approach in the development of CPGs, which can affect the quality of the published recommendations kredo_guide_2016. In the clinical context, several challenges make creating high-quality CPGs difficult and can limit their widespread acceptance and use. First, while the policymaker’s objective function may align with the decision-maker’s, principal-agent settings may also arise where the policymaker (principal) uses guidelines to align the decision-maker’s (agent’s) behavior with their own objectives. The extent to which the objective sought by the policymaker (e.g., increase expected welfare, minimize the risk of serious adverse events, cost control, etc.) is acceptable to the decision-maker can therefore play a role in whether and how guidelines are used. Throughout the paper, I assume that the agent and the principal both share the same objective of maximizing expected welfare up to certain constraints defined by the policymaker, where guidelines are one of the many tools policymakers can use to influence the decision-makers' choices. Although strict compliance with guideline recommendations could be enforced in practice, uncertainty in health outcomes and information asymmetry is likely to make such a regulation inefficient in several cases boone_discretion_2025,einav_limitations_2024. Second are a set of concerns related to the quality of the underlying clinical evidence used to formulate CPGs graham_clinical_2011. Evidence coming from meta-analysis and randomized controlled trials (RCTs) is often considered the most internally valid among all types of evidence chloros_has_2023. However, it is well known that RCTs often lack external validity : results are often not generalizable to important subsets of the target population due to restrictive eligibility criteria when conducting the trials lawlor_commentary_2004, atkins_creating_2007. When the characteristics of the trial participants in studies that inform CPGs do not match the target population, the applicability of the findings and treatment recommendations is questioned. While studies based on observational data may have the benefit of increased representativeness, concerns about biased treatment effect estimates are prominent. Although quasi-experimental designs that can address confounding are generally not well understood or accepted in clinical contexts, interest in the use of real-world data by industry, health care payers, and regulatory bodies has been increasing in recent years cadth_guidance_2023, inesss_state_2022, us_food_and_drug_administration_use_2023. The policy learning (PL) literature has also considered policymakers using observational (offline) data to create guidelines (policy rules) to inform decision-makers’ choices from among a set of treatment options. The use of observational data with rigorous research designs has also been advocated in the clinical literature in cases where typical RCTs cannot answer the research question, costs are prohibitive, or when it is impossible to randomize interventions due to legal or ethical concerns atkins_creating_2007,greenfield_making_2017,reeves_including_2019. Although results from observational studies are often biased, several methods have been and are still being developed to reduce biases and better inform policymakers on the effectiveness of medical interventions in the real world wager_estimation_2018,basu_2sls_2018,chernozhukov_doubledebiased_2018. In this paper, I show how to derive a feasible policy rule from non-randomized, observational data. Specifically, I consider the choice of pharmaceutical treatment options for patients co-infected with HepC and HIV, a setting where at least 3 treatment options exist, and there are no clear CPGs in the Canadian context. I use weaker and more practical assumptions than the ones typically invoked in the corresponding PL literature, namely that treatment assignments are unconfounded conditional on the value of a latent state, which is used to form groups of homogeneous patients. I then estimate conditional average treatment effects (CATEs) for each group and use these estimates to derive the corresponding policy rule, where the conditionality is with respect to the group membership. A feasible version of this rule uses decision trees to predict group membership with a small subset of covariates, which could then be leveraged to improve existing CPGs or build new CPGs. The first contribution of this paper is the development of a more practical and parsimonious approach than the ones typically employed for policy learning with observational data. This approach is described below and contrasts with the traditional PL literature, which usually assumes that treatment assignment is independent of potential outcomes for every treatment choice conditional on observed characteristics, an assumption known as unconfoundedness (selection on observables; kitagawa_who_2018,kitagawa_equality-minded_2021, chernozhukov_semi-parametric_2019,mbakop_model_2021,athey_policy_2021, zhou_offline_2023). However, unconfoundedness is often unrealistic in real-world problems where important confounders may remain unobserved due to practical or ethical considerations kallus_confounding-robust_2019. More importantly, unconfoundedness does not tell applied researchers which variables should be conditioned on (and at what values) to overcome selection biases. The choice of conditioning variable can be a serious dilemma even in low-dimensional settings where the number of observed covariates is less than the number of observations, but still relatively large chernozhukov_fisher-schultz_2025. If the use of machine learning (ML) techniques somewhat alleviates this concern by selecting subsets of $X$ that best predict the observed outcome, it remains unclear which technique offers the best performance for a given context.\footnote{Specifically, athey_policy_2021 argue that their approach, which relies on unconfoundedness, “allows for the use of black-box machine learning tools provided we can verify that they are accurate in mean-squared error". When moving from a binary to a multi-action PL problem, zhou_offline_2023 state that “[t]his leaves open a question of both theoretical and practical significance, namely which algorithm should be used [...] to achieve maximum statistical efficiency."} To address the latter problem, I replace unconfoundedness by a different, more practical assumption that I call “group-wise" unconfoundedness : assignment to treatment is unconfounded conditional on the value of a latent state, unobserved by both the econometrician and the decision-maker, rather than on the value of all observed covariates. In the context of a physician prescribing from among a choice of drugs or medical procedures, it is often reasonable to assume that treatment response heterogeneity will depend on a discrete and finite number of health states rather than on each possible value taken by the entire set of covariates. The health state can be a function of different observed (e.g., age, chronic conditions, income) and unobserved (e.g., frailty, social support) characteristics. Consequently, group-wise unconfoundedness is more parsimonious than standard unconfoundedness since it implies that the conditioning set can be reduced from the space of all possible values taken by the covariates to a discrete state space without introducing any selection bias. This leads to CATE estimates that are less prone to overfitting and that are equally or more precise than the ones produced by standard ML techniques such as random forests, provided that the true number of groups is small relative to the sample size. The estimated CATEs will also be robust to the presence of unobserved confounders, as long as those confounders are constant within each group. Under group-wise unconfoundedness, consistent estimation of all CATEs will be possible if the true group memberships are known for all observations in the sample. It is therefore necessary to rely on a consistent estimator of group memberships to make the entire estimation procedure globally consistent. Clustering algorithms such as the K-means algorithm are a practical and computationally efficient way to estimate such group memberships. The second contribution of this paper relates directly to the literature on clustering algorithms; I show how a weighted K-means algorithm can be used to consistently estimate all within-group parameters (including CATEs) under a set of relatively weak assumptions formulated at the group level. I then show how to use those estimated CATEs to build an optimal policy rule that lies within the class of all policy rules that are deemed feasible by the policymaker.\footnote{Since all policy rules are by definition optimal for a given class of feasible policy rules, the adjective “optimal" becomes superfluous. This is why I omit it for the rest of the paper unless relevant.} One advantage of using the proposed weighted K-means algorithm is that it uses all available covariates to recover a small number of groups where the outcome of each untreated individual is used as a counterfactual for every treated individual in the same group. Another advantage is that all estimated CATEs remain constant within each group if all outcome models are separable in $X$ (i.e., every element in $X$ affects the outcome only marginally). In this case, predicting every CATE for a given patient reduces to predicting his individual group membership. This contrasts with other ML techniques where selecting a subset of $X$ may discard useful information, -- hence leading to regularization bias chernozhukov_doubledebiased_2018 -- or where the large number of estimated CATEs makes the interpretation of the results challenging. Although the weighted K-means algorithm leads to interpretable results within each group, estimation of group memberships relies on a (potentially very) large number of covariates. Since group membership estimation is essential for correctly predicting every CATE for a given patient, this can lead to practical issues in clinical settings where decision-makers may not be able to correctly generate or use all the information required by the algorithm. In order to translate this policy rule into a feasible one that is implementable by decision-makers, I consider the policy class of classification/decision trees, which are known for their flexibility and interpretability, and are also widely used in healthcare to guide decision-making processes podgorelec_decision_2002. More precisely, I use classification trees to predict the estimated group memberships based on a clinically relevant subset of $X$ that is observable by both the econometrician and the decision-maker. The resulting decision tree might then be used to complement or replace existing CPGs once its predictive performance has been empirically validated wei_zhao_construction_2020. However, the use of classification trees for clinical decision-making entails a trade-off between feasibility and prediction accuracy : the simpler the tree, the greater the classification error and vice versa. If both the sampled and the target populations have a similar treatment response distribution, then the policymaker could implement a policy rule by allocating treatment to anyone in the target population based on their predicted group membership. In the context of the weighted K-means algorithm, a similar treatment response distribution between the sampled and target populations means that both populations contain the same set of latent groups, with sufficiently large sample sizes such that every treatment option is observed within each group. Therefore, differences between the sampled and the target population do not lead to a lack of generalizability as long as the two populations are composed of the same groups with sufficient group sizes in the sampled population. This constitutes another benefit of the suggested approach compared to ML techniques, where significant deviations from perfect representativeness may lead to poor external validity bay_machine_2024. If the population of interest consists of a single group, the corresponding policy rule would be to assign the welfare-maximizing treatment option to everyone in this population (provided that the policymaker has access to a consistent estimate of the average treatment effect for every option). If such a situation is verified, then it is unnecessary to rely on the weighted K-means algorithm to guide clinical decision-making. However, it is often highly unrealistic to assume treatment effect homogeneity for all treatment options, thus advocating for the development of policy rules under assumptions that can adapt to a wide array of distinct patterns in terms of treatment effect heterogeneity. In the medical literature, such policy rules are often characterized as “personalized medicine" and are increasingly used to improve both quality of care and patients' quality of life chakraborty_statistical_2013,einav_limitations_2024. So far, most of the literature on PL has focused on the asymptotic and statistical properties of different estimators/algorithms whose objectives are to learn a policy rule $\pi^{*}$ within a predetermined class of feasible policies $\Pi$ under unconfoundedness or bounded violations of unconfoundedness kitagawa_who_2018, chernozhukov_semi-parametric_2019, athey_policy_2021, kallus_minimax-optimal_2021, mbakop_model_2021, zhou_offline_2023,ek_off-policy_2023, kitagawa_constrained_2023, sahoo_learning_2024, viviano_policy_2024,zhan_policy_2024,hess_efficient_2025. To the best of my knowledge, these estimators/algorithms have not yet been used in any applied work except for those presented in the papers cited above. Therefore, this paper also serves to bridge the gap between theory and practice via a concrete example of how a policy rule can be derived using standard classification/regression procedures with observational data from a large-scale prospective cohort. The paper's second contribution relates to the application of the proposed methodology. According to saeed_frequent_2024, 20% to 30% of Canadians infected with the human immunodeficiency virus (HIV) are also infected with the hepatitis C virus (HCV). This rate goes above 50% among people who inject drugs (PWID). According to the CDC, about 21% of people with HIV in the United States are also infected with HCV, while this rate lies between 62 and 80% for PWID with HIV cdc_viral_2025. Even though the evolution of HIV/HCV co-infection is similar in the two countries, the U.S. guidelines recommend prolonged treatments for co-infected individuals (compared to HCV mono-infected ones), whereas Canadian guidelines make no distinction between treatments for HCV mono-infected and HIV/HCV co-infected individuals martel-laferriere_prise_2022-1. Moreover, potential drug-drug interactions between antiretroviral therapy (ART) and the recently approved direct-acting antiviral (DAA) agents used to treat HCV create an additional unknown risk that needs to be assessed before formulating recommendations or developing new guidelines american_association_for_the_study_of_liver_diseases_patients_2022. Even though DAAs are very effective in the HCV mono-infected population, it is less clear in HIV/HCV co-infected patients since most co-infected individuals were not eligible to participate in DAA trials saeed_how_2016. Observational studies looking at the effectiveness of DAAs administered to co-infected persons obtain significantly lower \textit{sustained virological response} (SVR) rates compared to HCV mono-infected persons alam_real-world_2017, wilton_real-world_2020. However, it is believed that such lower responses to treatment are due to social and behavioral rather than biological factors, especially among PWID wilton_real-world_2020. Reaching SVR is equivalent to being cured from the disease and implies an undetectable viral load in the blood 12 weeks after the end of the treatment. Note also that, unlike Hepatitis B, there is no vaccine against HCV, which makes reinfection always possible for any previously cured patient. Taking on the policymaker's perspective, I use data from the Canadian Co-Infection Cohort (CCC) to learn a policy rule that can be implemented for routine decision-making processes. The objective of this policy rule is to maximize the probability of reaching SVR for all observations in the sample under a set of prespecified constraints. Four treatment options are considered, including the absence of treatment (the “null" treatment option). Decision trees are then used to derive the feasible version of the policy rule. Such a feasible rule could then be included in CPGs to guide future decision-making processes by prescribing physicians, provided that the sampled and the target populations share a similar treatment response distribution. This is believed to be the case with the CCC data given the cohort's very broad eligibility criteria, and given that the CCC has enrolled around 23% of all HIV/HCV co-infected patients engaged in care in Canada saeed_frequent_2024. Assuming that the objective of the policymaker is to maximize social welfare, an “optimal" policy rule from the perspective of the policymaker balances the cost of each treatment option with its expected health benefits. As suggested by kitagawa_who_2018, it is possible to do this by reformulating the outcome $Y_i$ such that it corresponds to the health benefits net of the costs for each treatment option. However, such an approach implies that the net outcome becomes a function of the policymaker's willingness-to-pay (WTP). A more realistic and practical approach is to adjust the unconstrained policy rule -- that is, a policy rule that does not consider any cost or budget constraint -- so that the recommended treatment option maximizes welfare while satisfying the policymaker's WTP. These constraints can be incorporated into the infeasible policy rule based on all available information, or into its feasible version based on the tree's predicted group memberships. Another particularity of the empirical application is that patients do not necessarily fully adhere to treatment over time. Given that adherence to treatment can be endogenous with respect to potential outcomes, this requires group-wise unconfoundedness to hold for any level of adherence to treatment.\footnote{Adherence to treatment corresponds to the proportion of days covered (PDC, see, for instance, prieto-merino_estimating_2021 for more details) and is defined as the ratio between the time during which the patient has taken the treatment and the prescribed duration of the treatment. Therefore, low adherence can be caused by either a lack of patient compliance with the prescribed regimen or a switch to another treatment regimen. The idea is that adherence should be directly correlated with the degree of exposure to the prescribed treatment.} Moreover, imperfect adherence to treatment adds another step to the learning problem since adherence is likely to vary across treatment options for a given patient. For example, it might be preferable to prescribe a treatment option with a shorter duration to individuals who are less likely to fully adhere to treatment, even if the prescribed treatment is less effective than all other options under perfect adherence. This is why predicted adherence is incorporated as another “layer" of the learning policy problem under the hypothesis that the probability of reaching SVR never decreases with higher adherence. To fix ideas, Figure (ref) draws a schematic representation of the different objects used in the empirical application to solve the problem of interest. First, I estimate all group memberships via the weighted K-means algorithm. Estimated group memberships are then used to identify within-group CATEs using standard linear probability models with observed SVR as the outcome. The selection of hyperparameter values and variables to include in the outcome models is assessed through 5-fold cross-validation with several different data splittings. Inference on the within-group CATEs is performed during the cross-validation procedure following the general logic of chernozhukov_fisher-schultz_2025. The estimated CATEs are then used to derive both infeasible and feasible policy rules. Finally, adherence to treatment is modeled via a two-part model where predicted adherence is then incorporated into both policy rules to account for variations in adherence across individuals and treatment options. All four combinations of policy rules are represented graphically in Section (ref) as a function of the policymaker's WTP. The resulting policy rules can also be interpreted as four different cost-effectiveness analyses where the in-sample optimal treatment allocation is derived for various levels of the policymaker's WTP, conditional on the cost of each treatment option.

Estimation results show that out-of-sample (OOS) prediction errors are minimized when the total number of groups in the empirical application is set to five. This result is achieved using simple linear probability models within each group, where the outcome corresponds to whether a patient has reached SVR or not ten weeks after the end of the prescribed treatment. All estimated CATEs are strictly positive, but important treatment effect heterogeneity is identified across groups. More precisely, one of the groups identifies CATEs that are not significantly different from zero for two out of four treatment options, and where the probability of being “spontaneously" cured from the HepC infection is very high ($\approxeq 80\%$).\footnote{Spontaneous clearance is associated with an acute infection of HepC, and is characterized by an absence of detectable viral load in the patient's blood up to 6 months after the onset of the infection. On the other hand, patients with a chronic HepC infection do not usually experience spontaneous clearance within 6 months. More details on spontaneous clearance are presented in Section (ref)} According to the simplest decision tree that distinguishes between all five groups, belonging to this group can be predicted with pre-treatment characteristics such as the degree of liver cirrhosis/fibrosis, the blood concentration of certain enzymes or hormones, and whether the patient is indigenous or not. Simple descriptive statistics show that the average patient within this group is slightly younger and has a higher income compared to the average patients in all other groups. The four other groups feature median CATEs that are all significantly different from zero (ranging from 0.7 to 0.94), with similar orderings between the effectiveness of all treatment options. In addition to the literature on PL that was cited earlier, this paper also relates to the literature on statistical decision theory in econometrics hirano_asymptotics_2009, hirano_asymptotic_2020, stoye_minimax_2009,stoye_minimax_2012, bhattacharya_inferring_2012, tetenov_statistical_2012, manski_statistical_2004, manski_econometrics_2021, to the literature on clinical practices with HIV/HCV co-infected patients laniece_delaunay_gaps_2022,peters_retrospective-prospective_2021,roman_lopez_effectiveness_2024,saeed_eliminating_2020,saeed_how_2016,wilton_real-world_2020, and also to the literature on individualized treatment rules in biostatistics and epidemiology chakraborty_statistical_2013, wallace_doubly-robust_2015, zhao_estimating_2012, zhao_new_2015, zhang_robust_2013, li_optimal_2023. The first strand of literature focuses on the statistical properties of different decision criteria when the treatment effect depends on an unknown “state of nature". When applied to the problem of interest, the possible states of nature are analogous to the group memberships recovered by the K-means algorithm. More details on this topic are given in Section (ref). This paper contributes to the second strand of literature by comparing the effectiveness of three DAA regimens used with HIV/HCV co-infected patients : Mavyret, Epclusa, and Harvoni. Our results show that Harvoni weakly dominates Epclusa in terms of median treatment effect for all five identified groups. In other words, it is never strictly optimal to prescribe Epclusa to any given patient under infinite maximum WTP and perfect adherence to treatment. However, Epclusa becomes the optimal treatment choice under finite WTP and perfect adherence since it is less expensive but as effective as Harvoni for all patients in one of the five identified groups.\footnote{Epclusa is also the treatment that is minimax-regret optimal under perfect adherence; see Section (ref) for more details.} To summarize, my results show that the optimal treatment choice among the considered options is not trivial and depends on a set of pre-treatment covariates, the policymaker's maximum WTP, and the chosen decision criterion. Finally, the last strand of literature is mainly concerned with dynamic treatment regimes where treatment decisions need to be made at different points in time in the future. Such a setup goes beyond the scope of this paper, but several insights derived from this literature are used in this paper. The remainder of the paper is structured as follows. Section (ref) introduces the notation and the general framework of the problem of interest. Section (ref) describes the assumptions and the proposed methodology used to solve this problem. Section (ref) presents the characteristics of the cohort and the associated health conditions (HCV and HIV). Section (ref) also details how various policy rules can be derived from this dataset when cost constraints are applied to the class of feasible policies. Section (ref) presents the results, while Section (ref) briefly discusses the results and concludes.

Notation and General Framework

Notation

Suppose we have access to a random sample of $N$ independent individuals, with each individual represented by the triple $(Y_i,X_i,A_i)$, where $Y_i \in \mathcal{Y} \subset \mathbb{R}$ is the $i^{th}$ individual's observed outcome, $X_i = (X_{i,1},...,X_{i,p})^{\top} \in \mathcal{X} \subseteq \mathbb{R}^{p} $ is a vector of pre-treatment covariates of size $p$, and where $A_i = (A_{i,0},...,A_{i,D}) \in [0,1]^{D+1}$ corresponds to the set of observed adherence to treatment. The observationally assigned treatment for individual $i$ is denoted by $D_i \in \mathcal{D} = \{0,1...,D\}$. Absence of treatment is represented by the number “0", which automatically leads to $A_{i,0}=0$ if no treatment was ever given to individual $i$. For now, I consider the setup where $A_{i,d} = 0$ for all $d \in \mathcal{D}\backslash D_i$, which implies that I exclude cases where patients can switch treatments over time and cases where patients have been treated in the past. This restriction can be easily relaxed provided that information about past treatment history is included in $X_i$ and that the subscript $i$ refers to individual-period rather than individual alone. Cases where the adherence is always either null or perfect can be easily obtained by imposing the constraint $\sum_{d=0}^D A_{i,d} = \{0,1\}$ to this general setup. Given that $A_{i,d} \in [0,1]$ for each value of $d \in \mathcal{D}$, I define the set of potential outcomes as $\{Y_i(d,a)\}_{d=0}^D$, where $Y_i(d,a)$ corresponds to the value we would have observed had the adherence $A_{i,d}$ for treatment $d$ been equal to $a$. Consequently, the following equivalence holds only for the observed outcome : $Y_i = Y_i(D_i, A_{i,D_i})$. All other potential outcomes have to be estimated from the data for any value of adherence to treatment between 0 and 1. It is also possible to modify the notation so that potential outcomes depend also on the individual's current “health state"; this modified notation is introduced in the next section.

Policy Rule and Decision Criteria

A policy rule $\pi(X_i)$ is a mapping from the space of covariates $\mathcal{X}$ to a decision $d \in \mathcal{D}$, where $\pi(X_i) = d$ implies that the policy rule $\pi$, when applied to a patient with covariates $X_i$, selects the $d^{th}$ treatment option as the “optimal" option. Different notions of optimality are explicated below. For simplicity, this setup excludes probabilistic assignments where the resulting policy rule would randomize treatment assignments for a given value of $X_i$. If we assume a utilitarian social welfare function and that utility is a known and strictly increasing function of $Y_i$, then the in-sample policy rule is defined as follows

where adherence is assumed, for now, to be constant across treatment options and individuals. A more general version of this problem treats adherence as a function of observed variables (see Section (ref)). Note that the maximization is performed over the sum of individual outcomes, which is important if $\Pi$ is globally constrained, and also that $Y_i$ does not incorporate any notion of costs, thus implying that the above policy rule is “optimal" from the patients' perspective.\footnote{It is also optimal from the physicians' perspective if the financial incentives are aligned with each patient's welfare, which is implicit throughout the paper.} Note also that since $Y_i$ is a scalar, there is no need to account for the patient's preferences in the above optimization. This logic remains valid if higher values of $Y_i$ have an impact on other outcomes that are beneficial, on average, to the patient's welfare (e.g., higher productivity, social effects, etc.). Following the literature on SDT, I now define the state variable $s_i \in \mathcal{S} = \{1,2...,S\}$ as an unobserved, discrete, and finite random variable that governs the distribution of potential outcomes for individual $i$ conditional on both treatment assignment and adherence. The state variable $s_i$ represents the unobserved “true" state value for individual $i$. Such a variable allows the introduction of fundamental uncertainty in the relationship between the treatment option $d$ (and the related adherence) and the outcome $Y_i$.\footnote{In the literature on SDT, the state value is assumed to be constant across individuals, hence representing the “state of nature" and how it affects potential outcomes. I depart from this literature by treating the state variable as a latent state variable that characterizes the individual's (unobserved) heterogeneity. As a result, the state variable is similar to a sampling indicator, where each possible value taken by $s_i$ corresponds to a different sampling distribution, as described in Section 2.2 of manski_econometrics_2021. This also explains why individuals are not assumed to be identically distributed.} If the value of $s_i$ is known for all $i \in [N]$, the in-sample policy rule can be written as follows

The subscript $C$ refers to the fact that this simple setup assumes no fundamental uncertainty about the relationship between the preferred treatment option $\pi(X_i)$ and the outcome $Y_i$. In practice, the value of $s_i$ is usually unobserved. In this case, we can take on a Bayesian perspective and write the in-sample policy rule as follows

where $\mathbb{E}_S$ is the expected value with respect to $s_i$, and $\Pr[s_i=g]$ is a subjective prior distribution over $s_i$. The notation $Y_i(d,a,g)$ makes it clear that potential outcomes are now treated as a function of the chosen treatment, the degree of adherence to this treatment, and the value of the state variable. Note that this fundamental uncertainty adds to the idiosyncratic uncertainty and needs to be properly taken into account for inference. Since $s_i$ is unobserved for all $i \in [N]$, it might be convenient for the policymaker to build a policy rule that yields good performance over all possible states for every individual. This leads to two additional decision criteria : the maximin criterion and the minimax regret criterion. The former criterion maximizes the minimum welfare attainable across all states, whereas the latter minimizes the maximum welfare loss incurred by the chosen policy across all states. The in-sample versions of the policy rules that can be derived from those criteria are defined as follows, respectively

which denotes the regret of applying policy $\pi$ to individual $i$ in state $g$. This corresponds to the difference between the maximum welfare attainable in state $g$ for individual $i$ and the welfare generated by policy $\pi$ in the same state for the same individual. The policy rules $\pi^{*}_{M}$ and $\pi^{*}_{MR}$ will generally differ from each other, except in special cases where the maximum welfare is invariant across all states manski_econometrics_2021. Theoretically speaking, choosing a specific criterion over another is still an open question. The statistical properties of the minimax-regret criterion have been widely studied in econometrics since it is argued to yield less pessimistic results than the maximin criterion, especially when applied to results from randomized experiments manski_adaptive_2007, manski_econometrics_2021,stoye_minimax_2009,stoye_minimax_2012,tetenov_statistical_2012, hu_minimax-regret_2024. On the other hand, the Bayesian criterion is easy to interpret and implement, provided that one has access to a prior distribution over the different states. However, it is well known that different priors can lead to different posterior estimates, hence different policy rules. Other decision criteria and policy rules that account for estimation uncertainty have also been recently proposed in the literature chernozhukov_policy_2025. However, such policy rules depend on the degree of risk aversion of the decision-maker, and are more useful in cases where larger estimated marginal treatment effects are associated with larger variances and vice-versa. This situation does not correspond to what is observed in the empirical application. This is why this paper focuses on the “Certainty" and the Bayesian decision criteria.

Treatment Effect Heterogeneity

A central feature of PL is treatment effect heterogeneity. Suppose each treatment option features a single, homogeneous treatment effect across all states possibly present in the population of interest. Then the corresponding policy rule would be to assign the treatment option with the highest true treatment effect to all individuals in the population. Introducing treatment effect heterogeneity in the population of interest is not only more realistic, but it also allows the policymaker to design more or less complex policy rules that will ultimately lead to higher welfare in the population for a given cost level. Following the notation developed in the previous section, I generally define the true CATE as follows

where the expected value is taken over $i$ conditional on $X_i=x$ for a given value of adherence $a$, of treatment option $d$ and of state/group variable $g$. It is also possible to obtain the definition of the ATE from Eq. ((ref)) by integrating over the distribution of $X_i$ as follows

where $p_{X}(x)$ is the (multivariate) density of the covariates when $X$ is treated as random. I also define the conditional average treatment effect on the treated (CATT) as follows

where the expectation is now conditional on both $X_i$ and the observationally assigned treatment option $D_i$. This quantity will differ from $\tau_{d,g}(a,x)$ if $D_i$ is correlated with $\mathbb{E}[Y_i(d,a,g)]$ when treated as a function of treatment choice $d$ and adherence $a$. Average treatment effect on the treated (ATT) can be obtained by substituting $\tau_{d,g}(a,x)$ by $\eta_{d,g}(a,x)$ in Eq. ((ref)). One of the challenges in estimating either CATEs or CATTs is that it requires a very large sample size to obtain precise estimates for each relevant partition of $X_i$, especially when some of the covariates in $\mathcal{X}$ are continuous. Obtaining precise values for each estimated CATE generally implies that a large number of individuals share identical or similar values of covariates for each treatment option in $\mathcal{D}$, which might not be feasible even in large samples. One way to deal with this issue is to use non-parametric, kernel-based methods or machine learning techniques such as random forests and neural networks so that observations with sufficiently close values of $X_i$ will be used to produce an estimate for each CATE. Another way to deal with this issue is to assume that the state variable $s_i$ captures treatment effect heterogeneity, and that the observed covariates are an unknown function of the state's value for every individual in the target population. This assumption is formalized in Assumption (ref).

Assumption (ref) is considered easy to satisfy in practice since a sufficiently large number of groups $\mathcal{S}$ will incorporate most of the true unobserved heterogeneity present in the target population. Furthermore, assuming that $\Pr[s_i = g|X_i = x] \ne \Pr[s_i = g]$ implies that the state variable and the covariates are not independent of each other, which is generally less restrictive than the opposite. However, this assumption might be difficult to verify in nonlinear models where marginal effects are functions of the covariates. A nonlinear version of Assumption (ref) could be formulated at the level of the linear predictor when using a nonlinear parametric specification. The former assumption does not imply that covariates have no role in estimating potential outcomes since some elements in $X_i$ might be important predictive variables of the outcome. Predictive variables correspond to variables that affect the outcome of interest only marginally and feature no qualitative interaction with any treatment option zhang_variable_2018,zhang_subgroup_2022. Omitted-variable biases are likely to occur if predictive variables correlated with the assignment mechanism (or the adherence level) are not properly “controlled for" when estimating CATEs from predicted outcome values. This is not an issue in randomized controlled trials (RCTs) where treatment assignment is randomized and compliance with the assigned treatment is (almost) perfect. This is however an issue with observational data where treatment assignment is not randomized and adherence may be far from perfect. Nonetheless, not properly controlling for correlated predictive variables (i.e., confounders) is a second-order issue in the context of PL problems if the resulting biases do not alter the ordering of all CATEs within each group. For instance, it is often convenient to assume a partially linear relationship of the form $Y_i(d,a_i,g) = a_i\tau_{d,g} + h_{d,g}(X_i) + \epsilon_{i}$ robinson_root-n-consistent_1988,chernozhukov_doubledebiased_2018, where $h_{d,g}(\cdot)$ is an unknown function of $X_i$ for treatment option $d$ in group $g$, and where the treatment effect $\tau_{d,g}$ is assumed to be constant for any $a_i \in [0,1]$. In this setup, $X_i$ is a set of predictive variables that has to be included in the model to obtain unbiased predictions of the outcome $Y_i(d,a_i,g)$, but not of the true CATEs $\tau_{d,g}$ if $a_i$ is uncorrelated with $X_i$. Conversely, prescriptive variables exhibit a qualitative interaction with the prescribed treatment and are required to consistently recover the distribution of the true CATEs. Prescriptive variables would therefore correspond to the set of variables in $X_i$ that can be used to predict the true group membership $s_i$ for all $i \in [N]$, as shown in the next section.

Solving the PL problem

Consistent Estimation of the CATEs

It is well known that selection into treatment might bias treatment effect estimates when treatment assignments are not randomized. Even if treatment assignments are randomized, imperfect adherence to treatment can also lead to biased estimates if adherence is not independent of potential outcomes conditionally on the observed covariates bowden_connecting_2021, bowden_instrumental_2024. When using real-world data with potentially imperfect adherence to treatment, it becomes useful to assume that both treatment assignment and adherence are independent of potential outcomes, conditional on observed covariates for any degree of adherence. This assumption, known as unconfoundedness, is formalized below in Assumption (ref).

Assumption (ref) differs from the typical unconfoundedness assumption since it requires both treatment assignment and adherence to be independent of all potential outcomes for any possible value of adherence $a$. This assumption must be satisfied for any level of adherence; otherwise, individuals would be able to select the treatment option with the highest reward for a given value of adherence below unity. However, it only needs to be satisfied for the true value of the state variable $s_i$ since it is implicitly assumed throughout the paper that the state variable is exogenously determined.\footnote{This assumption might however be violated if the state variable becomes a function of the assigned treatment over its duration. Because the treatment options I consider in the empirical application are relatively short (up to three months), I rule out this possibility for all observations in the sample. It is also more realistic to assume the exogeneity of the state variable rather than the exogeneity of all elements in $X_i$, which is implicit under the typical unconfoundedness assumption (see, for instance, hunermund_double_2023 for more details).} Instead of relying on standard unconfoundedness, this paper employs a slightly different and more practical assumption that is formulated below in Assumption (ref).

Assumption (ref) supposes that the observed treatment choice and adherence are exogenous within each group formed by the true value of the state variable. This assumption complements Assumption (ref) and is necessary to obtain consistent estimates of all within-group CATEs. Note that Assumption (ref) corresponds, in fact, to a special case of Assumption (ref) : the two assumptions become identical if $S = |\mathcal{X}|$, where $|\mathcal{X}|$ represents the total number of possible values taken by any $x \in \mathcal{X}$. If at least one element in $\mathcal{X}$ is continuous, then $|\mathcal{X}|$ is infinite and local approximations of Assumption (ref) are instead required to consistently estimate all CATEs in finite samples. Therefore, Assumption (ref) is more practical than Assumption (ref) without being stronger than the latter since it automatically adjusts the conditioning set to the minimum amount of information that is required for conditional exogeneity to hold. The fact that Assumption (ref) corresponds to a special case of Assumption (ref) implies that both assumptions will be violated under the same circumstances, namely when treatment selection and/or adherence to treatment is driven by the patient's idiosyncratic response to treatment after conditioning on $X_i$ or $s_i$. This will occur, for instance, when patient-specific response to treatment is unbiasedly anticipated (by the patient and/or the treating physician), and adherence or treatment choice is then adjusted accordingly. For instance, every patient who responds negatively to a given treatment may exhibit, during the course of the treatment, side effects that are not predictable by any set of observable characteristics. If these side effects generally induce the treating physician to stop the current treatment regimen or switch to another one, then Assumption (ref) will be violated, and the estimated coefficient relating adherence to observed outcomes will be positively biased lamiraud_therapeutic_2007. This means that Assumption (ref) is equivalent to assuming that either patient-specific responses to treatment are predictable by observables, or that anticipated responses to treatment do not systematically influence the patient's or the physician's behavior. A formal example showing under which conditions Assumption (ref) will and will not hold under additive group-wise heterogeneity is presented in Appendix (ref). The main practical difficulty concerning Assumption (ref) is that each $s_i$ is unobserved. Therefore, it is necessary to obtain consistent estimates of the true state values $\mathbf{s} = (s_1,...,s_i,...,s_N)^{\top}$ to be able to recover all CATEs. To do so, I rely on the next six assumptions.

Assumption (ref) is standard in the PL literature and supposes that all potential outcomes are uniformly bounded both from above and below mbakop_model_2021, zhou_offline_2023. Assumption (ref) states that all treatment options must feature a non-zero probability of being assigned to anyone in the sampled population given the true value of the state variable. In other words, no formal contraindication must prevent the assignment of any treatment option to all observations in a given group. This assumption is weaker than the typical overlap assumption since it is conditional on the state value rather than on the value of the observed covariates. It thus allows for the presence of formal contraindications at the individual level that do not translate at the group level (although this might be rare in practice). Assumption (ref) states that each estimated within-group outcome model is $\sqrt{\tilde{N}_g}$-consistent if and only if every observation is correctly classified as the sample size tends to infinity, where $\tilde{N}_g$ corresponds to the true number of observations in group $g$. This assumption is also weaker than the typical $\sqrt{N}$-consistency assumption made in the semiparametric estimation literature due to its conditional nature, but also because it does not assume uniform consistency of $\hat{m}_{d,g}(a,x)$ as $N \to \infty$ across all data-generating processes. For instance, if $\tilde{N}_g = N^{\lambda_g}$ with $0 < \lambda_g < 1$ for a given $g \in \mathcal{S}$, Assumption (ref) still holds even if the rate of convergence of $\hat{m}_{d,(z_i)}(a,x)$ is slower than $N^{-1/2}$. In the current form, Assumption (ref) however excludes estimators of $\hat{m}_{d,(z_i)}(a,x)$ that converge at a slower rate than $\tilde{N}_g^{-1/2}$. Relaxing such an assumption could be done, for instance, to allow for the nonparametric estimation of the regression function $m_{d,g}(a,x)$ for every $g \in \mathcal{S}$, but at the cost of a slower convergence rate for the estimated CATEs (see Corollary (ref)). Assumption (ref) may also be too restrictive in high-dimensional settings where each $\hat{m}_{d,z_i}(a,x)$ is modeled semiparametrically. As shown by chernozhukov_doubledebiased_2018, such an approach will lead to regularization bias if the estimated nuisance part converges at a slower rate than the standard parametric rate. To overcome this bias, the authors suggest estimating all parameters of interest using orthogonalized score functions and cross-fitted estimators.\footnote{Cross-fitting is a sample splitting approach that limits the incidence of regularization biases on the parameters of interest. See Section 3 of chernozhukov_doubledebiased_2018 for more details.} In contrast, Assumption (ref) directly implies that all nuisance parameters are $\sqrt{\tilde{N}_g}$-consistent, which implies the absence of any regularization bias. This assumption is realistic in settings where the space of relevant predictive variables remains low-dimensional and linear within each group, which is a credible assumption for several drug-health outcome relationships. The realism of this assumption is confirmed in the empirical application, where the set of nuisance parameters that maximizes OOS prediction accuracy corresponds to the empty set within each group. Extending the empirical application to other outcomes might however invalidate Assumption (ref), and therefore require the use of the double/debiased approach with cross-fitting. Nonetheless, sample splitting is used in the empirical application to obtain valid inference on all estimated CATEs (see Section (ref)). If the selected classifier can correctly classify all observations in the sample as $N \to \infty$, then it is said to be uniformly consistent. Uniformly consistent classifiers are typically hard to obtain without making strong assumptions on the data-generating process su_identifying_2016,dzemski_convergence_2021, langevin_bias-reduced_2026. In practice, it is also more frequent to use a consistent classifier, which implies that the proportion of misclassified observations goes to zero in the limit. In this case, the rate of convergence of the difference $\hat{m}_{d,(z_i)}(a,x) - m_{d,g}(a,x)$ will typically be slower than $\tilde{N}_g^{-1/2}$, but the estimates will remain consistent as $\tilde{N}_g \to \infty$. This case is believed to be more realistic for most empirical applications, including the one presented in this paper.\footnote{Modifying Assumption (ref) such that $\hat{m}_{d,(z_i)}(a,x) - m_{d,g}(a,x)$ is $\sqrt{\tilde{N}_g}$-consistent if and only if $z_i = (s_i)$ for almost every observation in the sample implies a faster rate of convergence than the standard parametric rate for each outcome model taken separately, which is too strong an assumption to make in most empirical settings.} Assumptions (ref) and (ref) are similar to Assumptions 2$(iv)$ and 2$(v)$ of langevin_bias-reduced_2026 respectively, and correspond to sufficient requirements on the groups' structure to guarantee the consistency of the K-means procedure under correct specification of the total number of groups. Assumption (ref) may be relaxed in favor of finite nested conditional moments, where both $\mathbb{E}[X_{i}] = \boldsymbol{\mu}^0_g$ and $\text{Var}[X_{i}] = \Sigma^0_g \Leftrightarrow g = s_i^* \in \mathcal{S}^*$, where $\mathcal{S}^* = \{1,2...,S^*\}$ with $S^* > S$ and $\mathbbm{1}[s_i^*=g] \Rightarrow \mathbbm{1}[(s_i^*)^*=s_i]$, and where $(g)^*:\mathcal{S}^* \rightarrow \mathcal{S}$ is a surjective permutation in the group labels. Specifically, the unobserved heterogeneity in the covariates can be more important than the one described by Assumption (ref), as long as the former can be mapped into the latter by merging some groups of covariates. The inverse is however not possible without violating Assumption (ref), meaning that the heterogeneity in the treatment effects cannot be larger than the heterogeneity in the covariates. Note that Assumption (ref) requires each true covariance matrix $\Sigma_g^0$ to have only a few different elements across groups to be satisfied. However, the overall misclassification rate is likely to be smaller as the between-group differences in $\Sigma_g^0$ become larger for a given value of $p$. Finally, Assumption (ref) supposes that the size of each group will go to infinity as the total sample size goes to infinity, but leaves the growth rate of each group unspecified. It is worth noting here that the implicit assumption that $S$ is a known and fixed population parameter can be problematic in small samples, especially if some groups represent a very small share of the whole population. In practice, this will lead to generalizability issues only if the sample size is too small such that the parameters associated with one or several groups in the population cannot be identified from the sample. The probability that such an issue arises in practice is highly context-dependent and is believed not to occur in the considered empirical application given the representativeness of the cohort and the relatively large sample size (i.e., $N=1,232$). On the other hand, it is sometimes realistic to consider the true number of groups to be a function of the sample size if more nuanced differences across observations are likely to be identified as $N$ increases. Such an assumption can practically coexist with the former one if all analyses are performed conditional on the sample size. Although there might exist a (much) larger number of groups in the population, it becomes approximated by the “true" number of in-sample groups $S$, which in fact combines true population groups that are similar enough not to lead to significant biases in the estimates. The main theoretical results of the paper are based on the following definitions.

I can now state the main theoretical results of this paper. All proofs are provided in Appendix (ref). Note that all theoretical results implicitly assume that the econometrician knows the true number of groups, $S$, which can be estimated using $K$-fold cross-validation, or another technique from the relevant literature tibshirani_estimating_2001,sugar_finding_2003.

Lemma (ref) shows that using the weighted K-means algorithm never increases the TSMD between two consecutive iterations. This lemma is very practical for searching the global minimum of the TSMD across the entire parameter space, $\mathcal{Z} \times \Theta^{S} \times \Omega^S$. The addition of the logarithm of the determinant of $\Sigma_g$ in the definition of the squared Mahalanobis distance ensures that the TSMD will never increase when going from $\Sigma^{(k)}_g$ to $\Sigma^{(k+1)}_g$. However, this lemma requires every $\Sigma^{(k)}_g$ to be positive-definite at each iteration of the algorithm. This may pose a problem in finite samples if $p > N^{(k)}_g$, where $N^{(k)}_g$ denotes the number of observations assigned to the $g^{th}$ group at the $k^{th}$ iteration. The matrix $\Sigma^{(k)}_g$ will also be singular if some covariates show no within-group variation or near multicollinearity among each other, which is likely to occur in small groups and/or if $X_i$ contains binary or discrete elements. To solve this type of problem, I use the eigenvalue decomposition of $\Sigma^{(k)}_g$ to enforce a minimum, strictly positive eigenvalue for every group at every iteration of the weighted K-means algorithm. Doing so guarantees that the resulting covariance matrix will be invertible and positive-definite. This can however invalidate the conclusion of Lemma (ref) as the resulting set of covariance matrices will no longer correspond to the global minimizer of the TSMD. In the empirical application, the objective function decreases between two consecutive iterations for all sets of initial parameter values, except when the enforced minimum eigenvalue, denoted by $\lambda_{min}$, is set at a relatively high level (i.e., more than 1.5). This is why I recommend to set $\lambda_{min}$ as close to zero as possible in order to limit the impact of this restriction on the convergence behavior of the algorithm.\footnote{Note also that the weighted K-means algorithm reduces to the standard K-means algorithm when $\lambda_{min} \ge \max \{\lambda_{1,max},...,\lambda_{S,max}\}$ at each iteration, where $\lambda_{g,max}$ is the maximum eigenvalue of $\Sigma_g$ for a given $g \in \mathcal{S}$. In this case, it is easy to show that $d_i^{2M}(\boldsymbol{\mu}_{g}, \Sigma_g) = \lambda_{min}^{-1} (X_i - \boldsymbol{\mu}_{g})^{\top} (X_i - \boldsymbol{\mu}_{g}) + p\log (\lambda_{min})$, and that minimizing $d^{2TM}(\mathbf{z},\boldsymbol{\mu}, \Sigma)$ with respect to $\mathbf{z}$, $\boldsymbol{\mu}$, and $\Sigma$ will yield the same results as the standard K-means algorithm. Such an equivalence also highlights the fact that $\lambda_{min}$ determines the maximum weight that can be given to any covariate introduced in the weighted K-means algorithm.} The main benefit of using the minimum eigenvalue approach is that it allows the use of all available covariates to estimate all group memberships. As shown in Theorem (ref), the probability of misclassifying any observation in the sampled population goes to zero as the number of covariates goes to infinity. This theorem is similar to Theorem 1 of langevin_bias-reduced_2026 for a finite number of groups and offers a clear advantage to the Mahalanobis distance over the Euclidean distance for classification tasks under mild regularity conditions. Consequently, it becomes crucial to include as many covariates as possible in the weighted K-means algorithm to bring down the misclassification rate to zero in finite samples, even though this might lead to $p > N^{(k)}_g$. This is why the minimum eigenvalue approach is preferred over other approaches based on l$_1$/l$_2$-regularization, or on direct estimation of $\{\Sigma^{(k)}_g\}^{-1}(X_i - \boldsymbol{\mu}^{(k)}_{g})$ cai_chime_2019. Note also that $p \to \infty$ does not necessarily violate Assumption (ref) if $p/N \to 0$ as $N,p \to \infty$, or if the set of relevant predictive variables is different from the set of relevant prescriptive variables. Theorem (ref) states that the asymptotic location of the global minimum of the TSMD is located at the true mean values and covariance matrices if almost every observation in the sample is correctly classified asymptotically. Note that the asymptotic global minimum of the TSMD is unique under Assumptions $\ref{ass:7}$-$\ref{ass:8}$, unless two groups (or more) share the same vector of mean values and the same covariance matrix, which is (very) unlikely when $p$ is (very) large. Combining the results of both theorems implies that globally minimizing the TSMD leads to consistent estimates of all mean values, covariance matrices, and group memberships (Corollary (ref)). Note also that Corollary (ref) does not require $p$ to increase at a faster rate than $N$, which limits the occurrence of the issues mentioned above under Assumption $\ref{ass:8}$. In fact, $p$ can be much smaller than the number of observations within each group without sacrificing consistency, even though some observations may remain misclassified asymptotically. Corollary (ref) requires that $p$ increases at a faster rate than $N$ to obtain $\sqrt{\tilde{N}_g}$-consistent estimates of all CATEs under Assumption (ref). If this is not the case, the rate of convergence of the CATEs in the $g^{th}$ group might be slower than $\tilde{N}_g^{-1/2}$, but the convergence of the weighted K-means algorithm will be facilitated if the number of covariates remains relatively small. This shows that the proposed approach entails a trade-off between efficiency and practicability under Assumption (ref) and Assumptions (ref)-(ref). It is however possible to have both $p < N$ and $z_i = (s_i)$ for every $i \in [N]$ in finite samples under the same set of assumptions. Such a situation will occur if $\boldsymbol{\mu}^0_g$ and $\Sigma^0_g$ vary sufficiently across groups such that it will bring classification error down to zero, and therefore lead to efficient estimation of all CATEs under Assumption (ref). Overall, the above results show that minimizing the TSMD through the weighted K-means algorithm leads to the recovery of almost every true group membership under Assumptions (ref)-(ref). Predicted group memberships can then be used to obtain consistent estimates of all CATEs. The precision of the estimates will depend on both the “quality and quantity" of the covariates introduced in the weighted K-means algorithm. To note, there is no guarantee that the global minimum of the TSMD lies close to the set of true group memberships in finite samples with a relatively small number of covariates, or when most covariates do not satisfy Assumption (ref). This is why I use 5-fold cross-validation with several data splittings to validate the results obtained in the empirical application presented in Section (ref). The proposed methodology does not require the estimation of the treatment selection model $\Pr[A_{i,d}>0|s_i=g,X_i=x], \ \forall (d,g,x) \in \mathcal{D} \times \mathcal{S}\times \mathcal{X}$, as long as Assumptions (ref)-(ref) are satisfied and all treatment options are observed within each group. Estimating $\Pr[A_{i,d}>0|s_i=g,X_i=x]$ is useful when Assumption (ref) is violated and orthogonalizaton of the adherence/assigned treatment variable is required (as shown in chernozhukov_doubledebiased_2018), or when a doubly-robust approach is considered athey_policy_2021,zhang_subgroup_2022,li_optimal_2023. The latter is usually performed by plugging the different CATEs into a “score" vector that is based on augmented inverse probability weights (AIPW). However, it is unclear how the double-robustness property is maintained in the multi-action framework when the vector of estimated scores is defined as follows zhou_offline_2023,zhan_policy_2024

for each $d \in \mathcal{D}$, where adherence is set equal to $A_{i,D_i}$ across all treatment options for simplicity. When $d \ne D_i$, the $d^{th}$ element of $\hat{\Gamma}_{i,g}$ is equal to $\hat{m}_{d,g}(A_{i,D_i},X_i)$, whose consistency will depend on the chosen outcome model regardless of how one models the conditional selection process $\Pr[A_{i,d}>0|s_i=g,X_i=x]$ for each $d \in \mathcal{D}$. Moreover, the fact that the realized outcome $Y_i$ is included in the definition of the estimated score $\hat{\Gamma}_{i,g}$ necessarily requires an additional prediction step to identify the policy rule (e.g., Equation (38) of athey_policy_2021 and step 7 in Algorithm 1 of zhou_offline_2023) when the class of feasible policies $\Pi$ is left unconstrained. If the first term on the RHS of Eq. ((ref)) is dropped, then the policy rule can be written as $\hat{\pi}^*(X_i) = \arg \max_{d\in \mathcal{D}} \hat{m}_{d,g}(A_{i,D_i},X_i)$ for a given value of $g \in \mathcal{S}$. In contrast, the methodology proposed in this paper does not rely on the inclusion of the observed or predicted outcome value in the weighted K-means algorithm, which is believed to reduce generalization error and improve OOS prediction accuracy. Using the estimated residuals from the outcome models $\hat{m}_{d,g}(A_{i,D_i},X_i)$ to identify the CATEs also creates a form of redundancy in the weighted K-means algorithm, where the estimated CATEs are present in both the inputs and the output of the algorithm, and therefore have to be removed for test samples where the outcome is not yet revealed. Sensitivity analyses also indicate that doing so leads to worse OOS prediction accuracy in the empirical application (results not shown). Other approaches based on nonparametric estimators and statistical learning, such as kernel regression, sieve-based methods, random forests, and neural networks, can certainly be used to estimate CATEs foster_subgroup_2011,athey_policy_2021. However, those algorithms/estimators mostly rely on the standard unconfoundedness assumption and/or might not offer good OOS performance if the size of the training dataset is not large enough. Moreover, most of the literature on PL focuses on the discrete treatment case, whereas the empirical application considered in this paper is a multi-action PL problem with a semicontinuous treatment variable for each action (i.e., the observed adherence). Despite this, empirical results obtained with the proposed approach are compared to the ones obtained by an “off-the-shelf" approach based on random forests under discrete enforcement of the treatment variable (see Section (ref) for more details). Additional comparisons between the proposed approach and other ML techniques are however left to further research.

Incorporating Imperfect Adherence to Treatment

If adherence to treatment can be imperfect in practice, it is important to account for this feature when learning a policy rule from observational data. Adherence is modeled as a function of expected treatment effectiveness and of a set of observed individual and environmental features (not necessarily a subset of $\mathcal{X}$). Taking on the Bayesian perspective, it is thus possible to modify the policy rule as follows

where $a(\cdot,\cdot)$ denotes the adherence function, $W_i \in \mathcal{W} \subseteq \mathbb{R}^{p'}$ with $p'>0$, and where $\mathbb{E}_S[\tau_{\pi(X_i),s_i}(1)]$ refers to the expected effectiveness of the selected treatment with perfect adherence. Under Assumption (ref), it is possible to rewrite this equation as follows

where it is now explicit that the policy rule $\pi^{*}_{B}$ depends directly on the CATE $\tau_{d,g}(\cdot)$. The choice of elements contained in $W_i$ can be determined using, for instance, LASSO-based regularization procedures or theoretical insights from the relevant literature. In practice, this results in a two-step estimation procedure where the first step uses observed adherence $A_{i,D_i}$ to estimate all CATEs via the approach described in the previous section. The second step models the observed adherence $A_{i,D_i}$ as a function of $W_i$ and predicted treatment effectiveness under perfect adherence for each treatment option. This last quantity can be estimated for each $i \in [N]$ with

for each $d \in \mathcal{D}$ using the estimated CATEs $\hat{\tau}_{d,g}(\cdot)$ from the first step and a subjective prior over $s_i$. Introducing predicted treatment effectiveness in the model allows adherence to depend on the patients' or the physicians' anticipations regarding treatment effectiveness, which is a natural choice if adherence is assumed to be unconditionally endogenous.\footnote{As opposed to being conditionally endogenous, which would violate Assumptions (ref)-(ref).} Indeed, a poor response to treatment is likely to create negative anticipations over time, which will entail lower adherence values, hence leading to a positive, endogenous correlation between observed adherence and potential outcomes. Predicted treatment effectiveness thus acts as a proxy for anticipated treatment effectiveness and should not be omitted when modeling adherence. Given that observed adherence is typically a semi-continuous variable whose distribution is bounded between 0 and 1 with a probability mass at 1 exactly, I recommend modeling adherence using a two-part model leung_choice_1996. The first part models whether or not the $i^{th}$ individual completely adheres to treatment, whereas the second part models the degree of adherence in case the $i^{th}$ individual does not completely adhere to treatment. Once a proper adherence model $\hat{a}(W_i,\hat{\tau}_{d}(1))$ is estimated, predicted adherence values can be plugged into the outcome model in order to obtain the Bayesian version of the estimated policy rule

where $\hat{\tau}_{d,g}(\cdot)$ is defined as in Corollary (ref), and where $\hat{\tau}_{d}(1)$ is defined as in Eq. ((ref)). For simplicity, the adherence model is assumed to be independent of the state variable. This can be relaxed via the estimation of within-group adherence models provided sufficient variation in observed adherence across groups.

Identifying a Feasible Policy Rule

Under Assumptions (ref) and (ref), and no global constraint over the class of feasible policies $\Pi$, the infeasible (and globally optimal) policy rule is the one we would obtain if $s_i$ were known for all observations. In this case, we can define the estimated version of the infeasible policy rule for observation $i$ as follows

where I now denote the estimated policy rule for the $i^{th}$ observation as a function of $W_i$, and where $\hat{\tau}_{d}(1)$ is defined as above. If predicted adherence does not significantly differ across treatment options in the sampled population, then optimal treatment choices depend only on the value of $s_i$ for all $i \in [N]$. Since $s_i$ is unobserved, replacing $s_i$ by $\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ for each observation in the sample leads to an efficient policy rule as a result of Theorem (ref), where $(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ correspond to the global minimizers of the TSMD. Because $\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ is a consistent estimator of $s_i$ for any $i \in [N]$, we can treat $\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ as the true group membership value for every observation in the sample. However, replacing $s_i$ by $\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ typically requires the use of a large number of covariates to ensure that the number of misclassified observations is close to zero in finite samples. This is why I still call this policy rule “infeasible" given that it is based on too large an amount of information to be practical in most clinical contexts. I define the feasible policy rule as the rule that optimally predicts $\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ using only a subset $V \in \mathcal{V}_K \subset \mathcal{X}$, where $\mathcal{V}_K$ is the space of all subsets of $\mathcal{X}$ such that $|V| = K < |X|$ with $|A|$ denoting the size of the $A$, and with $K$ being a small positive discrete number (often less than 10). If we consider $\mathcal{H}$ as the class of all real-valued functions $h: \mathcal{V}_K \to \mathcal{S}$, then we can define the following misclassification rate

where $V = (V_1,...,V_i,...,V_N)^{\top}$, with $V_i$ corresponding to the $K$-sized vector of selected covariates for individual $i$. It is then possible to define the optimal function $h^*(\cdot) \in \mathcal{H}$ with optimal subset $V^*$ as follows

Combining this definition with the definition of $\hat{\pi}^{I}_i$ gives us the estimated feasible policy rule for the $i^{th}$ observation

where the subscript $C$ refers to “Certainty", as in Section (ref). It is also possible to apply the Bayesian decision criterion to obtain the Bayesian version of this feasible policy rule for observation $i$, which is written as follows

where $\Pr[\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma}) = g|V^*_i=v]$ is a data-dependent prior that defines the probability of belonging to group $g$ conditional on the value taken by $V^*_i$ for any pair $(i,g) \in [N]\times \mathcal{S}$. This last policy rule has the benefit of accounting for the uncertainty coming from the fact that only the subset of characteristics $V^*$ is used to identify group membership. Even though the function $h^*(\cdot)$ and the subset $V^*$ are optimal in terms of the minimization of the misclassification rate, they will yield suboptimal results compared to using $\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ in the policy rule, unless the true group memberships are entirely determined by the subset $V^*$ and there exists a function $h^*(\cdot)$ such that the misclassification rate $\text{Mis}(h^*(V^*),V^*)$ is equal to zero in general. Finally, it is possible to (re)introduce global constraints on the class of feasible policies $\Pi$ by summing $\hat{\pi}^{F}_{B,i}$ over all observations as follows

where $d_i$ is the selected treatment option for observation $i$, and where $\Pi \subset \mathcal{D}^N$ such that all global constraints are satisfied. Note that Eq. ((ref)) and Eq. ((ref)) can be easily modified to include other decision criteria (i.e., maximin or minimax-regret). In practice, the most difficult part is to determine which function $h(\cdot) \in \mathcal{H}$ and subset $V \subset X$ are such that the misclassification rate Mis$(h,V)$ is minimized. To do so, several authors have relied on decision trees given their flexibility and interpretability sverdrup_policytree_2020,athey_policy_2021,zhou_offline_2023. Other methods, such as simple linear regression or multinomial logistic regression, can also be used, but they tend to be less flexible than decision trees and are not necessarily easy to use for routine decision-making processes.

Empirical Application

This section is divided into four subsections. The first subsection briefly describes the cohort used for the empirical application. The second describes the health condition shared by the cohort's participants and the different treatment options currently used to treat this condition. The third subsection briefly describes the clinical guidelines concerning HepC treatment in both mono- and HIV/HCV co-infected patients. Finally, the fourth subsection then describes the details of the employed PL strategy.

The Canadian Co-Infection Cohort

The Canadian Co-Infection Cohort (CCC) prospectively follows around 2,200 patients with documented HIV infection and chronic Hepatitis C (HepC) infection or evidence of HCV exposure across 6 Canadian provinces between April 2003 and December 2023.\footnote{A chronic HepC infection, as opposed to an acute HepC infection, refers to an infection that does not spontaneously clear 6 months after a first positive HCV blood test.} Inclusion criteria into the cohort are very broad : $>$ 16 years old, documented HIV seropositive infection, documented HepC infection or evidence of exposure, and capacity to provide informed consent. All eligible patients were approached to participate in order to avoid selection bias klein_cohort_2010. A total of 19 clinical sites are participating in the CCC as of the end of 2023. Figure (ref) shows the geographical distribution of all participants enrolled in the CCC between 2003 and 2023.

The primary objective of the CCC is to evaluate the effect of highly active antiretroviral therapy (HAART) on the progression of end-stage liver disease (ESLD) in HIV/HCV co-infected individuals, while the secondary objective is to evaluate the “role of HCV treatment in the evolution of liver disease with a particular emphasis on evaluating access to treatment" klein_cohort_2010. To meet these objectives, socio-demographic, medical, behavioral, and patient-reported quality of life information is collected for all participants through questionnaires at baseline and during follow-up visits. Follow-up visits are planned every 6 months to update information on risk behaviors, medical treatments, and diagnoses. Blood tests are also performed at each follow-up visit, which includes liver profile (i.e., concentration of various liver enzymes), blood cell count, hepatitis A and B serology, plasma HIV and HCV RNA, and HCV genotype. Causes of death, withdrawals from the study, and losses to follow-up (defined as missing more than 2 consecutive follow-up visits) are also systematically recorded. It is estimated that the CCC has enrolled 23% of the total HIV/HCV co-infected population in care in Canada saeed_frequent_2024. The CCC is representative of the co-infected population engaged in clinical care young_rate_2023. Figure (ref) presents the socioeconomic characteristics of CCC participants at enrollment. Note the large proportion of Indigenous people and low-income individuals (less than CAN\$1,500 per month) compared to the overall Canadian population.\footnote{According to Statistics Canada's Canadian Income Survey, around 15% of the Canadian population lives on less than CAN\$1,500 per month statistics_canada_table_2024. The 2021 Census data also revealed that Indigenous people comprise about 5% of the overall Canadian population government_of_canada_canadas_2023.} Because HCV and HIV are blood-borne diseases that spread mostly through injection drug use, low-income Indigenous people are at increased risk of being co-infected. In Canada, evidence suggests that Indigenous peoples account for 70% to 80% of new HepC infections among individuals who inject drugs saeed_frequent_2024.

HIV, HCV, and Direct-Acting Antiviral Agents

Around 6 million people are believed to be co-infected with HIV and HCV worldwide soriano_changing_2013. Since the approval of the first antiretroviral therapy (ART) in 1987, several other antiretroviral drugs have been approved and used to treat HIV infection. This has led to a dramatic reduction in HIV morbidity and mortality in the developed world tseng_evolution_2015. HepC is a viral infection that causes liver inflammation, which can lead to lifelong illnesses such as liver cirrhosis and cancer. In Canada, it is estimated that more than 200,000 persons were living with chronic Hepatitis C infection as of the end of 2021, with an incidence of nearly 8,000 new cases per year. HepC-related mortality has been decreasing steadily since the approval of the new direct-acting antiviral (DAA) agents by Health Canada in 2014 public_health_agency_of_canada_canadas_2024. Although DAAs have been covered by public health insurance plans in all Canadian provinces since 2018, recent evidence suggests that treatment initiation rates are too low to reach the HCV reduction targets in Canada by 2030 young_rate_2023, feld_timing_2022, tian_feasibility_2024.\footnote{The Canadian HCV elimination targets aim to reduce chronic hepatitis C infections by 80% and related deaths by 65% between 2016 and 2030.} Systemic barriers such as incarceration and discrimination toward Indigenous people are believed to represent some of the main causes for current low DAA uptake despite universal health insurance coverage in Canada jiang_hcv_2023, hosein_increased_2024. Therefore, extensive scale-up in screening, treatment, and harm-reduction strategies may be necessary to reach the Canadian targets tian_feasibility_2024. Note that, unlike Hepatitis B, there is no vaccine against HepC, which implies that individuals can be reinfected after being completely cleared of HepC viral load. In HCV mono-infected patients, clinical trials have shown more than 90% efficacy of DAAs in achieving SVR 12 weeks after treatment, regardless of HCV genotype esteban_efficacy_2018,laurain_safety_2019. As mentioned in the introduction, reaching SVR is equivalent to being cured from the disease and is the health outcome of interest in most studies on the efficacy of HCV treatments. Similar results have been obtained in trials with both HCV mono-infected and HIV/HCV co-infected patients ford_expanding_2012, rockstroh_efficacy_2015. However, the generalizability of those results outside the clinical trial setting is not obvious due to the very restrictive eligibility criteria used in most DAA trials. Specifically, uptake of certain antiretroviral therapy regimens, active illicit drug use, low CD4 cell count, and detectable HIV RNA are the reasons that lead to the exclusion of co-infected patients from most DAA trials saeed_how_2016. In parallel, observational studies seeking to estimate the effectiveness of DAAs in co-infected patients obtain SVR rates that are between 10 and 17 percentage points lower for co-infected persons compared to HCV mono-infected persons alam_real-world_2017, wilton_real-world_2020. The lower effectiveness of DAAs in co-infected individuals can be attributed to different factors. A recent study looking at 4098 co-infected individuals across 6 different countries found that unsuccessful DAA treatment was more likely among individuals with evidence of higher immunodeficiency and among those who reported recently injecting drugs harney_unsuccessful_2025. Other evidence suggests that the non-achievement of SVR for persons who initiate DAA therapies is mainly driven by nonvirological factors such as treatment discontinuation, nonattendance for viral load testing after treatment initiation, or death darvishian_loss_2020, wilton_real-world_2020. Other elements, such as drug-drug interactions (DDIs) with current and/or past ART, high degree of hepatic fibrosis, and HCV genotype, might also affect DAA effectiveness in co-infected individuals chen_treatment_2015, kwo_regimens_2015, martel-laferriere_prise_2022.

DAA Regimens and Clinical Guidelines

Three main DAAs are currently used to treat HepC infection in co-infected individuals, with the following commercial names : Epclusa (EPC), Mavyret (MAV), and Harvoni (HAR). Both the Epclusa and Harvoni treatment regimens consist of a single pill to be taken daily for 12 weeks. The Mavyret treatment regimen consists of three pills to be taken daily for 8 weeks. The American Association for the Study of Liver Diseases (AASLD), which provides guidelines and practical recommendations to health practitioners in the United States regarding HCV treatment, indicates that Harvoni features fewer suspected DDIs with ARTs than Epclusa, which in turn features fewer suspected DDIs with ARTs than Mavyret american_association_for_the_study_of_liver_diseases_patients_2022. This may be part of the reason why the price negotiated by pan-Canadian Pharmaceutical Alliance for Harvoni ($\approx$ CAN\$798/pill) is higher than Epclusa's ($\approx$ CAN\$714/pill), which is in turn higher than Mavyret's ($\approx$ CAN\$238/pill).

For patients who do not have any one of the following : a current Hepatitis B infection, an ongoing pregnancy, a current hepatocellular carcinoma (HCC) diagnosis, a history of liver transplant, and who were never treated with DAAs or older interferon-based treatment regimens, Canadian guidelines recommend prescribing either Epclusa or Mavyret for their respective length, regardless of the HCV genotype martel-laferriere_prise_2022-1. On the other hand, the AASLD guidelines recommend extending Mavyret's regimen length to 12 weeks in case of a co-infection with HIV. Unlike the AASLD guidelines, Canadian guidelines do not distinguish between HCV mono-infected and HIV/HCV co-infected individuals american_association_for_the_study_of_liver_diseases_patients_2022. Given the general lower SVR rates for HIV/HCV co-infected patients, the absence of specific Canadian guidelines for co-infected patients might lead to suboptimal health outcomes within the Canadian co-infected population. Finally, the AASLD recommends using the Harvoni treatment regimen only against certain HCV genotypes (1a, 1b, 4, 5, and 6), whereas Canadian guidelines recommend using Harvoni only in patients where interferon-based treatment regimens have not been successful in the past martel-laferriere_prise_2022-1. Since there is no vaccine for HepC, it is not possible to prevent the spread of the virus through mass vaccination campaigns. Moreover, it is estimated that 20% to 30% of Canadians infected with HIV are co-infected with HCV. This percentage is likely to rise in the future due to higher reinfection rates in co-infected individuals, especially among those who inject drugs young_rate_2023,hosein_increased_2024. Therefore, it is essential to treat every co-infected patient with the most effective option in order to maintain the HCV transmission rate as low as possible in the future. Given that the marginal effect of each DAA regimen on the probability of reaching SVR in co-infected patients is uncertain, the proposed estimation strategy can be used to derive a policy rule that will identify which DAA among the three listed above is the most effective for any given patient with observed characteristics $X_i$ subject to a set of predetermined (cost) constraints. The next section details such a strategy.

Estimation Strategy

Sampled Population

The sampled population consists of all HCV-infected CCC participants (regardless of the enrollment date) who were still alive and not lost to follow-up before January 1st, 2014, and who were prescribed either Epclusa, Mavyret, or Harvoni after January 1st, 2014, or who were never prescribed any treatment against HCV between January 1st, 2010, and September 1st, 2024. These never-treated individuals will act as the set of “untreated individuals" to which every treatment option will be compared, and allows for the fact that any individual infected with HepC can spontaneously clear the virus. The outcome $Y_i$ corresponds to whether an individual has reached SVR between January 1st, 2014, and September 1st, 2024. Individuals who were treated multiple times with DAAs or with older interferon-based treatment regimens and at least one DAA regimen between January 1st, 2014, and September 1st, 2024, are not excluded from the sample. This implies that the same treated individual can be observed multiple times in the sampled population.\footnote{Because this situation concerns only 39 patients out of 1,187, all within-patient outcome values were assumed to be independent from each other after conditioning on a set of dummy variables indicating all previously assigned treatment options. Observations with missing or unknown outcome values were also discarded from the analysis (less than 20 observations out of 745 remaining treated observations).} SVR outcomes for both treated and untreated patients are provided by the participating centers and are confirmed via a careful examination of the laboratory test results. Additional details concerning the identification of spontaneous clearance in untreated patients and SVR outcomes in treated patients are described in Appendix (ref).

Classification Procedure and Outcome Models

I use a total of 146 pre-treatment variables to perform the classification procedure, including standard individual characteristics (e.g., age, gender, weight, etc.), socioeconomic indicators (e.g., income brackets, educational attainment, etc.), medical diagnoses, past and current ART medications, previous HCV treatment regimens (and the respective SVR outcome), and other health-related information. Detailed information about the construction of the set of variables entering the classification procedure is provided in Appendix (ref). I then introduce this set of variables into the weighted K-means algorithm to minimize the TSMD, as shown in Definition (ref)$(iii)$ and Algorithm (ref). Because the TSMD is a non-convex function of the estimated parameters and group memberships, I use a thousand initial random assignment values (i.e., $\mathbf{z}^{(0)}$) to search for the global minimum of the TSMD for each combination of hyperparameters in a large grid of realistic hyperparameter values. To limit the influence of missing values on the classification procedure and given the large number of missing values in the CCC data, I impute missing covariates' values with the group mean values at each iteration of the algorithm, which allows me not to drop any observations during the classification procedure.\footnote{One caveat of such an imputation scheme is that it underestimates (co)variance parameters, which is however compensated through the enforcement of the minimum eigenvalue $\lambda_{min}$. Direct adjustment to appropriately “inflate" the estimated (co)variance parameters led to several convergence issues, which is why I did not consider this approach in the empirical application. Note also that the number of missing values varies across covariates, from a few missing values to several hundreds.} The weighted K-means algorithm is governed by two hyperparameters : the total number of groups $S$, and the minimum eigenvalue of the variance-covariance matrix $\Sigma_g$ for all $g \in \mathcal{S}$, $\lambda_{min}$. As mentioned earlier, enforcing a given value of $\lambda_{min} > 0$ guarantees that every $\Sigma_g$ remains positive-definite at each iteration of the algorithm, which is necessary for the algorithm to converge. Increasing the value of $\lambda_{min}$ reduces the weight associated with variables that feature a relatively small degree of within-group variation. Since some groups might be defined by a specific value for a single or a set of covariates, it is important not to automatically discard variables that show little or no within-group variation. Enforcement of positive-definiteness through a common minimum eigenvalue allows me to do so. I model the SVR outcomes $Y_i$ via a linear probability model (LPM) within each group identified by the weighted K-means algorithm. LPMs are convenient in this context given that : 1) they are more likely to meet Assumption (ref) than nonlinear binary choice models; 2) they easily deal with imperfect adherence since $\hat{\tau}_{d,g}(a) \equiv \hat{\tau}_{d,g}(1)\times a$ for any $a \in [0,1]$, where $\hat{\tau}_{d,g}(1)$ is estimated using observed adherence $A_{i,D_i}$ for any pair $(d,g) \in \mathcal{D} \times \mathcal{S}$; and 3) they can identify treatment effects even if all treated individuals with non-null adherence are associated with a non-null outcome value in a given group.\footnote{This is not the case for nonlinear binary choice models where identification of treatment effects requires the presence of at least one null (i.e., zero) outcome for each treatment option in a given group.} Moreover, information criteria (e.g., AIC and BIC) indicate better in-sample goodness-of-fit with LPMs compared to Probit binary choice models in most cases (results not shown). Observed adherence $A_{i,d} \in [0,1]$ was computed for each treated individual using the recorded DAA start and end dates and the standard treatment length for each treatment option (8 or 12 weeks). Complete adherence to the assigned treatment was also confirmed by the participating centers for each patient. I use 5-fold cross-validation with (at least) four different data splittings to validate the choice of the total number of groups $S$, the choice of the minimum eigenvalue $\lambda_{min}$, and the specification of each within-group outcome model. Both the weighted K-means algorithm and the estimation of all outcome models were performed using only the “training set" of observations for each fold and each data splitting. This leads to the estimation of OOS prediction errors over an OOS size (at least) four times the size of the sampled population. I use the area under the receiver operating characteristic (AUC-ROC) curve and the maximum sum of sensitivity and specificity obtained by the OOS predicted outcome values to assess the OOS accuracy of each combination of hyperparameter values and choice of outcome models. To facilitate the cross-validation procedure, I use the optimal assignment values obtained with the entire sample, $\hat{\mathbf{z}}(\hat{\boldsymbol{\mu}},\hat{\Sigma})$, as the initial assignment values for each training set. Finally, inference on the selected within-group CATEs is performed by increasing the number of total splits to 200 (which leads to 1,000 values for all CATEs) and by using quantile aggregation to generate median values and confidence intervals with 95% coverage, as shown in chernozhukov_fisher-schultz_2025.

Adherence Model

As mentioned in Section (ref), I model adherence with a two-part model where the first part models the probability of fully adhering to treatment, whereas the second part models the continuous adherence value if the observed adherence is below unity. For simplicity, both parts are modeled using a linear model with a relatively small set of not fully overlapping covariates for each part of the model. This namely includes : predicted treatment effectiveness under full adherence ($\hat{\tau}_d(1)$; first part only), clinical center fixed effects (first part only), the degree of hepatic cirrhosis (first part only), recent injection drug use (both parts), ethnicity indicators (both parts), and a fixed address indicator (both parts). The final choice of covariates in each part of the model is based on statistical significance and theoretical importance of the variables. To compute $\hat{\tau}_d(1)$, I use a subjective prior distribution where the group that minimizes the squared Mahalanobis distance $d_i^{2M}(\hat{\boldsymbol{\mu}}_g,\hat{\Sigma}_g)$ is assigned a weight of one for all $i \in [N]$. This choice of prior distribution is based on the idea that the estimated group memberships are converging to the true ones as $p$ tends to infinity according to Theorem (ref). The number of covariates included in the second part is much lower than in the first part since the majority of patients perfectly adhered to their prescribed treatment. I then incorporate predicted adherence values in each outcome model to yield the estimated CATEs with the predicted adherence $\hat{\tau}_{d,g}(\hat{a}(W_i,\hat{\tau}_{d}(1)))$ for every observation in the sampled population and every treatment option. However, I set $\hat{a}(W_i,\hat{\tau}_{d}(1)) = A_{i,d}$ if $d=D_i$ such that the adherence model features no prediction error when the selected treatment option corresponds to the observationally assigned one.

Feasible Policy Rule

Once I obtain the set of estimated group memberships $\hat{\mathbf{z}}(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ from the weighted K-means algorithm, I then use the rpart command from the rpart package in R to predict $\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ for each $i \in [N]$ in order to derive a feasible policy rule. The rpart command allows the econometrician to fix the maximum depth of the decision tree and the degree of complexity of the tree. For simplicity, I set the degree of complexity of the decision tree to the largest degree (i.e., less complex) such that each group membership is predicted by at least one terminal leaf. Missing/unknown values in the matrix of covariates $X$ are handled by surrogate rules following the recommendation of breiman_classification_2017. The rpart command also directly provides the conditional probability $\text{Pr}[\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma}) =g |L_i = k]$, where $L_i$ is the terminal leaf associated with the $i^{th}$ observation, and where $k \in \{1,2...,|L|\}$ with $|L|$ corresponding to the number of terminal leaves in the tree. Those probabilities are used to compute the Bayesian version of the feasible policy rule, as described in Section (ref).

Cost-Effectiveness Analysis

Cost-effectiveness analyses (CEAs) are then realized for both perfect and imperfect (i.e., predicted) adherence to treatment and for each “type" of estimated group memberships, which leads to four different CEAs. The goal of each CEA is to show how the optimal treatment allocation changes in the sampled population according to the policymaker's WTP. WTP is expressed in terms of costs (\$) per additional percentage point of probability of reaching SVR in each CEA. I define the incremental cost-effectiveness ratio (ICER) of treatment option $d$ for any individual $i$ with adherence level $a_i \in [0,1]$ as follows

where $c_d$ refers to the cost of prescribing treatment option $d$ to any individual in the sample, $\hat{s}_i = \{\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma}), \hat{h}(V^*_i)\}$ is an estimator of the true group membership $s_i$, and where the denominator is expressed in percentage points. Optimal treatment allocation is then obtained by selecting the treatment option that maximizes the chance of reaching SVR for every individual in the sample, and for which the corresponding ICER value remains below the policymaker's WTP. Formally, this can be represented by the following constrained maximization program

where the feasible policy rule is now treated as a function of the adherence vector $\mathbf{a} = (a_1,...,a_i,...,a_N)^{\top}$, and the policymaker's WTP, denoted by $\omega$, and where the $C$ subscript refers to the fact that $\hat{s}_i$ is treated as the “true" group membership in the above maximization program. It is then possible to illustrate how the in-sample optimal treatment allocation changes as a function of adherence and/or the policymaker's WTP. It is also possible to use the predicted adherence $\hat{a}_j(\cdot,\cdot)$ to identify the optimal treatment choice for a new patient $j$ and her corresponding set of covariates $V^*_j$ for a given value of $\omega$. It is also straightforward to adapt such a CEA to any decision criterion by modifying the maximand in Eq. ((ref)) accordingly, as shown in Section (ref). Finally, I also compare the total costs and predicted outcome values provided by any treatment allocation $\hat{\pi}(\mathbf{a},\omega)$ to a “benchmark" allocation. This comparison is performed through the computation of an aggregate ICER value, which is allowed to vary with adherence, WTP, and the choice of benchmark allocation. This aggregate ICER value is defined as follows

where $\hat{\pi}_i(\mathbf{a},\omega)$ stands as the treatment choice for individual $i$ at WTP $\omega$ and adherence vector $\mathbf{a}$, $\mathbf{b} = \{b_i\}_{i \in \mathcal{B}}$ is the set of benchmark treatment choices for all patients in the subset $\mathcal{B} \subseteq [N]$ with $|\mathcal{B}| = N_b$, and where $\mathbf{\tilde{a}} = (\tilde{a}_1,...,\tilde{a}_i,...,\tilde{a}_N)^{\top}$ corresponds to the set of adherence values used for benchmarking (which may be equal to $\mathbf{a}$). Note that $\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ is used as the true group membership for any $i \in [N]$. The value of $\text{ICER}_A(\mathbf{a},\mathbf{\tilde{a}},\mathbf{b},\omega)$ can then be plotted against WTP and/or adherence values for a given benchmark allocation $\mathbf{b}$ to identify the situations in which the treatment allocation $\hat{\pi}(\mathbf{a},\omega)$ dominates (or is dominated) by the benchmark at the aggregate level. Note that if $\text{ICER}_A(\mathbf{a},\mathbf{\tilde{a}},\mathbf{b},\omega)$ is positive, this implies that there is a trade-off between total health benefits and total costs, and neither $\hat{\pi}(\mathbf{a},\omega)$ nor the benchmark allocation dominates the other one. If $\text{ICER}_A^*(\mathbf{a},\mathbf{\tilde{a}},\mathbf{b},\omega)$ is negative, then one of the two allocations dominate the other at the aggregate level, which does not imply that one of the two allocations is pareto optimal. However, Pareto optimality can be achieved by enforcing an additional positivity constraint such that every term in the sum of the denominator of Eq. ((ref)) is positive for any $(\mathbf{a},\mathbf{\tilde{a}},\mathbf{b},\omega) \in [0,1]^{2N} \times \mathcal{D}^{N_b} \times \mathbb{R}_{>0}$.

Results

Descriptive Statistics

Table (ref) shows various descriptive statistics by treatment regimen in the sampled population. All the values presented in Table (ref) are mean values with their respective standard errors in parentheses. Line (1) indicates that around 80% of the patients prescribed Harvoni in the sample are men, compared to less than 70% for all other treatment regimens. Line (2) shows that never-treated patients are younger, on average, than all treated patients in the sample, although this difference is not statistically significant. Line (3) indicates that there exists a slight gradient in terms of liver stiffness across treatment regimen on average, where the average patient prescribed Harvoni has a higher liver stiffness than the one prescribed Epclusa, which has a higher liver stiffness than the one prescribed Mavyret.\footnote{A normal liver stiffness corresponds to a FibroScan measure that lies between 2 and 7 kilopascals (kPa). A FibroScan measure higher than 19 kPa is usually associated with liver cirrhosis.} Line (4) also shows that untreated patients have somewhat lower income levels, which can be explained by financial constraints preventing access to DAA treatments before the advent of public coverage by all Canadian provinces in 2018 (except Newfoundland). Typically, such differences in observed covariates likely indicate the presence of endogenous selection into treatment, but only unconditionally. Under Assumption (ref), conditioning on the unobserved group membership leads to plausible exogenous treatment assignments for every individual in the sample. Lines (5) and (6) of Table (ref) show the average adherence level to the corresponding regimen and the relative occurrence of SVR by treatment regimen, respectively. Line (6) indicates that SVR is achieved slightly more frequently under Harvoni compared to the two other DAA regimens, even if patients who are prescribed Harvoni are somewhat older and have worse liver conditions, on average, than the average patient associated with any other treatment option. Line (5) of Table (ref) also shows that adherence levels are very high for all treatment regimens. This is explained by the fact that adherence is perfect for more than 90% of all treated individuals for each treatment regimen, separately. Note that for the Mavyret treatment regimen, only one observation presented imperfect adherence with an adherence value virtually equal to zero (i.e., the patient was prescribed Mavyret but did not initiate the treatment). Note also that the first column of Table (ref) indicates a spontaneous clearance rate of 18% in the never-treated population.

Classification Procedure and Outcome Models

Table (ref) shows the results of the classification and regression procedures both in- and out-of-sample, with OOS referring to the results obtained with the 5-fold cross-validation procedure. The first line of Table (ref) presents the results of the regression procedure when there is only one group, hence without applying the weighted K-means algorithm. The first column of Table (ref) shows the value $\lambda^*_{min}$, which corresponds to the minimum eigenvalue $\lambda_{min}$ that is associated with the largest OOS maximum sum of specificity and sensitivity (column (6)).\footnote{Here, I consider that the OOS maximum sum of specificity and sensitivity is better than the OOS area under the ROC curve to identify the general performance of a binary classifier. This is because the former indicator relies exclusively on the optimal point in the ROC curve, which is the point that is closest to the top left corner of the ROC curve (i.e., the “perfect classifier"). In fact, the area under the ROC curve has often been criticized for considering regions of the ROC curve that are not relevant regarding the predictive accuracy of a binary classifier halligan_disadvantages_2015,janssens_reflection_2020.} The second column of Table (ref) shows the value of global minimum of the TSMD, $d^{2TM}(\hat{\mathbf{z}}(\hat{\boldsymbol{\mu}},\hat{\Sigma}),\hat{\boldsymbol{\mu}},\hat{\Sigma})$, obtained by the weighted K-means algorithm for each value of $S \in \{2,...,6\}$, where $\hat{\mathbf{z}}(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ corresponds to the group memberships associated with the global minimizers $\hat{\boldsymbol{\mu}}$ and $\hat{\Sigma}$ identified by the algorithm. To note, the total distance $d^{2TM}(\hat{\mathbf{z}}(\hat{\boldsymbol{\mu}},\hat{\Sigma}),\hat{\boldsymbol{\mu}},\hat{\Sigma})$ can be negative if $\lambda_{min}$ is between zero and one due to the presence of the term $\log \det(\Sigma_g)$ in the objective function. Columns (3) to (6) in Table (ref) show that the maximum OOS sum of specificity and sensitivity is obtained when $S=5$ and $\lambda_{min} = 0.70$. The OOS area under the ROC curve associated with this combination of hyperparameters is also large (0.8196), but not the largest among all results. Taken together, I consider that $S=5$ and $\lambda_{min} = 0.70$ identify group memberships and outcome models with the best OOS predictive accuracy among all results. A more thorough search across the hyperparameter space could eventually yield marginally better hyperparameters, but at the cost of greater computational overhead. Note also that small values of $\lambda_{min}$ tend to lead to high levels of predictive accuracy in-sample, but to lower predictive accuracy out-of-sample (e.g., when $S=3$ in Table (ref)). This is mainly caused by in-sample overfitting, where small values of $\lambda_{min}$ assign very large weights to a small number of covariates during the classification procedure. The results shown in columns (3) to (6) of Table (ref) are all based on the following LPM

where $\beta_{0,g}$ is the intercept value for the $g^{th}$ group, $\tau_{d,g}$ is the CATE of treatment option $d$ for the $g^{th}$ group, and where $\epsilon_{i} \sim N(0,\sigma^2_{\epsilon})$ with $\sigma^2_{\epsilon} > 0$. Those simple models do not include any of the covariates $X_i$ used for classification due to systematically lower OOS performances compared to the performances of simple LPMs with no additional covariates (results not shown). The absence of necessity to control for additional covariates in the LPMs is also a sign that the estimated group memberships with $S=5$ appropriately account for most of the unobserved heterogeneity across patients. The intercept value for the $g^{th}$ group, $\beta_{0,g}$, denotes the probability of observing spontaneous clearance for any untreated individual in the $g^{th}$ group.

The estimated coefficients of each outcome model for the best combination of hyperparameters (i.e., $S=5$ and $\lambda_{min}= 0.7$) and the full sample are presented in Figure (ref). These coefficients are obtained by taking the median over the 1,000 values obtained from the sample splitting procedure, explaining why some bars are not centered within the confidence interval. The first set of bars on the left indicates that individuals belonging to group 1 have a probability of around 80% to spontaneously clear the HepC virus. This result could be leveraged to reduce total prescription costs in the future if there are no substantial costs of waiting to confirm spontaneous clearance for individuals within this group in the target population. This might not be the case in practice, as there are significant costs to society associated with not treating a patient who may spread the infection while waiting. The descriptive statistics within each group (shown in Appendix (ref)) also show that patients in Group 1 tend to be younger and have a higher income than patients in all other groups. However, these differences are not statistically significant between groups. Results from Figure (ref) also show that Epclusa is weakly dominated when there is no cost constraint due to the lower or equal effectiveness of Epclusa in all groups compared to the other DAA regimens. Those results also imply that 58.6% of all treated individuals in the sample would have observed an equal or higher probability of reaching SVR if their prescription had been switched to another treatment option, as shown in Table (ref).

The last set of bars in Figure (ref) shows the estimated coefficients of the LPM when the full sampled population is not divided into different groups. This last set of results produces an effectiveness ranking that is slightly different from the ranking portrayed in Table (ref) due to imperfect adherence in the LPMs. Harvoni remains the most effective option, but is followed by Epclusa and then Mavyret. This global effectiveness ranking is preserved in all groups, except in Group 1, where the ranking is completely reversed, and in Group 5, where Harvoni and Epclusa share the same median effectiveness estimate. Note, however, that no DAA treatment regimen (apart from the no-treatment regimen) is significantly more effective than any other regimen within each group. Although this precludes us from concluding that a given DAA regimen is truly more effective than the two other ones for any given patient belonging to a given group, such results can still be used to guide future decisions. If there is no constraint imposed on treatment choice and estimates of marginal effects are unbiased, then it is always better to prescribe a superior treatment in terms of the conditional probability of reaching a desirable outcome. For comparison purposes, Figure (ref) shows the estimated marginal effects on SVR of each treatment option when the hyperparameter values yield the best in-sample predictive performance (i.e., $G=3$ and $\lambda_{min}=0.04$). The most remarkable difference between Figure (ref) and Figure (ref) lies in the first group, where reducing the number of groups to $G=3$ now leads to negative and statistically significant marginal effects of Epclusa and Harvoni on SVR, while Mavyret does not feature any significant effect on SVR. Although the results shown in (ref) are more parsimonious and feature better in-sample performances than those shown in Figure (ref), using these results to guide prescription behaviors in the general population would be misguided if the OOS maximum sum of specificity and sensitivity is a better indicator of external validity. This example also illustrates how misspecification of the number of groups and in-sample overfitting can lead to biased estimates and subsequent suboptimal choices.

For further comparison, I use the multi_arm_causal_forest command from the grf package to generate estimates of CATEs in a multi-action PL setting using causal random forests, which is one of the most popular approaches for policy learning in the econometric literature tibshirani_generalized_2024,wager_estimation_2018,athey_policy_2021. However, two elements limit the comparison between the proposed approach and the one using causal random forests. First, the multi_arm_causal_forest command does not handle semi-continuous treatment effect indicators (due to imperfect adherence). To bypass this issue, I enforced $A_{i,D_i} = 1$ if and only if the observed adherence for assigned treatment $D_i$ was higher than 0.5. Second, the same command does not provide the predicted outcome values for the left-out observations that were not used for training the different trees. Consequently, it is impossible to compare the indicators of OOS predictive accuracy between the two methods. Despite these two caveats, the grf package yielded an AUC of the in-sample ROC curve of 0.7383, with an in-sample maximum sum of specificity and sensitivity of 1.6395, both much lower than what is shown in columns (3) and (5) of Table (ref). Given that indicators of OOS performance are usually lower than in-sample indicators, likely, the multi_arm_causal_forest command would not perform as well as the proposed approach in the target population.

Adherence to Treatment

Table (ref) shows the estimated coefficients for each part of the adherence model described in Section (ref) when applied to treated individuals only. The predicted effectiveness of the Epclusa treatment regimen was removed from the binary part of the model due to the very high collinearity with the predicted effectiveness of the Harvoni treatment regimen. The absence of strong collinearity among the selected variables was assessed using the variance inflation factor (VIF $<$ 5). The first column of Table (ref) shows that predicted treatment effectiveness of Mavryet, liver stiffness, and quality of life are significantly and positively associated with the probability of completely adhering to treatment.\footnote{The statistically significant coefficient on predicted treatment effectiveness of Mavyret is because the only patient who did not completely adhere to treatment features a relatively low predicted treatment effectiveness.} This is also the case for people who have a fixed address (+14.1 percentage points), who are either Asian (+7.6 percentage points) or Metis (+7.2 percentage points), and who are living in a long-term care facility (+4.8 percentage points). The fixed effects associated with clinical centers 3, 4, 11, 16, and 19 are also positive and statistically significant at the 95% confidence level (results not shown). Such results warrant further investigation to better understand why patients treated in those clinical centers are more likely to completely adhere to treatment. On the other hand, living alone or sharing a room/apartment is negatively and significantly associated with the probability of fully adhering to treatment. Apart from the intercept, being a part-time worker is the only variable significantly associated with a lower level of adherence among observations with imperfect adherence (63 observations only). The other variables are included in the continuous part of the model since their respective p-value is lower than 0.3 and also because those variables are suspected to affect the level of adherence according to the literature on DAAs harney_unsuccessful_2025, darvishian_loss_2020, wilton_real-world_2020, jiang_hcv_2023. Note that the signs of the estimated coefficients align with what is expected according to this same literature. Note also that the Asian and Metis ethnicity binary indicators were not included in the continuous part due to a lack of variation among the 63 observations with imperfect adherence.

Feasible Policy Rule and Cost-Effectiveness Analysis

The structure of the decision tree used to predict group memberships is presented in the appendix and was provided directly by the rpart.plot command of the rpart package in R. The selected decision tree contains 20 terminal leaves and has a depth of 10 (shown in Appendix (ref)). This tree shows that group memberships can be predicted using the results of various blood tests (e.g., blood levels of insulin, glucose, haemoglobin, total bilirubin, aspartate transferase (AST), gamma-glutamyl transferase (GGT), and some others), the degree of liver stiffness (based on FibroScan measures), whether if the patient is Indigenous, the degree of pain reported by the patient, and the patient's height.

Table (ref) shows the treatment choice for each terminal leaf of the selected decision tree based on the Bayesian version of the feasible policy rule when there is no cost or global constraint, as shown in Eq. ((ref)). Treatment choices presented in Table (ref) are computed using the conditional probabilities $\text{Pr}[\hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma}) =g |L_i = k]$ provided by the rpart function, as shown in Section (ref). Table (ref) shows that Mavyret and Epclusa are the optimal choices in one (no. 13) and two (no. 6 and no. 15) terminal leaves out of 20, respectively. Moreover, the Epclusa treatment regimen is the minimax-regret optimal treatment regimen based on the results shown in Figure (ref). This means that the Epclusa treatment regimen performs relatively well across all groups, with a maximum regret (or welfare loss) of 10.4 percentage points -- compared to 11.4 and 16.2 percentage points for Harvoni and Mavyret, respectively --, which occurs when the patient truly belongs to Group 1 and Epclusa is prescribed instead of Mavyret. However, minimax-regret optimality depends on the effectiveness of every considered alternative, implying that Epclusa may not remain minimax-regret optimal if other treatment regimens had been added in the analysis.

Figure (ref) shows the results of the CEA described in Section (ref) under perfect adherence for every treatment option and every individual in the sample when $\hat{s}_i = \hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma})$. Panel (a) of Figure (ref) shows the total health benefits on the left axis, represented by the sum of the marginal effects (in total percentage points, $pp$) across every individual's optimal treatment choice, and the corresponding total prescription costs on the right axis (in CAN\$M), both as a function of the payer's WTP. Panel (b) of Figure (ref) shows the corresponding optimal treatment allocation as the percentage of individuals in the sample that receive each treatment option as a function of WTP. Finally, panel (c) of Figure (ref) shows the value of the aggregate ICER when the benchmark corresponds to the observationally assigned treatment $D_i$ (the status quo allocation), as described by Eq. ((ref)). Note that Panel (c) of Figure (ref) is based exclusively on treated patients given the nature of the benchmark. Figure (ref) shows that most of the health benefits are generated when the WTP reaches CAN\$580 per additional percentage point. Specifically, 87.4% of the maximum health benefits are generated when the WTP is equal to CAN\$580/pp, with total costs corresponding to 56.4% of the maximum total prescription costs (CAN\$44.4 million vs CAN\$78.9 million). Panel (b) of Figure (ref) also shows that 90.2% of all individuals in the sample are optimally treated with Mavyret when WTP is between CAN\$580/pp and CAN\$650/pp. Panel (b) also shows that all 121 individuals belonging to Group 1 are left untreated when the WTP is between CAN\$580/pp and CAN\$2,010/pp. When WTP goes beyond CAN\$2,010/pp, then optimal allocation implies that 85.3%, 9.8%, and 4.9% of all individuals in the sample are prescribed Harvoni, Mavyret, and Epclusa, respectively, which leads to the maximum health benefits achievable with perfect adherence. Panel (c) of (ref) shows that the corresponding optimal treatment allocation leads to reduced costs (savings of CAN\$2.9 million) and higher health benefits (+21.3 pp) compared to the status quo when WTP lies between CAN\$900/pp and CAN\$2,010/pp. Here, the negative values of the aggregate ICER correspond to negative total costs (i.e., savings) divided by positive total health benefits, therefore corresponding to an allocation that dominates the status quo at the aggregate level. To give a better sense of the magnitude of the health benefits compared to the cost of treatment, it is possible to compare the total cost of a given treatment option per expected gains in quality-adjusted life years (QALYs). For instance, if we assume that Mavyret features a marginal effect of 75 percentage points on the probability of reaching SVR, and that being cured from HCV leads to an increase in quality of life between 1.3 and 2 QALYs mattingly_changes_2020, such a hypothetical scenario translates into a total cost that lies between CAN\$28,500 and CAN\$43,850 per additional QALY given the total cost of the Mavyret treatment regimen (CAN\$40,000.00 for 8 weeks). The upper bound is below the standard conservative threshold of CAN\$50,000 per QALY. In contrast, the list price of Harvoni is CAN\$67,000.00 for 12 weeks. This leads to a total cost that lies between CAN\$41,875 and CAN\$64,425 per additional QALY if we assume a marginal effect of 80 percentage points on the probability of reaching SVR, which is still relatively low.

Panel (a) of Figure (ref) shows the optimal treatment allocations when the tree-based predicted group memberships $\hat{h}(V^*)$ are used instead of the estimated group memberships $\hat{\mathbf{z}}(\hat{\boldsymbol{\mu}},\hat{\Sigma})$. Using the tree-based predicted group memberships leads to a misclassification rate $\text{Mis}(\hat{h}(V^*),V^*)$ of 23.1%. It is possible to quantify the impact of these classification errors by computing an aggregate ICER value using Eq. ((ref)) where $\hat{\pi}_i(\mathbf{a},\omega)$ corresponds to the treatment allocation shown just above in Panel (a), adherence values $\mathbf{a}$ and $\tilde{\mathbf{a}}$ are all set to one, and where the benchmark allocation $b_i$ is set to either the status quo (Panel (b)), or to the optimal treatment allocation obtained when $\hat{s}_i = \hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ for every $i\in [N]$ in Eq. ((ref)) (Panel (c)). Panel (b) of Figure (ref) shows that using $\hat{h}(V^*)$ to compute treatment choices eliminates the dominant cases where lower costs and higher health benefits would occur compared to the status quo. In fact, when WTP is between CAN\$840/pp and CAN\$890/pp, the status quo allocation provides more health benefits at lower costs (on aggregate) than the optimal treatment allocation based on $\hat{h}(V^*)$ due to misclassification errors. All other WTP values imply trade-offs between costs and health benefits. On the other hand, Panel (c) of Figure (ref) shows that, when the benchmark corresponds to the allocation obtained when $\hat{s}_i = \hat{z}_i(\hat{\boldsymbol{\mu}},\hat{\Sigma})$, treatment choices based on $\hat{h}(V^*)$ are dominated for every value of WTP above CAN\$2,020/pp. Note that the total health benefits and costs associated with this CEA are not shown for brevity.

Figure (ref) shows the results of the CEA with predicted adherence and estimated group memberships. All curves in each panel of Figure (ref) are now “smooth" since predicted adherence is a continuous variable that lies between 0 and 1 and each estimated CATE with predicted adherence $\hat{\tau}_{d,g}(\hat{a}(W_i,\hat{\tau}_d(1)))$ is obtained by multiplying $\hat{\tau}_{d,g}(1)$ with the predicted adherence obtained from the two-part model. The upper panel of (ref) shows the optimal treatment allocation obtained when predicted adherence is used instead of perfect adherence, which leads to slight increases in the share of patients optimally treated by Mavyret and Epclusa when WTP is high ($>$ CAN\$2,600/pp). Panel (b) of Figure (ref) shows that the status quo allocation under predicted adherence is dominated by the allocation shown in Panel (a) for every value of WTP that is between CAN\$900/pp and CAN\$2,010/pp, where the probability of reaching SVR increases by a total of 460.8 pp to 841.4 pp while total costs are reduced between CAN\$3.6M and CAN\$4.9M. The allocation shown in Panel (a) also dominates the status quo when WTP is between CAN\$2,010/pp and CAN\$2,080/pp, but with much smaller savings ($<$ CAN\$0.6M). Because using predicted adherence might change the optimal treatment allocation -- and therefore change total health benefits and costs -- it is possible to perform the CEA with optimal treatment allocation under perfect adherence as the benchmark. This is what is shown in Panel (c) of Figure (ref). Given that predicted adherence is always equal to or lower than perfect adherence, and that lower adherence always leads to lower health benefits, the aggregate ICER shown in Panel (c) of Figure (ref) is always positive. This implies that using predicted adherence instead of perfect adherence always implies a trade-off where savings are realized, but at the expense of lower health benefits (or everyone remains in the null treatment regimen). Nonetheless, the CEA with predicted adherence is believed to yield more credible results since observed adherence is indeed imperfect for a non-negligible number of individuals in the sample. Finally, Figure (ref) shows the result of the CEA with both predicted adherence and tree-based predicted group memberships, which is the analog of Figure (ref) but with predicted adherence. As in Figure (ref), the aggregate ICER associated with the treatment allocation shown in Panel (a) when the status quo allocation is used as the benchmark is shown in Panel (b) of Figure (ref). This comparative exercise shows that using the tree-based predicted group memberships with predicted adherence always leads to trade-offs between total costs and total health benefits when compared to the status quo. When WTP is below CAN\$900/pp, both total costs and health benefits are negative compared to the status quo allocation (with predicted adherence), whereas total costs and health benefits become both positive if WTP $>$ CAN\$900/pp. Finally, Panel (c) of Figure (ref) shows that the treatment allocation of Panel (a) is dominated by the treatment allocation obtained from the estimated group memberships under perfect adherence when WTP lies between CAN\$730/pp and CAN\$830/pp, and also when WTP $>$ CAN\$880/pp. These results show that, although using the tree-based predicted group memberships can be quite practical, it increases the probability of suboptimal treatment allocations compared to using the estimated group memberships, especially when imperfect adherence is taken into account.

Discussion and Conclusion

Decision-making processes are pervasive in daily life. In this paper, I developed a relatively simple methodology to help policymakers build policy rules that can guide decision-makers using observational data. Commonly known policy rules are often not mandatory; this can take the form of guidelines, evaluation grids, or even personal criteria for daily consumption choices. Even if they are not strictly enforced, policy rules remain important in practice if they serve as a baseline decision rule upon which agents can improve based on contextual information and personal judgment. In health care and health sciences in general, clinical judgment of health practitioners is often opposed to “evidence-based medicine", where the former is assumed to be often unreliable by the proponents of the latter, whereas the latter is often criticized for not taking into account the contextual elements of any single choice problem alpert_why_2010,knaapen_evidence-based_2014. As mentioned recently by several authors, those two visions are not dichotomous but should rather coexist for the benefit of the patients riviere_can_2020, subbiah_next_2023, einav_limitations_2024. This includes leveraging observational data to better identify the real-world elements that influence treatment effectiveness at the population level, and then using those (macro) elements to guide the decision-making process at the patient's (micro) level without abstracting away from contextual information. Better clinical guidelines that are based on real-world evidence and that acknowledge the usefulness of discretionary decisions under specific circumstances are therefore warranted and necessary to improve the quality of care and to reduce the incidence of medical errors in practice. The case of DAA treatment choice for HIV/HCV co-infected patients is one of the many examples where observational data can be used to inform policymakers about the trade-offs that exist in the real world when formulating policy rules, especially when those policy rules can be used in a legal context to defend the various choices made by the implied decision-makers. The methodology developed in this paper to design (optimal) policy rules has several benefits : it relies more practical assumptions than what is typically assume in the literature in policy learning, it advocates for the use of a specific clustering algorithm that leads to the consistent estimation of all CATEs in the sample, it easily deals with missing values in covariates, it can handle semi-continuous treatment variable in multi-action framework, it is more parsimonious than typical ML methods -- which leads to more precise estimates of all CATEs in the sample while reducing the probability of overfitting at the same time --, and it yields interpretable results that can be easily used to formulate updated guideline recommendations, as shown in the paper. Nonetheless, the proposed methodology has several limitations. First, the credibility and usefulness of the estimated policy rules will depend mostly on whether group-wise unconfoundedness and group-wise heterogeneity (Assumptions (ref) and (ref)) are both satisfied in practice. If they are not, the estimated policy rules might lead to suboptimal treatment allocations, which may then lead to worse health outcomes and/or higher costs than the status quo. Making the proposed estimation methodology robust to violations of these two assumptions is left to future research. Second, the computational burden of the proposed approach is relatively high, especially when a minimum eigenvalue $\lambda_{min}$ has to be enforced for the weighted K-means algorithm to yield reliable results. To reduce the computational burden, it might be a good idea to exclude all discrete variables from the set of covariates used to perform the weighted K-means algorithm. By doing so, all within-group estimated variances will be bounded away from zero, unless groups can be formed by a single observation. Discarding initial values that lead to such a degenerate result thus allows for the minimum eigenvalue $\lambda_{min}$ to be very close to zero and greatly alleviates the computation burden. Third, it requires enough within-group variation in order to identify all CATEs. This is, however, an issue with all ML techniques used to identify treatment effect heterogeneity, and the fact that the method is more parsimonious than most other ML techniques makes this issue less relevant in practice. Fourth and last, the method cannot identify contraindications that are being strictly followed by most prescribing physicians in their current practice. Although this limitation is not related to the proposed estimation strategy per se, it is important to keep in mind that contraindications that are well known by physicians are likely to be absent from observational datasets. This may then lead to estimated “optimal" treatment allocations that do not account for these contraindications, which may worsen health outcomes compared to the status quo. Results of the empirical application show that it is possible to improve upon the current practice when treating co-infected HIV/HCV patients in Canada. Specifically, estimation results show that it would have been possible to reduce prescription costs in the CCC by an amount lying between CAN\$3.6 million and CAN\$4.9 million while still increasing the total health benefits just by reallocating treatment options among treated individuals in the cohort. However, such important savings are mostly due to the existence of Group 1, where it is estimated that 121 individuals in the sample have approximately an 80% chance of spontaneously clearing HCV without taking any treatment. Because acute infection is temporary but still lasts around 6 months on average, not treating these patients may not be optimal if this leads to an important rise in the reinfection rate within the target population young_rate_2023. Results from the empirical application also show that treatment reallocation could have substantially increased the probability of reaching SVR for every treated individual in the sample (+2,849.2 pp among 745 treated individuals under predicted adherence), but this would have come at the expense of slightly larger prescription costs (+CAN\$0.57 million). Although this paper focused exclusively on SVR, other outcomes are also important for both mono-infected HCV and co-infected HIV/HCV patients: probability of reinfection, degree of liver stiffness, mortality, etc. Still, it is known that being cured of HCV improves patients' quality of life and other nonvirological outcomes simmons_long-term_2015, park_impact_2021. This is why SVR remains the main outcome when assessing the effectiveness of DAA regimens, although other outcomes should be included in future analyses. At last, it is important to recall that, according to the obtained results in the context of HIV/HCV co-infection, no DAA regimen is dominated by any other single DAA regimen or set of DAA regimens when the policymaker's WTP is below a certain level. In other words, the optimal treatment choice for any HIV/HCV co-infected patient is not trivial and depends on the characteristics of the patient, on the chosen decision criterion, and the policymaker's WTP. Such a result warrants caution when physicians prescribe an anti-HCV treatment for a never-treated, co-infected patient, especially when costs vary greatly across treatment options. Another interesting feature of the result is that optimal treatment choice for the physician does not seem to depend on the ART taken by the patient, although this could be a consequence of contraindications in the prescribing patterns of physicians regarding HIV/HCV co-infection, as mentioned above. More research is required to validate such an intuition.