{'fraction_non_alphanumeric': 0.03349316159167714, 'fraction_numerical': 0.01440223641607247, 'mean_word_length': 4.823802163833076, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Using Multiple Imputation to Classify Potential Outcomes Subgroups

With medical tests becoming increasingly available, concerns about over-testing and over-treatment dramatically increase. Hence, it is important to understand the influence of testing on treatment selection in general practice. Most statistical methods focus on average effects of testing on treatment decisions. However, this may be ill-advised, particularly for patient subgroups that tend not to benefit from such tests. Furthermore, missing data are common, representing large and often unaddressed threats to the validity of statistical methods. Finally, it is desirable to conduct analyses that can be interpreted causally. We propose to classify patients into four potential outcomes subgroups, defined by whether or not a patient's treatment selection is changed by the test result and by the direction of how the test result changes treatment selection. This subgroup classification naturally captures the differential influence of medical testing on treatment selections for different patients, which can suggest targets to improve the utilization of medical tests. We can then examine patient characteristics associated with patient potential outcomes subgroup memberships. We used multiple imputation methods to simultaneously impute the missing potential outcomes as well as regular missing values. This approach can also provide estimates of many traditional causal quantities. We find that explicitly incorporating causal inference assumptions into the multiple imputation process can improve the precision for some causal estimates of interest. We also find that bias can occur when the potential outcomes conditional independence assumption is violated; sensitivity analyses are proposed to assess the impact of this violation. We applied the proposed methods to examine the influence of 21-gene assay, the most commonly used genomic test, on chemotherapy selection among breast cancer patients.


Introduction
Newly diagnosed cancer patients and their physicians frequently face the challenges of deciding which tests and subsequent treatments to use. The testing itself can be both invasive and expensive, and the cost of medical tests reaches 70 billion dollars annually 1 . Furthermore, test results do not predictably inform treatment decisions, either because guidelines are not always available or they are not followed in practice. Thus, describing how testing impacts the choice of treatment in general practice is of fundamental importance. We will view this problem in a causal framework, in which taking the test is the exposure of interest, and the treatment received is the outcome of interest.
An example is the use of the 21-gene recurrence score (RS) assay in breast cancer care, the most commonly used genomic assay in the United States. The test is used to inform recommendations about the use of chemotherapy. Women who get a high score on the test are recommended to take chemotherapy to reduce their risk of recurrence, whereas women with a low score are not because they have a good prognosis and are expected to receive little additional benefit from chemotherapy. Prior studies have demonstrated that RS tests can provide more precise estimation of the risk of cancer recurrence and the need for chemotherapy than do traditional histopathologic features (e.g., tumor grade and size) [2][3][4][5][6][7][8] . Even though the current clinical guidelines only endorse the use of RS assay in selected patients with certain clinical characteristics 9,10 , more and more patients who do not meet these criteria are being selected for RS testing. Also, chemotherapy is highly morbid, and many believe it is being over-prescribed 11,12 . Thus, there is a pressing need to measure the influence of the RS assay on chemotherapy decisions to evaluate whether the test is being used in an efficient manner and acted upon appropriately in general practice.
Determining the influence of testing on treatment selection in a community practice setting is not well studied and is challenging to quantify using standard methods. First, selection bias may occur because testing is not offered randomly. Both clinical factors as well as patient and physician preferences influence whether testing is offered, and this confounds the effect of testing on treatment selection 13,14 . Second, while causal inference methods such as propensity-score based methods can address selection bias and give effect estimates with desirable causal interpretations, they are limited to estimating an overall effect of testing on treatment selection, not how testing changes treatment plans within patient subgroups. The overall effect can obscure the degree to which treatment plans are changed by testing. For example, the overall effect can be zero, if testing encourages and discourages treatment for an equal number of patients, even though testing did change treatment plans. Third, missing data are common, especially when multiple data sources are combined as in many studies. Missing data are often-unaddressed threats to the validity of statistical methods.
In this paper, we propose to use a potential outcomes framework (i.e., a causal inference framework) 33 to measure the influence of testing by comparing treatment use when everyone is tested with when no one is tested. This minimizes selection bias affecting test use that is attributable to observed confounders. Specifically, we propose to classify patients into four potential outcomes subgroups defined by their expected treatment decisions with a given test versus without the test: 1) Never-Treated: those who would not have been treated with or without the test; 2) Treatment-Encouraged: those who would have been treated with the test but would not have been treated without the test; 3) Treatment-Discouraged: those who would not have been treated with the test but would have been treated without the test; 4) Always-Treated: those who would have been treated with or without the test. These four subgroups (listed in Table 1) are called potential outcomes subgroups 33 . A patients membership in these subgroups may be driven by patient characteristics or by systemic factors (e.g., treating physician). Such classification can also suggest ways to improve the utilization of medical tests in treatment decisions. The classification recognizes that a medical test can have different influences on treatment decisions for different patients; i.e., there may be effect heterogeneity of testing on treatment. With the classification, we will be able to measure how many patients treatment selections are changed (or unchanged) by test results for overall and key clinical subgroups, and also quantify how much test results encourage (or discourage) treatment selections. The potential outcomes subgroups bear some resemblance to principal stratification 16 , which has been extensively used by us and others to study the causal effect of intermediate variables (e.g., degree of compliance, surrogate markers) [17][18][19] . However, unlike principal strata, which are defined by the joint distribution of intermediate variables under different exposure status, here we focus on the joint distribution of the outcome of interest (here treatment received) under different exposure status (here test received). We believe this framework provides important conceptual and practical advances in studying the influence of medical testing in treatment selection. We will also demonstrate its advantages of providing unique and highly relevant causal quantities of interest over alternative causal inference methods.
We will classify patients into these potential outcomes subgroups by exploring the close connection between missing data analysis and causal inference. Since two potential outcomes of the same patients are not simultaneously observed, we will use multiple imputation methods to impute the unobserved potential outcomes to enable us to make causal inferences. Multiple imputation methods have been used to impute missing potential outcomes in solving different problems in causal inference literature [22][23][24][25]27,28,30 . These methods have the advantage of simultaneously handling standard missing data in observable variables (e.g. due to non-response) and missing potential outcomes. However, since missing data analysis and causal inference have different objectives and data structures, there are important distinctions between them. First, causal inference focuses on comparing the potential outcomes of the same patient under different exposure levels; while missing data analyses do not. Second, the causal assumptions (e.g., ignorability of treatment assignment or potential outcomes conditional independence) are not usually made in standard missing data analyses. Third, causal inference distinguishes pre-exposure variables from post-exposure variables and usually does not use a post-exposure variable, while missing data analyses do not differentiate them and often utilize a post-exposure variable in the imputation models if it is predictive of other variables with missing values. Fourth, the missing potential outcomes problem is highly structured such that only one of two potential outcomes is observed for any patient when the exposure variable is binary. In this situation, the missing data indicator is identical to the binary exposure variable. Fifth, since no two potential outcomes of the same patient can be observed simultaneously, there is no information about the correlation between two potential outcomes of the same patient. This distinct data structure is similar to the missing data structure in the literature on "statistical matching" 20,21 . More detailed discussions about missing data analyses and causal inference can be found by Ding and Li (2018) and   29,30 . In this paper, we will study the performance of our proposed methods. We will also examine the efficiency gain of our causal estimates when we incorporate the knowledge of our data structure and causal inference assumptions into the multiple imputation process. We will propose sensitivity analyses to investigate the violation of one of the key causal inference assumptions on causal estimates (i.e., potential outcomes conditional independence). Additionally, we will explore the impact of incorporating a post-exposure variable in the multiple imputation process on the causal estimates. Furthermore, we will demonstrate that this approach is more robust towards certain types of model mis-specification than other causal methods (e.g., propensity-score based methods).
To illustrate our methods, we will examine the impact of the RS assay on the decision to undergo chemotherapy for breast cancer patients. The analyses will use combined data sets from the Surveillance, Epidemiology, and End Results databases (SEER) registry, genomic testing laboratories, and patient/physician surveys, recently collected from the ICanCare Study 31,32 . The data linkage was a unique survey-registry-industry collaboration and the ICanCare study was part of the only program project funded by National Cancer Institute that focused on cancer treatment decision making.

