数据预处理基本方法

发表于 | 分类于
本文介绍进行数据挖掘和机器学习之前需要对数据进行的基本预处理方法[1]，包括量纲标准化、有偏数据的检测和转换、离群点的检测和处理、降维、特征选择、缺失值处理、去除或者增加predictor、样本类别失衡处理，针对不同的数据挖掘任务和机器学习模型需要采用相应的预处理方法。采用UCI上的adult数据集作为例子分析。

量纲标准化

问题描述

如果一个数据集中包括了工资和年龄这两个预测特征(predictor)，如果分别用”元”和”周岁”量纲来记录这两个预测特征的数据，则前者该列的数据的平均值在几千甚至几万这个水平，数据跨度可能是4000~40000;而后者的平均值一般在几十这个水平，数据跨度可能是20～100.两厢比较，在特征空间来看工资和年龄这两个维度的话，前者的平均值远远大于后者，而且前者分布的方差远远大于后者的方差，这种预测特征的平均值和方差差异可能会造成模型的性能减弱，导致模型不稳定

处理方法

标准化(standardization)：centering + scaling;
其中centering是指预测特征的每个值减去该预测特征列的平均值;scaling是指预测特征的每个值除以该预测特征列的方差

优缺点

能够提升计算上的数值稳定性，改善模型(比如偏线性回归)的性能，但是scale操作改变了预测特征原来的量纲，因此会让模型的可解释性变差。

例子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
library(data.table)
# 将adult训练数据集加载到R的工作环境中
data <- as.data.frame(fread('./adult.data',na.strings='?',stringsAsFactors=T))
# 增加每一列的列名
colnames(data) <- c("age","workclass","fnlwgt","education","education.num","marital.status","occupation",
		"relationship","race","sex","capital.gain","capital.loss","hours.per.week","native.country","output")
# 划分为训练数据和输出类别两部分
train_output <- data[['output']]
train_data <- data
train_data$output <- NULL
# 区分数值型特征和类别型特征
col_numeric <- colnames(train_data)[sapply(train_data, is.numeric)]
col_category <- colnames(train_data)[sapply(train_data, is.factor)]
> col_numeric
[1] "age"            "fnlwgt"         "education.num"  "capital.gain"  
[5] "capital.loss"   "hours.per.week"
> col_category
[1] "workclass"      "education"      "marital.status" "occupation"    
[5] "relationship"   "race"           "sex"            "native.country"
# 查看每一列的大致分布情况，注意类别特征必须是因子类型才能够进行统计
> summary(train_data)
      age                   workclass         fnlwgt       
 Min.   :17.00   Private         :22696   Min.   :  12285  
 1st Qu.:28.00   Self-emp-not-inc: 2541   1st Qu.: 117827  
 Median :37.00   Local-gov       : 2093   Median : 178356  
 Mean   :38.58   NA              : 1836   Mean   : 189778  
 3rd Qu.:48.00   State-gov       : 1298   3rd Qu.: 237051  
 Max.   :90.00   Self-emp-inc    : 1116   Max.   :1484705  
                 (Other)         :  981                    
        education     education.num                 marital.status 
 HS-grad     :10501   Min.   : 1.00   Divorced             : 4443  
 Some-college: 7291   1st Qu.: 9.00   Married-AF-spouse    :   23  
 Bachelors   : 5355   Median :10.00   Married-civ-spouse   :14976  
 Masters     : 1723   Mean   :10.08   Married-spouse-absent:  418  
 Assoc-voc   : 1382   3rd Qu.:12.00   Never-married        :10683  
 11th        : 1175   Max.   :16.00   Separated            : 1025  
 (Other)     : 5134                   Widowed              :  993  
           occupation           relationship                   race      
 Prof-specialty :4140   Husband       :13193   Amer-Indian-Eskimo:  311  
 Craft-repair   :4099   Not-in-family : 8305   Asian-Pac-Islander: 1039  
 Exec-managerial:4066   Other-relative:  981   Black             : 3124  
 Adm-clerical   :3770   Own-child     : 5068   Other             :  271  
 Sales          :3650   Unmarried     : 3446   White             :27816  
 Other-service  :3295   Wife          : 1568                             
 (Other)        :9541                                                    
     sex         capital.gain    capital.loss    hours.per.week 
 Female:10771   Min.   :    0   Min.   :   0.0   Min.   : 1.00  
 Male  :21790   1st Qu.:    0   1st Qu.:   0.0   1st Qu.:40.00  
                Median :    0   Median :   0.0   Median :40.00  
                Mean   : 1078   Mean   :  87.3   Mean   :40.44  
                3rd Qu.:    0   3rd Qu.:   0.0   3rd Qu.:45.00  
                Max.   :99999   Max.   :4356.0   Max.   :99.00  
                                                                
       native.country 
 United-States:29170  
 Mexico       :  643  
 NA           :  583  
 Philippines  :  198  
 Germany      :  137  
 Canada       :  121  
 (Other)      : 1709

从结果可以看到，6个数值型特征中，”age”,”education.num”,”hours.per.week”这三个特征的平均分布水平与”fnlwgt”的平均分布相差巨大，另外后者与另外的两个特征”capital.loss”,”capital.gain”的分布离散程度远远大于前三个特征，因此可能在建模之前需要对这六个数值型特征进行标准化处理。

