特征依賴分析（Feature Dependency Analysis）對于模型分析而言至關(guān)重要，對于注重可解釋性的工業(yè)界應(yīng)用更是如此。對于一個模型，一般會有多個輸入，這些輸入對于整個模型輸出結(jié)果的影響程度是各不相同的，比如對于銀行客戶違約預(yù)測模型而言，可能該用戶的歷史違約行為輸入對于模型輸出違約的影響最大。對于某個模型而言，我們希望可以分析出哪些特征對其輸出的影響最大，而這種分析是一種全局的模型無關(guān)的分析（Global model-agnostic），也就是說

全局（Global）：意味著這個分析是以整個模型為單位進(jìn)行的，給定一定數(shù)量的目標(biāo)樣本集，我們從樣本集的宏觀角度，去評估某維特征對于模型的影響。如果期望給定一個樣本，然后評估這個樣本的某維特征的影響，這是全局分析做不到的。模型無關(guān)（Model-Agnostic）：意味著這種分析是適用于所有模型的，并不是某個模型（比如樹模型）所獨有的。以下介紹幾種常用的方法，可用于分析模型的特征依賴。

部分依賴曲線（Partial Dependency Plot）部分依賴曲線（Partial Dependency Plot, PDP）就是這樣一種用以評估某維特征對于模型輸出影響的可視化分析方法，而且這種方法也是全局的模型無關(guān)的。總體而言，這個方法會對模型輸入的特征的所有可能取值，都給出其平均的模型輸出打分，然后繪制成一個曲線。這種方法有一個基本假設(shè)

Fundamental Assumption： 模型的輸入都是獨立無關(guān)的。

在滿足這個假設(shè)的前提下，假如我們已經(jīng)訓(xùn)練好了某個模型f，如果我們需要對輸入 x 的第 S 維特征進(jìn)行依賴分析，將其記為 ，剩余的其他維度特征記為? 。對于n nn個給定的樣本，其中的 d 為特征維度，由于和X獨立無關(guān)，因此維持樣本原本的特征不變，將設(shè)置為待分析的目標(biāo)值v vv，此時模型可以得到打分（注意到我們有 n 個樣本，因此有n nn個模型輸出打分），求平均值，其中的表示第 i 個樣本的剩余特征值。顯然，這個平均值就是當(dāng)特征取值為v vv的時候，該模型的平均輸出值大小。整個過程同樣可以用公式(1-1)進(jìn)行表示。

注意到這個方法描述的是模型在整個數(shù)據(jù)集的維度，當(dāng)其他維度特征不變的情況下，待分析特征的平均輸入輸出響應(yīng)曲線。筆者準(zhǔn)備了一些例子，具體見[5]的partial_dependency_plot_demo.ipynb，這里簡單進(jìn)行說明，在本例子中采用的是sklearn.dataset中的boston房價預(yù)測數(shù)據(jù)集，boston['data']大小為506 × 13 ，表示506個樣本，13維特征，而boston['target']是這506個樣本的真實房價。一共有13維特征可供部分依賴曲線分析，我們優(yōu)先選擇哪個或者哪些進(jìn)行呢？我們可以通過特征權(quán)重（也稱之為特征重要性）進(jìn)行篩選，特征權(quán)重可以通過樹模型的分裂增益進(jìn)行統(tǒng)計，也可以通過改變某個特征的排序，計算指標(biāo)的diff大小進(jìn)行統(tǒng)計，具體的方法我們以后博文再聊。如Fig 1.1所示，其中的Feature Importance就是通過分裂增益統(tǒng)計的特征權(quán)重排序，而Permutation Importance就是亂序某維特征排序統(tǒng)計出來的特征權(quán)重排序。從中我們可以發(fā)現(xiàn)，第12和第5維特征是最為重要的特征，我們嘗試對這兩維特征的PDP進(jìn)行繪制。

Fig 1.1 兩種計算特征權(quán)重的方式，（1） 通過樹模型的節(jié)點分裂增益進(jìn)行統(tǒng)計，（2）通過亂序某維特征，統(tǒng)計亂序特征帶來的指標(biāo)變化從而統(tǒng)計重要性。前者是限定在樹模型中使用，后者可以是模型無關(guān)的。

如Fig 1.2所示，筆者在（a）和（b）分別繪制了第12維特征和第5維特征的PDP曲線，同時繪制了聯(lián)合第12維和第5維的2維PDP曲線，如（c）所示。2維PDP曲線將式子(1-1)中的 從一個變量變?yōu)閮蓚€變量，可以描述兩個變量對輸出的共同影響。從圖中我們能發(fā)現(xiàn)一些有意思的點：

