03 марта 2018

Работа с пакетом MuMIn. Часть 2: Объединение моделей для получения коллективного решения



Автор: Владимир Шитиков

Методы обобщения моделей и прогнозов

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

С использованием примеров из обучающей выборки построим новую «глобальную» модель и сформируем ансамбль из 20 наилучших моделей при delta <= 5:

 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))

     

На основе проведенных расчетов можно сделать ряд выводов:

  1. Эффективность коллектива моделей, как правило, выше одной оптимальной модели, поскольку комбинированное решение может быть эффективным инструментом стабилизации дисперсии прогнозных значений. Такое решение значительно более устойчиво к не вполне объяснимым "провалам", которые свойственны отдельным индивидуальным моделям. Кроме того коллектив может дополнительно учитывать часто весьма полезный вклад "слабых" предикторов, не включенных в оптимальную модель.  
  2. Сам по себе размер ансамбля не должен быть чрезмерно большим. Согласно известной аксиоме (Розенберг Г. (1987) "Тройка, семерка, туз...". Знание-Сила, №1, с. 97), этот  оптимум обычно находится в диапазоне от трех до семи членов  (предпочтительно при delta < 2).
  3. Необходима разработка модуля для тестирования ансамблей моделей и нахождения оптимальных гиперпараметров агрегированной модели (механизма усреднения коэффициентов и числа членов 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


1 комментарий :

kodulehe tegemine комментирует...

Очень информативно!

Отправить комментарий