Using Propensity Scores to Develop and Evaluate Treatment Rules with Observational Data

TL;DR

This paper presents a flexible method for estimating individualized treatment rules from observational data using propensity scores, addressing challenges of unmeasured confounders and providing tools for evaluation.

Contribution

It introduces a novel estimation framework that leverages observation weights based on propensity scores, applicable to various prediction methods, and offers a way to assess the population benefit of treatment rules.

Findings

Framework accommodates any prediction method with observation weights

Provides a method to estimate population benefit of treatment rules

Implementation available in R package DevTreatRules

Abstract

In this paper, we outline a principled approach to estimate an individualized treatment rule that is appropriate for data from observational studies where, in addition to treatment assignment not being independent of individual characteristics, some characteristics may affect treatment assignment in the current study but not be available in future clinical settings where the estimated rule would be applied. The estimation framework is quite flexible and accommodates any prediction method that uses observation weights, where the observation weights themselves are a ratio of two flexibly estimated propensity scores. We also discuss how to obtain a trustworthy estimate of the rule's population benefit based on simple propensity-score-based estimators of average treatment effect. We implement our approach in the R package DevTreatRules and share the code needed to reproduce our results on GitHub.

Figures

Tables

Table 1. Table 1: Proposed partitioning of individual characteristics in an observational study

Potentially influences treatment	Observed in
in current study	independent clinical settings
	Yes	No
Yes	$C^{\text{TI}}$	$C^{\text{TN}}$
No	$C^{\text{NI}}$	$C^{\text{NN}}$

Table 2. Table 2: Selected methods for estimating treatment rule 1 Holloway et al., ( 2019 ) 2 Huling et al., ( 2019 )

Framework	Type of approach	Accommodates observational data	Distinguishes between $(C^{\text{TI }},C^{\text{NI }})$ and $\left(C^{\text{TN }}\right)$	Accommodates a range of statistical learning methods	R package	Shares code to reproduce simulations or application
Zhao et al., (2012)	Direct	Yes	No	No	Yes1	No
Zhang et al., 2012a ; Zhang et al., 2012b	Direct	Yes	No	Yes	No	Yes
Zhang et al., (2015)	Direct	Yes	No	No	No	Yes
Qian and Murphy, (2011)	Direct	No	No	No	No	No
Chen et al., (2017)	Direct	Yes	No	Yes	Yes2	Yes
Cai et al., (2011)	Indirect	No	No	Yes	No	No
Lu et al., (2013)	Indirect	No	No	Yes	No	No
Kang et al., (2014)	Indirect	No	No	Yes	No	Yes
McKeague and Qian, (2014)	Indirect	No	No	Yes	No	No
Ciarleglio et al., (2015)	Indirect	Yes	No	No	No	Yes
This paper	Indirect	Yes	Yes	Yes	Yes	Yes

Table 3. Table 3: Mean outcome probability, as a function of rule type and development set sample size, averaged over 1000 replications and calculated in evaluation sets of size 10000

	Sample size in development set
Type of Rule	50	100	200	500	1000
Split-regression	0.543	0.553	0.562	0.57	0.572
Split-regression (naive, no weights)	0.552	0.554	0.553	0.55	0.549
Optimal rule	0.574	0.574	0.574	0.574	0.574
Treating all	0.479	0.479	0.479	0.479	0.479
Treating none	0.543	0.543	0.543	0.543	0.543

Table 4. Table 4: Summary of Selected Rules on the Validation and Evaluation Sets (Outcome: No breast cancer after 10 years). The selected propensity method/rule method are in parentheses

	Positives	Negatives	ATE in Positives	ATE in Negatives	ABR
Estimates on Validation Set
Split regression (ridge/lasso)	3544	15758	0.013	-0.051	0.044
Split regression (logistic/logistic)	3692	15610	0.021	-0.053	0.047
Treat no one (logistic/NA)	0	19302	NA	-0.045	0.045
Estimates on Evaluation Set
Split-regression (ridge/lasso)	3408	15894	0.015 (-0.027, 0.068)	-0.045 (-0.067, -0.031)	0.04
Split-regression (logistic/logistic)	3569	15733	0.002 (-0.049, 0.052)	-0.048 (-0.07, -0.035)	0.04
Treat no one (logistic/NA)	0	19302	NA	-0.05	0.05

Table 5. Table 5: Summary of Outcome and Treatment Variables

	Overall (93676)	Non-event (3063)	Event (90613)	N missing
Outcomes
No CHD within 10 years of enrollment, n (%)	90613 (97%)	-	-	0
No breast cancer within 10 years of enrollment, n (%)	88883 (95%)	-	-	0
Treatment
Currently using unopposed estrogen and/or estrogen plus progesterone, n (%)	41630 (44%)	1003 (33%)	40627 (45%)	85

Table 6. Table 6: Summary of Variables Influencing Only Treatment Assignment

	Overall (93676)	Non-event (3063)	Event (90613)	N missing
Highest grade completed				767
None, n (%)	84 (0%)	2 (0%)	82 (0%)
1-4, n (%)	356 (0%)	11 (0%)	345 (0%)
5-8, n (%)	1121 (1%)	51 (2%)	1070 (1%)
9-11, n (%)	3288 (4%)	184 (6%)	3104 (3%)
High school, n (%)	15122 (16%)	592 (19%)	14530 (16%)
Vocational, n (%)	9123 (10%)	369 (12%)	8754 (10%)
Some college, n (%)	24812 (27%)	828 (27%)	23984 (27%)
College, n (%)	10669 (11%)	277 (9%)	10392 (12%)
Some post-graduate, n (%)	11018 (12%)	314 (10%)	10704 (12%)
Master’s, n (%)	14732 (16%)	343 (11%)	14389 (16%)
Doctoral, n (%)	2584 (3%)	67 (2%)	2517 (3%)
Ethnicity				265
American Indian or Alaskan Native, n (%)	421 (0%)	18 (1%)	403 (0%)
Asian or Pacific Islander, n (%)	2671 (3%)	52 (2%)	2619 (3%)
Black or African-American, n (%)	7635 (8%)	283 (9%)	7352 (8%)
Hispanic/Latino, n (%)	3609 (4%)	56 (2%)	3553 (4%)
White (non-Hispanic), n (%)	78016 (84%)	2611 (86%)	75405 (83%)
Other, n (%)	0 (0%)	28 (1%)	1031 (1%)
Heard about study				1281
Mailed letter, n (%)	47623 (52%)	1722 (57%)	45901 (51%)
Brochure, n (%)	9789 (11%)	317 (10%)	9472 (11%)
TV, n (%)	2731 (3%)	93 (3%)	2638 (3%)
Radio, n (%)	1017 (1%)	24 (1%)	993 (1%)
Newspaper or magaize, n (%)	14610 (16%)	415 (14%)	14195 (16%)
Meeting, n (%)	1158 (1%)	32 (1%)	1126 (1%)
Friend or relative, n (%)	9408 (10%)	223 (7%)	9185 (10%)
Other, n (%)	6059 (7%)	198 (7%)	5861 (7%)
Family income				4119
Less than $10,000, n (%)	3917 (4%)	248 (8%)	3669 (4%)
$10,000 - $19,999, n (%)	10101 (11%)	504 (17%)	9597 (11%)
$20,000 - $34,999, n (%)	20226 (23%)	838 (29%)	19388 (22%)
$35,000 - $49,999, n (%)	17430 (19%)	536 (18%)	16894 (19%)
$50,000 - $74,999, n (%)	17487 (20%)	409 (14%)	17078 (20%)
$75,000 - $99,999, n (%)	8181 (9%)	169 (6%)	8012 (9%)
$100,000 - $149,999, n (%)	6034 (7%)	73 (3%)	5961 (7%)
$150,000 or more, n (%)	3393 (4%)	42 (1%)	3351 (4%)
Don’t know, n (%)	2788 (3%)	99 (3%)	2689 (3%)

Table 7. Table 7: Summary of Variables Influencing Treatment Assignment and Rule

	Overall (93676)	Non-event (3063)	Event (90613)	N missing
