1.5 Performance Evaluation

为了衡量不同模型的效果，我们通常需要在测试数据利用特定的评价标准进行评价。那么我们需要定义用什么数据进行测试，也就是如何进行数据划分；用什么指标进行模型评估，也就是模型评价指标。

1) 数据集划分

1a) Training, Test and Validation Sets

当数据量非常大时(在深度学习中比较常见)，常见的做法是将数据集划分为三部分:

• Training Set：在调整超参数的过程中用于训练模型的数据集

• 在调参的过程中用于评估超参数效果的数据集: 有的文献中称为validation set，有的文献中称为testing set，上图中标注为testing set

• 用来最终测试模型的泛化能力的数据集：有的文献中称为testing set，有的文献中称为validation set，上图中标注为validation set

• 在数据集划分的过程中，我们通常不希望某个label只在特定一个数据集中出现。一般而言这种情况出现的可能性不大，但是当样本的类别比较多，或者某些类别的样本数量特别少时，也是容易发生的。这是我们通常会采取分层划分的方式，即对于每个类别的样本，分别进行数据集划分，这样就能排除上述的可能性。

• 超参数是模型不能直接从数据中学习得到，而需要事先人为给定的参数。所以为了选择合适的超参数(即所谓的"调参")，需要用一个外部的数据集上的分类效果作为挑选的指标。调参的过程属于一种模型选择(model selection)。

• 不同领域的文献中对数据集的描述可能会有一些出入。在传统机器学习的文献中，测试集(test set)指的是在性能评估前不能见到label的样本(即最终用于性能评估的样本)，验证集(validation set)指的是用来调整超参数的数据集；但是有些临床研究中会把最终用于性能评估的样本称为"验证集"。这需要大家结合具体的语境进行区分，例如上图给出的数据集命名方式对应后者。

• 希望具体了解的同学可以参考维基百科的介绍Training, validation, and_test data sets

• 拟合最终用于评估模型泛化能力的模型时，所有可以用到的数据有时被称为"discovery set"。用来最终测试模型的泛化能力的数据集，剩下的数据都属于discovery set。在上图中，training set+testing set = discovery set

• 上图涉及了两类性能评估:

  • 一是评估超参数的分类效果:我们选择不同的超参数组合，用training set拟合模型，用上图的"testing set"的标签作为指标，挑选出最大化分类性能的超参数组合

  • 二是评估模型最终的泛化能力: 挑选出我们认为最合适的超参数组合后，我们用这组超参数在所有可用数据，即discovery set上拟合模型，再用一些从未见到过标签的数据(上图的"validation set")评估模型的泛化能力，作为分类效果的最终度量

由于我们在调参的过程中最大化的是上图中的testing set上的分类效果，所以在testing set上最好的一组超参数对应的分类效果是高估的，不能最终用于模型的性能评估。

如果我们不去调整模型的超参数，直接使用一个独立于数据的经验值，则不需要上图中的"testing set"。这时training set等价于discovery set，剩下的样本可全部用来评估模型的泛化能力

当数据量比较少时(在传统机器学习的应用场景中比较常见)，我们通常会多次重复这一过程，取平均值以降低最终汇报的结果受随机因素的影响。具体的做法可以大致分为三类，我们以下逐一进行介绍:

• 多次random split

• 交叉验证(Cross-Validation)

• bootstrapping

1b) 多次random split

相当于把我们在1a) 中介绍的方法重复很多次，每次都能得到一个评估泛化能力的指标，这样得到的分布就考虑到了随机性的影响

1c) Cross-Validation

• K-fold Cross-validation

• 在K折交叉验证中，数据集被均匀地划分为$k$个部分（folds）。在每轮验证中，把$\frac{k-1}{k}$个folds当做discovery set，在剩下的一个fold上评估模型分类效果。

• K折交叉验证确保训练样本和测试样本之间没有重叠，K轮结束后，每个样本会被设置为测试样品一次。最后，模型平均表现是在 $k$轮次中计算指标的平均值得到的。

• 如果希望在交叉验证的过程中同时对模型超参数进行调节，一种常见的做法被称为"nested cross validation"，即在每个划分出的discovery set上，对给定的每一组超参数组合，都进行一轮交叉验证，用平均的分类效果当做为这一个discovery set挑选超参数的指标，再用表现最好的超参数组合在整个discovery set上拟合模型。sklearn的文档提供了一个很好的例子，请大家参考plot nested cross validation iris

