Estimated vs. User Defined ITR
The package allows to compare the performance of estimated ITRs with
user defined ITRs. The estimate_itr function takes the
following arguments:
fit |
a fitted object from the estimate_itr
function |
user_itr |
a function defined by users that returns a unit-level
continuous score for treatment assignment (we assume those that have
score less than 0 should not have treatment) |
data |
a data frame |
treatment |
a character string specifying the treatment variable in
the data |
outcome |
a character string specifying the outcome variable in
the data |
budget |
a numeric value specifying the maximum percentage of
population that can be treated under the budget constraint |
The function returns an object that contains the estimated GATE, ATE,
and AUPEC for the user defined ITR.
# estimate ITR
fit <- estimate_itr(
treatment = "T",
form = user_formula,
data = star_data,
algorithms = c("causal_forest"),
budget = 0.2,
split_ratio = 0.7)
#> Evaluate ITR under sample splitting ...
# user's own ITR
score_function <- function(data){
data %>%
mutate(score = case_when(
school_urban == 1 ~ 0.1, # inner-city
school_urban == 2 ~ 0.2, # suburban
school_urban == 3 ~ 0.4, # rural
school_urban == 4 ~ 0.3, # urban
)) %>%
pull(score) -> score
return(score)
}
# evalutate ITR
compare_itr <- evaluate_itr(
fit = fit,
user_itr = score_function,
data = star_data,
treatment = "T",
outcome = outcomes,
budget = 0.2)
#> Cannot compute PAPDp
# summarize estimates
summary(compare_itr)
#> -- PAPE ------------------------------------------------------------------------
#> estimate std.deviation algorithm statistic p.value
#> 1 -0.37 1.4 causal_forest -0.26 0.8
#> 2 0.00 0.0 user_itr NaN NaN
#>
#> -- PAPEp -----------------------------------------------------------------------
#> estimate std.deviation algorithm statistic p.value
#> 1 -0.11 1.08 causal_forest -0.1 0.92
#> 2 1.06 0.67 user_itr 1.6 0.11
#>
#> -- PAPDp -----------------------------------------------------------------------
#> Cannot compute PAPDp
#>
#> -- AUPEC -----------------------------------------------------------------------
#> estimate std.deviation algorithm statistic p.value
#> 1 -0.63 1.07 causal_forest -0.59 0.555
#> 2 -0.91 0.42 <NA> -2.19 0.028
#>
#> -- GATE ------------------------------------------------------------------------
#> estimate std.deviation algorithm group statistic p.value upper lower
#> 1 -34 108 causal_forest 1 -0.32 0.750 -212.2 143
#> 2 -223 107 causal_forest 2 -2.08 0.038 -399.9 -47
#> 3 181 108 causal_forest 3 1.68 0.093 3.8 359
#> 4 -28 108 causal_forest 4 -0.26 0.795 -205.4 149
#> 5 126 106 causal_forest 5 1.19 0.235 -48.7 301
#> 6 126 58 user_itr 1 2.16 0.031 30.2 222
#> 7 96 59 user_itr 2 1.62 0.105 -1.4 194
#> 8 -33 59 user_itr 3 -0.56 0.579 -129.7 64
#> 9 -139 59 user_itr 4 -2.36 0.018 -236.5 -42
#> 10 -32 59 user_itr 5 -0.54 0.589 -129.4 65
We plot the estimated Area Under the Prescriptive Effect Curve
(AUPEC) for the writing score across a range of budget constraints for
user defined ITR and estimated ITRs. The plot shows that the estimated
ITRs have better performance than the user defined ITR.
# plot the AUPEC
plot(compare_itr)
Existing Model vs. User-Defined Model
The package also allows to compare the performance of estimated ITRs
of existing ML packages with user defined models. The following code
shows an example using causal forest from the grf package
with sample splitting. The estimate_itr function takes the
following arguments:
treatment |
a character string specifying the treatment variable in
the data |
form |
a formula specifying the outcome and covariates |
data |
a data frame |
algorithms |
a character vector specifying the ML algorithms to be
used |
budget |
a numeric value specifying the maximum percentage of
population that can be treated under the budget constraint |
split_ratio |
a character string specifying the outcome variable in
the data |
user_model |
a character string specifying the user defined
model |
The user_model input should be a function that takes two
arguments: training_data and test_data. The
function will make use of the training_data to fit a model
and then use the test_data to estimate CATE or other
metrics of interest. It should also specify the way to get the ITR,
based on the estimated effects.
In the following example, we fit a linear model with sample splitting
and use the estimated CATE. We compute the ITR by assigning treatment to
those with positive CATE and no treatment to those with negative CATE.
The function user_model takes in the training data and test
data and return a list that contains (1) an ITR; (2) a fitted model; and
(3) a continuous score with the same length as the input data.
# user-defined model
user_model <- function(training_data, test_data){
# model fit on training data
fit <- train_model(training_data)
# estimate CATE on test data
compute_hatf <- function(fit, test_data){
score <- fit_predict(fit, test_data)
itr <- score_function(score)
return(list(itr = itr, score = score))
}
hatf <- compute_hatf(fit, test_data)
return(list(
itr = hatf$itr,
fit = fit,
score = hatf$score))
}
Note that the user defined model can be any model that returns a
unit-level continuous score for treatment assignment. It does not have
to be a linear model or model that estimate CATE. We can specify custom
functions in the train_model function and the
fit_predict function to compute the score. If the model
does not have a default predict function, we need to write
up a custom function with fit_predict.
# train model
train_model <- function(data){
fit <- lm(
Y ~ T*(cov1 + cov1 + cov3),
data = data)
return(fit)
}
# predict function
fit_predict <- function(fit, data){
# need to change this function if
# the model does not have a default predict function
score <- predict(fit, data)
return(score)
}
In addition, we can also choose any scoring rule that maps the score
to a binary indicator of treatment assignment.
# score function
score_function <- function(score){
itr <- (score >= 0) * 1
return(itr)
}
If split_ratio is specified, the function will split the
data into training and test data. The split_ratio should be
a numeric value between 0 and 1. Alternatively, if n_folds
is specified, the function will use the entire data to fit the user
defined model via cross-validation.
# estimate ITR
compare_fit <- estimate_itr(
treatment = "T",
form = user_formula,
data = star_data,
algorithms = c("causal_forest"),
budget = 0.2,
split_ratio = 0.7,
user_model = "user_model")
#> Evaluate ITR under sample splitting ...
# evaluate ITR
compare_est <- evaluate_itr(compare_fit)
# summarize estimates
summary(compare_est)
#> -- PAPE ------------------------------------------------------------------------
#> estimate std.deviation algorithm statistic p.value
#> 1 2.9e-01 1.2 causal_forest 0.25 0.81
#> 2 -2.2e-16 0.0 user_model -Inf 0.00
#>
#> -- PAPEp -----------------------------------------------------------------------
#> estimate std.deviation algorithm statistic p.value
#> 1 2.56 1.1 causal_forest 2.31 0.021
#> 2 0.62 1.0 user_model 0.61 0.542
#>
#> -- PAPDp -----------------------------------------------------------------------
#> estimate std.deviation algorithm statistic p.value
#> 1 1.9 1.6 causal_forest x user_model 1.2 0.23
#>
#> -- AUPEC -----------------------------------------------------------------------
#> estimate std.deviation algorithm statistic p.value
#> 1 0.20 0.93 causal_forest 0.22 0.83
#> 2 -0.11 0.81 user_model -0.13 0.89
#>
#> -- GATE ------------------------------------------------------------------------
#> estimate std.deviation algorithm group statistic p.value upper lower
#> 1 -72 108 causal_forest 1 -0.66 5.1e-01 -250 106
#> 2 107 107 causal_forest 2 1.00 3.2e-01 -69 284
#> 3 -94 108 causal_forest 3 -0.87 3.9e-01 -272 84
#> 4 57 108 causal_forest 4 0.53 6.0e-01 -121 235
#> 5 17 107 causal_forest 5 0.16 8.7e-01 -159 193
#> 6 -137 106 user_model 1 -1.30 1.9e-01 -311 37
#> 7 -571 104 user_model 2 -5.47 4.5e-08 -743 -400
#> 8 -529 105 user_model 3 -5.02 5.2e-07 -702 -355
#> 9 436 107 user_model 4 4.06 4.9e-05 260 613
#> 10 817 104 user_model 5 7.83 5.1e-15 645 989
plot(compare_est)
