DALEX

book：https://ema.drwhy.ai/introduction.html

Which variables contribute to the selected prediction

加性归因的分解图

在试图理解模型对单个观察的预测时，最常见的问题可能是：哪些变量对这个结果的贡献最大？没有单一的最佳方法可以用来回答这个问题。我们介绍了分解 (BD) 图，它提供了一种可能的解决方案。这些图可用于呈现“变量归因”，即将模型的预测分解为可归因于不同解释变量的贡献。

首先看几幅图，这几幅图是理解分解图的关键：

预测的分布.

标记为“all data”的行中的小提琴图总结了数据集中所有 2207 个观测值的预测分布。红点表示预测的平均值，可以理解为对所有数据而言，期望的预测概率为0.235 .

为了评估单个解释变量对该特定单实例预测的贡献，固定连续变量值时模型预测的变化。图中面板 A 中标记为“age=8”的行中的小提琴图总结了当年龄变量取值“8（也就是将数据集中所有数据的年龄特征设置为8） 岁”时获得的预测分布，红点表示预测的均值，可以理解为当年龄取8的时候，模型对于数据的期望预测概率。“class=1st”行中的小提琴图用变量值描述了值age和class分别设置为“8 years”和“1st时候预测的分布和平均值。最后一行则表示这一条样本的的预测概率。

灰色细线显示了当特定解释变量的值被行名称中指示的值替换时，不同个体的预测变化。例如，第一行和第二行之间的线条表示将年龄变量的值固定为“8 岁”对不同的个体有不同的影响。对于某些人（很可能是 8 岁的乘客），模型的预测根本没有改变。对于其他人，预测值增加（可能对于 8 岁以上的乘客）或减少（最有可能对于 8 岁以下的乘客）。

B 和 C 中所示的简化图可能会引起人们的兴趣。请注意，在面板 C 中，标记为“截距”的行显示了整个数据集预测的总体平均值 (0.235)。连续的行表示通过固定特定解释变量的值引起的平均预测的变化。正变化用绿色条表示，而负差异用红色条表示。最后一行标记为“预测”，包含总体平均值和变化的总和。

这个图总结了特定解释变量对模型预测的影响。

从图中可以看到，泰坦尼克号数据集的随机森林模型的平均预测等于 23.5%。这是泰坦尼克号上所有人的平均生存预测概率。请注意，这不是幸存者的百分比，而是平均模型预测。因此，对于不同的模型，我们很可能会获得不同的平均值。

预测等于 42.2%。它远高于平均预测。对这一预测影响最大的两个解释变量是班级（值为“1st”）和年龄（值为 8）。通过固定这两个变量的值，我们将平均预测值提高了 35.6 个百分点。但是，性别将预测的生存概率降低了约 8.3 个百分点。

还值得一提的是，对于包含交互的模型，归因于变量的预测部分取决于解释变量值的设置顺序。

我们来看一个例子:

第一步构建模型，然后创建一个解释器：

library(DALEXtra)

## Loading required package: DALEX

## Welcome to DALEX (version: 2.4.2).
## Find examples and detailed introduction at: http://ema.drwhy.ai/

library(DALEX)
library(randomForest)

## randomForest 4.7-1.1

## Type rfNews() to see new features/changes/bug fixes.

titanic_imputed$survived <- as.factor(titanic_imputed$survived)
aps_lm_model4 <- randomForest(survived ~., data = titanic_imputed)
aps_lm_explainer4 <- DALEX::explain(aps_lm_model4, data = titanic_imputed, label = "survived")