有偏数据的检测和处理

问题描述

无偏的预测特征是指其取值的分布是大致对称的，一种理想的无偏分布的例子就是正态分布。与无偏相反的有偏有两种情况：大部分值集中在数值偏小的一端，只有小部分值分布在数值偏大的一端，称为正偏（右偏）;与之相反的分布称为负偏（左偏）。

检测与处理方法

一种简单的检测有偏的方法是查看该预测特征列的直方图，可以通过调整组距来探寻数据分布更多的细节信息，但是这种方法在准确度和效率方面较差。另一种方法是采用统计上的偏度值统计量来量化偏度的程度，偏度值越接近0,分布越对称;偏度值越大于0,正偏越严重;偏度值越小于0,负偏越严重。

两种处理方法：BoxCox转换，只适用于正数的预测特征列值; Yeo-Johnson方法，能够适用于一般的数值列

例子

1
2
3
4
5
library(e1071)
skewValues <- apply(train_data[col_numeric], 2, skewness)
> skewValues
           age         fnlwgt  education.num   capital.gain   capital.loss    hours.per.week 
     0.5586919      1.4468468     -0.3116472     11.9527463      4.5942058         0.2276216

可以看到，”fnlwgt”存在轻度正偏，”capital.gain”和”capital.loss”的正偏比较严重，其余三个数值型特征的分布偏度较弱

1
2
3
4
5
library(car)
yj_data <- yjPower(train_data[['capital.gain']],lambda=0)  #lambda一般取值[-2, 2] 
new_skew <- skewness(yj_data)
> new_skew
[1] 3.095858

我们选择偏度最严重的”capital.gain”进行去偏处理，因为该特征取值包含了0，而BoxCox要求转换的值列大于0，
因此只能用Yeo-Johnson方法。处理后的新的特征”yj_data”的偏度严重减弱不少。

离群点的检测和处理

问题描述

离群点的非正式定义是：datas that exceptionally far from the mainstream of the data[1]。一般需要针对具体的数据和任务对其中的”exceptionally”和”mainstream”进行界定。一些很明显的scientifically unvalid的预测特征值(比如年龄为负值)可以直接删除。在一些数据样本数量不够大且总体分布本身有偏时，处于尾部的疑似离群点的样本点可能也是来自于该总体分布的tail部分，这是因为总体分布的尾部部分采集样本时被采集到的概率很小，如果采样次数太少，则有可能来自于tail部分的样本点会比较少，从而看起来像是离群点，对于这些样本不应该删除，而是减弱它们的影响同时参与模型的构建。

检测方法

最基本的离群点侦测方法包括利用直方图、箱线图、DBSCAN聚类来定位可能的离群点。箱线图假设数据符合高斯分布，通过计算$Q_3 + 1.5 \times IQR$和$Q_1 - 1.5 \times IQR$(其中$IQR = Q_3 - Q_1$表示第三个四分位点与第一个四分位点的差)，大于前者或者小于后者的数据均被列为可能的离群点，如果用在不同类中，则为了得到离群点的排名，需要去除不同类的分布差异所带来的影响，因此除以对应类下的IQR值。箱线图中定位离群点时选择了分位点来衡量属性的分布集中情况(第二个四分位点，即中位数)和散步情况(IQR)，是因为它们相对于平均值和方差来说对于异常点具有较好的鲁棒性。

上面所说的基本的离群点的检测方法只能够用在特征较少的情况下。在特征数较大的数据中常用的离群点的检测方法包括Robust Covariance、One_Class SVM、局部离群点检测LOF算法以及Isolation Forest。

Robust Convariance

这种检测异常点的大概思路是通过一种具备鲁棒性的方式求出符合inlier数据分布的参数，然后计算所有样本点到该分布的Mahalanobis距离，作为异常点的衡量标准。Robust Convariance方法在高维特征的情况下效果会变差，实际应用中一般需要满足：样本数大于两倍的特征数。scikit-learn实现了高斯分布情况下的Robust Convariance方法：sklearn.covariance.EllipticEnvelope。

One_Class SVM

严格来说One_class SVM属于Novelty Detection(异常点检测：训练集不包含离群点，检测一个新的样本是否异常)。该方法在embedding的p-维空间中学习一个闭环边界，确定一个包含初始样本点(不含离群点)的子空间，对于新的样本确定其是否在这个闭环所包的子空间内，如果在，则被认为是与初始样本集同分布的;否则，如果落在子空间外面，则我们以一定的概率认为它是异常点。sklearn.svm.OneClassSVM实现这一异常点检测算法。

Isolation Forest