Age, mean (IQR)	63.6 (58, 69)	68.3 (64, 73)	63.5 (57, 69)	0
Angina ever, n (%)	5547 (6%)	584 (19%)	4963 (6%)	708
Aortic aneurysm ever, n (%)	187 (0%)	30 (1%)	157 (0%)	1523
Breast cancer ever, n (%)	5299 (6%)	208 (7%)	5091 (6%)	879
Coronary bypass surgery ever, n (%)	881 (1%)	204 (7%)	677 (1%)	1513
Cancer ever, n (%)	12075 (13%)	481 (16%)	11594 (13%)	752
Cardiac catheterization ever, n (%)	3837 (4%)	453 (15%)	3384 (4%)	1513
Carotid endarterectomy/angioplasty ever, n (%)	344 (0%)	59 (2%)	285 (0%)	1510
Cervix cancer ever, n (%)	1205 (1%)	44 (1%)	1161 (1%)	916
Heart failure ever, n (%)	893 (1%)	134 (4%)	759 (1%)	7
Cadiovascular disease ever, n (%)	17523 (19%)	1206 (40%)	16317 (18%)	2045
Hysterectomy ever, n (%)	39149 (42%)	1415 (46%)	37734 (42%)	87
Diabetes ever, n (%)	5318 (6%)	544 (18%)	4774 (5%)	96
Stroke ever, n (%)	1415 (2%)	142 (5%)	1273 (1%)	56
Have a lot of energy?				772
All the time, n (%)	4740 (5%)	108 (4%)	4632 (5%)
Most the time, n (%)	33633 (36%)	766 (25%)	32867 (37%)
A good bit, n (%)	20668 (22%)	642 (21%)	20026 (22%)
Some times, n (%)	19397 (21%)	780 (26%)	18617 (21%)
A little bit, n (%)	9956 (11%)	481 (16%)	9475 (11%)
Never, n (%)	4510 (5%)	258 (9%)	4252 (5%)
General health				655
Excellent, n (%)	16576 (18%)	263 (9%)	16313 (18%)
Very good, n (%)	37684 (41%)	861 (28%)	36823 (41%)
Good, n (%)	29669 (32%)	1280 (42%)	28389 (32%)
Fair, n (%)	8210 (9%)	550 (18%)	7660 (9%)
Poor, n (%)	882 (1%)	82 (3%)	800 (1%)
High cholesterol requiring pills ever, n (%)	13773 (15%)	748 (25%)	13025 (15%)	2071
Hot flash in past 4 weeks				733
No, n (%)	15158 (16%)	386 (13%)	14772 (16%)
Mild, n (%)	4593 (5%)	139 (5%)	4454 (5%)
Moderate, n (%)	1267 (1%)	38 (1%)	1229 (1%)
Severe, n (%)	0 (0%)	0 (0%)	0 (0%)
Hypertension				1699
Never, n (%)	61196 (67%)	1270 (42%)	59926 (67%)
Yes, untreated, n (%)	7317 (8%)	338 (11%)	6979 (8%)
Yes, treated, n (%)	23464 (26%)	1392 (46%)	22072 (25%)
Recent physical/emotional problems socially				747
Not at all, n (%)	68565 (74%)	2004 (66%)	66561 (74%)
Slightly, n (%)	14249 (15%)	549 (18%)	13700 (15%)
Moderately, n (%)	6121 (7%)	279 (9%)	5842 (6%)
Quite a bit, n (%)	3217 (3%)	169 (6%)	3048 (3%)
Extremely, n (%)	777 (1%)	34 (1%)	743 (1%)
Quality of life (1-10), mean (IQR)	8.3 (8, 9)	8.1 (7, 9)	8.3 (8, 9)	724
Menopause before age 40, n (%)	8352 (9%)	339 (12%)	8013 (9%)	3951
MENPSYMP, n (%)	64608 (71%)	1922 (65%)	62686 (71%)	2425
Limited in daily activities?				726
Yes, limited a lot, n (%)	6263 (7%)	442 (15%)	5821 (6%)
Yes, limited a little, n (%)	23110 (25%)	1132 (37%)	21978 (24%)
No, not limited at all, n (%)	63577 (68%)	1459 (48%)	62118 (69%)
One or both ovaries removed				552
No, n (%)	65240 (70%)	2019 (66%)	63221 (70%)
Yes, one taken out, n (%)	6583 (7%)	217 (7%)	6366 (7%)
Yes, both taken out, n (%)	18890 (20%)	713 (23%)	18177 (20%)
Yes, unknown number taken out, n (%)	738 (1%)	29 (1%)	709 (1%)
Yes, part of ovary taken out, n (%)	893 (1%)	30 (1%)	863 (1%)
Don’t know, n (%)	780 (1%)	36 (1%)	744 (1%)
Osteoporosis ever, n (%)	8282 (9%)	385 (13%)	7897 (9%)	1240
Any part of ovaries removed before age 40, n (%)	8279 (9%)	301 (10%)	7978 (9%)	1120
Peripheral arterial disease ever, n (%)	2084 (2%)	249 (8%)	1835 (2%)	784
Pregnant ever, n (%)	84005 (90%)	2760 (90%)	81245 (90%)	315
Angioplasty of coronary arteries ever, n (%)	1128 (1%)	177 (6%)	951 (1%)	1509
Stroke ever, n (%)	1415 (2%)	142 (5%)	1273 (1%)	56
Health limits vigorous activities				812
Yes, limited a lot, n (%)	30022 (32%)	1542 (51%)	28480 (32%)
Yes, limited a little, n (%)	41367 (45%)	1162 (38%)	40205 (45%)
No, not limited at all, n (%)	21475 (23%)	328 (11%)	21147 (24%)

Equations

$$\psi \equiv {\int }_{{\mathbf{C}}^{\text{T }}}\left\{E\left[Y\mid {\mathbf{C}}^{\text{T }},T=1\right]-E\left[Y\mid {\mathbf{C}}^{\text{T }},T=0\right]\right\}dP({\mathbf{C}}^{\text{T }}),$$

$$E[Y^1-Y^0\mid \mathbf{R}=\mathbf{r}].$$

$$\psi (\mathbf{r})\equiv {\int }_{{\mathbf{C}}^{\text{T }}}\left\{E\left[Y\mid {\mathbf{C}}^{\text{T }},\mathbf{R}=\mathbf{r},T=1\right]-E\left[Y\mid {\mathbf{C}}^{\text{T }},\mathbf{R}=\mathbf{r},T=0\right]\right\}dP({\mathbf{C}}^{\text{T }}|\mathbf{R}=\mathbf{r})$$

$$f_t(\mathbf{r})={\int }_{{\mathbf{C}}^{\text{T }}}E\left[Y\mid {\mathbf{C}}^{\text{T }},\mathbf{R}=\mathbf{r},T=t\right]dP({\mathbf{C}}^{\text{T }}|\mathbf{R}=\mathbf{r}).$$

$$f_t(\mathbf{r})\equiv \underset{f\in \mathcal{F}}{\mathrm{argmin}}{\int }_{\mathbf{R}}{\int }_{{\mathbf{C}}^{\text{T }}}{w}_t(\mathbf{r},{\mathbf{c}}^{\text{T}})\mathcal{L}(y,f(\mathbf{r}))dP(y,{\mathbf{c}}^{\text{T}},\mathbf{r}|T=t),$$

$$w_t(\mathbf{r},{\mathbf{c}}^{\text{T}})=\frac{P(T=t|\mathbf{R}=\mathbf{r})}{P(T=t|\mathbf{R}=\mathbf{r},{\mathbf{C}}^{\text{T }}={\mathbf{c}}^{\text{T}})}.$$

$${\tilde{f}}_t\equiv {\mathrm{argmin}}_{f\in \mathcal{F}}\frac{1}{n_t}\sum _{T_i=t}\left(\frac{\tilde{P}(T=t\mid {\mathbf{r}}_i)}{\tilde{P}(T=t\mid {\mathbf{r}}_i,{{\mathbf{c}}^{\text{T}}}_i)}\right)\mathcal{L}(y_i,f({\mathbf{r}}_i)),$$

$${\tilde{\beta }}_0^{\text{D1}}$$

$${\tilde{\beta }}_1^{\text{D1}}$$

$${\widehat{\text{ATE}}}^{+}\equiv \frac{1}{N^{+}}\sum _{\{j\mid \tilde{B}({\mathbf{r}}_j)=1\}}\frac{t_j{y}_j}{\tilde{P}(T=1\mid {\mathbf{C}}^{\text{T, eval }}=c^{\text{T, eval}})}-\frac{1}{N^{+}}\sum _{\{j\mid \tilde{B}({\mathbf{r}}_j)=1\}}\frac{(1-t_j)y_j}{\tilde{P}(T=0\mid {\mathbf{C}}^{\text{T, eval }}=c^{\text{T, eval}})}$$