權(quán)重高的特征，其輸入對于輸出的影響都是較為遞增/遞減的，整體趨勢比較統(tǒng)一。權(quán)重高的特征，的輸出響應(yīng)范圍大，比如（a）的能從-6到8，峰值差達(dá)到了14，（b）的范圍從-2到8，峰值差達(dá)到了10。權(quán)重高低與否，與PDP曲線是遞增趨勢或者遞減趨勢無關(guān)。部分依賴曲線的趨勢只能代表該特征從平均的角度，給模型的輸出帶來的影響是正向的亦或是負(fù)向的，并不能決定這個權(quán)重的重要與否。權(quán)重的重要性高低從PDP曲線中，更應(yīng)該從打分峰值差距中判斷，能打出更大峰值差距的特征會傾向于更重要。如Fig 1.2 （a）所示，第12維特征是遞減趨勢的特征，但是其特征權(quán)重確是最高的。我們理解這維特征的時候，應(yīng)該理解為該特征在對整個模型的輸出上有著顯著的影響，這體現(xiàn)在對模型的輸出能產(chǎn)生更大范圍的影響，但是這個影響從預(yù)測打分的角度而言是負(fù)向的。只要我們將第12維特征取個相反數(shù)，如Fig 1.2 （d）所示，我們就能發(fā)現(xiàn)第12維特征的趨勢會變成遞增，但是此時特征權(quán)重沒有任何變化。同理的，第5維特征也是如此，如Fig 1.2（b）所示，只是該特征對于模型輸出而言是正向的，意味著特征值越大，模型輸出傾向于越大。

Fig 1.2 第12和第5維特征的部分依賴曲線圖，（a）為第12維特征的PDP曲線，（b）為第5維特征的PDP曲線，（c）為聯(lián)合第12維和第5維的2維PDP曲線。

為了繼續(xù)驗證這個觀點，我們對其他維度（0，2，6，9，10，11維）的特征進(jìn)行了PDP曲線繪制，如Fig 1.3 （a）-（f）所示，我們發(fā)現(xiàn)這些權(quán)重不高的特征，其打分峰值差約為1.4，0.8，2，2，2.5，1.4，的確比第12和第5維的特征打分峰值差低很多。同時，我們也能發(fā)現(xiàn)，在以下的PDP曲線中存在很多特征的模型輸出趨勢是較為動蕩的，如（a）、（b）、（d）。這類型的特征對于模型輸出來說，通常具備較高的非線性，由于PDP曲線表示的平均值的關(guān)系，如果某個區(qū)段中平均值之間有著距離的波動，那么很可能其樣本中，該特征在這個區(qū)段的影響方差也很大，這樣在線上很容易導(dǎo)致一些bad case，需要著重分析和考慮。這類型的區(qū)段見Fig 1.3中的紅線虛框。

Fig 1.3 （a）-（f）列舉了其他6維特征的PDP曲線，其中的紅線虛框表示PDP變化劇烈的區(qū)段。PDP曲線雖然能描述某維模型特征，對于模型輸出的全局、模型無關(guān)的依賴分析，方法直觀而且簡單，但是具有以下缺點。

1、PDP曲線最多只能分析到2維的聯(lián)合特征分析。如Fig 1.2（c）所示，PDP通常只能對一維特征和兩維特征進(jìn)行特征依賴分析（因為高于2維將無法進(jìn)行可視化）；2、PDP曲線計算代價高。對于某維特征而言，其中的某個值v vv的對應(yīng)1需要進(jìn)行 n 次模型打分計算，對于一個具有m mm中可能取值的特征而言，其模型計算復(fù)雜度就是，對于 k 個特征而言就是，其復(fù)雜度還是很高的。3、PDP曲線依賴于特征間獨立無關(guān)的假設(shè)。通常數(shù)據(jù)的特征之間都不會是獨立無關(guān)的，這意味著我們在遍歷某特征的可能取值的時候，可能會出現(xiàn)“不合理”的情況。舉個例子，假設(shè)某個特征 x 和某個特征 y 的關(guān)系為 y = x ，顯然這兩個特征是相關(guān)的。如果此時對特征y yy進(jìn)行PDP繪制，那么會遍歷 y  的所有可能取值范圍，比如是range(0,100,0.5)，y 作為被看成獨立無關(guān)的特征是需要被固定的，此時會出現(xiàn)個“詭異”的情況，當(dāng) x 遍歷到20的時候， y  可能為10（因為固定y yy特征），這違背了y = x 。顯然，對于非獨立無關(guān)的數(shù)據(jù)而言，PDP的繪制策略將會遍歷到一些不可能存在的數(shù)據(jù)點，這會導(dǎo)致繪制出來的點不置信。4、PDP曲線可能隱藏了數(shù)據(jù)里面的隱藏模式。由于有PDP曲線只是對平均值進(jìn)行計算，這種全局方法可能導(dǎo)致一些數(shù)據(jù)中模式被隱藏。舉個例子，假如對于某個特征而言，數(shù)據(jù)集中有一半數(shù)據(jù)是模型輸出和輸入特征正相關(guān)，一半數(shù)據(jù)是模型輸出和輸入特征負(fù)相關(guān)，這平均來看，最終PDP的結(jié)果可能是繪制出一個直線！顯然這不是我們期望的。解決方案是采用獨立條件期望曲線（Individual Conditional Expectation curve， ICE）[6]進(jìn)行分析，這個屬于局部模型無關(guān)模型，我們以后再聊。以上的第3點可以用Fig 1.4舉例，假設(shè)total_bill為隨機(jī)變量 X , tip為隨機(jī)變量 Y ，那么Fig 1.4 （a）繪制了邊緣分布P ( X ) 和P ( Y ) ，可知，而Fig 1.4（b）表示的是當(dāng) X 取值在( 25 , 30 ) 之間時候的條件分布。可以發(fā)現(xiàn)，邊緣分布和某個區(qū)間內(nèi)的條件分布有很大的差別，如果采用部分依賴曲線，則會假設(shè)每一個區(qū)間內(nèi)的條件分布和全局的邊緣分布一致，這在兩個變量非獨立無關(guān)的情況下顯然是不合理的！而當(dāng)兩變量是獨立無關(guān)的時候，如Fig 1.4（c）和（d）所示，此時邊緣分布和條件分布的差別不大，這才是符合部分依賴曲線分析的假設(shè)的場景~ 此處相關(guān)的繪圖代碼見Appendix Code A.和Code B.。