## Preparation of a new explainer is initiated
##   -> model label       :  survived 
##   -> data              :  2207  rows  8  cols 
##   -> target variable   :  not specified! (  WARNING  )
##   -> predict function  :  yhat.randomForest  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package randomForest , ver. 4.7.1.1 , task classification (  default  ) 
##   -> model_info        :  Model info detected classification task but 'y' is a NULL .  (  WARNING  )
##   -> model_info        :  By deafult classification tasks supports only numercical 'y' parameter. 
##   -> model_info        :  Consider changing to numerical vector with 0 and 1 values.
##   -> model_info        :  Otherwise I will not be able to calculate residuals or loss function.
##   -> predicted values  :  numerical, min =  0 , mean =  0.2355061 , max =  1  
##   -> residual function :  difference between y and yhat (  default  )
##   A new explainer has been created!

aps_lm_explainer4

## Model label:  survived 
## Model class:  randomForest.formula,randomForest 
## Data head  :
##   gender age class    embarked  fare sibsp parch survived
## 1   male  42   3rd Southampton  7.11     0     0        0
## 2   male  13   3rd Southampton 20.05     0     2        0

解释器对象允许统一访问预测模型，而不管它们的内部结构如何。有了这个对象，我们就可以进行模型分析了。

DALEX::predict_parts()函数将模型预测分解为可归因于各个变量的部分。在最简单的调用中，该函数需要三个参数：

1. explainer- 解释器对象，使用函数创建DALEX::explain()；
2. new_observation- 需要解释的观察结果；它应该是一个数据框，其结构与用于拟合模型的数据集的结构相匹配；
3. ype- 变量归因的计算方法；可能的方法是”break_down”，可选包括”shap”、“oscillations”和”break_down_interactions”

bd_lm <- predict_parts(explainer = aps_lm_explainer4,
                 new_observation = titanic_imputed[1,],
                            type = "break_down")
bd_lm

##                                  contribution
## survived: intercept                     0.236
## survived: gender = male                -0.127
## survived: class = 3rd                  -0.009
## survived: age = 42                     -0.031
## survived: fare = 7.11                  -0.048
## survived: embarked = Southampton       -0.003
## survived: sibsp = 0                    -0.001
## survived: parch = 0                    -0.010
## survived: survived = 0                  0.000
## survived: prediction                    0.006

plot(bd_lm)

其他常用参数还包括：

1. order- 字符（列名）或整数（列索引）的向量，指定用于计算可变重要性度量的解释变量的顺序；
2. keep_distributions- 一个逻辑值（FALSE默认）；如果TRUE，则有关预测的条件分布的附加诊断信息存储在结果对象中，并且可以使用通用plot()函数绘制。

bd_lm <- predict_parts(explainer = aps_lm_explainer4,
                 new_observation = titanic_imputed[1,],
                            type = "break_down",
                 keep_distributions = TRUE,
                 order = c("class", "age", "gender", "fare", 
                             "parch", "sibsp", "embarked")
                 )
bd_lm

##                                  contribution
## survived: intercept                     0.236
## survived: class = 3rd                  -0.063
## survived: age = 42                     -0.045
## survived: gender = male                -0.059
## survived: fare = 7.11                  -0.048
## survived: parch = 0                    -0.008
## survived: sibsp = 0                    -0.005
## survived: embarked = Southampton       -0.002
## survived: prediction                    0.006

plot(bd_lm, plot_distributions = TRUE)

## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
## ℹ The deprecated feature was likely used in the iBreakDown package.
##   Please report the issue at
##   <]8;;https://github.com/ModelOriented/iBreakDown/issueshttps://github.com/ModelOriented/iBreakDown/issues]8;;>.

相互作用分解图

对于某些模型，例如具有交互作用的模型，相互作用分解图该算法识别变量对之间的相互作用，并在构建分解 (BD) 图时将它们考虑在内。

相互作用（偏离可加性）意味着解释变量的影响取决于其他变量的值，交互作用使解释变量相对于模型预测的重要性的评估变得复杂。

bd_lm <- predict_parts(explainer = aps_lm_explainer4,
                 new_observation = titanic_imputed[1,],
                            type = "break_down_interactions")
bd_lm