Notation, the Causal Framework and Quantities of Interest
Let Y represent the treatment receipt (i.e., the outcome of interest, 0/1 for no/yes), Z the status of being tested (i.e., the exposure of interest, 0/1 for no/yes), X a set of covariates and R the test result. Naturally, R only exists for testers and is a post-exposure variable. We assume some observations are missing in Y , R and in some variables of X. Let M be the associated missing-data indicator matrix for these observable variables. The observed and missing parts of Y are denoted by Y obs and Y mis respectively, and the same for X.
For each patient, we assume there are two potential outcomes, Y 1 and Y 0 , corresponding to being tested or not (Z = 1, 0). Let Y 1 indicate a patients treatment status that would be observed if the patient is tested and Y 0 otherwise. If we know both Y 1 and Y 0 for every patient, we can infer whether a test causes a change in treatment selection by comparing Y 1 to Y 0 , and calculate many causal quantities that are free of selection bias. In this manuscript, we choose several representative causal quantities including: (1) average testing effect (ATE), E(Y 1 − Y 0 ), for the overall population and patient subgroups; (2) marginal odds ratio (MOR), [P 1 (1 − P 0 )]/[P 0 (1 − P 1 )] where P 0 = P (Y 0 = 1) and P 1 = P (Y 1 = 1); and (3) the fractions of patients in each of the four potential outcome subgroups: (Y 0 = 0, Y 1 = 0) for the never-treated, (Y 0 = 0, Y 1 = 1) for the treatment-encouraged, (Y 0 = 1, Y 1 = 0) for the treatment-discouraged, and (Y 0 = 1, Y 1 = 1) for the always-treated group, which are denoted by P (00), P (01), P (10) and P (11) respectively. The MOR defined in (2) compares the odds of treatment use if everyone were tested to the odds if no one were tested. Odds ratios in observational studies are usually association measures and lack causal interpretations because of confounding. In contrast, MORs do have causal interpretations. Based on the fractions defined in (3), we can calculate many quantities that are functions of them. For example, we can calculate the proportion of patients whose treatment plans are changed by the use of testing (i.e., the proportion of patients who were either treatment encouraged or discouraged), P (01) + P (10), and the proportion of patients whose treatment plans are not changed (i.e., the proportion of patients who were either nevertreated or always-treated), P (00) + P (11). In fact, ATE is also a function of these fractions, which is equal to the difference in the proportions in the "treatment encouraged" and "treatment discouraged" subgroups, i.e., P (01) − P (10). We are also interested in estimating the associations between patient characteristics and patient membership in these potential outcomes subgroups.

