特征工程（上）

特征選擇 (feature_selection)

Filter

1. 移除低方差的特征 (Removing features with low variance)
2. 單變量特征選擇 (Univariate feature selection)

Wrapper

3. 遞歸特征消除 (Recursive Feature Elimination)

Embedded

4. 使用SelectFromModel選擇特征 (Feature selection using SelectFromModel)
5. 將特征選擇過程融入pipeline (Feature selection as part of a pipeline)

當數據預處理完成后，我們需要選擇有意義的特征輸入機器學習的算法和模型進行訓練。

通常來說，從兩個方面考慮來選擇特征：

特征是否發散

如果一個特征不發散，例如方差接近於0，也就是說樣本在這個特征上基本上沒有差異，這個特征對於樣本的區分並沒有什么用。

特征與目標的相關性

這點比較顯見，與目標相關性高的特征，應當優選選擇。除移除低方差法外，本文介紹的其他方法均從相關性考慮。

根據特征選擇的形式又可以將特征選擇方法分為3種：

• Filter：過濾法，按照發散性或者相關性對各個特征進行評分，設定閾值或者待選擇閾值的個數，選擇特征。
• Wrapper：包裝法，根據目標函數（通常是預測效果評分），每次選擇若干特征，或者排除若干特征。
• Embedded：嵌入法，先使用某些機器學習的算法和模型進行訓練，得到各個特征的權值系數，根據系數從大到小選擇特征。類似於Filter方法，但是是通過訓練來確定特征的優劣。

特征選擇主要有兩個目的：

• 減少特征數量、降維，使模型泛化能力更強，減少過擬合；

• 增強對特征和特征值之間的理解。

拿到數據集，一個特征選擇方法，往往很難同時完成這兩個目的。通常情況下，選擇一種自己最熟悉或者最方便的特征選擇方法（往往目的是降維，而忽略了對特征和數據理解的目的）。接下來將結合 Scikit-learn提供的例子 介紹幾種常用的特征選擇方法，它們各自的優缺點和問題。

Filter

1）移除低方差的特征 (Removing features with low variance)

假設某特征的特征值只有0和1，並且在所有輸入樣本中，95%的實例的該特征取值都是1，那就可以認為這個特征作用不大。如果100%都是1，那這個特征就沒意義了。當特征值都是離散型變量的時候這種方法才能用，如果是連續型變量，就需要將連續變量離散化之后才能用。而且實際當中，一般不太會有95%以上都取某個值的特征存在，所以這種方法雖然簡單但是不太好用。可以把它作為特征選擇的預處理，先去掉那些取值變化小的特征，然后再從接下來提到的的特征選擇方法中選擇合適的進行進一步的特征選擇。

from sklearn.feature_selection import VarianceThreshold
X = [[0, 0, 1], [0, 1, 0], [1, 0, 0], [0, 1, 1], [0, 1, 0], [0, 1, 1]]
# 方差低於此閾值的特征將被刪除
sel = VarianceThreshold(threshold=(.8 * (1 - .8)))
sel.fit_transform(X)

array([[0, 1],
       [1, 0],
       [0, 0],
       [1, 1],
       [1, 0],
       [1, 1]])

果然, VarianceThreshold 移除了第一列特征，第一列中特征值為0的概率達到了5/6.

2）單變量特征選擇 (Univariate feature selection)

單變量特征選擇的原理是分別單獨的計算每個變量的某個統計指標，根據該指標來判斷哪些變量重要，剔除那些不重要的變量。

對於分類問題(y離散)，可采用：

• 卡方檢驗
• f_classif
• mutual_info_classif
• 互信息

對於回歸問題(y連續)，可采用：

• 皮爾森相關系數
• f_regression,
• mutual_info_regression
• 最大信息系數

這種方法比較簡單，易於運行，易於理解，通常對於理解數據有較好的效果（但對特征優化、提高泛化能力來說不一定有效）。

• SelectKBest 移除得分前 k 名以外的所有特征(取top k)
• SelectPercentile 移除得分在用戶指定百分比以后的特征(取top k%)
• 對每個特征使用通用的單變量統計檢驗： 假正率(false positive rate) SelectFpr, 偽發現率(false discovery rate) SelectFdr, 或族系誤差率 SelectFwe.
• GenericUnivariateSelect 可以設置不同的策略來進行單變量特征選擇。同時不同的選擇策略也能夠使用超參數尋優，從而讓我們找到最佳的單變量特征選擇策略。