##                                  contribution
## survived: intercept                     0.236
## survived: gender = male                -0.127
## survived: class = 3rd                  -0.009
## survived: age = 42                     -0.031
## survived: fare = 7.11                  -0.048
## survived: embarked = Southampton       -0.003
## survived: sibsp = 0                    -0.001
## survived: parch = 0                    -0.010
## survived: survived = 0                  0.000
## survived: prediction                    0.006

plot(bd_lm)

平均归因的 Shapley 加法解释 (SHAP)

在存在交互作用的情况下，归因的计算值取决于计算中使用的解释性协变量的顺序。SHAP基于对所有（或大量）可能排序的变量属性值进行平均的想法。

图中10 个随机排序（由每个图中的行的顺序表示）的解释变量的 BD 图，用于预测统一样本。这些图显示了各种变量对不同排序的贡献的明显差异。对于变量fare和class可以观察到最显着的差异，贡献根据顺序改变符号。

为了消除变量排序的影响，我们可以计算属性的平均值。下图中显示的十个排序计算的平均值。分别呈现红色和绿色条，负平均值和正平均值紫色箱形图总结了不同排序中每个解释变量的属性分布。从结果来看，最重要的变量是age、class和gender。

bd_lm <- predict_parts(explainer = aps_lm_explainer4,
                 new_observation = titanic_imputed[1,],
                            type = "shap",B=25)
bd_lm

##                                           min           q1       median
## survived: age = 42               -0.054181242 -0.045599456 -0.038365202
## survived: class = 3rd            -0.094495696 -0.063498867 -0.053847757
## survived: embarked = Southampton -0.012328953 -0.009799728 -0.007578614
## survived: fare = 7.11            -0.069949252 -0.028183280 -0.009963752
## survived: gender = male          -0.139885818 -0.127352968 -0.121020843
## survived: parch = 0              -0.015545990 -0.009640689 -0.005661531
## survived: sibsp = 0              -0.008888083 -0.001809017  0.007355686
## survived: survived = 0            0.000000000  0.000000000  0.000000000
##                                          mean           q3          max
## survived: age = 42               -0.038905374 -0.033879021 -0.024763933
## survived: class = 3rd            -0.050668709 -0.047488899  0.018955143
## survived: embarked = Southampton -0.006297454 -0.001517898  0.002366108
## survived: fare = 7.11            -0.020580444 -0.001175351  0.007544178
## survived: gender = male          -0.112069343 -0.091645446 -0.059080199
## survived: parch = 0              -0.005513584 -0.001940643  0.011640236
## survived: sibsp = 0               0.004528790  0.008734708  0.016188491
## survived: survived = 0            0.000000000  0.000000000  0.000000000

plot(bd_lm)

How dose a variable affect the prediction

Ceteris-paribus Profiles

这种方法的基本思想是：根据变量值变化引起的模型预测变化来评估所选解释变量的影响。该方法通过假设所有其他变量的值不变来检查解释变量的影响。主要目标是了解变量值的变化如何影响模型的预测。

这个方法用于对单个探索性变量的值发生变化，模型的预测将如何变化！

本质上，CP Profiles 显示因变量（响应）的条件期望对特定解释变量值的依赖性

该方法根据变量值变化引起的模型预测变化来评估所选解释变量的影响。该方法基于其他条件不变的原则。如果我们可以通过分别调查解释变量的影响并一次更改一个解释变量来探索模型，那么似乎更容易理解黑盒模型的工作原理。

cp_titanic_rf <- predict_profile(explainer = (aps_lm_explainer4), 
                           new_observation = titanic_imputed[1,])
cp_titanic_rf