Estimating the Distributions of Potential Treatment Selection
To measure the causal influence of a test, ideally, we wish to know the treatment plan for each patient when she is tested and when she is not. That is, we wish to know both Y 0 and Y 1 . However, the treatment plan is observable for a given patient in only one of these two scenarios because she may not get tested and; even if she does, her pre-test treatment plan is not usually documented. Formally, this is a missing data problem. We use sequential regression multiple imputation (or sometimes called the fully conditional specification method) 34,35 to impute the unobserved (counterfactual) treatment plan for each patient, which we term the potential outcomes multiple imputation method (POMI). That is, for patients who received the test, we impute what would have been their treatment plans had they not received the test. And for patients who did not receive the test, we impute what would have been their treatment plans had they received the test from the predictive distribution of Y given X and/or R. Once we have complete information for Y 0 and Y 1 after imputation, we can calculate the causal quantities of interest. Since we have incomplete information in observable variables Y , X and R, we can simultaneously impute missing potential outcomes and missing values in these observable variables during the multiple imputation process using the same imputation methods.
Although we can view causal inference as a missing data problem, there are distinct features of causal inference different from standard missing data problems. As stated in the Introduction section, the objectives of causal inference are to obtain the causal effect of an exposure on the outcome, while missing data analyses are to handle missing values in variables, utilize incomplete information and reduce bias due to non-response. Hence, some assumptions are closely related but others are different. Below we list four assumptions used in missing data and/or causal inference literature: 1. Missing at random 36 , f (M |Z, Y, X, R) = f (M |Z, Y obs , X obs , R obs ). That is, missingness can be explained by the observed data. This is a standard assumption in the missing data literature. Note that f () denotes a probability density or mass function for the distribution of either a continuous or discrete variable in this manuscript.
2. Stable unit treatment value assumption (SUTVA) 22,37 . The potential outcomes of one patient are unaffected by the potential outcomes or testing status of other patients. This is a standard assumption in the causal inference literature. This assumption is not usually explicitly made, but it is implied in the missing data literature.
X, that is, there are no unmeasured confounders between Z and Y after adjustment for X. It is a standard assumption in the causal inference literature, and a special case of the missing-at-random assumption (see Assumption 1). 4. Potential outcomes conditional independence, (Y 0 ⊥Y 1 )|X. This assumption states that two potential outcomes are conditionally independent. However, since the two potential outcomes are never observed together, we have no information about the correlation, independence or joint distribution of Y 0 and Y 1 from the observed data. This assumption is frequently made in the causal inference literature [39][40][41] . It is also made in the missing data and survey methodology literature when applying the less-known statistical matching technique to conduct a joint analysis on variables that are never jointly observed 20,21 . We will explore robustness of our methods to this conditional independence assumption using sensitivity analyses.
These assumptions will be essential to facilitating identification and estimation of causal quantities of interest. We propose to impute the missing values in X and R as well as the missing potential outcomes (Y 0 and Y 1 ) iteratively variable-by-variable using a sequence of conditional regression models until convergence: for any j 34,35,43 . Note that X = (X 1 , · · · , X k ) and X (−j) = (X 1 , · · · , X (j−1) , X (j+1) , · · · , X k ) denoting the collection of the k − 1 variables in X except X j . This method is the sequential regression multiple imputation method and has its advantages in handling nonnormal data 34 . We can evaluate how well the missing variables are imputed as in regular missing data imputation problems.
We impute the missing values D times to give a range for each imputed value to reflect prediction uncertainty. We then obtain the estimates of the quantities of interest for each of the D imputed data sets before combining these estimates and their variances across D imputed data sets using Rubin's formula 38 . Let ∆ represent any quantity of interest, ∆ d the point estimate and W d the withinimputation variance estimate of ∆ d for the dth imputed data set. The POMI estimate of ∆ is obtained by averaging across D imputed data sets, such that . Note that W represents the average within-imputation variance and B represents the between-imputation variance. The estimated ∆ follows a t distribution, with the degree of freedom (D − 1)[1 + D D+1 (W /B) 2 ] 42 . Specifically, the within-imputation variance for the estimated causal quantity of interest, ATE, for each imputed data set is calculated as follows: 1 and N is the sample size. The calculation of the withinimputation variance of the MOR estimate is more complicated than ATE because it is a non-linear function of P 0 and P 1 . Since MOR is skewed and positive, we calculate the variance for the MOR estimate on a log scale. Specifically, for each imputed data set, we calculate the variance of log( M OR) using the delta method as follows: . From the variance of log( M OR), we can then calculate the variance of the estimated MOR using the delta method again. The details for the variance calculation of estimated MOR are provided in the appendix. The within-imputation variances for estimated P (00), P (01), P (11) and P (10) are straightforward to calculate as multinomial proportions.
What we described so far are standard multiple imputation procedures. We will explore several features of causal inference that are distinct from standard missing data problems and see how the incorporation of them in the multiple imputation process impacts the causal estimates. These features will be incorporated through the four different POMI methods investigated in the following Simulation section. First, we will explicitly incorporate the causal inference assumptions (listed above as Assumptions 2, 3 and 4) and distinct data structures into the multiple imputations of missing potential outcomes. Under Assumption 2 and 3, we can simplify our imputation models for Y 0 and Y 1 such that . Under Assumptions 2, 3 and 4, we can further simplify our imputation models for Y 0 and Considering R only exists for testers, we can further simplify the models such that f (Y 0 |Y 1 , Z, X, R) = f (Y 0 |X). Second, we will explore the use of Z in the multiple imputation process. Z is a missing data indicator for Y 0 and Y 1 . With an MAR assumption, multiple imputation does not need to include the missing data indicator 35 . However, here Z also serves the dual role as the exposure of interest. Third, we will investigate the use of test results R, a correlate of Y and a post-exposure variable. In causal inference, post-exposure variables are not usually utilized because they may be in the causal pathway, and adjusting for them can bias results 44 . However, in standard missing data problems, postexposure variables are usually included in the imputation models to increase the predictibility of variables with missing values, when post-exposure variables are correlated with the variables with missing values and hence contain additional information about them.

