Pavel Polishchuk, 2014
Консенсусное моделирование.
1. Усреднение результатов прогнозов нескольких моделей.
2. Стекинг результатов прогнозов нескольких моделей.
Загрузим данные для первой выборки по растворимости, построим несколько моделей методами pls, randomForest, knn, gbm и svm и сохраним их. На сайте можно найти уже построенные модели.
x <- local.load("data/sol_x1.RData")
y <- local.load("data/sol_y1.RData")
require(caret)
require(doParallel)
registerDoParallel(3)
preProc <- preProcess(x, method=c("scale", "center"))
x <- predict(preProc, x)
set.seed(42)
cv <- createFolds(y, 5, returnTrain=TRUE)
trControl <- trainControl(method="LGOCV",
index=cv,
savePredictions=TRUE,
preProcOptions=NULL)
plsGrid <- data.frame(ncomp=1:8)
m.pls <- train(x, y,
method="pls",
trControl=trControl,
tuneGrid=plsGrid)
knnGrid <- data.frame(k=seq(1,20,2))
m.knn <- train(x, y,
method="knn",
trControl=trControl,
tuneGrid=knnGrid)
rfGrid <- data.frame(mtry=c(100,500,1000))
m.rf <- train(x, y,
method="rf",
trControl=trControl,
tuneGrid=rfGrid,
importance=TRUE)
m.gbm <- train(x, y,
method="gbm",
trControl=trControl)
svmGrid <- expand.grid(C=c(10, 20, 50, 100), sigma=10^-(3:7))
m.svm <- train(x, y,
method="svmRadial",
trControl=trControl,
tuneGrid=svmGrid)
save(m.gbm, file="data/sol1_gbm.RData")
save(m.knn, file="data/sol1_knn.RData")
save(m.pls, file="data/sol1_pls.RData")
save(m.rf, file="data/sol1_rf.RData")
save(m.svm, file="data/sol1_svm.RData")
Загрузим построенные модели и значение свойства
m.gbm <- local.load("data/sol1_gbm.RData")
m.knn <- local.load("data/sol1_knn.RData")
m.pls <- local.load("data/sol1_pls.RData")
m.rf <- local.load("data/sol1_rf.RData")
m.svm <- local.load("data/sol1_svm.RData")
y <- local.load("data/sol_y1.RData")
Теперь необходимо из каждой модели извлечь результаты прогноза для кросс-валидации для лучшей модели, чтобы затем усреднить их и сопоставить с наблюдаемыми значениями.
Создадим функцию, которая будет извлекать результаты прогноза для кросс-валидации из модели, полученной пакетом caret
getCV <- function(caret.model) {
df <- caret.model$pred
best <- caret.model$bestTune
ids <- apply(df[ ,names(best), drop=FALSE], 1, function(r) all(r == best[1,]) )
df <- df[ids, ]
df <- df[order(df$rowIndex), c("pred")]
return(df)
}
Создадим еще одну вспомогательную функцию, которая при создании списка из объектов, названия обектов будет присваивать элементам списка.
named.list <- function(...) {
names <- as.list(substitute(list(...)))[-1]
setNames(list(...), names)
}
Объединим все модели в виде списка и извлечем из них значения прогноза для кросс-валидации
models.list <- named.list(m.gbm, m.knn, m.pls, m.rf, m.svm)
cv.pred <- as.data.frame(sapply(models.list, getCV))
head(cv.pred)
m.gbm m.knn m.pls m.rf m.svm
1 -2.799118 -3.626667 -2.148206 -3.000497 -2.677521
2 -2.683472 -3.553333 -2.214603 -2.639931 -2.540740
3 -2.589948 -3.603333 -2.378869 -3.280679 -3.100913
4 -2.639194 -3.853333 -2.297653 -3.249204 -3.106995
5 -3.001957 -2.823333 -1.795756 -3.036543 -2.766391
6 -3.231845 -3.556667 -2.535705 -3.353757 -3.070914
Добавим колонку со средним значением прогноза по всем моделям
cv.pred$mean <- rowMeans(cv.pred)
head(cv.pred)
m.gbm m.knn m.pls m.rf m.svm mean
1 -2.799118 -3.626667 -2.148206 -3.000497 -2.677521 -2.850402
2 -2.683472 -3.553333 -2.214603 -2.639931 -2.540740 -2.726416
3 -2.589948 -3.603333 -2.378869 -3.280679 -3.100913 -2.990748
4 -2.639194 -3.853333 -2.297653 -3.249204 -3.106995 -3.029276
5 -3.001957 -2.823333 -1.795756 -3.036543 -2.766391 -2.684796
6 -3.231845 -3.556667 -2.535705 -3.353757 -3.070914 -3.149778
Теперь посчитаем величину \( RMSE \) и \( R^2 \) для каждой модели и усредненного прогноза. Для этого создадим функцию, которая будет считать \( R^2 \) как \( 1-\frac{PRESS}{SS} \)
R2test <- function(ts.pred, ts.obs, ws.obs.mean = mean(ts.obs)) {
return(1 - sum((ts.pred - ts.obs)**2)/sum((ts.obs - ws.obs.mean)**2))
}
Применим функцию к каждой колонке и округлим получившийся результат до 3 знаков
cv.rmse <- sapply(cv.pred, Metrics::rmse, y)
cv.r2 <- sapply(cv.pred, cor, y) ^ 2
cv.r2test <- sapply(cv.pred, R2test, y)
cv.rmse <- round(cv.rmse, 3)
cv.r2 <- round(cv.r2, 3)
cv.r2test <- round(cv.r2test, 3)
Объединим результаты в таблицу
cv.res <- as.data.frame(rbind(cv.rmse, cv.r2, cv.r2test))
cv.res
m.gbm m.knn m.pls m.rf m.svm mean
cv.rmse 0.697 1.140 0.866 0.719 0.738 0.680
cv.r2 0.896 0.734 0.841 0.893 0.884 0.907
cv.r2test 0.896 0.722 0.839 0.889 0.883 0.901
Видно, что модель knn заметно уступает в качестве прогноза другим моделям.
Стекинг это способ объединения прогнозов разных моделей используя линейную модель, в которой коээфициенты играют роль весов, определеяющих вклад каждой из объединяемых моделей. Поскольку результаты прогноза различных моделей сильно коррелируют между собой, то для их объединения в виде линейной модели используется метод неотрицательных наименьших квадратов (non-negative least squares), т.е. все коэффициенты в уравнении больше либо равны нулю и по сути отражают вклад каждой модели в консенсусный прогноз.
Для построения таких моделей можно воспользоваться пакетом nnls. Загрузим его и построим модель. В качестве параметров в модель передается матрица значений предсказанных разными моделями и вектор наблюдаемых значений.
require(nnls)
m.nnls <- nnls(as.matrix(cv.pred[, -ncol(cv.pred)]), y)
Посчитаем статистические характеристики и добавим их к ранее созданной таблице результатов.
cv.res$nnls <- round(c(Metrics::rmse(m.nnls$fitted, y),
cor(m.nnls$fitted, y) ^ 2,
R2test(m.nnls$fitted, y)), 3)
cv.res
m.gbm m.knn m.pls m.rf m.svm mean nnls
cv.rmse 0.697 1.140 0.866 0.719 0.738 0.680 0.622
cv.r2 0.896 0.734 0.841 0.893 0.884 0.907 0.919
cv.r2test 0.896 0.722 0.839 0.889 0.883 0.901 0.917
В результате видно, что стекинг незначительно улучшает результаты прогноза, если ориентироваться на результаты полученные на кросс-валидации.
Посмотрим коеффициенты в полученной модели
coef(m.nnls)
[1] 0.38438567 0.00000000 0.05772164 0.20288088 0.37905466
Видно, что модель knn не включена в итоговую модель, а модель pls вносит очень незначительный вклад. svm и gbm модели имеют практически одинаковый вклад, превосходящий вклад rf модели. Хотя модель rf сама по себе имеет чуть лучшую прогностическую способность по сравнению с svm.
Загрузим данные для второй выборки по растворимости.
x.test <- local.load("data/sol_x2.RData")
y.test <- local.load("data/sol_y2.RData")
Поскольку это разные датасеты, то сперва надо добавить отсутствующие дескрипторы, удалить лишние дескрипторы и нормировать их. Поскольку мы не сохранили объект preProcess для обучающей выборки, загрузим первый набор дескрипторов и произведем все необходимые манипуляции.
x <- local.load("data/sol_x1.RData")
preProc <- preProcess(x, method=c("scale", "center"))
x <- predict(preProc, x)
x.test[setdiff(names(preProc$mean), colnames(x.test))] <- 0
x.test <- x.test[,names(preProc$mean)]
x.test <- predict(preProc, x.test)
Предскажем для новой выборки значения растворимости всеми моделями из списка models.list
test.pred <- as.data.frame(sapply(models.list, predict, x.test))
head(test.pred)
m.gbm m.knn m.pls m.rf m.svm
1 -1.1778582 -0.440000 -0.3168113 -0.8795913 -0.6089708
2 -1.0579281 -0.440000 -0.2862568 -1.0169123 -1.0142728
3 -0.7143699 -0.440000 -0.2963245 -1.1513807 -1.7784464
4 -0.8748695 -0.440000 -0.0164618 -1.1532090 -2.2388806
5 -0.9348645 -1.346667 -1.4939263 -1.1176407 -0.8709753
6 -1.7129724 -1.473333 -1.3859660 -1.2913827 -1.3519756
Вычислим среднее значение по всем моделям
test.pred$mean <- rowMeans(test.pred)
head(test.pred)
m.gbm m.knn m.pls m.rf m.svm mean
1 -1.1778582 -0.440000 -0.3168113 -0.8795913 -0.6089708 -0.6846463
2 -1.0579281 -0.440000 -0.2862568 -1.0169123 -1.0142728 -0.7630740
3 -0.7143699 -0.440000 -0.2963245 -1.1513807 -1.7784464 -0.8761043
4 -0.8748695 -0.440000 -0.0164618 -1.1532090 -2.2388806 -0.9446842
5 -0.9348645 -1.346667 -1.4939263 -1.1176407 -0.8709753 -1.1528147
6 -1.7129724 -1.473333 -1.3859660 -1.2913827 -1.3519756 -1.4431260
Рассчитаем статистические характеристики прогноза для внешней тестовой выборки
test.rmse <- sapply(test.pred, Metrics::rmse, y.test)
test.r2 <- sapply(test.pred, cor, y.test) ^ 2
test.r2test <- sapply(test.pred, R2test, y.test, mean(y))
test.rmse <- round(test.rmse, 3)
test.r2 <- round(test.r2, 3)
test.r2test <- round(test.r2test, 3)
Объединим результаты в таблицу
test.res <- as.data.frame(rbind(test.rmse, test.r2, test.r2test))
test.res
m.gbm m.knn m.pls m.rf m.svm mean
test.rmse 0.844 1.342 1.234 0.828 0.969 0.722
test.r2 0.685 0.530 0.617 0.702 0.611 0.779
test.r2test 0.709 0.264 0.377 0.720 0.616 0.787
Предскажем значение растворимости в соответствии с полученной ранее моделью для стекинга. Поскольку стандартная функция predict не реализована для таких моделей, то прогноз выполним самостоятельно.
stack.pred <- apply(test.pred[, -ncol(test.pred)], 1, function(i) sum(i * coef(m.nnls)) )
Расчитаем статистику и добавим ее к таблице результатов
test.res$nnls <- round(c(Metrics::rmse(stack.pred, y.test),
cor(stack.pred, y.test) ^ 2,
R2test(stack.pred, y.test, mean(y))),
3)
Сравним полученные результаты
cv.res
m.gbm m.knn m.pls m.rf m.svm mean nnls
cv.rmse 0.697 1.140 0.866 0.719 0.738 0.680 0.622
cv.r2 0.896 0.734 0.841 0.893 0.884 0.907 0.919
cv.r2test 0.896 0.722 0.839 0.889 0.883 0.901 0.917
test.res
m.gbm m.knn m.pls m.rf m.svm mean nnls
test.rmse 0.844 1.342 1.234 0.828 0.969 0.722 0.753
test.r2 0.685 0.530 0.617 0.702 0.611 0.779 0.748
test.r2test 0.709 0.264 0.377 0.720 0.616 0.787 0.768
Видно, что усреднение моделей, даже включая довольно слабые модели, такие как knn в нашем случае, дает заметное улучшение точности прогноза по сравнению с отдельными моделями. При чем для тестовой выборки это более заметно. Результаты стекинга же в данном случае мало отличаются от результатов простого усреднения моделей.
test.pred и cv.pred этого урока + те результаты которые вы получите в ходе выполнения домашнего задания. Формат текстовых файлов. Первая строка - имена моделей, все последующие строки результаты прогноза соединений этими моделями. Первая колонка - имена соединений.