Summary

Имеется набор данных “German Credit” — обезличенная информация о людях, бравших кредит в банке и вернувших, либо не вернувших его.

Необходимо было построить модель, прогнозирующую кредитный дефолт, и сделать краткое описание методики построения модели, самой модели, также дать численную оценку качества модели.

Было построено несколько моделей, предсказывающих дефолт (дерево решений, Байесовский классификатор с выбором параметров и выбором точки отсечения, Random Forest с выбором точки отсечения), произведена оценка их точности с учетом повышенной “стоимости” ошибок II рода (False Negative, предсказано отсутствие дефолта при том, что клиент совершил дефолт). Лучшей моделью оказался Баейсовский классификатор с правильным выбором параметров и оптимально подобранным порогом вероятности отнесения клиента к дефолтным.

Решение задачи

Первичный анализ данных

Загрузим наши данные и посмотрим на них.

setwd('~/R/DoubleData/')
german<-read.csv("germancredit.csv")
str(german)
## 'data.frame':    1000 obs. of  21 variables:
##  $ Default        : int  0 1 0 0 1 0 0 0 0 1 ...
##  $ checkingstatus1: Factor w/ 4 levels "A11","A12","A13",..: 1 2 4 1 1 4 4 2 4 2 ...
##  $ duration       : int  6 48 12 42 24 36 24 36 12 30 ...
##  $ history        : Factor w/ 5 levels "A30","A31","A32",..: 5 3 5 3 4 3 3 3 3 5 ...
##  $ purpose        : Factor w/ 10 levels "A40","A41","A410",..: 5 5 8 4 1 8 4 2 5 1 ...
##  $ amount         : int  1169 5951 2096 7882 4870 9055 2835 6948 3059 5234 ...
##  $ savings        : Factor w/ 5 levels "A61","A62","A63",..: 5 1 1 1 1 5 3 1 4 1 ...
##  $ employ         : Factor w/ 5 levels "A71","A72","A73",..: 5 3 4 4 3 3 5 3 4 1 ...
##  $ installment    : int  4 2 2 2 3 2 3 2 2 4 ...
##  $ status         : Factor w/ 4 levels "A91","A92","A93",..: 3 2 3 3 3 3 3 3 1 4 ...
##  $ others         : Factor w/ 3 levels "A101","A102",..: 1 1 1 3 1 1 1 1 1 1 ...
##  $ residence      : int  4 2 3 4 4 4 4 2 4 2 ...
##  $ property       : Factor w/ 4 levels "A121","A122",..: 1 1 1 2 4 4 2 3 1 3 ...
##  $ age            : int  67 22 49 45 53 35 53 35 61 28 ...
##  $ otherplans     : Factor w/ 3 levels "A141","A142",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ housing        : Factor w/ 3 levels "A151","A152",..: 2 2 2 3 3 3 2 1 2 2 ...
##  $ cards          : int  2 1 1 1 2 1 1 1 1 2 ...
##  $ job            : Factor w/ 4 levels "A171","A172",..: 3 3 2 3 3 2 3 4 2 4 ...
##  $ liable         : int  1 1 2 2 2 2 1 1 1 1 ...
##  $ tele           : Factor w/ 2 levels "A191","A192": 2 1 1 1 1 2 1 2 1 1 ...
##  $ foreign        : Factor w/ 2 levels "A201","A202": 1 1 1 1 1 1 1 1 1 1 ...
summary(german)
##     Default    checkingstatus1    duration    history      purpose   
##  Min.   :0.0   A11:274         Min.   : 4.0   A30: 40   A43    :280  
##  1st Qu.:0.0   A12:269         1st Qu.:12.0   A31: 49   A40    :234  
##  Median :0.0   A13: 63         Median :18.0   A32:530   A42    :181  
##  Mean   :0.3   A14:394         Mean   :20.9   A33: 88   A41    :103  
##  3rd Qu.:1.0                   3rd Qu.:24.0   A34:293   A49    : 97  
##  Max.   :1.0                   Max.   :72.0             A46    : 50  
##                                                         (Other): 55  
##      amount      savings   employ     installment    status     others   
##  Min.   :  250   A61:603   A71: 62   Min.   :1.000   A91: 50   A101:907  
##  1st Qu.: 1366   A62:103   A72:172   1st Qu.:2.000   A92:310   A102: 41  
##  Median : 2320   A63: 63   A73:339   Median :3.000   A93:548   A103: 52  
##  Mean   : 3271   A64: 48   A74:174   Mean   :2.973   A94: 92             
##  3rd Qu.: 3972   A65:183   A75:253   3rd Qu.:4.000                       
##  Max.   :18424                       Max.   :4.000                       
##                                                                          
##    residence     property        age        otherplans housing   
##  Min.   :1.000   A121:282   Min.   :19.00   A141:139   A151:179  
##  1st Qu.:2.000   A122:232   1st Qu.:27.00   A142: 47   A152:713  
##  Median :3.000   A123:332   Median :33.00   A143:814   A153:108  
##  Mean   :2.845   A124:154   Mean   :35.55                        
##  3rd Qu.:4.000              3rd Qu.:42.00                        
##  Max.   :4.000              Max.   :75.00                        
##                                                                  
##      cards         job          liable        tele     foreign   
##  Min.   :1.000   A171: 22   Min.   :1.000   A191:596   A201:963  
##  1st Qu.:1.000   A172:200   1st Qu.:1.000   A192:404   A202: 37  
##  Median :1.000   A173:630   Median :1.000                        
##  Mean   :1.407   A174:148   Mean   :1.155                        
##  3rd Qu.:2.000              3rd Qu.:1.000                        
##  Max.   :4.000              Max.   :2.000                        
##