Simulations
We conduct simulations to examine: 1) the impact of Z, R and causal inference assumptions on estimating the quantities of interest through comparisons of four POMI methods, which differ by whether the imputation models include Z, and/or R or make certain causal inference assumptions; 2) the properties of causal estimates of ATE and MOR estimated by these four POMI methods and the inverse probability weighting method (IPW); 3) the properties of the unique causal estimates that are obtainable from the four POMI methods; 4) the robustness of the POMI methods towards model misspecification in comparisons with IPW; 5) sensitivity analyses towards the potential outcomes conditional independence assumption (Assumption 4); 6) the properties of the estimated association between patient characteristics and potential outcomes subgroup memberships through the POMI methods.

Simulation Design
We assume that X includes two variables, X 1 and X 2 . The distributions of X 1 , We further categorize R * into three groups to mimic the distribution of RS assay results in our data example.
Note that U is an unmeasured predictor of Y , which is used to induce a non-zero correlation between Y 0 and Y 1 after adjusting for X and hence a violation of the potential outcomes conditional independence assumption. We assume U ∼ Normal(0,1). The default model parameter specifications in our simulations are: α = 0 (i.e., no interaction between X 1 and X 2 ), and β 1 = β 2 = 0 (i.e., conditional independence between potential outcomes). However, we will let α = 0 when conducting analyses to test the robustness of our methods towards model misspecifications. We will also allow β 1 and β 2 to be non-zero when conducting sensitivity analyses on the potential outcomes conditional independence assumption (Assumption 4). Each simulated data set consists of X 1 , X 2 , Z, R and Y . To be consistent with reality, the data set used in the data example does not have information on U , or on R when Z = 0. Further, we allow X 2 and Y be partially missing. Let M x and M y be 0 or 1, indicating X 2 and Y being available or missing, respectively. The models for the probability of missing X 2 and Y are specified as follows: logit[P (M x = 1)] = −1 + 0.2X 1 and logit[P (M y = 1)] = −2.5 + X 1 . We assume X 1 is completely measured. For each parameter specification, we simulate 500 data sets, each of which has N = 1, 000 subjects. Unless specified otherwise, the default scenarios in our data simulations satisfy Assumptions 1-4. Additional simulations with N = 500 and 250 subjects are also conducted to assess the impact of different sample sizes and summarized in Appendix.
We investigate the performance of four POMI methods. These imputation methods differ by whether the imputation models include Z, R or explicitly make causal assumptions of ignorabilty of exposure and potential outcomes independence (Assumptions 3 and 4). We will investigate the impact of these choices on estimates of the causal quantities of interest. These choices are reflected by different conditional regression models as described below. All four methods make the MAR and SUTVA assumptions (Assumptions 1 and 2), and are carried out by the IVEWARE software 34 .
1. Potential outcome imputation method 1 (POMI-Z): this is the default imputation procedure carried out by IVEWARE which consists of the following sequence of conditional regression models until convergence: Note that this imputation method does not make Assumptions 3 and 4, or utilize Z (the exposure variable and the missing indicator for Y 0 and Y 1 ), but it utilizes the post-exposure variable R.
2. Potential outcome imputation method 2 (POMI): different from POMI-Z, we now include Z in the imputation of missing X 2 , using the following sequential regression models: The rationale of using the test status Z in imputing missing X 2 is that, the probability of receiving the test Z is associated with X 2 ; hence, Z contains information about X 2 . Similar to POMI-Z, this method does not make Assumptions 3 and 4, but it utilizes R.
3. Potential outcome imputation method 3 (POMI+IND): different from POMI, we now incorporate Assumptions 3 and 4, and consider the fact that test results R are only available among testers. The sequential regression models are simplified as follows: . This is also consistent with the default setup in our simulations.
4. Potential outcome imputation method 4 (POMI+IND-R): different from POMI+IND, we now exclude the post-exposure variable R, using the following sequential regression models: Note that f (X 2 |.) represents a conditional probability density function modelled by a linear regression for X 2 , f (R|.) represents a conditional probability mass function modelled by a multinomial regression for R, f (Y 0 |.) and f (Y 1 |.) represent conditional probability mass functions modelled by logistic regression for Y 0 and Y 1 respectively. We iteratively impute the missing values in X 2 , R, Y 0 and Y 1 , with the sequential conditional regression models until convergence. To evaluate how well the missing variables are imputed, we compare the distributions of observed and imputed values conditional on the estimated propensity of being observed 45 . We repeat the imputation process 5 times to have 5 imputed data sets. We then conduct statistical inference on these 5 imputed data sets.
We compared these four methods with the inverse probability weighting method (IPW) 46 which is a propensity-score based method and commonly used to obtain causal inference. Before using IPW, we need to address the issue of missing data first. We use the same sequential regression multiple imputation method to impute missing values in the observable variables X 2 and Y 5 times with the following sequence of conditional regression models: f (X 2 |X 1 , Y, Z), f (Y |X 1 , X 2 , Z). For each of the 5 imputed data sets, we calculate the propensity score (denoted by π i = P (Z i = 1)) from a logistic regression of Z as a function of X 1 and X 2 . We then estimate P (Y 1 = 1) and P (Y 0 = 1) as Hence, the estimated ATE is AT E IP W = P 1,ipw − P 0,ipw . The marginal odds ratio is estimated as M OR IP W = P 1,ipw (1 − P 0,ipw ) / P 0,ipw (1 − P 1,ipw ) . The variance calculations of the estimated ATE and MOR use the delta method. For MOR, the log transformation is used to obtain the variance of log( M OR) first before obtaining the variance of M OR, similar to what has been described in the Appendix. Finally, Rubin's formulae are used to combine the estimates and their variances across 5 imputed data sets.
The true values for the causal quantities of interest cannot be obtained analytically but can be closely approximated by simulating 10 million patients with complete information on both Y 0 and Y 1 , based on which the true causal quantities can be obtained. This approach has often been used in causal inference research. For each set of parameter specifications, we simulate 500 data sets. For each quantity of interest, we estimate the mean of each estimate, bias (the difference between the estimate and the true value), coverage probability (CR) of the 95% confidence interval, standard error (SE), empirical standard deviation (ESD) and mean squared error (MSE) across 500 data sets.