孤立点检测算法iForest跟Random Forest类似，同样是由多颗子树iTree组成的森林，不同的是iTree的构建方法不同。iTree构建时从$n$个样本的数据集中均匀采样$\psi$个样本，在这个子数据集中，每次随机的选择一个属性、在该属性取值的$[min, max]$范围内随机取值分别作为切分变量和切分点，据此对$\psi$个样本划分到左右结点中：小于切分点的样本放入左结点，大于切分点的样本放入右结点。在左右结点中重复上述过程，直到满足终止条件(样本不可分或者iTree树的高度超过$\log \psi$).按照上面的步骤建好了$t$颗树之后，怎么判断一个样本是否有可能是异常点呢？对每个样本$x$计算它的anomaly score,这个得分的取值范围在$[0,1]$中，得分越靠近1,越有可能这个样本是异常点，样本x的anomaly score的具体计算公式为：
$$ s(x)=2^{-\frac{E[h(x)]}{c(\psi)}} $$
其中$h(x)$是样本x在第h棵树中检索到的结点的高度，$E[\cdot]$是对样本x在全部t棵树中检索到高度求平均值;$c(\psi)$是在$\psi$个样本上构建的t颗树的平均路径长度。在高维属性情况下，每个iTree分别采用Kurtosis test(峰度：描述属性取值分布形态的陡缓程度的统计量)对属性进行筛选，只考虑那些容易引起空间位置突变的属性(文本分类时还可以采用卡方检验，结合那些对分类有意义的特征[9])。

在应用上，可以用于欺诈检测，及时发现异常交易、异常用户和异常信息流，另外也可以用在为文本分类器剔除错误标记的语料，为聚类任务去噪等[9]。

iForest例子

选择UCI中的SMS Spam Collection[10]数据集，该数据集包括5574条文本，分两类，747条文本为spam文本(垃圾文本)，剩余的4827条文本为ham文本，即正常的文本，我希望iForest算法能够以较高的召回率找出这些spam文本。

为了方便后面load_files函数将各条文本即对应的类别读取到Bunch对象中进行TF-IDF处理，需要先将从UCI[10]下载的文件进行处理，原始数据的每一行代表一条文本以及对应的类别，中间以制表符分隔，制表符前面是类别，后面是对应的文本。我们将属于spam类的行保存为spam目录下的单个的文本，将属于ham类的行保存为ham目录下的单个文本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import re
dir_str = '/mnt/shared/Dataset/smsspamcollection'
f = open(dir_str+'/SMSSpamCollection','r')
i = 0
for line in f.readlines():
    split = re.split(r'\t', line)
    if 'spam' == split[0]:
        with open(dir_str+'/spam/'+bytes(i)+'.txt', 'w') as write_spam:
            write_spam.write(split[1])
    else:
        with open(dir_str+'/ham/'+bytes(i)+'.txt', 'w') as write_ham:
            write_ham.write(split[1])
    i += 1
f.close()
print 'split done!'

re.split函数用于对每一行文本根据制表符进行切分，第一部分为类别，第二部分为文本正文。对原始的数据处理完之后，需要将每个文本转换为TF-IDF的向量表示形式，在此基础上再利用iForest进行异常样本检测：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
from sklearn.datasets import load_files
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.feature_extraction.text import TfidfTransformer
from sklearn.ensemble import IsolationForest
from sklearn import metrics
# 将sms spam文本转换为TF-IDF矩阵形式
dir_str = '/mnt/shared/Dataset/smsspamcollection/'
spamdata = load_files(dir_str+'splited/',
                  categories = None,
                  shuffle = True,
                  random_state = 42)
count_vect = CountVectorizer(spamdata.data, stop_words = 'english', min_df=20, max_df=0.85)
train_counts = count_vect.fit_transform(spamdata.data)
train_tfidf = TfidfTransformer().fit_transform(train_counts)
# 在TF-IDF形式的训练集上训练iForest模型并对训练集数据进行预测
clf = IsolationForest(n_estimators=300, max_samples=200, contamination=0.1)
clf.fit(train_tfidf)
outlier = clf.predict(train_tfidf)
# 将预测结果进行转换方便打印混淆矩阵
for i in range(0, outlier.shape[0]):
    if outlier[i] == 1:
        outlier[i] = 0
    else:
        outlier[i] = 1
print('Print Confusion Matrix')
print(metrics.classification_report(spamdata.target, outlier,target_names=spamdata.target_names))
>>>>>
Print Confusion Matrix
             precision    recall  f1-score   support
        ham       0.92      0.96      0.94      4827
       spam       0.64      0.48      0.55       747
avg / total       0.88      0.89      0.89      5574

上面代码中有几个地方说明一下：CountVectorizer()函数中的参数stop_words=’english’表示采取scikit-learn中内置的英文停用词表对文本进行停用词过滤。实际上因为我们选择的sms spam collection数据集中的文本大部分(4827/5574)是人们之间正常沟通的短信，因此里面包含了大量的口语和具有个人风格的缩略词，下面是一部分文本样本：

1
2
3
4
5
ham	Is that seriously how you spell his name?
ham	I‘m going to try for 2 months ha ha only joking
ham	So ü pay first lar... Then when is da stock comin...
ham	Aft i finish my lunch then i go str down lor. Ard 3 smth lor. U finish ur lunch already?
ham	Ffffffffff. Alright no way I can meet up with you sooner?

因此设置min_df=20表示如果一个词在少于20个文本中出现，则过滤掉这个词;max_df=0.85表示如果一个词在超过85%的文本中出现过，则被认为不具有区分不同类别文本的信息，同样过滤掉。

