【R/English】Individual Conditional Expectation (ICE)

# Objective
To explore Individual Conditional Expectation (ICE) and run it in R script. This article is based on information in <a href = "https://www.amazon.co.jp/exec/obidos/ASIN/429712226X/hatena-blog-22/">「機械学習を解釈する技術 ~Techniques for Interpreting Machine Learning~」</a>by <a href = "https://dropout009.hatenablog.com/">Mitsunosuke Morishita</a>. In this book, the author does not go through all the methods by R, so I decided to make a brief note with an R script.  

# Individual Conditional Expectation (ICE)    

**Individual Conditional Expectation (ICE)** is a method of visualizing and understanding the relationship between the response variable and individual observation. The technique of calculating ICE is pretty simple if you are familiar with calculating <a href = "https://qiita.com/Shu25/items/e9c191c02ecd56730c7a">Partial Dependence (PD)</a>. PD is just a mean of ICE. PD is calculated by changing the value of an explanatory variable while keeping track of the prediction. For example, if we increase the unit of explanatory variable $X0$, and the response variable also increases, these variables have a positive relationship. In PD, we average out all observation orders to see the overview of the relationship, but in ICE, we plot all observations individually.  


![ICE1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/b091fb41-a71d-93d1-baeb-123b5b04c5f5.png)


### Formula for ICE  

ICE formula looks like a formula for PD. For your review, PD is  

$$
\large \hat{PD_j}(x_j) = \displaystyle \frac{1}{N} \sum_{i = 1}^{N} \hat{f}(x_j,x_{i,-j})
$$

Your set of explanatory variable as $X = (X_0,...,X_J)$, and define your trained model as $\hat{f}(X)$. Your target explanatory variable is $X_j$, and $X$ without $X_j$ is $X_{-j} = (X_0,...,X_{j-1},X_{j+1},...,X_J)$. Actual observation of $X_j$ at $i$th observation is defined as $x_{j,i}$. So actual observation without $x_{j,i}$ is $x_{i,-j} = (x_{i,0},...,x_{i,j-1},x_{i,j+1},...,x_{i,J})$.   
Explanatory variables are $X_j = x_j$, $\hat{PD_j}(x_j)$.   

For ICE, we take out the  averaging out part $\displaystyle \frac{1}{N} \sum_{i = 1}^{N}$ and add little *i* for individual observation.  

$$
\large \hat{ICE_{i,j}}(x_j) = \hat{f}(x_j,x_{i,-j})
$$



# Why? (The Limitation of PD)   

Let's say you are predicting customer satisfaction at a restaurant. The restaurant has a variety of customers of all ages. One of your explanatory variables is the amount of food per dish (g). You assume for a younger customer, would have higher satisfaction with more food on the dish. On the other hand, older people might not be interested in the amount of food. In the situation like that and if you only look at PD, you would miss the relationship between individual observation and the response variable. 


# Execution with Real Data   

Now, let's see how to run PFI with actual <a href = "http://math.furman.edu/~dcs/courses/math47/R/library/mlbench/html/BostonHousing.html">dataset</a>. 


### Get Dataset  

```r
# Set up
library(mlbench)
library(tidymodels)
library(DALEX)
library(ranger)
library(Rcpp)
library(corrplot)
library(ggplot2)
library(gridExtra)

data("BostonHousing")
df = BostonHousing
`%notin%` <- Negate(`%in%`)
set.seed(3)
```

### Overview of the Dataset  

Here is overview of the dataset

![dataset_overvirew.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/0c22fd86-df67-7dce-1101-46832c1345ea.png)





### Build a Model

We won't cover building a model in this article. I used XGBoost.  

```r
split = initial_split(df, 0.8)
train = training(split)
test = testing(split)

model = rand_forest(trees = 100, min_n = 1, mtry = 13) %>% 
  set_engine(engine = "ranger", seed(25)) %>% 
  set_mode("regression")

fit = model %>% 
  fit(medv ~., data=train)
fit
```


### Predict medv  

```r
result = test %>% 
  select(medv) %>% 
  bind_cols(predict(fit, test))

metrics = metric_set(rmse, rsq)

result %>% 
  metrics(medv, .pred)
```

# Interpre ICE   

Use the function *explain* to create an explainer object that helps us to interpret the model. 

```r
explainer = fit %>% 
  explain(
    data = test %>% select(-medv),
    y = test$medv
  )
```


Use *predict_profile* function to get ICE plot. The source code of *predict_profile* is <a href = "https://rdrr.io/cran/DALEX/src/R/model_profile.R">here</a>. You can designate which plot you like to plot by giving the variables method a vector of variable names in the plot function. (The blue dots are actual prediction and variable)    

```r
ice = explainer %>% 
  predict_profile(
    new_observation = test
  )
plot(ice,variables = c("lstat", "rm", "dis", "crim"), title = "ICE for lstat, rm, dis, crim")
```

![ice_plots.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/aef5b655-8f46-21a1-fe25-1a2169a7757b.png)

#### FYI 

|Method|Function|
|---|---|
|Permutation Feature Importance(PFI)|*model_parts()*|  
|Partial Dependence(PD)|*model_profile()*|  
|Individual Conditional Expectation(ICE)|*predict_profile()*|  
|SHAP|*predict_parts()*|  

# rm   

```r
plot(ice,variables = c("rm"), title = "ICE for rm")
```
![ice_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/423a1939-3bf2-9d44-db6c-25f5ce73c7d9.png)


If you look at the outlier on rm plot, as other observations follow the path, there is one unique observation with an opposite relationship with prediction. 


```r
cor  <- df[,unlist(lapply(df, is.numeric))]
pairs(cor, pch = 19, upper.panel = NULL, panel=panel.smooth, col = dplyr::case_when(row(cor) == 373 ~"red",T ~"grey"))
```
![pair_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/99e898c7-7e03-631f-9847-53baf3344b0a.png)



```r
# Merge prediction and test data.  
test_pred = cbind(test, explainer$y_hat)
# Get outlier
outlier = test_pred[(test_pred$`explainer$y_hat`>35)&(test_pred$rm<6),]
outlier
```
|index|crim|zn|indus|chas|nox|rm|age|dis|rad|tax|ptratio|b|lstat|medv|y_hat|  
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|  
373|8.26725|0|18.1|1|0.668|5.875|89.6|1.1296|24|666|20.2|347.88|8.88|50|40.461|   

With its extremely high rad, tax, and low dis, we can assume the houses are the center of the city. Maybe there ain't many houses with more than 6 rooms, so our model could not do a good job on that (extrapolation). or As the room number increases, the size of each room get smaller and smaller, so the house price would decrease due to the decrease of the house value.  


# Conclution  

While PD could not visualize individual relationships, ICE is a great way to visualize individual relationships to the response variable. It is a substance of PD, so if you understand PD, I bet ICE was straightforward to you. 

### References


https://www.amazon.co.jp/exec/obidos/ASIN/429712226X/hatena-blog-22/

https://dropout009.hatenablog.com/

https://qiita.com/Shu25/items/e9c191c02ecd56730c7a

http://math.furman.edu/~dcs/courses/math47/R/library/mlbench/html/BostonHousing.html

https://rdrr.io/cran/DALEX/src/R/model_profile.R

## Methods of Interpreting Machine Learning Qiita Links  

https://qiita.com/Shu25/items/446d30d88ac64919aa77

https://qiita.com/Shu25/items/e9c191c02ecd56730c7a

https://qiita.com/Shu25/items/f6320d8612de42119d6e

https://qiita.com/Shu25/items/344c852830ed88557aea