Simulation Results
First, we evaluate the properties of estimates for ATE and MOR from these five different methods (Table 2). These quantities of interest are obtainable from all methods. Based on the evaluation, POMI-Z has larger bias and lower coverage rates while all other methods give negligible bias and close to nominal 95% coverage rates. The IPW method tends to give slightly bigger variance estimates; while all the POMI methods give similar variance estimates. Additional simulations with N = 500 and 250 show that the bias of the estimates becomes larger as the sample size decreases in most scenarios; however, they demonstrate similar patterns across different methods as N=1000 (See Tables  8-11 in Appendix).
Second, we evaluate the properties of the estimated proportions of patients in each potential outcome subgroup by these five methods (Table 3). We estimate the fraction of patients that are either never-treated P (00), treatment-encouraged P (01), treatment-discouraged P (10), or always-treated P (11). Note that IPW Table 2. Comparison of five methods in estimating average testing effect (ATE), and marginal odds ratios (MOR). ESD: empirical standard deviation; SE: average standard error; CR: coverage rate of 95% confidence interval. Here we assume potential outcomes conditional independence and no interaction between X1 and X2 (i.e., α = β1 = β2 = 0).

Quantities Methods
True cannot provide estimates of these proportions; hence, no estimates from IPW are provided. The simulation results indicate that POMI+IND and POMI+IND-R give estimates with smaller bias, smaller ESD and SE than POMI-Z and POMI. The coverage rates from POMI+IND and POMI+IND-R are also closer to the nominal 95% coverage rates. Note that the conservative coverage rates from POMI-Z and POMI result from the variance being severely overestimated, which overwhelms the bias in these methods. Hence, the results demonstrate the benefit of taking advantage of the causal assumptions of ignorability and potential outcomes conditional independence when we estimate the proportions of potential outcome subgroups. The performances of POMI+IND and POMI+IND-R are very similar to each other, indicating that there is a minimal advantage in utilizing the information from the post-exposure variable R to estimate these quantities in our simulations. Additional simulations with N = 500 and 250 show that the bias of the estimates becomes larger as the sample size decreases in most scenarios; however, they demonstrate similar patterns across different methods as N=1000 (See Tables 8-11 in Appendix).
Third, we evaluate the performance of a generalized logit model for potential outcomes subgroup membership using POMI+IND ( Table 4). The model is specified as log[ P (ij) P (00) ] = γ 0ij + γ 1ij X 2 + γ 2ij X 1 , where ij = 01, 10, or 11 and γ 0ij , γ 1ij , γ 2ij are the corresponding coefficients. Since POMI+IND demonstrates advantages over POMI-Z and POMI and similarity as POMI+IND-R in performance based on prior simulations (Tables 2 and 3), we focus on the performance of POMI+IND in estimating the coefficients of the generalized logit model for the potential outcomes subgroup membership. Our simulations demonstrate that the POMI+IND method performs well in estimating the model Table 3. Comparison of five methods in estimating the proportion of patients in each potential outcome subgroup. ESD: empirical standard deviation; SE: standard error; CR: coverage rate of 95% confidence interval. Here we assume potential outcomes conditional independence and no interaction between X1 and X2 (i.e., α = β1 = β2 = 0).

Quantities Methods
True Fourth, we continue to focus on the POMI+IND method and examine the impact of an unmeasured predictor U of Y on the quantities of interest ( Figure  1). The presence of U introduces a correlation between Y 0 and Y 1 and subsequently a violation of Assumption 4 (e.g. potential outcome conditional independence). We let β 1 = β 2 and vary their sizes to be 0, 0.9, 1.06, 1.2, 1.38 and 1.49, which correspond to the proportions of variation of Y explained by U being 0%, 5%, 10%, 15%, 20% and 25%, respectively. As the correlations between Y 0 and Y 1 increase, the violation of Assumption 4 becomes more severe. Our simulations show that U has little impact on the bias of ATE and MOR. However, as U explains more of the variance in Y , Y 0 and Y 1 become increasingly more correlated; subsequently, the treatment-encouraged and treatment-discouraged subgroups become more rare and the bias of the estimated proportion in each potential outcome subgroup becomes more severe. Table 4. Properties of model coefficients from generalized logit model for potential outcomes subgroup membership using the POMI+IND method. ESD: empirical standard deviation; SE: standard error; CR: coverage rate of 95% confidence interval. Here we assume potential outcomes conditional independence and no interaction between X1 and X2 (i.e., α = β1 = β2 = 0).