• K折交叉验同样可以重复多次，以考虑随机性的影响。

• Leave-one-out Cross-validation

Leave-one-out Cross-validation (LOOCV)是一种极端的K折交叉验证，对应把k设成样本总数的情形。这就是说每次进行只留下一个样本用作测试集，其余m-1全为训练集，进行m次训练，取m次的评估结果的平均值进行模型选择。

LOO的好处在于有确定的结果，缺点在于计算量比较大，而且无法反映分类效果受随机因素的影响。

1d) Bootstraping

当数据集较小，难以有效划分Training/Test时，可以采用**Bootstraping(自助采样)的方法。给定包含m个样本的数据集D，我们对它进行采样产生数据集 D′：每次随机从D中挑选出一个样本，将其拷贝放入D′, 然后再将该样本放回初始数据集D中，使得该样本在下次采样时仍有可能被采样到；这个过程重复执行m次后，我们就得到可包含m个样本数据的数据集D′,这就是Bootstraping(自助采样)**的结果。样本在m次采样中始终不被采到到概率为：

$$lim_{m\to\infty} (1-\frac{1}{m})^m \to \frac{1}{e} \approx 0.368$$

初始数据集D中约有36.8%的样本未出现在采样数据集D′中。于是我们可将D′ 用于模型拟合，D-D′用于评估泛化能力。

由于bootstraping进行的是有放回抽样，同一个样本在discovery set中可能会出现多次。

1e) 实现方法

简单起见，我们暂时不考虑调参的问题。

python

[sklearn] package在sklearn.model_selection子模块实现了很多种数据集划分的方式。如果理解了上述的数据集划分具体在做什么事情，直接用numpy也可以非常容易的实现。

• random split: 可使用train_test_split函数

sklearn.model_selection.train_test_split(*arrays, test_size=None, train_size=None, random_state=None, shuffle=True, stratify=None)
# train_size, test_size定义了两个数据集所占的样本比例
# random_state: 如果给定一个固定的整数作为random seed，可以保证结果的可重复性
# stratify: 用于分层划分的变量，可以是样本标签

• 多次random split: 我们可以在一个循环里多次调用train_test_split函数。除此之外，sklearn还提供了ShuffleSplit和StratifiedShuffleSplit两个类，分别针对需要分层抽样课不需要分层抽样的情形，可以简化这一过程的实现

import numpy as np
from sklearn.model_selection import StratifiedShuffleSplit,ShuffleSplit
X = np.array([[1, 2], [3, 4], [1, 2], [3, 4], [1, 2], [3, 4]])
y = np.array([0, 0, 0, 1, 1, 1])

# n_splits控制了重复次数，test_size控制测试集大小
rs = ShuffleSplit(n_splits=5, test_size=.25, random_state=0)
for train_index, test_index in rs.split(X):
    print("TRAIN:", train_index, "TEST:", test_index)
    X_train, X_test = X[train_index], X[test_index]
    y_train, y_test = y[train_index], y[test_index]    

sss = StratifiedShuffleSplit(n_splits=5, test_size=0.5, random_state=0)
# 根据样本标签进行分层划分
for train_index, test_index in sss.split(X, y):
    print("TRAIN:", train_index, "TEST:", test_index)
    X_train, X_test = X[train_index], X[test_index]
    y_train, y_test = y[train_index], y[test_index]

• 交叉验证: sklearn 提供了多个类用于实现交叉验证，具体使用请大家自行参阅文档

  • KFold: 常规的K折交叉验证的数据集划分

  • StratifiedKFold: K折交叉验证，按label分层划分fold

  • RepeatedKFold: 常规的K折交叉验证,重复多次

  • RepeatedStratifiedKFold: K折交叉验证，按label分层划分fold，重复多次

• bootstrapping: 我们可以用sklearn.utils.resample函数来实现

import numpy as np
from sklearn.utils import resample
X = np.array([[1, 2], [3, 4], [1, 2], [3, 4], [1, 2], [3, 4]])
y = np.array([0, 0, 0, 1, 1, 1])
indices = np.arange(X.shape[0])
n_repeats = 100 #重复次数
for i in range(n_repeats):
    train_index = resample(indices, replace=True, n_samples=X.shape[0], random_state=None, stratify=y)
    # replace=True 意味着有放回的抽样
    # n_samples=X.shape[0]意味着一共抽出和原始数据集大小相当的样本
    # stratify=y: 按样本标签分层抽样
    test_index = np.setdiff1d(indices,train_index)
    X_train, X_test = X[train_index], X[test_index]
    y_train, y_test = y[train_index], y[test_index]

