[R - ml] 模型的评估

2018-11-16 本文已影响1人 XuXiaolu

主要的模型评估方式有三种：

1. hold out

以构建模型的训练数据来评估模型，显然不是最佳选择，因为我们的目标是对未来数据的分析。
而模型是以训练集来建立的。
所以要 train, test

然而我们不能重复的进行测试以选择最佳模型，因此常用的方式是
train, validataion, test
用validatation 数据集来帮助寻找最佳的模型，然后用test 完成最后一次测试，作为对未来数据的验证

random_ids = order(runif(1000))
train = credit[random_ids[1:500], ]
validate = credit[random_ids[501:750], ]
test = credit[random_ids[751:1000], ]

如果数据的class 是非均衡的，就可能有问题

分层随机采样

require(caret)
in_train = createDataPartition(credit$default, p = 0.75, list = FALSE) # 不以list保存
train = credit[in_train, ]
test = credit[-in_train, ]

2. Cross-validatation : CV

10-fold cv

100 行数据
20 行用于测试
80 行用于训练模型
例如对于KNN 模型，我们想找到最佳k
对于每个k，我们进行10-fold CV。即把训练数据分10组，取其中的9组（72行）构建模型，
然后计算剩余的 8 行数据的预测错误，之后我们对10次计算后的错误平均。最终我们选择最小的那个模型。
0.2， 0.3， 0.4， 0.5 取最小的那个k
利用剩余80行数据，以及前面选择的 k 来构建模型
利用剩余的20行数据验证我们的模型

caret for CV

require(caret)
folds = createFolds(credit$Good.Loan, k = 10)
str(folds)
credit01_train = credit[ -folds$Fold01, ]
credit01_test = credit[folds$Fold01, ]

信用卡的例子：

library(caret)
library(C50)
library(irr)
set.seed(123)
folds = createFolds(credit$Good.Loan, k = 10)
cv_results = lapply(folds, function(x) {
  credit_train = credit[-x, ]
  credit_test = credit[x, ]
  credit_model = C5.0(Good.Loan ~. , data = credit_train)
  credit_pred = predict(credit_model, credit_test)
  credit_actual = credit_test$Good.Loan
  kappa = kappa2(data.frame(credit_actual, credit_pred))$value
  return(kappa)
})

str(cv_results)
mean(unlist(cv_results))

说明前面我们关于信用的坏账模型不怎么样。
k-fold 选多大是一个bias 与 vias 的问题。
k 越小显然 bias（偏差）更大， variance（方差）小，
而k越大， bias 小，但是 variance大。

3. boostrap

boostrap 源于统计中利用随机抽样来估计总体的方法。
类似机器学习中我们用多个随机选取的训练与测试及来估计模型的性能统计。
最终利用均值来作为未来新数据下的性能估计。
k-fold 与 bootstrap 的区别在于 bootstrap是有放回的，也就是说某些样本可能重复出现。

bootstrap 错误估计同时考虑过于乐观的训练集与过于悲观的测试集。
$error = 0.632 × error_{test} + 0.368 × error_{train}$

回归问题

$MeanSqureErroe MSE = \frac{1}{N} {\sum{(y_{est}-y_{actual})^2}}$
$RootMeanSquareError RMSE = \sqrt{MSE}$
$Relative Squre Erroe RSE = \frac{MSE} {MSE_{baseline}} = \frac{ \sum{(y_{est}-y_{actual})^2}}{\sum{(y_{mean}-y_{actual})^2}}$
R^2 模型的解释性：
$R^2 = \frac{MSE_{baseline}-MSE}{MSE_{baseline}}$
$R^2 = 1- RSE$

demo

require(car)
?Prestige
Prestige_clean = Prestige[!is.na(Prestige$type), ]
model = lm(prestige~. , data = Prestige_clean)
actual = Prestige_clean$prestige
rmse = (mean(score - actual)^2)^0.5
rmse
mu = mean(actual)
rse = sum(score - actual)^2 /sum(mu - actual)^2
rse
rsquare = 1 - rse
rsquare

$MeanSquareLogError MSLE = \frac{1}{N} {\sum(logy_est - logy_actual)^2 }$
$RootMeanSquareLogError MSLE = \sqrt MSLE$

分类问题

如何判断分类器性能？

实际分类值
预测分类
估计的预测分类概率

实际与预测对比

actual_outcom = test_data$out_come
predicted_outcome = predict(model, test_data)

预测分类的概率

R语言不同模型，参数不同。prob, posterior, raw, probability 都有可能

# naive bayes
predict(model, test_data, type = 'raw')
# decision tree
predict(model, test_data, type = 'prob')

Confusion Matrix

$accuracy = \frac{TN + TN} {TP + TN + FP + FN}$
$erroratte =\frac{FP + FN} {TP + TN + FP + FN} = 1 - accuracy$
table()