Notice:
　　The methods based on F-test estimate the degree of linear dependency between two random variables. (F檢驗用於評估兩個隨機變量的線性相關性)On the other hand, mutual information methods can capture any kind of statistical dependency, but being nonparametric, they require more samples for accurate estimation.(另一方面，互信息的方法可以捕獲任何類型的統計依賴關系，但是作為一個非參數方法，估計准確需要更多的樣本)

卡方(Chi2)檢驗

經典的卡方檢驗是檢驗定性自變量對定性因變量的相關性。比如，我們可以對樣本進行一次chi2 測試來選擇最佳的兩項特征：

from sklearn.datasets import load_iris
from sklearn.feature_selection import SelectKBest
from sklearn.feature_selection import chi2
iris = load_iris()
X, y = iris.data, iris.target
X.shape

(150, 4)

X[:5, :]

array([[5.1, 3.5, 1.4, 0.2],
       [4.9, 3. , 1.4, 0.2],
       [4.7, 3.2, 1.3, 0.2],
       [4.6, 3.1, 1.5, 0.2],
       [5. , 3.6, 1.4, 0.2]])

X_new = SelectKBest(chi2, k=2).fit_transform(X, y)
X_new.shape

(150, 2)

X_new[:5, :]

array([[1.4, 0.2],
       [1.4, 0.2],
       [1.3, 0.2],
       [1.5, 0.2],
       [1.4, 0.2]])

• f_classif
  ANOVA F-value between label/feature for classification tasks.

• mutual_info_classif
  Mutual information for a discrete target.

• chi2
  Chi-squared stats of non-negative features for classification tasks.

• f_regression
  F-value between label/feature for regression tasks.

• mutual_info_regression
  Mutual information for a continuous target.

• SelectPercentile
  Select features based on percentile of the highest scores.

• SelectFpr
  Select features based on a false positive rate test.

• SelectFdr
  Select features based on an estimated false discovery rate.

• SelectFwe
  Select features based on family-wise error rate.

• GenericUnivariateSelect
  Univariate feature selector with configurable mode.

Pearson相關系數 (Pearson Correlation)

皮爾森相關系數是一種最簡單的，能幫助理解特征和響應變量之間關系的方法，該方法衡量的是變量之間的線性相關性，結果的取值區間為[-1，1]，-1表示完全的負相關，+1表示完全的正相關，0表示沒有線性相關。

import numpy as np
from scipy.stats import pearsonr
np.random.seed(0)
size = 300
x = np.random.normal(0, 1, size)
# pearsonr(x, y)的輸入為特征矩陣和目標向量，能夠同時計算 相關系數 和p-value.
print("Lower noise", pearsonr(x, x + np.random.normal(0, 1, size)))
print("Higher noise", pearsonr(x, x + np.random.normal(0, 10, size)))

Lower noise (0.7182483686213841, 7.32401731299835e-49)
Higher noise (0.057964292079338155, 0.3170099388532475)

這個例子中，我們比較了變量在加入噪音之前和之后的差異。當噪音比較小的時候，相關性很強，p-value很低。

我們使用Pearson相關系數主要是為了看特征之間的相關性，而不是和因變量之間的。

Wrapper

遞歸特征消除 (Recursive Feature Elimination)

遞歸消除特征法使用一個基模型來進行多輪訓練，每輪訓練后，移除若干權值系數的特征，再基於新的特征集進行下一輪訓練。

對特征含有權重的預測模型(例如，線性模型對應參數coefficients)，RFE通過遞歸減少考察的特征集規模來選擇特征。首先，預測模型在原始特征上訓練，每個特征指定一個權重。之后，那些擁有最小絕對值權重的特征被踢出特征集。如此往復遞歸，直至剩余的特征數量達到所需的特征數量。

RFECV 通過交叉驗證的方式執行RFE，以此來選擇最佳數量的特征：對於一個數量為d的feature的集合，他的所有的子集的個數是2的d次方減1(包含空集)。指定一個外部的學習算法，比如SVM之類的。通過該算法計算所有子集的validation error。選擇error最小的那個子集作為所挑選的特征。

from sklearn.feature_selection import RFE
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_iris

rf = RandomForestClassifier()
iris=load_iris()
X,y=iris.data,iris.target
print(X.shape)
print(X[:5, :])

(150, 4)
[[5.1 3.5 1.4 0.2]
 [4.9 3.  1.4 0.2]
 [4.7 3.2 1.3 0.2]
 [4.6 3.1 1.5 0.2]
 [5.  3.6 1.4 0.2]]

# 使用遞歸特征消除進行特征排序
rfe = RFE(estimator=rf, n_features_to_select=3)
X_rfe = rfe.fit_transform(X,y)
print(X_rfe.shape)
print(X_rfe[:5, :])

