【R语言 数据建模】模型验证武器

分类模型是数据挖掘中应用非常广泛的算法之一，常用的分类算法有Logistic模型、决策树、随机森林、神经网络、Boosting等。针对同一个数据集，可以有这么多的算法进行分析，那如何评估什么样的模型比较合理呢？本文就讲讲常用的模型验证武器，主要包括混淆矩阵、ROC曲线、提升度、增益法和KS统计量。

一、混淆矩阵

混淆矩阵就是如下图所示的那样，也是最简单的一种模型验证方法：

通过混淆矩阵可以算出模型预测精度((a+d)/(a+b+c+d))、正例覆盖率(d/(c+d))、负例覆盖率(a/(a+b))等。通过这么些指标综合考虑模型的预测准确率。

二、ROC曲线

在讲解ROC曲线之前，我们先看看几个定义：

Sensitivity:正确预测到的正例数/实际正例总数,即d/(c+d)
Specificity:正确预测到的负例数/实际负例总数,即a/(a+b)

ROC曲线就是根据这两个指标值绘制出来的，其中x轴为1-Specificity，y轴为Sensitivity。
通过比较ROC曲线与45°直线可以直观的反映模型的好坏，但并不能从定量的角度反馈模型好是好到什么程度或模型差是差到什么程度。那么就引申出了AUC的概念，即ROC曲线下的面积。当曲线偏离45°直线越远，则AUC越大，模型相应就会越好。一般认为AUC在0.75以上，模型就可以接受了。

三、提升度Lift

在讲解提升度曲线之前，我们先看看几个定义：

Pi：测试集中正例的比例，即(c+d)/(a+b+c+d)
Ptp：正确预测到的正例个数占总观测值的比例，即d/a+b+c+d=Pi1* Sensitivity
Pfp：把负例错误地预测成正例的个数占总数的比例，即b/a+b+c+d=(1-Pi1)*(1- Specificity)
Depth：预测成正例的比例，即b+d/a+b+c+d=Ptp+Pfp
PV_Plus:正确预测到的正例数/预测正例总数，即d/(b+d)=Ptp/depth
提升度Lift=(d/b+d)/(c+d/a+b+c+d)=PV_plus/Pi1

Lift曲线就是根据Depth和Lift两个指标绘制而成，它反映了预测正例的正真准确率。

四、增益法Gain

其实增益法Gain与提升度是一个事物的两种说法，从公式中就可以看出：
Gain=d/(b+d)=PV_plus
Gain与提升度相比并没有除以Pi值。

K-S统计量
统计学中，对于单样本的K-S检验就是利用样本数据来推断其是否服从某种分布，对于两样本的K-S检验主要推测的是两个样本是否具有相同的分布，对于模型的评估，希望正例的累积概率分布与负例的累积概率分布存在显著差异。

所以我们使用K-S统计量刻画模型的优劣，即使正例与负例的累积概率差达到最大。这是一个定量的判断规则，如下图所示，为传统的评价准则

ks值         预测表现
<0.2           差
0.2~0.4        一般
0.4~0.5        好
0.5~0.6        很好
0.6~0.75       非常好
0.75~1         完美

通常要求模型KS值在0.4以上。

废话不多说，下面我们看看如何使用R语言实现这些评估模型的方法。

R语言源代码学习：

#读取数据
dmagecr <- read.table(file = file.choose(), head = TRUE, sep = '')

#数据结构
str(dmagecr)

#指定感兴趣的水平为bad
dmagecr$good_bad <- factor(dmagecr$good_bad, levels = c('good','bad'),ordered = TRUE)

#创建训练集和测试集
set.seed(1234)
index <- sample(c(1,2), size = nrow(dmagecr), replace = TRUE, prob = c(0.7,0.3))
train <- dmagecr[index == 1,]
test <- dmagecr[index == 2,] 

#构建Logistic模型
model <- glm(formula = good_bad ~ checking+history+duration+savings+property, family = binomial(link = "logit"), data = train)

#模型结果查看
summary(model)

#模型的显著性检验
anova(object = model, test = 'Chisq')


#模型的拟合优度检验--HL检验
library(sjmisc)
HL_test <- hoslem_gof(x = model)
HL_test