## Top profiles    : 
##       gender        age class    embarked fare sibsp parch survived _yhat_
## 1     female 42.0000000   3rd Southampton 7.11     0     0        0  0.402
## 1.1     male 42.0000000   3rd Southampton 7.11     0     0        0  0.006
## 11      male  0.1666667   3rd Southampton 7.11     0     0        0  0.482
## 1.110   male  0.9050000   3rd Southampton 7.11     0     0        0  0.596
## 1.2     male  1.6433333   3rd Southampton 7.11     0     0        0  0.576
## 1.3     male  2.3816667   3rd Southampton 7.11     0     0        0  0.578
##       _vname_ _ids_  _label_
## 1      gender     1 survived
## 1.1    gender     1 survived
## 11        age     1 survived
## 1.110     age     1 survived
## 1.2       age     1 survived
## 1.3       age     1 survived
## 
## 
## Top observations:
##   gender age class    embarked fare sibsp parch survived _yhat_  _label_ _ids_
## 1   male  42   3rd Southampton 7.11     0     0        0  0.006 survived     1

plot(cp_titanic_rf,variables = "embarked")

## 'variable_type' changed to 'categorical' due to lack of numerical variables.
## 'variable_type' changed to 'categorical' due to lack of numerical variables.

这种方法提供了一种统一、易于交流且可扩展的模型探索方法。它们的图形表示易于理解和解释。

Ceteris-paribus Oscillations

这算是一种计算变量重要性的方法！

对其他条件不变 (CP) 剖面的目视检查是有见地的。但是，如果模型具有大量解释变量，我们最终可能会得到大量可能不堪重负的图。在这种情况下，选择最有趣或最重要变量很有用。

值得注意的是，解释变量对特定实例预测的影响越大，对应的 CP 剖面的波动就越大。对于对模型预测影响很小或没有影响的变量，曲线将是平坦的或几乎不会改变。换句话说，CP 配置文件的值应该接近模型对特定实例的预测值。

bd_lm <- predict_parts(explainer = aps_lm_explainer4,
                 new_observation = titanic_imputed[1,],
                            type = "oscillations_uni")
bd_lm

##    _vname_ _ids_ oscillations
## 5     fare     1   0.38138614
## 1   gender     1   0.19800000
## 3    class     1   0.12485714
## 2      age     1   0.11922772
## 7    parch     1   0.05550495
## 6    sibsp     1   0.05075248
## 4 embarked     1   0.00350000
## 8 survived     1   0.00000000

plot(bd_lm)

基于 ceteris-paribus 振荡的可变重要性度量。

Does the models fit well around the prediction

本地诊断图

可能会发生这样的情况，尽管模型的预测性能总体上令人满意，但模型对某些观察结果的预测却非常糟糕。在这种情况下，人们常说“模型没有很好地覆盖输入空间的某些区域”。

例如，适用于某家医院“典型”患者数据的模型可能不适用于另一家医院可能具有不同特征的患者。或者，为评估春季假期消费贷款的风险而开发的模型可能不适用于圣诞节假期礼物的秋季贷款。

出于这个原因，在做出重要决策的情况下，有必要检查模型在与感兴趣实例相似的观察中的局部行为。

种解决此问题的本地诊断技术。第一个是局部保真图（local-fidelity plots），用于评估模型围绕感兴趣的观察的局部预测性能。第二个是局部稳定性图（local-stability plots），用于评估围绕感兴趣的观察的预测的（局部）稳定性。

局部保真图（local-fidelity plots）

基本思想：假设，对于感兴趣的观察，我们已经从训练数据中识别出一组具有相似特征的观察。我们将这些类似的观察结果称为“邻居”。局部保真图背后的基本思想是将邻居的残差分布（即因变量的观测值和预测值之间的差异；如果分布图呈现合理的分布，则说明模型的整体性能合理。

rm(list = ls())
library(randomForest)
titanic_glm_model <- randomForest(survived ~ .,data = titanic_imputed)

## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?

explainer_glm <- DALEX::explain(titanic_glm_model,
                         data = titanic_imputed,
                         y = titanic_imputed$survived)

