Введение в R: часть 5

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 в нашем случае, дает заметное улучшение точности прогноза по сравнению с отдельными моделями. При чем для тестовой выборки это более заметно. Результаты стекинга же в данном случае мало отличаются от результатов простого усреднения моделей.

Домашнее задание

  1. Получить консенсусный прогноз на пяти описанных в занятии моделях, используя обычный метод наименьших квадратов (lm), по аналогии с nnls.
  2. Оценить возможности консенсусных подходов (усреднее, стекинг с использованием nnls и lm), используя три и четыре лучшие индивидуальные модели.
  3. Сохранить два текстовых файла с результатами кросс-валидации и результатами прогноза внешней тестовой выборки всеми построенными моделями, т.е. по сути это объекты test.pred и cv.pred этого урока + те результаты которые вы получите в ходе выполнения домашнего задания. Формат текстовых файлов. Первая строка - имена моделей, все последующие строки результаты прогноза соединений этими моделями. Первая колонка - имена соединений.