{"id":1753,"date":"2024-09-02T10:16:54","date_gmt":"2024-09-02T08:16:54","guid":{"rendered":"https:\/\/lorentzen.ch\/?p=1753"},"modified":"2024-09-02T12:17:15","modified_gmt":"2024-09-02T10:17:15","slug":"explaining-a-causal-forest","status":"publish","type":"post","link":"https:\/\/lorentzen.ch\/index.php\/2024\/09\/02\/explaining-a-causal-forest\/","title":{"rendered":"Explaining a Causal Forest"},"content":{"rendered":"\n

We use a causal forest [1] to model the treatment effect in a randomized controlled clinical trial. Then, we explain this black-box model with usual explainability tools. These will reveal segments where the treatment works better or worse, just like a forest plot, but multivariately.<\/p>\n\n\n\n

Data<\/h2>\n\n\n\n

For illustration, we use patient-level data of a 2-arm trial of rectal indomethacin against placebo to prevent post-ERCP pancreatitis (602 patients) [2]. The dataset is available in the package {medicaldata}.<\/p>\n\n\n\n

The data is in fantastic shape, so we don\u2019t need to spend a lot of time with data preparation.<\/p>\n\n\n\n

    \n
  1. We integer encode factors.<\/li>\n\n\n\n
  2. We select meaningful features, basically those shown in the forest plot of [2] (Figure 4) without low-information features and without hospital.<\/li>\n<\/ol>\n\n\n\n

    The marginal estimate of the treatment effect is -0.078, i.e., indomethacin reduces the probability of post-ERCP pancreatitis by 7.8 percentage points. Our aim is to develop and interpret a model to see if this value is associated with certain covariates.<\/p>\n\n\n\n

    library(medicaldata)\nsuppressPackageStartupMessages(library(dplyr))\nlibrary(grf)          #  causal_forest()\nlibrary(ggplot2)\nlibrary(patchwork)    #  Combine ggplots\nlibrary(hstats)       #  Friedman's H, PDP\nlibrary(kernelshap)   #  General SHAP\nlibrary(shapviz)      #  SHAP plots\n\nW <- as.integer(indo_rct$rx) - 1L      # 0=placebo, 1=treatment\ntable(W)\n#   0   1 \n# 307 295\n\nY <- as.numeric(indo_rct$outcome) - 1  # Y=1: post-ERCP pancreatitis (bad)\nmean(Y)  # 0.1312292\n\nmean(Y[W == 1]) - mean(Y[W == 0])      # -0.07785568\n\nxvars <- c(\n  "age",         # Age in years\n  "male",        # Male (1=yes)\n  "pep",         # Previous post-ERCP pancreatitis (1=yes)\n  "recpanc",     # History of recurrent Pancreatitis (1=yes)\n  "type",        # Sphincter of oddi dysfunction type\/level (0=no, to 3=type 3)\n  "difcan",      # Cannulation of the papilla was difficult (1=yes)\n  "psphinc",     # Pancreatic sphincterotomy performed (1=yes)\n  "bsphinc",     # Biliary sphincterotomy performed (1=yes)\n  "pdstent",     # Pancreatic stent (1=yes)\n  "train"        # Trainee involved in stenting (1=yes)\n)\n\nX <- indo_rct |>\n  mutate_if(is.factor, function(v) as.integer(v) - 1L) |> \n  rename(male = gender) |> \n  select_at(xvars)\n\nhead(X)\n            \n# age  male   pep recpanc  type difcan psphinc bsphinc pdstent train\n#  26     0     0       1     1      0       0       0       0     1\n#  24     1     1       0     0      0       0       1       0     0\n#  57     0     0       0     2      0       0       0       0     0\n#  29     0     0       0     1      0       0       1       1     1\n#  38     0     1       0     1      0       1       1       1     1\n#  59     0     0       0     1      1       0       1       1     0\n            \nsummary(X)\n            \n#     age             male             pep            recpanc     \n# Min.   :19.00   Min.   :0.0000   Min.   :0.0000   Min.   :0.000  \n# 1st Qu.:35.00   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.000  \n# Median :45.00   Median :0.0000   Median :0.0000   Median :0.000  \n# Mean   :45.27   Mean   :0.2093   Mean   :0.1595   Mean   :0.299  \n# 3rd Qu.:54.00   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:1.000  \n# Max.   :90.00   Max.   :1.0000   Max.   :1.0000   Max.   :1.000  \n#      type           difcan          psphinc          bsphinc      \n# Min.   :0.000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  \n# 1st Qu.:1.000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  \n# Median :2.000   Median :0.0000   Median :1.0000   Median :1.0000  \n# Mean   :1.743   Mean   :0.2608   Mean   :0.5698   Mean   :0.5714  \n# 3rd Qu.:2.000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000  \n# Max.   :3.000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  \n#    pdstent           train       \n# Min.   :0.0000   Min.   :0.0000  \n# 1st Qu.:1.0000   1st Qu.:0.0000  \n# Median :1.0000   Median :0.0000  \n# Mean   :0.8239   Mean   :0.4701  \n# 3rd Qu.:1.0000   3rd Qu.:1.0000  \n# Max.   :1.0000   Max.   :1.0000  <\/pre><\/div>\n\n\n\n
    The model<\/h2>\n\n\n\n
    We use the {grf} package to fit a causal forest<\/em> [1], a tree-ensemble trying to estimate conditional average treatment effects (CATE) E[Y(1) – Y(0) | X = x]. As such, it can be used to study treatment effect inhomogeneity<\/em>.<\/p>\n\n\n\n
    In contrast to a typical random forest:<\/p>\n\n\n\n