Семинар 8. Кросс-валидация
Ноябрь, 12, 2018
Шаманское заклинание для настройки глобальных опций отчёта:
library(tidyverse) # обработка данных, графики...
library(skimr) # описательные статистики
library(rio) # импорт фантастического количества форматов данных
library(ISLR) # ещё данные
library(caret) # пакет для подбора параметров разных моделей
library(FFTrees) # быстрые деревья
library(margins) # для подсчёта предельных эффектов
library(rpart.plot) # для картинок деревьев
library(plotROC) # визуализация ROC-кривой
library(ggeffects) # графики для предельных эффектов
library(MLmetrics) # метрики качества
library(ranger) # случайный лес
library(factoextra) # графики для кластеризации и pca
library(elasticnet) # LASSO
library(latex2exp) # формулы в подписях к графику
library(distances) # расчет различных расстоянийСнова импортируем набор данных по успеваемости студентов и смотрим на него :)
educ <- import('data/xAPI-Edu-Data.csv')
skim(educ)## Skim summary statistics
## n obs: 480
## n variables: 17
##
## Variable type: character
## variable missing complete n min max empty n_unique
## Class 0 480 480 1 1 0 3
## GradeID 0 480 480 4 4 0 10
## NationalITy 0 480 480 2 11 0 14
## ParentAnsweringSurvey 0 480 480 2 3 0 2
## ParentschoolSatisfaction 0 480 480 3 4 0 2
## PlaceofBirth 0 480 480 3 11 0 14
## Relation 0 480 480 3 6 0 2
## SectionID 0 480 480 1 1 0 3
## Semester 0 480 480 1 1 0 2
## StageID 0 480 480 10 12 0 3
## StudentAbsenceDays 0 480 480 7 7 0 2
## Topic 0 480 480 2 9 0 12
## gender 0 480 480 1 1 0 2
##
## Variable type: integer
## variable missing complete n mean sd p0 p25 p50 p75 p100
## AnnouncementsView 0 480 480 37.92 26.61 0 14 33 58 98
## Discussion 0 480 480 43.28 27.64 1 20 39 70 99
## VisITedResources 0 480 480 54.8 33.08 0 20 65 84 99
## raisedhands 0 480 480 46.77 30.78 0 15.75 50 75 100
## hist
## <U+2587><U+2586><U+2585><U+2585><U+2585><U+2583><U+2583><U+2582>
## <U+2587><U+2587><U+2587><U+2587><U+2583><U+2585><U+2585><U+2583>
## <U+2586><U+2583><U+2582><U+2581><U+2583><U+2583><U+2587><U+2587>
## <U+2587><U+2586><U+2583><U+2583><U+2583><U+2586><U+2586><U+2583>Данные взяты с kaggle.com. Целевая переменная — успеваемость студента, Class, принимает три значения, ‘H’ (high), ‘M’ (mid), ‘L’ (low). Остальные переменные — характеристики студента.
Повторяем шаги предыдущего семинара:
- Для целей бинарной классификации объединяем две верхних категории в одну
- Уберём старую переменную ‘Class’ с тремя категориями
- Объявляем все текстовые переменные факторными:
educ_logit <- mutate(educ, y = fct_collapse(Class, H = c('M', 'H'))) %>%
select(-Class)
educ_fct <- mutate_if(educ_logit, is.character, factor)По умолчанию, базовой категорией факторной переменной будет категория, название которой стоит раньше в алфавите. В нашем примере ‘H’ идёт в алфавите раньше ‘L’, поэтому при построении логистической регрессии категория ‘H’ будет соответствовать \(y_i=0\) в формулах, а категория ‘L’ — \(y_i=1\). Следовательно, положительный коэффициент при регрессоре будет означать его положительное влияение на вероятность ‘L’.
Конечно, можно не полагаться на алфавит и самому указать любую базовую категорию.
Создадим другую набор данных с альтернативной базовой категорией:
educ_fct_rel <- mutate(educ_fct, y = fct_relevel(y, 'L'))- Разбиваем выборку на две части: 80% в обучающей и 20% в тестовой
set.seed(777)
train_rows <- createDataPartition(educ_fct$y, p = 0.8, list = FALSE)
educ_train <- educ_fct[train_rows, ]
educ_test <- educ_fct[-train_rows, ]- Упражнение 1.
Разделите набор данных Default из пакета ISLR на обучающую и тестовую выборки так, чтобы в первую попало 70% всех наблюдений.
def <- Default
glimpse(def)## Observations: 10,000
## Variables: 4
## $ default <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, No...
## $ student <fct> No, Yes, No, No, No, Yes, No, Yes, No, No, Yes, Yes, N...
## $ balance <dbl> 729.5265, 817.1804, 1073.5492, 529.2506, 785.6559, 919...
## $ income <dbl> 44361.625, 12106.135, 31767.139, 35704.494, 38463.496,...set.seed(777)
train_rows <- createDataPartition(def$default, p = 0.7, list = FALSE)
def_train <- def[train_rows, ]
def_test <- def[-train_rows, ]1 Логистическая регрессия
Оцениваем логистическую регрессию с помощью пакета caret.
educ_lmodel <- train(data = educ_train,
y ~ gender + SectionID + raisedhands, family = binomial(link = 'logit'),
method = 'glm')
summary(educ_lmodel)##
## Call:
## NULL
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8009 -0.4881 -0.1803 0.6427 3.2122
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.745737 0.365406 2.041 0.04127 *
## genderM 1.004068 0.335521 2.993 0.00277 **
## SectionIDB -0.348347 0.319917 -1.089 0.27621
## SectionIDC -0.702561 0.601374 -1.168 0.24270
## raisedhands -0.069383 0.008099 -8.567 < 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: 445.18 on 384 degrees of freedom
## Residual deviance: 271.91 on 380 degrees of freedom
## AIC: 281.91
##
## Number of Fisher Scoring iterations: 6- Упражнение 2.
Оцените заново логистическую регессию по набору данных def с помощью функций пакета caret.
def_lmodel <- train(data = def_train, default ~ ., family = binomial(link = 'logit'), method = 'glm')
summary(def_lmodel)##
## Call:
## NULL
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5134 -0.1415 -0.0537 -0.0193 3.6837
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.101e+01 5.951e-01 -18.498 <2e-16 ***
## studentYes -6.193e-01 2.847e-01 -2.176 0.0296 *
## balance 5.786e-03 2.810e-04 20.587 <2e-16 ***
## income 5.782e-06 9.853e-06 0.587 0.5573
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2050.6 on 7000 degrees of freedom
## Residual deviance: 1098.9 on 6997 degrees of freedom
## AIC: 1106.9
##
## Number of Fisher Scoring iterations: 8Посторим прогнозы модели для тестовых данных educ_test. Для этого будем использовать функцию predict, которой передадим оценённую модель educ_lmodel. В переменной educ_predz уже будут лежать предсказанные классы. Чтобы получить их вероятности, нужно добавить аргумент type = 'prob'.
educ_pred <- predict(educ_lmodel, newdata = educ_test)
head(educ_pred)## [1] L H L H H L
## Levels: H Leduc_prob <- predict(educ_lmodel, newdata = educ_test, type = 'prob')
head(educ_prob)## H L
## 25 0.1975849 0.8024151
## 31 0.8432553 0.1567447
## 37 0.3001891 0.6998109
## 45 0.8240289 0.1759711
## 50 0.7801183 0.2198817
## 55 0.2855258 0.7144742- Упражнение 3.
Постройте прогнозы для тестовой выборки def_test. Сохраните как предсказанные классы, так и их вероятности.
def_pred <- predict(def_lmodel, newdata = def_test)
head(def_pred)## [1] No No No No No No
## Levels: No Yesdef_prob <- predict(def_lmodel, newdata = def_test, type = 'prob')
head(def_prob)## No Yes
## 2 0.9989207 1.079261e-03
## 3 0.9901872 9.812833e-03
## 4 0.9995653 4.346734e-04
## 5 0.9980561 1.943937e-03
## 6 0.9981011 1.898868e-03
## 11 0.9999899 1.011236e-05Теперь мы можем посмотреть на матрицу ошибок и узнать, насколько хорошую модель мы оценили. Для этого будем использовть функцию confusionMatrix из пакета caret. В качестве аргумента data нужно указать предсказанные значения, а в reference — правильные ответы.
confusionMatrix(data = educ_pred, reference = educ_test$y)## Confusion Matrix and Statistics
##
## Reference
## Prediction H L
## H 65 8
## L 5 17
##
## Accuracy : 0.8632
## 95% CI : (0.7774, 0.9251)
## No Information Rate : 0.7368
## P-Value [Acc > NIR] : 0.002303
##
## Kappa : 0.633
## Mcnemar's Test P-Value : 0.579100
##
## Sensitivity : 0.9286
## Specificity : 0.6800
## Pos Pred Value : 0.8904
## Neg Pred Value : 0.7727
## Prevalence : 0.7368
## Detection Rate : 0.6842
## Detection Prevalence : 0.7684
## Balanced Accuracy : 0.8043
##
## 'Positive' Class : H
## - Упражнение 4.
Составьте матрицу ошибок для прогнозов модели def_lmodel.
confusionMatrix(data = def_pred, reference = def_test$default)## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 2885 65
## Yes 15 34
##
## Accuracy : 0.9733
## 95% CI : (0.9669, 0.9788)
## No Information Rate : 0.967
## P-Value [Acc > NIR] : 0.02642
##
## Kappa : 0.4474
## Mcnemar's Test P-Value : 4.293e-08
##
## Sensitivity : 0.9948
## Specificity : 0.3434
## Pos Pred Value : 0.9780
## Neg Pred Value : 0.6939
## Prevalence : 0.9670
## Detection Rate : 0.9620
## Detection Prevalence : 0.9837
## Balanced Accuracy : 0.6691
##
## 'Positive' Class : No
## Чтобы получичть и другие метрики качества для нашей модели, будем использовать функции twoClassSummary() и prSummary(). Первая возвращает значения ROC, специфичности и чувствительности.
Чувствительность (полнота) показывает долю правильно распознанных объектов отрицательного класса. \[ TPR = TP / (TP + FN) = recall = sensitivity \]
Специфичность показывает долю правильно распознанных объектов положительного класса. \[ FPR = FP / (FP + TN) = 1 - specificity \]
ROC — это площадь под кривой, построенной в осях доли правильных положительных ответов и доли неправильных положительных ответов в зависимости от значения порога для попадания в положительный класс.
Вторая функция возвращает AUC, полноту, точность и F-меру.
AUC — это площадь под кривой, построенной в осях точность-полнота. Точность — это доля объектов на самом деле принадлежащих к положительному классу среди всех объектов, которых модель отнесла к положительному классу. \[ precision = TP / (TP + FP) \]
А F-мера — это среднее гармоническое между точностью и полнотой.
Однако обе функции требуют, чтобы в данных были столбцы с вероятностями классов под названиями самих этих классов, истинные ответы под названием obs, и бинарные предсказания, названные pred. Поэтому создадим отдельную таблицу educ_test_set со всеми результатами оценивания.
educ_test_set <- data.frame(H = educ_prob$H,
L = educ_prob$L,
pred = educ_pred,
obs = educ_test$y)
glimpse(educ_test_set)## Observations: 95
## Variables: 4
## $ H <dbl> 0.1975849, 0.8432553, 0.3001891, 0.8240289, 0.7801183, 0....
## $ L <dbl> 0.802415134, 0.156744738, 0.699810922, 0.175971069, 0.219...
## $ pred <fct> L, H, L, H, H, L, H, L, L, H, L, L, H, H, L, H, H, H, H, ...
## $ obs <fct> L, H, L, H, H, L, H, L, L, H, H, L, H, H, L, H, L, H, L, ...levs <- levels(educ_test_set$obs)
twoClassSummary(educ_test_set, lev = levs) # don't work with tibble## ROC Sens Spec
## 0.9037143 0.9285714 0.6800000prSummary(educ_test_set, lev = levs) # нужен пакет MLmetrics## AUC Precision Recall F
## 0.9300243 0.8904110 0.9285714 0.9090909Что находится в векторе
levs?Упражнение 5.
Выведите характеристики модели def_lmodel с помощью функций twoClassSummary() и prSummary().
def_test_set <- data.frame(No = def_prob$No,
Yes = def_prob$Yes,
pred = def_pred,
obs = def_test$default)
twoClassSummary(def_test_set, lev = levels(def_test_set$obs))## ROC Sens Spec
## 0.9495472 0.9948276 0.3434343prSummary(def_test_set, lev = levels(def_test_set$obs))## AUC Precision Recall F
## 0.9977132 0.9779661 0.9948276 0.9863248Почти все метрики можно визуализировать! Но для примера мы построим ROC-кривую. Для этого вдобавок к базовому слою ggplot мы будем использовать слой geom_roc из пакета plotROC. В эстетиках нужно указать аргументы d — истинные значения — и m — метки положительного класса. Если добавить аргумент color, то можно получить разные ROC-кривые по категориям какой-нибудь переменной.
ggplot(educ_test_set, aes(d = obs, m = L)) +
geom_roc(n.cuts = 0) +
labs(title = 'Кривая ROC',
x = 'False positive ratio = FP / (FP + TN)',
y = 'True positive ratio = TP / (TP + FN)')
educ_test_set <- mutate(educ_test_set, gender = educ_test$gender)
ggplot(educ_test_set, aes(d = obs, m = L, color = gender)) +
geom_roc(n.cuts = 0) +
labs(title = 'Кривые ROC в зависимости от пола',
x = 'False positive ratio = FP / (FP + TN)',
y = 'True positive ratio = TP / (TP + FN)',
color = 'Пол')
- Упражнение 6.
Визуализируйте ROC-кривую для прогнозов модели def_lmodel. Сначала общую, а затем отдельно для студентов и всех остальных.
ggplot(def_test_set, aes(d = obs, m = Yes)) +
geom_roc(n.cuts = 0)
def_test_set <- mutate(def_test_set, stud = def_test$student)
ggplot(def_test_set, aes(d = obs, m = Yes, color = stud)) +
geom_roc(n.cuts = 0)
2 На природу, к деревьям и в лес!
? Пакет caret представляет собой общий интерфейс ко многим моделям, однако некоторые модели он ещё не поддерживает. Например, быстрые и стройные, fast and frugal, деревья.
Быстрые деревья нужны, например, в медицине, чтобы создавать быстрые и простые рекомендации для спасения жизни.
Мы применим их для спасения неуспевающих студентов :) Увы, пакет принимает на вход только 0/1 в зависимой переменной, поэтому заменим названия категорий на числа:
educ_train2 <- mutate(educ_train, ybin = ifelse(y == 'H', 1, 0)) И построим картинку для быстрого и стройного дерева:
fftree_model <- FFTrees(formula = ybin ~ .,
data = educ_train2)
plot(fftree_model)
- Упражнение 7.
Используйте быстрые и стройные деревья для набора данных def_train. Помните, что целевая переменная в нём закодирована как ‘Yes’ и ‘No’.
def_train2 <- mutate(def_train, defaultbin = ifelse(default == 'Yes', 1, 0))
def_fftree <- FFTrees(formula = defaultbin ~ .,
data = def_train2)
plot(def_fftree)
3 Случайный лес
В алгоритме случайного леса мы
Выращиваем целый лес, скажем 500, деревьев.
Строим прогноз с помощью каждого дерева.
Агрегируем прогнозы деревьев. Можно в качестве итогового прогноза выбрать ту категорию, за которую проголосовало большинство деревьев. Можно оценить вероятности категорий, взяв процент деревьев, проголосовавших за ту или иную категорию.
Деревья оказываются не идеальными копиями друг друга по двум причинам:
Каждое дерево обучается на случайной выборке из исходной выборки. Обычно для каждого дерева берут подвыборку с повторениями из исходной выборки, так чтобы размер подвыборки равнялся размеру исходной выборки.
При каждом делении каждой ветки на две части происходит предварительный случайный отбор переменных. Скажем, из исходных 100 переменных, каждый раз случайно отбирается 10, а затем из этих 10 выбирается наилучшая, по которой ветвь и делится на две ветви.
У идеи есть куча вариантов исполнения, отличающихся деталями:
- критерием деления ветви на две;
- критерием остановки деления дерева;
- количеством предварительно отбираемых переменных перед каждым делением;
- количество деревьев;
Посмотрим на все вариации случайного леса, которые перебрал ranger.
ranger_model <- train(y ~ ., data = educ_train,
method = 'ranger',
na.action = na.omit,
importance = 'impurity')
ranger_model## Random Forest
##
## 385 samples
## 16 predictor
## 2 classes: 'H', 'L'
##
## No pre-processing
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 385, 385, 385, 385, 385, 385, ...
## Resampling results across tuning parameters:
##
## mtry splitrule Accuracy Kappa
## 2 gini 0.8599090 0.5731322
## 2 extratrees 0.8312404 0.4562407
## 31 gini 0.9135944 0.7696853
## 31 extratrees 0.9113253 0.7635196
## 60 gini 0.9058607 0.7472557
## 60 extratrees 0.9065539 0.7498327
##
## Tuning parameter 'min.node.size' was held constant at a value of 1
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were mtry = 31, splitrule = gini
## and min.node.size = 1.plot(ranger_model,
xlab = 'Количество случайно отбираемых регрессоров',
ylab = 'Точность (бутстрэп оценка)',
main = 'Зависимость точности от настроек дерева',
auto.key = list(title = 'Алгоритм разбиения ветки', cex.title = 1))
Как мы поняли, что за название легенды отвечает аргумент auto.key?
- Мы проверили класс созданного графика —
trellis, значит, этоlattice! - Потом загуглили ‘lattice change legend title’.
- Нашли полезную статью
- Профит!
И более подробно про наилучшую:
ranger_model$finalModel## Ranger result
##
## Call:
## ranger::ranger(dependent.variable.name = ".outcome", data = x, mtry = param$mtry, min.node.size = param$min.node.size, splitrule = as.character(param$splitrule), write.forest = TRUE, probability = classProbs, ...)
##
## Type: Classification
## Number of trees: 500
## Sample size: 385
## Number of independent variables: 60
## Mtry: 31
## Target node size: 1
## Variable importance mode: impurity
## Splitrule: gini
## OOB prediction error: 8.31 %- Упражнение 8.
Реализуйте случайный лес для данных про кредиты, выведите описание лучшей модели и постройте визуализацию.
def_ranger <- train(default ~ ., data = def_train,
method = 'ranger',
importance = 'impurity')## note: only 2 unique complexity parameters in default grid. Truncating the grid to 2 .def_ranger$finalModel## Ranger result
##
## Call:
## ranger::ranger(dependent.variable.name = ".outcome", data = x, mtry = param$mtry, min.node.size = param$min.node.size, splitrule = as.character(param$splitrule), write.forest = TRUE, probability = classProbs, ...)
##
## Type: Classification
## Number of trees: 500
## Sample size: 7001
## Number of independent variables: 3
## Mtry: 2
## Target node size: 1
## Variable importance mode: impurity
## Splitrule: extratrees
## OOB prediction error: 3.16 %plot(def_ranger,
xlab = 'Количество случайно отбираемых регрессоров',
ylab = 'Точность (бутстрэп оценка)',
main = 'Зависимость точности от настроек дерева',
auto.key = list(title = 'Алгоритм разбиения ветки', cex.title = 1))
Построить информативно на графике все сотни или тысячи деревьев невозможно.
Можно попытаться выделить важность переменных:
ranger_import <- varImp(ranger_model)
ranger_import## ranger variable importance
##
## only 20 most important variables shown (out of 60)
##
## Overall
## VisITedResources 100.000
## StudentAbsenceDaysUnder-7 76.711
## raisedhands 54.587
## AnnouncementsView 29.433
## Discussion 14.961
## ParentAnsweringSurveyYes 9.844
## ParentschoolSatisfactionGood 3.194
## GradeIDG-06 2.705
## SectionIDB 2.507
## NationalITySaudiArabia 2.435
## TopicChemistry 2.311
## genderM 2.295
## NationalITyJordan 1.854
## PlaceofBirthKuwaIT 1.434
## TopicIT 1.405
## TopicScience 1.395
## TopicEnglish 1.370
## NationalITyKW 1.280
## PlaceofBirthJordan 1.138
## RelationMum 1.075plot(ranger_import,
main = 'Важность переменных случайного леса',
xlab = 'Среднее падение индекса Джини')
TODO: или суммарное?
- Упражнение 9.
Выделите важность переменных в модели def_ranger.
def_import <- varImp(def_ranger)
def_import## ranger variable importance
##
## Overall
## balance 100.00
## income 47.56
## studentYes 0.00plot(def_import)
И, конечно, построить прогнозы:
educ_ranger <- mutate(educ_test,
yhat = predict(ranger_model, educ_test, na.action = na.pass))
confusionMatrix(educ_ranger$yhat, educ_ranger$y)## Confusion Matrix and Statistics
##
## Reference
## Prediction H L
## H 68 5
## L 2 20
##
## Accuracy : 0.9263
## 95% CI : (0.8541, 0.9699)
## No Information Rate : 0.7368
## P-Value [Acc > NIR] : 2.613e-06
##
## Kappa : 0.8024
## Mcnemar's Test P-Value : 0.4497
##
## Sensitivity : 0.9714
## Specificity : 0.8000
## Pos Pred Value : 0.9315
## Neg Pred Value : 0.9091
## Prevalence : 0.7368
## Detection Rate : 0.7158
## Detection Prevalence : 0.7684
## Balanced Accuracy : 0.8857
##
## 'Positive' Class : H
## По умолчанию, пакет caret сам решает, сколько значений параметров перебирать и какие конкретно.
Список перебираемых параметров:
modelLookup(model = 'ranger')## model parameter label forReg forClass
## 1 ranger mtry #Randomly Selected Predictors TRUE TRUE
## 2 ranger splitrule Splitting Rule TRUE TRUE
## 3 ranger min.node.size Minimal Node Size TRUE TRUE
## probModel
## 1 TRUE
## 2 TRUE
## 3 TRUEНо мы можем заказать перебор любых.
Можно заказать количество перебираемых значений для каждого параметра:
ranger_model <- train(y ~ ., data = educ_test,
method = "ranger",
na.action = na.omit,
importance = 'impurity',
tuneLength = 4)Или явно значения
grid <- expand.grid(mtry = c(5, 10), min.node.size = c(1, 5), splitrule = 'gini')
ranger_model <- train(y ~ ., data = educ_test,
tuneGrid = grid, method = "ranger", na.action = na.omit)
ranger_model## Random Forest
##
## 95 samples
## 16 predictors
## 2 classes: 'H', 'L'
##
## No pre-processing
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 95, 95, 95, 95, 95, 95, ...
## Resampling results across tuning parameters:
##
## mtry min.node.size Accuracy Kappa
## 5 1 0.9478084 0.8523469
## 5 5 0.9371202 0.8239818
## 10 1 0.9501749 0.8607567
## 10 5 0.9477034 0.8524988
##
## Tuning parameter 'splitrule' was held constant at a value of gini
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were mtry = 10, splitrule = gini
## and min.node.size = 1.4 Кросс-валидация для k ближайших соседей
Метод классификации k-ближайших соседей реализует пословицу «Скажи мне, кто твой друг, и я скажу, кто ты». Мы прогнозируем класс нового наблюдения, как наиболее часто встречающийся класс у ближайших соседей данного наблюдения.
Реализуем метод k ближайших соседей для набора данных про студентов. Подбирать число соседей будем с помощью кросс-валидации. Для этого у объекта trainControl зададим опции method = 'cv' и number = 5. Первая отвечает за тип кросс-валидации, а вторая за колиество разбиений выборки. В настройках обучения train укажем метод knn, добавим параметры кросс-валидции из объекта trxrl, а также попросим центрировать и отмасштабировать данные.
trctrl <- trainControl(method = 'cv', number = 5)
set.seed(777)
educ_knn_fit <- train(y ~., data = educ_train, method = 'knn',
trControl = trctrl,
preProcess = c('center', 'scale'))
educ_knn_fit## k-Nearest Neighbors
##
## 385 samples
## 16 predictor
## 2 classes: 'H', 'L'
##
## Pre-processing: centered (60), scaled (60)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 308, 308, 308, 309, 307
## Resampling results across tuning parameters:
##
## k Accuracy Kappa
## 5 0.8285960 0.5182881
## 7 0.8129081 0.4794467
## 9 0.8259644 0.5090629
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.Визуализируем долю правильных ответов в зависимости от числа соседей.
plot(educ_knn_fit, xlab = 'Число соседей',
ylab = 'Доля правильных ответов',
main = 'Метод k ближаших соседей')
- Упражнение 10.
Реализуйте метод k ближайших соседей на данных о кредитах. Не забудьте центрировать и отмасштабировать их! Затем визуализрйуте полученный результат.
def_knn_fit <- train(default ~ ., data = def_train, method = 'knn',
trControl = trctrl,
preProcess = c('center', 'scale'))
plot(def_knn_fit, xlab = 'Число соседей',
ylab = 'Доля правильных ответов',
main = 'Метод k ближаших соседей')
Теперь получим прогнозы для тестовой выборки и оценим их качество.
educ_knn_predict <- predict(educ_knn_fit, newdata = educ_test)
confusionMatrix(educ_knn_predict, educ_test$y)## Confusion Matrix and Statistics
##
## Reference
## Prediction H L
## H 65 7
## L 5 18
##
## Accuracy : 0.8737
## 95% CI : (0.7897, 0.933)
## No Information Rate : 0.7368
## P-Value [Acc > NIR] : 0.0009533
##
## Kappa : 0.6657
## Mcnemar's Test P-Value : 0.7728300
##
## Sensitivity : 0.9286
## Specificity : 0.7200
## Pos Pred Value : 0.9028
## Neg Pred Value : 0.7826
## Prevalence : 0.7368
## Detection Rate : 0.6842
## Detection Prevalence : 0.7579
## Balanced Accuracy : 0.8243
##
## 'Positive' Class : H
## - Упражнение 11.
Постройте прогнозы метода k ближайших соседей для данных по кредитам и оцените их качество.
def_knn_pred <- predict(def_knn_fit, newdata = def_test)
confusionMatrix(def_knn_pred, def_test$default)## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 2879 64
## Yes 21 35
##
## Accuracy : 0.9717
## 95% CI : (0.9651, 0.9773)
## No Information Rate : 0.967
## P-Value [Acc > NIR] : 0.08142
##
## Kappa : 0.4382
## Mcnemar's Test P-Value : 5.225e-06
##
## Sensitivity : 0.9928
## Specificity : 0.3535
## Pos Pred Value : 0.9783
## Neg Pred Value : 0.6250
## Prevalence : 0.9670
## Detection Rate : 0.9600
## Detection Prevalence : 0.9813
## Balanced Accuracy : 0.6731
##
## 'Positive' Class : No
## 5 Кросс-валидация для LASSO
Чтобы получить простую модель с хорошими предсказаниями, можно использовать LASSO-регрессию. От обычной она отличается тем, что в целевую функцию мы добавляем дополнительное слагаемое — сумму модулей всех коэффициентов кроме первого: \[ Q = \sum_{i=1}^n (y_i - \hat y_i)^2 + \lambda \sum_{i=2}^k |\beta_i|, \] где \(\sum_{i=1}^n (y_i - \hat y_i)^2\) — сумма квадратов остатков, \(\lambda\) — параметр регуляризации. Добавляя сумму модулей коэффициентов, мы штрафуем модель за слишком большие веса. При значении параметра \(\lambda = 0\) эта модель совпадает с обычным МНК, а при росте \(\lambda\) всё большее число коэффициентов зануляется.
Поcтроим LASSO-регрессию для данных о студентах. В качестве зависимой переменной возьмём количество раз, когда студент поднимал руку, а в качестве объясняющих — все остальные. В параметрах обучения укажем method = 'lasso' и зададим явно значения параметров, среди которых мы хотим найти наилучший, аргументом tuneGrid.
trctrl <- trainControl(method = 'cv', number = 5)
educ_lasso_fit <- train(raisedhands ~ ., data = educ_train,
method = 'lasso', trControl = trctrl,
tuneGrid = expand.grid(.fraction = c(0.1, 0.2, 0.3, 0.5, 0.7)))
educ_lasso_fit## The lasso
##
## 385 samples
## 16 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 308, 306, 309, 308, 309
## Resampling results across tuning parameters:
##
## fraction RMSE Rsquared MAE
## 0.1 20.03849 0.6132161 16.62564
## 0.2 19.10408 0.6245219 15.15029
## 0.3 19.37723 0.6115496 15.32103
## 0.5 19.96130 0.5897134 15.93492
## 0.7 20.11173 0.5839026 16.10837
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was fraction = 0.2.plot(educ_lasso_fit, xlab = TeX('Доля $lambda$ от $lambda_{max}$', output = 'text'), ylab = 'RMSE', main = 'LASSO регрессия')
- Упражнение 12.
Постройте прогнозы для тестовой выборки.
educ_lasso_pred <- predict(educ_knn_fit, newdata = educ_test)6 Махалонобис
Шикарное объяснение метрики Махалонобиса от whuber.
Вспомним данные по потреблению протеинов Европе из книги Practial Machine Learning Cookbook. Загрузим их и посмотрим описательные статистики.
protein <- import('data/Europenaprotein.csv')
skim(protein)## Skim summary statistics
## n obs: 25
## n variables: 10
##
## Variable type: character
## variable missing complete n min max empty n_unique
## Country 0 25 25 2 14 0 25
##
## Variable type: numeric
## variable missing complete n mean sd p0 p25 p50 p75 p100
## Cereals 0 25 25 32.25 10.97 18.6 24.3 28 40.1 56.7
## Eggs 0 25 25 2.94 1.12 0.5 2.7 2.9 3.7 4.7
## Fish 0 25 25 4.28 3.4 0.2 2.1 3.4 5.8 14.2
## Fr&Veg 0 25 25 4.14 1.8 1.4 2.9 3.8 4.9 7.9
## Milk 0 25 25 17.11 7.11 4.9 11.1 17.6 23.3 33.7
## Nuts 0 25 25 3.07 1.99 0.7 1.5 2.4 4.7 7.8
## RedMeat 0 25 25 9.83 3.35 4.4 7.8 9.5 10.6 18
## Starch 0 25 25 4.28 1.63 0.6 3.1 4.7 5.7 6.5
## WhiteMeat 0 25 25 7.9 3.69 1.4 4.9 7.8 10.8 14
## hist
## <U+2585><U+2587><U+2581><U+2582><U+2582><U+2581><U+2581><U+2582>
## <U+2581><U+2582><U+2581><U+2581><U+2587><U+2582><U+2583><U+2582>
## <U+2585><U+2587><U+2582><U+2583><U+2581><U+2582><U+2581><U+2581>
## <U+2585><U+2587><U+2587><U+2587><U+2582><U+2581><U+2586><U+2583>
## <U+2582><U+2587><U+2582><U+2587><U+2583><U+2587><U+2581><U+2581>
## <U+2587><U+2587><U+2581><U+2583><U+2582><U+2586><U+2581><U+2581>
## <U+2582><U+2586><U+2587><U+2586><U+2581><U+2583><U+2581><U+2582>
## <U+2583><U+2581><U+2585><U+2583><U+2586><U+2587><U+2586><U+2587>
## <U+2581><U+2583><U+2587><U+2581><U+2581><U+2587><U+2583><U+2583>Отмасштабируем все числовые переменные с помощью функции scale(). Затем спрячем текстовую переменную Country в названия строк:
protein_no_country <- protein %>%
mutate_if(is.numeric, ~ as.vector(scale(.))) %>%
column_to_rownames(var = 'Country')Дополнение в виде функции as.vector нужно потому, что функция scale возвращает матрицу, а каждый столбец должен быть вектором :)
Делаем кластеризацию в два шага: находим расстояние махалонобиса, а затем рисуем картинку:
protein_dist <- distances(protein_no_country,
normalize = "mahalanobize") %>% as.matrix()
rownames(protein_dist) <- rownames(protein_no_country)
protein_hcl <- hcut(protein_dist, k = 4)
fviz_dend(protein_hcl,
cex = 0.5, # размер подписи
color_labels_by_k = TRUE) # цвет подписей по группам
Ура :)