R

R的caret package提供了很多用于数据集划分的函数，例如:

• createDataPartition: 相当于sklearn的model_selection.train_test_split或ShuffleSplit

• createFolds: 类似sklearn的StratifiedKFold

• createMultiFolds: 类似sklearn的RepeatedStratifiedKFold

• createResample: 有放回的抽样/bootstraping

这些函数接受样本标签作为输入进行分层抽样，并返回一个列表，列表中是一些数组，对应一次划分的训练集或测试集的标号。详细的说明请大家参考caret的文档。

y <- factor(c(1,1,2,2,2,2,2,2,2,2))
> y[createDataPartition(y)[[1]]]
[1] 1 2 2 2 2
Levels: 1 2
# createDataPartition总是会把把2个1中的1个划分到训练集中

• 请大家注意，R package通常通过数据类型区分分类问题和回归问题。

• 例如,如果我们们希望对分类问题根据标签进行分层划分，如果希望caret有预期的行为，就应当保证输入的样本标签是因子类型而不是数值类型

2) Performance Evaluation

2a) Confusion matrix

最常用的评估分类模型表现的方法是构建一个混淆矩阵(confusion matrix)

Confusion matrix会总结模型正确和错误分类的样本数量，并将预测的样本分成如下四类：

	真实值	真实值	总数
	Condition Positive	Condition Negative
Predicted Condition Positive	True Positive（TP）	False Positive（FP）	P'
Predicted Condition Negative	False Negative（FN）	True Negative（TN）	N'
总数	P	N

真阳性（TP）:正确的肯定，诊断为有，实际上也有癌症。

伪阳性（FP）:错误的肯定，诊断为有，实际上却没有癌症。第一型错误。

真阴性（TN）:正确的否定，诊断为没有，实际上也没有癌症。

伪阴性（FN）:错误的否定，诊断为没有，实际上却有癌症。第二型错误。

我们可以通过python函数获得TP，FP,TN,FN。

from sklearn.metrics import confusion_matrix
y_true=[0, 1, 0, 1]
y_pred=[1, 1, 1, 0]
tn, fp, fn, tp = confusion_matrix(y_true,y_pred ).ravel()
# 输出
(tn, fp, fn, tp)
(0, 2, 1, 1)

2b) Evaluation Scores

• Accuracy (0 ~ 1)

summarizes both positive and negative predictions, but is biased if the classes are imbalanced:

$\text{Accuracy} = \frac{TP + TN}{TP + TN + FP + FN}$

from sklearn.metrics import accuracy_score
y_true=[0, 1, 0, 1]
y_pred=[1, 1, 1, 0]
accuracy_score(y_true, y_pred)
# 输出
0.25

• Sensitivity (0 ~ 1)

Also called Recall, True Positive Rate (TPR). It _**_summarizes how well the model finds out positive samples:

$\text{Sensitivity/Recall/TPR} = \frac{TP}{TP + FN}$

import numpy as np
from sklearn.metrics import recall_score
y_true=[0, 1, 0, 1]
y_pred=[1, 1, 1, 0]
recall=recall_score(y_true, y_pred, average='macro')
# 输出
recall
0.25

• Specificity (0 ~ 1)

Also called True Negative Rate, which equals 1 - False Positive Rate (FPR), FPR=FP/(TN+FP)。

$\text{Specificity} = \frac{TN}{FP + TN}$

$\text{FPR} = \frac{FP}{FP + TN}=\text{1-Specificity}$

from sklearn.metrics import confusion_matrix
y_true=[0, 1, 0, 1]
y_pred=[1, 1, 1, 0]
tn, fp, fn, tp = confusion_matrix(y_true,y_pred ).ravel()
specificity=tn/(fp+tn)
# 输出
specificity
0.0

• Precision/positive Predictive Value (PPV) (0 ~ 1)

summarizes how well the model finds out negative samples:

$\text{Precision/Positive Predictive Value} = \frac{TP}{TP + FP}$

import numpy as np
from sklearn.metrics import precision_score
y_true=[0, 1, 0, 1]
y_pred=[1, 1, 1, 0]
precision=precision_score(y_true, y_pred, average='macro')
# 输出
precision
0.16666

• F1 score (0 ~ 1)

balances between positive predictive value (PPV) and true positive rate (TPR) and is more suitable for imbalanced dataset:

$\text{F1 score} = 2 \frac{PPV \cdot TPR}{PPV + TPR}$

import numpy as np
from sklearn.metrics import f1_score
y_true=[0, 1, 0, 1]
y_pred=[1, 1, 1, 0]
f1_score=f1_score(y_true, y_pred, average='macro')
# 输出
f1_score
0.2

• Matthews correlation coefficient (MCC) (-1 ~ 1)

another metric that balances between recall and precision:

$\text{MCC} = \frac{TP \times TN - FP \times FN} {(TP + FN)(TP + FP)(TN + FP)(TN + FN)}$

2c) ROC Curve

有时一个固定的cutoff不足以评估模型性能。 Receiver Operating Characterisic（ROC）曲线可以评估模型的表现。 ROC曲线是根据在不同分类阈值下计算得到TPR和FPR，然后得到一系列的成对的TPR和FPR。ROC曲线对于类别不平衡问题也有比较好的评估。

如下图，描述了True Positive Rate (TPR)和False Positive Rate (FPR)之间的关系。TPR和FPR之间是成正比的，TPR高，FPR也高。x轴是FPR，y轴是TPR。ROC曲线下面积（AUROC）是一个单值，它总结了不同截止值下的模型平均表现，常常用于报告模型的分类表现。AUROC接近于1，可以认为模型的分类效果很好。

2d) Plot ROC with python

另外我们可以通过python函数完成这个过程。示例中是通过x预测y。最后用预测的y对照真实的y进行ROC绘制。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, auc
from sklearn.linear_model  import LogisticRegression
# 数据加载
x=np.array([[12, 87, 70, 74, 46, 58, 55, 33, 60, 50, 70, 22, 68, 90, 10],
[55, 44, 23, 35, 43, 31, 40, 50, 20, 22, 50, 60, 10, 30, 40],
[74, 46, 56, 69, 59, 90, 30, 76, 34, 11, 60, 90,  9, 20, 39]])
x_=np.array([[78, 50, 70, 81, 44, 20, 51, 40, 30, 81],
[50, 60, 33, 31, 11, 56, 31, 11, 60, 13],
[23, 82, 51, 12,  5, 44, 17,  4, 57, 24]])
y=np.array([1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1])
y_=np.array([0, 1, 0, 0, 0, 1, 0, 0, 0, 0])

X_train = x.T
X_test = x_.T
y_train = y
y_test = y_

#模型训练测试
log_reg = LogisticRegression()
log_reg.fit(X_train, y_train)
# 预测
y_predict=log_reg.predict(X_test)
proba=log_reg.predict_proba(X_test).T[1]
fpr, tpr, thresholds = roc_curve(y_test, proba)
roc_auc = auc(fpr, tpr)
# 输出AUC
roc_auc
0.9375
# 画ROC曲线
plt.plot(fpr, tpr, '-', color='b', label='ROC', lw=2)
plt.plot([0, 1], [0, 1], '--', color=(0.6, 0.6, 0.6), label='Random Chance')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.title('ROC curve')
plt.xlabel('FPR')
plt.ylabel('TPR')
plt.legend()
plt.show()
plt.close()

结果如下：

2e) Plot ROC with R

我们也可以通过R完成这个过程。示例中是通过x预测y。最后用预测的y对照真实的y进行ROC绘制。

library("ROCR")
#数据加载
X_train <- t(matrix(c(12, 87, 70, 74, 46, 58, 55, 33, 60, 50, 70, 22, 68, 90, 10,
                    55, 44, 23, 35, 43, 31, 40, 50, 20, 22, 50, 60, 10, 30, 40,
                    74, 46, 56, 69, 59, 90, 30, 76, 34, 11, 60, 90,  9, 20, 39),
                  byrow = TRUE, nrow = 3))
X_test <- t(matrix(c(78, 50, 70, 81, 44, 20, 51, 40, 30, 81,
                   50, 60, 33, 31, 11, 56, 31, 11, 60, 13,
                   23, 82, 51, 12,  5, 44, 17, 4, 57, 24),
                 byrow = TRUE, nrow = 3))
Y_train <- t(matrix(c(1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1),
                  byrow = TRUE, nrow = 1))
Y_test <- t(matrix(c(0, 1, 0, 0, 0, 1, 0, 0, 0, 0),
                   byrow = TRUE, nrow = 1))
