Quick Start Guide: Optimal Policy Trees with Numeric Treatment

In this guide we will give a demonstration of how to use Optimal Policy Trees with numeric treatment options. For this example, we will use the auto-mpg dataset, where the task is usually to predict a car's fuel efficiency (in miles-per-gallon, or MPG) based on other characteristics of the car. To apply a prescriptive lens to this case, we will instead treat the amount of acceleration as a treatment that can be controlled, and try to find the value that optimizes the MPG for a given car.

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. For a real-world case study using similar techniques, see the grocery pricing case study.

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

using CSV, DataFrames
df = CSV.read("auto-mpg.csv", DataFrame, missingstring="?")
dropmissing!(df)
392×9 DataFrame
 Row │ mpg      cylinders  displacement  horsepower  weight  acceleration  mod ⋯
     │ Float64  Int64      Float64       Int64       Int64   Float64       Int ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │    18.0          8         307.0         130    3504          12.0      ⋯
   2 │    15.0          8         350.0         165    3693          11.5
   3 │    18.0          8         318.0         150    3436          11.0
   4 │    16.0          8         304.0         150    3433          12.0
   5 │    17.0          8         302.0         140    3449          10.5      ⋯
   6 │    15.0          8         429.0         198    4341          10.0
   7 │    14.0          8         454.0         220    4354           9.0
   8 │    14.0          8         440.0         215    4312           8.5
  ⋮  │    ⋮         ⋮           ⋮            ⋮         ⋮          ⋮            ⋱
 386 │    36.0          4         135.0          84    2370          13.0      ⋯
 387 │    27.0          4         151.0          90    2950          17.3
 388 │    27.0          4         140.0          86    2790          15.6
 389 │    44.0          4          97.0          52    2130          24.6
 390 │    32.0          4         135.0          84    2295          11.6      ⋯
 391 │    28.0          4         120.0          79    2625          18.6
 392 │    31.0          4         119.0          82    2720          19.4
                                                  3 columns and 377 rows omitted

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

Warning

Please refer to the Reward Estimation documentation for a detailed description on how to perform reward estimation for numeric treatments properly. For simplicity and to keep the focus on Optimal Policy Trees, this quick start guide does not cover tuning the reward estimation parameters, but in real problems this tuning is an important step.

First, we split into training and testing:

X = select(df, Not(["mpg", "acceleration", "car name"]))
treatments = df.acceleration
outcomes = df.mpg

(train_X, train_treatments, train_outcomes), (test_X, test_treatments, test_outcomes) =
    IAI.split_data(:policy_maximize, X, treatments, outcomes, seed=123, train_proportion=0.5)

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

The treatment, acceleration, is a numeric value, so we follow the process for estimating rewards with numeric treatments. We will consider prescribing acceleration values from 11 to 20, in steps of 3:

treatment_candidates = 11:3:20

The outcome is continuous, so we use a NumericRegressionRewardEstimator to estimate the MPG under our candidate acceleration values using a Random Forest model:

reward_lnr = IAI.NumericRegressionRewardEstimator(
    propensity_estimator=IAI.RandomForestRegressor(),
    outcome_estimator=IAI.RandomForestRegressor(),
    reward_estimator=:doubly_robust,
    propensity_min_value=0.1,
    random_seed=1,
)

train_predictions, train_reward_score = IAI.fit_predict!(
    reward_lnr, train_X, train_treatments, train_outcomes, treatment_candidates)
train_rewards = train_predictions[:reward]
196×4 DataFrame
 Row │ 11       14       17       20
     │ Float64  Float64  Float64  Float64
─────┼────────────────────────────────────
   1 │ 15.5192  14.3211  15.9095  15.7161
   2 │ 17.9843  17.0684  16.1826  16.6844
   3 │ 15.1993  11.2804  14.1719  15.8209
   4 │ 13.6143  11.0216  14.1719  15.8209
   5 │ 14.0478  11.2476  14.5093  16.1168
   6 │ 15.0572  15.5953  16.1991  16.6844
   7 │ 15.0314  13.373   14.793   16.0183
   8 │ 13.3217  15.0754  18.3585  18.2202
  ⋮  │    ⋮        ⋮        ⋮        ⋮
 190 │ 21.6221  23.0741  22.5122  26.9959
 191 │ 26.0113  36.0119  34.4337  30.4533
 192 │ 26.4589  34.6858  42.8009  29.7727
 193 │ 26.3356  37.3628  32.1039  30.9279
 194 │ 26.4589  39.36    34.6957  31.9937
 195 │ 19.6138  18.9977  35.9689  19.9677
 196 │ 24.0463  26.3079  26.1804  39.972
                          181 rows omitted