## Preparation of a new explainer is initiated
##   -> model label       :  randomForest  (  default  )
##   -> data              :  2207  rows  8  cols 
##   -> target variable   :  2207  values 
##   -> predict function  :  yhat.randomForest  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package randomForest , ver. 4.7.1.1 , task regression (  default  ) 
##   -> predicted values  :  numerical, min =  0.01158472 , mean =  0.3215986 , max =  0.9911881  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -0.788393 , mean =  0.0005581408 , max =  0.9083918  
##   A new explainer has been created!

johny_d <- titanic_imputed[24,]

id_johny <- predict_diagnostics(explainer_glm, johny_d, variables = NULL)

## Warning in ks.test.default(residuals_other, residuals_sel): p-value will be
## approximate in the presence of ties

id_johny

## 
##  Asymptotic two-sample Kolmogorov-Smirnov test
## 
## data:  residuals_other and residuals_sel
## D = 0.27986, p-value = 0.0009476
## alternative hypothesis: two-sided

plot(id_johny)

局部稳定性图（local-stability plots）

局部稳定性图背后的想法是检查解释变量的微小变化（以邻居集内的变化为代表）是否对预测产生了很大影响。如果轮廓几乎平行并且彼此非常接近，则表明模型的预测围绕感兴趣的实例是稳定的。

当然，不同解释变量的 CP曲线可能非常不同，所以一个自然的问题是：我们应该检查哪些变量？显而易见的选择是根据可变重要性度量。

函数predict_diagnostics()也可用于构建局部稳定性图。为实现这一目标，我们必须选择我们要为其创建绘图的解释变量。

id_johny <- predict_diagnostics(explainer_glm, johny_d, variables = "class")
id_johny

## Top profiles    : 
##      gender age            class    embarked   fare sibsp parch survived
## 24   female   2              1st Southampton 151.16     1     2        0
## 24.1 female   2              2nd Southampton 151.16     1     2        0
## 24.2 female   2              3rd Southampton 151.16     1     2        0
## 24.3 female   2        deck crew Southampton 151.16     1     2        0
## 24.4 female   2 engineering crew Southampton 151.16     1     2        0
## 24.5 female   2 restaurant staff Southampton 151.16     1     2        0
##         _yhat_ _vname_ _ids_      _label_
## 24   0.4648416   class    24 randomForest
## 24.1 0.6369031   class    24 randomForest
## 24.2 0.2900404   class    24 randomForest
## 24.3 0.3468279   class    24 randomForest
## 24.4 0.2879101   class    24 randomForest
## 24.5 0.3404406   class    24 randomForest
## 
## 
## Top observations:
##    gender age class    embarked   fare sibsp parch survived    _yhat_
## 24 female   2   1st Southampton 151.16     1     2        0 0.4648416
##         _label_ _ids_
## 24 randomForest     1

plot(id_johny)

## 'variable_type' changed to 'categorical' due to lack of numerical variables.
## 'variable_type' changed to 'categorical' due to lack of numerical variables.

How good is the model

一些评估（预测）模型的整体性能有用的度量。

模型指标总结

(eva_rf <- DALEX::model_performance(explainer_glm))

## Measures for:  regression
## mse        : 0.1090639 
## rmse       : 0.3302482 
## r2         : 0.5005587 
## mad        : 0.1725322
## 
## Residuals:
##          0%         10%         20%         30%         40%         50% 
## -0.78839299 -0.25390932 -0.20797935 -0.16576801 -0.14602219 -0.10837962 
##         60%         70%         80%         90%        100% 
## -0.07594227  0.02264149  0.20604309  0.55445529  0.90839181

该DALEX::model_performance()函数的应用返回一个类“model_performance”的对象，其中包括几个模型性能度量的估计值，以及一个包含因变量的观察值和预测值及其差异（即残差）的数据框。可以通过将通用plot()函数应用于对象来构建 ROC 曲线或提升图。

