Pavel Polishchuk, 2014
caretcaretПакет caret является удобным интерфейсом к очень многим методам машинного обучения, что в значительной мере упрощает их использование. Однако при этом теряется полный контроль над процессом построения модели. Много подробной информации содержится на официальном сайте http://caret.r-forge.r-project.org/
require(caret)
Loading required package: caret
Loading required package: lattice
Loading required package: ggplot2
Для генерации на основе рабочей выборки обучающего и тестового наборов в пакете caret используется функция createDataPartition, которая возвращает индексы обучающей выборки. Функция подходит как для классификационных задач так и для регрессионных. При этом разделение на обучающий и тестовый наборы осуществляется сбалансировано с точки зрения значения свойства.
Загрузим значения растворимости для первой выборки.
y1 <- local.load("data/sol_y1.RData")
Выделим 40% в обучающую, а остальное в тестовую выборку.
set.seed(42)
train.ids <- createDataPartition(y1, p=0.4, times=1)
Результатом будет список будет список, каждый элемент списка это вектор индексов соединений обучающей выборки.
train.ids
$Resample1
[1] 3 8 17 47 52 57 60 61 62 64 67 70 72 73 83 85 87
[18] 89 145 147 150 162 164 177 179 182 186 187 195 196 200 201 203 204
[35] 207 208 211 215 219 220 224 269 315 393 394 395 415 417 418 419 493
[52] 500 502 509 513 548 554 563 575 631 632 633 636 643 646 652 660 666
[69] 668 669 670 672 674 678 692 699 712 725 738 739 741 1 5 21 34
[86] 35 37 38 40 42 107 123 126 127 129 135 153 161 165 170 174 175
[103] 272 286 308 364 366 368 376 404 405 407 412 416 426 458 460 477 488
[120] 501 510 511 514 519 520 532 545 546 547 549 550 556 562 566 567 606
[137] 610 623 638 641 673 689 694 696 697 700 702 704 707 708 714 716 722
[154] 730 740 744 770 771 782 785 795 23 31 94 95 96 99 105 110 111
[171] 112 113 114 115 139 140 144 172 234 243 249 253 256 262 263 283 294
[188] 305 306 307 310 311 328 330 334 337 340 341 358 361 422 433 439 446
[205] 447 448 455 483 485 486 494 505 518 521 523 527 529 536 542 558 559
[222] 565 572 607 608 613 618 624 626 650 654 656 709 718 724 746 754 758
[239] 759 769 774 783 92 226 229 230 233 236 237 245 246 247 251 259 271
[256] 275 279 281 293 297 298 303 309 313 318 320 321 323 336 346 348 351
[273] 352 353 355 367 382 383 384 386 390 401 403 430 442 449 461 465 466
[290] 467 472 479 499 576 577 579 580 583 584 586 600 601 602 604 605 619
[307] 653 661 682 715 731 752 755 761 762 763 773 777 779 786 790 800
В нашем случае в списке только один элемент - одна обучающая выборка. Если необходимо создать несколько пар обучающих/тестовых выборок, меняем значение times на нужное число.
Поскольку у нас только одна тестовая выборка, то для упрощения дальнейшего кода преобразуем train.ids из списка в вектор индексов.
train.ids <- train.ids[[1]]
Загрузим дескрипторы для первой выборки и разделим всю выборку на обучающую и тестовую
x1 <- local.load("data/sol_x1.RData")
x1.train <- x1[train.ids,]
x1.test <- x1[-train.ids,]
y1.train <- y1[train.ids]
y1.test <- y1[-train.ids]
Основная идея - создание максимально представительной тестовой выборки. Это достигается выбором в тестовую выборку соединений максимально отличающихся друг от друга. Процедура итеративна, на каждом шаге к тестовому набору добавляется одно соединение, которое максимально отличается от соединений, которые уже включены в тестовую выборку.
В качестве меры подобия по умолчанию используется эвклидово расстояние, поэтому чтобы получить адекватные результаты предварительно необходимо шкалировать дескрипторы всей рабочей выборки для чего используем стандартную функцию scale.
x1.scaled <- scale(x1)
Выберем случайным образом несколько соединений и поместим их в начальную тестовую выборку, а остальные в обучающую
set.seed(24)
startSet <- sample(1:nrow(x1.scaled), 5)
x1.test.diss <- x1.scaled[startSet, ]
x1.train.diss <- x1.scaled[-startSet, ]
Вызовем функцию maxDissim (требует установленного пакета proxy), в которой передадим тестовый и обучающий наборы соединений и также укажем количество соединений, которое мы хотим переместить из обучающего набора в тестовый, например 20% (\( ^1/_5 \) от 800 = 160, с учетом того что 5 соединений мы уже выбрали следует указать 155). Чтобы ускорить выполнение функции укажем параметр randomFrac, который определяет какая доля обучающей выборки будет рассматриваться на каждом шаге поиска подобных соединений.
newSamples <- maxDissim(x1.test.diss, x1.train.diss, n = 155, randomFrac = 0.1)
Attaching package: 'proxy'
Следующий объект скрыт from 'package:stats':
as.dist, dist
newSamples
[1] 530 776 562 561 508 539 769 694 687 664 782 703 692 699 524 673 560
[18] 672 630 145 698 655 749 559 515 511 529 516 505 333 643 527 504 146
[35] 507 543 657 691 552 777 476 701 765 514 541 554 707 149 615 750 535
[52] 778 665 705 781 150 661 411 575 81 660 761 760 422 674 621 489 675
[69] 790 501 219 685 547 538 647 666 741 794 546 537 717 606 406 548 753
[86] 767 582 658 479 580 745 570 708 627 677 791 147 343 667 686 220 779
[103] 437 754 563 551 668 773 417 702 404 542 774 639 80 566 733 413 472
[120] 642 795 136 579 662 626 500 482 595 216 793 488 525 376 496 495 484
[137] 713 399 390 764 521 771 469 629 477 415 768 209 735 736 467 576 549
[154] 493 789
Перенесем эти соединения из обучающей выборки в тестовую
x1.test.diss <- rbind(x1.test.diss, x1.train.diss[newSamples,])
x1.train.diss <- x1.train.diss[-newSamples,]
Кроме евклидового расстояния в качестве меры подобия можно использовать другие метрики реализованные в пакете proxy, например манхэттенское расстояние (полный список можно посмотреть выполнив summary(proxy::pr_DB))
# Not run
newSamples_Manhattan <- maxDissim(x1.test.diss, x1.train.diss, n = 20, method="Manhattan", randomFrac = 0.1)
Наличие в выборке дескрипторов с нулевой вариацией может не совсем корректно обрабатываться некоторыми методами, поэтому такие переменные лучше сразу исключить из рассмотрения. Также переменные, имеющие малое число уникальных значений, могут рассматриваться как малоинформативные и тоже могут быть исключены из рассмотрения.
Воспользуемся функцией nearZeroVar, с помощью которой проверим набор данных на наличие дескрипторов с постоянными значениями
nzv <- nearZeroVar(x1.train, freqCut=0, uniqueCut=0)
freqCut - отношение количества наиболее часто встречающегося значения переменной к количеству второго по частоте встречаемости значения переменной \( \frac{N(x_1)}{N(x_2)} \)
uniqueCut - отношение числа уникальных значений переменной к общему числу объектов в выборке в процентах \( \frac{N_{unique}}{N_{total}} \times 100 \)
Если значение freqCut больше заданной границы И значение uniqueCut меньше, то переменная будет отмечена как имеющая малый разброс значений. Поэтому, чтобы отметить только переменные с постоянными значениями, достаточно установить uniqueCut = 0.
nzv
[1] 7 15 38 45 46 51 55 56 57 65 83 86 88 91
[15] 92 104 117 128 131 132 133 135 139 140 149 154 166 177
[29] 183 184 185 190 196 204 209 211 238 248 250 258 260 266
[43] 276 294 296 297 326 328 330 334 343 345 353 363 373 379
[57] 381 394 398 399 405 409 431 433 443 444 446 458 461 470
[71] 472 491 494 536 540 544 552 558 559 563 570 582 616 621
[85] 658 674 675 677 682 689 693 696 697 704 721 729 732 733
[99] 751 756 757 759 762 766 769 775 784 797 811 819 825 840
[113] 846 857 860 869 880 884 899 907 909 915 940 948 950 952
[127] 954 964 967 972 973 975 997 1008 1014 1016 1030 1035 1037 1039
[141] 1041 1046 1047 1052 1053 1059 1060 1068 1069 1086 1088 1104 1106 1122
[155] 1138 1139 1140 1155 1158 1171 1178 1190 1205 1219 1240 1255 1263 1266
[169] 1277 1297 1314 1330 1340 1341 1342 1351 1355 1356 1359 1369 1373 1376
[183] 1411 1412 1418 1419 1421 1426 1434 1437 1450 1476 1487 1506 1516 1517
[197] 1527 1528 1549 1565 1582 1596 1598 1601 1608 1614 1622 1633 1636 1643
[211] 1645 1649 1669 1670 1676 1678 1686 1689 1690 1696 1709 1721 1731 1735
[225] 1745 1750 1752 1754 1756 1759 1768 1772 1774 1775 1782 1791 1795 1802
[239] 1803 1805 1810 1814 1816 1818 1820 1828 1831 1832 1836 1840 1842 1845
[253] 1848 1852 1861 1870 1885 1893 1911 1947 1948 1949 1957 1966 1969 1976
[267] 1978 1979 1984 2001 2021 2035 2055 2059 2063 2074 2077 2080 2094 2097
[281] 2102 2108 2129 2132 2143 2152 2161 2170 2179 2183 2188 2198 2199 2208
[295] 2211 2223 2224 2228 2234 2238 2254 2260 2262 2263 2264 2273 2277 2291
[309] 2292 2295 2297 2302 2303 2304 2307 2310 2312 2314 2319 2336 2362 2366
[323] 2367 2382 2393 2409 2420 2428 2433 2435 2436 2439 2455 2471 2492 2498
[337] 2514 2524 2528 2530 2531 2535 2536 2537 2550 2552 2560 2594 2597 2602
[351] 2603 2613 2617 2645 2646 2647 2648 2652 2656 2680 2681 2684 2707 2711
[365] 2714 2717 2722 2737 2743 2757 2758 2763 2779 2783 2789 2790 2806 2807
[379] 2819 2845 2888 2898 2901 2917 2920 2925 2945 2947 2951 2970 2982 2989
[393] 2990 3006 3007 3041 3048 3061 3113 3114 3116 3123 3138 3150 3173 3191
[407] 3233 3237 3259 3271 3281 3287 3290 3309 3320 3326 3330 3341 3342 3349
[421] 3409 3413 3416 3431 3433 3435 3437 3442 3458 3461 3473 3478 3498 3500
[435] 3501 3517 3523 3570 3577 3608 3638 3644 3650 3673 3677 3679 3693 3695
[449] 3701 3711 3739 3743 3749 3800 3803 3804 3805 3810 3820 3821 3822 3836
[463] 3837 3838 3839 3840 3841 3849 3851 3866 3877 3879 3884 3899 3907 3909
[477] 3913 3927 3933 3947 3956 3959 3968 3971 3972 3977 3980 3982 3983 3986
[491] 3990 3997 3998 4000 4001 4011 4015 4016 4021 4022 4023 4024 4025 4028
[505] 4032 4033 4034 4050 4051
Т.е. в нашем случае для соединений обучающей выборки есть 509 дескрипторов с постоянными значениями, которые нужно удалить.
x1.train <- x1.train[,-nzv]
Иногда для получения надежной модели необходимо удалить переменные с высокой взаимной корреляцией. Для чего можно воспользоваться функцией findCorrelation из пакета caret.
highCor <- findCorrelation(cor(x1.train), cutoff=0.9)
head(highCor, 10)
[1] 1836 1795 3418 1825 3303 1827 857 3398 3320 1706
length(highCor)
[1] 2154
highCor содержит индексы переменных с высокой взаимной корреляцией, которые можно исключить из рассмотрения при построении модели.
x1.train <- x1.train[,-highCor]
dim(x1.train)
[1] 322 1395
Если в выборке, например, есть переменная, которая может быть выражена через сумму двух других (или нескольких) переменных, то такую переменную можно рассматривать как малоинформативную и несущую дублирующую информацию и ее можно исключить из построения модели. Для того, чтобы найти переменные связанные такими линейными зависимостями используется функция findLinearCombos.
linCombo <- findLinearCombos(x1.train)
linCombo это список из двух элементов:
linCombo$linearCombos - список найденных линейных комбинаций.
linCombo$remove - вектор индексов переменных, которые можно выразить через линейную комбинацию остальных переменных.
Т.е. всего было найдено
length(linCombo$remove)
[1] 1094
переменных, которые можно исключить из построения модели, воспользовавшись свойством числовых индексов
Шкалирование (или масштабирование) переменных обучающей выборки часто используется перед построением модели, т.к. если абсолютные значения переменных сильно отличаются между собой, то переменные с большими абсолютными значениями искусственно получат больший вес и значимость при построении некоторых моделей.
Наиболее распространенным преобразованием является центрирование и нормирование на величину стандартного отклонения
\[ x_{i,scaled}=\frac{x_i-\overline{x}}{sd(x)} \]
С этой целью можно использовать функцию scale из базового пакета, или функцию preProcess из пакета caret. Удобнее использовать второй вариант, т.к. объект создаваемый функцией preProcess содержит информацию о параметрах шкалирования (например, величину среднего значения и величину стандартного отклонения для каждой переменной) и его можно легко применять к другим наборам данных, для которых необходимо спрогнозировать значения свойства по построенным моделям.
Зададим способ нормировки переменных
preProc <- preProcess(x1.train, method=c("center", "scale"))
Чтобы не вызывать каждый раз эту функцию создавая объект класса preProcess, его лучше сохранить и потом загружать при необходимости.
save(preProc, file="data/sol_x1_preprocess.RData")
Масштабируем переменные обучающей выборки и сохраним их в файл
x1.train <- predict(preProc, x1.train)
save(x1.train, file="data/sol_x1_train_scaled.RData")
Опция method функции preProcesss может принимать и другие значения, которые позволяют преобразовывать данные. В частности значение method = "range" приводит значения переменных к диапазону [0,1].
caretОсновные возможности caret: поиск и выбор оптимальных параметров модели и построение конечной модели. На сегодня в пакете caret можно строить модели, используя около 150 различных методов. Помимо них можно самостоятельно реализовать поддержку любых других методов, не вошедших в это число.
Основной алгоритм
Сперва определяются параметры, которые следует оптимизировать. После чего выполняется построение моделей и на основе оценок, получаемых в ходе валидации моделей, определяют оптимальные параметры построения модели, которые уже используются для построения конечной модели на всем наборе данных.
Сперва указываются параметры процедуры обучения модели. В данном примере это 5-кратная кросс-валидация, повторенная дважды.
trControl1 <- trainControl(method = "repeatedcv",
number = 5,
repeats = 2,
preProcOptions = NULL,
savePredictions = TRUE)
Опция savePredictions позволяет сохранить значения, спрогнозированные в ходе кросс-валидации.
Теперь с использованием функции train строится модель, где указываются наборы значений X и Y для обучающей выборки, указывается метод построения модели и передается параметр trControl.
m1.pls <- train(x1.train,
y1[train.ids],
method = 'pls',
trControl = trControl1)
m1.pls
Partial Least Squares
322 samples
1395 predictors
No pre-processing
Resampling: Cross-Validated (5 fold, repeated 2 times)
Summary of sample sizes: 257, 258, 257, 258, 258, 258, ...
Resampling results across tuning parameters:
ncomp RMSE Rsquared RMSE SD Rsquared SD
1 1.72 0.387 0.235 0.0802
2 1.49 0.54 0.172 0.0714
3 1.31 0.642 0.182 0.0781
RMSE was used to select the optimal model using the smallest value.
The final value used for the model was ncomp = 3.
В результате кросс-валидации показано как меняются статистические характеристики модели с изменением числа компонент в PLS уравнении. Оптимальным числом компонент в данном случае является 3. Именно для этого числа компонент и построена итоговая модель, которая хранится в объекте m1.pls и которую можно в дальнейшем использовать для прогноза новых выборок.
Результат можно вывести в виде графика, показывающего зависимости ошибки модели от настроечных параметров модели
plot(m1.pls)
Чтобы спрогнозировать значения растворимости для тестовой выборки дескрипторы тестового набора необходимо масштабировать с использованием параметров, полученных для обучающей выборки. Предварительно нужно оставить только те дескрипторы, которые присутствовали в обучающей выборке, и расположить их в той же последовательности. Для этого можно воспользоваться тем, что объект класса preProcess хранит имена дескрипторов.
preProc <- local.load("data/sol_x1_preprocess.RData")
x1.test <- x1.test[,names(preProc$mean)]
x1.test <- predict(preProc, x1.test)
В общем случае необходимо делать проверку на
NA, т.к. какие-то дескрипторы возможно отсутствуют в тестовом наборе.
Сам прогноз выполняется функцией predict
pred <- predict(m1.pls, x1.test)
Определим среднеквадратичное отклонение значения прогноза для внешней тестовой выборки
Metrics::rmse(y1.test, pred)
[1] 1.298963
trainПод коэффициентом детерминации, возвращаемом функцией train подразумевается обычный квадрат коэффициента корреляции Пирсона, а не рассчитываемый по формуле
\[ R^2 = 1 - \frac{PRESS}{SS} = 1 - \frac{\sum_{i=1}^{n}(\hat{y}_i-y_i)^2}{\sum_{i=1}^{n}(y_i-\overline{y}_{train})^2} \]
Если необходимо вычислить коэффициент корреляции по этой формуле, это надо делать вручную, т.е. написать функцию. Значения, спрогнозированные в ходе кросс-валидации хранятся в объекте m1.pls$pred.
head(m1.pls$pred)
pred obs rowIndex ncomp Resample
1 -1.771647 -3.84 1 3 Fold1.Rep1
2 -1.995828 -3.82 3 3 Fold1.Rep1
3 -2.946894 -3.96 7 3 Fold1.Rep1
4 -2.893327 -4.89 9 3 Fold1.Rep1
5 -2.686042 -4.22 15 3 Fold1.Rep1
6 -3.202500 -6.89 16 3 Fold1.Rep1
Написать функцию, которая бы принимала в качестве единственного параметра объект класса train (m1.pls), и возвращала для выбранной модели значение \( R^2 \), рассчитанное по формуле приведенной выше.
При валидации модели можно задавать наборы соединений в каждом из фодов. Сгенерировать соединения, которые войдут в каждый из фолдов можно либо с использованием встроенных функций createFolds и createMultiFolds. Или, если валидации требует более сложной процедуры генерации фолдов, создав собственную функцию, которая бы возвращала список индексов соединений обучающей выборки каждого фолда.
set.seed(42)
# создание набора фолдов для 5-кратной кросс-валидации
cv <- createFolds(y1.train, k=5, returnTrain=TRUE)
# или создание набора фолдов для 5-кратной кросс-валидации повторенной трижды
cv2 <- createMultiFolds(y1.train, 5, 3)
Фолды это списки с индексами соединений обучающей выборки для каждого фолда
summary(cv)
Length Class Mode
Fold1 257 -none- numeric
Fold2 258 -none- numeric
Fold3 258 -none- numeric
Fold4 258 -none- numeric
Fold5 257 -none- numeric
summary(cv2)
Length Class Mode
Fold1.Rep1 258 -none- numeric
Fold2.Rep1 258 -none- numeric
Fold3.Rep1 257 -none- numeric
Fold4.Rep1 257 -none- numeric
Fold5.Rep1 258 -none- numeric
Fold1.Rep2 258 -none- numeric
Fold2.Rep2 257 -none- numeric
Fold3.Rep2 258 -none- numeric
Fold4.Rep2 258 -none- numeric
Fold5.Rep2 257 -none- numeric
Fold1.Rep3 258 -none- numeric
Fold2.Rep3 258 -none- numeric
Fold3.Rep3 258 -none- numeric
Fold4.Rep3 257 -none- numeric
Fold5.Rep3 257 -none- numeric
Созданные один раз фолды можно сохранить.
Еще одна возможность повлиять на процедуру моделирования заключается в указании значений настроечных параметров, которые будут использоваться при оптимизации. Для этого необходимо создать data.frame, колонки которого называются также как и настроечные параметры модели. Список параметров, которые можно оптимизировать для каждой модели, приведен на странице http://caret.r-forge.r-project.org/modelList.html
Для pls метода возможно оптимизация только числа компонент в уравнении ncomp
plsGrid <- data.frame(ncomp=1:7)
plsGrid
ncomp
1 1
2 2
3 3
4 4
5 5
6 6
7 7
Теперь задаем параметры обучения модели через функцию trainControl, куда передаем индексы фолдов, а в качестве метода указываем LGOCV. Последнее не принципиально, т.к. при передаче параметра index, параметр method уже не учитывается. Но из соображений удобства лучше использовать LGOCV.
trControl2 <- trainControl(method="LGOCV",
index=cv,
preProcOptions=NULL,
savePredictions = TRUE)
m2.pls <- train(x1.train,
y1.train,
method = 'pls',
trControl = trControl2,
tuneGrid = plsGrid)
m2.pls
Partial Least Squares
322 samples
1395 predictors
No pre-processing
Resampling: Repeated Train/Test Splits Estimated (5 reps, 0.75%)
Summary of sample sizes: 257, 258, 258, 258, 257
Resampling results across tuning parameters:
ncomp RMSE Rsquared RMSE SD Rsquared SD
1 1.71 0.397 0.162 0.0858
2 1.5 0.528 0.143 0.0756
3 1.33 0.635 0.121 0.0468
4 1.21 0.698 0.132 0.0437
5 1.16 0.721 0.0965 0.044
6 1.12 0.742 0.0858 0.0393
7 1.08 0.76 0.112 0.0521
RMSE was used to select the optimal model using the smallest value.
The final value used for the model was ncomp = 7.
plot(m2.pls)
Для ускорения процедуры кросс-валидации, которая может выполняться независимо для каждого фолда, следует использовать возможности многопоточности. Для этого предварительно инициализируем кластер из заданного числа узлов.
# Windows / Unix
require(doParallel)
registerDoParallel(2)
# Unix
library(doMC)
registerDoMC(2)
Параллельное выполнение кода приводит к дублированию в памяти используемых данных, поэтому при запуске построения модели в многопоточном режиме надо учитывать этот факт. Использование пакета doMC возможно только в Unix-системах и дает возможность избежать многократного дублирования данных в памяти.
Теперь процесс построения модели по умолчанию будет использовать несколько потоков. Сравним два подхода.
trControl2 <- trainControl(method="LGOCV",
index=cv2,
preProcOptions=NULL,
savePredictions=TRUE,
allowParallel=FALSE)
system.time(m2.pls <- train(x1.train,
y1.train,
method = 'pls',
trControl = trControl2,
tuneGrid = plsGrid))
user system elapsed
21.81 2.89 24.74
trControl3 <- trainControl(method="LGOCV",
index=cv2,
preProcOptions=NULL,
savePredictions=TRUE,
allowParallel=TRUE)
system.time(m3.pls <- train(x1.train,
y1.train,
method = 'pls',
trControl = trControl3,
tuneGrid = plsGrid))
user system elapsed
1.45 0.16 20.37
Выигрыш по времени не так велик, как можно было ожидать. Этому есть несколько причин, которые необходимо учитывать при построении моделей. Основные это расход времени на создание потоков и копирование в них необходимых данных. Даже если памяти в компьютере достаточно эта процедура может быть сравнительно медленной. Поэтому, если модели строятся быстро сами по себе, то польза от их построения в многопоточном режиме не будеть очень заметной. Наибольшего эффекта можно достигнуть, когда построение одной модели достаточно длительно. Тогда затраты на создание потоков и копирование данных будут не так заметны.
all.equal(m2.pls, m3.pls)
[1] "Component 7: target, current do not match when deparsed"
[2] "Component 10: Component 21: 1 element mismatch"
[3] "Component 11: Component 14: Mean relative difference: 0.1"
[4] "Component 19: Component 1: Mean relative difference: 0.5611394"
[5] "Component 19: Component 2: Mean relative difference: 0.04477612"
Запущенный кластер можно остановить в любой момент, используя функцию stopImplicitCluster (только в случае использования пакета doParallel).
# if cluster was initialized via registerDoParallel
stopImplicitCluster()
В пакете caret предусмотрена функция varImp, которая возвращает значения важности переменных. Эта функция содержит несколько параметров.
varImp(object, useModel = TRUE, nonpara = TRUE, scale = TRUE, ...)
object - объект класса train, который возвращает одноименная функция из пакета caret
useModel - следует ли использовать оценку важности переменных используя построенную модель. Возможно для ряда моделей, которые имеют собственный способ оценки важности переменных (randomForest, pls, lm, gbm и т.д., полный список доступен в справке по описываемой функции).
nonpara - если для модели нет собственноо алгоритма расчета важности переменных, или useModel = FALSE, то использовать непераметрическую или параметрическую оценку важности переменных.
scale - надо ли шкалировать значения важности переменных в диапазоне от 0 до 100 или нет.
Определим важности переменных в модели m1.pls
vi1.pls <- varImp(m1.pls)
vi1.pls
pls variable importance
only 20 most important variables shown (out of 1395)
Overall
`S_A(d_a)/I_I_I_I/1_3a,2_4a,3_4a/6` 100.00
`S_A(d_a)/I_I_I_I/1_2s,1_3a,2_4a/6` 82.94
`S_A(chg)/B_B_B_B/1_4a,2_4a,3_4a/5` 73.21
`S_A(lip)/B_C_D_D/1_2a,1_4s/4` 73.18
`S_A(rf)/B_B_B_C/1_2s,2_3a,3_4s/6` 69.44
`S_A(rf)/B_B_B_B/1_2s,1_3a,2_4a/6` 67.44
`S_A(lip)/B_B_C_C/1_3a,2_4a,3_4s/6` 64.77
`S_A(chg)/B_B_B_B/2_3s,3_4a/4` 64.62
`S_A(lip)/B_B_C_D/1_2a,1_4s,2_3a/6` 61.19
`S_A(chg)/B_B_B_B/1_2s,3_4a/3` 59.99
`S_A(rf)/A_B_B_C/1_2s,2_3a/4` 59.69
`S_A(chg)/D_D_D_D/1_3a,2_4a,3_4a/6` 59.00
`S_A(chg)/C_C_D_D/1_2s,1_3a,2_4a/6` 58.54
`S_A(lip)/B_B_C_D/1_3a,1_4s,2_3a/6` 56.34
`S_A(rf)/B_B_B_B/1_2s,2_4a,3_4a/6` 56.30
`S_A(lip)/C_C_C_D/1_2s,2_3a/4` 53.47
`S_A(rf)/A_B_B_C/1_2s,2_3a,3_4s/6` 52.82
`S_A(rf)/A_B_B_C/2_3a,2_4s/4` 52.26
`S_A(lip)/B_B_B_C/1_2a,2_3a,3_4a/6` 51.78
`S_A(lip)/B_B_C_C/1_2a,2_4a,3_4s/6` 50.95
Получившиеся значения щкалированы в диапазоне [0-100], чтобы вывести оригинальные значения зададим параметр scale = FALSE
vi2.pls <- varImp(m1.pls, scale = FALSE)
vi2.pls
pls variable importance
only 20 most important variables shown (out of 1395)
Overall
`S_A(d_a)/I_I_I_I/1_3a,2_4a,3_4a/6` 0.03006
`S_A(d_a)/I_I_I_I/1_2s,1_3a,2_4a/6` 0.02496
`S_A(chg)/B_B_B_B/1_4a,2_4a,3_4a/5` 0.02204
`S_A(lip)/B_C_D_D/1_2a,1_4s/4` 0.02203
`S_A(rf)/B_B_B_C/1_2s,2_3a,3_4s/6` 0.02091
`S_A(rf)/B_B_B_B/1_2s,1_3a,2_4a/6` 0.02032
`S_A(lip)/B_B_C_C/1_3a,2_4a,3_4s/6` 0.01951
`S_A(chg)/B_B_B_B/2_3s,3_4a/4` 0.01947
`S_A(lip)/B_B_C_D/1_2a,1_4s,2_3a/6` 0.01844
`S_A(chg)/B_B_B_B/1_2s,3_4a/3` 0.01808
`S_A(rf)/A_B_B_C/1_2s,2_3a/4` 0.01800
`S_A(chg)/D_D_D_D/1_3a,2_4a,3_4a/6` 0.01779
`S_A(chg)/C_C_D_D/1_2s,1_3a,2_4a/6` 0.01765
`S_A(lip)/B_B_C_D/1_3a,1_4s,2_3a/6` 0.01699
`S_A(rf)/B_B_B_B/1_2s,2_4a,3_4a/6` 0.01698
`S_A(lip)/C_C_C_D/1_2s,2_3a/4` 0.01613
`S_A(rf)/A_B_B_C/1_2s,2_3a,3_4s/6` 0.01594
`S_A(rf)/A_B_B_C/2_3a,2_4s/4` 0.01577
`S_A(lip)/B_B_B_C/1_2a,2_3a,3_4a/6` 0.01563
`S_A(lip)/B_B_C_C/1_2a,2_4a,3_4s/6` 0.01538
Универсальные значения важности переменных, можно получить если задать параметр useModel = FALSE или применять функцию varImp к модели, для которой нет собственной реализации расчета важности переменных. В этом случае результаты уже не будут зависеть от самой модели.
vi3.pls <- varImp(m1.pls, useModel = FALSE, scale = FALSE, nonpara = FALSE)
vi3.pls
Linear model variable importance
only 20 most important variables shown (out of 1395)
Overall
S_A(d_a)/I_I_I_I/1_3a,2_4a,3_4a/6 13.045
S_A(d_a)/I_I_I_I/1_2s,1_3a,2_4a/6 11.980
S_A(lip)/B_C_D_D/1_2a,1_4s/4 11.798
S_A(rf)/B_B_B_C/1_2s,2_3a,3_4s/6 11.193
S_A(rf)/A_B_B_C/1_2s,2_3a/4 10.697
S_A(rf)/B_B_B_B/1_2s,1_3a,2_4a/6 9.712
S_A(lip)/B_B_C_D/1_2a,1_4s,2_3a/6 9.587
S_A(lip)/C_C_C_D/1_2s,2_3a/4 9.251
S_A(rf)/A_B_B_C/1_2s,2_3a,3_4s/6 9.185
S_A(rf)/A_B_B_C/2_3a,2_4s/4 9.126
S_A(lip)/B_B_C_C/2_4a,3_4s/4 8.875
S_A(lip)/B_B_C_C/1_3a,2_4a,3_4s/6 8.780
S_A(lip)/B_C_C_D/1_3a,2_3a/4 8.498
S_A(lip)/B_B_C_C/1_2a,2_4a,3_4s/6 8.465
S_A(lip)/B_C_C_D/1_3a,1_4s/4 8.461
S_A(lip)/B_B_C_D/1_3a,1_4s,2_3a/6 8.360
S_A(lip)/B_C_C_D/1_3a,1_4s,2_3a/6 8.200
S_A(rf)/B_B_B_B/1_2s,2_4a,3_4a/6 8.086
S_A(lip)/B_B_B_C/2_3a,3_4a/4 7.912
S_A(lip)/B_B_C_C/1_3a,2_4a/3 7.672
В данном случае использовался параметрический подход, который заключается в том, что строится линейная модель, и важность переменной выражается как абсолютное значение t-критерия (коэффициента Стьюдента) для коэффициента перед переменной в линейном уравнении.
Если использовать непараметрический подход, то важность переменных в этом случае представляется как величина коэффициента детерминации для локальной регрессии (loess - locally weighted scatterplot smoothing), включающей одну переменную.
vi4.pls <- varImp(m1.pls, useModel = FALSE, scale = FALSE, nonpara = TRUE)
vi4.pls
loess r-squared variable importance
only 20 most important variables shown (out of 1395)
Overall
S_A(d_a)/I_I_I_I/1_3a,2_4a,3_4a/6 0.3472
S_A(d_a)/I_I_I_I/1_2s,1_3a,2_4a/6 0.3096
S_A(lip)/B_C_D_D/1_2a,1_4s/4 0.3031
S_A(rf)/B_B_B_C/1_2s,2_3a,3_4s/6 0.2813
S_A(rf)/A_B_B_C/1_2s,2_3a/4 0.2634
S_A(rf)/B_B_B_B/1_2s,1_3a,2_4a/6 0.2277
S_A(lip)/B_B_C_D/1_2a,1_4s,2_3a/6 0.2231
S_A(lip)/C_C_C_D/1_2s,2_3a/4 0.2110
S_A(rf)/A_B_B_C/1_2s,2_3a,3_4s/6 0.2086
S_A(rf)/A_B_B_C/2_3a,2_4s/4 0.2065
S_A(lip)/B_B_C_C/2_4a,3_4s/4 0.1975
S_A(lip)/B_B_C_C/1_3a,2_4a,3_4s/6 0.1941
S_A(lip)/B_C_C_D/1_3a,2_3a/4 0.1841
S_A(lip)/B_B_C_C/1_2a,2_4a,3_4s/6 0.1829
S_A(lip)/B_C_C_D/1_3a,1_4s/4 0.1828
S_A(lip)/B_B_C_D/1_3a,1_4s,2_3a/6 0.1792
S_A(lip)/B_C_C_D/1_3a,1_4s,2_3a/6 0.1736
S_A(d_a)/A_I_I_I/1_4s,3_4s/4 0.1729
S_A(rf)/B_B_B_B/1_2s,2_4a,3_4a/6 0.1696
S_A(lip)/B_B_B_C/2_3a,3_4a/4 0.1636
Как видно во всех случаях некоторые переменные ранжируются приблизительно одинаково с точки зрения их важности, независимо от выбранного метода оценки важности.
По умолчанию при построении и выборе оптимальной модели качество регрессионных моделей оценивается по величине RMSE, а классификационных по точности прогноза. Однако в ряде случаев может потребоваться другая оценочная функция, например MAE в случае регрессионных задач.
Функция trainControl среди прочих имеет параметр summaryFunction, который по умолчанию равен defaultSummary - функция, которая возвращает статистику для кросс-валидации. Можно создать собственную функцию расчета статистических показателей и передать в качестве параметра summaryFunction.
Создадим функцию, которая будет возвращать в качестве статистических показателей величины \( MAE \), \( RMSE \) и \( R^2 \). Функция должна содержать три обязательных параметра. data представляет собой data.frame , который содержит колонки obs и pred, для наблюдаемых значений и значений предсказанных в ходе процедура кросс-валидации. Результатом функции должен быть вектор значений с именами.
summaryMAE <- function(data, lev=NULL, model=NULL) {
r2 <- cor(data$obs, data$pred) ^ 2
rmse <- sum((data$obs - data$pred) ^ 2) / nrow(data)
mae <- sum(abs(data$obs - data$pred)) / nrow(data)
out <- c(mae, rmse, r2)
names(out) <- c("MAE", "RMSE", "R2")
out
}
Создадим новый объект trainControl, включающий созданную функцию, и построим модель.
trControl.mae <- trainControl(method="LGOCV",
index=cv,
savePredictions=TRUE,
preProcOptions=NULL,
summaryFunction=summaryMAE)
m4.pls <- train(x1.train,
y1.train,
method="pls",
trControl=trControl.mae,
tuneGrid=plsGrid)
m4.pls
Partial Least Squares
322 samples
1395 predictors
No pre-processing
Resampling: Repeated Train/Test Splits Estimated (5 reps, 0.75%)
Summary of sample sizes: 257, 258, 258, 258, 257
Resampling results across tuning parameters:
ncomp MAE RMSE R2 MAE SD RMSE SD R2 SD
1 1.29 2.95 0.397 0.0961 0.576 0.0858
2 1.15 2.28 0.528 0.099 0.437 0.0756
3 1 1.79 0.635 0.0704 0.323 0.0468
4 0.914 1.49 0.698 0.0979 0.318 0.0437
5 0.862 1.35 0.721 0.101 0.219 0.044
6 0.837 1.25 0.742 0.0918 0.186 0.0393
7 0.815 1.17 0.76 0.114 0.241 0.0521
RMSE was used to select the optimal model using the smallest value.
The final value used for the model was ncomp = 7.
В списке результатов появилось новое значение, однако выбор оптимальной модели происходит на основании величины RMSE. Если необходимо использовать для выбора оптимальной модели величину MAE, укажем это в качестве параметра функции train
m5.pls <- train(x1.train,
y1.train,
method="pls",
trControl=trControl.mae,
tuneGrid=plsGrid,
metric="MAE",
maximize=FALSE)
m5.pls
Partial Least Squares
322 samples
1395 predictors
No pre-processing
Resampling: Repeated Train/Test Splits Estimated (5 reps, 0.75%)
Summary of sample sizes: 257, 258, 258, 258, 257
Resampling results across tuning parameters:
ncomp MAE RMSE R2 MAE SD RMSE SD R2 SD
1 1.29 2.95 0.397 0.0961 0.576 0.0858
2 1.15 2.28 0.528 0.099 0.437 0.0756
3 1 1.79 0.635 0.0704 0.323 0.0468
4 0.914 1.49 0.698 0.0979 0.318 0.0437
5 0.862 1.35 0.721 0.101 0.219 0.044
6 0.837 1.25 0.742 0.0918 0.186 0.0393
7 0.815 1.17 0.76 0.114 0.241 0.0521
MAE was used to select the optimal model using the smallest value.
The final value used for the model was ncomp = 7.
В этом случае алгоритм выбрал ту же модель с числом компонент 7, но уже на основании величины MAE
sol_x1.RData построить модели методами pls, knn, rf, gbm, svmRadial, avNNet.sol_x2.RData.