Итак, у нас 1000 записей о клиентах, которые совершили либо не совершили кредитный дефолт (столбец Default). По каждому клиенту есть значения двадцати различных параметров, на основании которых мы можем пытаться строить модели. Незаполненных значений нет.

Параметр Default сейчас не факторный, исправим это.

german$Default=as.factor(german$Default)
summary(german$Default)
##   0   1 
## 700 300

Попробуем построить несколько моделей и оценить их с помощью перекрестной проверки, разбивая данные на 10 частей (10-fold cross-validation). Метрика задана в описании данных, с помощью cost-матрицы:

costmatrix=matrix(data = c(0,1,5,0), nrow = 2, dimnames = list(c(0,1), c(0,1)))
costmatrix
##   0 1
## 0 0 5
## 1 1 0

В отличие от описания, в нашей реализации ряды матрицы соответствуют предсказанным классам, а столбцы - настоящим значениям. Это cоответствует принятым конвенциям в R.

Заметим, что с такой cost-матрицей на наших тренировочных данных в качестве benchmark можно использовать модель, которая всегда будет предсказывать дефолт. Средняя “стоимость” такой модели будет равна \(\frac{700\cdot1 + 300\cdot0}{1000} = 0.7\).

Подготовим разбиение данных для кросс-валидации:

library(cvTools)
## Loading required package: lattice
## Loading required package: robustbase
folds = cvFolds(1000, K = 10)

Деревья принятия решений

Построим дерево с помощью алгоритма C5.0. В качестве входного параметра можно передавать cost-матрицу, дерево строится с ее учетом и не требует последующих доработок.

library(C50)
costs<-vector()
# Цикл кросс-валидации
for (i in 1:10){
  # строим модель по данным, тренировочным для данной итерации кросс-валидации
  model <- C5.0(Default ~ ., data = german[folds$which!=i,], costs = costmatrix)
  # предсказываем на оставшихся данных
  pred = predict(model, german[folds$which==i,])
  # запоминаем стоимость этой итерации
  costs[i]<-sum(table(pred, german[folds$which==i,]$Default)*costmatrix)/100
}
# оценка стоимости модели
mean(costs)
## [1] 0.603

Выбор параметров с помощью decision tree для Байесовского классификатора

В статье Scaling up the Naive Bayesian Classifier: Using Decision Trees for Feature Selection рассматривается алгоритм, который показывает на наших данных очень высокую точность классификации. Правда, он применялся в ситуации, когда в качестве метрики использовалась обычная misclassification error, но попробуем применить его и в нашей ситуации с соответствующими доработками. Смысл алгоритма в следующем: сначала с помощью деревьев принятия решений, построенных с помощью алгоритма C4.5 отбирались важные параметры, а затем строился Байесовский классификатор, для построения которого использовались только отобранные параметры.

Мы будем использовать доработанную версию алгоритма построения деревьев, C5.0.

Выбор параметров с помощью деревьев принятия решений

