# Quick Start Guide: Optimal Policy Trees with Categorical Treatment

*This is an R version of the corresponding OptimalTrees quick start guide.*

In this guide we will give a demonstration of how to use Optimal Policy Trees with categoric treatment options. For this example, we will use the Credit Approval dataset, where the task is to predict the approval outcome of credit card application. Because the features have been anonymized for confidentiality reasons, we will arbitrarily select one of the categoric variables `A1`

to be the treatment.

*Note: this case is not intended to serve as a practical application of policy trees, but rather to serve as an illustration of the training and evaluation process.*

First we load in the data and drop 37 rows with missing values:

```
df <- read.csv("crx.data", na.strings = "?", stringsAsFactors = T, header = F,
col.names = paste("A", as.character(seq(1, 16)), sep = ""))
df <- na.omit(df)
```

```
A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 A14 A15 A16
1 b 30.83 0.00 u g w v 1.25 t t 1 f g 202 0 +
2 a 58.67 4.46 u g q h 3.04 t t 6 f g 43 560 +
3 a 24.50 0.50 u g q h 1.50 t f 0 f g 280 824 +
[ reached 'max' / getOption("max.print") -- omitted 650 rows ]
```

Policy trees are trained using a features matrix/dataframe `X`

as usual and a rewards matrix that has one column for each potential treatment that contains the outcome for each sample under that treatment.

There are two ways to get this rewards matrix:

- in rare cases, the problem may have full information about the outcome associated with each treatment for each sample
- more commonly, we have observational data, and use this partial data to train models to estimate the outcome associated with each treatment

Refer to the documentation on data preparation for more information on the data format.

In this case, the dataset is observational, and so we will use RewardEstimation to estimate our rewards matrix.

## Reward Estimation

First, we split into training and testing:

```
X <- subset(df, select = -c(A1, A16))
treatments <- df$A1
outcomes <- df$A16 == "+"
split <- iai::split_data("policy_maximize", X, treatments, outcomes, seed = 123,
train_proportion = 0.5)
train_X <- split$train$X
train_treatments <- split$train$treatments
train_outcomes <- split$train$outcomes
test_X <- split$test$X
test_treatments <- split$test$treatments
test_outcomes <- split$test$outcomes
```

Note that we have used a training/test split of so that we save more data for testing to ensure high-quality reward estimation on the test set.

The treatment, `A1`

, is a categoric variable, so we follow the process for estimating rewards with categorical treatments.

Our outcome is binary, so we use a `categorical_classification_reward_estimator`

to estimate the approval outcome under each option with a doubly-robust reward estimation method using random forests to estimate both propensity scores and outcomes:

```
reward_lnr <- iai::categorical_classification_reward_estimator(
propensity_estimator = iai::random_forest_classifier(),
outcome_estimator = iai::random_forest_classifier(),
reward_estimator = "doubly_robust",
random_seed = 123,
)
train_rewards <- iai::fit_predict(
reward_lnr, train_X, train_treatments, train_outcomes,
propensity_score_criterion = "auc", outcome_score_criterion = "auc")
train_rewards$rewards
```

```
a b
1 1.1569342 0.9800000
2 0.4936364 1.0034353
3 0.4000000 1.2678540
4 0.6075000 1.0675597
5 1.2689371 0.7750000
6 0.3000000 1.0688216
7 0.6925000 1.0499854
8 0.0300000 1.1949500
9 2.6753089 0.5726993
10 0.3116667 1.8181659
11 1.4972260 0.9050000
12 1.5646903 0.9700000
13 0.5500000 1.2753524
14 1.1846420 0.8300000
15 1.5822727 0.9800000
16 1.4536463 0.8735714
17 0.8036364 1.0453690
18 0.7436364 1.0838822
19 2.3086239 0.8416667
20 0.6087500 1.0484107
21 0.5700000 1.1226228
22 0.8600000 1.0084987
23 0.8687500 1.0154274
24 0.6636364 1.0085941
25 1.2897594 0.8700000
26 0.8300000 1.0548315
27 0.1400000 1.3806278
28 0.6200000 1.1009086
29 0.4700000 1.2016472
30 0.5473333 1.0571861
[ reached 'max' / getOption("max.print") -- omitted 296 rows ]
```

```
train_rewards$score
```

```
$propensity
[1] 0.6089214
$outcome
$outcome$b
[1] 0.9225966
$outcome$a
[1] 0.9566708
```