Fig 1.4 當(dāng)特征之間（如tip和total_bill之間）不獨立無關(guān)的時候，（a）total_bill和tip各自的邊緣分布，（b）當(dāng)total_bill值在(25,30)之間時候，tip的條件分布。對比（a）和（b）可以發(fā)現(xiàn)tip的邊緣分布和條件分布之間有較大差別。 當(dāng)特征之間獨立無關(guān)的時候，如（c）和（d）所示，邊緣分布和條件分布差別不大。

條件依賴曲線（Marginal Plot）

部分依賴曲線強(qiáng)烈依賴于特征之間獨立無關(guān)的假設(shè)，然而現(xiàn)實生活中的模型輸入特征間無法保證都是獨立無關(guān)的，怎么解決這個問題呢？我們可以對條件分布進(jìn)行求平均，而不是對邊緣分布進(jìn)行求平均，這個方法稱之為條件依賴曲線（Marginal-Plot，M-Plot），通過這種方法我們可以避免對現(xiàn)實中不可能存在的數(shù)據(jù)點進(jìn)行求平均。當(dāng)然，在實際中我們通常是對某個區(qū)間內(nèi)的所有數(shù)據(jù)點進(jìn)行求平均，用公式表示如式子(2-1)所示。

然而，采用M-Plot雖然可以避免對不可能存在的數(shù)據(jù)點進(jìn)行求平均，但是仍然會受到特征之間非獨立無關(guān)的影響。如Fig 2.1所示，假設(shè)特征X對于模型輸出來說，在真實情況下其實并沒有太大影響，但是和存在強(qiáng)相關(guān)，而對于模型輸出來說影響重要。我們對處進(jìn)行M-Plot計算，由于 和之間存在強(qiáng)相關(guān)關(guān)系，即便只考慮對條件分布進(jìn)行求平均，仍然很可能會導(dǎo)致這種情況，而這個分析是不準(zhǔn)確的，因為這個其實并不是對于模型的影響，而是由和的相關(guān)關(guān)系而被傳遞的影響。M-Plot會混合一些相關(guān)特征的影響，此時我們就被誤導(dǎo)了。可以舉個栗子，此處的是天空中的烏云數(shù)量，是天空中降雨量，而我們預(yù)估的Y是街道上行人打傘的數(shù)量，可知降雨量和天空中烏云數(shù)量強(qiáng)相關(guān)，但是從我們直觀理解上，天空中的烏云數(shù)量不應(yīng)該會影響行人是否打傘，而降雨量才會。如果我們從條件依賴曲線中得到“天空中烏云越多，行人打傘的可能越大”的結(jié)論，其實是不符合邏輯的。

Fig 2.1 條件依賴曲線，對條件分布進(jìn)行求平均，從而解決PDP中會取得現(xiàn)實中不存在的樣本點的問題。然而M-Plot仍然會受到特征間非獨立無關(guān)的影響，導(dǎo)致不能單獨反應(yīng)某個特征對輸出的真實影響。根據(jù)本章對Marginal-Plot的介紹，我們嘗試對上一章的Partial Dependency Plot的圖示進(jìn)行對應(yīng)的Marginal-Plot的繪制和對比，采用的代碼見Appendix Code C.，繪制結(jié)果如Fig 2.2所示。我們可以發(fā)現(xiàn)幾點區(qū)別：