iForest对训练集预测得到的原始outlier值只包括1和-1两个取值，1代表iForest对该文本预测为非异常的，-1则代表iForest认为该文本是异常的。事实上，iForest对于每个文本样本计算的原始的normal score值可以通过clf.decision_function(train_tfidf)得到，这个normal score值的取值在[0,1]区间，那么predict()函数根据什么来判断一个样本是否是异常的呢？很显然需要一个判断normal score的阈值来进行比较。这个阈值的计算通过IsolationForest()函数的contamination参数确定，这个参数指定的是数据集中异常点的比例，从而通过这个比例值确定decision_function()返回的normal score列表中对应的百分位值，即为比较的阈值，normal score值小于这个阈值的样本被标记为-1，对应spam类，通过查看spamdata.target_names可得知spam在下标为1的位置，ham在下标为0的位置。因此在上面的代码中将iForest的预测值还原为可以跟真实的类别spamdata.target进行比较的形式。下面的代码来自于IsolationForest的scikit-learn实现，是计算阈值和输出预测值1和-1的关键部分：

1
2
3
4
5
6
# from fit() function
self.threshold_ = -sp.stats.scoreatpercentile(
            -self.decision_function(X), 100. * (1. - self.contamination))
# from predict() function
is_inlier = np.ones(X.shape[0], dtype=int)
is_inlier[self.decision_function(X) <= self.threshold_] = -1

所以增大参数contamination的值能够提高spam类的召回率，但是精确率会下降，下面的结果是contamination取0.3的结果：

1
2
3
4
5
6
7
Print Confusion Matrix
             precision    recall  f1-score   support
        ham       0.96      0.77      0.85      4827
       spam       0.34      0.77      0.48       747
avg / total       0.87      0.77      0.80      5574

我们看一下iForest对于两个类别的文本预测的normal score的分布范围，scikit-learn中的decision_function用于计算样本的normal score值，这个值越低表示该样本越异常，代码和对应的值分布图如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import matplotlib.pyplot as plt
clf.fit(train_tfidf)
score = clf.decision_function(train_tfidf)
spam_index = [i for i,x in enumerate(spamdata.target) if 1==x]
ham_index = [i for i,x in enumerate(spamdata.target) if 0==x]
spam_score = score[spam_index]
ham_score = score[ham_index]
plt.hist(spam_score,facecolor='green',alpha=0.5,label='spam')
plt.hist(ham_score,facecolor='blue',alpha=0.5,label='ham')
plt.xlabel('score of spam text')
plt.ylabel('Count')
plt.legend(loc='upper left')
plt.grid(True)
plt.show()


从normal score值的分布来看，iForest并没有很好的把spam类文本和ham类文本区分开来，中间存在较大的重叠部分，从而在我们增大参数contamination的值使得spam类的召回率改善时会让spam类的精确率极低，因为此时重叠部分大量的ham类样本被错分为spam类。

Local Outlier Factor

LOF算法是基于密度的离群点检测算法，它基于这样一个观察：The density around an outlier object is significantly different from the density around its neighbors[11]，即离群点附近的密度与其邻居的邻近密度差距很大。因此LOF计算一个数据点对它的邻居数据点的相对密度(LOF值)来作为它是否为离群点的度量方式，这个相对密度越大，说明该数据点越有可能是离群点。

LOF算法描述如下：
输入：输入数据D，k(选择的邻居数)，距离度量方式T；
输出：根据LOF值从大到小排序的样本标号；
过程：

  1. 根据距离度量方式T计算D中两两数据点之间的距离$dist(d_i, d_j)$;
  2. 根据设定的k值计算数据点$o$的k-distance:$dist_k (o)$，这是该数据点到其第k个最近邻居之间的距离；
  3. 计算数据点$o$的k-distance邻居$N_k (0)=\{ \hat o | \hat o in D, dist(o,\hat o) \le dist_k(o)\}$，即找到那些距离数据点不大于其与第k个最近邻的距离的其他数据点(可能不止k个，因为其中可能有多个数据点到o的距离相同)；
  4. 计算数据点$o$的local reachability density: $lrd_k (o)={\|N_k(o) \|} \over {\sum_{\hat o \in N_k(o)} reachdist_k (\hat o \leftarrow o)}$，其中分母部分称为从$\hat o$到$o$的reachability distance：$max \{dist_k(o), dist(o,\hat o) \}$;
  5. 计算数据点$o$的LOF值：$LOF_k(o)={{\sum_{\hat o \in N_k(o) {{lrd_k(\hat o)} \over {lrd_k(o)}}} \over {\| N_k(o)\|}}$，可以看到数据点$o$LOF值的计算实际为它的k-distance邻居的local reachability对$o$的比值的平均值；
  6. 对于数据集D中的所有数据点重复过程2-5,然后根据LOF值从大到小对数据点进行排序输出。

处理方法

显式的离群点处理方法spatial sign，它的思想是对于样本的每一个属性除以该样本在P-维属性空间中距离原点的欧式距离，公式如下。几何上来看，相当于把所有的样本点映射到p-维属性空间的单位球面上，缓解离群点的影响。
$$ x_{ij}^\star = {x_{ij} \over \sqrt{\sum_{j=1}^P x_{ij}^2}} \tag{0} $$

某些模型本身具有对离群点较好的鲁棒性，包括基于回归树的模型、软间隔SVM(只选择少数的支撑向量建模)、逻辑回归模型(sigmoid非线性变换缓解离群点的影响)。

降维和特征选择

问题描述

降维和特征选择是有区别的，降维是指将p-维的特征空间减小到k-维同时尽量保留原来特征空间的某种信息，而特征选择是指去除一些无用的或者信息量太少的特征，保留具有高信息量的强特征。特征选择往往会造成降维的效果，但是降维不意味着进行了特征选择。

方法

PCA

原理

PCA是降维方法，严格来说不具有特征选择的作用，因为它的主要思想是通过原始的预测特征的线性组合构造新的主成分)，构造的原则是找到一组两两正交的方向形成一个新的特征坐标系，然后将原始的特征向量在这个新的坐标系下进行投影，从而得到一组新的主成分，如果新的坐标系的维度小于原始的特征空间坐标系，则此时PCA具有降维的作用。关键就在于如何找到这一组两两正交的方向。