(150, 3)
[[5.1 1.4 0.2]
 [4.9 1.4 0.2]
 [4.7 1.3 0.2]
 [4.6 1.5 0.2]
 [5.  1.4 0.2]]


E:\Anaconda3\envs\sklearn\lib\site-packages\sklearn\ensemble\forest.py:248: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.
  "10 in version 0.20 to 100 in 0.22.", FutureWarning)
E:\Anaconda3\envs\sklearn\lib\site-packages\sklearn\ensemble\forest.py:248: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.
  "10 in version 0.20 to 100 in 0.22.", FutureWarning)

Embedded

使用SelectFromModel選擇特征 (Feature selection using SelectFromModel)

from sklearn.feature_selection import SelectFromModel
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_iris

clf = RandomForestClassifier()
iris=load_iris()
X,y=iris.data,iris.target
print(X.shape)
print(X[:5, :])
sfm = SelectFromModel(clf, threshold=0.25)
X_sfm = sfm.fit_transform(X,y)
print(X_sfm.shape)
print(X_sfm[:5, :])

(150, 4)
[[5.1 3.5 1.4 0.2]
 [4.9 3.  1.4 0.2]
 [4.7 3.2 1.3 0.2]
 [4.6 3.1 1.5 0.2]
 [5.  3.6 1.4 0.2]]
(150, 2)
[[1.4 0.2]
 [1.4 0.2]
 [1.3 0.2]
 [1.5 0.2]
 [1.4 0.2]]


E:\Anaconda3\envs\sklearn\lib\site-packages\sklearn\ensemble\forest.py:248: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.
  "10 in version 0.20 to 100 in 0.22.", FutureWarning)

基於L1的特征選擇 (L1-based feature selection)

使用L1范數作為懲罰項的線性模型(Linear models)會得到稀疏解：大部分特征對應的系數為0。當你希望減少特征的維度以用於其它分類器時，可以通過 feature_selection.SelectFromModel 來選擇不為0的系數。

特別指出，常用於此目的的稀疏預測模型有 linear_model.Lasso（回歸）， linear_model.LogisticRegression 和 svm.LinearSVC（分類）

from sklearn.feature_selection import SelectFromModel
from sklearn.svm import LinearSVC
print(X.shape)
lsvc = LinearSVC(C=0.01, penalty="l1", dual=False).fit(X,y)
model = SelectFromModel(lsvc, prefit=True)
X_embed = model.transform(X)
X_embed.shape

(150, 4)





(150, 3)

X_embed[:5,:]

array([[5.1, 3.5, 1.4],
       [4.9, 3. , 1.4],
       [4.7, 3.2, 1.3],
       [4.6, 3.1, 1.5],
       [5. , 3.6, 1.4]])

那么工作中我們更傾向於使用什么方法呢？

首先來回顧一下我們在業務中的模型會遇到什么問題。

• 模型效果不好
• 訓練集效果好，跨時間測試效果不好
• 跨時間測試效果也好，上線之后效果不好
• 上線之后效果還好，幾周之后分數分布開始下滑
• 一兩個月內都比較穩定，突然分數分布驟降
• 沒有明顯問題，但模型每個月逐步失效

然后我們來考慮一下業務所需要的變量是什么。

• 變量必須對模型有貢獻，也就是說必須能對客群加以區分
• 邏輯回歸要求變量之間線性無關
• 邏輯回歸評分卡也希望變量呈現單調趨勢 （有一部分也是業務原因，但從模型角度來看，單調變量未必一定比有轉折的變量好）
• 客群在每個變量上的分布穩定，分布遷移無可避免，但不能波動太大

為此我們從上述方法中找到最貼合當前使用場景的幾種方法。

import pandas as pd
import numpy as np
df_train = pd.read_csv('train.csv')
df_train.head()

	PassengerId	label	Pclass	Name	Sex	Age	SibSp	Parch	Ticket	Fare	Cabin	Embarked