$${\widehat{\text{ATE}}}^{-}\equiv \frac{1}{N^{-}}\sum _{\{j\mid \tilde{B}({\mathbf{r}}_j)=0\}}\frac{t_j{y}_j}{\tilde{P}(T=1\mid {\mathbf{C}}^{\text{T, eval }}=c^{\text{T, eval}})}-\frac{1}{N^{-}}\sum _{\{j\mid \tilde{B}({\mathbf{r}}_j)=0\}}\frac{(1-t_j)y_j}{\tilde{P}(T=0\mid {\mathbf{C}}^{\text{T, eval }}=c^{\text{T, eval}})}$$

$$\widehat{\text{ABR}}\equiv \left(\frac{N^{+}}{N^{+}+N^{-}}\right){\widehat{\text{ATE}}}^{+}+\left(\frac{N^{-}}{N^{+}+N^{-}}\right)\left(-{\widehat{\text{ATE}}}^{-}\right),$$

$$X_i\sim \text{Uniform}(0,2),L_i\sim \text{Bernoulli}(0.5),G_i\sim \text{Normal}(0,1)$$

$$T_i\mid L_i$$

$$P(Y_i=1\mid X_i,L_i,T_i)$$

$$\Uppsi_{\mathbf{R}}\equiv\int_{\mathbf{C^{\text{T }}}}\left\{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}E(Y\mid\mathbf{C^{\text{T }}},\mathbf{R},T=1)}-{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}E(Y\mid\mathbf{C^{\text{T }}},\mathbf{R},T=0)}\right\}\>dP^{\text{theoretical}}(\mathbf{C^{\text{T }}}\mid\mathbf{R}),$$

$$\text{Uppsi}(T=1\mid \mathbf{R})=$$

$$\displaystyle=\int_{\mathbf{C^{\text{T }}}}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\min_{f\in\mathcal{F}}\left\{\int_{\mathcal{Y}}\left(y-f(\mathbf{R})\right)^{2}\>dP^{\text{actual}}(Y\mid\mathbf{C^{\text{T }}},\mathbf{R},T=1)\right\}}\>dP^{\text{theoretical}}(\mathbf{C^{\text{T }}}\mid\mathbf{R}),$$

$$=\underset{f\in \mathcal{F}}{\min }{\int }_{{\mathbf{C}}^{\text{T }}}{\int }_{\mathcal{Y}}{\left(y-f(\mathbf{R})\right)}^{2}d{P}^{\text{actual}}(Y\mid {\mathbf{C}}^{\text{T }},\mathbf{R},T=1)d{P}^{\text{theoretical}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R}).$$

$$\text{Uppsi}(T=1,\mathbf{R})\equiv$$

$$\displaystyle=\int_{\mathbf{R}}\int_{\mathbf{C^{\text{T }}}}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\min_{f\in\mathcal{F}}\bigg{\{}\int_{\mathcal{Y}}\left(y-f(\mathbf{R})\right)^{2}\>dP^{\text{actual}}(Y\mid\mathbf{C^{\text{T }}},\mathbf{R},T=1)\bigg{\}}}\>dP^{\text{theoretical}}(\mathbf{C^{\text{T }}}\mid\mathbf{R})\>dP^{\text{actual}}(\mathbf{R}\mid T=1)$$

$$={\int }_{\mathbf{R}}\underset{f\in \mathcal{F}}{\min }{\int }_{{\mathbf{C}}^{\text{T }}}{\int }_{\mathcal{Y}}{\left(y-f(\mathbf{R})\right)}^{2}d{P}^{\text{theoretical}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R})d{P}^{\text{actual}}(Y\mid {\mathbf{C}}^{\text{T }},\mathbf{R},T=1)d{P}^{\text{actual}}(\mathbf{R}\mid T=1)$$

$$={\int }_{\mathbf{R}}\underset{f\in \mathcal{F}}{\min }{\int }_{{\mathbf{C}}^{\text{T }}}{\int }_{\mathcal{Y}}{\left(y-f(\mathbf{R})\right)}^{2}d{P}^{\text{theoretical}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R})d{P}^{\text{actual}}(Y\mid {\mathbf{C}}^{\text{T }},\mathbf{R},T=1)d{P}^{\text{actual}}(\mathbf{R}\mid T=1)\frac{d{P}^{\text{actual}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R},T=1)}{d{P}^{\text{actual}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R},T=1)}$$

$$={\int }_{\mathbf{R}}\underset{f\in \mathcal{F}}{\min }{\int }_{{\mathbf{C}}^{\text{T }}}{\int }_{\mathcal{Y}}{\left(y-f(\mathbf{R})\right)}^2\frac{d{P}^{\text{theoretical}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R})}{d{P}^{\text{actual}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R},T=1)}d{P}^{\text{actual}}(Y,{\mathbf{C}}^{\text{T }},\mathbf{R}\mid T=1).$$

$$P^{\text{actual}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R},T=1)$$

$$=\frac{P^{\text{actual}}(T=1\mid {\mathbf{C}}^{\text{T }},\mathbf{R})P^{\text{actual}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R})P^{\text{actual}}(\mathbf{R})}{P^{\text{actual}}(T=1\mid \mathbf{R})P^{\text{actual}}(\mathbf{R})},$$

$$=\frac{P^{\text{actual}}(T=1\mid {\mathbf{C}}^{\text{T }},\mathbf{R})P^{\text{actual}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R})}{P^{\text{actual}}(T=1\mid \mathbf{R})},$$

$$\equiv \frac{\pi ({\mathbf{C}}^{\text{T }},\mathbf{R})P^{\text{actual}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R})}{\pi (\mathbf{R})},$$

$$\text{Uppsi}(T=1,\mathbf{R})$$

$$={\int }_{\mathbf{R}}\underset{f\in \mathcal{F}}{\min }{\int }_{{\mathbf{C}}^{\text{T }}}{\int }_{\mathcal{Y}}{\left(y-f(\mathbf{R})\right)}^2\frac{d{P}^{\text{theoretical}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R})\pi (\mathbf{R})}{\pi ({\mathbf{C}}^{\text{T }},\mathbf{R})d{P}^{\text{actual}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R})}d{P}^{\text{actual}}(Y,{\mathbf{C}}^{\text{T }},\mathbf{R}\mid T=1)$$

$$={\int }_{\mathbf{R}}\underset{f\in \mathcal{F}}{\min }{\int }_{{\mathbf{C}}^{\text{T }}}{\int }_{\mathcal{Y}}{\left(y-f(\mathbf{R})\right)}^2\frac{d{P}^{\text{theoretical}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R})}{d{P}^{\text{actual}}({\mathbf{C}}^{\text{T }}\mid \mathbf{R})}\cdot \frac{\pi (\mathbf{R})}{\pi ({\mathbf{C}}^{\text{T }},\mathbf{R})}d{P}^{\text{actual}}(Y,{\mathbf{C}}^{\text{T }},\mathbf{R}\mid T=1)$$

$$\text{Uppsi}(T=1,\mathbf{R})={\int }_{\mathbf{R}}\underset{f\in \mathcal{F}}{\min }{\int }_{{\mathbf{C}}^{\text{T }}}{\int }_{\mathcal{Y}}{\left(y-f(\mathbf{R})\right)}^2\frac{\pi (\mathbf{R})}{\pi ({\mathbf{C}}^{\text{T }},\mathbf{R})}d{P}^{\text{actual}}(Y,{\mathbf{C}}^{\text{T }},\mathbf{R}\mid T=1).$$

$$\underset{f\in \mathcal{F}}{\min }\frac{1}{N_1}\sum _{i=1}^{N}I(T_i=1)\frac{\pi ({\mathbf{R}}_i)}{\pi ({{\mathbf{C}}^{\text{T }}}_i,{\mathbf{R}}_i)}{\left[Y_i-f({\mathbf{R}}_i)\right]}^2,$$

$$\text{Uppsi}(T=0,\mathbf{R})$$

$$\displaystyle=\int_{\mathbf{R}}\int_{\mathbf{C^{\text{T }}}}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\min_{f\in\mathcal{F}}\bigg{\{}\int_{\mathcal{Y}}\left(y-f(\mathbf{R})\right)^{2}\>dP^{\text{actual}}(Y\mid\mathbf{C^{\text{T }}},\mathbf{R},T=0)\bigg{\}}}\>dP^{\text{theoretical}}(\mathbf{C^{\text{T }}}\mid\mathbf{R})\>dP^{\text{actual}}(\mathbf{R}\mid T=0)$$