PCA寻找这组方向的标准是在与已经找到的方向正交的前提下尽可能多的capture到原始预测特征中的方差(variability)，示意图[2]如下：

可以看到，The largest principal component is the direction that maximazes the variance of the projected data[2].为了寻找到这样一组方向，需要借助SVD分解。

SVD(奇异值分解)是将一个矩阵分解为其他三个矩阵相乘的形式，核心表达式为：

其中矩阵$U,V$的列向量分别由下面两个正定矩阵的特征向量组成：

这里重要的一点是这些特征向量在$U,V$中的排列顺序，它们的排列顺序取决于对应特征值的大小。对应(2)中的两个正定矩阵，分别存在一组特征值$\lambda_1^1,\lambda_2^1,\ldots,\lambda_m^1$和$\lambda_1^2,\lambda_2^2,\ldots,\lambda_n^2$，如果将这两组特征值从大到小排序，则前$min(m,n)$个特征值是一样的且均大于零，后面的$max(m,n)-min(m,n)$个特征值均为0,因此只需要按照其中一组特征值的大小来排列特征向量在$U,V$中的顺序。剩下的矩阵$D_{m \times n}$的主对角线上的值就是从大到小排列后的$min(m,n)$个特征值的平方根(称为A的奇异值)。

那么如何从SVD联系到PCA呢？按照奇异值从大到小的顺序我们取前$r \lt min(m,n)$个奇异值和对应的奇异向量，根据式(1)有下面的等式近似成立：

约等式两边同时右乘$V_{n \times r}$，又因为$V_{r \times n}^T V_{n \times r} = I_{r \times r}$，因此有：

如果$A_{m \times n}$表示有m个样本n个预测特征的数据集，则通过式(4)的变换将预测特征空间的维度从n减小到r。r的选择可以通过交叉验证来实现或者通过scree plot来选择。

对比PCA的降维思想和SVD，可以发现，PCA所希望寻找那一组两两正交的方向就是A矩阵的一组右奇异向量，而A的右奇异值用来衡量原始数据集在对应的右奇异向量上的方差。

优劣势

PCA创建的主成分特征之间是线性无关的，适合于对特征弱相关性有要求的学习模型(比如朴素贝叶斯、线性回归模型)；同时PCA也能够用于去除数据集中无关的冗余信息，在图像存储中可以用较小的存储空间来保存图片同时保留较好的图像还原度。

从PCA降维的过程可以发现，PCA是个无监督的降维过程，也就是说它在构建这些主成分的时候并没有考虑到它们与输出变量之间的关系;另外主成分的构建标准是投影方差最大化，但是一个预测特征所包含的信息可能无法通过方差来很好的衡量，如果一个与输出变量关系不大的预测特征的方差非常大，则经过PCA降维处理后，该预测特征可能会在多个主成分中占据较大的分量(loading)，为了克服这种情况的出现，最好对有偏的预测特征进行去偏处理，并对所有的特征进行标准化处理。关于PCA更具体和详细的分析和推导参见参考文献[4].

PLS

偏最小二乘跟PCA类似的地方在于它所寻找的PLS方向同样是原来预测特征的线性组合，不同的地方在于确定PLS方向时的准则既做到尽可能的提取原来预测特征的方差信息，又做到所确定的PLS方向是与输出变量的相关性最大。同样，PLS需要对数据集(包括输出变量)进行标准化处理，具体的算法步骤参考了文献[3],下面对算法过程进行简单梳理。

  • 输入：数据集预测变量部分$E_0$,数据集输出变量部分$Y$，假设包括m个样本，n个预测变量，则矩阵$E_0$维度为$m \times n$，只考虑一个输出变量，即$Y_{m \times 1}$，预测特征集$X=(x_1,\ldots,x_n)$.
  • 输出：r个PLS成分$(t_1,t_2,\ldots,t_r)$
  1. 求矩阵$E_0^T Y Y^T E_0$最大特征值所对应的特征向量$w_1$,求得成分$t_1 = w_1^T X$,计算成分得分向量$\hat t_1 = E_0 w_1$和残差矩阵$E_1 = E_0 - \hat {t_1} \alpha_1^T$，其中$\alpha_1 = {{E_0^T \hat t_1} \over {\| \hat t_1 \|}^2}$;
  2. 求矩阵$E_1^T Y Y^T E_1$最大特征值所对应到特征向量$w_2$，求得成分$t_2 = w_2^T X$，计算成分得分向量$\hat t_2 = E_1 w_2$和残差矩阵$E_2 = E_1 - \hat {t_2} \alpha_2^T$，其中$\alpha_2 = {{E_1^T \hat t_2} \over {\| \hat t_2 \|}^2}$;
  3. 重复上述步骤，直到通过交叉有效性确定共抽取r个成分$t_1,t_2,\ldots,t_r$为止。