Model
Parameter  Table 5. Comparing the robustness of POMI+IND and IPW towards model mis-specifications in estimating average testing effect (ATE) and marginal odds ratio (MOR). ESD: empirical standard deviation; SE: standard error; CR: coverage rate of 95% confidence interval. Here we assume potential outcomes conditional independence (i.e., β1 = β2 = 0) and there is an interaction between X1 and X2 (i.e., α = 0). Lastly, we examine the robustness of the POMI+IND method towards model misspecification. We compare the method with the IPW method. For IPW to work properly, the model for the test assignment Z needs to be correctly specified. Here, we specify the true model for Z as logit[P (Z = 1)] = −1.5 − X 1 + 0.5X 2 + 2X 1 × X 2 . However, in data analyses, neither IPW nor POMI+IND methods include the interaction term X 1 × X 2 . The simulation results are summarized in Table 5. It shows that IPW gives biased estimates of ATE and MOR and much below 95% coverage rates for their confidence intervals; however, the model misspecification for Z has little impact on estimates from the POMI+IND method. This highlights the need for IPW to specify correct models for Z; however, this is not required for the imputation method.  Figure 1. Relative bias of estimated average testing effect (ATE), marginal odds ratio (MOR), and proportion of patients in each potential outcomes subgroup, as the variance of Y explained by an unmeasured predictor of Y , U , increases and subsequently violation of potential outcomes conditional independence increases (i.e., β1 = 0 and β2 = 0). Here we assume no interaction between X1 and X2 (i.e., α = 0). The true values for ATE is -0.186, MOR is 0.458 and the proportions of patients in four potential outcomes subgroups, P (00), P (01), P (10) and P (11), are 0.212, 0.104, 0.290, and 0.394, respectively.

Application: the Impact of the RS Assay on Breast Cancer Chemotherapy Selection
The analysis data set includes 2,931 incident breast cancer patients in Georgia and Los Angeles who were diagnosed between 2013-2014 with favorable prognosis (estrogen positive, human epidermal growth factor receptor 2 negative and invasive tumor behavior) 18,31 . These patients were chosen because the RS assay has the potential to make a difference in their chemotherapy treatment decisions. The data set comes from three data sources: 1) SEER databases for clinical features; 2) patient survey data from the ICanCare study; 3) the RS assay use and its testing results from the laboratory, Genomic Health, Inc. The outcome of interest is the receipt of chemotherapy (yes/no). The exposure of interest is the receipt of the RS assay (yes/no). The covariates include age (year), tumor grade (1, 2, or 3), tumor size (≤1 cm, 1 − 2 cm or ≥2 cm), node status (positive or negative), the number of major co-morbidities, menopausal status, and the risk of genetic mutation (low or high), education level (high school or less, some college, graduate degree), race (White, Black, Latino, Asian, other), insurance type (medicaid, medicare, private, other) and family income.
Overall, 35% of the patients received the RS assay and 27% received chemotherapy. Among these RS assay testers, the test results showed that 65%, 27%, and 8% of patients received low-, intermediate-and high-risk scores respectively. We first fit a logistic regression with above covariates without addressing missing data. We find that the factors that are significantly associated with the use of chemotherapy include: the receipt of the RS assay (odds ratio  Table 12 in Appendix).
We then apply the POMI+IND method to obtain causal estimates. We estimate that the average causal effect of the RS assay testing on chemotherapy (i.e., ATE) is −10% with its 95% CI of (−14%, −7%) and the MOR is 0.59 with its 95% CI of (0.49, 0.70). Table 6 presents the proportion of patients in each of the four potential outcome subgroups by overall and by covariates. Overall, it is estimated that the receipt of the RS assay encouraged a change of treatment plan in 9%(8%, 11%) of patients from no chemotherapy to chemotherapy, and 20%(17%, 22%) of patients from chemotherapy to no chemotherapy. Hence, the receipt of RS assay has reduced the overall chemotherapy use by 10% and changed the treatment plan for a total of 30% of patients. Regardless of whether they were tested or not, 59% of patients would never receive chemotherapy and 12% would always receive chemotherapy.
In addition, we find that, independent of the use of the RS assay, older patients (≥65 years old) are more likely to reject chemotherapy than younger patients (73% vs. 49% estimated to be in the "never-treatment" subgroup). We also find that, with the use of the RS assay, patients were more likely to be discouraged from chemotherapy if they had a larger tumor (36% vs. 19%, 8% in the "treatmentdiscouraged" subgroup for patients with tumor size of ≥ 2cm vs. ≤ 1cm, 1-2cm) or if they had positive nodes (40% vs. 15% for patients with positive vs. negative nodes). It is particularly noteworthy that even though the guidelines did not make any recommendations about the use of the RS assay for node-positive patients, 14% of them were offered the RS assay in practice. And in fact these patients were more likely to change their chemotherapy plans after testing with the RS assay than node-negative patients (46% vs. 25% changed for patients with positive vs. negative nodes).
It is also evident that the RS assay played an influential role in chemotherapy decisions based on how much the results of the assay changed the treatment selections: for patients received a low-risk score, 2% were encouraged to receive chemotherapy vs. 22% discouraged from chemotherapy; for patients received a high-risk score, 34% were encouraged vs. 5% discouraged. In contrast, for patients received an intermediate-risk score, 19% were encouraged vs. 17% discouraged.
We also fit a multivariable multinomial logistic regression to examine the relationship between several key patient clinical predictors and the potential outcomes subgroup membership. These clinical predictors include tumor node status, tumor grade, tumor size, and age. They are the most important predictors for our study population and have been endorsed by the American Joint Committee on Cancer 47 . Table 7 listed the estimated probability of patients' potential outcome subgroup membership as a function of these clinical factors. We estimate that, for patients with negative node, tumor grade 1, tumor size ≤ 1cm, and age<65 years old, there was a 7% chance to be encouraged to receive chemotherapy and 5% chance to be discouraged from chemotherapy by using the RS assay, an 88% chance of staying with the decision of no chemotherapy and a 1% of chance of staying with chemotherapy decision regardless of the use of the RS assay. We also estimate that, for patients with negative node, tumor grade of 3, tumor size of 1-2cm and age < 65 years old, the chances of changing to chemotherapy from no chemotherapy, from no chemotherapy to chemotherapy or staying with chemotherapy or no chemotherapy after the use of the RS assay are 12%, 28%, 52% and 8% respectively. We also found that, generally, patients with a mix of low-risk factors (e.g., smaller tumor size, lower tumor grade, negative node) and high-risk factors (e.g., larger tumor size, higher tumor grade, positive node) for recurrence are more likely to be affected by using the RS assay than patients with either all low-risk or high-risk factors. For example, node-positive patients with small and low-grade tumors are most likely to be discouraged from chemotherapy by using the RS assay. Similarly, node-negative patients with small but high-grade tumor are most likely to be encouraged to take chemotherapy by using the RS assay. In contrast, patients with consistently low (high) risk clinical factors are more likely to remain in the "never treated" ("always treated") subgroup.
We further examine how robust our analysis results are towards the violation of the potential outcomes conditional independence assumption (Assumption 4). Since Y 0 and Y 1 are not observed simultaneously in our data, we cannot evaluate this conditional independence assumption nor estimate the correlation between Y 0 and Y 1 based on our data. We artificially induce a correlation by simulating an additional predictor U of Y that is not in the data set, and examine how sensitive our analyses are towards the correlation between Y 0 and Y 1 . Since U is a predictor of Y , U is a common predictor of both Y 0 and Y 1 . Since U is not observed or Table 6. Estimated proportion of patients in each of the four potential outcomes subgroups defined by the impact of the RS assay on chemotherapy according to patient characteristics, test results and insurance status.