$$={\int }_{\mathbf{R}}\underset{f\in \mathcal{F}}{\min }{\int }_{{\mathbf{C}}^{\text{T }}}{\int }_{\mathcal{Y}}{\left(y-f(\mathbf{R})\right)}^2\frac{1-\pi (\mathbf{R})}{1-\pi ({\mathbf{C}}^{\text{T }},\mathbf{R})}d{P}^{\text{actual}}(Y,{\mathbf{C}}^{\text{T }},\mathbf{R}\mid T=0),$$

$$\underset{f\in \mathcal{F}}{\min }\frac{1}{N_0}\sum _{i=1}^{N}I(T_i=0)\frac{1-\pi ({\mathbf{R}}_i)}{1-\pi ({{\mathbf{C}}^{\text{T }}}_i,{\mathbf{R}}_i)}{\left[Y_i-f({\mathbf{R}}_i)\right]}^2$$

Full text

1 Introduction

Precision medicine strives to leverage an individual’s specific characteristics to determine the most beneficial course of treatment for that individual. In this paper, we consider the practice of precision medicine in settings that present two particular sets of challenges: 1) data come from observational studies where individual characteristics may influence treatment assignment; and 2) the data at hand may measure individual characteristics that will not be available in future clinical settings. As we summarize in Table 2, there are statistical methods to account for the observational study design in 1) but they do not account for the subtlety regarding variable roles in 2), and they often lack a user-friendly software implementation that would allow practitioners to reliably apply them.

Some characteristics that affect treatment recommendation in a particular study may be unavailable in future clinical settings, while other characteristics that did not directly affect treatment assignment may nonetheless be informative about treatment response in future clinical settings. For example, some information that influences treatment assignment (e.g. prognosis) may be measured subjectively in a given observational study and not be measured in that same manner in clinical settings outside the scope of the study; such a variable directly affects treatment assignment in the observational study but cannot directly influence treatment recommendations in future clinical settings. On the other hand, some available characteristics may not be interpretable in a clinical setting if there is a lack of scientific knowledge about their role in a disease pathway (e.g. gene expression levels); such characteristics do not directly affect the treatment decisions of clinicians but nonetheless may be predictive of treatment benefit. Each of these variable types should be handled differently by statistical methods that require distinct prediction of treatment assignment and prediction of outcome, but past statistical methods have not pointed out the distinction.

We propose a principled framework (along with a user-friendly implementation in the R package DevTreatRules, available on CRAN) that appropriately handles these distinct variable when developing and evaluating a treatment rule based on data from an observational study where treatment is not independent of individual characteristics. A treatment rule is a function that recommends treatment based on individual characteristics; to be useful to practitioners, a treatment rule must appropriately handle individual characteristics as they are actually observed in clinical settings. Our framework yields a treatment rule that accounts for the clinically distinct roles of individual characteristics in an interpretable way by asking two simple questions – 1) Which characteristics will be measured in the same manner in future settings where the estimated rule would be applied? 2) Which characteristics potentially influence treatment assignment and thus must be accounted for as confounders? – that reflect clinical rather than statistical expertise. Our framework also places a clear emphasis on developing and evaluating the treatment rule on independent datasets which is also absent from other papers in the treatment rule literature, as discussed in Section 4.

2 Previous Work

The framework presented in this paper draws on two previous bodies of work: 1) literature on estimating treatment rules;and 2) tools from the causal inference literature that describe how to estimate the average treatment effect (ATE) in study designs where treatment is not randomized. There is notable overlap between these two topics even if the connection is not always made explicit in the literature; nearly all methods for estimating treatment effect that are appropriate in settings with non-randomized treatment rely on the inverse-probability-of-treatment weighting (IPW) approach to balance observed confounders across the treatment groups. The propensity score (the probability of treatment, conditional on individual characteristics) has played a vital role in facilitating reliable comparisons of outcomes across treatment groups in non-randomized studies at least back to Rosenbaum and Rubin, (1983). Austin, (2011) presents an excellent and accessible overview of the ATE derived from the causal inference literature and how it informs propensity-score-based estimation strategies for non-randomized treatment assignment (e.g. the IPW method).

In this section, we highlight previous work on estimating treatment rules and we defer discussion of the causal inference literature to Section 6.1, where we can more clearly show how it shapes the target parameter of interest that is critical to our framework.

There is active statistical research on methods to predict the most beneficial treatment option for a specific individual, in line with the objectives of precision medicine. An existing approach generally belongs to one of two categories, either indirect or direct.

Indirect approaches (e.g. Kang et al., (2014); Cai et al., (2011); Lu et al., (2013); McKeague and Qian, (2014); Ciarleglio et al., (2015)) typically assume structure on the regression function linking the conditional mean outcome to individual outcomes and a treatment indicator. The ideal treatment assignment for a individual is then inferred by predicting his or her expected outcomes across possible values of the treatment variable using the regression model, and choosing the treatment option with the most desirable predicted outcome (e.g. a larger mean time until relapse or a smaller probability of 5-year relapse). In Lipkovich et al., (2017), indirect approaches are given the label “global outcome modeling”. As detailed in Section 6, our framework is an indirect approach.

On the other hand, direct methods are motivated by optimizing performance of the treatment rule itself in a population of interest rather then optimizing the accuracy of predicted outcomes and then using each individual-level predicted outcome to decide which treatment to assign (which would be the case with an indirect method). That is, these methods directly estimate an optimal treatment rule rather than prioritizing estimation of the expected outcome conditional on individual characteristics and treatment assignment (i.e. a regression function), and then taking the additional step of assigning an individual to the treatment with the most desirable predicted outcome.

Direct methods generally seek the treatment rule within a particular class of allowable rules that maximizes an estimate of clinical benefit, and have the appealing motivation of not being vulnerable to mis-specification of the regression function predicting outcome based on individual characteristics (though direct methods must still make other modeling assumptions to which they are sensitive). Lipkovich et al., (2017) present an outstanding recent survey of statistical methods for estimating treatment rules in RCTs; direct methods are part of that review’s “optimal treatment regimes” category. Zhang et al., (2015) pre-specify that the treatment rule must be a nested sequence of “if, then” statements with one or two individual characteristics involved in each statement. Zhang et al., 2012b take a similar approach but instead use a regression model to dictate the rule’s structure.

Other direct methods make less restrictive assumptions about the form of the treatment rule at the cost of interpretability of the rule. Zhao et al., (2012) propose outcome-weighted learning (OWL), which defines the optimal treatment rule as the solution to a weighted classification problem whose solution can be approximated using a modified form of support vector machines (SVM), a well-established statistical learning tool (Hastie et al.,, 2008). As noted by Zhang et al., 2012a , OWL can be viewed as part of a general weighted-classification framework for estimating a treatment rule, so any classification method that accommodates observation weights (e.g. classification trees or penalized regression (Hastie et al.,, 2008)) is a viable alternative to SVM in the estimation stage of OWL. The interactions-based procedure from Chen et al., (2017) based on the earlier work of Tian et al., (2014) is a recent example of a promising direct approach to estimating treatment rules with even greater flexibility.

Importantly, as mentioned in Lipkovich et al., (2017), some existing direct and indirect methods for estimating treatment rules are adaptable to observational study designs where treatment assignment is not independent of individual characteristics by accommodating the IPW approach (Austin,, 2011), which re-weights observations by the inverse of their estimated propensity scores so clinically observed confounders are roughly balanced between the treatment groups. The IPW adjustment thus allows for a sensible direct comparison of mean re-weighted outcomes across treatment groups that is reflective of the underlying population where the rule would be applied in the future. There is no consensus among researchers on whether the direct or indirect approach to estimating treatment rules yields superior results.

3 Categorization of Individual Characteristics

We propose partitioning each individual characteristic collected in an observational study (aside from outcome and treatment variables) into the four clinically distinct categories shown in Table 1: $C^{\text{TI }}$, $C^{\text{TN }}$, $C^{\text{NI }}$, or $C^{\text{NN }}$, where the abbreviations in each superscript tell us whether each characteristic (C) affects treatment assignment in the current study (TN), is expected to be observed in independent studies (NI), both (TI), or neither (NN). This helps make explicit (as is not done in previous work) that only the $C^{\text{TI }}$and $C^{\text{NI }}$variables are viable candidates for inclusion in a treatment rule because the remaining $C^{\text{TN }}$and $C^{\text{NN }}$will not be available in future clinical situations where a proposed treatment rule would be implemented. Here is a brief example of how each variable type may present itself in a hypothetical observational study:

Example of $C^{\text{TI }}$: Age. In the study population, physicians might have been more likely to recommend treatment to older individuals. In independent clinical settings, we are confident that age can be reliably measured on the same scale. 2. 2.

Example of $C^{\text{TN }}$: Center-specific measure of prognosis. In the study population, clinicians may have been more likely to recommend treatment to individuals with a poor prognosis as estimated by a center-specific set of guidelines (e.g. a hospital’s standard rule-of-thumb procedure based on their specific doctors’ prior experiences with individuals from the population). In independent clinical settings taking place in different centers, clinicians would not estimate an individual’s prognosis with that same center-specific approach. 3. 3.

Example of $C^{\text{NI }}$: Gene expression levels. In the study population, individuals’ gene expression levels may have been measured but clinicians did not consider the information when recommending treatment due to a lack of scientific knowledge about the genes’ roles in disease progression and response to treatment. In independent clinical settings, gene expression levels might still be reliably measured and would be eligible to inform treatment decisions if a newly developed treatment rule suggests their importance or if other scientific knowledge becomes available in the interim. 4. 4.

Discussion of $C^{\text{NN }}$: $C^{\text{NN }}$consists of variables in a dataset that are believed to have no role in influencing treatment and cannot be reliably collected in future clinical settings; the variables in $C^{\text{NN }}$are not of interest to development or evaluation of treatment rules. An example could be a study ID variable.

It will be useful in later sections to define: $\mathbf{C^{\text{I }}}\equiv(C^{\text{TI }},C^{\text{NI }})$ as all observed individual characteristics that are also expected to be observed in independent clinical settings. Also, we define $\mathbf{R}$ as a subset of the variables contained in $\mathbf{C^{\text{I }}}$ that the researcher believes may affect response to treatment and thus are viable candidates for a treatment rule. We also define $\mathbf{C^{\text{T }}}\equiv(C^{\text{TI }},C^{\text{TN }})$ as all individual characteristics that potentially affect treatment in the current observational study. In the BuildRule() and EvaluateRule() functions from DevTreatRules, users are required to provide the characteristics in $\mathbf{C^{\text{T }}}$ with the argument names.influencing.treatment and the characteristics in $\mathbf{R}$ with the argument names.influencing.rule. That is, individual characteristics are never entered into the package without first being categorized by the user.

4 Gaps in the Literature

To help situate this paper in the literature, Table 2 presents a selection of direct and indirect methods for estimating treatment rules and whether each satisfies five criteria. As seen in Table 2, none of these selected previous methods (and, to our knowledge, no other method in the literature) distinguishes between observed individual characteristics using the clinically meaningful categorization of whether they will be available in future clinical settings (in which case they are sensible candidates for building the treatment rule) or are not expected to be available in future settings (in which case they should not be used to build the rule). The lack of available R packages implementing previous approaches is also re-enforced by Table 2.

We believe that our work fills a gap in the literature by providing practitioners with a principled approach to appropriately classify variable types as in Table 1 and by providing the user-friendly DevTreatRules, which actually requires users to make their clinically informed variable categorizations when they apply the work to real data (using the arguments names.influencing.treatment and names.influencing.rule in its BuildRule() and EvaluateRule() functions). In addition, although it is not considered in the table, none of those previous methods explicitly integrates data-splitting to ensure that the stages of development/evaluation (or development/validation/evaluation) of a treatment rule are conducted on independent datasets to yield a trustworthy estimate of the rule’s population impact; we emphasize data-splitting throughout this paper and in the formalized procedure in Section 6.2.

5 Motivating Example

Suppose an observational study recruits individuals with a particular type of cancer from a single hospital and collects baseline information at the time of recruitment. Further suppose that, for each individual, this observational dataset measures: months until relapse ($Y$); an indicator of receiving standard-of-care ($T=0$) or additional chemotherapy in addition to the standard-of-care option ($T=1$); age; a measure of day-to-day life functioning; and expression levels of $p$ genes. For simplicity, we assume that clinicians decide to treat individuals based only on age and day-to-day life functioning. In contrast, the gene expression levels are unobservable by clinicians due to high cost and are also uninterpretable due to a lack of scientific knowledge about the genes’ roles in disease progression and response to treatment. In this illustrative example, we assume no unmeasured confounding.111In practice, however, gene expression levels may serve as surrogates for unobserved confounders. For instance, it is estimated that 20%-30% of breast cancer cases consist of an over-expression of human epidermal growth factor receptor type 2 (HER2), which can be targeted by additional treatment that inhibits HER2 generation (Joensuu et al.,, 2006; Hudis,, 2007). For types of cancer with unknown subtypes the gene expressions conducted in an observational study may be associated with (but not explicitly define, as with HER2) disease subtypes that may be more or less vulnerable to disruption by the treatment mechanism.

In practice, there are two alternative scores used to an individual’s day-to-day functioning. One is the Eastern Cooperative Oncology Group (ECOG) Performance Status, a grade ranging from 0 to 5 where a lower grade represents fewer restrictions in day-to-day life (Oken et al.,, 1982). The other is the Karnofsky Performance Status (KPS), an 11-point scale ranging from 0 to 100 where a lower value represents more day-to-day restrictions (Karnofsky,, 1949). As an added complication, we suppose that in this particular observational study day-to-day functioning is measured using the ECOG score, but that KPS is used instead of ECOG score in future clinical settings where we would like to apply our developed treatment rule in the future.222Although in practice there are mappings between ECOG and KPS, for illustrative purposes here we treat them as variables representing the same individual characteristic using qualitatively distinct scales (e.g. higher scores mean lower quality of life for KPS but a higher quality of life for ECOG) and suppose they are non-conformable. In contrast, we assume an individual’s age and $p$ gene expression measurements will be reliably collected in future clinical settings.

Figure 1 shows our hypothesized data-generating mechanism. We have the following variable types in this example: C^{\text{NI }}$$= (gene1, …, genep), C^{\text{TI }}$$= age, C^{\text{TN }}$$= ECOG, \mathbf{C^{\text{T }}}$$= (age, ECOG), and \mathbf{R}$$=(gene1, …, genep). The KPS variable is a potential confounder only in independent clinical settings (but not the current dataset) and as such KPS plays an important role in evaluation of the rule that we discuss in Section 6.3.

We are interested in building and evaluating a classifier, for future individuals, that indicates their optimal treatment based on their gene expression values. However, in constructing this classifier, we must account for a) the confounding influence of age and day-to-day functioning on the relationship between treatment and months until relapse; and b) the distinct approaches to measuring day-to-day functioning across different clinical settings (based on use of ECOG or KPS).

6 Method

Here, we provide some theoretical justification for how the individual characteristics $C^{\text{TN }}$, $C^{\text{TI }}$, and $C^{\text{NI }}$should be used to develop a treatment rule and evaluate the rule’s benefit. We will refer to this method as the split-regression approach to developing a treatment rule. As in Section 5, we will interpret $T$ as an indicator of receiving a new treatment in addition to standard-of-care ($T=1$) or standard-of-care alone ($T=0$) and outcome $Y$ as months until relapse.

6.1 Estimating the Rule

6.1.1 The Target Parameter

Our goals are 1) to identify a subset of the population – in terms of $\mathbf{R}$, the researcher-chosen subset of individual characteristics that are valid candidates for inclusion in a treatment rule – we expect to benefit from treatment and 2) to estimate the extent of benefit in this subpopulation.

Our estimation strategy begins with a foundational parameter of interest, $E[Y^{1}-Y^{0}]$, often called the average treatment effect (ATE) in the causal inference literature, where in potential outcomes notation, $Y^{1}$ is the months until relapse that would have been observed had an individual received treatment and $Y^{0}$ is the months until relapse that would have been observed had the individual received standard-of-care (see Kennedy, (2015) for an excellent review). So the ATE gives the average difference in months until relapse when an individual in the population of interest receives treatment instead of standard-of-care. The ATE is identifiable under three assumptions: 1) consistency: $T=t$ implies $Y=Y^{t}$; 2) no unmeasured confounding: conditional on observing $\mathbf{C^{\text{T }}}$, $T$ is independent of $Y^{t}$; 3) positivity: if $P(\mathbf{C^{\text{T }}}>0)$ then $P(T=t\mid\mathbf{C^{\text{T }}}=\mathbf{c^{\text{T}}})>0$, for $t=0,1$.