對于第12維和5維特征而言，采用Mplot和PDP的方法繪制出來的曲線趨勢一致，但是對于第2和9維特征而言則有比較大的區(qū)別。Mplot和PDP的打分范圍差別很大。這兩點很好理解，由于PDP會網(wǎng)格遍歷某維特征取值范圍內(nèi)的所有取值，那么就可能會取得一些真實數(shù)據(jù)中不可能存在的數(shù)據(jù)點，導(dǎo)致取平均后打分有偏差。這個會導(dǎo)致在PDP曲線中出現(xiàn)一些突變，或者打分范圍和真實的打分范圍有所偏差。Mplot對比PDP有很多優(yōu)勢，比如其計算速度更快（不需要進(jìn)行網(wǎng)格遍歷），避免對不可能存在的數(shù)據(jù)點進(jìn)行求平均導(dǎo)致繪圖不置信等等，但是其還是會受到特征之間存在依賴關(guān)系導(dǎo)致的“依賴傳遞”的影響。

Fig 2.2 Marginal-Plot 和 對應(yīng)的 Partial Dependency Plot的對比，(a)-(d)分別是第12，5，2，9維特征的對比。為了解決這個問題，我們可以繼續(xù)探索Accumulate Local Estimate（ALE）累計局部估計，我們后文見。

Reference

[1]. https://christophm.github.io/interpretable-ml-book/pdp.html, Partial Dependence Plot (PDP)

[2]. https://zhuanlan.zhihu.com/p/428466235，

[3]. https://christophm.github.io/interpretable-ml-book/ale.html， Accumulated Local Effects (ALE) Plot

[4]. https://scikit-learn.org/stable/modules/partial_dependence.html, Partial Dependence and Individual Conditional Expectation plots

[5]. https://github.com/FesianXu/FeatureDependencyAnalysisDemo

[6]. https://christophm.github.io/interpretable-ml-book/ice.html#ice，Individual Conditional Expectation (ICE)

AppendixCode A. 非獨立無關(guān)下，條件分布和邊緣分布的分布差別

import seaborn as sns
sns.set_theme(style="darkgrid")
# marginal distribution
tips = sns.load_dataset("tips")
g = sns.jointplot(x="total_bill", y="tip", data=tips,
                  kind="reg", truncate=False,
                  xlim=(0, 60), ylim=(0, 12),
                  color="m", height=7)
# conditional distribution
tip_range = tips[(tips['total_bill'] < 30) & (tips['total_bill'] > 25)]
g = sns.jointplot(x="total_bill", y="tip", data=tip_range,
                  kind="reg", truncate=False,
                  xlim=(0, 60), ylim=(0, 12),
                  color="m", height=7)

Code B. 獨立無關(guān)下，條件分布和邊緣分布的分布差別

import numpy as np
import seaborn as sns
x = np.random.uniform(0,60,size=1000)
y = np.random.uniform(0,60,size=1000)
inps = {"x":x, "y":y}
inps = pd.DataFrame(inps)
g = sns.jointplot(x="x", y="y", data=inps,
                  kind="reg", truncate=False,
                  xlim=(0, 80), ylim=(0, 80),
                  color="m", height=7)
inps_ = inps[(inps['x'] < 30) & (inps['x'] > 25)]
g = sns.jointplot(x="x", y="y", data=inps_,
                  kind="reg", truncate=False,
                  xlim=(0, 80), ylim=(0, 80),
                  color="m", height=7)

Code C. 繪制Marginal-Plot的代碼

def draw_marginal_plot(data, feat_axis, num_parts):
    max_v = np.max(data[:, feat_axis])
    min_v = np.min(data[:, feat_axis])
    step = (max_v - min_v) / num_parts
    df = pd.DataFrame(data)
    data_record = []
    for each_part in range(num_parts):
        low_bound = min_v + each_part * step
        upper_bound = min_v + (each_part + 1) * step
        data_select = df[(df[feat_axis] > low_bound) & (df[feat_axis] < upper_bound)]
        if len(data_select) > 0:
            data_pred = GBDTreg.predict(data_select)
            data_mean = np.array(data_pred).mean()
            data_std = np.std(np.array(data_pred))
            data_record.append((each_part, data_mean, data_std))
    return data_record, step

feat_id = 12
data_record, step = draw_marginal_plot(X_train, feat_id, 50)
data_record = np.array(data_record)
plt.plot(data_record[:, 0] * step, data_record[:, 1])
plt.xlabel("feat {}".format(feat_id))
plt.ylabel("avg predict")
plt.title("Marginal Plot")