Fisher’s Discriminant Analysis

关于费雪判别分析的内容主要参考自[5]，本文只概述其原理和梳理其主要的过程，具体的数学推导在[5]中有详细和具体的论述。最原始的费雪判别分析是专用于二类别分类问题的降维方法，它希望能够找到一根过原点的直线，使得两个类别的数据样本在这条直线上的投影尽可能的分离开，同时使同一组内的投影点之间尽可能集中。为此，费雪提出的准则为：
$$ \max J(w) =\max {{m_2 - m_1}^2 \over {s_1^2 + s_2^2}} $$
其中分子部分称为组间散布，用于衡量投影后两个类别样本之间的分离程度;分母部分称为整体组内散布，分别代表两个类别样本投影之后各类之间的集中程度，它们的定义分别为：

将原始的$J(w)$进行改写后变为$J(w)=\frac{w^T S_B w}{w^T S_W w}$，此时最大化费雪准则等价与下列约束优化问题：

通过求解这个优化问题得到唯一的最佳投影直线的方向向量为$w \propto S_W^{-1}(m_2 - m_1)$,其中：

上式表示整体组内散布矩阵。如果需要将费雪的判别分析引入分类功能，还需要在最佳投影直线上确定一个门槛值$w_0$来进行类别判断。在确定$w_0$时候，如果假设分类样本的条件概率$p(x | c_1),p(x | c_2)$服从正态分布并且两类别的数据点具有相同的散布形态(即具有相同的协方差矩阵)，则对应的分类方法称为线性判别分析LDA，关于门槛值$w_0$的具体求解过程在[5]中有详细的推导和说明。

前面所讲的费雪判别分析限定在二分类问题上，实际上可以将其推广到多分类的问题，具体的过程在[5]中末尾处有详细的推导，多分类下的费雪判别分析算法过程如下[6]:

  • 输入：N个D维向量$x_1,\ldots,x_N$，数据能够被分成d个类别;
  • 输出：投影矩阵$W = (w_1,\ldots,w_{d-1})$,其中每一个$w_i$都是D维列向量
  1. 求出类内散度矩阵$S_w$和类间散度矩阵$S_b$;
  2. 对$S_w = UDV^T$做奇异值分解，求得$S_w^{-1} = VD^{-1}U^T$;
  3. 对矩阵$S_w^{-1} S_b$做特征值分解;
  4. 取最大的前d-1个特征值对应的特征向量$w_1,\ldots,w_{d-1}$.

上述算法之所以选择取前面d-1个特征向量，是因为k类别判别分析所能产生的有效线性特征总数至多为k-1[5]。利用这k-1个D维特征向量，就可以得到k-1个新的预测特征值列。注意利用费雪判别分析做降维处理时不需要对数据进行去中心化，避免去中心化后每个类的中心会重叠[6].

自编码器

自编码器是一种特殊的神经网络，以单个隐藏层的自编码器为例，输出层尝试使得输出$\hat x$接近于输入$x$，从而通过隐藏层学习到输入数据的压缩表示，当给定一个数据的压缩表示，结合隐藏层和输出层之间的连接权重矩阵就可以重构出最原始未压缩的数据[7].为了学习到输入数据的压缩表示，隐藏层神经元个数要求小于输入层神经元个数，或者在隐藏层神经元个数大于输入层个数时加入稀疏性限制(稀疏自编码器)，从而学习到输入数据的低维结构。如果输入数据的维度较大或者希望找到数据维度更低的压缩表示，可以采用栈式自编码器，它以单隐藏层的自编码器为基础进行连接：前一个自编码器学习到的低维结构作为下一个自编码器的输入层。

其他

其他的降维方法还包括多维缩放(MDS)(用于图片的降维[8])，基于流形学习的降维方法[6].

例子

我们利用PCA来对一张图片进行实验。

1
2
3
4
5
6
7
library(jpeg)
greyflower <- readJPEG('greyflower.jpg')  # 将灰度图像中像素点的值以矩阵的形式保存起来
svdflower <- svd(greyflower)  # 对图像矩阵进行SVD分解
sumD <- sum((svdflower$d)^2) # 计算所有奇异值平方的总和
percentD <- (svdflower$d)^2/sumD*100  # 求出每一个奇异值平方值在总和中所占的比例
plot(x=c(1:min(dim(greyflower))),y=percentD,
        xlab='主成分标号',ylab='所占奇异值总和百分比(%)',type='b',col='red')

这里的图片是灰度图，对其SVD分解后计算每个奇异值平方的占比，绘制陡坡图如下;如果是RGB彩色图，需要分别对三个颜色通道进行SVD分解，然后再还原为彩色图。

可以看到，原始图片的variability信息几乎可以完全被前面几个主成分解释。我们一次选择5、10、20、50、100和全部主成分用于灰度图还原：

1
2
3
4
singularNum <- 5   # 用于选择还原灰度图时的主成分个数
greyflowerPca <- svdflower$u[,c(1:singularNum)]%*%
        diag(svdflower$d[1:singularNum])%*%t(svdflower$v[,c(1:singularNum)])   # 根据SVD表达式而来