0	1	0	3	Braund, Mr. Owen Harris	male	22.0	1	0	A/5 21171	7.2500	NaN	S
1	2	1	1	Cumings, Mrs. John Bradley (Florence Briggs Th...	female	38.0	1	0	PC 17599	71.2833	C85	C
2	3	1	3	Heikkinen, Miss. Laina	female	26.0	0	0	STON/O2. 3101282	7.9250	NaN	S
3	4	1	1	Futrelle, Mrs. Jacques Heath (Lily May Peel)	female	35.0	1	0	113803	53.1000	C123	S
4	5	0	3	Allen, Mr. William Henry	male	35.0	0	0	373450	8.0500	NaN	S

1）變量重要性

• IV值
• 卡方檢驗
• 模型篩選

這里我們使用IV值或者模型篩選多一點

IV其實就是在WOE前面加上一項。

• \(p_{y_i}=\frac{y_i}{y_T}\)
• \(p_{n_i}=\frac{n_i}{n_T}\)
• \(woe_i = ln(\frac{p_{y_i}}{p_{n_i}})\)
• \(iv_i = (p_{y_i} - p_{n_i}) \times woe_i\)

最后只需要將每一個區間的iv加起來就得到總的iv值：

\[IV = \sum iv_i \]

import math
a = 0.4  
b = 0.6
iv = (a - b) * math.log(a / b)
iv

0.08109302162163284

或者集成模型輸出特征重要性：

# lightGBM中的特征重要性
feature = pd.DataFrame(
            {'name' : model.booster_.feature_name(),
            'importance' : model.feature_importances_
          }).sort_values(by =  ['importance'],ascending = False)

2）共線性

• 相關系數 COR
• 方差膨脹系數 VIF

在做很多基於空間划分思想的模型的時候，我們必須關注變量之間的相關性。單獨看兩個變量的時候我們會使用皮爾遜相關系數。

df_train.corr()

	PassengerId	label	Pclass	Age	SibSp	Parch	Fare
PassengerId	1.000000	-0.005007	-0.035144	0.036847	-0.057527	-0.001652	0.012658
label	-0.005007	1.000000	-0.338481	-0.077221	-0.035322	0.081629	0.257307
Pclass	-0.035144	-0.338481	1.000000	-0.369226	0.083081	0.018443	-0.549500
Age	0.036847	-0.077221	-0.369226	1.000000	-0.308247	-0.189119	0.096067
SibSp	-0.057527	-0.035322	0.083081	-0.308247	1.000000	0.414838	0.159651
Parch	-0.001652	0.081629	0.018443	-0.189119	0.414838	1.000000	0.216225
Fare	0.012658	0.257307	-0.549500	0.096067	0.159651	0.216225	1.000000

import seaborn as sns
sns.set(color_codes=True)
np.random.seed(sum(map(ord, "distributions")))
# 在數據集中繪制成對關系
sns.pairplot(df_train) # 對角線上是單維度分布

<seaborn.axisgrid.PairGrid at 0xe7abfa160>

在多元回歸中，我們可以通過計算方差膨脹系數VIF來檢驗回歸模型是否存在嚴重的多重共線性問題。定義：

\[VIF = \frac{1}{1-R^2} \]

其中，\(R_i\)為自變量 對其余自變量作回歸分析的負相關系數。方差膨脹系數是容忍度\(1-R^2\)的倒數。

方差膨脹系數VIF越大，說明自變量之間存在共線性的可能性越大。一般來講，如果方差膨脹因子超過10，則回歸模型存在嚴重的多重共線性。又根據Hair(1995)的共線性診斷標准，當自變量的容忍度大於0.1，方差膨脹系數小於10的范圍是可以接受的，表明白變量之間沒有共線性問題存在。

VIF函數詳細使用方法可以看statsmodels官方文檔.

from statsmodels.stats.outliers_influence import variance_inflation_factor
import numpy as np

data = [[1,2,3,4,5],
        [2,4,6,8,9],
        [1,1,1,1,1],
       [2,4,6,4,7]]
X = np.array(data).T

variance_inflation_factor(X,0)

98.33333333333381

3）單調性

• bivar圖

# 等頻切分
df_train.loc[:,'fare_qcut'] = pd.qcut(df_train['Fare'], 10)
print(df_train.head())
df_train = df_train.sort_values('Fare')
alist = list(set(df_train['fare_qcut']))
badrate = {}
for x in alist:
    
    a = df_train[df_train.fare_qcut == x]
    
    bad = a[a.label == 1]['label'].count()
    good = a[a.label == 0]['label'].count()
    
    badrate[x] = bad/(bad+good)
badrate

     PassengerId  label  Pclass                             Name   Sex   Age  \
271          272      1       3     Tornquist, Mr. William Henry  male  25.0   
277          278      0       2      Parkes, Mr. Francis "Frank"  male   NaN   
263          264      0       1            Harrison, Mr. William  male  40.0   
597          598      0       3              Johnson, Mr. Alfred  male  49.0   
302          303      0       3  Johnson, Mr. William Cahoone Jr  male  19.0   

     SibSp  Parch  Ticket  Fare Cabin Embarked       fare_qcut  
271      0      0    LINE   0.0   NaN        S  (-0.001, 7.55]  
277      0      0  239853   0.0   NaN        S  (-0.001, 7.55]  
263      0      0  112059   0.0   B94        S  (-0.001, 7.55]  
597      0      0    LINE   0.0   NaN        S  (-0.001, 7.55]  
302      0      0    LINE   0.0   NaN        S  (-0.001, 7.55]  





{Interval(39.688, 77.958, closed='right'): 0.5280898876404494,
 Interval(14.454, 21.679, closed='right'): 0.42045454545454547,
 Interval(7.55, 7.854, closed='right'): 0.2988505747126437,
 Interval(8.05, 10.5, closed='right'): 0.23076923076923078,
 Interval(10.5, 14.454, closed='right'): 0.42857142857142855,
 Interval(77.958, 512.329, closed='right'): 0.7586206896551724,
 Interval(-0.001, 7.55, closed='right'): 0.14130434782608695,
 Interval(27.0, 39.688, closed='right'): 0.37362637362637363,
 Interval(7.854, 8.05, closed='right'): 0.1792452830188679,
 Interval(21.679, 27.0, closed='right'): 0.5168539325842697}

f = zip(badrate.keys(),badrate.values())
f = sorted(f,key = lambda x : x[1],reverse = True )
badrate = pd.DataFrame(f)
badrate

	0	1
0	(77.958, 512.329]	0.758621
1	(39.688, 77.958]	0.528090
2	(21.679, 27.0]	0.516854
3	(10.5, 14.454]	0.428571
4	(14.454, 21.679]	0.420455
5	(27.0, 39.688]	0.373626
6	(7.55, 7.854]	0.298851
7	(8.05, 10.5]	0.230769
8	(7.854, 8.05]	0.179245
9	(-0.001, 7.55]	0.141304

badrate.columns = pd.Series(['cut','badrate'])
badrate = badrate.sort_values('cut')
print('===============================================')
print(badrate)
badrate.plot('cut','badrate')

===============================================
                 cut   badrate
9     (-0.001, 7.55]  0.141304
6      (7.55, 7.854]  0.298851
8      (7.854, 8.05]  0.179245
7       (8.05, 10.5]  0.230769
3     (10.5, 14.454]  0.428571
4   (14.454, 21.679]  0.420455
2     (21.679, 27.0]  0.516854
5     (27.0, 39.688]  0.373626
1   (39.688, 77.958]  0.528090
0  (77.958, 512.329]  0.758621





<matplotlib.axes._subplots.AxesSubplot at 0xe018a4550>

4）穩定性

• PSI
• 跨時間交叉檢驗

跨時間交叉檢驗

就是將樣本按照月份切割，一次作為訓練集和測試集來訓練模型，取進入模型的變量之間的交集，但是要小心共線特征！

解決方法

• 不需要每次都進入模型，大部分都在即可
• 先去除共線性（這也是為什么集成模型我們也會去除共線性）

群體穩定性指標(population stability index)

公式：

\[PSI = \sum{(實際占比-預期占比)*{\ln(\frac{實際占比}{預期占比})}} \]

來自知乎的例子：
比如訓練一個logistic回歸模型，預測時候會有個概率輸出p。
你測試集上的輸出設定為p1吧，將它從小到大排序后10等分，如0-0.1,0.1-0.2,......。
現在你用這個模型去對新的樣本進行預測，預測結果叫p2,按p1的區間也划分為10等分。
實際占比就是p2上在各區間的用戶占比，預期占比就是p1上各區間的用戶占比。
意義就是如果模型跟穩定，那么p1和p2上各區間的用戶應該是相近的，占比不會變動很大，也就是預測出來的概率不會差距很大。
一般認為psi小於0.1時候模型穩定性很高，0.1-0.25一般，大於0.25模型穩定性差，建議重做。

def var_PSI(dev_data, val_data):
    dev_cnt, val_cnt = sum(dev_data), sum(val_data)
    if dev_cnt * val_cnt == 0:
        return None
    PSI = 0
    for i in range(len(dev_data)):
        dev_ratio = dev_data[i] / dev_cnt
        val_ratio = val_data[i] / val_cnt + 1e-10
        psi = (dev_ratio - val_ratio) * math.log(dev_ratio/val_ratio)
        PSI += psi
    return PSI

注意分箱的數量將會影響着變量的PSI值。

PSI並不只可以對模型來求，對變量來求也一樣。只需要對跨時間分箱的數據分別求PSI即可。