p1 <- plot(eva_rf, geom = "histogram") 
p2 <- plot(eva_rf, geom = "prc") 
p3 <- plot(eva_rf, geom = "boxplot") 
 



library("patchwork")
p1 + p2 + p3

Which variable are important to the model

变量重要性

计算变量的重要性可用于多种目的。

1. 模型简化：不影响模型预测的变量可能会被排除在模型之外。
2. 模型探索：比较变量在不同模型中的重要性可能有助于发现变量之间的相互关系。此外，变量在其重要性函数中的排序有助于决定我们应该以何种顺序进行进一步的模型探索。
3. 基于领域知识的模型验证：识别最重要的变量可能有助于评估基于领域知识的模型的有效性。
4. 知识生成：识别最重要的变量可能会导致发现特定机制中涉及的新因素。 评估变量重要性的方法通常可以分为两组：模型特定的和模型不可知的。

对于线性模型和许多其他类型的模型，有一些方法可以利用模型结构的特定元素来评估解释变量的重要性。这些是特定于模型的方法。例如，对于线性模型，可以使用归一化回归系数的值或其对应的 p 值作为变量重要性度量。对于基于树的集合，这样的度量可以基于特定树中特定变量的使用。这方面的一个很好的例子是基于袋外数据的随机森林模型的可变重要性度量（Leo Breiman 2001 a），但也有其他方法，如XgboostExplainer包中实现的方法（Foster 2017）用于梯度提升和randomForestExplainer （Paluszynska 和 Biecek 2017）用于随机森林。

DALEXR 包中基于排列的可变重要性度量的实现。关键功能是model_parts()允许计算度量。为了计算的目的，可以在包括loss_sum_of_squares()、loss_root_mean_square()、loss_accuracy()、loss_cross_entropy()和在内的几个损失函数中进行选择loss_one_minus_auc()。

为了计算基于排列的可变重要性度量，我们应用该model_parts()函数。唯一需要的参数是explainer。其他（可选）参数是：

1. loss_function，要使用的损失函数（默认情况下，它是loss_root_mean_square函数）。
2. type，变量重要性度量的形式，raw，difference，ratio
3. variables，一个字符向量，提供解释变量的名称，要为其计算变量重要性度量。
4. variable_groups，解释变量名称的特征向量列表。对于每个向量，为向量中提供名称的变量的联合效应计算单个变量重要性度量。
5. B，用于计算（平均）可变重要性度量的排列数，B = 10默认使用。要获得基于单排列的度量，请使用B = 1.
6. N，为了计算可变重要性度量而从解释器对象中可用的数据中抽样的观察次数；默认情况下，N = 1000使用；如果N = NULL，则使用整个数据集。

set.seed(1980)
x <- model_parts(explainer = explainer_glm, 
        loss_function = loss_root_mean_square,
                    B = 1)

plot(x)

# set.seed(1980)
# model_parts(explainer = explainer_rf, 
#         loss_function = loss_root_mean_square,
#                     B = 1,
#             variables = colnames(explainer_rf$data))

# set.seed(1980)
# (vip.50 <- model_parts(explainer = explainer_rf, 
#                    loss_function = loss_root_mean_square,
#                                B = 50,
#                             type = "difference"))

# plot(vip_rf, vip_svm, vip_lm) 同时绘制三种模型结果

How dose a variable affect the average prediction

Partial-dependence Profiles

部分依赖 (PD) 图，有时也称为 PD 图。是由Friedman ( 2000 )在梯度提升机 (GBM) 的背景下引入的。

基本思想是展示模型预测的预期值如何作为选定解释变量的函数。PD是通过数据集中所有实例（观察）的 CP的平均值来估计的。

rm(list = ls())
library(randomForest)
titanic_glm_model <- randomForest(survived ~ .,data = titanic_imputed)

## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?