train_reward_score[:propensity]
Dict{String, Float64} with 4 entries:
  "17" => 0.120533
  "20" => 0.200617
  "11" => 0.287145
  "14" => 0.165352
train_reward_score[:outcome]
Dict{String, Float64} with 4 entries:
  "17" => 0.864211
  "20" => 0.250276
  "11" => 0.653976
  "14" => 0.833718

We can see that the R2 of the internal outcome models for each candidate treatment are between 0.25 and 0.86, whereas the propensity models have R2 between 0.12 and 0.29. The non-zero propensity scores tell us that there is a small amount of treatment assignment bias in this data (keeping in mind that we arbitrarily selected a feature to use as the treatment column). It seems we can predict the outcomes reasonably well, so we have a good base for the reward estimates.

Optimal Policy Trees

Now that we have a complete rewards matrix, we can train a tree to learn an optimal prescription policy that maximizes MPG. We will use a GridSearch to fit an OptimalTreePolicyMaximizer (note that if we were trying to minimize the outcomes, we would use OptimalTreePolicyMinimizer):

grid = IAI.GridSearch(
    IAI.OptimalTreePolicyMaximizer(
        random_seed=1,
        minbucket=15,
    ),
    max_depth=4:5,
)
IAI.fit!(grid, train_X, train_rewards)
Optimal Trees Visualization

The resulting tree recommends different accelerations based on the characteristics of the car. The intensity of the color in each leaf shows the difference in quality between the best and second-best acceleration values.

We can see that a variety of our candidate acceleration values are prescribed by the tree. For instance, older cars with higher horsepower receive the highest acceleration.

We can make treatment prescriptions using predict:

prescriptions = IAI.predict(grid, train_X)
196-element Vector{String}:
 "20"
 "20"
 "20"
 "20"
 "20"
 "20"
 "20"
 "20"
 "20"
 "20"
 ⋮
 "14"
 "14"
 "20"
 "14"
 "14"
 "14"
 "14"
 "14"
 "20"

The prescriptions are always returned as strings matching the column names of the input rewards matrix. In our case the treatments are numeric values, and if we want them in numeric form to use later we can convert them to numeric treatments using convert_treatments_to_numeric:

IAI.convert_treatments_to_numeric(prescriptions)
196-element Vector{Int64}:
 20
 20
 20
 20
 20
 20
 20
 20
 20
 20
  ⋮
 14
 14
 20
 14
 14
 14
 14
 14
 20

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:

rank = IAI.predict_treatment_rank(grid, train_X)
196×4 Matrix{String}:
 "20"  "11"  "17"  "14"
 "20"  "11"  "17"  "14"
 "20"  "11"  "17"  "14"
 "20"  "11"  "17"  "14"
 "20"  "11"  "17"  "14"
 "20"  "11"  "17"  "14"
 "20"  "11"  "17"  "14"
 "20"  "11"  "17"  "14"
 "20"  "11"  "17"  "14"
 "20"  "11"  "17"  "14"
 ⋮
 "14"  "17"  "20"  "11"
 "14"  "17"  "20"  "11"
 "20"  "17"  "11"  "14"
 "14"  "17"  "20"  "11"
 "14"  "17"  "20"  "11"
 "14"  "17"  "20"  "11"
 "14"  "17"  "20"  "11"
 "14"  "11"  "17"  "20"
 "20"  "17"  "11"  "14"

For each point in the data, this gives the treatments in order of effectiveness. As before, this are returned as strings, but we can convert the treatments to numeric values with convert_treatments_to_numeric:

IAI.convert_treatments_to_numeric(rank)
196×4 Matrix{Int64}:
 20  11  17  14
 20  11  17  14
 20  11  17  14
 20  11  17  14
 20  11  17  14
 20  11  17  14
 20  11  17  14
 20  11  17  14
 20  11  17  14
 20  11  17  14
  ⋮
 14  17  20  11
 14  17  20  11
 20  17  11  14
 14  17  20  11
 14  17  20  11
 14  17  20  11
 14  17  20  11
 14  11  17  20
 20  17  11  14

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)
196×4 DataFrame
 Row │ 11       14       17       20
     │ Float64  Float64  Float64  Float64
─────┼────────────────────────────────────
   1 │ 17.021   14.7505  16.9222  18.2866
   2 │ 17.021   14.7505  16.9222  18.2866
   3 │ 17.021   14.7505  16.9222  18.2866
   4 │ 17.021   14.7505  16.9222  18.2866
   5 │ 17.021   14.7505  16.9222  18.2866
   6 │ 17.021   14.7505  16.9222  18.2866
   7 │ 17.021   14.7505  16.9222  18.2866
   8 │ 17.021   14.7505  16.9222  18.2866
  ⋮  │    ⋮        ⋮        ⋮        ⋮
 190 │ 22.1578  21.3878  22.3794  27.5812
 191 │ 26.3143  31.3847  29.363   28.3271
 192 │ 26.3143  31.3847  29.363   28.3271
 193 │ 26.3143  31.3847  29.363   28.3271
 194 │ 26.3143  31.3847  29.363   28.3271
 195 │ 18.6563  19.3162  18.4529  17.3021
 196 │ 22.1578  21.3878  22.3794  27.5812
                          181 rows omitted