writeJPEG(greyflowerPca,target='greyflowerPca.jpg',quality=1)   # 在当前工作目录下生成灰度图greyflowerPca.jpg

得到的对比图如下：

我们对比一下用于还原灰度图的主成分的总特征值之和的占比分别是多少：

缺失值处理

缺失值如果集中在少数样本上或者少数预测特征上时，可以考虑将这些样本或者预测特征去掉，但是预测特征的移除需要谨慎考虑。更为一般的处理方法是进行缺失值的填充(imputation)，有两种常用的填充方法：

  1. 基于样本的K-近邻方法
    找到距离缺失值所在样本最近的k个邻居，利用这k个邻居的预测特征值的平均值或者加权平均值作为缺失值填充值。这种方法同样需要确定k的取值和选择距离的度量方式，好处是缺失值的填充值必定在该预测特征的range范围内。适用于缺失值集中在样本行上的数据集。

  2. 基于预测特征的局部建模方法
    如果缺失值集中在某些预测特征A上，且存在其他预测特征B与A有线性相关关系或非线性的关系(可以从已有的数据判断或者从数据集的实际背景得到)，则以A为输出变量，以B为预测变量建立相应的拟合模型$f:B \rightarrow A$，利用建立的拟合模型$f$对A中的缺失值进行预测性填充。

例子

1
2
3
4
5
col_NA <- colnames(train_data)[apply(train_data,2,anyNA)]
num_NA <- apply(train_data[col_NA],2,function(x){sum(is.na(x))})
> num_NA
    workclass     occupation   native.country 
     1836           1843            583

可以看到，所有的缺失值集中在以上三个类别型特征上，考虑到样本数(32561)较大，可以考虑对这些列进行值填充。caret包中的preProcess方法支持三种缺失值填充方法:基于KNN填充的、利用Bagging树填充的、利用中值填充的，我们选择用KNN填充作为例子,利用KNN进行填充时，它只会选择其中的数值型特征作为邻居的衡量方法。