explainer_glm <- DALEX::explain(titanic_glm_model,
                         data = titanic_imputed,
                         y = titanic_imputed$survived)

## Preparation of a new explainer is initiated
##   -> model label       :  randomForest  (  default  )
##   -> data              :  2207  rows  8  cols 
##   -> target variable   :  2207  values 
##   -> predict function  :  yhat.randomForest  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package randomForest , ver. 4.7.1.1 , task regression (  default  ) 
##   -> predicted values  :  numerical, min =  0.01431641 , mean =  0.3217734 , max =  0.9894182  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -0.7811962 , mean =  0.0003833328 , max =  0.908653  
##   A new explainer has been created!

DALEX包中计算 PD 配置文件的函数是model_profile(). 唯一需要的参数是explainer。其他有用参数包括：

1. variables, 一个字符向量，提供要计算轮廓的解释变量的名称；默认情况下，variables = NULL，在这种情况下，将对模型中包含的所有数值变量执行计算。
2. N, 用于计算 PD 配置文件的（随机采样）观测值的数量（N = 100默认情况下）；N = NULL暗示使用解释器对象中包含的整个数据集。
3. type，PD 配置文件的类型，具有值”partial”（默认值）“conditional”、 和”accumulated”。
4. variable_type, 一个字符串，指示是否应仅对”numerical”（连续）解释变量（默认）或仅对”categorical”变量执行计算。
5. groups，将用于对配置文件进行分组的解释变量的名称，groups = NULL默认情况下为（在这种情况下不应用配置文件分组）。
6. k,hclust(), 在函数的帮助下要创建的集群的数量，k = NULL默认情况下使用并且意味着没有集群。

pdp_rf <- model_profile(explainer = explainer_glm, variables = "age")
pdp_rf

## Top profiles    : 
##   _vname_      _label_        _x_    _yhat_ _ids_
## 1     age randomForest  0.1666667 0.5501549     0
## 2     age randomForest  2.0000000 0.5759993     0
## 3     age randomForest  4.0000000 0.5939621     0
## 4     age randomForest  7.0000000 0.5573200     0
## 5     age randomForest  9.0000000 0.5287923     0
## 6     age randomForest 13.0000000 0.4777920     0

使用plot函数，获得了 PD 曲线图。

library("ggplot2")

## 
## Attaching package: 'ggplot2'

## The following object is masked from 'package:randomForest':
## 
##     margin

plot(pdp_rf) +  ggtitle("Partial-dependence profile for age") 

PD 曲线可以绘制在 CP 曲线之上。

plot(pdp_rf, geom = "profiles") + 
    ggtitle("Ceteris-paribus and partial-dependence profiles for age") 

聚类pd

要计算聚类 PD ，我们必须对 CP 进行聚类。为了这个目标，我们使用函数的k参数

pdp_rf_clust <- model_profile(explainer = explainer_glm, 
                              variables = "age", k = 3)

pdp_rf_clust

## Top profiles    : 
##   _vname_        _label_        _x_ _cluster_    _yhat_ _ids_
## 1     age randomForest_1  0.1666667         1 0.5913137     0
## 2     age randomForest_1  2.0000000         1 0.6337502     0
## 3     age randomForest_1  4.0000000         1 0.6499235     0
## 4     age randomForest_1  7.0000000         1 0.6120171     0
## 5     age randomForest_1  9.0000000         1 0.5733029     0
## 6     age randomForest_1 13.0000000         1 0.4988271     0

plot(pdp_rf_clust)

通过使用geom = “profiles”函数中的参数，可以在 CP 配置文件的顶部绘制集群 PD 配置文件plot()。

plot(pdp_rf_clust, geom = "profiles") + 
    ggtitle("Clustered partial-dependence profiles for age") 

分组 pd

设置groups 参数

pdp_rf_gender <- model_profile(explainer = explainer_glm, 
                               variables = "age", groups = "gender")

pdp_rf_gender