table(sms_raw_test$type, sms_test_pred)

require(gmodels)
CrossTable(sms_test_pred, sms_raw_test$type,
                     prop.chisq = FALSE, prop.t = FALSE,
                     dnn = c('predicted', 'actual'))

caret::ConfusionMatrix()

confusionMatrix(sms_test_pred, sms_raw_test$type, positive = 'spam')

次函数将输出Kappa
kappa statistic
.< 0.2 预测比差
.> 0.8 好
k = [Pr(a) - Pr(e)] / [1 - Pr(e)]

Kappa 统计是相对预期一致的准确性调整
Pr(e) 是预期值与实际值一致的概率
$Pr(e) = Pr(actual postive)Pr(predicted positive) + Pr(actual negative)Pr(predicted negative)$

require(vcd)
table(sms_raw_test$type)
table(sms_raw_test$type, sms_test_pred)
kappa(table(sms_raw_test$type, sms_test_pred))
pr.ap = 186 / (186 + 1206)
pr.pp = (157 + 4) / (4 + 157 + 1202 + 29)
pr.an = 1206 / (1206 + 186)
pr.pn = (1202 + 29) / (1202 + 29 + 4 + 157)
pr.e = pr.ap*pr.pp+pr.an*pr.pn
pr.a = (1202 + 157) / (1202 + 29 + 4 + 157)
kappa.c = (pr.a - pr.e) / (1 - pr.e)
kappa.c

sensitivity 和 specificity

sensitivity 度量 postive 样本被正确分类的比率
specificity 度量 negative 样本被正确分类的比率

# manual calculation
sensitivity(sms_test_pred, sms_raw_test$type, positive = 'spam')
specificity(sms_test_pred, sms_raw_test$type, negative = 'ham')

0~1 越接近 1 越好
同时考虑的方式 trade-off

precision 和 recall

$正确性: precision =\frac{TP}{TP + FP}$
$完备性: recall =\frac{TP}{TP + FN}$
实际上 recall = sensitivity
posPredValue 由caret提供函数计算

posPredValue(sms_test_pred, sms_raw_test$type, positive = 'spam')

F - measure

F - measure 是以调和均值将precision和recall 以一个指标来衡量。
$F - measure = \frac{ 2 * precision * recall} {precision + recall} = \frac{2TP} {2TP + FN + FP}$

ROC 曲线

receiver operating characteristic curve

AUC

Area Under Curve 被定义为ROC曲线下的面积，显然这个面积的数值不会大于1

require(ROCR)
sms_test_pred = predict(sms_classifier, sms_test, type = 'raw')
pred = prediction(sms_test_pred[, 2], sms_raw_test$type)
perf = performance(pred, measure = 'tpr', x.measure = 'fpr')
plot(perf, main = "ROC curve", col = 'blue', lwd = 2)
abline(a = 0, b = 1, lwd = 2, lty = 2) # 两步要一起操作

AUC 可以通过以下方式获得：

perf.auc = performance(pred, measure = 'auc')
# performance 是s4 的一个对象
str(perf.auc)
unlist(perf.auc@y.values)

library(ROCR)
library(e1071)
inTrain = createDataPartition(y = iris$Species, p = 0.75, list = FALSE)
iristrain = iris[inTrain, ]
iristest = iris[-inTrain, ]

nb_model = naiveBayes(Species~. , data = iristrain)
nb_prediction = predict(nb_model, iristest[, -5], type = 'raw')
score = nb_prediction[, c('virginica')]
actual_class = iristest$Species == 'virginica'
# Mix some noise to the score
# Make the score less precise for illustration
score = (score + runif(length(score))) / 2
pred = prediction(score, actual_class)
perf = performance(pred, 'prec', 'rec')
plot(perf)

library(ROCR)
library(e1071)
nb_model = naiveBayes(Species~. , data = iristrain)
nb_prediction = predict(nb_model, iristest[, -5], type = 'raw')
score = nb_prediction[, c('virginica')]
actual_class = iristest$Species == 'virginica'
# Mix some noise to the score
# Make the score less precise for illustration
score = (score + runif(length(score))) / 2
pred = prediction(score, actual_class)
perf = performance(pred, 'prec', 'rec')
auc = performance(pred, 'auc')
auc = unlist(slot(auc, 'y.values'))
plot(perf)
legend(0.6, 0.3, c(c(paste('AUC is', auc),'\n'), border = 'white', cex = 1.0,
                   box.col = 'white'))

ROC - cost

根据FP， FN 的cost来选择最佳的cutoff点
cost.fp, cost.fn 惩罚因子

# Plot the cost curve to find the best cutoff
# Assign the cost for False + ve and False - ve
perf = performance(pred, 'cost', cost.fp = 4, cost.fn = 1) 
plot(perf)