Discussion
We proposed a potential outcomes framework to measure the influence of medical testing on treatment selection, stratifying patients into four potential outcomes subgroups based on whether the test has any impact on treatments election (either undergoing the treatment or not), and whether the test encourages or discourages treatment selection. This framework captures well the differential causal effect of testing on treatment selection for different patient subgroups, and allows us to examine the distribution of these subgroups and characteristics associated with these subgroup memberships. The classification suggests targets for more efficient utilization of medical tests. With advanced tests designed to support precision medicine becoming increasingly available and increasing concern for over-testing 3 , the causal estimates produced from this framework become more relevant. The framework also allows us to estimate many traditional causal quantities of interest (e.g., average causal effects, marginal odds ratios) on which existing causal inference methods more often focus. By taking advantage of the connection between missing data and causal inference, we can simultaneously implement the causal framework as well as handle missing data in the observable variables (i.e., due to non-response). We also explore the advantages of incorporating unique features of causal inference assumptions in the data imputation process. These assumptions are often made in causal inference literature, but cannot be estimated or verified from the data. Simulation studies suggested that the proposed methods provide reasonably accurate inference. Additionally, the explicit incorporation of the causal assumptions in the data imputation process can improve precision for some causal estimates. We also find that bias can occur for some causal estimates when the potential outcomes conditional independence assumption is violated, which can be assessed through sensitivity analyses. The proposed framework was then used to examine how the most commonly used genomic test in the United States, 21-gene assay, impacts the use of chemotherapy among breast cancer patients.
Existing research has often focused on estimating the overall relationship between medical testing and treatment options in a population. For example, breast magnetic resonance imaging has been found to be associated with receiving a mastectomy instead of breast conserving surgery (adjusted odds ratio = 1.21) 49 . In another example, it is found that 55% of women with invasive breast cancer and positive HER2 genetic testing results received trastuzumab treatment 50 . Although this information is of interest, the average effect of testing on treatment choice may not apply equally across subgroups of patients who may not benefit from (or, may be harmed by) testing. Hence, it is important to recognize that a medical test can have different influences on treatment decisions for different patients; i.e., there may be effect heterogeneity of testing on treatment. We proposed the use of potential outcomes subgroups to classify patients by the presence (vs. absence) and direction of the effect of testing on treatment decisions. Patients are classified into one of the four groups: 1) never treated; 2) treatment encouraged; 3) treatment discouraged; 4) always treated. These subgroups naturally capture the effect heterogeneity of medical testing on treatment decisions. Additionally, we investigated how patients' membership in these subgroups may be driven by patient characteristics or systemic factors. Such classification would suggest targets for improving the utilization of medical tests and the allocation of limited resources. For example, in the data analyses, we found that patients with positive nodes were more influenced by the testing and more of these patients were discouraged from receiving chemotherapy than patients with negative nodes, even though guidelines make no recommendations for how to use the 21-gene assay in node-positive patients. We also find that patients with a mix of low and high risk factors were more likely to be influenced by the use of 21-gene assay, and the tests were more likely to either encourage or discourage these patients for chemotherapy use. Given that the test was recommended to aid decision-making in people with intermediate risk of recurrence, our results show that efforts to ensure that the patients with a mix of low and high risk factors are offered testing are likely to produce substantial changes in plans and maximize the utility of the test.
Our methods differ from traditional latent class models, which usually deal with multiple or longitudinal outcomes and do not classify patients by potential outcomes. Subgroup analysis or adding interaction terms can also investigate effect heterogeneity; however, these approaches require us to pre-specify variables that modify effect and suffer from false discovery rate 51 . In summary, standard analyses cannot provide the classification targeted by our methods. Our methods of analyzing heterogeneous testing effect using the potential outcomes subgroups approach are at the individual level. For each patient, the testing effect on treatment is obtained by comparing her underlying tendency of receiving treatment with vs. without a given test. Our approach may be complementary to the standard analyses based on known effect modifiers. Further discussions of individual heterogeneous treatment effects are discussed by others 39,52,53 .
We proposed the use of multiple imputation methods to impute the missing potential outcomes for each patient (Y 0 or Y 1 ) since we do not observe two potential outcomes of the same patient simultaneously. Subsequently, we can classify patients into four potential outcomes subgroups according to individual heterogeneous treatment effects. We chose the sequential regression multiple imputation method. It has the advantages of accommodating many different data types and structure compared with other missing data methods 54 . The main uncertainty in theoretical justification for this method is the likely incompatibility between the assumed conditional models and the true conditional models derived from the joint conditional distributions 34 . However, many simulation studies have demonstrated the robustness of this method towards possible misspecifications of the assumed conditional models 54 .
Our methods have the advantage of imputing missing potential outcomes and regular missing values in the observable variables simultaneously during the data imputation process. This is more convenient and faster than most of other causal inference methods in the presence of missing data, which will likely need to handle the missing data problem first using multiple imputation methods before applying causal inference methods and combining analysis results across multiple imputed data sets.
The use of multiple imputation methods utilizes the close connection between causal inference and missing data framework, since the causal inference problems here are essentially a missing data problem. However, there are some unique features of causal inference that are different from the regular missing data problems in our context. First, causal assumptions are unique. Even though there is no information from the data to test the validity of the assumptions regarding ignorability or potential outcomes conditional independence, we demonstrate that the explicit incorporation of these causal assumptions into the data imputation process can increase the precision of some of our causal estimates such as the proportion of patients in each of the potential outcomes subgroups. Second, we find that the violation of the potential outcomes conditional independence assumption does not have much impact on traditional causal quantities such as the average causal effect or marginal odds ratio. However, it does impact the causal quantities regarding the potential outcomes subgroups, which concern the joint distribution of the potential outcomes. As the correlation between two potential outcomes increases, patients in the treatment-encouraged and treatment-discouraged subgroups become fewer. Since there is no information from the empirical data about the validity of this assumption, we proposed a method to conduct sensitivity analyses to relax this assumption and evaluate the impact of violation of this assumption. Our proposed sensitivity analyses are closely related to the random effect models proposed by Zhang et al. 39 , which use random effects to induce the correlation between Y 0 and Y 1 . Here we used an unmeasured predictor of Y to induce the correlation. Third, Z is the missing data indicator for Y 0 and Y 1 . With an MAR assumption, multiple imputation does not usually need to include the missing indicator 35 . However, Z is also the exposure of interest. And it is a function of the covariate X 2 , which is partially missing. Hence, for multiple imputation of X 2 to work properly, the imputation model for X 2 needs to include Z to satisfy the MAR assumption. This explains why POMI-Z performs more poorly in our simulations compared with methods including Z in the imputation models. Fourth, a post-exposure variable (e.g., testing results) that is correlated with the outcome is not usually utilized in the causal inference methods because it can be in the causal pathway between the exposure and outcome. However, it is usually used in the multiple imputation process and not treated differently from other variables. Our study finds that the incorporation of a post-exposure variable in the multiple imputation process still generates valid causal estimates. In fact, it did not make much difference whether or not we included the post-exposure variable in our simulations. However, we believe that if the post-exposure variable is correlated with variables that have more missingness, the utilization of the post-exposure variable will increase precision. In our simulations, there is only 10% missingness in the variable Y , a correlate of the post-exposure variable R; hence, the incorporation of R in the imputation process did not provide noticeable efficiency gain. Lastly, we find that our POMI methods do not require the need to correctly specify the model for the exposure to obtain valid causal estimates; however, a propensity-score based causal method usually does. This feature occurs because the POMI methods do not require one to model the exposure but only to provide a correct predictive model for the outcome. Our methods can be implemented through IVEware in SAS, or through the MICE (Multivariate Imputation by Chained Equations) package in R. We have provided code examples in GitHub at https://github.com/yunliyunli/POMI. In summary, we proposed a potential outcomes framework to classify patients according to the differential causal effect of medical testing on treatment selection for different individual patients. It has the potential to be very useful in examining effect heterogeneity of medical testing on treatment selection and in proposing targets for the optimal use of medical tests. Future research includes conducting more sensitivity analyses on the causal assumptions and improving the efficiency for these causal quantities of interest.