As originally developed by Robins, (1986), if the consistency, no unmeasured confounding, and positivity assumptions hold, then the ATE is equivalent to

[TABLE]

which is known as a “g-computation” formula in the causal inference literature. To estimate whether treatment offers a superior outcome for an individual with characteristics $\mathbf{R}=\mathbf{r}$, we simply consider the subgroup-specific ATE

[TABLE]

The subgroup-specific ATE in (2) is also a parameter of interest in Cai et al., (2011), the previous method that we believe is most similar to the one we discuss in this section. The practical limitations of Cai et al., (2011) on which this paper expands are that it was designed only for the RCT setting, did not explicitly distinguish between individual characteristics $\mathbf{C^{\text{T }}}$ and $\mathbf{C^{\text{I }}}$, was not implemented in an R package, and did not share code to guide users in implementation.

To estimate (2) in the setting of non-randomized treatment assignment, the g-computation formula in (3) similarly changes to

[TABLE]

which estimates the average treatment effect for individuals with characteristics $\mathbf{R}=\mathbf{r}$ that will be observable in future clinical settings.

From (3), we see that we should treat the individuals defined by $\Omega^{+}=\left\{\mathbf{r}\mid\psi(\mathbf{r})>C\right\}$ for $C\geq 0$. If $\mathbf{R}$ were a set of gene expression levels and $C=0$, for example, then $\Omega^{+}$ would be the subset of gene expression levels for which treatment increases the expected number of months until relapse. The expected improvement in months until relapse among the treated subpopulation is $\int_{\mathbf{r}\in\Omega^{+}}\psi(\mathbf{r})dP(\mathbf{r}).$

6.1.2 Estimating the Target Parameter

One possible approach for estimating the modified ATE in (3) would be to separately estimate $E\left[Y\mid\mathbf{C^{\text{T }}},\mathbf{R},T=t\right]$ for $t=0,1$ and estimate $dP(\mathbf{C^{\text{T }}}|\mathbf{R})$, then plug these estimates into (3); however, this approach requires estimation of a conditional density function that is only practical when the individual characteristics in $\mathbf{R}$ are perhaps one or two categorical variables with very few levels while in other situations the approach would prescribe a very complicated and highly variable average (as described in greater detail in Chapter 3 of Varadhan and Seeger, (2013), for example).

Our goal is now to re-write our estimation target (3) as a simple minimization problem that does not involve estimation of the conditional density $dP(\mathbf{C^{\text{T }}}|\mathbf{R})$. We begin by re-stating (3) as $\psi(\mathbf{r})=f_{1}(\mathbf{r})-f_{0}(\mathbf{r})$, where, for $t=0,1$,

[TABLE]

We emphasize that (4) is exactly equivalent to (3), just with altered notation. By taking the perspective in (4) we only need to estimate $f_{t}(\mathbf{r})$ for $t\in\{0,1\}$ to obtain an estimate of the $\psi(\mathbf{r})$ in (3). As derived in the Appendix, it turns out that $f_{t}(\mathbf{r})$ can be written as the minimizer

[TABLE]

where $\mathcal{L}(y,f(\mathbf{r}))$ is any function of y and $\mathbf{r}$  such that its conditional expectation $E[y\mid f(\mathbf{r})]$ is the minimizer – which in practice we may think of as a “canonical” loss function such as squared-error loss for a continuous $y$ or logistic loss for a binary $y$ – and $\mathcal{F}$ is a function class in which the rule is known to lie (one possibility could be $\ell_{1}(P)$, the space of all absolutely integrable functions over $P$).

A natural weight function $w_{t}(\mathbf{r},l)$ that can be used turns out to be

[TABLE]

We note that the standard IPW observation weight would replace the numerator of (6) with a 1. In fact, the weight function in (6) is closely related to the “stabilized weights” proposed by Robins et al., (2000), which in this case would be $P(T=t)/P(T=t|\mathbf{C^{\text{T }}}=\mathbf{c^{\text{T}}},\mathbf{R}=\mathbf{r})$. The weight function implied by Robins et al., (2000) differs from (6) by the latter’s extra conditioning on $\mathbf{R}=\mathbf{r}$, the subset of individual characteristics that are potentially informative inputs to the treatment rule.

The derivation in (5) implies a natural estimate of $f_{t}$ that is not complicated by the dimensionality of $\mathbf{R}$: the minimizer of the weighted sample average over individuals in the $T=t$ group

[TABLE]

where $n_{t}$ is the number of individuals in the group $T=t$, $\tilde{P}(T=t\mid\mathbf{r}_{i})$ is an estimate of $P(T=t\mid\mathbf{R}=\mathbf{r})$, and $\tilde{P}(T=t|\mathbf{r}_{i},\mathbf{c^{\text{T}}}_{i})$ is an estimate of $P(T=t\mid\mathbf{R}=\mathbf{r}_{i},\mathbf{C^{\text{T }}}=\mathbf{c^{\text{T}}}_{i})$. We note that the estimate $\tilde{f}_{t}$ is only reasonable if $\mathcal{F}$ is a suitably constrained class (e.g. a class with smoothness constraints, like a Sobolev class or class with bounded total variation (van de Geer,, 2000)).

However, the formula (7) suggests another approach for estimation of $f_{t}$: Instead of necessarily solving a formal empirical minimization as in (7), one might use any predictive modeling method that accommodates observation weights (e.g. generalized linear models, lasso, ridge regression, boosted trees, neural nets, and many others). All of these methods can be written as, either exactly or approximately, minimizing a weighted least-squares-like loss over a, potentially complicated, function class. Many of the methods indicated in Table 2 as being flexible in fact only support penalized or non-penalized weighted regression; our framework supports much more flexibility by moving beyond regression-based methods.

Now that we can estimate $\tilde{\psi}(\mathbf{r})=\tilde{f}_{1}(\mathbf{r})-\tilde{f}_{0}(\mathbf{r}),$ we can simply form the treatment rule as $\tilde{B}(\mathbf{r})\equiv I\left[\tilde{f}_{1}(\mathbf{r})-\tilde{f}_{0}(\mathbf{r})>0\right]$, which recommends treatment to an individual with characteristics $\mathbf{R}=\mathbf{r}$ if the estimated months until relapse is higher under $T=1$ than $T=0$.

6.2 The Recipe

The work in Sections 6.1 yields the following procedure applied to a development dataset (D1), which must be independent of the evaluation dataset (D2); we use the superscripts D1 and D2 to emphasize the dataset on which the accompanying estimate is formed.

Use the scientific knowledge underlying D1 to partition observed individual characteristics into the four categories presented in Table 1: $C^{\text{TI }}$, $C^{\text{TN }}$, $C^{\text{NI }}$, and $C^{\text{NN }}$. Also form $\mathbf{C^{\text{T }}}=(C^{\text{TI }},C^{\text{TN }})$, $\mathbf{C^{\text{I }}}=(C^{\text{TI }},C^{\text{NI }})$, and form the potential inputs for the treatment rule $\mathbf{R}\subseteq\mathbf{C^{\text{I }}}$. 2. 2.

For observations $i=1,\ldots,n$ on D1:

1. (a)

Choose a prediction method and estimate the propensity scores $\tilde{P}^{\text{D1}}(T=1\mid\mathbf{R}=\mathbf{r}_{i})$ and $\tilde{P}^{\text{D1}}(T=1\mid\mathbf{R}=\mathbf{r}_{i},\mathbf{C^{\text{T }}}=\mathbf{c^{\text{T}}}_{i})$.333It can be useful to truncate estimated propensity scores so they are not too close to 0 or 1 (which can lead to very large observation weights); a default setting in our software implementation truncates estimated propensity scores to stay between 0.05 and 0.95, but this choice can be overwritten by the user. 2. (b)

Compute the weights $\tilde{W}_{t}(\mathbf{R}=\mathbf{r}_{i},\mathbf{C^{\text{T }}}=\mathbf{c^{\text{T}}}_{i})=\frac{\tilde{P}^{\text{D1}}(T=t\mid\mathbf{R}=\mathbf{r}_{i})}{\tilde{P}^{\text{D1}}(T=t\mid\mathbf{R}=\mathbf{r}_{i},\mathbf{C^{\text{T }}}=\mathbf{c^{\text{T}}}_{i})}$, for $t=0,1$. 3. (c)