We can see that the internal outcome estimation models have AUCs of 0.92 and 0.96, which gives us confidence that the reward estimates are of decent quality, and good to base our training on. The AUC for the propensity model is lower at 0.61, which is not particularly high, and suggests difficulty in reliably estimating the propensity. The doubly-robust estimation method should help to alleviate this problem, as it is designed to deliver good results if either propensity scores or outcomes are estimated well. However, we should always pay attention to these scores and proceed with caution if the estimation quality is low.

## Optimal Policy Trees

Now that we have a complete rewards matrix, we can train a tree to learn an optimal prescription policy that maximizes credit approval outcome. We will use a `grid_search`

to fit an `optimal_tree_policy_maximizer`

(note that if we were trying to minimize the outcomes, we would use `optimal_tree_policy_minimizer`

):

```
grid <- iai::grid_search(
iai::optimal_tree_policy_maximizer(
random_seed = 121,
max_categoric_levels_before_warning = 20,
),
max_depth = 1:5,
)
iai::fit(grid, train_X, train_rewards$rewards)
iai::get_learner(grid)
```

```
Fitted OptimalTreePolicyMaximizer:
1) Split: A6 in [d,e,ff,i,k,m,q,r]
2) Split: A11 < 3.5
3) Prescribe: b, 128 points, error 5.416
4) Prescribe: a, 37 points, error 4.886
5) Split: A8 < 1.438
6) Prescribe: a, 83 points, error 5.231
7) Prescribe: b, 78 points, error 5.022
```

The resulting tree recommends different values for `A1`

based on other characteristics: namely, the values of `A6`

, `A8`

, and `A11`

. The intensity of the color in each leaf shows the difference in quality between the best and second-best values. We can see that both `a`

and `b`

values are prescribed by the tree, implying each treatment is better suited to certain subgroups.

We can make treatment prescriptions using `predict`

:

```
iai::predict(grid, train_X)
```

```
[1] "a" "b" "b" "b" "a" "b" "b" "b" "b" "b" "a" "a" "b" "a" "b" "b" "a" "b" "b"
[20] "b" "a" "b" "b" "b" "a" "b" "b" "a" "a" "b" "b" "b" "a" "a" "b" "b" "b" "b"
[39] "b" "b" "a" "b" "b" "b" "b" "a" "b" "b" "a" "a" "b" "b" "a" "b" "b" "a" "b"
[58] "b" "b" "b"
[ reached getOption("max.print") -- omitted 266 entries ]
```

If we want more information about the relative performance of treatments for these points, we can predict the full treatment ranking with `predict_treatment_rank`

:

```
iai::predict_treatment_rank(grid, train_X)
```

```
[,1] [,2]
[1,] "a" "b"
[2,] "b" "a"
[3,] "b" "a"
[4,] "b" "a"
[5,] "a" "b"
[6,] "b" "a"
[7,] "b" "a"
[8,] "b" "a"
[9,] "b" "a"
[10,] "b" "a"
[11,] "a" "b"
[12,] "a" "b"
[13,] "b" "a"
[14,] "a" "b"
[15,] "b" "a"
[16,] "b" "a"
[17,] "a" "b"
[18,] "b" "a"
[19,] "b" "a"
[20,] "b" "a"
[21,] "a" "b"
[22,] "b" "a"
[23,] "b" "a"
[24,] "b" "a"
[25,] "a" "b"
[26,] "b" "a"
[27,] "b" "a"
[28,] "a" "b"
[29,] "a" "b"
[30,] "b" "a"
[ reached getOption("max.print") -- omitted 296 rows ]
```

To quantify the difference in performance behind the treatment rankings, we can use `predict_treatment_outcome`

to extract the estimated quality of each treatment for each point:

```
iai::predict_treatment_outcome(grid, train_X)
```

```
a b
1 0.9036333 0.7599537
2 0.4797862 0.7681123
3 0.1111449 0.3740325
4 0.1111449 0.3740325
5 0.5584478 0.3209701
6 0.1111449 0.3740325
7 0.4797862 0.7681123
8 0.4797862 0.7681123
9 0.1111449 0.3740325
10 0.1111449 0.3740325
11 0.9036333 0.7599537
12 0.9036333 0.7599537
13 0.1111449 0.3740325
14 0.5584478 0.3209701
15 0.4797862 0.7681123
16 0.4797862 0.7681123
17 0.9036333 0.7599537
18 0.4797862 0.7681123
19 0.1111449 0.3740325
20 0.4797862 0.7681123
21 0.5584478 0.3209701
22 0.4797862 0.7681123
23 0.4797862 0.7681123
24 0.4797862 0.7681123
25 0.9036333 0.7599537
26 0.4797862 0.7681123
27 0.1111449 0.3740325
28 0.5584478 0.3209701
29 0.5584478 0.3209701
30 0.4797862 0.7681123
[ reached 'max' / getOption("max.print") -- omitted 296 rows ]
```