#H-L的P值显著大于0.05，即接受实际值与预测值相吻合的原假设，再次说明模型是比较理想的
#模型预测
probility <- predict(object = model, newdata = test[,-21], type = 'response')
predict <- ifelse(probility > 0.5, 'bad', 'good')
#转型为因子
predict <- factor(predict, levels = c('good','bad'), order = TRUE)
#模型评估混淆矩阵
Freq <- table(test[,21], predict)
#预测精度
Accuracy <- sum(diag(Freq))/sum(Freq)
Freq;Accuracy


#还可以使用ROC曲线下的面积AUC、提升度Lift、增益法Gain和K-S统计量。下面就深入介绍这几种方法：

#ROC曲线
library(pROC)
roc_curve <- roc(test[,21],probility)
names(roc_curve)
Specificity <- roc_curve$specificities
Sensitivity <- roc_curve$sensitivities
library(ggplot2)
p <- ggplot(data = NULL, mapping = aes(x= 1-Specificity, y = Sensitivity))
p + geom_line(colour = 'red') +geom_abline(intercept = 0, slope = 1)+ annotate('text', x = 0.4, y = 0.5, label=paste('AUC=',round(roc_curve$auc,2)))+ labs(x = '1-Specificity',y = 'Sensitivity', title = 'ROC Curve')



#Lift曲线
Pi <- table(test$good_bad)[2]/sum(table(test$good_bad))
Ptp <- Pi*Sensitivity
Pfp <- (1-Pi)*(1-Specificity)
Depth <- Ptp + Pfp
PV_Plus <- Ptp/Depth
Lift <- PV_Plus/Pi

p <- ggplot(data = NULL, mapping = aes(x= Depth, y = Lift))
p + geom_line(colour = 'blue') + labs(x = 'Depth',y = 'Lift', title = 'Lift Curve')


#Gain曲线
p <- ggplot(data = NULL, mapping = aes(x= Depth, y = PV_Plus))
p + geom_line(colour = 'blue')  + labs(x = 'Depth',y = 'PV_Plus', title = 'Gain Curve')

#实际上，Lift曲线与Gain曲线长的一模一样，只不过是纵坐标不同而已。


#准备K-S数据
ks_data <- as.data.frame(cbind(good_bad=test[,21], probility))
good_ks <- ks_data[which(ks_data$good_bad==1),'probility']
bad_ks <- ks_data[which(ks_data$good_bad==2),'probility']

#自定义计算累计分布函数值
KS_Data <- function(x, y){
  gaps_x <- seq(min(x), max(x), length=1000)
  cauculate_x <- numeric()
  for(i in 1:length(gaps_x)){
    cauculate_x[i] <- sum(x<=gaps_x[i])/length(x)
  }
  gaps_x <- sort((gaps_x-min(gaps_x))/(max(gaps_x)-min(gaps_x)))

  gaps_y <- seq(min(y), max(y), length=1000)
  cauculate_y <- numeric()
  for(i in 1:length(gaps_y)){
    cauculate_y[i] <- sum(y<=gaps_y[i])/length(y)
  }
  gaps_y <- sort((gaps_y-min(gaps_y))/(max(gaps_y)-min(gaps_y)))
  return(list(df = data.frame(rbind(data.frame(Gaps = gaps_x,Cauculate = cauculate_x,Type = 'Positive'),data.frame(Gaps = gaps_y,Cauculate = cauculate_y,Type = 'Negtive'))), KS = max(abs(cauculate_y-cauculate_x)), x = gaps_y[which.max(abs(cauculate_y-cauculate_x))],y = abs(cauculate_x[which.max(abs(cauculate_y-cauculate_x))]-cauculate_y[which.max(abs(cauculate_y+cauculate_x))])/2))
}

#绘制K-S曲线
ggplot(data = KS_Data(bad_ks,good_ks)$df, mapping = aes(x = Gaps, y = Cauculate, colour = Type)) + geom_line() + theme(legend.position='none') + annotate(geom = 'text', x = KS_Data(bad_ks,good_ks)$x, y = KS_Data(bad_ks,good_ks)$y, label = paste('K-S Value: ', round(KS_Data(bad_ks,good_ks)$KS,2))) + labs(x = 'Probility', y = 'CDF')