Declaration of conflicting interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
This work was supported by the National Cancer Institute through grants P01CA163233, CA46592, CA129102 and U057141. Table 8. Additional simulations: comparison of five methods in estimating average testing effect (ATE), and marginal odds ratios (MOR) when N=250. ESD: empirical standard deviation; SE: average standard error; CR: coverage rate of 95% confidence interval. Here we assume potential outcomes conditional independence and no interaction between X1 and X2 (i.e., α = β1 = β2 = 0).

Quantities Methods
True  Table 9. Additional simulation: comparison of five methods in estimating the proportion of patients in each potential outcome subgroup when N=250. ESD: empirical standard deviation; SE: standard error; CR: coverage rate of 95% confidence interval. Here we assume potential outcomes conditional independence and no interaction between X1 and X2 (i.e., α = β1 = β2 = 0).  Table 10. Additional simulation: comparison of five methods in estimating average testing effect (ATE), and marginal odds ratios (MOR) when N=500. ESD: empirical standard deviation; SE: average standard error; CR: coverage rate of 95% confidence interval. Here we assume potential outcomes conditional independence and no interaction between X1 and X2 (i.e., α = β1 = β2 = 0).  Table 11. Additional simulation: comparison of five methods in estimating the proportion of patients in each potential outcome subgroup when N=500. ESD: empirical standard deviation; SE: standard error; CR: coverage rate of 95% confidence interval. Here we assume potential outcomes conditional independence and no interaction between X1 and X2 (i.e., α = β1 = β2 = 0).

Variance Estimates of Marginal Odds Ratio
The MOR is defined as MOR= [P 1 /(1 − P 1 )]/[P 0 /(1 − P 0 )]. Since MOR is skewed and bounded to be positive, we use log transformation to make it more symmetric so that we can use the delta method to obtain the variance estimate of log( M OR)