Choose a prediction method that accommodates observation weights (e.g. generalized linear regression, lasso, boosted trees, and many others) and estimate $\tilde{f}_{0}$ and $\tilde{f}_{1}$ as suggested by (7). For example, weighted linear regression with a continuous response would yield, for observations $i=1,\ldots,n$ on D1,

[TABLE]

where $\tilde{W}_{t}(\mathbf{R}=\mathbf{r}_{i},\mathbf{C^{\text{T }}}=\mathbf{c^{\text{T}}}_{i})=\frac{\tilde{P}^{\text{D1}}(T=t\mid\mathbf{R}=\mathbf{r}_{i})}{\tilde{P}^{\text{D1}}(T=t\mid\mathbf{R}=\mathbf{r}_{i},\mathbf{C^{\text{T }}}=\mathbf{c^{\text{T}}}_{i})}$ for $t=0,1$ and where we define $\tilde{f}_{0}^{\text{D1}}(\mathbf{r})\equiv\mathbf{r}^{\top}\tilde{\beta}_{0}^{\text{D1}}$ and $\tilde{f}_{1}^{\text{D1}}(\mathbf{r})\equiv\mathbf{r}^{\top}\tilde{\beta}_{1}^{\text{D1}}$. 3. 3.

Form the treatment rule $\tilde{B}(\mathbf{r})\equiv I\left[\tilde{f}_{1}^{\text{D1}}(\mathbf{r})-\tilde{f}_{0}^{\text{D1}}(\mathbf{r})>0\right]$, where $I(\cdot)$ is the indicator function. 4. 4.

Use scientific knowledge underlying D2 to select the potential confounders $\mathbf{C^{\text{T, eval }}}$. 444Ideally, D1 and D2 would be datasets from separate observational studies where D2 independently samples from the population where future intervention would take place, but D1/D2 may also be a random partition of data from a single observational study; in the latter case, we will have $\mathbf{C^{\text{T, eval }}}=\mathbf{C^{\text{T }}}$. The DevTreatRules package supports developing and evaluating rules in either situation. 5. 5.

For observations $j=1,\ldots,m$ on D2:

1. (a)

Assign the recommended treatment with $\tilde{B}_{j}^{\text{D2}}\equiv\tilde{B}^{\text{D2}}(\mathbf{r}_{j})$. 2. (b)

As discussed next in Section 6.3, form the IPW-based estimators of the ATE in the test-positives and in the test-negatives using (10) and (11), respectively.

6.3 Evaluating the Rule

Now we define $\mathbf{C^{\text{T, eval }}}$ as the set of potential confounders of the association between treatment and response in the evaluation dataset. We note that if the development and evaluation datasets are partitions of the same observational dataset then the variables in $\mathbf{C^{\text{T, eval }}}$ will be identical to those in $\mathbf{C^{\text{T }}}$, but this need not be the case when development and evaluation are carried out using data from separate studies. For example, in the motivating example from Section 5, we had \mathbf{C^{\text{T }}}$$=(age, ECOG) but \mathbf{C^{\text{T, eval }}}$$= (age, KPS).

To evaluate the the developed rule $\tilde{B}$ using Section 6.2, we can use an estimate of the ATE in the test-positives population with the IPW-based estimator (see e.g. Austin, (2011))

[TABLE]

where $N^{+}$ is the number of test-positives. A small modification to (10) also estimates the effect of avoiding treatment among the test-negatives:

[TABLE]

where $N^{-}$ is the number of test-negatives. We note that we would expect (11) to be a negative number for a rule that accurately identifies the test-negatives.

We can also form a simple estimator of the average benefit of the rule (ABR) in the population from which our the evaluation dataset is a representative sample with

[TABLE]

a weighted average of the benefit of receiving treatment the test-positives from (10) and the benefit of avoiding treatment in the test-negatives from (11), respectively.

Alternatively one could estimate ATE with the estimator developed by Robins et al., (1994) that is sometimes called the doubly-robust or augmented analog of (10); Lunceford and Davidian, (2004) present an excellent derivation and simulation study comparing different estimators of the ATE. We chose to present only the IPW-based estimators for ease of exposition.

6.4 R Implementation

The R package DevTreatRules implements this paper’s split-regression approach. In particular, the functions SplitData(), BuildRule(), and EvaluateRule() handle, respectively: the development/evaluation partitioning of a dataset (or development/validation/evaluation partitioning if model selection is also performed); the development of the treatment rule (in dataset D1) as in steps 2-4 of Section 6.2; and the evaluation of the rule (in dataset D2) as in steps 5-6 of Section 6.2. The vignette accompanying DevTreatRules walks through an example of building and evaluating a treatment rule using the package, in a situation where model selection is also performed using the CompareRulesOnValidation() function.

7 Simulations

We simulate data for a development sample of size $n$ with the individual characteristics

[TABLE]

the treatment indicator

[TABLE]

and binary outcome

[TABLE]

for $i=1,\ldots,n$, where $\beta_{0,T=t},\beta_{1,T=t},\gamma_{T=t}\in\mathbb{R}$ for $t=0,1$. Using the notation from Section 1, we would categorize $\mathbf{C^{\text{T }}}=\mathbf{C^{\text{T, eval }}}=L$ and $\mathbf{C^{\text{I }}}=\mathbf{R}=(X,G)$.

Figure 2 depicts the relationship between $P(Y\mid X,L,T)$ and $X$, where the true response probability for the $L=1$ group is shown with circles and for the $L=0$ group with triangles. The standard-of-care group is shown in blue and the treatment group is shown in orange. As seen in Figure 2, the outcome ($Y$) is more likely under treatment than under standard-of-care for individuals with a value of $X$ about 1.3 in both the $L=1$ and $L=0$ groups. Thus, the optimal treatment rule would recommend treatment to an individuals with $X>1.3$. However, due to the confounding effect of $L$, an empirical average of the response-curves for each treatment group (solid blue line and dotted orange line in the left panel of Figure 2) incorrectly suggests that there is no subset of individuals for whom treatment makes the outcome more likely. On the other hand, the IPW approach (solid blue line and dotted orange line in the right panel of Figure 2) uses a weighted average of response-curves, where each weight is the IPW weight (based on $L$), to correctly identify the value of $X$ (about 1.3) where the response-curves cross.

In Table 3, we present the mean probability of the (desirable) outcome for a range of development set sample sizes, specifying logistic regression (with $\mathbf{R}$ as the predictors) for steps 2 in Section 6.2. The first row shows the mean outcome probability for rules built using the split-regression approach described in Section 6.2. The second row reports the mean outcome probability for a modification to the split-regression approach that “naively” uses the incorrect sample averaging shown in the left panel of Figure 2 (i.e. it uses identical observation weights in step 2 of Section 6.2).

The estimated outcome probabilities in the first three rows of Table 3 can be compared to the benchmark value of $0.574$ for the optimal rule (known from the data-generating mechanism) that perfectly assigns treatment to only those who benefit and withholds it from those who do not benefit. We also compare to a rule that recommends everyone receive treatment (which may be the prevailing policy when an available treatment is believed to be uniformly effective) and to a rule that recommends no one receive treatment (which might be the preferred strategy when the effectiveness of a proposed treatment has not yet been established).

Table 3 shows the advantage of using split-regression with IPW weights relative to the “naive” uniform weighting: The bias of the naive approach forces the estimation of the incorrect non-crossing response-curves shown with the solid lines in the left panel of Figure 2. In contrast, split-regression with the correct IPW weighting approaches the optimal treatment rule with large enough sample size because it is estimating the crossing response-curves in the right panel of Figure 2. The code used to perform the simulations is shared at github.com/jhroth/simulations-split-regression.

8 Data Example: WHI-OS

We also illustrate the split-regression approach and compare it to alternatives by applying the R package DevTreatRules to the Women’s Health Initiative Observational Study component (WHI-OS). A detailed description of the study design and summaries of baseline measurements for participants (postmenopausal women between the ages of 49 and 81 who were recruited for the WHI clinical trial but either declined to participate or were later deemed ineligible) are available in Langer et al., (2003). We also present summary tables of the variables we retained for analysis in the appendix. The GitHub page github.com/jhroth/data-examplesplit-regression contains the R code needed to go start-to-finish from loading the raw WHI-OS datasets (access to which requires additional permission) to replicating our estimates in the evaluation subset.

