Практика 3. Параметрические классификаторы для бинарной зависимой переменной (Y)
Задание
На наборе данных из своего варианта построить указанные модели для прогноза бинарной зависимой переменной. Доля обучающей выборки – 75%.
Построить три графика:
1. Матричный график взаимного разброса переменных модели (ggpairs).
1. Две ROC-кривые на одних осях: сравнение качества прогноза сравниваемых моделей на обучающей выборке.
1. Две ROC-кривые на одних осях: сравнение качества прогноза сравниваемых моделей на тестовой выборке.
В конце файла с кодом в комментарии сравнить модели по качеству с помощью ROC-кривых. Сделать предположение о том, что в данном случае повлияло на преимущество одного метода над другим.
Скрипт решения и графики разместить в репозитории на github, прислать ссылку на почту преподавателя.
Номер варианта – номер студента в списке. Студент под номером 19 берёт вариант 1, под номером 20 – 2, и т.д.
Вариант №5
| Вариант | Ядро для set.seed() |
Данные | Зависимая переменная | Объясняющие переменные | Методы |
|---|---|---|---|---|---|
| 5 | 234 | Default{ISLR} – долги по кредитным картам | default (Yes – наличие признака, No – отсутствие) |
все остальные | Логистическая регрессия, QDA |
Подключим пакеты функций:
library('ISLR')
library('GGally')## Loading required package: ggplot2## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library('MASS')
library('ROCR') # ROC-кривыеЗададим ядро генератора случайных чисел (по номеру варианта) и объём обучающей выборки.
my.seed <- 234
train.percent <- 0.75
options("ggmatrix.progress.bar" = FALSE)Подключим исходные данные и исследуем их структуру
## 'data.frame': 10000 obs. of 4 variables:
## $ default: Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ student: Factor w/ 2 levels "No","Yes": 1 2 1 1 1 2 1 2 1 1 ...
## $ balance: num 730 817 1074 529 786 ...
## $ income : num 44362 12106 31767 35704 38463 ...## default student balance income
## 1 No No 729.5265 44361.625
## 2 No Yes 817.1804 12106.135
## 3 No No 1073.5492 31767.139
## 4 No No 529.2506 35704.494
## 5 No No 785.6559 38463.496
## 6 No Yes 919.5885 7491.559Построим матричный график взаимного разброса переменных модели На графиках класс неплательщиков (default = ‘Yes’) показан красным.
# графики разброса
ggp <- ggpairs(Default, aes(colour = default, alpha = 0.4)) +
scale_colour_manual(values = c('darkgreen', 'red')) +
scale_fill_manual(values = c('darkgreen', 'red'))
print(ggp, progress = FALSE)## Warning in ggmatrix_gtable(x, ...): Please use the 'progress' parameter in
## your ggmatrix-like function call. See ?ggmatrix_progress for a few examples.
## ggmatrix_gtable 'progress' and 'progress_format' will soon be deprecated.TRUE## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Априорные вероятности классов, которые можно оценить как доли наблюдений каждого класса в выборке, неодинаковы. Неплательщиков (default == 1) намного меньше:
# доли наблюдений в столбце default
table(Default$default) / sum(table(Default$default))##
## No Yes
## 0.9667 0.0333Доля наименьшего класса, в данном случае 0.0333, это ошибка нулевого классификатора: если бы мы прогнозировали default = Yes для всех наблюдений, ровно в такой доле случаев мы бы ошиблись. Точность моделей целесообразно будет сравнивать с этой величиной. Построим модели, используя в качестве объясняющего фактора все остальные переменные.
Делим наблюдения на обучающую и тестовую выборку
set.seed(my.seed)
inTrain <- sample(seq_along(Default$default), nrow(Default)*train.percent)
df_train <- Default[inTrain, ]
df_test <- Default[-inTrain, ]
# фактические значения на обучающей выборке
Факт <- df_train$defaultСтроим модели для прогнозирования default
1.Логистическая регрессия
model.logit <- glm(default ~ ., data = df_train, family = 'binomial')
summary(model.logit)##
## Call:
## glm(formula = default ~ ., family = "binomial", data = df_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1867 -0.1277 -0.0489 -0.0173 3.8071
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.159e+01 6.089e-01 -19.035 <2e-16 ***
## studentYes -4.828e-01 2.842e-01 -1.699 0.0894 .
## balance 5.960e-03 2.838e-04 21.006 <2e-16 ***
## income 9.880e-06 9.772e-06 1.011 0.3120
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2103.9 on 7499 degrees of freedom
## Residual deviance: 1092.0 on 7496 degrees of freedom
## AIC: 1100
##
## Number of Fisher Scoring iterations: 8Параметр balance модели логистической регрессии значим с вероятностью 0.99.
# прогноз: вероятности принадлежности классу 'Yes' (дефолт)
forecast_logit_train <- predict(model.logit, df_train, type = 'response')Вектор forecast_logit_train состоит из вероятностей принадлежности наблюдений к классам, а не из меток этих классов. Поэтому для прогноза нужно сделать разделение на два класса вручную, используя какую-то границу отсечения. В данном случае это 0.5 – значение по умолчанию.
Прогноз <- factor(ifelse(forecast_logit_train > 0.5, 2, 1),levels = c(1, 2),labels = c('No', 'Yes'))
# матрица неточностей
conf.m <- table(Факт, Прогноз)
conf.m## Прогноз
## Факт No Yes
## No 7237 26
## Yes 148 89У этой модели низкая чувствительность:
# чувствительность
conf.m[2, 2] / sum(conf.m[2, ])## [1] 0.3755274# специфичность
conf.m[1, 1] / sum(conf.m[1, ])## [1] 0.9964202# верность
sum(diag(conf.m)) / sum(conf.m)## [1] 0.97682.QDA
Построим модель:
model.qda <- qda(default ~ ., data = Default[inTrain, ])
model.qda## Call:
## qda(default ~ ., data = Default[inTrain, ])
##
## Prior probabilities of groups:
## No Yes
## 0.9684 0.0316
##
## Group means:
## studentYes balance income
## No 0.2899628 802.8205 33627.38
## Yes 0.3628692 1766.1912 32944.40# прогноз: вероятности принадлежности классу 'Yes' (дефолт)
forecast_qda_train <- predict(model.qda, df_train, type = 'response')
Прогноз <- factor(ifelse(forecast_qda_train$posterior[, 'Yes'] > 0.5, 2, 1),levels = c(1, 2), labels = c('No', 'Yes'))
# матрица неточностей
conf.m <- table(Факт, Прогноз)
conf.m## Прогноз
## Факт No Yes
## No 7245 18
## Yes 161 76Разделяющая граница квадратичного дисперсионного анализа нелинейна, поэтому коэффициентов в отчёте мы не видим. Чувствительность немного выше, чем у логистичесокй регрессии.
# чувствительность
conf.m[2, 2] / sum(conf.m[2, ])## [1] 0.3206751# специфичность
conf.m[1, 1] / sum(conf.m[1, ])## [1] 0.9975217# верность
sum(diag(conf.m)) / sum(conf.m)## [1] 0.9761333Такая чувствительность не может нас устраивать, поскольку высокое значение верности модели (accuracy) обусловлено исключительно большой долей одного из классов в выборке. В такой ситуации надо пожертвовать небольшой частью специфичности, чтобы подтянуть чувствительность. Сделаем это, изменив границу отсечения классов.
ROC-кривые
Построим следующие ROC-кривые:
Две ROC-кривые на одних осях: сравнение качества прогноза сравниваемых моделей на обучающей выборке.
Две ROC-кривые на одних осях: сравнение качества прогноза сравниваемых моделей на тестовой выборке.
# функция построения ROC-кривой: pred -- прогноз, truth -- факт
rocplot <- function(pred, truth, ...){
require(ROCR)
# вызов функции непосредственно из пакета ROCR
predob = ROCR::prediction(pred, truth)
perf = performance(predob, "tpr", "fpr")
plot(perf,...)}
# Прогноз для обучающей выборки
p_logit_train <- predict(model.logit, df_train, type = 'response')
p_qda_train <- predict(model.qda, df_train, type = 'response')
forecast_logit_train <- factor(ifelse(forecast_logit_train > 0.5, 2, 1),levels = c(1, 2), labels = c('No', 'Yes'))
forecast_qda_train <- factor(ifelse(forecast_qda_train$posterior[, 'Yes'] > 0.5, 2, 1),levels = c(1, 2), labels = c('No', 'Yes'))
# Прогноз для тестовой выборки
p_logit_test <- predict(model.logit, df_test, type = 'response')
p_qda_test <- predict(model.qda, df_test, type = 'response')
forecast_logit_test <- factor(ifelse(p_logit_test > 0.5, 2, 1),levels = c(1, 2), labels = c('No', 'Yes'))
forecast_qda_test <- factor(ifelse(p_qda_test$posterior[, 'Yes'] > 0.5, 2, 1),levels = c(1, 2), labels = c('No', 'Yes'))
par(mfrow = c(1, 2))
# график для обучающей выборки
rocplot(p_logit_train, df_train[, "default"], main = "Обучающий набор")
rocplot(p_qda_train[["posterior"]][,1], df_train[, "default"], add = T, col = "red")
# график для тестовой выборки
rocplot(p_logit_test, df_test[, "default"], main = "Тестовый набор")
rocplot(p_qda_test[["posterior"]][,1], df_test[, "default"], add = T, col = "red")