R Notebook
Математическое моделирование
Практика 3
Параметрические классификаторы для бинарной зависимой переменной (\(Y\))
В практических примерах ниже показано:
- как построить модель логистической регрессии
- как построить модель линейного дискриминантного анализа (LDA)
- как построить модель квадратичного дискриминантного анализа (QDA)
- как изменять границу отсечения вероятности для классов
- как построить ROC-кривую
Модели: логистическая регрессия, LDA, QDA.
Данные: сгенерированный набор по кредитным картам Default{ISLR}. В наборе 10000 наблюдений и 4 показателя:
* default – бинарная переменная: если Yes, держатель карты не вернул долг по кредитке;
* student – бинарная переменная: если Yes, держатель карты является студентом;
* balance – средний баланс, который оставался на балансе кредитной карты, после ежемесячного платежа;
* income – доход держателя карты.
Пакеты:
library('ISLR')
library('GGally')
library('MASS')Зададим ядро генератора случайных чисел и объём обучающей выборки.
my.seed <- 12345
train.percent <- 0.85
options("ggmatrix.progress.bar" = FALSE)Исходные данные: набор Default
# ?Default # справка по набору данных
head(Default)str(Default)## '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 = '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)
Обратите внимание: априорные вероятности классов, которые можно оценить как доли наблюдений каждого класса в выборке, неодинаковы. Неплательщиков (default == 1) намного меньше:
# доли наблюдений в столбце default
table(Default$default) / sum(table(Default$default))##
## No Yes
## 0.9667 0.0333Доля наименьшего класса, в данном случае 0.0333, это ошибка нулевого классификатора: если бы мы прогнозировали default = Yes для всех наблюдений, ровно в такой доле случаев мы бы ошиблись. Точность моделей целесообразно будет сравнивать с этой величиной.
Мы построим простые модели, используя всего один объясняющий фактор: balance.
Отбираем наблюдения в обучающую выборку
set.seed(my.seed)
inTrain <- sample(seq_along(Default$default),
nrow(Default)*train.percent)
df <- Default[inTrain, ]
# фактические значения на обучающей выборке
Факт <- df$defaultСтроим модели, чтобы спрогнозировать default
Логистическая регрессия
model.logit <- glm(default ~ balance, data = df, family = 'binomial')
summary(model.logit)##
## Call:
## glm(formula = default ~ balance, family = "binomial", data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2319 -0.1474 -0.0602 -0.0229 3.7485
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.056e+01 3.915e-01 -26.98 <2e-16 ***
## balance 5.423e-03 2.391e-04 22.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2428.1 on 8499 degrees of freedom
## Residual deviance: 1355.8 on 8498 degrees of freedom
## AIC: 1359.8
##
## Number of Fisher Scoring iterations: 8Параметры модели логистической регрессии значимы с вероятностью 0.99.
# прогноз: вероятности принадлежности классу 'Yes' (дефолт)
p.logit <- predict(model.logit, df, type = 'response')
Прогноз <- factor(ifelse(p.logit > 0.5, 2, 1),
levels = c(1, 2),
labels = c('No', 'Yes'))
# матрица неточностей
conf.m <- table(Факт, Прогноз)
conf.m## Прогноз
## Факт No Yes
## No 8195 30
## Yes 200 75Обратите внимание: вектор p.logit состоит из вероятностей принадлежности наблюдений к классам, а не из меток этих классов. Поэтому для прогноза нужно сделать разделение на два класса вручную, используя какую-то границу отсечения. В данном случае это 0.5 – значение по умолчанию. У этой модели низкая чувствительность:
# чувствительность
conf.m[2, 2] / sum(conf.m[2, ])## [1] 0.2727273# специфичность
conf.m[1, 1] / sum(conf.m[1, ])## [1] 0.9963526# верность
sum(diag(conf.m)) / sum(conf.m)## [1] 0.9729412LDA
model.lda <- lda(default ~ balance, data = Default[inTrain, ])
model.lda## Call:
## lda(default ~ balance, data = Default[inTrain, ])
##
## Prior probabilities of groups:
## No Yes
## 0.96764706 0.03235294
##
## Group means:
## balance
## No 806.2403
## Yes 1738.4972
##
## Coefficients of linear discriminants:
## LD1
## balance 0.002201984# прогноз: вероятности принадлежности классу 'Yes' (дефолт)
p.lda <- predict(model.lda, df, type = 'response')
Прогноз <- factor(ifelse(p.lda$posterior[, 'Yes'] > 0.5,
2, 1),
levels = c(1, 2),
labels = c('No', 'Yes'))
# матрица неточностей
conf.m <- table(Факт, Прогноз)
conf.m## Прогноз
## Факт No Yes
## No 8207 18
## Yes 222 53Отчёт по модели LDA содержит три раздела: априарные вероятности классов (Prior probabilities of groups), групповые средние объясняющих переменных (Group means) и коэффициенты линейной разделяющей границы (Coefficients of linear discriminants).
У этой модели чувствительность ещё меньше.
# чувствительность
conf.m[2, 2] / sum(conf.m[2, ])## [1] 0.1927273# специфичность
conf.m[1, 1] / sum(conf.m[1, ])## [1] 0.9978116# верность
sum(diag(conf.m)) / sum(conf.m)## [1] 0.9717647QDA
model.qda <- qda(default ~ balance, data = Default[inTrain, ])
model.qda## Call:
## qda(default ~ balance, data = Default[inTrain, ])
##
## Prior probabilities of groups:
## No Yes
## 0.96764706 0.03235294
##
## Group means:
## balance
## No 806.2403
## Yes 1738.4972# прогноз: вероятности принадлежности классу 'Yes' (дефолт)
p.qda <- predict(model.qda, df, type = 'response')
Прогноз <- factor(ifelse(p.qda$posterior[, 'Yes'] > 0.5,
2, 1),
levels = c(1, 2),
labels = c('No', 'Yes'))
# матрица неточностей
conf.m <- table(Факт, Прогноз)
conf.m## Прогноз
## Факт No Yes
## No 8203 22
## Yes 215 60Разделяющая граница квадратичного дисперсионного анализа нелинейна, поэтому коэффициентов в отчёте мы не видим. Чувсивительность чуть хуже, чем у логистичесокй регрессии.
# чувствительность
conf.m[2, 2] / sum(conf.m[2, ])## [1] 0.2181818# специфичность
conf.m[1, 1] / sum(conf.m[1, ])## [1] 0.9973252# верность
sum(diag(conf.m)) / sum(conf.m)## [1] 0.9721176Очевидно, такая ситуация с чувствительностью не может нас устраивать, поскольку высокое значение верности модели (accuracy) обусловлено исключительно большой долей одного из классов в выборке. В такой ситуации надо пожертвовать небольшой частью специфичности, чтобы подтянуть чувствительность. Сделаем это, изменив границу отсечения классов.
Подбор границы отсечения вероятностей классов
ROC-кривая для LDA
Для начала построим график совместного изменения чувствительности и специфичности с изменением вероятности отсечения от 0 до 1 – ROC-кривую. Для примера возьмём модель LDA.
# считаем 1-SPC и TPR для всех вариантов границы отсечения
x <- NULL # для (1 - SPC)
y <- NULL # для TPR
# заготовка под матрицу неточностей
tbl <- as.data.frame(matrix(rep(0, 4), 2, 2))
rownames(tbl) <- c('fact.No', 'fact.Yes')
colnames(tbl) <- c('predict.No', 'predict.Yes')
# вектор вероятностей для перебора
p.vector <- seq(0, 1, length = 501)
# цикл по вероятностям отсечения
for (p in p.vector){
# прогноз
Прогноз <- factor(ifelse(p.lda$posterior[, 'Yes'] > p,
2, 1),
levels = c(1, 2),
labels = c('No', 'Yes'))
# фрейм со сравнением факта и прогноза
df.compare <- data.frame(Факт = Факт, Прогноз = Прогноз)
# заполняем матрицу неточностей
tbl[1, 1] <- nrow(df.compare[df.compare$Факт == 'No' & df.compare$Прогноз == 'No', ])
tbl[2, 2] <- nrow(df.compare[df.compare$Факт == 'Yes' & df.compare$Прогноз == 'Yes', ])
tbl[1, 2] <- nrow(df.compare[df.compare$Факт == 'No' & df.compare$Прогноз == 'Yes', ])
tbl[2, 1] <- nrow(df.compare[df.compare$Факт == 'Yes' & df.compare$Прогноз == 'No', ])
# считаем характеристики
TPR <- tbl[2, 2] / sum(tbl[2, 2] + tbl[2, 1])
y <- c(y, TPR)
SPC <- tbl[1, 1] / sum(tbl[1, 1] + tbl[1, 2])
x <- c(x, 1 - SPC)
}
# строим ROC-кривую
par(mar = c(5, 5, 1, 1))
# кривая
plot(x, y, type = 'l', col = 'blue', lwd = 3,
xlab = '(1 - SPC)', ylab = 'TPR',
xlim = c(0, 1), ylim = c(0, 1))
# прямая случайного классификатора
abline(a = 0, b = 1, lty = 3, lwd = 2)
# точка для вероятности 0.5
points(x[p.vector == 0.5], y[p.vector == 0.5], pch = 16)
text(x[p.vector == 0.5], y[p.vector == 0.5], 'p = 0.5', pos = 4)
# точка для вероятности 0.2
points(x[p.vector == 0.2], y[p.vector == 0.2], pch = 16)
text(x[p.vector == 0.2], y[p.vector == 0.2], 'p = 0.2', pos = 4)
Видно, что изменение границы отсечения с 0.5 до 0.2 увеличивает чувствительность модели почти в три раза, в то время как специфичность ухудшается незначительно. Матрица неточностей и её характеристики для LDA с p = 0.2:
Прогноз <- factor(ifelse(p.lda$posterior[, 'Yes'] > 0.2,
2, 1),
levels = c(1, 2),
labels = c('No', 'Yes'))
conf.m <- table(Факт, Прогноз)
conf.m## Прогноз
## Факт No Yes
## No 8028 197
## Yes 120 155# чувствительность
conf.m[2, 2] / sum(conf.m[2, ])## [1] 0.5636364# специфичность
conf.m[1, 1] / sum(conf.m[1, ])## [1] 0.9760486# верность
sum(diag(conf.m)) / sum(conf.m)## [1] 0.9627059Упражнение 3
На наборе данных из своего варианта построить указанные модели для прогноза бинарной зависимой переменной. Доля обучающей выборки – 75%.
Построить три графика:
1. Матричный график взаимного разброса переменных модели (ggpairs).
1. Две ROC-кривые на одних осях: сравнение качества прогноза сравниваемых моделей на обучающей выборке.
1. Две ROC-кривые на одних осях: сравнение качества прогноза сравниваемых моделей на тестовой выборке.
В конце файла с кодом в комментарии сравнить модели по качеству с помощью ROC-кривых. Сделать предположение о том, что в данном случае повлияло на преимущество одного метода над другим.
Скрипт решения и графики разместить в репозитории на github, прислать ссылку на почту преподавателя.
Номер варианта – номер студента в списке. Студент под номером 19 берёт вариант 1, под номером 20 – 2, и т.д.
Варианты
| Вариант |
Ядро для set.seed()
|
Данные | Зависимая переменная | Объясняющие переменные | Методы |
| 1 | 123 | Default{ISLR} – долги по кредитным картам | default (Yes – наличие признака, No – отсутствие) | все остальные | Логистическая регрессия, LDA |
| 2 | 123 | titanic_train{titanic} – выжившие в катастрофе Титаника | survived | все остальные, кроме name | Логистическая регрессия, LDA |
| 3 | 123 | PimaIndiansDiabetes{mlbench} – случаи диабета у женщин индейского племени Пима | diabetes (pos – наличие признака, neg – отсутствие) | все остальные | Логистическая регрессия, LDA |
| 4 | 123 | Glass{mlbench} – химический состав разных типов стекла | Type 1(1 – наличие признака, все остальные – отсутствие) | все остальные | Логистическая регрессия, LDA |
| 5 | 234 | Default{ISLR} – долги по кредитным картам | default (Yes – наличие признака, No – отсутствие) | все остальные | Логистическая регрессия, QDA |
| 6 | 234 | titanic_train{titanic} – выжившие в катастрофе Титаника | survived | все остальные, кроме name | Логистическая регрессия, QDA |
| 7 | 234 | PimaIndiansDiabetes{mlbench} – случаи диабета у женщин индейского племени Пима | diabetes (pos – наличие признака, neg – отсутствие) | все остальные | Логистическая регрессия, QDA |
| 8 | 234 | Glass{mlbench} – химический состав разных типов стекла | Type 1(1 – наличие признака, все остальные – отсутствие) | все остальные | Логистическая регрессия, QDA |
| 9 | 345 | Default{ISLR} – долги по кредитным картам | default (Yes – наличие признака, No – отсутствие) | все остальные | LDA, QDA |
| 10 | 345 | titanic_train{titanic} – выжившие в катастрофе Титаника | survived | все остальные, кроме name | LDA, QDA |
| 11 | 345 | PimaIndiansDiabetes{mlbench} – случаи диабета у женщин индейского племени Пима | diabetes (pos – наличие признака, neg – отсутствие) | все остальные | LDA, QDA |
| 12 | 345 | Glass{mlbench} – химический состав разных типов стекла | Type 1(1 – наличие признака, все остальные – отсутствие) | все остальные | LDA, QDA |
| 13 | 678 | Glass{mlbench} – химический состав разных типов стекла | Type 2(1 – наличие признака, все остальные – отсутствие) | все остальные | Логистическая регрессия, LDA |
| 14 | 678 | Glass{mlbench} – химический состав разных типов стекла | Type 2(1 – наличие признака, все остальные – отсутствие) | все остальные | Логистическая регрессия, QDA |
| 15 | 678 | Glass{mlbench} – химический состав разных типов стекла | Type 2(1 – наличие признака, все остальные – отсутствие) | все остальные | LDA, QDA |
| 16 | 678 | Glass{mlbench} – химический состав разных типов стекла | Type 7(1 – наличие признака, все остальные – отсутствие) | все остальные | Логистическая регрессия, LDA |
| 17 | 678 | Glass{mlbench} – химический состав разных типов стекла | Type 7(1 – наличие признака, все остальные – отсутствие) | все остальные | Логистическая регрессия, QDA |
| 18 | 678 | Glass{mlbench} – химический состав разных типов стекла | Type 7(1 – наличие признака, все остальные – отсутствие) | все остальные | LDA, QDA |
Как загрузить данные Default
# install.packages('ISLR')
library('ISLR')
head(Default)dim(Default)## [1] 10000 4Как загрузить данные titanic_train
# install.packages('titanic')
library('titanic')
head(titanic_train)dim(titanic_train)## [1] 891 12Как загрузить данные PimaIndiansDiabetes
# install.packages('mlbench')
library('mlbench')
data(PimaIndiansDiabetes)
head(PimaIndiansDiabetes)dim(PimaIndiansDiabetes)## [1] 768 9Как загрузить данные Glass
# install.packages('mlbench')
library('mlbench')
data(Glass)
head(Glass)dim(Glass)## [1] 214 10Источники
- James G., Witten D., Hastie T. and Tibshirani R. An Introduction to Statistical Learning with Applications in R. URL: http://www-bcf.usc.edu/~gareth/ISL/ISLR%20First%20Printing.pdf