We can also extract the standard errors of the outcome estimates with predict_treatment_outcome_standard_error:

IAI.predict_treatment_outcome_standard_error(grid, train_X)
196×4 DataFrame
 Row │ 11        14        17        20
     │ Float64   Float64   Float64   Float64
─────┼────────────────────────────────────────
   1 │ 0.80353   1.06319   0.770092  0.710896
   2 │ 0.80353   1.06319   0.770092  0.710896
   3 │ 0.80353   1.06319   0.770092  0.710896
   4 │ 0.80353   1.06319   0.770092  0.710896
   5 │ 0.80353   1.06319   0.770092  0.710896
   6 │ 0.80353   1.06319   0.770092  0.710896
   7 │ 0.80353   1.06319   0.770092  0.710896
   8 │ 0.80353   1.06319   0.770092  0.710896
  ⋮  │    ⋮         ⋮         ⋮         ⋮
 190 │ 0.513389  0.882271  1.02996   1.77642
 191 │ 0.948884  1.26322   0.739374  0.448817
 192 │ 0.948884  1.26322   0.739374  0.448817
 193 │ 0.948884  1.26322   0.739374  0.448817
 194 │ 0.948884  1.26322   0.739374  0.448817
 195 │ 0.23138   0.407594  0.625738  0.466373
 196 │ 0.513389  0.882271  1.02996   1.77642
                              181 rows omitted

These standard errors can be combined with the outcome estimates to construct confidence intervals in the usual way.

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_predictions, test_reward_score = IAI.fit_predict!(
    reward_lnr, test_X, test_treatments, test_outcomes, treatment_candidates)
test_rewards = test_predictions[:reward]
196×4 DataFrame
 Row │ 11       14       17       20
     │ Float64  Float64  Float64  Float64
─────┼────────────────────────────────────
   1 │ 18.8265  17.9006  15.2741  16.5961
   2 │ 18.5837  17.3373  15.3651  16.8543
   3 │ 16.8317  17.7591  15.4516  17.111
   4 │ 14.8371  12.4136  14.3317  16.1539
   5 │ 15.8857  13.5001  14.0596  16.2479
   6 │ 14.9192  16.8369  14.9     16.2122
   7 │ 28.6542  23.7261  23.4555  23.7447
   8 │ 28.7612  26.9357  25.2422  24.9468
  ⋮  │    ⋮        ⋮        ⋮        ⋮
 190 │ 30.2826  38.5807  28.1592  27.4065
 191 │ 27.5855  46.7469  31.719   25.8422
 192 │ 28.4471  27.5669  26.3109  25.6221
 193 │ 27.1386  22.2712  26.9411  27.5965
 194 │ 30.6563  35.0577  35.2216  30.77
 195 │ 70.4636  32.0715  31.719   25.7882
 196 │ 29.963   31.6384  25.1198  30.191
                          181 rows omitted
test_reward_score[:propensity]
Dict{String, Float64} with 4 entries:
  "17" => 0.158079
  "20" => 0.267284
  "11" => 0.393047
  "14" => 0.217413
test_reward_score[:outcome]
Dict{String, Float64} with 4 entries:
  "17" => 0.705786
  "20" => 0.770381
  "11" => 0.470776
  "14" => 0.886658

We see the internal models of our test reward estimator have scores that are similar to those on the test set. As with the training set, this gives us confidence that the estimated rewards are a fair reflection of reality, and will serve as a good basis for evaluation.

We can now evaluate the quality using these new estimated rewards. First, we will calculate the average predicted MPG 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)
196-element Vector{Float64}:
 16.59605211
 16.85433253
 17.11098041
 16.15388911
 16.24788885
 16.21219682
 23.74468928
 25.2421742
 30.65150663
 24.82638225
  ⋮
 19.79697195
 24.61787132
 27.40654131
 46.74687468
 25.62213452
 27.59645482
 35.05765523
 32.07152412
 30.19096294

We can then get the average estimated MPG under our treatments:

using Statistics
mean(policy_outcomes)
23.93172935

We can compare this number to the average estimated MPG under a baseline policy that assigns the same treatment to all observations regardless of their features. Let us compare to a policy that assigns 17 everywhere:

mean(test_rewards[:, "17"])
23.32036357

We can see that the personalization offered by the tree policy indeed improves upon a baseline where each observation receives the same treatment.