## Top profiles    : 
##   _vname_             _label_        _x_ _groups_    _yhat_ _ids_
## 1     age randomForest_female  0.1666667   female 0.7308659     0
## 2     age randomForest_female  2.0000000   female 0.7154967     0
## 3     age randomForest_female  4.0000000   female 0.7418063     0
## 4     age randomForest_female  7.0000000   female 0.7306238     0
## 5     age randomForest_female  9.0000000   female 0.7405401     0
## 6     age randomForest_female 13.0000000   female 0.7608780     0

plot(pdp_rf_gender)

plot(pdp_rf_gender, geom = "profiles") + 
    ggtitle("Partial-dependence profiles for age, grouped by gender") 

# 如果要比较两个模型的pd，直接调用plot函数
# plot(pdp_rf, pdp_lmr) 

Local-dependence and Accumulated-local Profiles LD 和LA

PD曲线是CP曲线的平均值。

LD

rm(list = ls())
library(randomForest)
titanic_glm_model <- randomForest(survived ~ .,data = titanic_imputed)

## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?

explainer_glm <- DALEX::explain(titanic_glm_model,
                         data = titanic_imputed,
                         y = titanic_imputed$survived)

## Preparation of a new explainer is initiated
##   -> model label       :  randomForest  (  default  )
##   -> data              :  2207  rows  8  cols 
##   -> target variable   :  2207  values 
##   -> predict function  :  yhat.randomForest  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package randomForest , ver. 4.7.1.1 , task regression (  default  ) 
##   -> predicted values  :  numerical, min =  0.01086255 , mean =  0.3228053 , max =  0.9911796  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -0.7861884 , mean =  -0.0006484816 , max =  0.9012987  
##   A new explainer has been created!

ld_rf <- model_profile(explainer = explainer_glm,
                      type       = "conditional",
                      variables  = c("age", "sibsp"))

ld_rf

## Top profiles    : 
##                               _vname_      _label_        _x_    _yhat_ _ids_
## age.randomForest.0.1666666667     age randomForest  0.1666667 0.5541816     0
## age.randomForest.2                age randomForest  2.0000000 0.5905590     0
## age.randomForest.4                age randomForest  4.0000000 0.5991920     0
## age.randomForest.7                age randomForest  7.0000000 0.5639429     0
## age.randomForest.9                age randomForest  9.0000000 0.5472621     0
## age.randomForest.13               age randomForest 13.0000000 0.4651306     0

plot(ld_rf,geom = "profiles")

LA

ld_rf <- model_profile(explainer = explainer_glm,
                      type       = "accumulated",
                      variables  = c("age", "sibsp"))

ld_rf

## Top profiles    : 
##   _vname_      _label_        _x_    _yhat_ _ids_
## 1     age randomForest  0.1666667 0.5238453     0
## 2     age randomForest  2.0000000 0.5586784     0
## 3     age randomForest  4.0000000 0.5675521     0
## 4     age randomForest  7.0000000 0.5249199     0
## 5     age randomForest  9.0000000 0.5077423     0
## 6     age randomForest 13.0000000 0.4309084     0

plot(ld_rf, geom = "profiles")

# plot(pd_rf, ld_rf, al_rf)  同时显示所有图

Does the model fit well in general

残差诊断

dm <- model_diagnostics(explainer_glm)

默认情况下，将plot()函数应用于model_diagnostics-class 对象会根据因变量的预测值（在水平轴上）生成残差散点图（在垂直轴上）。通过使用参数variable和yvariable，可以用其他变量分别指定用于水平轴和垂直轴的图。除了解释变量的名称外，这两个参数接受以下值：

1. “y”对于因变量，
2. “y_hat”对于因变量的预测值，
3. “obs”对于观察的标识符，
4. “residuals”对于残差，
5. “abs_residuals”对于残差的绝对值。

plot(dm, variable = "y_hat", yvariable = "y") 

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'