#模型训练测试
X <- as.data.frame(cbind(X_train,Y_train))
logistic <- glm(V4 ~ ., data = X, family='binomial', control=list(maxit=100))
summary(logistic)
#预测
X_test <- as.data.frame(X_test)
predicted_probs <- predict(logistic, X_test, type = 'response')
pred <- prediction(predicted_probs, Y_test)
auc <- performance(pred, 'auc')@y.values
#输出auc
auc
0.9375
#画ROC曲线
roc <- performance(pred, "tpr", "fpr")
plot(roc, main = 'ROC Curve')

结果如下：

3) Tips：Confidence interval of sensitivity

根据混淆矩阵可以计算得到Sensitivity和Specificity：

$\text{Sensitivity/Recall/TPR} = \frac{TP}{TP + FN}$

$\text{Specificity} = \frac{TN}{FP + TN}$

Sensitivity是指是指实际为阳性的样本中，判断为阳性的比例，可以认为是伯努利实验成功的概率p的点估计。

因此Sensitivity的置信区间可等价于努利实验成功的概率p的置信区间估计。

3a) 大样本量计算 - 正态分布近似

当样本数较大(N>30)时，根据中心极限定律可以利用正太分布近似成功概率p。

$\hat p \pm 1.96 \sqrt{\frac{\hat p \left(1 - \hat p\right)}{n}}$

• 在线工具

通过下面链接完成计算， http://vassarstats.net/clin1.html。

• R实现

可以通过R软件binom包实现，需要已知TP，FN。

library(binom)
TP <-100
FN<-50
binom.confint(x = TP, n = FN+TP, method ="as")
'''
      method   x   n      mean     lower     upper
1 asymptotic 100 150 0.6666667 0.5912276 0.7421057
'''

3b) 大样本/中等样本量计算 - Bootstrap方法

Bootstrap估计是利用重复抽样的方法对参数进行估计的。我们可以通过R软件pROC包中Bootstrap估计进行Sensitivity的区间估计，这里为非参数估计。当样本量较大或者中等时，可采用这种方式。

可以利用R软件pROC包实现：

library(pROC)
data(aSAH)
r=plot.roc(aSAH$outcome, aSAH$s100b, ci=TRUE, of="thresholds")
ci(r)
'''
95% CI: 0.6301-0.8326 (DeLong)
'''
ci.se(r)
'''
95% CI (2000 stratified bootstrap replicates):
  sp se.low se.median se.high
 0.0 1.0000    1.0000  1.0000
 0.1 0.9067    0.9756  1.0000
 0.2 0.8195    0.9268  1.0000
 0.3 0.7555    0.8780  0.9683
 0.4 0.6780    0.8262  0.9356
 0.5 0.6000    0.7707  0.9024
 0.6 0.5366    0.6829  0.8488
 0.7 0.4878    0.6585  0.8000
 0.8 0.3510    0.5854  0.7805
 0.9 0.2195    0.3902  0.6098
 1.0 0.1463    0.2927  0.4390
 '''

3c) 小样本量计算 - 其他精确近似方法

当样本量较小时，我们可以威尔逊区间法或者Pearson-Klopper方法等。我们可以通过R的binom包实现，需要已知TP，FN。

library(binom)
TP <-100
FN<-50
binom.confint(x = TP, n = FN+TP, method ="
exact")
'''
  method   x   n      mean     lower     upper
1  exact 100 150 0.6666667 0.5851570 0.7414436
'''

binom.confint(x = TP, n = FN+TP, method ="
wilson")
'''
  method   x   n      mean     lower     upper
1 wilson 100 150 0.6666667 0.5878976 0.7371124
'''

• exact - Pearson-Klopper method. See also binom.test

• wilson - Wilson method.

其他的近似方法可以调节binom.confint中method参数。

4) Homework

用你熟悉的脚本语言对BreastCancer数据集的预测结果 BreastCancer_predict.txt 进行评估，画出ROC curve，提交代码和ROC curve的图片。BreastCancer_predict.txt 中第一列为真实的数据标签，包括良性(benign,0)和恶性(malignant,1)，第二列为预测为正类(malignant,1)的概率值，用制表符"\t"分隔。

对机器学习比较熟悉的同学可以选择自己训练二分类模型(使用任意分类器均可)并绘制ROC曲线。输入为 BreastCancer.csv 数据集，第一列为样本编号，最后一列是样本标签，剩下9列对应9个特征，用逗号","分隔。 如果自己拟合模型，需要进行数据空缺值填充（用特征的样本中位数补全缺失值），必要时做data scaling，并自行划分数据集。