Для выбора параметров, с помощью которых мы будем делать предсказания, воспользуемся следующим алгоритмом:

  1. Случайно выберем 10% данных.
  2. Запустим по этим данным алгоритм C5.0 построения дерева решений. В качестве дополнительного параметра будем передавать алгоритму нашу cost-матрицу.
  3. Выберем те параметры, которые использовались на трех верхних уровнях дерева.
  4. Повторим шаги 1-3 пять раз.
  5. В качестве окончательного списка параметров, которые мы будем использовать в дальнейшем, возьмем объединение параметров из всех пяти итераций.

Ниже приведен код для одной итерации.

set.seed(1234)
hundred = sample (nrow(german), 100)
germantree = C5.0(Default ~ ., data = german[hundred,], costs = costmatrix)
plot(germantree)

Интересующие нас параметры: checkingstatus1, age.

После пяти итераций полный список выглядит так: savings, checkingstatus1, duration, others, liable, residence, age, history, employ, installment, otherplans.

Байесовский классификатор

Для того, чтобы учесть cost-матрицу, будем строить классификатор следующим образом:

  • На выходе модели будем получать не точную классификацию, а вероятности отнесения наблюдения в один из классов.
  • Для каждой итерации механизма кросс-валидации будем строить ROC-кривую, соответствующую нашей cost-матрице: по оси x будем откладывать точки отсечения, по оси y — стоимость модели при выборе соответствующей точки отсечения.
  • Сложим стоимости модели поточечно для всех десяти итераций кросс-валидации и усредним.
  • В получившемся графике найдем минимальную стоимость и соответствующую ей точку отсечения.
features<-c('Default','savings', 'checkingstatus1', 'duration','others', 'liable', 'residence', 'age', 'history', 'employ', 'installment', 'otherplans')
library(e1071)
library(ROCR)
## Loading required package: gplots
## 
## Attaching package: 'gplots'
## 
## Следующий объект скрыт от 'package:stats':
## 
##     lowess
#векторы для хранения точек отсечения и средней стоимости
cv_cutoff_y<-vector(length = 101)
cv_cutoff_x<-seq(from = 0, to = 1, by = .01)
for (i in 1:10){
  # строим модель по данным, тренировочным для данной итерации кросс-валидации
  model <- naiveBayes(Default ~., data = german[folds$which!=i,features])
  # предсказываем на оставшихся данных, тип "raw" дает нам вероятности вместо классификации
  pred = predict(model, german[folds$which==i,features], type = "raw")
  # упаковываем предсказания и истинные значения в объект prediction для работы с ROC-кривыми
  rocpred<-prediction(pred[,2], german[folds$which==i,]$Default)
  # генерим ROC-кривую с нашей cost-матрицей
  perf <- performance(rocpred,"cost", cost.fp = 1, cost.fn = 5)
  # линейно интерполируем ROC-кривую в равномерно взятых на отрезке [0, 1] точках, чтобы можно было осмысленно складывать ROC-кривые разных итераций.
  cv_cutoff_y<-cv_cutoff_y+approx(x = perf@x.values[[1]], y = perf@y.values[[1]], xout = cv_cutoff_x, yleft = Inf, yright = Inf)$y
  }
cv_cutoff_y<-cv_cutoff_y/10
#график "суммарной" ROC-кривой
plot(cv_cutoff_x, cv_cutoff_y, type = "l", xlab = "cutoff", ylab = "cost")

#минимальная стоимость модели
min(cv_cutoff_y)
## [1] 0.5398536
#значение порога для минимальной стоимости.
cv_cutoff_x[which.min(cv_cutoff_y)]
## [1] 0.14

Другие способы выбора переменных для Байесовского классификатора

В качестве альтернативного варианта попробуем выбирать переменные для Байесовского классификатора следующим образом:

  • Построим дерево решений C5.0 по всем данным
  • Отсортируем переменные в порядке их важности
  • Начнем с модели, содержащей только самую важную переменную, и будем добавлять следующие переменные по одной и строить классификатор для каждого поднабора из первых самых важных \(j\) переменных. Стоимость моделей для каждой такой итерации будем сохранять.
  • Выберем наиболее “дешевую” модель из построенных.