## Evaluating Optimal Policy Trees

It is critical for a fair evaluation that we do not evaluate the quality of the policy using rewards from our existing reward estimator trained on the training set. This is to avoid any information from the training set leaking through to the out-of-sample evaluation.

Instead, what we need to do is to estimate a new set of rewards using only the test set, and evaluate the policy against these rewards:

```
test_rewards <- iai::fit_predict(
reward_lnr, test_X, test_treatments, test_outcomes, propensity_score_criterion = "auc", outcome_score_criterion = "auc")
test_rewards$rewards
```

```
a b
1 0.6027619 1.1714377
2 1.1994289 0.7066667
3 0.7400000 1.1499867
4 0.8175000 1.0664644
5 0.7155000 1.0679181
6 0.6233333 1.1766472
7 0.8200000 1.0171730
8 0.7014286 1.3139850
9 1.1087486 0.5662857
10 1.5083790 0.6450000
11 1.3935686 0.7800000
12 0.9542857 1.0130374
13 0.8700000 1.0399533
14 0.9000000 1.0104069
15 0.7700000 1.0877628
16 1.0471866 0.7500000
17 0.8066667 1.0557882
18 1.6507572 0.7400000
19 0.9000000 1.0866270
20 1.0000000 1.0529102
21 0.6100000 1.1498051
22 0.8600000 1.0820338
23 0.1466667 1.0968707
24 0.0500000 1.2888514
25 2.0438714 0.5386667
26 0.5900000 1.1883981
27 0.6425641 1.0793621
28 0.8200000 1.1932089
29 1.2323624 0.7083333
30 1.1865466 0.6650000
[ reached 'max' / getOption("max.print") -- omitted 297 rows ]
```

```
test_rewards$score
```

```
$propensity
[1] 0.613628
$outcome
$outcome$b
[1] 0.9246304
$outcome$a
[1] 0.9197715
```

We see the scores are similar to those on the training set, giving us confidence that the estimated rewards are a fair reflection of reality, and will serve as a good basis for evaluation. The AUC for the propensity model is again on the low side, but should be compensated for by means of the doubly-robust estimator.

We can now evaluate the quality using these new estimated rewards. First, we will calculate the average predicted credit approval probability under the treatments prescribed by the tree for the test set. To do this, we use `predict_outcomes`

which uses the model to make prescriptions and looks up the predicted outcomes under these prescriptions:

```
policy_outcomes <- iai::predict_outcomes(grid, test_X, test_rewards$rewards)
```

```
[1] 0.6027619 0.7066667 1.1499867 1.0664644 0.7155000 1.1766472
[7] 1.0171730 0.7014286 0.5662857 0.6450000 1.3935686 1.0130374
[13] 1.0399533 1.0104069 0.7700000 1.0471866 1.0557882 1.6507572
[19] 0.9000000 1.0529102 0.6100000 1.0820338 1.0968707 0.0500000
[25] 2.0438714 0.5900000 0.6425641 0.8200000 0.7083333 0.6650000
[31] 1.0499801 0.9800000 1.1012993 -0.1767049 -0.6512686 -1.0381208
[37] 0.6192308 -0.2816803 -0.2016309 -0.1904316 -0.1088713 -0.7966556
[43] -0.2980677 0.7800000 -0.9520040 -0.2536429 0.5800000 0.6791667
[49] 0.7016667 0.7049524 0.7800000 0.8600000 0.8866667 1.1272797
[55] 1.0012934 0.8400000 0.8800000 3.4358133 1.1467204 0.7864286
[ reached getOption("max.print") -- omitted 267 entries ]
```

We can then get the average estimated approval probability under our treatments:

```
mean(policy_outcomes)
```

`[1] 0.4071621`

We can compare this number to the average approval probability under the treatment assignments that were actually observed:

```
evaluate_reward_mean <- function(treatments, rewards) {
total <- 0.0
for (i in seq(treatments)) {
total <- total + rewards[i, treatments[i]]
}
total / length(treatments)
}
evaluate_reward_mean(test_treatments, test_rewards$rewards)
```

`[1] 0.4065159`

We see a small improvement over the real treatments.