Briefly, we aim to build a treatment rule to assign baseline hormone therapy (HRT) – defined as currently using unopposed estrogen and/or estrogen plus progesterone at baseline – to postmenopausal women if it will increase a woman’s probability of remaining free of coronary heart disease (CHD) after 10 years and, in a separate analysis, if it will increase a woman’s probability of remaining free of breast cancer after 10 years. All variables included in our analysis besides the outcomes were measured at baseline. We used the “adjudicated” outcome variables in the WHI-OS as described in Curb et al., (2003).

We classified $4$ categorical variables as belonging to $C^{\text{TN }}$: education level, ethnicity, family income, and how each participant heard about the study. We identified $31$ self-reported variables to make up $C^{\text{TI }}$and we chose $\mathbf{R}=C^{\text{TI }}$, so there are $31$ candidates to be used as inputs in our treatment rule. We did not specify any variables as belonging to $C^{\text{NI }}$. In the Appendix, we present the complete list of these variables and simple summaries including counts of missing values. We intend for this to be primarily an illustrative example rather than a definitive claim about appropriate variable classifications in this study.

About 17.6% of the $94140$ observations in the initial dataset had a missing value of at least one variable in $\mathbf{C^{\text{T }}}$ or $\mathbf{R}$. Instead of conducting a complete-case analysis that would drop observations with missing values, we used an IPW-based adjustment for missingness (Seaman and White,, 2013). Our shared code shows the details of our adjustment for missingness using the additional.weights argument of the BuildRule() and EvaluateRule() functions in DevTreatRules.

The top-half of Table 4 presents estimated ATEs and ABR in the validation set for two split-regression specifications using the outcome of no breast cancer after 10 years: one that used ridge regression for the propensity score and lasso for the rule, and another specification that used logistic regression for both the propensity score and rule models. On the validation set, none of the split-regression specifications for the outcome of no CHD after 10 years appeared to be an improvement over the naive strategy of treating no one and thus none of those specifications were chosen. Since the estimated ATEs and ABR in Table 4 informed our model selection – in particular, we chose these models because their ABRs in the validation set were relatively high – they do not serve as trustworthy estimates of ATE and ABR in independent samples drawn from this population in the future.

The bottom-half of Table 4 presents the estimated ATEs and ABR in the evaluation set, which did not inform model selection. Unfortunately, the 95% CI for estimated ATE among the treated population contains [math] for both rules in the evaluation set; as a result, we do not find evidence that either rule has identified a subpopulation of individuals who appear to benefit from treatment.

9 Discussion

We outlined a principled approach to classify the roles of variables collected in an observational study into clinically meaningful categories and, using that knowledge, to develop a treatment rule along with a trustworthy estimate of the rule’s population benefit. Since this paper is intended to be a practical guide to help practitioners go from start to finish in estimating and evaluating treatment rules without getting bogged down by the onerous and error-prone tasks of coding the method from scratch in statistical software, we implemented our approach in the R package DevTreatRules and shared the code needed to reproduce our simulations and data example on GitHub.

In a simple simulation study, we saw the benefit of estimating a treatment rule using this paper’s preferred approach with IPW weighting compared to using uniform observational weights that ignore the observational study design. In the WHI-OS data example, we used the split-regression approach to develop a treatment rule that assigns baseline hormone therapy to postmenopausal women if it is expected to increase their probability of remaining free of coronary heart disease after 10 years or free of breast cancer after 10 years. With the 10-year CHD outcome, split-regression did not estimate a rule with a positive estimate of ATE in the treated subgroup on the validation set. With the 10-year breast cancer outcome, however, split-regression did identify a rule that recommends HRT to about 18% of women and, among this treated subpopulation, had an estimated 1.5 percentage-point decrease in the probability of breast cancer. However, this 1.5 percentage-point decrease lacks statistical significance (95% CI: $-0.027,0.068$) and, as a result, we would simply recommend not assigning HRT to any women in this population if the goal is to reduce 10-year breast cancer incidence or 10-year CHD incidence. This null finding is, unfortunately, often a typical result in the search for informative treatment rules with observational data.

One notable limitation of this work is that, while the algorithm outlined in Section 6.2 is fairly general, the accompanying R implementation does not have as much flexibility (e.g. the currently available estimation methods are linear/logistic regression and its lasso/ridge counterparts); in future work we hope to expand the package to support more estimation methods. Another limitation is that our data example is not publicly reproducible because we are unable to share the underlying WHI-OS dataset; however, we do share the code that will reproduce the data example for users who have access to the raw WHI-OS data files.

Appendix A Motivation for IPW Estimator of Average Treatment Effect

We are starting with

[TABLE]

and, focusing on the $T=1$ group for a fixed $\mathbf{R}$, we can write

[TABLE]

For random $\mathbf{R}$, we have

[TABLE]

Since

[TABLE]

we can rewrite (16) as

[TABLE]

If we assume $dP^{\text{theoretical}}(\mathbf{C^{\text{T }}}\mid\mathbf{R})=dP^{\text{actual}}(\mathbf{C^{\text{T }}}\mid\mathbf{R})$ (i.e. representative sampling of covariates, conditional on values of the biomarkers), then this becomes

[TABLE]

If we assume the propensity scores are known, then the plug-in estimator of (17) is

[TABLE]

where $N_{1}=\sum_{i=1}^{N}I(T_{i}=1)$. We note that (18) is weighted squared-error loss among $T=1$ group with weights $w_{1}(\mathbf{C^{\text{T }}},\mathbf{R})\equiv\pi(\mathbf{R})/\pi(\mathbf{C^{\text{T }}},\mathbf{R})$.

Similarly in the $T=0$ group, the target parameter is

[TABLE]

whose plug-in estimator, assuming the propensity scores are known, is

[TABLE]

where $N_{0}=\sum_{i=1}^{N}I(T_{i}=0)$. We note that (18) is weighted squared-error loss among $T=0$ group with weights $w_{0}(\mathbf{C^{\text{T }}},\mathbf{R})\equiv(1-\pi(\mathbf{R}))/(1-\pi(\mathbf{C^{\text{T }}},\mathbf{R}))$.

Appendix B Evaluating the Rule

As described in the paper, a foundational target parameter in our framework is the average treatment effect (ATE) in a subpopulation of individuals with characteristics $\mathbf{R}=\mathbf{r}$.

[TABLE]

which under the assumptions of consistency, no unmeasured confounding, and positivity, as detailed in Kennedy, (2015) can be re-written as

[TABLE]

which implies we should recommend treatment to individuals in the subgroup

[TABLE]

which is known as the test-positives group. If $\mathbf{R}$ were a set of gene expression levels, for example, then $\Omega^{+}$ would be the subset of gene expression levels for which treatment increases the expected number of months until relapse. Again, the expected improvement in months until relapse among this treated subpopulation is

[TABLE]

Similarly, treatment should not be recommended to individuals in the test-negatives subpopulation defined by

[TABLE]

and

[TABLE]

would yield the average increase in months until relapse for the subpopulation that avoids treatment.

To obtain a trustworthy estimate of how a developed treatment rule will perform for individuals seen in future clinical settings, it is absolutely essential that development of the rule and evaluation of the selected rule are performed independently. Chapter 7 (Model Assessment and Selection) of Hastie et al., (2008) provides an excellent discussion of this topic with regard to the predictive performance of model-based estimates. In brief, if the development and evaluation datasets coincide then our evaluation of the treatment rule’s benefit will be overly optimistic because it will reward the rule for incorrectly classifying noise in the development dataset as signal that would persist when we apply the rule to future individuals.

We can take the treatment rule $\tilde{B}(\mathbf{r})$ – which gives us a mapping from an individual’s particular characteristics $\mathbf{R}=\mathbf{r}$ to a treatment recommendation – that was estimated on a development dataset and apply it to the independent evaluation dataset to yield, for the $n_{\text{eval}}$ individuals in the evaluation dataset indexed by $i=1,\ldots,n_{\text{eval}}$,

[TABLE]

the test-positives subset of observations in the evaluation dataset, based on $\mathbf{R}$, who are expected to have more months until relapse under treatment than under standard-of-care. Similarly,

[TABLE]

yields the test-negatives subset of observations in the evaluation dataset, based on $\mathbf{R}$, who are expected to have fewer months until relapse under treatment than under standard-of-care.

Appendix C Data Example: Summary of Dataset