Автор: Владимир Шитиков
Методы обобщения моделей и прогнозов
По аналогии с коллективными методами принятия решений, столь эффективно используемыми в человеческом обществе, принято считать, что суммарная эффективность любой мультимодельной системы распознавания или прогнозирования теоретически будет в среднем выше отдельных ее членов. Поэтому в последние несколько десятилетий активно разрабатывались возможные подходы к тому, как построить на одних и тех же исходных данных некоторый "коллектив" (ensemble) частных одно- или разнотипных моделей и выполнить их обобщение (averaging) с целью получить более обоснованное комбинированное решение (forecast combinations, или multimodel inference).
В предыдущем сообщении мы с помощью пакета MuMIn выполняли ранжирование построенных моделей по их относительной подтверждающей силе (strength of evidence), которая количественно оценивалась с использованием информационных критериев. Однако часто существует обоснованное сомнение в том, что та или иная модель, отобранная как оптимальная, сможет в полной мере объяснить все механизмы моделируемых процессов и выполнять устойчивые прогнозы на всей области определения исходных данных. Не исключено, что существенную долю полезной информации несут в себе вторая, третья и последующие модели из рассматриваемого набора. Учет этой информации часто позволяет скомпенсировать стохастические флуктуации прогнозов единственной модели, выбранной в качестве оптимальной.
Выделяют два возможных подхода к получению коллективного решения (Burnham & Anderson 2002): объединение прогнозов и усреднение параметров моделей. В первом случае рассматриваются расчетные значения \(\hat{y}_1, \hat{y}_2, \dots, \hat{y}_r\) отклика \(Y\), полученные с помощью \(r\) различных моделей. Тогда коллективным мультимодельным прогнозом \(\hat{y}_g\) будет расчетное значение анализируемой переменной \(Y\), вычисленное как некоторая функция индивидуальных прогнозов \( \hat{y}_g = f(\hat{y}_1, \hat{y}_2, \dots, \hat{y}_r) \). На практике чаще всего применяют простейшие методы взвешенного усреднения, которые представляют решение в виде линейной комбинации наиболее информативного множества частных прогнозов \[\hat{y}_g = \sum_{k=1}^{r} \hat{y}_k \omega_{k},\] где \(\omega_{k}\) – весовые коэффициенты для каждого \(\hat{y}_k\), значения которых неизменны на всем диапазоне варьирования исходных данных \(\boldsymbol{X}\).
Второй путь основан на использовании различных способов объединения оценок параметров \(\hat{\beta}_{k}\) для отобранного подмножества регрессионных моделей \(g_1, g_2, \dots, g_r\). В результате этого рассчитываются коэффициенты новой агрегированной модели регрессии \(\hat{\beta}_{g} = f(\hat{\beta}_{k} | g_k, \omega_{k})\), обеспечивающие необходимую точность прогнозирования. Функция model.avg() из пакета MuMIn использует при этом два простых механизма усреднения моделей (model averaging). В обоих случаях для вычисления средневзвешенных коэффициентов \(\hat{\beta}_{g}\) чаще всего применяются веса \(w_k\) (weight), полученные на основе разностей информационных критериев (см. первое сообщение).
Агрегированная модель с "полными" коэффициентами ("full", или "shrinkage model-averaged coeffcients") основана на формуле простого взвешенного среднего, предполагающей, что веса \(\omega_{k}\) для коэффициентов каждой объединяемой модели одинаковы для всех предикторов:
\[ \hat{\beta}_{ig} = \sum_{k=1}^{r} \hat{\beta}_{ik} \omega_{k}, \qquad \omega_{k} = w_k / \sum _{k=1}^{r} w_k , \]
где \(\hat{\beta}_{ik}\) - коэффициент при \(i\)-й переменной в \(k\)-й модели; \(\hat{\beta}_{ik}\) = 0, если переменная \(i\) не включена в модель \(g_k\). Отметим такую важную особенность, что ошибки коэффициентов такой агрегированной регрессионной модели включают дополнительный член (Lukacs et al. 2009): \[ var(\hat{\beta}_{ig}) = \sum_{k=1}^{r} \omega_{k} [var(\hat{\beta}_{ik}) + (\hat{y}_g - \hat{y}_k)^2 ], \] что обычно трактуется как некая особая форма регуляризации (сходная с гребневой регрессией или лассо), а сами коэффициенты называют сжатыми ("coefficients with shrinkage").
Для агрегированной модели с "условными" коэффициентами ("conditional", или "subset model-averaged coefficients") используется механизм взвешенного усреднения с учетом коэффициентов только тех переменных, которые включены в конкретные частные модели. Для этого составляется матрица индикаторных значений \(\Gamma_{ik}\), принимающих значения 1, если предиктор i входит в k-ю модель, и 0 в противном случае. Тогда коэффициенты агрегированной модели рассчитываются как \[ \hat{\beta}_{ig} = \sum_{k=1}^{r} \hat{\beta}_{ik} \omega_{k} \Gamma_{ik} / \sum_{k=1}^{r} \omega_{k} \Gamma_{ik}.\] Условные коэффициенты по своей величине всегда превышают полные и равны им, если \(i\)-я переменная включена во все \(r\) моделей.
Существуют разные мнения относительно того, какие коэффициенты обобщенной модели ("условные" или "полные") являются более адекватными. Мы полагаем, что эта дискуссия не имеет особого смысла, поскольку оба подхода основаны на двух теоретически корректных способах нахождения средних (т.е. с учетом отсутствующих значений или без). Многое зависит от того, какова структура причинно-следственных связей в конкретном наборе данных: использование "условных" коэффициентов, придающих несколько большее значение относительно "слабым" предикторам, чаще всего оказывается предпочтительным в сложных случаях статистических взаимодействий.
Для оценки весов \(\omega_{k}\) кроме использования разностей информационных критериев \(\Delta\), можно применять и иные методы, иногда более подходящие для решения конкретной задачи. В пакет MuMIn включены функции расчета вектора \(\omega\), основанные на создании повторных выборок (бутстреп, кросс-проверка, метод «складного ножа»), алгоритме Бейтса-Гренджера, косинусных функциях, адаптивной регрессии, методах Монте-Карло и т.д.
Теплота гидратации цемента
В предыдущем сообщении мы рассмотрели набор данных Cement и построили для него 16 возможных частных моделей, которые были ранжированы по возрастанию информационного критерия AICc. Сформируем подмножество наилучших моделей, покрывающих 95% вычисленных весов, т.е. cumsum(weight) <= .95. Дополнительно с использованием параметра extra рассчитаем для каждой модели такие общепринятые показатели, как обычный и скорректированный коэффициенты детерминации, а также F-критерий:
library(MuMIn) options(na.action = "na.fail") # Пример из Burnham and Anderson (2002), page 100: fm1 <- lm(y ~ X1 + X2 + X3 + X4, data = Cement) ms1 <- dredge(fm1, extra = list( "R^2", "*" = function(x) { s <- summary(x) c(adjRsq = s$adj.r.squared, F = s$fstatistic[[1]]) }) ) ms1[1:4,6:13] R^2 *.adjRsq *.F df logLik AICc delta weight 4 0.9786784 0.9744140 229.5037 4 -28.15620 69.31239 0.000000 0.5657106 12 0.9823355 0.9764473 166.8317 5 -26.93314 72.43771 3.125321 0.1185603 8 0.9822847 0.9763796 166.3449 5 -26.95180 72.47503 3.162633 0.1163690 10 0.9724710 0.9669653 176.6270 4 -29.81705 72.63411 3.321714 0.1074715 confset.95p <- get.models(ms1, subset = cumsum(weight) <= .95)
Обратим внимание на неоднозначный выбор оптимальной модели с использованием разных критериев адекватности. Выполним построение агрегированной модели, обобщающей 4 наилучшие частные модели:
avgm <- model.avg(confset.95p) summary(avgm) Component model call: lm(formula = y ~ <4 unique rhs>, data = Cement) Component models: df logLik AICc delta weight 12 4 -28.16 69.31 0.00 0.62 124 5 -26.93 72.44 3.13 0.13 123 5 -26.95 72.48 3.16 0.13 14 4 -29.82 72.63 3.32 0.12 Term codes: X1 X2 X3 X4 1 2 3 4 Model-averaged coefficients: Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 60.4843 17.9260 18.2147 3.321 0.000898 *** X1 1.4920 0.1575 0.1747 8.542 < 2e-16 *** X2 0.6250 0.1203 0.1292 4.839 1.31e-06 *** X4 -0.4160 0.2289 0.2408 1.728 0.084009 . X3 0.2500 0.1847 0.2132 1.173 0.240898 Full model-averaged coefficients (with shrinkage): Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 60.48430 17.92604 18.21471 3.321 0.000898 *** X1 1.49198 0.15745 0.17467 8.542 < 2e-16 *** X2 0.55106 0.23133 0.23552 2.340 0.019299 * X4 -0.10354 0.21306 0.21628 0.479 0.632129 X3 0.03204 0.10656 0.11317 0.283 0.777105 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Relative variable importance: X1 X2 X4 X3 Importance: 1.00 0.88 0.25 0.13 N containing models: 4 3 2 1
Как видим, значения условного (subset) коэффициента для переменной X1, включенной во все частные модели, не отличается от полного (full with shrinkage), тогда как для переменных X3 и X4 разности коэффициентов весьма значительны.
Поскольку ошибки условных коэффициентов можно найти как взвешенные суммы ошибок коэффициентов частных моделей \( \hat{\beta}_{ig} = \sum_{k=1}^{r} \omega_{k} var(\hat{\beta}_{ik}) \), то легко вычислить соответствующие доверительные интервалы для задаваемых доверительных вероятностей:
CI <- cbind(confint(avgm), confint(avgm, level=0.85)) CI[ , c(1,3,4,2)] 2.5 % 7.5 % 92.5 % 97.5 % (Intercept) 24.7841330 34.43403292 86.53457235 96.18447223 X1 1.1496363 1.25038159 1.73357355 1.83431889 X2 0.3718565 0.44425065 0.80580132 0.87819546 X4 -0.8878522 -0.75563682 -0.07634164 0.05587375 X3 -0.1678276 -0.04066862 0.54070384 0.66786280
Как было показано выше, при расчете ошибок полных коэффициентов необходимо учесть дополнительный "штрафной" член (Lukacs et al., 2009). После этого доверительные интервалы можно вычислить обычным способом как \( \hat{\beta_{ig}} \mp t_{f, \alpha / 2} [var(\hat{\beta_{ig}})]^{0.5} \), где число степеней свободы \(f = 35\) можно считать вполне допустимым для практических расчетов.
# объект model.averaging включает оценки "сжатых" коэффициентов shrinkage.coef <- avgm$coef.shrinkage # и значения бета-коэффициентов для всех частных моделей coef.array <- avgm$coefArray coef.array <- replace(coef.array, is.na(coef.array), 0) # создадим пустую таблицу для хранения результатов расчетов shrinkage.estimates <- data.frame(shrinkage.coef,variance=NA) # вычислим shrinkage-скорректированные ошибки for(i in 1:dim(coef.array)[3]){ input <- data.frame(coef.array[,,i], weight = Weights(avgm)) variance <- rep(NA, dim(input)[1]) for (j in 1:dim(input)[1]){ variance[j] <- input$weight[j] * (input$Std..Err[j]^2 + (input$Estimate[j] - shrinkage.estimates$shrinkage.coef[i])^2) } shrinkage.estimates$variance[i] <- sum(variance) } # вычислим доверительные интервалы shrinkage.estimates$lci <- shrinkage.estimates$shrinkage.coef - 2.03*sqrt(shrinkage.estimates$variance) shrinkage.estimates$uci <- shrinkage.estimates$shrinkage.coef + 2.03*sqrt(shrinkage.estimates$variance) print(shrinkage.estimates) shrinkage.coef variance lci uci (Intercept) 60.48430263 321.34280013 24.09444766 96.8741576 X1 1.49197757 0.02479082 1.17235203 1.8116031 X2 0.55105655 0.05351159 0.08146534 1.0206478 X4 -0.10354105 0.04539387 -0.53604956 0.3289675 X3 0.03203825 0.01135571 -0.18428500 0.2483615
Заметим, что как регуляризованные, так и условные коэффициенты при Х3 и Х4 оказались статистически незначимыми, поскольку их доверительные интервалы включают значение 0.
Совокупность расчетных значений \(\hat{y_g}\) можно оценить обычным образом с помощью функции predict(), которая в данном случае имеет специальный параметр full, принимающий значения TRUE или FALSE в зависимости от того, по какому набору коэффициентов следует вести расчет. Оценим точность ансамблей из усредненных моделей, основанных на полных и условных коэффициентах, по трем показателям: среднему абсолютному отклонению (MAE), корню из среднеквадратичного отклонения (RSME) и квадрату коэффициента детерминации Rsq = 1 - NSME, где NSME - относительная ошибка, равная отношению средних квадратов отклонений от регрессии и от общего среднего:
# Функция, выводящая вектор критериев ModCrit <- function (pred, fact) { mae <- mean(abs(pred-fact)) rmse <- sqrt(mean((pred-fact)^2)) Rsq <- 1-sum((fact-pred)^2)/sum((mean(fact)-fact)^2) c(MAE=mae, RSME=rmse, Rsq = Rsq ) } rbind( averaged.subset = ModCrit(predict(avgm, full = FALSE),Cement$y), averaged.full = ModCrit(predict(avgm, full = TRUE),Cement$y)) MAE RSME Rsq averaged.subset 5.906004 7.120158 0.7573218 averaged.full 1.709184 1.958280 0.9816430
Чтобы понять смысл полученных расхождений в оценках качества моделей, отобразим имеющиеся статистические закономерности на графике. Для этого составим новую таблицу, обрабатываемую функцией predict(), в которой значения X1 "пробегают" по всему диапазону исходных данных, а остальные три переменные зафиксированы на уровне значений их средних. Покажем на графике прямые, связанные с каждой из четырех частных моделей, и две прямые, соответствующие усредненным моделям (полным - "full averaged", и условным - "subset averaged").
nseq <- function(x, len = length(x)) seq(min(x, na.rm = TRUE), max(x, na.rm=TRUE), length = len) newdata <- as.data.frame(lapply(lapply(Cement[, -1], mean), rep, 25)) newdata$X1 <- nseq(Cement$X1, nrow(newdata)) n <- length(confset.95p) pred <- data.frame( model = sapply(confset.95p, predict, newdata = newdata), averaged.subset = predict(avgm, newdata, full = FALSE), averaged.full = predict(avgm, newdata, full = TRUE) ) matplot(newdata$X1, pred, type = "l", lwd = c(rep(2,n),3,3), lty = 1, xlab = "X1", ylab = "y", col=1:7)
Поскольку функция predict() с параметром se.fit = TRUE возвращает стандартную ошибку регрессии, то для полной модели построим также линии доверительных огибающих:
pred.se <- predict(avgm, newdata, se.fit = TRUE) y <- pred.se$fit ci <- pred.se$se.fit * 2 matplot(newdata$X1, cbind(y, y - ci, y + ci), add = TRUE, type="l", lty = 2, col=c(1,1), lwd = 1) legend("topleft", legend=c(lapply(confset.95p, formula), paste(c("subset", "full"), "averaged"), "averaged predictions + CI"), lty = c(rep(1,6),2), lwd = c(rep(2,n),3,3,1), cex = .75, col=c(1:6,1))
На графике видно, что модель с условными коэффициентами склонна к очевидному смещению в сторону уменьшения значений теплоты гидратации Y.
Вероятность заболеть диабетом
Используем в качестве второго примера набор данных pima, обсуждаемый в известной книге Дж. Фаравея (Faraway, 2002) и входящий в пакет faraway. Взрослых женщин-индианок из племени пима протестировали на наличие признаков диабета (test = 0 если заболевание отсутствует, и 1 если диагноз положительный). Одновременно учитывали следующие показатели: число беременностей (pregnant), концентрация глюкозы в плазме (glucose), диастолическое кровяное давление (diastolic), толщина кожно-жировой складки трицепса (triceps), инсулин в сыворотке крови (insulin), индекс массы тела (bmi), индекс родственно-генетической предрасположенности к диабету (diabetes) и возраст пациента (age).
Следуя рекомендациям Дж. Фаравея, обратим внимание на некоторые недопустимые "небрежности" при формировании набора данных. Ряд медицинских показателей равны 0, что, вероятно, означает пропуск данных, поскольку трудно себе представить живого пациента с нулевым кровяным давлением. Удалим подобные строки из набора данных. Вызывает также сомнение максимум рожденных детей (pregnant = 17), но сочтем такое значение показателя реально возможным.
library(faraway) data(pima) str(pima) pima$diastolic [pima$diastolic == 0] <-NA pima$glucose [pima$glucose == 0] <-NA pima$triceps [pima$triceps == 0] <-NA pima$insulin [pima$insulin == 0] <-NA pima$bmi [pima$bmi == 0] <-NA pima <- pima[complete.cases(pima),] summary(pima$test <- as.factor(pima$test)) 0 1 262 130
Построим "глобальную" модель логистической регрессии, включающую все наблюдаемые переменные:
fmlr <- glm(test ~ .,data = pima, family = binomial()) summary(fmlr) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.004e+01 1.218e+00 -8.246 < 2e-16 *** pregnant 8.216e-02 5.543e-02 1.482 0.13825 glucose 3.827e-02 5.768e-03 6.635 3.24e-11 *** diastolic -1.420e-03 1.183e-02 -0.120 0.90446 triceps 1.122e-02 1.708e-02 0.657 0.51128 insulin -8.253e-04 1.306e-03 -0.632 0.52757 bmi 7.054e-02 2.734e-02 2.580 0.00989 ** diabetes 1.141e+00 4.274e-01 2.669 0.00760 ** age 3.395e-02 1.838e-02 1.847 0.06474 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Null deviance: 498.10 on 391 degrees of freedom Residual deviance: 344.02 on 383 degrees of freedom AIC: 362.02
Создадим набор из 256 всех возможных частных моделей и ограничим его условием cumsum(weight) < 0.85. После этого создадим усредненную модель, объединяющую 17 отобранных частных моделей:
allmod <- dredge(fmlr) avgmod <-model.avg(allmod,cumsum(weight) < 0.85) summary(avgmod) Component model call: glm(formula = test ~ <17 unique rhs>, family = binomial(), data = pima) Component models: df logLik AICc delta weight 12357 6 -172.44 357.10 0.00 0.17 1235 5 -173.62 357.39 0.29 0.14 123578 7 -172.21 358.72 1.61 0.07 123567 7 -172.23 358.76 1.65 0.07 2357 5 -174.36 358.88 1.77 0.07 12356 6 -173.35 358.92 1.82 0.07 12358 6 -173.37 358.96 1.85 0.07 123457 7 -172.44 359.17 2.07 0.06 12345 6 -173.62 359.45 2.35 0.05 23578 6 -173.93 360.09 2.98 0.04 1235678 8 -172.02 360.41 3.31 0.03 123568 7 -173.12 360.54 3.43 0.03 23567 6 -174.21 360.65 3.54 0.03 1234578 8 -172.21 360.79 3.69 0.03 1234567 8 -172.23 360.83 3.72 0.03 23457 6 -174.31 360.85 3.74 0.03 123456 7 -173.35 360.99 3.89 0.02 Term codes: age bmi diabetes diastolic glucose insulin pregnant triceps 1 2 3 4 5 6 7 8 Model-averaged coefficients: Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) -9.9727948 1.1355744 1.1389894 8.756 < 2e-16 *** age 0.0426278 0.0185999 0.0186439 2.286 0.02223 * bmi 0.0744313 0.0230849 0.0231533 3.215 0.00131 ** diabetes 1.1378544 0.4259816 0.4273032 2.663 0.00775 ** glucose 0.0371798 0.0052995 0.0053156 6.994 < 2e-16 *** pregnant 0.1018981 0.0606241 0.0607644 1.677 0.09355 . triceps 0.0122192 0.0170973 0.0171509 0.712 0.47619 insulin -0.0008761 0.0013114 0.0013155 0.666 0.50543 diastolic -0.0002998 0.0118512 0.0118880 0.025 0.97988 Full model-averaged coefficients (with shrinkage): Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) -9.973e+00 1.136e+00 1.139e+00 8.756 < 2e-16 *** age 3.580e-02 2.313e-02 2.316e-02 1.546 0.12214 bmi 7.443e-02 2.308e-02 2.315e-02 3.215 0.00131 ** diabetes 1.138e+00 4.260e-01 4.273e-01 2.663 0.00775 ** glucose 3.718e-02 5.300e-03 5.316e-03 6.994 < 2e-16 *** pregnant 6.291e-02 6.872e-02 6.879e-02 0.915 0.36045 triceps 3.250e-03 1.034e-02 1.036e-02 0.314 0.75381 insulin -2.452e-04 7.975e-04 7.994e-04 0.307 0.75907 diastolic -6.369e-05 5.464e-03 5.481e-03 0.012 0.99073 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Relative variable importance: bmi diabetes glucose age pregnant insulin triceps diastolic Importance: 1.00 1.00 1.00 0.84 0.62 0.28 0.27 0.21 N containing models: 17 17 17 13 11 7 6 6
Оценим ошибки предсказания модели с использованием обоих наборов коэффициентов ("full" и "subset"). Предварительно определим функцию logit2prob() для преобразования значений логита в соответствующие вероятности:
logit2prob <- function(logit){ odds <- exp(logit); prob <- odds / (1 + odds) return(prob) } Pred <- ifelse(logit2prob(predict(avgmod, full=TRUE))>0.5,1,0) table(Факт=pima$test, Прогноз=Pred) Acc <- mean(pima$test == Pred) paste("Точность=", round(100*Acc, 2), "%", sep="") # Для модели с «полными» коэффициентами Прогноз Факт 0 1 0 233 29 1 55 75 [1] "Точность=78.57%" Pred <- ifelse(logit2prob(predict(avgmod, full=FALSE))>0.5,1,0) table(Факт=pima$test, Прогноз=Pred) Acc <- mean(pima$test == Pred) paste("Точность=", round(100*Acc, 2), "%", sep="") # Для модели с «условными» коэффициентами Прогноз Факт 0 1 0 218 44 1 37 93 [1] "Точность=79.34%"
Формально модель с условными коэффициентами оказалась в целом чуть точнее, хотя у нее наблюдается отчетливая тенденция к гипердиагностике, тогда как для модели с полными коэффициентами характерен более частый пропуск анализируемого признака.
Мы выполнили обобщение 17 частных моделей. Возникают правомерные вопросы: численность какого коллектива является оптимальной и не удастся ли получить более эффективный ансамбль из меньшего числа моделей?
Чтобы сделать наш анализ независимым, воспользуемся внешним критерием - точностью предсказания на отдельном проверочном подмножестве, не участвовавшем в построении моделей. Проведем разделение имеющейся выборки на обучающую и проверочную. Использование предсказываемой качественной переменной в качестве аргумента функции createDataPartition() обеспечивает необходимый баланс уровней этой переменной при разделении исходных данных. Для обучения моделей выделим 75% исходной выборки (295 пациентов), а для независимого тестирования оставим 25% (97 пациентов):
library(caret) set.seed(7) train <- unlist(createDataPartition(pima$test, p=0.75)) c(nrow(pima[train,]), nrow(pima[-train,])) [1] 295 97
fmlrt <- glm(test ~ .,data = pima[train,], family = binomial()) allmodt <- dredge(fmlrt) mod20 <- get.models(allmod, subset = delta <= 5) # Последовательность значений delta model.avg(mod20)$msTable$delta [1] 0.0000000 0.2871189 1.6125558 1.6531455 1.7733025 1.8214838 1.8543884 [8] 2.0688914 2.3480942 2.9831460 3.3082866 3.4323749 3.5422478 3.6909647 [15] 3.7247098 3.7420870 3.8874319 3.9254308 4.7904076 4.9851575
Составим таблицу, включающую сравнительные показатели эффективности коллективов моделей с численностью MCount от 2 до 20. Столбцы таблицы - точность прогнозирования переменной test различными моделями (с условными и полными коэффициентами) на тестовой выборке (subsetCP и fullCP) и обучающей последовательности (subsetCT и fullCT соответственно). В первой строке таблицы (MCount = 1) - логистическая регрессия с оптимальным составом предикторов.
predtab <- pima[-train,-9] traintab <- pima[train,-9] lmbest <- glm(test ~ bmi+diabetes+glucose+pregnant, data = pima[train,], family = binomial()) summary(lmbest) Ptrain <- ifelse(logit2prob(predict(lmbest))>0.5,1,0) Pred <- ifelse(logit2prob(predict(lmbest,predtab))>0.5,1,0) Res <- data.frame(MCount=1, subsetCP = mean(pima[-train,9]==Pred), fullCP = mean(pima[-train,9]==Pred), subsetCT = mean(pima[train,9]==Ptrain), fullCT = mean(pima[train,9]==Ptrain)) for(k in 2:length(mod20)) { avgmodt <- model.avg(mod20[1:k]) Pred.subset <- ifelse(logit2prob(predict(avgmodt,predtab, full=FALSE))>0.5,1,0) Pred.full <- ifelse(logit2prob(predict(avgmodt,predtab, full=TRUE))>0.5,1,0) Ptrain.subset <- ifelse(logit2prob(predict(avgmodt, traintab, full=FALSE))>0.5,1,0) Ptrain.full <- ifelse(logit2prob(predict(avgmodt, traintab, full=TRUE))>0.5,1,0) Res <- rbind(Res, c(k, mean(pima[-train,9]==Pred.subset), mean(pima[-train,9]==Pred.full), mean(pima[train,9]==Ptrain.subset), mean(pima[train,9]==Ptrain.full))) }
Представим данные из полученной таблицы на графике:
plot ( Res$MCount,Res$subsetCP, type="b", lwd=2, pch =17, col=4, xlab="Число моделей", ylab="Ошибка предсказания", ylim=c(0.76, 0.81)) lines ( Res$MCount,Res$fullCP, type="b", lwd=2, pch =16, col=2) lines ( Res$MCount,Res$subsetCT, type="b", lwd=2, pch =2, col=3) lines ( Res$MCount,Res$fullCT, type="b", lwd=2, pch =1, col=5) legend("bottomright", c("Test Subset", "Test Full","Train Subset", "Train Full"), col= c(4,2,3,5), pch = c(17,16,2,1))
На основе проведенных расчетов можно сделать ряд выводов:
- Эффективность коллектива моделей, как правило, выше одной оптимальной модели, поскольку комбинированное решение может быть эффективным инструментом стабилизации дисперсии прогнозных значений. Такое решение значительно более устойчиво к не вполне объяснимым "провалам", которые свойственны отдельным индивидуальным моделям. Кроме того коллектив может дополнительно учитывать часто весьма полезный вклад "слабых" предикторов, не включенных в оптимальную модель.
- Сам по себе размер ансамбля не должен быть чрезмерно большим. Согласно известной аксиоме (Розенберг Г. (1987) "Тройка, семерка, туз...". Знание-Сила, №1, с. 97), этот оптимум обычно находится в диапазоне от трех до семи членов (предпочтительно при delta < 2).
- Необходима разработка модуля для тестирования ансамблей моделей и нахождения оптимальных гиперпараметров агрегированной модели (механизма усреднения коэффициентов и числа членов r), который был бы подобен функции train() из пакета caret.
Сравнение с другими типами моделей
Используем набор данных pima для построения моделей бинарной классификации другими различными методами, подробно рассмотренными нами ранее (Шитиков, Мастицкий 2017, глава 6). Оптимизацию гиперпараметров и подгонку коэффициентов будем выполнять с применением функции train() из пакета caret. В качестве сравниваемых показателей эффективности применим точность прогнозирования исхода test различными моделями на тестовой выборке (Acc.test) и обучающей выборке (Acc.train). Первые две строки итоговой таблицы заполним значениями из ранее сформированной таблицы Res:
MethodRes <- data.frame(Acc.test = Res$subsetCP[1],Acc.train = Res$subsetCT[1]) MethodRes <- rbind(MethodRes, c(max(c(Res$subsetCP[-1], Res$fullCP[-1])), max(c(Res$subsetCT[-1], Res$fullCT[-1]))) )
# определение схемы тестирования control <- trainControl(method="repeatedcv", number=10, repeats=3, classProbs = TRUE)
# LDA - линейный дискриминантный анализ set.seed(7) fit.lda <- train(pima[train, 1:8], pima$test[train], method="lda", trControl=control) MethodRes <- rbind(MethodRes, c(mean(pima[-train, 9]==predict(fit.lda, predtab)), mean(pima[train, 9]==predict(fit.lda, traintab)))) # SVML - метод опорных векторов с линейным ядром set.seed(7) fit.svL <- train(pima[train, 1:8], pima$test[train], method="svmLinear", trControl=control) MethodRes <- rbind(MethodRes, c(mean(pima[-train, 9]==predict(fit.svL, predtab)), mean(pima[train, 9]==predict(fit.svL, traintab)))) # SVMR - метод опорных векторов с радиальным ядром set.seed(7) fit.svR <- train(pima[train, 1:8], pima$test[train], method="svmRadial", trControl=control, tuneGrid = expand.grid(sigma = 0.4, C = 2)) MethodRes <- rbind(MethodRes, c(mean(pima[-train, 9]==predict(fit.svR, predtab)), mean(pima[train, 9]==predict(fit.svR, traintab)))) # CART - дерево классификации set.seed(7) fit.cart <- train(pima[train, 1:8], pima$test[train], method="rpart", trControl=control) MethodRes <- rbind(MethodRes, c(mean(pima[-train, 9]==predict(fit.cart, predtab)), mean(pima[train, 9]==predict(fit.cart, traintab)))) # RF - случайный лес set.seed(7) fit.rf <- train(pima[train, 1:8], pima$test[train], method="rf", trControl=control) MethodRes <- rbind(MethodRes, c(mean(pima[-train, 9]==predict(fit.rf, predtab)), mean(pima[train, 9]==predict(fit.rf, traintab)))) rownames(MethodRes) <- c( "Логистическая регрессия","Усредненная модель ALR", "Дискриминантный анализ LDA","Опорные векторы SVM (лин.ядро)", "Опорные векторы SVM (рад.ядро)","CART - дерево классификации", "RF - случайный лес") MethodRes Acc.test Acc.train Логистическая регрессия 0.7628866 0.7830508 Усредненная модель ALR 0.7938144 0.8033898 Дискриминантный анализ LDA 0.7938144 0.7898305 Опорные векторы SVM (лин.ядро) 0.7731959 0.7898305 Опорные векторы SVM (рад.ядро) 0.7319588 0.9525424 CART - дерево классификации 0.7731959 0.8033898 RF - случайный лес 0.8041237 1.0000000
Снова (и это становится уже очевидной тенденцией) вне конкуренции оказалась модель случайного леса (Random Forrest), также использующая коллектив моделей, основанных на деревьях классификации. Интересны также результаты, показанные моделью опорных векторов с радиальным ядром: найденная этой моделью гиперплоскость почти идеально разделила объекты обучающей выборки, но это привело к многим ошибкам на тестовой выборке (т.е. произошло типичное "переобучение" модели). Обсуждаемая нами агрегированная модель на основе 6 частных моделей логистической регрессии ALR показала вполне достойные результаты на фоне остальных 6 моделей, построенных другими современными методами.
Рекомендуемая литература
- Burnham KP, Anderson DR (2002) Model selection and multimodel inference: a practical information-theoretic approach, 2nd ed. Springer, New York
- Faraway JJ (2002) Practical Regression and Anova using R. Электронная книга, адрес доступа: www.stat.lsa.umich.edu/~faraway/book
- Lukacs PM, Burnham KP, Anderson DR (2009) Model selection bias and Freedman's paradox. Annals of the Institute of Statistical Mathematics 62(1): 117-125
- Шитиков ВК, Мастицкий СЭ (2017) Классификация, регрессия и другие алгоритмы Data Mining с использованием R. Электронная книга, адрес доступа: https://github.com/ranalytics/data-mining
Отправить комментарий