1
2
knn_impute <- preProcess(train_data, method=c("knnImpute")
new_train <- predict(knn_impute, train_data)

预测特征的移除

问题说明

预测特征的移除一般是因为它没有提供有效的信息或者它提供了冗余的信息，而移除之后能够使得模型的训练更快，模型本身更简单，也能够提高模型的可解释性和稳定性。存在退化分布(degenerate distribution)的预测特征对于所有的样本提供了几乎相同的信息，因此对于模型而言没有提供具有区分度的信息，也就是说该预测特征没有提供有效信息；另一种情况是两个或者多个预测特征之间存在(共/多重共)线性关系时，相当于向模型提供了相同的信息，造成信息冗余，而且导致模型复杂度增加，稳定性变差，难以准确确定预测变量的相对重要性。下面分别对这两种情况的检测方法[1]进行说明。

共线性检测

  1. VIF(variance inflation factor)方差膨胀因子只适用于线性回归，而且要求样本数大于预测变量的个数;
  2. PCA主成分提取到的方差比例和对应的loading值(这个方法的说明后续补充);
  3. 预测变量之间的协相关系数矩阵：这是比较常用的方法

在实际中用的比较多的是第三种方法，尽可能少的移除预测变量使得剩余的预测变量之间的相关性低于某个设定的阈值，步骤如下：

  1. 计算所有预测变量之间的相关系数;
  2. 找出相关系数绝对值最大的系数所对应的两个预测变量A和B;
  3. 计算预测变量A与其他预测变量之间相关系数绝对值的平均值$\overline {cor}_A$,对B进行同样的计算得到$\overline {cor}_B$;
  4. 如果$\overline {cor}_A \gt \overline {cor}_B$，则移除预测变量A;否则移除B;
  5. 重复步骤2-4直到所有剩余预测变量之间的相关系数绝对值小于设定的阈值。

退化分布的预测特征的检测

一个预测变量产生退化分布是指该预测变量样本取值相同或者绝大部分相同，前者称为Zero Variance Predictor，后者称为Near-zero Variance Predictor。

判断一个预测变量是否是Near-zero Variance Predictor，一般有两个条件：

  1. 不同值的个数占所有取值个数的比例很低(比如10%);
  2. 出现次数最多的那个值的频数除以出现次数第二多的那个值的频数，这个比值很大(比如大于20)。

例子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
index_col_category <- sapply(col_category,simplify=T,function(x){grep(x,colnames(train_data))})
index_category <- unlist(index_col_category,use.names=F)
> cor(train_data[-index_category])
                       age        fnlwgt  capital.gain capital.loss
age             1.00000000 -0.0766458679  0.0776744982   0.05777454
fnlwgt         -0.07664587  1.0000000000  0.0004318858  -0.01025171
capital.gain    0.07767450  0.0004318858  1.0000000000  -0.03161506
capital.loss    0.05777454 -0.0102517117 -0.0316150630   1.00000000
hours.per.week  0.06875571 -0.0187684906  0.0784086154   0.05425636
               hours.per.week
age                0.06875571
fnlwgt            -0.01876849
capital.gain       0.07840862
capital.loss       0.05425636
hours.per.week     1.00000000

上面的代码是查看数值型特征之间的相关性，从结果可以看出，基本上不存在线性相关关系，而且特征维度不高，因此没必要用PCA来构建线性无关的主成分特征。我们来检查一下是否存在呈现退化分布的特征。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 10表示不同值对所有取值的占比阈值；30表示the ratio of the most common value to the second most common value
nzv <- nearZeroVar(train_data,uniqueCut=10,freqCut=30,saveMetrics=T)
> nzv
                freqRatio percentUnique zeroVar   nzv
age              1.011261   0.224194589   FALSE FALSE
workclass        8.931917   0.024569270   FALSE FALSE
fnlwgt           1.000000  66.484444581   FALSE FALSE
education        1.440269   0.049138540   FALSE FALSE
education.num    1.440269   0.049138540   FALSE FALSE
marital.status   1.401853   0.021498111   FALSE FALSE
occupation       1.010002   0.042996222   FALSE FALSE
relationship     1.588561   0.018426952   FALSE FALSE
race             8.903969   0.015355794   FALSE FALSE
sex              2.023025   0.006142317   FALSE FALSE
capital.gain    86.020173   0.365467891   FALSE  TRUE
capital.loss   153.673267   0.282546605   FALSE  TRUE
hours.per.week   5.398013   0.288688922   FALSE FALSE
native.country  45.365474   0.125917509   FALSE  TRUE

从结果可以看到，”capital.gain”,”capital.loss”,”native.country”这三个特征的freqRatio值非常大，意味着出现次数最多的取值远远大于出现次数第二多的取值，其中”capital.gain”,”capital.loss”的情况更加严重，我们观察一下这两个特征的直方图：

1
2
3
4
5
6
7
8
9
library(ggplot2)
library(gridExtra)
p1 <- ggplot(train_data,aes(capital.gain))+
	  geom_histogram(aes(y=..density..))+
	  geom_density(color='red')
p2 <- ggplot(train_data,aes(capital.loss))+
		geom_histogram(aes(y=..density..))+
		geom_density(color='green')
grid.arrange(p1,p2)


可以发现，这两个特征的取值几乎全部集中在0这个值上，因此可以考虑从数据集中去除这两个特征。

类失衡问题处理

训练数据出现类失衡问题时，模型也会倾向于预测多数类，除了决策树和基于决策树的Bagging方法、Boosting方法受到类失衡数据的影响较小之外，一般需要对模型或者数据进行相应的处理。在能够进一步获取训练数据的前提下，采集更多的数据看类失衡的情况是否能够得到缓解。另外就是需要更换模型的评判方式，一般的准确率和误差率不再适用，而应该选择混淆矩阵、精确率、召回率和F1值等评价指标。

样本权重

增加少数类的样本权重，减少多数类的样本权重，类似于重复使用少数类样本。这种方法不对数据集进行抽样，而是在模型中指定样本的权重向量。

抽样方法

抽样方法的思路不是让模型取应对类失衡问题，而是通过抽样方法手动平衡各个类别的样本比例。一般有向下抽样和向上抽样两种方法。

向下抽样

从多数类别样本中抽取样本使得不同类别样本量大致相同。一种方法是从多数类样本中随机抽取直到样本量大致相同;另一种方法是对所有类使用bootstrap抽样使得bootstrap样本集中的类分布大致平衡。

向上抽样

对少数类的样本进行有放回抽样直至不同类的样本量大致相同，类似于对样本赋予不一样的权重。

截点选择

通过ROC曲线选择模型分类的最佳截点值(即预测门槛值)，避免模型在少数类上的过高误分类率。

参考文献

  1. Max Kuhn. Applied Predictive Modeling
  2. Trevor Hastie. The Elements of Statistical Learning, Data Mining, Inference, and Prediction, second Edition
  3. http://cs.ananas.chaoxing.com/download/554f2029e4b017d277eb4bf7
  4. https://ccjou.wordpress.com/2013/04/15/%E4%B8%BB%E6%88%90%E5%88%86%E5%88%86%E6%9E%90/
  5. https://ccjou.wordpress.com/2014/03/14/%E8%B2%BB%E9%9B%AA%E7%9A%84%E5%88%A4%E5%88%A5%E5%88%86%E6%9E%90%E8%88%87%E7%B7%9A%E6%80%A7%E5%88%A4%E5%88%A5%E5%88%86%E6%9E%90/
  6. http://chenrudan.github.io/blog/2016/04/01/dimensionalityreduction.html#3.3
  7. http://ufldl.stanford.edu/wiki/index.php/%E8%87%AA%E7%BC%96%E7%A0%81%E7%AE%97%E6%B3%95%E4%B8%8E%E7%A8%80%E7%96%8F%E6%80%A7
  8. https://ccjou.wordpress.com/2013/05/29/%E5%8F%A4%E5%85%B8%E5%A4%9A%E7%B6%AD%E6%A8%99%E5%BA%A6%E6%B3%95-mds/
  9. http://qf6101.github.io/machine%20learning/2015/08/01/Isolation-Forest
  10. SMS Spam Dataset:https://archive.ics.uci.edu/ml/datasets/SMS+Spam+Collection
  11. LOF Example：http://www.cse.ust.hk/~leichen/courses/comp5331/lectures/LOF_Example.pdf