Efficiency augmentation

To demonstrate how to use efficiency augmentation and the propensity score utilities available in the personalized package, we simulate data with two treatments. The treatment assignments are based on covariates and hence mimic an observational setting with no unmeasured confounders.

library(personalized)

In this simulation, the treatment assignment depends on covariates and hence we must model the propensity score \(\pi(x) = Pr(T = 1 | X = x)\). In this simulation we will assume that larger values of the outcome are better.

Propensity score utilities

Estimation of the propensity score is a fundamental aspect of the estimation of individualized treatment rules (ITRs). The personalized package offers support tools for construction of the propensity score function used by the fit.subgroup() function. The support is via the create.propensity.function() function. This tool allows for estimation of the propensity score in high dimensions via glmnet. In high dimensions it can be important to account for regularization bias via cross-fitting (https://arxiv.org/abs/1608.00060); the create.propensity.function() offers a cross-fitting approach for high-dimensional propensity score estimation. A basic usage of this function with cross-fitting (with 10 folds) is as follows:

prop.func can then be passed to fit.subgroup() as follows:

## family:    gaussian 
## loss:      sq_loss_lasso 
## method:    weighting 
## cutpoint:  0 
## propensity 
## function:  propensity.func 
## 
## benefit score: f(x), 
## Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint'
## 
## Average Outcomes:
##                 Recommended 0      Recommended 1
## Received 0 -10.5061 (n = 116) -17.5969 (n = 284)
## Received 1 -15.1773 (n = 371) -11.0091 (n = 229)
## 
## Treatment effects conditional on subgroups:
## Est of E[Y|T=0,Recom=0]-E[Y|T=/=0,Recom=0] 
##                           4.6712 (n = 487) 
## Est of E[Y|T=1,Recom=1]-E[Y|T=/=1,Recom=1] 
##                           6.5878 (n = 513) 
## 
## NOTE: The above average outcomes are biased estimates of
##       the expected outcomes conditional on subgroups. 
##       Use 'validate.subgroup()' to obtain unbiased estimates.
## 
## ---------------------------------------------------
## 
## Benefit score quantiles (f(X) for 1 vs 0): 
##       0%      25%      50%      75%     100% 
## -16.5151  -3.0487   0.2011   3.6675  14.8049 
## 
## ---------------------------------------------------
## 
## Summary of individual treatment effects: 
## E[Y|T=1, X] - E[Y|T=0, X]
## 
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -33.0302  -6.0975   0.4022   0.7844   7.3349  29.6098 
## 
## ---------------------------------------------------
## 
## 9 out of 50 interactions selected in total by the lasso (cross validation criterion).
## 
## The first estimate is the treatment main effect, which is always selected. 
## Any other variables selected represent treatment-covariate interactions.
## 
##            Trt1     V2      V3     V4     V11     V12     V21     V29
## Estimate 0.6624 0.6976 -0.3349 -0.003 -0.0212 -0.0435 -1.1087 -0.0102
##             V36    V41
## Estimate 0.2934 0.7411

Augmentation utilities

Efficiency in estimating ITRs can be improved by including an augmentation term. The optimal augmentation term generally is a function of the main effects of the full outcome regression model marginalized over the treatment. Especially in high dimensions, regularization bias can potentially have a negative impact on performance. Cross-fitting is again another reasonable approach to circumventing this issue. Augmentation functions can be constructed (with cross-fitting as an option) via the create.augmentation.function() function, which works similarly as the create.propensity.function() function. The basic usage is as follows:

aug.func can be used for augmentation by passing it to fit.subgroup() like:

## family:    gaussian 
## loss:      sq_loss_lasso 
## method:    weighting 
## cutpoint:  0 
## augmentation 
## function: augment.func 
## propensity 
## function:  propensity.func 
## 
## benefit score: f(x), 
## Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint'
## 
## Average Outcomes:
##                 Recommended 0      Recommended 1
## Received 0  -6.4022 (n = 135) -19.8831 (n = 265)
## Received 1 -17.8593 (n = 202) -11.3133 (n = 398)
## 
## Treatment effects conditional on subgroups:
## Est of E[Y|T=0,Recom=0]-E[Y|T=/=0,Recom=0] 
##                          11.4572 (n = 337) 
## Est of E[Y|T=1,Recom=1]-E[Y|T=/=1,Recom=1] 
##                           8.5698 (n = 663) 
## 
## NOTE: The above average outcomes are biased estimates of
##       the expected outcomes conditional on subgroups. 
##       Use 'validate.subgroup()' to obtain unbiased estimates.
## 
## ---------------------------------------------------
## 
## Benefit score quantiles (f(X) for 1 vs 0): 
##      0%     25%     50%     75%    100% 
## -6.3155 -0.5672  0.9020  2.6326  9.1404 
## 
## ---------------------------------------------------
## 
## Summary of individual treatment effects: 
## E[Y|T=1, X] - E[Y|T=0, X]
## 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -12.631  -1.134   1.804   1.886   5.265  18.281 
## 
## ---------------------------------------------------
## 
## 1 out of 50 interactions selected in total by the lasso (cross validation criterion).
## 
## The first estimate is the treatment main effect, which is always selected. 
## Any other variables selected represent treatment-covariate interactions.
## 
##            Trt1     V2
## Estimate 0.9793 0.7474

Comparing performance with augmentation

We first run the training/testing procedure to assess the performance of the non-augmented estimator:

## family:  gaussian 
## loss:    sq_loss_lasso 
## method:  weighting 
## 
## validation method:  training_test_replication 
## cutpoint:           0 
## replications:       5 
## 
## benefit score: f(x), 
## Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint'
## 
## Average Test Set Outcomes:
##                               Recommended 0
## Received 0 -14.3377 (SE = 4.1142, n = 29.6)
## Received 1 -13.2299 (SE = 1.3556, n = 93.2)
##                               Recommended 1
## Received 0 -15.5863 (SE = 3.2318, n = 67.8)
## Received 1 -10.9153 (SE = 4.6092, n = 59.4)
## 
## Treatment effects conditional on subgroups:
## Est of E[Y|T=0,Recom=0]-E[Y|T=/=0,Recom=0] 
##           -1.1079 (SE = 4.6909, n = 122.8) 
## Est of E[Y|T=1,Recom=1]-E[Y|T=/=1,Recom=1] 
##             4.671 (SE = 6.9031, n = 127.2) 
## 
## Est of 
## E[Y|Trt received = Trt recom] - E[Y|Trt received =/= Trt recom]:                     
## 1.9651 (SE = 5.1967)

Then we compare with the augmented estimator. Although this is based on just 5 replications, we can see that the augmented estimator is much better at discriminating between benefitting and non-benefitting patients, as evidenced by the large treatment effect among those predicted to benefit by the augmented estimator versus the smaller conditional treatment effect above.

## family:  gaussian 
## loss:    sq_loss_lasso 
## method:  weighting 
## 
## validation method:  training_test_replication 
## cutpoint:           0 
## replications:       5 
## 
## benefit score: f(x), 
## Trt recom = 1*I(f(x)>c)+0*I(f(x)<=c) where c is 'cutpoint'
## 
## Average Test Set Outcomes:
##                               Recommended 0
## Received 0 -10.2154 (SE = 4.8772, n = 39.8)
## Received 1 -16.7413 (SE = 2.5499, n = 54.6)
##                               Recommended 1
## Received 0 -18.1801 (SE = 2.9749, n = 59.8)
## Received 1 -11.1547 (SE = 1.9546, n = 95.8)
## 
## Treatment effects conditional on subgroups:
## Est of E[Y|T=0,Recom=0]-E[Y|T=/=0,Recom=0] 
##             6.5259 (SE = 6.4483, n = 94.4) 
## Est of E[Y|T=1,Recom=1]-E[Y|T=/=1,Recom=1] 
##            7.0254 (SE = 4.2785, n = 155.6) 
## 
## Est of 
## E[Y|Trt received = Trt recom] - E[Y|Trt received =/= Trt recom]:                     
## 6.3519 (SE = 2.8961)