【R/English】Permutation Feature Importance (PFI)


# Objective 
To explore Permutation Feature Importance (PFI) 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. 

# Permutation Feature Importance  

**Permutation Feature Importance (PFI)** *is defined to be the decrease in a model score when a single feature value is randomly shuffled 1. This procedure breaks the relationship between the feature and the target, thus the drop in the model score is indicative of how much the model depends on the feature.* -<a href = "https://scikit-learn.org/stable/modules/permutation_importance.html#:~:text=The%20permutation%20feature%20importance%20is,model%20depends%20on%20the%20feature.">scikit-learn </a> Here are simple 5 steps of PFI   

1. Predict the target with **ALL** explanatory variables and calculate prediction error, which is going to be the baseline   
2. Pick one explanatory variable and permeate/shuffle it on the debatable. Predict the target and calculate prediction error  
3. Calculate the difference of prediction errors from steps 1 and 2. Make the difference the **importance** of variables piked at step 2  
4. Repeat steps for all explanatory variables   
5. See the importance of all variables and analyze   

![PFI](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/2db6e750-498c-83d8-3ff5-f9aae2634043.png)


# 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)

data("BostonHousing")
df = BostonHousing
`%notin%` <- Negate(`%in%`)
```

### Obserview of the Dataset  

Here are overview of the dataset
![dataset_overvirew.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/f9f1ab11-c82d-74ed-c771-4e1c2147bb24.png)


```r
head(df)
```

![head().png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/9088d3aa-dddf-1e29-e7ee-fe496a41c242.png)


medv is our response variable, We predict this.  

```r
hist(df$medv,breaks = 20, main = 'Histgram of medv', xlab = 'medv ($ in 1,000)') 
```
![hist.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/89d857e8-00a6-03d1-6b7e-a78d7b66ea3b.png)


### Build a Model

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

```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
```
![fit1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/b8679be2-5fb9-c2ef-38e2-951db941437a.png)


### Predict *medv*  

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

metrics = metric_set(rmse, rsq)

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

|.metric|.estimator|.estimate|
|---|---|---| 
|rmse|standard|3.8857|  
|rsq|standard|0.8627|  



### Interpre Feature Importance  
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
  )
```
![expmainer.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/92779bfd-7475-898e-662c-30721c62fd49.png)


Use *model_parts* function to get PFI. Here you can see *rm* and *lstat* are the top 2 important variables to predict *medv*. The bark blue box chart show distribution of error loss since we calculate it multiple times.  

- *loss_function*: Evaluation metrics  
- *B*: # of shuffles   
- *type*: method of calculating importance "difference" or "ratio" are applicable 


```r
pfi = explainer %>% 
  model_parts(
    loss_function = loss_root_mean_square,
    B = 10,
    type = "difference"
  )

plot(pfi)
```

![pfi1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/f0b79717-97b4-51bf-a6dc-3a441f8df738.png)


#### FYI 

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

## Grouped Permutation Feature Importance (GPFI)  

If some explanatory variables are correlated with each other, PFI won't work well. Let's say $X0$ and $X1$ are correlated. While calculating the importance of $X0$, the model still uses $X1$ on prediction. The performance of the model would not decrease much because $X0$ and $X1$ are correlated. Thus, PFI will underestimate the importance of $X1$. In the plot below, rad; index of accessibility to radial highway, and tax;full-value property-tax rate per $10,000. In a situation like this, we should shuffle both variables together. In addition to that, we should use this GPFI when the variables are encoded by one-hot encoding. Or you can use it when you are dealing with data like latitudes and longitudes.  

```r
cor  <- df[,unlist(lapply(df, is.numeric))  ]
corrplot(cor(cor), method = 'number', order = 'alphabet')
```
![corrplot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/48f750e9-d815-0370-a833-c2cfac8738a5.png)



### Rad and Tax  

So let's run GPFI on our dataset. model_parts function have *variable_groups* method. It takes list objects. So make a list that contains name of explanatory variables in this case *rad* and *tax*[^1]. The source code of *feature_importance* is <a href = "https://rdrr.io/cran/DALEX/src/R/model_parts.R">here</a>. 

![GPFI.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/aa5d4a38-822f-fe63-eff1-93544028b7e1.png)



```r
# Make list
paired_var = list(c("tax","rad"))

# Male vector of explanatory variables Do not forget to take out your response variable  

all_vars = c(colnames(df)[colnames(df) %notin% c("tax","rad","medv")])

# Gather 
var_list = c(all_vars, paired_var)
                  
pfi = explainer %>% 
  model_parts(
    loss_function = loss_root_mean_square,
    B = 10,
    type = "difference",
    variable_groups = var_list
  )

# Gather explaniner object 
plot(pfi)
```

![pfi3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/6674843f-f229-fa9a-1da4-a981781ffced.png)



If you keep *tax* and *rad* in the plot, you can see that the importance of *tax* and *rad* are dispersed. 


```r
# Make list
paired_var = list(c("tax","rad"))

# Make vector of explanatory variables Do not forget to take out your response variable  

all_vars = c(colnames(df)[colnames(df) %notin% c("medv")])

# Gather 
var_list = c(all_vars, paired_var)
                  
pfi = explainer %>% 
  model_parts(
    loss_function = loss_root_mean_square,
    B = 10,
    type = "difference",
    variable_groups = var_list
  )

# Gather explaniner object 
plot(pfi)
```
![pfi2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/58dc6369-87e9-3b1f-843b-c076dfc2ad82.png)


# Conclution  

PFI and GPFI are very sufficient models to calculate the importance of explanatory variables in the model. On the other hand, PFI does not explain how each variable affects the prediction of the model. This could be done by Partial Dependence (PD). 



### References

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

https://dropout009.hatenablog.com/

https://scikit-learn.org/stable/modules/permutation_importance.html#:~:text=The%20permutation%20feature%20importance%20is,model%20depends%20on%20the%20feature. 

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

https://stats.stackexchange.com/questions/92419/relative-importance-of-a-set-of-predictors-in-a-random-forests-classification-in

https://rdrr.io/cran/DALEX/src/R/model_parts.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

[^1]: It may not be right to pair up tax and rad variables without decent causal inference.