germantree = C5.0(Default ~ ., data = german, costs = costmatrix)
results<-data.frame(cutoff=vector(), cost = vector())
for (j in 1:length(rownames(C5imp(germantree)))){
  features<-c("Default", rownames(C5imp(germantree))[1:j])
  # векторы для хранения точек отсечения и средней стоимости
  cv_cutoff_y<-vector(length = 101)
  cv_cutoff_x<-seq(from = 0, to = 1, by = .01)
  for (i in 1:10){
    # строим модель по данным, тренировочным для данной итерации кросс-валидации
    model <- naiveBayes(Default ~., data = german[folds$which!=i,features])
    # предсказываем на оставшихся данных, тип "raw" дает нам вероятности вместо классификации
    pred = predict(model, german[folds$which==i,features], type = "raw")
    # упаковываем предсказания и истинные значения в объект prediction для работы с ROC-кривыми
    rocpred<-prediction(pred[,2], german[folds$which==i,]$Default)
    # генерим ROC-кривую с нашей cost-матрицей
    perf <- performance(rocpred,"cost", cost.fp = 1, cost.fn = 5)
    # линейно интерполируем ROC-кривую в равномерно взятых на отрезке [0, 1] точках, чтобы можно было осмысленно складывать ROC-кривые разных итераций.
    cv_cutoff_y<-cv_cutoff_y+approx(x = perf@x.values[[1]], y = perf@y.values[[1]], xout = cv_cutoff_x, yleft = Inf, yright = Inf)$y
  }
  cv_cutoff_y<-cv_cutoff_y/10
  # минимальная стоимость модели
  results<-rbind(results, c(cv_cutoff_x[which.min(cv_cutoff_y)], min(cv_cutoff_y)))
  # значение порога для минимальной стоимости.
}
colnames(results)<-c('cutoff','cost')
print(results)
##    cutoff      cost
## 1    0.23 0.5859169
## 2    0.21 0.5916352
## 3    0.20 0.5882811
## 4    0.16 0.5248223
## 5    0.18 0.5448823
## 6    0.09 0.5361921
## 7    0.16 0.5165513
## 8    0.13 0.5017926
## 9    0.11 0.5033886
## 10   0.10 0.5130845
## 11   0.12 0.5055973
## 12   0.10 0.5055682
## 13   0.13 0.5017194
## 14   0.11 0.5120122
## 15   0.14 0.5083248
## 16   0.16 0.5137367
## 17   0.14 0.5160599
## 18   0.12 0.5142644
## 19   0.12 0.5094448
## 20   0.13 0.5116233

Минимальная стоимость лучше, чем у предыдущих моделей и достигается при использовании следующих тринадцати переменных:

rownames(C5imp(germantree))[1:13]
##  [1] "checkingstatus1" "others"          "foreign"        
##  [4] "otherplans"      "purpose"         "history"        
##  [7] "employ"          "amount"          "job"            
## [10] "liable"          "duration"        "age"            
## [13] "savings"

Random forest

Построим модель с помощью RandomForest. Поскольку, в отличие от C5.0, алгоритму для построения random forest нельзя передать в качестве метрики cost-матрицу, построение и оценка модели делаются аналогично модели Байесовского классификатора:

library(randomForest)
## randomForest 4.6-10
## Type rfNews() to see new features/changes/bug fixes.
cv_cutoff_y<-vector(length = 101)
cv_cutoff_x<-seq(from = 0, to = 1, by = .01)
for (i in 1:10){
  model <- randomForest(Default ~ ., data = german[folds$which!=i,])
  pred = predict(model, german[folds$which==i,], type = "prob")
  rocpred<-prediction(pred[,2], german[folds$which==i,]$Default)
  perf <- performance(rocpred,"cost", cost.fp = 1, cost.fn = 5)
  # линейно интерполируем ROC-кривую в равномерно взятых на отрезке [0, 1] точках, чтобы можно было осмысленно складывать ROC-кривые разных итераций.
  cv_cutoff_y<-cv_cutoff_y+approx(x = perf@x.values[[1]], y = perf@y.values[[1]], xout = cv_cutoff_x, yleft = Inf, yright = Inf)$y
  }
cv_cutoff_y<-cv_cutoff_y/10
# график "суммарной" ROC-кривой
plot(cv_cutoff_x, cv_cutoff_y, type = "l", xlab = "cutoff", ylab = "cost")

# минимальная стоимость модели
min(cv_cutoff_y)
## [1] 0.5323428
# значение порога для минимальной стоимости.
cv_cutoff_x[which.min(cv_cutoff_y)]
## [1] 0.22

Мы видим, что результат оказался хуже, чем у Байесовского классификатора.

Заключение

Среди построенных моделей оптимальной оказался Байесовский классификатор, построенный на параметрах checkingstatus1, others, foreign, otherplans, purpose, history, employ, amount, job, liable, duration, age, savings и с точкой отсечения 0.13. Стоимость такой модели по оценке равна 0.5017194.