R: Анализ и визуализация данных

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

2018-03-03T12:44:00.000+00:00

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

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

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

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

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

Работа с пакетом MuMIn. Часть 1: Селекция моделей

2018-01-14T14:02:00.000+00:00

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

Введение

Современные исследования приобретают все более и более обобщающий и стратегический характер, а глубокая стратегия никогда не ограничивается рассмотрением какой-то одной идеи, гипотезы или модели. Принцип "множественности моделей", сформулированный еще в 1890 г. Т. Чемберленом, предполагает формирование набора альтернативных научных гипотез $H_1, H_2, \dots, H_r$, для каждой из которых подбирается адекватная математическая модель. В итоге вместо того, чтобы находить по Фишеру соотношение вероятностей основной и нулевой гипотез $H_0 | H_1$, оценивается относительная сила обоснованности (strength of evidence) каждой из рассматриваемых гипотетических моделей $g_1, g_2, \dots, g_r$. Эта новая парадигма оформилась в современную методологию "Model selection and Multimodel inference" (Burnham, Anderson, 2002), которая базируется на основных принципах теории информации Кульбака-Лейблера (Kullback-Leibler, 1951) и включает ранжирование моделей с последующим формированием статистического вывода на основе этих нескольких моделей.

Пакет MuMIn (от "Multi-Model Inference"), разработанный К. Бартоном, содержит набор функций, которые, используя информационные критерии, реализуют ранжирование и отбор статистических моделей различного типа и их последующее обобщение (model averaging) с целью получения коллективного решения. Ансамбль $G_r$ моделей, включаемых в анализ, может формироваться либо автоматически (например, как все возможные комбинации подмножеств предикторов заданной "глобальной" модели), либо задается исследователем в виде набора конкретных моделей (т.е. математически выраженных гипотез, нуждающихся в проверке). Для всех анализируемых моделей выполняется подгонка их коэффициентов по эмпирическим наборам данных и рассчитываются основные статистики и информационные критерии $IC$ качества аппроксимации. Полный список из $r$ моделей сортируют по уменьшению адекватности на основе заданного $IC$ и устанавливают порог, согласно которому некоторое количество "оптимальных" моделей далее будет использовано для формирования окончательного коллективного решения. Такие коллективные решения могут быть получены с использованием как традиционных взвешенных средних, так и других специальных алгоритмов (метод Бейтса-Гренджера, бутстреп, метод "складного ножа", адаптивная регрессия и др.).

В этой статье мы рассмотрим первую часть описанной процедуры, т.е. построение ранжированного списка моделей на основе информационных критериев с использованием пакета MuMIn.

Информационные критерии

Энтропийная мера Кульбака-Лейблера (K-L) трактуется как доля потерянной информации при использовании модели $g_i$ по сравнению с отображаемой этой моделью "полной реальностью" (full reality). Таким образом, для выбора оптимальной модели из набора $G_r$ необходимо подобрать такой критерий, который соответствует минимуму информационных потерь. В своей ключевой работе Хиротугу Акаике (Akaike, 1973) объединил понятия информационной и статистической теории и пришел к выводу фундаментальной важности: поиск оптимальной модели можно выполнить путем минимизации величины
\[\log(L(\hat{\theta}|data)) - Const,\] где $L(\hat{\theta}|data)$ - максимум функции правдоподобия модели $g_i$ с коэффициентами $\hat{\theta}$, подогнанными по выборочным данным $data$, а $Const$ - некоторая константа, корректирующая смещение и зависящая от числа степеней свободы. Акаике предложил конкретное определение соответствующего информационного критерия: \[ AIC = -2 \log(L(\hat{\theta}|data)) + 2K, \] где K - количество параметров модели.

Далее последовали достаточно многочисленные попытки "подправить" выражение для $Const$. В частности, для малых выборок при $n/K < 40$ рекомендуется использовать скорректированный $AIC$, который имеет вид (Sugiura, 1978): \[ AIC_{c} = AIC + 2K(K + 1) / (n - K - 1). \] Если объем выборки $n \gg K$, то $AIC_{c}$ сводится к $AIC$. В случае обычной регрессии, решаемой методом наименьших квадратов, $\log(L(\hat{\theta}|data)) = -(n/2) \log(RSS/n)$ и \[AIC_{c} = n \log(RSS/n) + 2K + 2K(K + 1) / (n - K - 1), \] где $RSS$ - сумма квадратов остатков подогнанной модели.

Весьма популярен байесовский информационный критерий, который позволяет выбирать модели на основе принципа максимального апостериорного правдоподобия (a posteriori most likely). При аппроксимации формулы Байеса, выполненной Шварцем (Schwarz, 1978), этот критерий имеет вид: \[BIC = -2 \log(L(\hat{\theta}|data)) + K \log(n).\] В литературе (Claeskens, Hjort, 2008) подробно описаны и другие версии информационных критериев: Takeuchi's IC ($TIC$), Deviance IC ($DIC$), Hannan and Quinn criterion ($HQ$), Focused IC ($FIC$), информационные критерии для счетных данных с высокой дисперсией ($QAIC, QAIC_{c}$), $CAICF$ (где $C$ соответствует слову "Consistent", а $F$ означает, что используется информационная матрица Фишера), а также близкие по сути статистики $C_p$ Мэллоуза (Mallows) и $ICOMP$ (Informational Complexity).

Для вычисления перечисленных критериев в пакетах stats, MuMIn, arm и GEE имеются соответствующие функции. Подлежащий использованию тип критерия выбирается из сложно формализуемых соображений и задается в представленных ниже примерах параметром rank = <имя IC-функции>, причем по умолчанию используется $AIC_{c}$.

Рассмотрим процесс построения и селекции моделей на трех примерах с разной степенью их детализации.

Теплота гидратации цемента

Сначала обратимся к хрестоматийному примеру, который основан на наборе данных Cement и подробно комментируется в основной литературе (см. Burnham, Anderson, 2002, стр. 100). Задача заключается в нахождении зависимости теплоты гидратации y (кал/г), выделяющейся при затвердевании различных марок цемента, от химического состава этих марок. Четыре независимых переменных соответствуют содержанию различных оксидных составляющих: X1: $3 CaO \times Al_2O_3$, X2: $3CaO \times SiO_2$, X3: $4CaO \times Al_2O_3 \times Fe_2O_3$, и X4: $2CaO \times SiO_2$. Построим линейную регрессионную модель:

library(MuMIn)
# можно изменить значение по умолчанию "na.omit", чтобы
# не учитывать пропущенные значения
options(na.action = "na.fail") 
# осуществляем построение глобальной модели
fm1 <- lm(y ~ ., data = Cement)
summary(fm1)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  62.4054    70.0710   0.891   0.3991  
X1            1.5511     0.7448   2.083   0.0708 .
X2            0.5102     0.7238   0.705   0.5009  
X3            0.1019     0.7547   0.135   0.8959  
X4           -0.1441     0.7091  -0.203   0.8441  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.446 on 8 degrees of freedom
Multiple R-squared:  0.9824,    Adjusted R-squared:  0.9736 
F-statistic: 111.5 on 4 and 8 DF,  p-value: 4.756e-07

Функция dredge() осуществляет построение всех возможных моделей из различных комбинаций независимых переменных полной линейной модели. В данном примере общее количество таких моделей составляет $r = 2^K = 16.$

 (ms1 <- dredge(fm1))
      Model selection table 
          (Intrc)    X1      X2      X3      X4 df  logLik  AICc delta weight
       4    52.58 1.468  0.6623                  4 -28.156  69.3  0.00  0.566
       12   71.65 1.452  0.4161         -0.2365  5 -26.933  72.4  3.13  0.119
       8    48.19 1.696  0.6569  0.2500          5 -26.952  72.5  3.16  0.116
       10  103.10 1.440                 -0.6140  4 -29.817  72.6  3.32  0.107
       14  111.70 1.052         -0.4100 -0.6428  5 -27.310  73.2  3.88  0.081
       15  203.60       -0.9234 -1.4480 -1.5570  5 -29.734  78.0  8.73  0.007
       16   62.41 1.551  0.5102  0.1019 -0.1441  6 -26.918  79.8 10.52  0.003
       13  131.30               -1.2000 -0.7246  4 -35.372  83.7 14.43  0.000
       7    72.07        0.7313 -1.0080          4 -40.965  94.9 25.62  0.000
       9   117.60                       -0.7382  3 -45.872 100.4 31.10  0.000
       3    57.42        0.7891                  3 -46.035 100.7 31.42  0.000
       11   94.16        0.3109         -0.4569  4 -45.761 104.5 35.21  0.000
       2    81.48 1.869                          3 -48.206 105.1 35.77  0.000
       6    72.35 2.312          0.4945          4 -48.005 109.0 39.70  0.000
       5   110.20               -1.2560          3 -50.980 110.6 41.31  0.000
       1    95.42                                2 -53.168 111.5 42.22  0.000
      Models ranked by AICc(x)

В полученной таблице представлены коэффициенты каждой из построенных моделей, значения логарифмов максимального правдоподобия и скорректированные значения информационного критерия Акаике, по убыванию которых выполнялась сортировка моделей. Важным показателем является величина $\Delta_{i} = AIC_{c_i} - AIC_{c_{min}}$, которая соответствует потере информации Кульбака-Лейблера или, эквивалентно, увеличению расстояния от текущей до оптимальной модели. Шкала $\Delta$ дает возможность исследователю выбрать допустимый порог информационных потерь и ограничить подмножество отбираемых моделей-претендентов по некоторому уровню адекватности.

Применительно к имеющимся данным, относительное правдоподобие $l_i$, оцененное для каждой модели $g_i$, соответствует формализованному уровню ее обоснованности (strength of evidence) по сравнению с наилучшей моделью и легко вычисляется как $l_i = \exp(-0.5 \Delta_{i})$. Это выражение позволяет сформировать хорошо интерпретируемые утверждения, такие как "Модель $g_1$, основанная на X1 и X2, в 192 раза лучше объясняет результаты эксперимента, чем полная модель $g_{16}$" (или "гипотеза $H_4$ в 192 раза более вероятна, чем $H_{16}$"). Нетрудно вычислить, что соотношение обоснованности гипотез для $\Delta = 2$ равно 2.7, а при $\Delta = 4$ оно достигает 7.4. Отметим, что в рамках обсуждаемой методологии не рекомендуется использовать термин "статистическая значимость" в том смысле, в каком он используется в рамках классической дихотомической процедуры Фишера.

Относительная вероятностная мера каждой модели $g_i$ в ансамбле из $G_r$ моделей, или ее вес (weight), может быть вычислена как нормированное значение силы обоснованности \[w_i = \exp(-0.5 \Delta_{i}) / \sum \exp(-0.5 \Delta_{i}) .\] Найденные веса позволяют также оценить относительную важность каждой независимой переменной как сумму весов тех моделей, в которых соответствующий предиктор присутствует:

      # Вычисляем степень важности отдельных предикторов
      importance(ms1)
                           X1   X2   X4   X3  
      Importance:          0.99 0.81 0.32 0.21
      N containing models:    8    8    8    8
      # Выводим на экран график таблицы сгенерированных моделей:
      par(mar = c(3,5,6,4))
      plot(ms1, labAsExpr = TRUE)

Влияние курения на заболеваемость раком

Пакет MuMIn, разумеется, не ограничивается построением моделей обычной линейной регрессии. Список поддерживаемых моделей включает более 30 типов функций из различных пакетов (glm, gam, lmer, polr, MCMCglmm, survreg, zeroinfl, betareg и проч.). Рассмотрим пример селекции моделей логистической регрессии с учетом парных взаимодействий факторов.

Набор данных esoph включает результаты обследования различных групп населения одного из департаментов Франции на предмет заболеваемости раком пищевода. Три порядковых переменных определяли интенсивность приема алкоголя (alcgp) и табакокурения (tobgp) в баллах от 1 до 4, а также градацию возраста обследуемых (agegp) от 1 до 6. Переменные ncases, ncontrols определяли число случаев заболевания раком и объем контрольной группы соответственно.

Хорошим способом познакомиться со структурой данных является отображение ее в виде мозаичной диаграммы:

   ## Реорганизация данных для мозаичной диаграммы
   ttt <- table(esoph$agegp, esoph$alcgp, esoph$tobgp)
   ttt[ttt == 1] <- esoph$ncases
   tt1 <- table(esoph$agegp, esoph$alcgp, esoph$tobgp)
   tt1[tt1 == 1] <- esoph$ncontrols
   tt <- array(c(ttt, tt1), c(dim(ttt),2), 
      c(dimnames(ttt), list(c("Болен", "Здоров"))))
   mosaicplot(tt, main="Заболеваемость раком", color = TRUE)

Особенностью этого примера является использование такого достаточно специфического типа данных, как порядковые категориальные переменные:

   str(esoph)
   'data.frame':   88 obs. of  5 variables:
     $ agegp    : Ord.factor w/ 6 levels "25-34"<"35-44"<..: 1 1 1 1 1 1 1 1 1 1 ...
     $ alcgp    : Ord.factor w/ 4 levels "0-39g/day"<"40-79"<..: 1 1 1 1 2 2 2 2 3 3 ...
     $ tobgp    : Ord.factor w/ 4 levels "0-9g/day"<"10-19"<..: 1 2 3 4 1 2 3 4 1 2 ...
    $ ncases   : num  0 0 0 0 0 0 0 0 0 0 ...
    $ ncontrols: num  40 10 6 5 27 7 4 7 2 1 ...
   fm2 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp + alcgp,
         data = esoph, family = binomial())
   summary(fm2)
   Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
   (Intercept) -1.77997    0.19796  -8.992  < 2e-16 ***
   agegp.L      3.00534    0.65215   4.608 4.06e-06 ***
   agegp.Q     -1.33787    0.59111  -2.263  0.02362 *  
   agegp.C      0.15307    0.44854   0.341  0.73291    
   agegp^4      0.06410    0.30881   0.208  0.83556    
   agegp^5     -0.19363    0.19537  -0.991  0.32164    
   tobgp.L      0.59448    0.19422   3.061  0.00221 ** 
   tobgp.Q      0.06537    0.18811   0.347  0.72823    
   tobgp.C      0.15679    0.18658   0.840  0.40071    
   alcgp.L      1.49185    0.19935   7.484 7.23e-14 ***
   alcgp.Q     -0.22663    0.17952  -1.262  0.20680    
   alcgp.C      0.25463    0.15906   1.601  0.10942    
   ---
   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
       Null deviance: 227.241  on 87  degrees of freedom
   Residual deviance:  53.973  on 76  degrees of freedom
   AIC: 225.45

Модель показывает, насколько возрастает риск заболеть раком при перемещении в следующую по порядку градацию каждого фактора. Построим серию моделей с разными комбинациями состава потенциальных предикторов:

  ms2 <- dredge(fm2))
  Model selection table 
    (Intrc) agegp alcgp tobgp df logLik   AICc  delta  weight
  8 -1.780  +     +     +     12 -100.727 229.6   0.00 0.798 
  4 -1.911  +     +            9 -106.026 232.4   2.75 0.202 
  6 -1.951  +           +      9 -133.754 287.8  58.20 0.000 
  7 -1.297        +     +      7 -139.482 294.4  64.75 0.000 
  3 -1.402        +            4 -143.135 294.8  65.14 0.000 
  2 -2.139  +                  6 -143.297 299.6  70.02 0.000 
  5 -1.417              +      4 -178.506 365.5 135.88 0.000 
  1 -1.584                     1 -187.361 376.8 147.15 0.000 
  Models ranked by AICc(x)

Согласно полученной сводной таблице, можно установить, что оптимальной является модель со всеми тремя факторами, которая идентична исходной полной модели (это, впрочем, вполне соответствует медицинским представлениям). Здесь знаки "+" лишь показывают, какие факторы следует использовать в моделях, а конкретные значения подогнанных коэффициентов для каждой градации фактора скрыты из соображений компактности представления результатов. Полную информацию о коэффициентах моделей можно получить, выполнив команду (attr(ms2,"coefTables").

Проведем небольшой статистический эксперимент, превратив порядковые категориальные переменные alcgp и tobgp в действительные числа и построив модели, включающие полученные целочисленные переменные, а также их парные взаимодействия:

  fm3 <- glm(cbind(ncases, ncontrols) ~ agegp + unclass(tobgp) *
    unclass(alcgp),data = esoph, family = binomial())
  summary(fm3)
  Coefficients:
                                Estimate Std. Error z value Pr(>|z|)    
  (Intercept)                   -4.70284    0.49575  -9.486  < 2e-16 ***
  agegp.L                        2.92868    0.65126   4.497 6.89e-06 ***
  agegp.Q                       -1.34088    0.58964  -2.274  0.02296 *  
  agegp.C                        0.12258    0.44850   0.273  0.78462    
  agegp^4                        0.06880    0.30787   0.223  0.82317    
  agegp^5                       -0.22286    0.19530  -1.141  0.25382    
  unclass(tobgp)                 0.60346    0.20014   3.015  0.00257 ** 
  unclass(alcgp)                 0.94297    0.17916   5.263 1.41e-07 ***
  unclass(tobgp):unclass(alcgp) -0.14099    0.07572  -1.862  0.06262 .  
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
      Null deviance: 227.241  on 87  degrees of freedom
  Residual deviance:  55.833  on 79  degrees of freedom
  AIC: 221.31
  (ms3 <- dredge(fm3))
  Model selection table 
     (Int)  agg unc(alc) unc(tbg) unc(alc):unc(tbg) df logLik   AICc  delta  weight
  16 -4.703 +   0.9430   0.6035   -0.1410           9  -101.657 223.6   0.00 0.612 
  8  -4.011 +   0.6531   0.2616                     8  -103.379 224.6   0.96 0.379 
  4  -3.570 +   0.6824                              7  -108.395 232.2   8.57 0.008 
  6  -2.766 +            0.3293                     7  -134.536 284.5  60.85 0.000 
  15 -4.200     0.9903   0.5544   -0.1445           4  -139.416 287.3  63.69 0.000 
  7  -3.491     0.6920   0.2048                     3  -141.396 289.1  65.46 0.000 
  3  -3.165     0.7221                              2  -144.846 293.8  70.21 0.000 
  2  -2.139 +                                       6  -143.297 299.6  76.01 0.000 
  5  -2.152              0.2957                     2  -179.350 362.8 139.22 0.000 
  1  -1.584                                         1  -187.361 376.8 153.15 0.000 
  Models ranked by AICc(x)

Приведенные расчеты имеют отношение к обсуждению такой нетривиальной проблемы, как статистическая природа повсеместно используемых балльных показателей, которые следует интерпретировать как порядковые категориальные переменные. В частности, по типологии Стивенса операции суммы или разности к ранговым переменным неприменимы, и поэтому некорректно вычислять средний балл, применяемый, например, в фигурном катании. Тем не менее, в нашем примере модели с "деклассированными" переменными оказались существенно более эффективными по AICc-критерию, чем модели, основанных на исходных данных. Кроме того, добавление парных взаимодействий показывает, что много пьющие и курящие люди имеют более низкую вероятность заболеть раком, нежели это предсказывается моделью, включающей только аддитивные эффекты.

Содержание озона в атмосферном воздухе

Набор данных airquality, входящий в базовую комплектацию R, содержит результаты ежедневных метеорологических наблюдений с 1 мая по 31 октября в окрестностях Нью-Йорка. Он включает такие показатели, как солнечная радиация Solar.R, средняя скорость ветра Wind, температура атмосферного воздуха Temp и содержание озона Ozone, которое интерпретируется как зависимая переменная. Чтобы оценить степень нелинейности связи между переменными, выполним визуализацию основных зависимостей с использованием функций удобного пакета visreg:

   data(airquality)
   library(visreg)
   airquality$MonthF <- as.factor(airquality$Month)
   fit <- lm(Ozone ~ Solar.R * MonthF, data=airquality)
   visreg(fit, "MonthF", by="Solar.R", 
      breaks=1, points=list(cex=1, pch=16))
   fit <- lm(Ozone ~ Wind * MonthF, data=airquality)
   visreg(fit, "MonthF", by="Wind", 
      breaks=1, points=list(cex=1, pch=16))

  library(splines)
  fit <- lm(Ozone ~ Solar.R +ns(Wind, df=2)*ns(Temp, df=2), data=airquality)
  visreg2d(fit, "Wind", "Temp", plot.type="persp", col = 3)

Нетрудно заметить, что имеющиеся предикторы не являются независимыми: между ними существуют выраженные взаимодействия, имеющие нелинейный характер. Поэтому построим две исходные регрессионные модели для прогнозирования концентрации озона, первая из которых будет учитывать неаддитивность парных взаимодействий предикторов, а вторая - нелинейность их влияния на отклик. Конечно, можно было бы построить и единую модель, включающую сразу все необходимые трансформации базовых переменных, но мы хотим показать, как работает функция merge(), выполняющая объединение двух списков частных моделей.

Перед началом расчетов выполним очистку набора данных от строк, содержащих пропущенные значения:

   airquality <- airquality[complete.cases(airquality),]
   airquality$Month.Day <- airquality$Month + airquality$Day/31
   mod1 <- lm(Ozone ~ Month.Day + (Solar.R + Wind + Temp )^2,
       data=airquality)
   summary(mod1)
   Coefficients:
                  Estimate Std. Error t value Pr(>|t|)   
   (Intercept)  -1.421e+02  6.404e+01  -2.218  0.02872 * 
   Month.Day    -1.690e+00  1.368e+00  -1.235  0.21963   
   Solar.R      -1.978e-01  2.114e-01  -0.936  0.35158   
   Wind          1.056e+01  4.280e+00   2.467  0.01527 * 
   Temp          2.527e+00  8.472e-01   2.982  0.00357 **
   Solar.R:Wind -7.228e-03  6.671e-03  -1.083  0.28113   
   Solar.R:Temp  4.594e-03  2.468e-03   1.862  0.06549 . 
   Wind:Temp    -1.612e-01  5.881e-02  -2.741  0.00723 **
   ---
   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
   
   Residual standard error: 19.12 on 103 degrees of freedom
   Multiple R-squared:  0.6909,    Adjusted R-squared:  0.6699 
   F-statistic: 32.89 on 7 and 103 DF,  p-value: < 2.2e-16

На основе этой "полной" модели выполним построение ансамбля моделей с использоваением всех возможных комбинаций сформированных переменных. Поскольку число исходных параметров $m = 7$, то наш ансамбль будет содержать $2^7 = 128$ моделей. Отсортируем его по убыванию BIC-критерия и с использованием функции subset() выберем подмножество наиболее адекватных моделей с delta < 4:

   ms1 <- dredge(mod1, rank=BIC)
   subset(ms1,delta < 4)
   Global model call: lm(formula = Ozone ~ Month.Day + 
            (Solar.R + Wind + Temp)^2, data = airquality)
   ---
   Model selection table 
     (Int) Mnt.Day    Slr.R   Tmp   Wnd Slr.R:Tmp Slr.R:Wnd Tmp:Wnd df   logLik    BIC delta weight
95  -136.8         -0.35310 2.451 11.15  0.005717           -0.1863  7 -482.309  997.6  0.00  0.435
79  -245.1          0.06599 3.914 14.38                     -0.2280  6 -485.303  998.9  1.28  0.230
96  -138.1  -1.690 -0.32480 2.656 11.15  0.005250           -0.1862  8 -481.502 1000.7  3.10  0.093
111 -232.6          0.18250 3.476 12.96           -0.010670 -0.1839  7 -483.930 1000.8  3.24  0.086
80  -235.4  -2.165  0.05826 4.022 14.05                     -0.2235  7 -484.010 1001.0  3.40  0.079
127 -140.8         -0.22600 2.322 10.55  0.005061 -0.007231 -0.1613  8 -481.689 1001.1  3.47  0.077
Models ranked by BIC(x)

Для построения второй модели сформируем дополнительные переменные как квадраты исходных предикторов, а для отражения сезонности в данных воспользуемся функцией синусоиды для верхнего полупериода.

   airquality$Solar.R.2 <- airquality$Solar.R^2
   airquality$Wind.2 <- airquality$Wind^2
   airquality$Temp.2 <- airquality$Temp^2
   norm.0_1 <- function(x) (x-min(x))/(max(x)-min(x))
   airquality$sin.period <- sin(pi*norm.0_1(airquality$Month.Day))
   mod2 <- lm(Ozone ~ Solar.R + Wind + Temp +  Solar.R.2 + Wind.2
       + Temp.2 + sin.period, data=airquality)
   summary(mod2)
   Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
   (Intercept)  3.541e+02  1.046e+02   3.384  0.00101 ** 
   Solar.R      1.810e-01  9.041e-02   2.001  0.04798 *  
   Wind        -1.330e+01  2.280e+00  -5.831 6.37e-08 ***
   Temp        -8.007e+00  2.812e+00  -2.847  0.00532 ** 
   Solar.R.2   -3.316e-04  2.608e-04  -1.272  0.20637    
   Wind.2       4.599e-01  9.986e-02   4.605 1.18e-05 ***
   Temp.2       5.936e-02  1.805e-02   3.289  0.00137 ** 
   sin.period   1.326e+01  7.046e+00   1.882  0.06264 .  
   ---
   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
   
   Residual standard error: 18.08 on 103 degrees of freedom
   Multiple R-squared:  0.7236,    Adjusted R-squared:  0.7048 
   F-statistic: 38.52 on 7 and 103 DF,  p-value: < 2.2e-16

На основе второй "полной" модели проведем построение еще одного ансамбля, также включающего 128 моделей. С использованием функции merge() объединим оба списка моделей и получим объединенный ансамбль из 164 моделей (в обоих списках модели, состоящие из базовых переменных Solar.R + Wind + Temp, совпадают и в объединенном списке их дубли исключаются). С помощью функции subset() выберем подмножество из 17 наиболее адекватных моделей с накопленной суммой весов cumsum(weight) < 0.95 (в следующем сообщении мы покажем, что эта процедура связана с оценкой доверительных интервалов).

   allmod <- merge(ms1, dredge(mod2, rank=BIC))
   dim(allmod)
   [1] 164  17
   ansmod <- subset(allmod, cumsum(weight) < 0.95)  
   dim(ansmod)
   [1] 17 17

Функции dredge(), merge(), subset() и используемая ниже функция model.sel() возвращают специальный объект класса model.selection, весьма похожий на обычную таблицу класса data.frame. Мы может использовать компоненты этого объекта для вывода нужных нам блоков информации. Предварительно, разумеется, стоит ознакомиться с "внутренностями" объекта, выполнив команду str():

 # Извлекаем список формул моделей
   ModList <- as.character(attr(ansmod,"model.calls"))
   # Из списка формул выделяем только предикторную часть
   ModList <- substring(ModList, regexpr("~", ModList)+1,
       regexpr(",", ModList) - 5)
   # Для компактности удаляем пробелы
   ModList <-gsub(" ", "", ModList, fixed = TRUE)
   # Объединяем данные
   data.frame(logLik=ansmod$logLik, BIC=ansmod$BIC,
       delta=ansmod$delta, weight=ansmod$weight, Formula=ModList)
  logLik   BIC  delta weight                                              Formula
1  -476.8 986.7 0.0000 0.2688                      Solar.R+Temp+Temp.2+Wind+Wind.2
2  -479.6 987.6 0.8917 0.1721                           Solar.R+Temp.2+Wind+Wind.2
3  -477.7 988.5 1.7847 0.1101                    Solar.R.2+Temp+Temp.2+Wind+Wind.2
4  -475.5 988.7 2.0011 0.0988           sin.period+Solar.R+Temp+Temp.2+Wind+Wind.2
5  -480.2 988.8 2.0989 0.0941                         Solar.R.2+Temp.2+Wind+Wind.2
6  -481.0 990.4 3.6889 0.0425                             Solar.R+Temp+Wind+Wind.2
7  -476.5 990.7 4.0253 0.0359            Solar.R+Solar.R.2+Temp+Temp.2+Wind+Wind.2
8  -476.7 991.2 4.5070 0.0282         sin.period+Solar.R.2+Temp+Temp.2+Wind+Wind.2
9  -479.1 991.3 4.5954 0.0270                sin.period+Solar.R+Temp.2+Wind+Wind.2
10 -481.6 991.6 4.8897 0.0233                           Solar.R.2+Temp+Wind+Wind.2
11 -474.6 991.7 4.9815 0.0222 sin.period+Solar.R+Solar.R.2+Temp+Temp.2+Wind+Wind.2
12 -484.1 991.9 5.1737 0.0202                                   Temp.2+Wind+Wind.2
13 -479.5 992.0 5.3085 0.0189                 Solar.R+Solar.R.2+Temp.2+Wind+Wind.2
14 -479.9 992.8 6.1169 0.0126              sin.period+Solar.R.2+Temp.2+Wind+Wind.2
15 -482.3 992.8 6.1473 0.0124                              Temp+Temp.2+Wind+Wind.2
16 -480.5 994.0 7.3394 0.0068                  sin.period+Solar.R+Temp+Wind+Wind.2
17 -485.4 994.4 7.7390 0.0056                                     Temp+Wind+Wind.2

Отметим, что, рассматривая информационные потери по Кульбаку, полученная выше оптимальная модель в $\exp((1012 - 986.7)/2) = 347000 $ раз лучше, чем линейная модель с тремя базовыми предикторами, и в $ exp((1101 - 986.7)/2) = 8 \times 10^{24}$ раз эффективнее по сравнению с нуль-моделью без параметров (К. Бёрнхэм и Д. Андерсон считают, что в том же соотношении находятся и вероятности соответствующих статистических гипотез). Далее заметим, что в окончательный список вошли только предискторы из второй модели, а относительная важность отдельных предикторов приобрела следующие значения:

importance(ansmod)
                         Wind Wind.2 Temp.2 Solar.R Temp Solar.R.2 sin.period
    Importance:          1.00 1.00   0.92   0.69    0.65 0.35      0.20      
    N containing models:   17   17     13      9      11    8         6

Для полноты изложения покажем, как к сформированному ансамблю ansmod добавить обобщенную аддитивную модель, в которой варьирование трех основных предикторов аппроксимируется в виде сглаживающих p-сплайнов:

   library(mgcv)
   mgam <- gam(Ozone ~ s(Solar.R) + s(Wind) + s(Temp), 
      data = airquality)
   summary(mgam)
   Parametric coefficients:
               Estimate Std. Error t value Pr(>|t|)    
   (Intercept)   42.099      1.663   25.32   <2e-16 ***
   ---
   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
   
   Approximate significance of smooth terms:
                edf Ref.df      F  p-value    
   s(Solar.R) 2.760  3.447  3.967  0.00769 ** 
   s(Wind)    2.910  3.657 13.768 1.52e-08 ***
   s(Temp)    3.833  4.753 11.769 7.40e-09 ***
   ---
   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
   
   R-sq.(adj) =  0.723   Deviance explained = 74.7%
   GCV score =  338.9  Scale est. = 306.83    n = 111
   BIC(mgam)
   [1] 993.7625
   gansmod <- merge(ansmod, model.sel(mgam, rank=BIC))

К сожалению, без тщательного подбора характеристик сглаживающих функций модель GAM продемонстрировала мало впечатляющую эффективность и в списке моделей заняла скромное 16 место. Однако наша задача заключалась лишь в иллюстрации работы функции model.sel(), которая создает объект model.selection, объединяя произвольное количество частных моделей.

Рекомендуемая литература

Akaike H (1973) Information theory as an extension of the maximum likelihood principle. In: Petrov BN, Csaki F (eds) Second international symposium on information theory. Akademiai Kiado, Budapest, pp 267-281
Burnham KP, Anderson DR (2002) Model selection and multimodel inference: a practical information-theoretic approach, 2nd edn. Springer, New York
Claeskens G, Hjort NL (2008) Model Selection and Model Averaging. Cambridge University Press, New York
Kullback S, Leibler RA (1951) On information and sufficiency. Ann Math Stat 22: 79-86

Обзор интересных R-пакетов за ноябрь-декабрь 2017 г.

2017-12-25T15:06:00.000+00:00

CVXR - утилиты для решения задач выпуклой оптимизации. Подробнее см. здесь и здесь.
dbplot - позволяет выполнять вычисления, необходимые для построения ggplot-графиков, на стороне базы данных.
ESTER - последовательная проверка статистических гипотез с использованием информационных критериев. Подробное описание подходов, реализованных в этом пакете, а также примеры кода можно найти здесь.
geoops - набор утилит для работы с данными в формате GeoJSON.
hdf5r - R-интерфейс для HDF5. См. примеры здесь.
iheatmapr - создание сложных интерактивных диаграмм типа "тепловая карта".
imbalance - реализация различных методов обработки несбалансированных наборов данных, необходимой перед построением моделей-классификаторов.
MlBayesOpt - байесовский подход к нахождению оптимальных гиперпараметров для таких алгоритмов, как "машина опорных векторов", "случайный лес" и "XGboost".
nnfor - прогнозирование временных рядов с использованием нейронных сетей.
ppclust - реализация различных методов нечеткой кластеризации.
quokar - набор утилит для диагностики качества моделей квантильной регрессии.
RMariaDB - R-интерфейс для работы с базами данных MariaDB и MySQL.
shinyaframe - визуализация трехмерных диаграмм рассеяния с использованием технологии виртуальной реальности A-Frame.
textrank - реализация алгоритма textrank (аналог алгоритма PageRank для текстовых данных).
usethis - набор утилит для автоматизации разработки R-пакетов.

С наступающим Новым Годом!

Добра и мира вам и вашим близким!

Обзор интересных R-пакетов за октябрь 2017 г.

2017-11-05T20:57:00.000+00:00

fastrtext - R-интерфейс к библиотеке fastText, разработанной исследователями из Facebook. Эта библиотека предназначена для работы с текстовыми данными и построения моделей путем обучения с учителем и без учителя.
googleLanguageR - пакет для работы с несколькими API от Google для обработки естественных языков (в частности, Natural Language API, Cloud Speech API и Cloud Translation API).
JuliaCall - R-интерфейс к языку Julia.
OutliersO3 - набор методов для обнаружения выбросов в наборах данных с большим количеством переменных.
powerlmm - анализ мощности для панельных данных с двумя или тремя уровнями случайных эффектов.
re2r - пакет для работы с регулярными выражениями, в основе которого лежит алгоритм RE2 от Google. См. примеры здесь.
slowraker - реализация алгоритма RAKE, предназначенного для автоматического извлечения ключевых слов из отдельных текстовых документов. См. примеры здесь.
tscount - построение моделей для временных рядов, представленных счетными переменными.
udpipe - пакет для выполнения распространенных задач, возникающих при обработке естественных языков (выделение токенов, лемматизация, определение части речи и др.).
Здесь можно ознакомиться со списком из 44 пакетов для R, обеспечивающих доступ к источникам данных разной тематики.

Обзор интересных R-пакетов за сентябрь 2017 г. (и другие новости)

2017-10-08T19:50:00.000+01:00

billboarder - R-интерфейс к billboard.js - JavaScript-библиотеке для интерактивной визуализации данных. См. примеры здесь и здесь.
blastula - хороший инструмент для формирования и отправки электронных писем с HTML-содержимым, которое будет правильно отображаться на любых устройствах. См. примеры здесь.
blogdown - пакет от RStudio, позволяющий использовать синтаксис R Markdown для подготовки статей для блогов. С полной документацией и примерами можно ознакомиться здесь.
cetcolor - набор хорошо воспринимаемых цветовых палитр, разработанных в соответствии с принципами, которые были описаны в работе Kovesi (2015).
esvis - набор функций для оценивания и визуализации "размера эффекта", вызванного тем или иным фактором (факторами) в экспериментальных и обсервационных исследованиях.
googledrive - набор утилит для работы с Google Drive из среды R.
Knoema - R-интерфейс к сервису Knoema, предоставляющему доступ к огромной коллекции публичных данных разной тематики (~2.5 миллиарда временных рядов из тысяч источников).
missRanger - реализация одного из эффективных алгоритмов заполнения пропущенных значений в наборах данных с разнотипными переменными.
skpr - инструмент для нахождения оптимального плана эксперимента. Имеется графический интерфейс, реализованный в виде Shiny-приложения.
spm - набор функций для моделирования пространственных данных с помощью как традиционных (random forests, GBM), так и гибридных методов машинного обучения.

Из других хороших новостей: в московском издательстве ДМК Пресс вышла книга "Машинное обучение с использованием библиотеки H2O" (ориг. назв. "Practical Machine Learning with H2O") в переводе, выполненным Андреем Огурцовым (автор блога "Биостатистика и язык R"). В аннотации к книге говорится следующее:

"Машинное обучение наконец-то достигло стадии зрелости. При помощи программного обеспечения H2O вы можете решать задачи машинного обучения и анализа данных с использованием простого в использовании и открытого (open source) фреймворка, который поддерживает большое количество операционный систем и языков программирования, а также масштабируется для обработки больших данных. Это практическое руководство научит вас использовать алгоритмы машинного обучения, реализованные в H2O, с упором на наиболее важные для продуктивной работы аспекты.

Если вы умеете программировать на R или Python, хотя бы немного знаете статистику и имеете опыт обработки данных, эта книга Даррена Кука познакомит вас с основами использования H2O и поможет вам поэкспериментировать с машинным обучением на наборах данных разного размера. Вы изучите несколько современных алгоритмов машинного обучения: глубокое обучение, «случайный лес», обучение на неразмеченных данных и ансамбли моделей.

Прочтя эту книгу, вы:

• узнаете, как импортировать данные в H2O, преобразовывать их и экспортировать их из H2O;

• изучите основные концепции машинного обучения, такие как перекрестная проверка и проверочные наборы данных;

• поработаете с тремя разными наборами данных, решая задачи регрессии, бинарной и многоклассовой классификации;

• используете H2O для анализа каждого набора данных при помощи четырех алгоритмов машинного обучения;

• поймете, как работает кластерный анализ и другие алгоритмы обучения на неразмеченных данных."

Следует добавить, что А. Огурцов выполнил не только отличный перевод книги, но и подготовил некоторые дополнительные и очень полезные материалы по H2O (в частности, по Deep Water), которые можно найти в соответствующем GitHub-репозитории.

Отчет с конференции "Enterprise Applications of the R Language"

2017-09-23T22:48:00.000+01:00

12-14 сентября в Лондоне прошла ежегодная конференция "Enterprise Applications of the R Language" (EARL) - один из наиболее значимых форумов, где пользователи R делятся своим опытом применения этой системы статистических вычислений для решения бизнес-задач.

Конференция традиционно началась с серии практических семинаров, которые в этом году были посвящены следующим темам:

Организация работы со Spark'ом в среде R с помощью пакета sparklyr
Работа с пакетом MicrosoftML
Веб-скрейпинг и анализ текстовых данных с помощью R
Работа с GitHub
Принципы написания функций в R
Введение в Shiny

В последующие два дня было представлено около 60 докладов с примерами использования R в самых разнообразных индустриях - страхование, банковское дело, маркетинг, программное обеспечение, спорт, медицина, автомобилестроение и др. Слайды этих сообщений можно найти на сайте конференции (см. раздел "Speakers"). Поделюсь некоторыми общими наблюдениями:

Два-три года назад было много обсуждений того, как основанные на R приложения можно запустить в производство. В этот раз таких обсуждений почти не было - то ли народ уже определился с "оптимальными" способами развертывания подобных приложений в своих организациях, то ли просто так получилось, что доклады были в основном посвящены другим темам. Одно из исключений - доклад вашего покорного слуги, в котором вместе с моим соавтором и коллегой Дмитрием Дерябиным мы привели детальный пример создания автоматизированного прогнозного приложения, написанного на R и развернутого в "облаке" Amazon с использованием целого ряда AWS-сервисов.
R начинают широко использовать в таких традиционно консервативных в отношении новых технологий индустриях, как банковское дело и страхование, а также в государственных организациях (во всяком случае, в Великобритании).
Фреймворк Shiny стал де-факто стандартом в мире R для создания интерактивных аналитических веб-приложений.

В представленных на конференции докладах было упомянуто большое количество интересных и полезных пакетов для R, на которые стоит обратить внимание, в частности: prophet, shiny.semantic, plumber, mlogit, ChoiceModelR, markovchain, keras, here, purrr, ChannelAttribution, dataCompareR.

Автор: Сергей Мастицкий

Обзор интересных R-пакетов за август 2017 г.

2017-09-03T17:53:00.000+01:00

brms - очень удобный в использовании пакет для построения регрессионных моделей, параметры которых оцениваются с использованием принципов байесовской статистики. По сути, это "обертка" вокруг известного фреймворка для создания байесовских моделей Stan.
blandr - пакет для выполнения анализа Бленда-Альтмана (Bland & Altman 1986). См. примеры здесь.
ggjoy - расширение для ggplot2, с помощью которого можно наглядно изображать распределения количественных переменных в нескольких группах в духе "joy-диаграмм" (диаграммы, которые своим внешним видом напоминают обложку известного альбома "Unknown Pleasures" группы Joy Division).
keras - высокоуровневый R-интерфейс к Keras - известному фреймворку для построения предсказательных моделей на основе нейронных сетей.
metaBMA - реализация метода "байесовского усреднения моделей" (Bayesian model averaging), который можно применять при выполнении мета-анализа.
mize - набор функций, реализующих некоторые распространенные алгоритмы оптимизации (CG, BFGS, L-BFGS) для параметров с неограниченным множеством возможных значений.
randomForestExplainer - набор функций для расчета и визуализации нескольких показателей важности предикторов в моделях, построенных с использованием алгоритма "случайный лес".
rsample - пакет, который окажется полезным при необходимости создания и анализа многократных повторных выборок из некоторого исходного набора данных (например, при выполнении бутстреп-анализа, перекрестной проверки, и т.п.).
sweep - одновременное построение нескольких прогнозных моделей для временных рядов. Пакет основан на принципах "опрятного" (tidy) анализа данных.
timetk - набор полезных функций для работы с временными рядами.

Классический труд по глубокому обучению теперь доступен на русском языке

2017-08-17T18:54:00.002+01:00

В издательстве ДМК Пресс вышла книга Я. Гудфеллоу и др. "Глубокое обучение", которая представляет собой один из наиболее полных и широко цитируемых трудов по этой (очень "горячей" сегодня) теме. В аннотации к книге говорится следующее:

"Глубокое обучение - это вид машинного обучения, наделяющий компьютеры способностью учиться на опыте и понимать мир в терминах иерархии концепций. Поскольку компьютер приобретает знания из опыта, отпадает нужда в человеке-операторе, который формально описывает необходимые компьютеру знания. Иерархическая организация позволяет компьютеру обучаться сложным концепциям, конструируя их из более простых; граф такой иерархии может содержать много уровней. В этой книге читатель найдет широкий обзор тем, изучаемых в глубоком обучении.

Книга содержит математические и концептуальные основы линейной алгебры, теории вероятностей и теории информации, численных расчетов и машинного обучения в том объеме, который необходим для понимания материала. Описываются приемы глубокого обучения, применяемые на практике, в том числе глубокие сети прямого распространения, регуляризация, алгоритмы оптимизации, сверточные сети, моделирование последовательностей, и др. Рассматриваются такие приложения, как обработка естественных языков, распознавание речи, компьютерное зрение, онлайновые рекомендательные системы, биоинформатика и видеоигры. Наконец, описываются перспективные направления исследований: линейные факторные модели, автокодировщики, обучение представлений, структурные вероятностные модели, методы Монте-Карло, статистическая сумма, приближенный вывод и глубокие порождающие модели.

Издание будет полезно студентами и аспирантам, а также опытным программистам, которые хотели бы применить глубокое обучение в составе своих продуктов или платформ."

P.S.: В качестве бонусной информации: на этой неделе на платформе начался новый курс-специализация по глубокому обучению. Курс разработан и преподается известным Эндрю Нг (Andrew Ng) в рамках его нового проекта deeplearning.ai

Обзор интересных R-пакетов за июль 2017 г.

2017-08-05T15:00:00.001+01:00

charlatan: пакет для создания "фейковых" наборов данных, которые могут включать адреса, имена людей, географические координаты, названия должностей, и т.п.
colordistance: пакет для работы с изображениями. В частности, пользователи могут избирательно маскировать определенные пиксели, вычислять количественные меры сходства между несколькими изображениями по присутствующим в них цветам, находить группы сходных изображений (кластерный анализ на основе доминирующих цветов), и т.п. См. примеры здесь, здесь и здесь.
diceR: пакет для выполнения кластерного анализа, в ходе которого одновременно используется несколько алгоритмов кластеризации. Такой подход позволяет выяснить, насколько разные алгоритмы "согласны" (проявляют "консенсус", англ. "consensus clustering") в отношении принадлежности того или иного наблюдения к определенному кластеру. См. примеры здесь.
ggformula: ggplot2-графики, спецификация которых задается c использованием стандартного для R "формульного" синтаксиса. См. многочисленные примеры здесь.
parallelDist: отличный инструмент для выполнения параллельных вычислений матриц сходства/различий ("параллельная" версия базовой R-функции dist()). Реализовано большое количество стандартных метрик сходства/различий.
RStudioConnect теперь поддерживает функционал plumber'а - пакета, с помощью которого можно легко создавать REST API для практически любых приложений, написанных на R.
secrete: позволяет зашифровывать пароли, API-ключи и т.п. важные элементы, используемые в R-пакетах и приложениях.
SentimentAnalysis: удобный пакет для анализа тональности текстов. Содержит несколько встроенных словарей, как общих, так и тематических (например, финансы). Имеется возможность работать не только с англоязычными текстами.
SimMultiCorrData: набор функций для создания "искусственных" наборов данных из нескольких переменных (непрерывные, бинарные, счетные) с заданной пользователем ковариационной матрицей. Такие наборы данных, в частности, являются ключевым компонентом имитационного моделирования.
walker: удобный пакет для построения байесовских регрессионных моделей с динамическими коэффициентами (т.е. коэффициентами, изменяющимися во "времени").

Обзор интересных R-пакетов за июнь 2017 г.

2017-07-08T11:20:00.000+01:00

В июне вышла новая версия R (v3.4.1) под кодовым названием "Single Candle".
dbplyr: этот новый пакет содержит весь код, который ранее отвечал за работу с удаленными базами данных в пакете dplyr. См. пояснения в статье Хэдли Уикхема.
desctable: формирование таблиц с описательными статистиками, как стандартными, так и пользовательскими. Пакет обладает простым синтаксисом в духе dplyr. См. примеры здесь.
dplyr: опубликована новая версия (v0.7) этого популярного пакета. С перечнем основных нововведений можно ознакомиться в блоге компании RStudio.
replyr: удобный dplyr-подобный пакет для работы с "большими данными" (например, с данными, обрабатываемыми на Spark-кластере). См. примеры здесь.
RJSplot: очередной пакет для создания интерактивных графиков и диаграмм средствами R и JavaScript. См. многочисленные примеры здесь.
simglm: позволяет создавать искусственные наборы данных на основе обобщенных линейных моделей с заданной пользователем структурой (включая модели со смешанными эффектами).
simstudy: еще один пакет для создания искусственных наборов данных с заданными пользователем свойствами. См. примеры здесь.
sparklyr: на платформе DataCamp появился новый обучающий курс по работе с этим пакетом.
tidygraph: визуализация графов и дендрограмм с использованием dplyr-подобного синтаксиса. См. примеры здесь.

Визуализация пространственно-распределенных данных с помощью пакетов ggmap и ggplot2

2017-07-05T22:42:00.000+01:00

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

Введение

Целью многих исследований является анализ конфигурации пространственных объектов и отображение ее структуры на картосхемах. В среде R для решения этой задачи часто используют методы и функции пакетов sp, maps, RgoogleMaps и связанных с ними ресурсов. При этом для создания основного слоя карты часто применяются свободно распространяемые shape-файлы, содержащие точечные и контурные графические примитивы, соответствующие отдельным географическим пунктам или регионам (см. пример здесь и здесь). Однако, по сравнению с картами, полученными на основе специализированных геоинформационных систем, таких как ArcGIS ESRI, подобная визуализация может показаться не столь симпатичной.

Ниже рассматриваются некоторые "продвинутые" методы для быстрой визуализации пространственных данных в среде R, основанные на двух идеях: (1) формирование "на лету" статических карт необходимого качества и масштаба с использованием актуальной информации серверов GoogleMap, OpenStreetMap, Stamen Maps или CloudMade и (2) широкое использование грамматики создания графических слоев для отображения необходимой информации на основе функций пакета ggplot2 (Wickham 2009, 2016; Мастицкий 2016). В результате развития этих концепций был разработан удобный пакет R, названный ggmap (Kahle & Wickham 2013). Продемонстрируем некоторые его возможности с использованием собственных данных.

Создание основного картографического слоя

Формирование основого слоя карты в ggmap состоит из двух этапов. На первом этапе с использованием функции get_map() осуществляется загрузка из Интернета необходимых картографических данных, их форматирование для создания наиболее подходящего изображения и формирование объекта ggmap, включающего полученный растр и связанные с ним атрибуты. На втором этапе с использованием функции qmap() осуществляется вывод на задаваемое графическое устройство основы карты, описанной объектом ggmap.

Источник данных, из которого осуществляется загрузка картоосновы, определяется при вызове функции get_map() параметром source = c("google", "osm", "stamen", "cloudmade"), значения которого соответствуют вышеперечисленным серверам. Кроме того, для загрузки из каждого ресурса можно использовать также специальные функции get_googlemap(), get_openstreetmap(), get_stamenmap() или get_cloudmademap(), для которых список используемых параметров становится более точно определенным. Все они возвращают объект ggmap стандартизованной структуры, однако далее мы будем рассматривать только карты, загруженные с сервисов GoogleMap и Stamen Maps, поскольку использование OpenStreetMap и CloudMade требуют установки дополнительных модулей или предварительной регистрации на сервере и получения доступа к соответствующим API.

Географическая локализация карты выполняется с использованием двух параметров: location и zoom. Для get_googlemap() стандартом location является пара координат "долгота/широта", определяющих центр карты. Величина радиуса пространственного охвата вокруг этого центра определяется аргументом масштабирования zoom, который является целым числом от 3 до 20. При zoom = 3 получаем примерно уровень континента, а при zoom = 20 имеем масштаб отдельно стоящих зданий.

Третьим важным аргументом является maptype, который определяет набор загружаемых с сервера графических компонентов и их эстетику. Для GoogleMap возможны следующие аргументы maptype = c("terrain", "satellite", "roadmap", "hybrid"). Рассмотрим пример создания основного слоя карты, выводящей фрагмент территории г. Саратова, на улицах которого, как уверяют, "огней так много золотых". Координаты населенного пункта определим с использованием функции geocode():

library(ggmap)
(reg_center = geocode("Saratov"))
      lon      lat
1 45.9608 51.59237
ggmap(get_map(location = as.numeric(reg_center), source = "google", zoom = 14, 
      maptype = "roadmap"))

После подготовки изображения основного картографического слоя можно приступить к нанесению дополнительных слоев данных, а также различных статистик или модельных объектов, полученных с использованием функций ggplot2.

Отображение точечных данных

В качестве примера будем использовать базу гидробиологических данных Института экологии Волжского бассейна РАН. В ходе экспедиционных исследований было взято 1158 проб донных сообществ беспозвоночных организмов в 100 малых и средних реках Среднего и Нижнего Поволжья, после чего результаты гидробиологической съемки были загружены в таблицы базы Microsoft Access. Покажем, как можно работать в среде R с данными, представленными в подобном формате. К сожалению, предоставление этой базы данных в публичный доступ невозможно, что не позволит читателям воспроизвести приведенные ниже примеры. Тем не менее, будем надеятся, что эти примеры окажутся полезными, и читатель легко сможет "спроецировать" их на собственные задачи.

Соединение с базой данных выполним через стандартный API-интерфейс Open DatabaseConnectivity, для чего первым делом создадим новый объект ODBC, используя панель администрирования операционной системы Windows. Присвоим объекту ODBC имя "База", свяжем его с актуальным файлом MDB и определим необходимый драйвер. Соединение с источником данных выполним при помощи функции odbcConnect() из пакета RODBC:

library(RODBC)
con <- odbcConnect("База")

Таблица "Координаты" нашей базы данных содержит географические координаты рек и их отдельных участков (для средних рек выделены отдельно истоки, среднее и нижнее течения и устье). Выведем на картосхему названия обследованных рек. Для создания основного слоя карты используем опцию maptype = "satellite":

SQL = "SELECT Наименование, Avg(Широта) AS lat, Avg(Долгота) AS lon 
       FROM Координаты GROUP BY Наименование"
koord <- sqlQuery(con, SQL)
myMap_gs <- ggmap(get_googlemap(center = c(50, 51.4), zoom = 6,
                  maptype = "satellite", extent = "device"))
myMap_gs + geom_text(aes(x = lon, y = lat, label = Наименование),
                     size = 2.7, colour = "green", data = koord) +
           ggtitle("Обследованные реки")

Теперь выведем точками на карте координаты мест отбора проб из таблицы "Координаты", взяв за основу карту из сервиса Stamen Maps. Отметим, что при использовании карт этого типа необходимо установить пакет ggmap не из репозитория CRAN, а из Github-репозитория его разработчиков (подробности см. здесь):

library(devtools)
install_github("dkahle/ggmap")
# --------------------------------
#  Участки отбора проб   
SQL =  "SELECT * FROM [Координаты]"
koord <- sqlQuery(con, SQL)
bbox <- c(left = 45, bottom = 48.5, right = 55, top = 55)
myMap_st <- ggmap(get_stamenmap(bbox, maptype = "terrain"))
myMap_st + geom_point(aes(x = Долгота, y = Широта), data = koord) +
           ggtitle("Точки отбора проб")

Обратим внимание, что локализация карты в этом случае может выполняться иным, часто более удобным способом: с указанием точных географических координат углов выделяемого прямоугольника поверхности Земли. Кроме того, для карт Stamen Maps имеется возможность использовать различные опции параметра maptype = c("terrain", "terrain-background", "terrain-labels", "terrain-lines", "toner", "toner-2010", "toner-2011", "toner-background", "toner-hybrid", "toner-labels", "toner-lines", "toner-lite", "watercolor"). Использование стиля toner приводит к резко контрастным черно-белым картам "газетного" вида, а карты watercolor, предназначенные для визуализации озерных систем, изображаются двумя группами оттенков цвета (желто-коричневого для суши и голубого для водных объектов). Мы предпочитаем использовать terrain или terrain-background, где подавляются все надписи, часто только мешающие анализу карт.

В ходе гидробиологической съемки было обнаружено более 700 видов и более крупных таксонов водных беспозвоночных. Выгрузим из базы показатель Count_Sample, определяющий, сколько раз тот или иной вид встретился в процессе исследований на каждом из участков рек. Например, рассмотрим, как распространены в реках Поволжья личинки кровососущих мокрецов рода Culicoides. Используем для отображения на карте кружки, диаметр которых пропорционален Count_Sample:

Species = 'CeCul.sp'
SQL =  paste("SELECT * FROM [_Численность с координатами] 
             WHERE Code = '", Species, "'", sep = "")
sumdat <- sqlQuery(con, SQL)
name_sp <- as.character(sumdat$NAME.1[1])
head(sumdat[, c(1, 10, 11, 6, 8, 9)])
         NAME   Широта  Долгота         NAME.1 Count_Sample Sum_CHISL
1        Анлы 53.93778 52.33430 Culicoides sp.            1        50
2 Б. Саморода 49.11437 46.84341 Culicoides sp.           19     22140
3  Б.Черемшан 54.30345 52.02079 Culicoides sp.            1        60
4  Б.Черемшан 54.50150 50.39124 Culicoides sp.            3       100
5    Байтуган 54.16788 52.28791 Culicoides sp.           11      2886
6     Домашка 52.93876 50.60469 Culicoides sp.            1        10

myMap_st + geom_point(aes(x = Долгота, y = Широта, size = Count_Sample), data = sumdat) +
           ggtitle(name_sp)

Распределение плотности вероятностей встречаемости видов

Для построения 2D-диаграммы плотности вероятностей выгрузим из базы данных полный список и координаты проб, в которых встретились виды Culicoides:

SQL =  paste("SELECT * FROM [_Проба с координатами] WHERE Code = '", 
             Species, "'", sep = "")
pundat <- sqlQuery(con, SQL)
name_sp <- as.character(pundat$NAME[1])
nrow(pundat)
[1] 222

Графически плотность двухмерного распределения этих точек на карте можно отобразить, используя контурные линии, которые ограничивают участки с одинаковой плотностью вероятности. В пакете ggplot2 для нахождения границ контуров плотности служит функция stat_density2d(), которая в свою очередь основана на функции kde2d() из базового R-пакета MASS. Функция scale_fill_gradient() определяет диапазон цветов заливки каждого из bins = 5 контуров, а функция scale_alpha() управляет уровнем прозрачности в зависимости от градации встречаемости вида:

myMap_st +  ggtitle(name_sp) +
stat_density2d(data = pundat,  aes(x = Долгота, y = Широта, 
               fill = ..level.., alpha = ..level..), bins = 5, geom = "polygon") +
scale_fill_gradient(low = "green", high = "red") +
scale_alpha(range = c(0, 0.3), guide = FALSE)

Предположим теперь, что мы хотим внимательнее рассмотреть ситуацию с распространением Culicoides в районе оз. Эльтон (сгущение в нижней части карты). В дополнение к предыдущему скрипту в отдельном слое расставим местоположения точек отбора проб и добавим легенду с градациями встречаемости личинок мокреца.

ggmap(get_googlemap(center = c(46.64, 49.21), zoom = 10,  
      maptype = "hybrid", extent = "device")) +
ggtitle(name_sp) +
geom_point(data = pundat, aes (x = Долгота, y = Широта), colour = "red") + 
stat_density2d(data = pundat,  aes(x = Долгота, y = Широта, 
               fill = ..level.., alpha = ..level..), bins = 5, geom = "polygon") +
scale_fill_gradient("Плотность\nвстречаемости", low = "green", high = "red") + 
guides(fill = guide_colorbar(barwidth = 1.5, barheight = 10)) +
scale_alpha(range = c(0, 0.3), guide = FALSE)

Литературные источники:

Wickham H. (2009, 2016) ggplot2: Elegant Graphics for Data Analysis. Springer, New York. URL: http://ggplot2.org/book/
Мастицкий С. Э. (2016) Визуализация данных с помощью ggplot2. М.: ДМК Пресс, 2016
Kahle D., Wickham H. (2013) ggmap: Spatial Visualization with ggplot2. The R Journal 5: 144–162. URL: http://vita.had.co.nz/papers/ggmap.html

Обзор интересных R-пакетов за май 2017 г.

2017-06-04T14:03:00.000+01:00

CausalImpact: великолепный инструмент от аналитиков из Google, который позволяет оценить эффект экспериментального воздействия на изменяющуюся во времени переменную при отсутствии контрольной группы.

defer: пакет, с помощью которого можно удобно "обернуть" несколько функций в одну высокоуровневую функцию. Вот небольшой пример:

C <- 10

f <- function(x) x*x + C
g <- function(y) (f(y) + 10) * f(y+3)
h <- function(z) mean(c(g(z), g(z+1), g(z+2)))

wrapper <- defer(h)
#> Found functions:
#>   g, f
#> variables:
#>   C
#> library calls:
#>   base::mean

rm(C, f, g, h)

wrapper(10)
#> [1] 29688.67

d3Tree: R-интерфейс к JavaScript-библиотеке D3 для создания "схлопывающихся" дендрограмм. См. примеры здесь.
MultiVarSel: реализация метода Perrot-Dockes et al. (2017) для селекции информативных предикторов для общих линейных моделей. См. примеры здесь.
networktools: набор функций для анализа графов (сетей). В частности, в этом пакете реализована новая метрика для оценки "влияния" (impact) отдельных узлов и ребер графа на устойчивость его структуры.

philentropy: позволяет вычислить 46 различных метрик сходства/расстояния между многомерными объектами. Вот список реализованных в этом пакете метрик:

library(philentorpy)

getDistMethods()

##  [1] "euclidean"         "manhattan"         "minkowski"         "chebyshev"        
##  [5] "sorensen"          "gower"             "soergel"           "kulczynski_d"     
##  [9] "canberra"          "lorentzian"        "intersection"      "non-intersection" 
## [13] "wavehedges"        "czekanowski"       "motyka"            "kulczynski_s"     
## [17] "tanimoto"          "ruzicka"           "inner_product"     "harmonic_mean"    
## [21] "cosine"            "hassebrook"        "jaccard"           "dice"             
## [25] "fidelity"          "bhattacharyya"     "hellinger"         "matusita"         
## [29] "squared_chord"     "squared_euclidean" "pearson"           "neyman"           
## [33] "squared_chi"       "prob_symm"         "divergence"        "clark"            
## [37] "additive_symm"     "kullback-leibler"  "jeffreys"          "k_divergence"     
## [41] "topsoe"            "jensen-shannon"    "jensen_difference" "taneja"           
## [45] "kumar-johnson"     "avg"

R, а также Python и F# теперь интегрированы в IDE Visual Studio, что, согласно задумке Microsoft, добавит удобства в процесс разработки "умных приложений" на основе этих языков благодаря использованию одной IDE.
shinymaterial: позволяет оформлять Shiny-приложения в духе "материального дизайна" (material design), предложенного Google. Подробнее см. здесь.
simmer: набор функций для моделирования дискретных событий.
Syberia: новый фреймворк для разработки предсказательных моделей средствами R и последующего внедрения этих моделей в бизнес-процессы. См. также видео доклада по Syberia, сделанного недавно на конференции R/Finance.

Обзор интересных R-пакетов за март-апрель 2017 г.

2017-05-10T22:34:00.000+01:00

21 апреля 2017 г. вышла новая версия R (v3.4.0), в которую вошло несколько изменений, связанных с увеличением скорости вычислений. Подробнее можно узнать здесь.
choroplethr: предназначен для создания картограмм. Кроме того, с помощью этого пакета можно выполнять соединение с несколькими удаленными источниками картографических данных.
cronR: позволяет задавать расписание выполнения R-скриптов на Unix/Linux машинах. Есть и аналогичный пакет для Windows-машин - taskscheduleR (см. примеры здесь).
errors: в физике и других научных дисциплинах любое измерение всегда сопровождается определенной погрешностью. В R одновременно работать с результатами подобных измерений и их погрешностями не трудно, но утомительно. Новый пакет errors предлагает удобное решение для работы с такими данными.
HIBPwned: R-интерфейс к API сервиса "Have I Been Pwned?". С помощью этого сервиса можно выяснить, находится ли тот или иной адрес электронной почты в списках адресов, учетные данные которых стали по тем или иным причинам публично доступными.
kerasR: R-интерфейс к популярной Python-библиотеке для "глубинного обучения" Keras. См. примеры здесь.
Компания Microsoft объявила о выпуске новой версии Microsoft R Server (v9.1).
modeleval: набор утилит для анализа качества предсказательных моделей, построенных с помощью пакета caret.
readxl: как следует из названия, этот пакет предназначен для организации удобной работы с файлами Excel из среды R.
RInside: с помощью этого пакета можно встраивать функционал R в приложения, написанные на C++.
Sinew: позволяет автоматически создать "скелет" справочного файла для той или иной функции в формате roxygen2, что значительно сэкономит время при разработке новых пакетов для R.
shinyWidgets: набор дополнительных виджетов для разработки Shiny-приложений.
shinyjqui: с помощью этого пакета можно добавлять jQuery UI-элементы в Shiny-приложения. Примеры приведены здесь и здесь.
Исправлены небольшие ошибки в популярном пакете quantmod, предназначенном для разработки трейдинговых моделей.

Делимся опытом: визуализация данных из лог-файлов

2017-04-30T20:53:00.001+01:00

Сегодня в рамках рубрики "Делимся опытом" мы представляем вашему вниманию небольшую заметку Дмитрия Берга, в которой он рассказывает о визуализации данных из лог-файлов средствами сервиса plot.ly.

Автор: Дмитрий Берг

Анализ логов веб-сайтов и других похожих систем - распространенная задача. В Linux это можно сделать из консоли терминала буквально одной строчкой кода, но возвращаемый при этом массив данных только в текстовом формате не позволяет сформировать полную картину о том, что происходит на сервере. Наглядная визуализация данных гораздо лучше помогла бы увидеть динамику изменений на сайте. В этой статье приводятся примеры подобной визуализации.

Прежде чем продолжить, сделаю несколько комментариев:

Для визуализациии я использую сервис plot.ly, и вся полученная инфографика будет экспортироваться на их платформу. Безусловно, существуют альтернативные варианты размещения и сохранения созданных графиков, но, как мне кажется, этот сервис и соответствующий R-интерфейс к нему (plotly) особенно удобны.
Для подготовки данных я применяю пакет data.table.
Спецификация моей системы: R - версия 3.3.2 (2016-10-31) - "Sincere Pumpkin Patch", RStudio - версия 1.0.136, ОС - Fеdora x86_64.

1. В качестве первого шага подключаем необходимые пакеты:

library("plotly")
library("devtools")
library("data.table")
library("RCurl")
library("XML")

2. Данные для воспроизведения рассматриваемого примера можно скачать по следующей ссылке: https://yadi.sk/d/TCGocq0K3ECZfH. В общем же случае, когда сервер сайта имеет адрес формата "ftp://server/hosting/", для скачивания логов можно выполнить примерно следующие команды:

url <- 'ftp://server/hosting/' # Через ftp подключаемся к серверу сайта
userpwd <- "Логин:Пароль"      # Вводим логин и пароль
filenames <- getURL(url, userpwd = userpwd, ftp.use.epsv = TRUE, 
           dirlistonly = TRUE) # Отправляем запрос на наличие файлов
filenames                      # Просматриваем список найденных файлов

Далее необходимо скачать обнаруженные лог-файлы с расширением .log. Из терминала Linux или OSX это можно сделать с помощью следующих команд:

curl -u логиг:пароль -O ftp://server/hosting/access_file1.log
curl -u логиг:пароль -O ftp://server/hosting/access_file2.log

3. Переходим к обработке данных. Как видите, моя стратегия заключается в том, чтобы сначала обработать каждую таблицу по отдельности, и только на заключительном этапе обьединить их в одну большую таблицу. Эта стратегия вызвана тем, что дальше при желании можно построить отдельные графики для каждой таблицы.

Сохраняем данные в таблицы:

f1 <- read.table("access.file1.log") 
f2 <- read.table("access.file2.log")

При необходимости, удаляем пустые столбцы (в моем случае оказалось два столбца):

f1$V3 <- NULL
f1$V4 <- NULL

В некоторых столбцах пришлось изменить формат:

f1$V5 <- as.POSIXct(strptime("2017-01-27", "%Y-%m-%d"))
f2$V5 <- as.POSIXct(strptime("2017-01-26", "%Y-%m-%d"))
f1$V8 <- as.numeric(df$V8)
f1$V2 <- as.character.numeric_version(df$V2)

Обьединяем таблицы в одну общую:

df <- merge(f1,f2, all = TRUE)

Переименуем названия столбцов на более осмысленные:

colnames(df) <- c('host', 'ip', 'date', 'method',
                  'bit', 'referer', 'useragent')

Просмотрим структуру таблицы:

str(df)

На этом подготовка данных завершена. Переходим к их визуализации.

4. Создадим два графика одновременно. На левом графике мы увидим, с каких IP адресов было выполнено обращение к тем или иным внутренним ресурсам, а на правом - распределение объема данных (в байтах), переданных клиентам/посетителям сайта.

p1 <- plot_ly(x = ~df$ip, y = ~df$referer) %>%
      add_markers(name = "refer/ip")
p2 <- plot_ly(x = ~df$useragent, y = ~df$bit) %>%
      add_markers(name = ~"useragent/bit")
p3 <- subplot(p1, p2)

Отправляем график в свой профиль на сервисе plot.ly (вводить логин и пароль достаточно только один раз - для этого задаются соответствующие глобальные системные переменные с помощью команды Sys.setenv()):

Sys.setenv("plotly_username" = "Логин")          # Вводим логин (имя) plot.ly
Sys.setenv("plotly_api_key" = "API KEY")         # Вводим ключ API plot.ly
plotly_POST(p3,

            filename = "ip-refer-useragent-bit") # Отправляем в личный профиль на plot.ly

Полученная интерактивная диаграмма выглядит следующим образом (см. также здесь):

Теперь создадим 3D график, который покажет нам, какие клиенты (посетители сайта) c уникальным IP адресом получили ответ от сервера, и каков был объем переданных данных:

p <- plot_ly(df, x = ~ip, y = ~useragent, z = ~bit, 
             marker = list(color = ~bit, 
                           colorscale = c('#ffc2c6', '#480607'), 
                           size = 3, showscale = TRUE)) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'IP'),
                      yaxis = list(title = 'User Agent'),
                      zaxis = list(title = 'Bite')),
         annotations = list(
           x = 1.13,
           y = 1.05,
           text = 'Bite',
           xref = 'paper',
           yref = 'paper',
           showarrow = FALSE
         ))

Отправляем график на платформу plot.ly:

plotly_POST(p, filename = "ip-useragent-bit")

Результат выглядит так (см. также здесь):

Скачать файл с полным кодом для воспроизведения приведенных выше примеров можно из моего GitHub-репозитория.

Новая книга по R: "Классификация, регрессия и другие алгоритмы Data Mining с использованием R"

2017-04-07T11:09:00.000+01:00

Друзья!

Владимир Кириллович и я рады представить вашему вниманию новый результат нашей совместной работы - книгу "Классификация, регрессия и другие алгоритмы Data Mining с использованием R".

В книге рассмотрена широкая совокупность методов построения статистических моделей классификации и регрессии для откликов, измеренных в альтернативной, категориальной и метрической шкалах. Подробно рассматриваются деревья решений, машины опорных векторов с различными разделяющими поверхностями, нелинейные формы дискриминантного анализа, искусственные нейронные сети и т.д. Показана технология применения таких методов бутстреп-агрегирования деревьев решений, как бэггинг, случайный лес и бустинг. Представлены различные методы построения ансамблей моделей для коллективного прогнозирования. Особое внимание уделяется сравнительной оценке эффективности и поиску оптимальных областей значений гиперпараметров моделей с использованием пакета caret. Рассматриваются такие алгоритмы Data Mining, как генерация ассоциативных правил и анализ последовательностей. Отдельная глава посвящена методам многомерной ординации данных и различным алгоритмам кластерного анализа.

Описание указанных методов сопровождается многочисленными примерами из различных областей на основе общедоступных исходных данных. Представлены несложные скрипты на языке R, дающие возможность читателю легко воспроизвести все расчеты.

Книга может быть использована в качестве учебного пособия по статистическим методам для студентов и аспирантов высших учебных заведений.

Книга распространяется совершенно бесплатно. Ее PDF-версию и соответствующие приложения можно скачать по следующим двум официальным адресам:

GitHub-репозиторий: https://github.com/ranalytics/data-mining
Сайт Института экологии Волжского бассейна: http://www.ievbras.ru/ecostat/Kiril/R/DM/

Кроме того, книга опубликована в виде отдельного сайта по адресу:

https://ranalytics.github.io/data-mining/

В ближайшее время мы планируем выложить также EPUB- и MOBI-версии, предназначенные для чтения на соответствующих электронных устройствах.

Мы будем благодарны за любые ваши замечания и пожелания касательно этой работы - отправляйте их, пожалуйста, по электронной почте stok1946["собака"]gmail.com и/или rtutorialsbook["собака"]gmail.com

Сергей Мастицкий, Владимир Шитиков

P.S.: Как отмечено выше, книга распространяется бесплатно. Однако если она окажется вам полезной и вы сочтете уместным отблагодарить авторов за их работу, вы можете перечислить любую сумму, воспользовавшись следующей кнопкой (все транзакции выполняются в безопасном режиме через систему электронных платежей PayPal; наличие у вас аккаунта в этой системе необязательно):

Построение деревьев регрессии

2017-02-12T14:03:00.000+00:00

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

В одном из предыдущих сообщений были определены основные понятия, преимущества и недостатки деревьев классификации и регрессии, которые являются одним из наиболее популярных методов решения многих практических задач (Breiman at al., 1984; Quinlan, 1986). По своей сути деревья используют "наивный подход" (naive approach) в том смысле, что исходят из предположения о взаимной независимости признаков. Поэтому модели регрессионных деревьев статистически наиболее работоспособны, когда мультиколлинеарность в комплексе анализируемых переменных выражена незначительно.

Алгоритм CART (Classification and Regression Trees) рекурсивно делит исходный набор данных на подмножества, которые становятся все более и более гомогенными относительно определенных признаков, в результате чего формируется древовидная иерархическая структура. Деление осуществляется на основе традиционных логических правил вида "ЕСЛИ (А) ТО (В)", где А - некоторое логическое условие, а В - процедура разбиения подмножества на две части, для одной из которых условие А истинно, а для другой - ложно. Примеры условий: X_i == F, X_i <= V; X_i >= V и др., где X_i - один из предикторов, F - выбранное значение категориальной переменной, V - специально подобранное опорное значение (порог) количественной переменной.

На первой итерации корневой узел дерева связывается с наиболее оптимальным условным суждением, и все множество объектов разбивается на две группы. От каждого последующего узла-родителя к узлам-потомкам также может отходить по две ветви, в свою очередь связанные c граничными значениями других наиболее подходящих переменных и определяющие правила дальнейшего разбиения (splitting criteria). Конечные узлы (terminal nodes) дерева (также известны как "листья", leaves) представлены совокупностями наблюдений из обучающей выборки, полученными в результате проверки всех предыдущих условных суждений. Общее правило выбора опорного значения для каждого узла построенного дерева можно сформулировать следующим образом: "выбранный признак должен разбить множество Х* так, чтобы получаемые в итоге подмножества Х*k, k = 1, 2, ..., p, состояли из наблюдений, принадлежащих к одному классу, или были максимально приближены к этому".

Описанный процесс относится к так называемым "жадным" алгоритмам, стремящимся построить максимально "кустистое" дерево (также "глубокое дерево", deep tree). Естественно, чем "кустистее" дерево, тем лучше будут результаты его тестирования на обучающей выборке, и, скорее всего, хуже - на незнакомых наблюдениях. Поэтому построенная модель должна быть еще и оптимальной по размерам, т.е. содержать информацию, улучшающую качество прогноза, и игнорировать ту информацию, которая его не улучшает. Для этого обычно проводят "обрезание" дерева (tree pruning) - отсечение ветвей там, где эта процедура не приводит к серьезному возрастанию ошибки.

Невозможно подобрать объективный внутренний критерий, приводящий к хорошему компромиссу между безошибочностью и компактностью модели, поэтому стандартный механизм оптимизации деревьев основан на перекрестной проверке (Loh & Shih, 1997). Для этого обучающая выборка разделяется, например, на 10 равных частей: 9 частей используются для построения дерева, а оставшаяся часть играет роль проверочной выборки. После многократного повторения этой процедуры из некоторого набора деревьев-претендентов с приемлемым (для решаемой задачи) разбросом значений того или иного критерия качества модели, выбирается дерево, показавшее наилучший результат в ходе перекрестной проверки.

Существует несколько алгоритмов построения деревьев, основанных на различной логике и критериях оптимизации. Функция rpart() из одноименного пакета выполняет рекурсивный выбор для каждого следующего узла таких опорных значений, которые приводят к минимальной сумме квадратов внутригрупповых отклонений D_t для всех t узлов дерева. Для оценки качества построенного дерева T в ходе его оптимизации используется следующая совокупность критериев:

СС(T) = Σ(D_t + λt) - штраф на сложность модели (cost complexity), включающий штрафной множитель λ за каждую неотсечённую ветвь;
D₀ - девианс для нулевого дерева (т.е. оценка разброса в исходных данных);
Cp = λ/D₀ - относительный штраф на сложность модели (также известен как "статистика Мэллоу", Mallow's Cp);
REL_er = ΣD_t /D₀ - относительная ошибка обучения для дерева из t узлов;
CV_er - ошибка перекрестной проверки (например, с разбиением на 10 блоков), также отнесенная к девиансу нуль-модели D₀. Как правило CV_erбольше, чем REL_er;
SE - стандартное отклонение для ошибки перекрестной проверки.

Оптимальным обычно считается дерево, состоящее из такого количества ветвей t, для которого сумма (CV_er + SE) является минимальной.

В качестве примера рассмотрим построение дерева CART, прогнозирующего обилие водорослей группы a1 в зависимости от гидрохимических показателей воды и условий отбора проб в различных водотоках (см. описание этих данных в сообщении, посвященном алгоритмам восстановления пропущенных значений). Сначала воспользуемся пакетом rpart. Обычно при построении деревьев регрессии применяется двухшаговая процедура: функция rpart() устанавливает связи между зависимой и независимыми переменными и формирует бинарное дерево, а функция prun() выполняет обрезание лишних ветвей:

library(DMwR)
library(rpart)
data(algae)
algae <- algae[-manyNAs(algae), ]  # Удаляем две записи с NA
(rt.a1 <- rpart(a1 ~ ., data = algae[, 1:12]))

n = 198
node), split, n, deviance, yval
      * denotes terminal node
 1) root 198 90401.290 16.996460  
   2) PO4>=43.818 147 31279.120  8.979592  
     4) Cl>=7.8065 140 21622.830  7.492857  
       8) oPO4>=51.118 84  3441.149  3.846429 *
       9) oPO4< 51.118 56 15389.430 12.962500  
        18) mnO2>=10.05 24  1248.673  6.716667 *
        19) mnO2< 10.05 32 12502.320 17.646880  
          38) NO3>=3.1875 9   257.080  7.866667 *
          39) NO3< 3.1875 23 11047.500 21.473910  
            78) mnO2< 8 13  2919.549 13.807690 *
            79) mnO2>=8 10  6370.704 31.440000 *
     5) Cl< 7.8065 7  3157.769 38.714290 *
   3) PO4< 43.818 51 22442.760 40.103920  
     6) mxPH< 7.87 28 11452.770 33.450000  
      12) mxPH>=7.045 18  5146.169 26.394440 *
      13) mxPH< 7.045 10  3797.645 46.150000 *
     7) mxPH>=7.87 23  8241.110 48.204350  
      14) PO4>=15.177 12  3047.517 38.183330 *
      15) PO4< 15.177 11  2673.945 59.136360 *

Приведенной выше командой мы построили полное дерево из 9 узлов и 10 листьев, обозначенных символом *. В каждой строке представлены по порядку: условие разбиения, число наблюдений, соответствующих этому условию, девианс (в данном случае - это эквивалент суммы квадратов отклонений от группового среднего) и среднее значение отклика для выделенной ветви. Например, перед первым разбиением общее множество из 198 наблюдений имеет среднее значение m = 16.99 при девиансе D = 90401. При PO4 >= 43.8 это множество делится на две части: 2) 147 наблюдений (m = 8.98, D = 31279) и 3) 51 наблюдение с высоким уровнем обилия водорослей (m = 40.1, D = 22442). Дальнейшие разбиения каждой из этих двух частей аналогичны.

Разумеется, лучший вариант - представить дерево графически. Популярны три варианта визуализации с использованием различных функций: plot(), prettyTree() из пакета DMwR и prp() из чрезвычайно продвинутого пакета rpart.plot:

prettyTree(rt.a1, compress = TRUE)

Полезно также проследить изменение перечисленных выше статистических критериев по мере выращивания дерева:

 printcp(rt.a1)
   Variables actually used in tree construction:
   [1] Cl   mnO2 mxPH NO3  oPO4 PO4 
   Root node error: 90401/198 = 456.57
   n= 198 
           CP nsplit rel error  xerror    xstd
   1 0.405740      0   1.00000 1.01753 0.13131
   2 0.071885      1   0.59426 0.76097 0.12507
   3 0.030887      2   0.52237 0.73959 0.12512
   4 0.030408      3   0.49149 0.68751 0.11861
   5 0.027872      4   0.46108 0.67415 0.11851
   6 0.027754      5   0.43321 0.67136 0.11856
   7 0.018124      6   0.40545 0.64653 0.11017
   8 0.016344      7   0.38733 0.66502 0.11160
   9 0.010000      9   0.35464 0.70171 0.11799

Функция rpart() и другие функции из этого пакета имеют собственные возможности выполнить перекрестную проверку и оценить ее ошибку при различных значениях штрафа на сложность модели cp (внимание: ваши результаты могут незначительно отличаться от приведенных в силу эффекта генератора случайных чисел!):

# Снижаем пороговые значения штрафа с шагом .005
rtp.a1 <- rpart(a1 ~ ., data = algae[, 1:12], control = rpart.control(cp = .005)) 
#  График изменения относительных ошибок от числа узлов дерева
plotcp(rtp.a1) 
with(rtp.a1, {lines(cptable[, 2] + 1, cptable[, 3],
     type = "b", col = "red")
     legend(locator(1), c("Ошибка обучения",
     "Ошибка крос-проверки (CV)", "min(CV ошибка)+SE"),
     lty = c(1, 1, 2), col = c("red", "black", "black"), bty = "n") })

На графике видно, что минимум относительной ошибки при перекрестной проверке примерно приходится на значение cp = 0.029. Выполним обрезку дерева при этом значении:

rtp.a1 <- prune(rtp.a1, cp = 0.029)
prettyTree(rtp.a1)

Выполним теперь дополнительную оптимизацию параметра ср с использованием функции train() из пакета caret (подробнее см. здесь). Будем тестировать деревья регрессии при 30 значениях критерия ср, применяя к каждому из них 10-кратную перекрестную проверку с 3 повторностями (как и выше, ваши результаты могут немного отличаться от приведенных):

library(caret)
cvCtrl <- trainControl(method = "repeatedcv", repeats = 3)
rt.a1.train <- train(a1 ~ ., data = na.omit(algae[, 1:12]), 
                     method = "rpart", tuneLength = 30, trControl = cvCtrl)
   184 samples
    15 predictors
   Resampling: Cross-Validation (10 fold, repeated 3 times) 
     cp      RMSE  Rsquared  RMSE SD  Rsquared SD
     0       15.8  0.45      3.19     0.18       
     0.0139  15.8  0.446     3.23     0.189      
     0.0279  15.6  0.459     3.5      0.208      
     0.0418  15.1  0.474     3.64     0.216      
     0.0557  15.2  0.468     3.79     0.216      
     0.0697  15.4  0.454     4.01     0.234      
     0.0836  15.7  0.451     4.18     0.242      
     0.0976  16.2  0.428     4.08     0.214      
     0.111   16.4  0.403     3.86     0.211      
   ...
     0.362   16.4  0.411     3.84     0.213      
     0.376   17    0.388     4.04     0.203      
     0.39    17.5  0.374     3.84     0.201      
     0.404   19.4  0.236     2.93     0.166      
   RMSE was used to select optimal model using smallest value.
   The final value used for the model was cp = 0.0418.
   
plot(rt.a1.train)

rtt.a1 <- rt.a1.train$finalModel
prettyTree(rtt.a1)

При cp = 0.0418 было получено существенно урезанное дерево, которое, правда, значительно потеряло в своей объясняющей ценности (с т.з. биологических механизмов, определяющих обилие водорослей).

Обратимся теперь к принципиально другим методам рекурсивного разбиения, представленным в пакете party. Стандартный механизм проверки статистических гипотез, который предотвращает переусложнение модели, реализован в функции ctree(), использующей метод построения деревьев на основе "условного вывода" (англ. conditional inference). Алгоритм принимает во внимание характер распределения независимых переменных и на каждом шаге рекурсивного разбиения данных осуществляет отбор несмещенного набора предикторов, используя формальный тест на основе статистического критерия Crit(t_j, µ_j, Σ_j), j = 1, ..., m, где, µ, Σ - соответственно среднее и ковариация (Hothorn et al., 2006). Оценка статистической значимости этого критерия выполняется на основе перестановочного теста, в результате чего формируются компактные деревья, не требующие процедуры обрезания.

library(party)  # Построение дерева методом "условного вывода"
(ctree.a1 <- ctree(a1 ~ ., data = algae[, 1:12]))
            
Conditional inference tree with 4 terminal nodes
   Response:  a1 
   Inputs:  season, size, speed, mxPH, mnO2, Cl, NO3, NH4, oPO4, PO4, Chla 
   Number of observations:  198 
   1) PO4 <= 43.5; criterion = 1, statistic = 46.49
     2)*  weights = 50 
   1) PO4 > 43.5
     3) oPO4 <= 51.111; criterion = 0.989, statistic = 10.861
       4) size == {small}; criterion = 0.995, statistic = 15.311
         5)*  weights = 15 
       4) size == {large, medium}
         6)*  weights = 49 
     3) oPO4 > 51.111
       7)*  weights = 84

plot(ctree.a1)

Оптимизацию параметра mincriterion выполним с использованием функции train() при тех же условиях перекрестной проверки:

ctree.a1.train <- train(a1 ~ ., data = algae[, 1:12], 
                        method = "ctree", tuneLength = 10, trControl = cvCtrl)
   Resampling: Cross-Validation (10 fold, repeated 3 times)
     mincriterion  RMSE  Rsquared  RMSE SD  Rsquared SD
     0.01          16    0.422     3.69     0.191      
     0.119         15.9  0.427     3.71     0.193      
     0.228         16    0.421     3.61     0.189      
     0.337         16    0.417     3.54     0.187      
     0.446         15.9  0.426     3.59     0.187      
     0.554         15.9  0.425     3.43     0.172      
     0.663         15.8  0.437     3.3      0.178      
     0.772         15.7  0.442     3.37     0.183      
     0.881         15.6  0.444     3.67     0.18       
     0.99          15.9  0.428     4.1      0.246      
   RMSE was used to select optimal model using smallest value.
   The final value used for the model was mincriterion = 0.881.

plot(ctree.a1.train)
ctreet.a1 <- ctree.a1.train$finalModel
plot(ctreet.a1)

Здесь имел место обратный процесс: число узлов дерева было предложено увеличить с 7 до 9. Обратим также внимание на то, что в дереве появились категориальные переменные (размер и скорость течения реки), которые ранее были проигнорированы rpart-деревьями.

Возникает естественный вопрос: а какому из полученных четырех деревьев следует отдать предпочтение при прогнозировании? Хорошую возможность ответить на этот вопрос предоставляет нам L. Torgo, подготовивший на сайте своей книги набор данных из 140 проверочных наблюдений (см. файл Eval.txt с предикторами и Sols.txt со значениями отклика). Пропущенные значения заполним с использованием алгоритма бэггинга:

Eval <- read.table('Eval.txt', header = F, dec = '.',
                   col.names = c('season', 'size', 'speed', 'mxPH', 'mnO2', 'Cl',
                                 'NO3', 'NH4', 'oPO4', 'PO4', 'Chla'),
                                 na.strings = c('XXXXXXX'))
Sols <- read.table('Sols.txt', header = FALSE, dec = '.',          
                   col.names = c('a1', 'a2', 'a3', 'a4', 'a5', 'a6', 'a7'),
                   na.strings = c('XXXXXXX'))
ImpEval <- preProcess(Eval[, 4:11], method = 'bagImpute')
Eval[, 4:11] <- predict(ImpEval, Eval[, 4:11])

Выполним прогноз для проверочной выборки и оценим точность каждой модели по трем показателям: среднему абсолютному отклонению (MAE), корню из среднеквадратичного отклонения (RSME) и коэффициенту детерминации (Rsq = 1 - NSME , где NSME – относительная ошибка, равная отношению средних квадратов относительно регрессии и общего среднего):

# Функция, выводящая вектор критериев
ModCrit <- function (pred, fact) {
     mae <- mean(abs(pred-fact))
     rmse <- sqrt(mean((pred-fact)^2))
     Rsq <- 1-mean((pred-fact)^2)/mean((mean(fact)-fact)^2)
     c(MAE=mae, RSME=rmse, Rsq=Rsq)
   } 

Result <- rbind(
    rpart_prune = ModCrit(predict(rtp.a1, Eval), Sols[, 1]),
    rpart_train = ModCrit(predict(rt.a1.train, Eval), Sols[, 1]),
    ctree_party = ModCrit(predict(ctree.a1, Eval), Sols[, 1]),
    ctree_train = ModCrit(predict(ctree.a1.train, Eval), Sols[, 1])
   )

Result
                    MAE     RSME      Rsq
   rpart_prune 11.16546 16.09485 0.3828278
   rpart_train 10.72834 15.36578 0.4374751
   ctree_party 11.32286 16.53470 0.3486336
   ctree_train 11.25551 16.40532 0.3587876

Можно с разумной осторожностью сделать вывод о том, что прогнозирующая ценность деревьев после дополнительной оптимизации функцией train() несколько возрастает, а деревья, построенные с помощью rpart(), немного точнее, чем деревья "условного вывода" ctree().

Использованные источники:

Breiman L., Friedman J.H., Olshen R.A. et al. (1984) Classifcation and Regression Trees. Belmont (CA): Wadsworth Int. Group, 368 p.
Quinlan J. R. (1986) Induction of Decision Trees. Machine Learning 1: 81-106
Loh W.-Y, Shih Y.-S. (1997) Split selection methods for classification trees. Statistica Sinica 7: 815-840
Torgo L. (2011) Data mining with R : learning with case studies. Chapman & Hall/CRC, 272 p.
Kuhn M. (2013) Predictive Modeling with R and the caret Package
Hothorn T, Hornik K, Zeileis A. (2006) Unbiased Recursive Partitioning: A Conditional Inference Framework. Journal of Computational and Graphical Statistics 15(3): 651-674
Hothorn T, Hornik K, Zeileis A. ctree: Conditional Inference Trees

Обзор интересных R-пакетов за январь 2017 г.

2017-02-04T18:59:00.001+00:00

Автор: Сергей Мастицкий

Перед тем как привести обзор новых (и обновленных) R-пакетов за январь, стоит сказать о важном событии: в этом месяце общее количество пакетов в хранилище CRAN перевалило за 10000 (эта цифра не включает в себя примерно 1300 пакетов, опубликованных в хранилище BioConductor, а также сотни пакетов, доступных только в репозиториях GitHub). Наличие такого разнообразия пакетов в распоряжении пользователей R - дело, несомненно, хорошее. Однако со временем пользователям становится все сложнее ориентироваться в этом разнообразии и находить нужные пакеты. В этой связи, возможно, стоит напомнить о некоторых инструментах, которые могут оказаться полезными при выполнении соответствующего поиска:

Обзор пакетов для решения ряда распространенных задач, приведенный на сайте CRAN - CRAN Task Views;
Список популярных пакетов на сайте RDocumentation.org;
Список популярных пакетов, опубликованных на плафторме GitHub;
CRANberries - обновляющийся в реальном времени список пакетов, добавляемых в хранилище CRAN;
Инструмент для поиска по хранилищу CRAN, доступный на сайте MRAN (Microsoft R Application Network);
Shiny-приложение для поиска пакетов по категориям CRAN Task Views, созданное Михаилом Поповым;
Функция findFn() из пакета sos.

А теперь, собственно, сам список новых/обновленных пакетов, показавшихся мне особенно интересными:

fst - пакет для выполнения (сверх-)быстрой сериализации и десериализации таблиц данных R.
intrval - набор утилит, облегчающих разбиение количественных переменных на интервалы.
newsflash - R-интерфейс к GDELT Television Explorer API, с помощью которого можно выяснить, сколько раз то и или иное ключевое слово встречалось в новостях телевизионных каналов США за определенное время.
roxygen2 - вышла новая версия этого популярного пакета, предназначенного для написания документации в ходе разработки R-пакетов. Теперь такую документацию можно писать и в формате Markdown.
RPushbullet - новая версия пакета, представляющего собой R-интерфейс к Pushbullet API - сервису для обмена сообщениями между устройствами.
RQGIS - новая версия пакета, представляющего собой R-интерфейс к QGIS API. QGIS - это геоинформационная система с открытым кодом, в которой реализовано более 1000 различных алгоритмов обработки геоданных.
sparklyr - новая версия пакета, предназначенного для работы со Spark-кластерами из среды R.
superheat - новый пакет для построения диаграмм типа "тепловая карта" (heatmap). Имеется возмжность создавать составные рисунки (например, путем добавления к тепловой карте диаграмм рассеяния, диаграмм размахов и т.п.) и выполнять тонкую настройку их элементов.
tensorflow - R-интерфейс к TensorFlow - известной библиотеке от Google для создания предсказательных моделей.
xml2 - новая версия пакета от RStudio для работы с данными в формате XML.

Особенности работы с функцией train() из пакета caret

2017-01-22T14:22:00.003+00:00

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

Как обсуждалось нами ранее, пакет caret (сокращение от Classification and Regression Training) был разработан как эффективная надстройка, позволяющая унифицировать и интегрировать использование множества различных функций и методов построения предсказательных моделей, реализованных в других пакетах R. При этом происходит всестороннее тестирование и оптимизация настраиваемых параметров и гиперпараметров (tuning parameters) моделей. Разработанная единая технология настройки моделей основана на использовании полуавтоматических интеллектуальных подходов и ряда широко применяемых критериев качества, рассчитываемых с применением алгоритмов создания повторных выборок (resampling).

Процесс поиска оптимальных значений параметров моделей с использованием функции train() в общем виде реализуется по следующей схеме:

1. Определение наборов данных и их предварительная обработка (при необходимости)
2. Определение спецификации параметров модели
3. Цикл для каждого параметра модели:
4. | Итеративная оптимизация параметров:
5. | | Выделение и предобработка подвыборок для обучения и тестирования модели
6. | | Подгонка модели по обучающим объектам
7. | | Прогнозирование отклика для тестовых объектов
8. | end
9. | Вычисление показателей средней эффективности прогноза
10. end
11. Установление оптимальных параметров
12. Подгонка итоговой модели по всей выборке с использованием оптимальных параметров

Перед обучением модели при помощи функции train() необходимо задать соответствующий алгоритм и всю совокупность условий процесса оптимизации, для чего функцией trainControl() создается специальный объект. Вызов этой функции с исходными настройками имеет следующий вид:

trainControl(method = "boot",

number = ifelse(grepl("cv", method), 10, 25),

p = 0.75,

repeats = ifelse(grepl("cv", method), 1, number),

search = "grid",

initialWindow = NULL,

horizon = 1,

fixedWindow = TRUE,

verboseIter = FALSE,

returnData = TRUE,

returnResamp = "final",

savePredictions = FALSE,

classProbs = FALSE,

summaryFunction = defaultSummary,

selectionFunction = "best",

seeds = NA,

preProcOptions = list(thresh = 0.95,

ICAcomp = 3,

k = 5),

sampling = NULL,

index = NULL,

indexOut = NULL,

timingSamps = 0,

predictionBounds = rep(FALSE, 2),

adaptive = list(min = 5,

alpha = 0.05,

method = "gls",

complete = TRUE),

trim = FALSE,

allowParallel = TRUE)

Остановимся на описании наиболее важных аргументов этой функции:

method - метод создания повторных выборок: "boot", "boot632", "cv", "repeatedcv", "LOOCV", "LGOCV" (для повторяющихся разбиений на обучающую и контрольную выборки), "none" (проверка качества модели выполняется только на обучающей выборке), "oob" (для таких алгоритмов, как случайные леса, бэггинг и др.), "adaptive_cv", "adaptive_boot" или "adaptive_LGOCV";
number - задает число итераций при создании повторных выборок, в частности - количество таких выборок (k-folds) при перекрестной проверке;
repeats - число повторностей для выполнения k-кратной перекрестной проверки;
p - доля обучающей выборки от общего объема данных (при выполнении k-кратной перекрестной проверки);
verboseIter - TRUE означает, что пользователь будет получать сообщения о ходе вычислений (это удобно для оценки оставшегося времени вычислений, которое зачастую бывает очень большим);
search - задает способ перебора параметров модели при ее настройке - по предварительно заданной сетке ("grid") или случайным образом ("random");
returnResamp и savePredictions - определяют условия сохранения результатов вычислений и предсказанных значений (возможные варианты: без сохранения - "none", сохранение только для итоговой модели - "final", или сохранение всех результатов - "all");
classProbs - при выполнении классификации TRUE означает, что в процессе вычислений будут сохраняться не только конечные метки предсказанного класса, но и значения вероятности принадлежности того или иного наблюдения к каждому из имеющихся классов;
summaryFunction - определяет функцию, которая вычисляет сводную метрику качества модели на основе повторных выборок;
selectionFunction - определяет функцию выбора оптимального значения настраиваемого параметра;
preProcOptions - список опций, который передается на функцию предварительной обработки данных preprocess().

Например, при создании объекта ctrl

ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 10)

параметры перекрестной проверки будут иметь следующий смысл:

method = "repeatedcv" означает, что необходимо выполнить повторную перекрестную проверку (также возможна перекрестная проверка без повторов, проверка по одному наблюдению и др.);
number = 10 означает, что в процессе перекрестной проверки исходные данные необходимо разбить на 10 (примерно) равных частей.
repeats = 10 означает, что перекрестная проверка будет запущена 10 раз.

Теперь перейдем непосредственно к описанию функции train(), которая имеет следующий формат вызова:

train(x, y,

rf",preProcess = NULL,

weights = NULL,

metric = ifelse(is.factor(y), "Accuracy", "RMSE"),

maximize = ifelse(metric %in%

c("RMSE", "logLoss"), FALSE, TRUE),

trControl = trainControl(), tuneGrid = NULL, tuneLength = 3)

Как всегда, исходные данные задаются либо матрицей предикторов и вектором со значениями отклика, либо объектом formula с одновременным указанием таблицы данных data. Следующий аргумент - method - в сущности, определяет модель классификации или регрессии, которую необходимо построить и протестировать. Если выполнить команду names(getModelInfo()), то можно увидеть список из 233 доступных методов (количество этих методов постоянно растет). Тот же список, но с различными возможностями поиска и сортировки можно найти по следующим ссылкам:

http://topepo.github.io/caret/modelList.html или http://topepo.github.io/caret/bytag.html

Пакет caret предоставляет пользователю унифицированный интерфейс доступа к пакетам, содержащим все эти функции. Необходимые компоненты автоматически подгружаются по мере их использования (предполагается, что соответствующие пакеты уже инсталлированы).

Например, следующим образом можно ознакомиться со списком всех моделей, имеющих отношение к линейной регрессии:

ls(getModelInfo(model = "lm"))
[1] "bayesglm"   "elm"        "glm"        "glmboost"   "glmnet"
[6] "glmStepAIC" "lm"         "lmStepAIC"  "plsRglm"    "rlm"

С каждым методом связан набор подлежащих оптимизации (гипер-)параметров. Например, легко убедиться в том, что простая линейная регрессия (method = "lm") не имеет таких параметров:

modelLookup("lm")
   model parameter     label forReg forClass probModel
 1    lm parameter parameter   TRUE    FALSE     FALSE

В свою очередь модель rpart (мы познакомимся с ней подробнее в следующем сообщении, в котором продолжим рассмотрение деревьев решений), имеет один параметр - Complexity Parameter (cp в сокращенном виде):

modelLookup("rpart")
   model parameter                label forReg forClass probModel
 1 rpart        cp Complexity Parameter   TRUE     TRUE      TRUE

Из вывода функции modelLookup() можно также увидеть, что линейная регрессия не используется для классификации (forClass = FALSE), тогда как rpart (от "Recursive Partitioning and Regression Trees") можно применять для построения деревьев как регрессии, так и классификации. В последнем случае модель осуществляет не только предсказание класса, но и оценивает апостериорные вероятности (probModel = TRUE).

Метод перекрестной проверки, заданный объектом trControl = trainControl(), хранит список опций, используемых на каждой итерации в ходе настройки параметров модели и оценки ее качества по определенным критериям. При построении каждой частной модели может осуществляться предварительная обработка данных с использованием методов, перечисленных в preProcess (и с учетом опций preProcOptions объекта trControl).

По умолчанию аргумент metric использует в качестве критерия качества точность предсказания ("Accuracy") в случае классификации и корень из среднеквадратичного отклонения прогнозных значений от наблюдаемых ("RMSE") в случае регрессии. Логический аргумент maximize уточняет, должен ли этот критерий быть максимизирован или минимизирован. Другие значения metric в совокупности с различными аргументами summaryFunction и selectionFunction объекта trControl обеспечивают широкие возможности для определения критериев поиска оптимальных моделей.

Количество перебираемых значений того или иного настраиваемого параметра задается аргументом tuneLength. Например, чтобы задать 30 повторов оценки параметра ср модели rpart, необходимо указать tuneLength = 30. Другой вариант - сохранить последовательность этих значений в отдельной таблице данных и затем подать эту таблицу на аргумент tuneGrid (например, tuneGrid = expand.grid(.cp = 0.5^(1:10))). Последний подход особенно полезен, когда диапазон возможных значений настраиваемого параметра известен заранее.

В результате спецификации всех положенных аргументов функции train() пользователь получает объект класса train, соответствующие элементы которого можно извлечь с помощью оператора $:

ls(mytrain)
[1]  "bestTune"     "call"         "coefnames"    "control"     
[5]  "dots"         "finalModel"   "maximize"     "method"      
[9]  "metric"       "modelInfo"    "modelType"    "perfNames"   
[13] "pred"         "preProcess"   "resample"     "resampledCM" 
[17] "results"      "terms"        "times"        "trainingData"
[21] "xlevels"      "yLimits"

Приведем краткий пример, иллюстрирующий нахождение оптимальной степени полинома для модели зависимости электрического сопротивления (Ом) мякоти фруктов киви от процентного содержания в ней сока (эти данные были подробно рассмотрены нами ранее). В одном из предыдущих сообщений мы выяснили, как найти оптимальную степень полинома d = 5 с использованием самостоятельно написанной функции скользящего контроля. К сожалению, функция train() не позволяет выполнить отбор предикторов для метода "lm" и тем самым выбрать оптимальную степень полинома. Однако мы можем выполнить любую из доступных типов перекрестной проверки и найти оптимальную степень полинома опосредованно, оценив характер изменения таких критериев качества, как например, RMSE или среднего коэффициента детерминации RSquared:

library(caret)
library(DAAG) 
data("fruitohms")

set.seed(123)
max.poly <- 7
degree <- 1:max.poly
RSquared <- rep(0, max.poly)
RMSE <- rep(0, max.poly)

# Выполним 10-кратную перекрестную проверку 10 раз
fitControl <- trainControl(method = "repeatedcv",
                           number = 10, repeats = 10)

# Использование функции train() для полиномиальной регрессии
for (d in degree)  {
     
     f <- bquote(juice ~ poly(ohms, .(d)))
     LinearRegressor <- train(as.formula(f),
     data = fruitohms,
     method = "lm", 
     trControl = fitControl)
     
     RSquared[d] <- LinearRegressor$results$Rsquared
     RMSE[d]<- LinearRegressor$results$RMSE
}

library(ggplot2)
Degree.RegParams = data.frame(degree, RSquared, RMSE)
ggplot(aes(x = degree, y = RSquared), data = Degree.RegParams) + geom_line()
ggplot(aes(x = degree, y = RMSE), data = Degree.RegParams) + geom_line()

Как и ранее, минимум ошибки и максимум коэффициента детерминации имеют место при d = 5. Выполним проверку качества итоговой модели, не указывая непосредственно объект trControl:

Poly5 <- train(juice ~ poly(ohms, 5), data = fruitohms, method = "lm")

summary(Poly5$finalModel)
 Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
 (Intercept)        35.152      0.766  45.893  < 2e-16 ***
 `poly(ohms, 5)1` -148.943      8.666 -17.187  < 2e-16 ***
 `poly(ohms, 5)2`   -5.540      8.666  -0.639  0.52383    
 `poly(ohms, 5)3`   51.078      8.666   5.894 3.43e-08 ***
 `poly(ohms, 5)4`  -13.905      8.666  -1.605  0.11118    
 `poly(ohms, 5)5`  -23.528      8.666  -2.715  0.00759 ** 
 ---
 Signif. codes:  0 С***Т 0.001 С**Т 0.01 С*Т 0.05 С.Т 0.1 С Т 1 
 Residual standard error: 8.666 on 122 degrees of freedom
 Multiple R-squared: 0.7362,     Adjusted R-squared: 0.7254 
 F-statistic:  68.1 on 5 and 122 DF,  p-value: < 2.2e-16 

Poly5
 Resampling: Bootstrap (25 reps) 
   RMSE  Rsquared  RMSE SD  Rsquared SD
   10    0.663     1.87     0.0745

Мы получили в точности те же коэффициенты модели, что и при использовании обычной функции lm(), однако для критериев качества RMSE и RSquared здесь были найдены стандартные ошибки. Эти статистики оценивались на основе 25 бутстреп-выборок (Шитиков, Розенберг 2014), формируемых функцией train() по умолчанию. Обратите внимание, что несмещенное бутстреп-значение коэффициента детерминации (0.663) несколько меньше полученного по полной модели (0.7362).

В заключение приведем сокращенную таблицу со списком моделей (Khun 2008), доступных в пакете caret при использовании функции train() (указаны также наименования оптимизируемых параметров этих моделей). Эта таблица полезна также тем, что является своеобразным путеводителем по пакетам R и реализованным в них статистическим методам:

Модели	Значение method	Пакет	Оптимизируемые параметры
Деревья на основе рекурсивного деления (recursive partitioning)	rpart	rpart	maxdepth
	ctree	party	mincriterion
Деревья на основе бустинга (boosted trees)	gbm	gbm	interaction.depth, n.trees, shrinkage
	blackboost	mboost	maxdepth, mstop
	ada	ada	maxdepth, iter, nu
Другие модели на основе бустинга	glmboost	mboost	mstop
	gamboost	mboost	mstop
	logitboost	caTools	nIter
Случайные леса (random forests)	rf	randomForest	mtry
Случайные леса (random forests)	cforest	party	mtry
Деревья на основе бэггинга (bagged trees)	treebag	ipred
Нейронные сети (neural networks)	nnet	nnet	decay, size
Модели на основе частных наименьших квадратов (partial least squares)	pls	pls, caret	ncomp
Машины опорных векторов с RBF ядром (support vector machines, RBF kernel)	svmRadial	kernlab	sigma, C
Машины опорных векторов с полиномиальным ядром (support vector machines, polynomial kernel)	svmPoly	kernlab	scale, degree, C
Гауссовы процессы с RBF ядром (Gaussian processes, RBF kernel)	gaussprRadial	kernlab	sigma
Гауссовы процессы с полиномиальным ядром (Gaussian processes, polynomial kernel)	gaussprPoly	kernlab	scale, degree
Линейные модели наименьших квадратов (linear least squares)	lm	stats
Многомерные адаптивные регрессионные сплайны (multivariate adaptive regression splines, MARS)	earth, mars	earth	degree, nprune
MARS на основе бэггинга (bagged MARS)	bagEarth	caret, earth	degree, nprune
Эластичные сети (elastic net)	enet	elasticnet	lambda, fraction
Модели лассо (lasso)	lasso	elasticnet	fraction
Машины релевантных векторов с RBF ядром (relevance vector machines, RBF kernel)	rvmRadial	kernlab	sigma
Машины релевантных векторов с полиномиальным ядром (relevance vector machines, polynomial kernel)	rvmPoly	kernlab	scale, degree
Линейный дискриминантный анализ (linear discriminant analysis)	lda	MASS
Пошаговый диагональный дискриминантный анализ (stepwise diagonal discriminant analysis)	sddaLDA, sddaQDA	SDDA
Логистическая регрессия для двух или более классов (logistic/multinomial regression)	multinom	nnet	decay
Регуляризованный дискриминантный анализ (regularized discriminant analysis)	rda	klaR	lambda, gamma
Гибкий дискриминантный анализ (flexible discriminant analysis, FDA)	fda	mda, earth	degree, nprune
FDA на основе бэггинга (Bagged FDA)	bagFDA*	caret, earth	degree, nprune
Машины опорных векторов на основе метода наименьших квадратов (least squares support vector machines)	ssvmRadial	kernlab	sigma
Метод k-ближайших соседей (k nearest neighbours)	knnЗ	caret	k
Разделение по центроидам (nearest shrunken centroids)	pam	pamr	threshold
Наивный байесовский классификатор (naive Bayes classifier)	nb	klaR	usekernel
Обобщенный метод частных наименьших квадратов (generalized partial least squares)	gpls	gpls	K.prov
Сети с квантованием обучающего вектора (learned vector quantization)	lvq	class	k

В будущем мы планируем описать большинство из перечисленных методов и продемонстрировать процесс оптимизации их параметров с использованием функции train().

Использованные источники:

Kuhn M. (2008) Building Predictive Models in R Using the caret Package. Journal of Statistical Software 5:113-142
Шитиков В. К., Розенберг Г. С. (2014) Рандомизация и бутстреп: статистический анализ в биологии и экологии с использованием R. Тольятти: Кассандра, 314 с. (PDF, данные и скрипты доступны на сайте авторов ievbras.ru/ecostat)

Обзор интересных R-пакетов за декабрь 2016 г.

2017-01-08T15:27:00.000+00:00

Автор: Сергей Мастицкий

aws.ses - R-клиент для работы с сервисом AWS Simple Email Service.
censys - R-клиент для поискового сервиса Censys.
clustRcompaR - пакет для выполнения кластеризации документов (например, для обнаружения групп документов, схожих по своей тематике). Позволяет автоматически выделять оптимальное число кластеров и оценивать качество получаемых решений путем сравнения с документами, чья групповая принадлежность известна.
ggmosaic - пакет для создания мозаичных диаграмм в стиле ggplot2.
gofastr - пакет для быстрого формирования терм-документных матриц.
lexicon - набор лексических таблиц, тематических списков слов и словарей (например, наиболее часто встречающиеся имена и фамилии в США, коллекция эмотиконов, список распространенных сокращений и др.). Может оказаться полезным при выполнении анализа текстовых данных.
nanotime - пакет для работы с метками времени с разрешением на уровне наносекунд.
rtrie - набор функций для визуализации префиксных деревьев (англ. trie). Подобные структуры данных широко используются для решения задач, связанных с выполнением быстрого поиска, автокоррекцией слов, рекомендацией следующего слова в предложении и др.
skeletor - создание "скелета" будущего пакета для R. В отличие от стандартных утилит (базовая R-функция package.skeleton() или create() из пакета devtools), skeletor создает "болванки" с более богатым (и полезным) содержимым.
textTinyR - набор функций для выполнения анализа текстовых данных большого объема (создание словарей, терм-документных матриц, подсчет частоты встречаемости слов и т.п. типичные задачи). Обработка выполняется по частям, что позволяет анализировать данные, объем которых превышает имеющийся объем памяти. Примеры можно найти в блоге автора.

Заполнение пропущенных значений в данных

2017-01-04T21:33:00.000+00:00

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

К сожалению, на практике в ходе сбора данных далеко не всегда удается получить полностью укомплектованные их наборы. Пропуски отдельных значений являются повсеместным явлением и поэтому, прежде чем начать применять статистические методы, обрабатываемые данные следует привести к "каноническому" виду. Для этого необходимо, либо удалить фрагменты объектов с недостающими элементами, либо заменить имеющиеся пропуски на некоторые разумные значения.

Проблема "борьбы с пропусками" столь же сложна, как и сама статистика, поскольку в этой области существует впечатляющее множество подходов. В русскоязычных книгах по использованию R (Кабаков, 2014; Мастицкий, Шитиков, 2015) бегло представлены только некоторые функции пакета mice, который, несмотря на свою "продвинутость", мало удобен для практической работы. Хорошей альтернативой являются методы "knnImpute", "bagImpute" и "medianImpute" функции preProcess() из пакета caret, которую мы рассмотрели ранее как инструмент для трансформации данных.

Используем для дальнейших примеров таблицу algae, включенную в пакет DMwR и содержащую данные гидробиологических исследований обилия водорослей в различных реках. Каждое из 200 наблюдений содержит информацию о 18 переменных, в том числе:

три номинальных переменных, описывающих размеры size = c("large", "medium", "small") и скорость течения реки speed = c("high", "low", "medium"), а также время года season = c("autumn", "spring", "summer", "winter"), сопряженное с моментом взятия проб;
8 переменных, составляющих комплекс наблюдаемых гидрохимических показателей: максимальное значение рН mxPH (1), минимальное содержание кислорода mnO2 (2), хлориды Cl (10), нитраты NO3 (2), ионы аммония NH4 (2), орто-фосфаты oPO4 (2), общий минеральный фосфор PO4 (2) и число клеток хлорофилла Chla (12); в скобках приведено число пропущенных значений;
средняя численность каждой из 7 групп водорослей a1-a7 (видовой состав не идентифицировался).

Читатель может самостоятельно воспользоваться функциями описательной статистики summary() или describe() из пакета Hmisc, а мы постараемся поддержать добрую традицию и привести парочку примеров диаграмм, построенных с использованием пакета ggplot2 (подробнее на русском языке см. Мастицкий, 2016):

library(DMwR)
library(ggplot2)

summary(algae) # вывод не приводится

ggplot(data=algae[!is.na(algae$mnO2),], aes(speed , mnO2)) + 
       geom_violin(aes(fill = speed), trim = FALSE, alpha = 0.3) + 
       geom_boxplot(aes(fill = speed), width = 0.2, outlier.colour = NA) + 
       theme(legend.position = "NA")

Мы получили, так называемую "скрипичную" диаграмму (англ. violin plots), которая объединяет в себе идеи диаграмм размахов и кривых распределения вероятности. Суть достаточно проста: продольные края "ящиков с усами" (для сравнения приведены тоже) замещаются кривыми плотности вероятности. В итоге, например, легко можно выяснить не только тот факт, что в потоках с быстрым течением (high) содержание кислорода выше, но и ознакомиться с характером распределения соответствующих значений.

Другой пример - категоризованные графики, удобные для визуализации данных, разбитых на отдельные подмножества (категории), каждое из которых отображается в отдельной диаграмме подходящего типа. Такие диаграммы (или "панели", от англ. panels, facets или multiples) определенным образом упорядочиваются и размещаются на одной странице:

qplot(PO4, a1, data = algae[!is.na(algae$PO4), ]) +
      facet_grid(facets = ~ season)+
      geom_smooth(color="red", se = FALSE)

Из представленных графиков легко сделать предположение о том, что численность водорослей группы а1 падает с увеличением концентрации фосфатов.

Однако в контексте темы этой статьи важно обратить внимание, что мы все время старались блокировать появление пропущенных значений: algae[!is.na(algae$PO4), ]. Если в обрабатываемой таблице обнаружены недостающие данные, то в общих чертах можно избрать одну из следующих возможных стратегий:

удалить строки с неопределенностями;
заполнить неизвестные значения выборочными статистиками соответствующей переменной (среднее, медиана и т.д.), полагая, что взаимосвязь между переменными в имеющемся наборе данных отсутствует (это соответствует известному "наивному" подходу);
заполнить неизвестные значения с учетом корреляции между переменными или меры близости между наблюдениями;
постараться обходить эту неприятную ситуацию, используя, например, формальный параметр na.rm некоторых функций.

Последняя альтернатива является самой ограничивающей, поскольку она далеко не во всех случаях позволяет осуществить необходимый анализ. В свою очередь, удаление строк данных более радикально, но может привести к серьезным потерям важной информации:

# Число строк с пропущенными значениями
nrow(algae[!complete.cases(algae),])
[1] 16

# Их удаление
algae <- na.omit(algae)

Можно удалить не все строки, а только те, в которых число пропущенных значений превышает, например, 20% от общего числа переменных, для чего существует специальная функция из пакета DMwR:

data(algae)

manyNAs(algae, 0.2)
[1]  62 199

algae <- algae[-manyNAs(algae, 0.2), ]

В результате мы удалили только две строки (62-ю и 199-ю), где число пропущенных значений больше одного. Обратите внимание, что выполняя команду data(algae), мы обновляем в памяти содержимое этого набора данных.

Если мы готовы принять гипотезу о том, что зависимости между переменными нет, то простым и часто весьма эффективным способом заполнения пропусков является использование средних значений. В том случае, если есть сомнения в нормальности распределения данных, предпочтительнее использовать медиану. Покажем, как это можно сделать с использованием функции preProcess() из пакета caret:

data(algae)
ind<- apply(algae, 1, function(x) sum(is.na(x))) > 0
algae[ind, 4:11]
    mxPH mnO2    Cl   NO3 NH4    oPO4     PO4  Chla
28  6.80 11.1 9.000 0.630  20   4.000      NA  2.70
38  8.00   NA 1.450 0.810  10   2.500   3.000  0.30
48    NA 12.6 9.000 0.230  10   5.000   6.000  1.10
55  6.60 10.8    NA 3.245  10   1.000   6.500    NA
56  5.60 11.8    NA 2.220   5   1.000   1.000    NA
57  5.70 10.8    NA 2.550  10   1.000   4.000    NA
58  6.60  9.5    NA 1.320  20   1.000   6.000    NA
59  6.60 10.8    NA 2.640  10   2.000  11.000    NA
60  6.60 11.3    NA 4.170  10   1.000   6.000    NA
61  6.50 10.4    NA 5.970  10   2.000  14.000    NA
62  6.40   NA    NA    NA  NA      NA  14.000    NA
63  7.83 11.7 4.083 1.328  18   3.333   6.667    NA
116 9.70 10.8 0.222 0.406  10  22.444  10.111    NA
161 9.00  5.8    NA 0.900 142 102.000 186.000 68.05
184 8.00 10.9 9.055 0.825  40  21.083  56.091    NA
199 8.00  7.6    NA    NA  NA      NA      NA    NA

library(caret)
pPmI <- preProcess(algae[, 4:11], method = 'medianImpute')
algae[, 4:11] <- predict(pPmI, algae[, 4:11])
(Imp.Med <- algae[ind, 4:11])
    mxPH mnO2     Cl   NO3      NH4    oPO4      PO4   Chla
28  6.80 11.1  9.000 0.630  20.0000   4.000 103.2855  2.700
38  8.00  9.8  1.450 0.810  10.0000   2.500   3.0000  0.300
48  8.06 12.6  9.000 0.230  10.0000   5.000   6.0000  1.100
55  6.60 10.8 32.730 3.245  10.0000   1.000   6.5000  5.475
56  5.60 11.8 32.730 2.220   5.0000   1.000   1.0000  5.475
57  5.70 10.8 32.730 2.550  10.0000   1.000   4.0000  5.475
58  6.60  9.5 32.730 1.320  20.0000   1.000   6.0000  5.475
59  6.60 10.8 32.730 2.640  10.0000   2.000  11.0000  5.475
60  6.60 11.3 32.730 4.170  10.0000   1.000   6.0000  5.475
61  6.50 10.4 32.730 5.970  10.0000   2.000  14.0000  5.475
62  6.40  9.8 32.730 2.675 103.1665  40.150  14.0000  5.475
63  7.83 11.7  4.083 1.328  18.0000   3.333   6.6670  5.475
116 9.70 10.8  0.222 0.406  10.0000  22.444  10.1110  5.475
161 9.00  5.8 32.730 0.900 142.0000 102.000 186.0000 68.050
184 8.00 10.9  9.055 0.825  40.0000  21.083  56.0910  5.475
199 8.00  7.6 32.730 2.675 103.1665  40.150 103.2855  5.475

Альтернативой "наивному" подходу является учет структуры связей между переменными. Например, можно воспользоваться тем, что между двумя формами фосфора oPO4 и PO4 существует тесная корреляционная связь. Тогда, например, можно избавиться от некоторых неопределенностей в показателе PO4, вычислив его пропущенные значения по уравнению простой регрессии:

data(algae)
lm(PO4 ~ oPO4, data = algae)

Coefficients:
(Intercept)         oPO4  
     42.897        1.293  

# Функция  вывода значений PO4 в зависимости от оPO4 
fillPO4 <- function(oP) {if (is.na(oP)) return(NA)
      else return(42.897 + 1.293 * oP)
}

# Восстановление пропущенных значений PO4
algae[is.na(algae$PO4), 'PO4'] <- sapply(algae[is.na(algae$PO4), 'oPO4'], fillPO4)
algae[ind, 10]
 [1]  48.069   3.000   6.000   6.500   1.000   4.000   6.000  11.000   6.000
[10]  14.000  14.000   6.667  10.111 186.000  56.091      NA

Одно из пропущенных значений удалось восстановить.

Разумеется, легко придти к мысли не утруждать себя перебором всех возможных корреляций, а учесть все связи скопом. Использование метода "bagImpute" осуществляет для каждой из имеющихся переменных построение множественной бутстреп-агрегированной модели или бэггинг-модели (англ. bagging) на основе деревьев регрессии, принимая все остальные переменные в качестве предикторов. Этот метод мудр и точен, но требует значительных затрат времени на вычисление, особенно при работе с данными большого объема:

data(algae)
pPbI <- preProcess(algae[, 4:11], method = 'bagImpute')
algae[, 4:11] <- predict(pPbI, algae[, 4:11])
Imp.Bag <- algae[ind, 4:11]

Наконец, третий метод функции preProcess() для заполнения пропусков - "knnImpute" - основан на простейшем, но чрезвычайно эффективном алгоритме k ближайших соседей (англ. k-nearest neighbours) или kNN. В основе метода kNN лежит гипотеза о том, что тестируемый объект d будет иметь примерно тот же набор признаков, как и обучающие объекты в локальной области его ближайшего окружения:

Если речь идет о классификации, то неизвестный класс объекта определяется голосованием k его ближайших соседей (на рис. k = 5). kNN-регрессия оценивает неопределенное значение неизвестной координаты Y, усредняя известные ее величины для тех же k соседних точек.

Одна из важных проблем kNN - выбор метрики, на основе которой оценивается близость объектов. Наиболее общей формулой для подсчета расстояния в m-мерном пространстве между объектами X₁ и X₂ является мера Минковского:

D_S(X₁, X₂) = [ Σ |x_1i – x_2i|^p]^1/r,

где i изменяется от 1 до m, а r и p - задаваемые исследователем параметры, с помощью которых можно осуществить нелинейное масштабирование расстояний между объектами. Мера расстояния по Евклиду получается, если принять в метрике Минковского r = p = 2, и является, по-видимому, наиболее общей мерой расстояния, знакомой всем по школьной теореме Пифагора. При r = p = 1 имеем "манхеттенское расстояние", не столь контрастно оценивающее большие разности координат x.

Вторая проблема с kNN заключается в решении вопроса о том, на мнение скольких конкретно соседей k нам целесообразно положиться? В свое время мы обсудим этот вопрос детально, а сейчас будем ориентироваться на значение k = 5, используемое по умолчанию:

data(algae)
pPkI <- preProcess(algae[, 4:11], method = 'knnImpute')
alg.stand <- predict(pPkI, algae[, 4:11])

Получив в результате применения predict() матрицу переменных с пропущенными значениями, заполненными этим методом, мы с удивлением обнаруживаем, что данные оказались стандартизованными (т.е. центрированными и нормированными на дисперсию). Но, а как иначе можно было посчитать меры близости с переменными, измеренными в разных шкалах? Пришлось для возвращения в исходное состояние применить операцию "решкалирования":

m <- pPkI$mean
sd <- pPkI$std
algae[, 4:11] <- t(apply(alg.stand, 1, function (r) m + r * sd))
(Imp.Knn <- algae[ind, 4:11])

Наконец, зададимся следующим закономерным вопросом: а какой метод лучше? Обычно эта проблема не имеет теоретического решения и исследователь полагается на собственную интуицию и опыт. Но мы можем оценить, насколько расходятся между собой результаты, полученные каждым способом заполнения. Для этого сформируем блок данных из 3 * 16 = 48 строк исходной таблицы с заполненными пропусками тремя методами ("Med", "Bag", "Knn") и выполним редукцию переменных методом главных компонент из многомерного пространства в двухмерное. Посмотрим, как "лягут карты" на плоскости:

ImpVal <- rbind(Imp.Med, Imp.Knn)
ImpVal <- rbind(ImpVal, Imp.Bag)
Imp.Metod <- as.factor(c(rep("Med", 16), rep("Knn", 16), rep("Bag", 16)))

library(vegan)
Imp.M <- rda(ImpVal ~ Imp.Metod, ImpVal)
plot(Imp.M, display = "sites", type = "p")
ordihull(Imp.M, Imp.Metod, draw = "polygon", alpha = 67, 
         lty = 2, col = c(1, 2, 3), label = TRUE)

Мы выделили контуром (англ. hull), проведенным через крайние точки, области каждого из трех блоков данных и поместили метку метода в центры тяжести полученных многоугольников. Понятно, что медианное заполнение характеризуется меньшей вариацией результатов, поскольку игнорирует специфичность свойств каждого объекта. Оба других метода, учитывающих внутреннюю структуру данных, дали приблизительно похожие результаты.

Использованные источники:

Кабаков Р.И. R в действии: Анализ и визуализация данных на языке R. М.: ДМК Пресс, 2014. 580 с.
Мастицкий С.Э., Шитиков В.К. Статистический анализ и визуализация данных с помощью R . М.: ДМК Пресс, 2015. 496 с.
Мастицкий С.Э. Визуализация данных с помощью ggplot2. - М.: ДМК Пресс, 2016. 222 с.
Torgo L. Data mining with R : learning with case studies. Chapman & Hall/CRC, 2011. 272 p.

Пять наиболее важных событий в мире R по итогам 2016 г.

2016-12-29T15:56:00.000+00:00

Автор: Сергей Мастицкий

Рисунок заимствован с сайта
https://www.r-consortium.org

Традиционно попытаюсь подвести итоги 2016 г., перечислив пять наиболее важных (на мой взгляд) событий и тенденций, имевших место в этом году в мире R. Если упустил что-то интересное, пожалуйста, сообщите об этом в комментариях.

IDE RStudio v1.0. В ноябре компания RStudio объявила о выходе новой версии (1.0) их одноименного IDE для R. Это был 10-й крупный релиз этого продукта, начиная с февраля 2011 г. В состав RStudio v1.0 вошло большое количество нововведений, среди которых стоит отметить следующие:

работа с R Notebooks - инструментом для интерактивного исполнения R-кода и одновременного документирования получаемых результатов;
интеграция с функционалом пакета sparkyr (см. ниже);
интеграция с функционалом пакета profvis, предназначенного для выполнения профилирования R-кода;
улучшенные средства импорта данных (на основе пакетов readr, readxl и haven);
возможность создавать несложные веб-сайты и документы размером с книгу при помощи языка разметки R Markdown.

Создание пакета sparklyr. Apache Spark - это чрезвычайно популярный инструмент для выполнения вычислений над "большими данными". Разработанный компанией RStudio пакет sparklyr представляет собой R-интерфейс для работы с этим инструментом. Важными особенностями sparklyr (в отличие от, например, встроенного в Spark пакета sparkR) являются следующие:

возможность написания легко читаемого кода с использованием синтаксиса, принятого в пакете dplyr;
пользователи R могут работать с большинством алгоритмов машинного обучения, реализованных в Spark'овской библиотеке MLlib;
как было отмечено выше, sparklyr интегрирован в RStudio v1.0, что делает работу с ним еще удобнее.

Релиз платформы Steam. Компания h20.ai выпустила Steam (под лицензией AGPL) - платформу, которая позволяет "...от начала и до конца выстроить весь процесс создания и разворачивания "умных" приложений". Если говорить "немаркетинговым" языком, то эта платформа помогает наладить совместную работу аналитиков, разрабатывающих предсказательные модели на базе движка h2o (поддерживаются такие языки программирования, как R, Python и Java), и разработчиков приложений (последние "потребляют" предсказания моделей посредством запросов к соответствующим API, которые можно легко опубликовать с помощью Steam). Выстраивание подобной связи между аналитиками и разработчиками эффективным образом - большая головная боль многих компаний, которые пытаются внедрять предсказательные модели в свои бизнес-процессы. В этом смысле Steam - платформа с большим потенциалом, особенно если учесть, что в состав h2o входят одни из наиболее эффективных алгоритмов машинного обучения среди существующих. Подробнее о Steam можно узнать из официальной документации, а также из приведенного ниже видео.

Проект RL10N. Создание Концсорциума R (R Consortium) начинает приносить первые плоды. Одним из восьми проектов, профинансированных этой организацией в 2016 г., стал "RL10N" - проект, целью которого является создание инфраструктуры, облегчающей локализацию R-приложений на языки, отличные от английского. Первым результатом работы над этим проектом стала разработка пакета poio, с помощью которого можно легче выполнять перевод сообщений, генерируемых R-функциями (поддерживаются все языки из списка ISO 639).
Книги по R на русском языке. В московском издательстве ДМК Пресс в 2016 г. вышло три книги по R:

Джеймс Г., Уиттон Д., Хасти Т., Тибширани Р. (2016) Введение в статистическое обучение с примерами на языке R. - М.: ДМК Пресс, 460 с. - пер. с англ. С. Э. Мастицкого.
К сожалению, в первое издание этой книги закралось большое количество опечаток и ошибок технического характера. В начале 2017 г. выйдет второе издание с исправлениями обнаруженных недочетов.
Мастицкий С. Э. (2016) Визуализация данных с помощью ggplot2. - М.: ДМК Пресс, 222 с. См. примеры кода из книги здесь.
Храмов Д. А. (2016) Сбор данных в Интернете на языке R. - М.: ДМК Пресс, 283 с.

С наступающим Новым Годом! Добра и мира вам и вашим близким!

Методы и критерии верификации регрессионных моделей

2016-12-18T22:59:00.002+00:00

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

В предыдущих сообщениях неоднократно поднимались важные проблемы, связанные с диагностикой регрессионных моделей: проверка статистических допущений, лежащих в основе используемого метода построения модели, оценка адекватности структуры систематической части модели, тестирование чувствительности модели к аномалиям в структуре исходных данных и др. (см. часть 1, часть 2, часть 3). В данном сообщения попробуем в сжатой форме дать сводку методов и критериев верификации регрессионных моделей и подвести некоторые итоги того, что было частично рассмотрено ранее.

Действительно, первое, что следует сделать сразу после построения модели, так это убедиться в том, насколько она полезна и эффективна. Проверка работоспособности полученных моделей (англ. model verification) представляет собой сложный и многостадийный процесс, по итогам которого может быть достигнут приемлемый уровень уверенности исследователя в том, что результаты, полученные при практическом использовании модели, окажутся правильными. Верификация включает оценку качества модели (англ. model evaluation) и ее адекватности выборочным данным (англ. model validation) по совокупности объективных критериев, анализ вклада отдельных предикторов и отбор их оптимальной комбинации (англ. model selection) и, наконец, независимое тестирование модели с целью настройки ее параметров (англ. model tuning).

Тестирование регрессионных моделей ставит своей задачей выявление следующих основных проблем:

Смещение (англ. bias), или систематическая ошибка модели (например, сдвиг предсказываемых значений на какую-то трудно объяснимую величину);
Высокая дисперсия прогноза, определяемая чаще всего излишней чувствительностью модели к небольшим изменениям в распределении обучающих данных;
Неадекватность - свойство модели, которая не отражает основных закономерностей в генеральной совокупности и в силу случайного характера обучающей выборки;
Переусложнение, или переобучение модели (англ. overfitting), которое «...так же вредно, как и ее недоусложнение» (Ивахненко, 1982). Переобученная модель дает хорошие результаты подгонки на обучающей совокупности, но часто ведет себя непредсказуемо при прогнозировании на "свежих" данных, не участвовавших в построении модели. Построение модели оптимальной сложности, в которой сбалансирована точность и устойчивость прогнозных значений, является фундаментальной задачей совершенствования методов моделирования. Тема переобучения моделей уже поднималась нами ранее, однако, безусловно, заслуживает дополнительного обсуждения.

Данные, используемые в примерах

В одном из предыдущих сообщений в качестве простого примера мы использовали данные по электрическому сопротивлению (Ом) мякоти фруктов киви в зависимости от процентного содержания в ней сока. Одним из наиболее эффективных способов выявления характера связи между подлежащей предсказанию переменной и предикторами является графическое изображение этой зависимости (глаз человека - все еще наиболее совершенный инструмент экспертизы нелинейности связи!).

Здесь и далее для построения графиков воспользуемся одним из наиболее популярных графических пакетов для R – ggplot2, который позволяет строить как всевозможные простые диаграммы, так и гораздо более сложные графики, такие как двухмерные диаграммы распределения плотности (см., например, Мастицкий, 2016). Легко показать линию простой регрессии с ее доверительными интервалами, которая в рассматриваемом нами примере не вполне адекватна имеющимся данным:

library(ggplot2)
library(DAAG)
data("fruitohms")
ggplot(fruitohms, aes(x = juice, y = ohms)) + geom_point() +
       stat_smooth(method = "lm") + 
       xlab("Содержание сока, %") + 
       ylab("Сопротивление, Ом")

Сглаживание данных c использованием базовой R-функции loess() лучше подходит для описания этой зависимости (в качестве примера, вместо обычной диаграммы рассеяния данные будут изображены в виде "сотовой диаграммы", где наблюдения распределены по ячейкам и число наблюдений в каждой ячейке представлено цветом соответствующего оттенка):

library(hexbin)
ggplot(fruitohms, aes(x = juice, y = ohms)) +
       geom_hex(binwidth = c(3, 500)) +
       geom_smooth(color = "red", se = FALSE) + 
       xlab("Содержание сока, %") + 
       ylab("Сопротивление, Ом")

Оценка адекватности (валидности) модели

При выполнении поиска оптимальной модели для тех или иных данных полезно иметь набор "эталонов" для сравнения. Если мы зададимся каким-нибудь подходящим критерием качества, то сможем оценить, насколько тестируемая модель отличается по этому критерию от эталона. Такими эталонными моделями могут быть следующие:

Нулевая модель (англ. null model) – модель, которая не содержит никаких предикторов. В случае с количественной зависимой переменной нулевая модель описывает ни что иное, как среднее значение этой переменной в имеющейся обучающей выборке. В случае же с качественной зависимой переменной (задача классификации) нулевая модель предсказывает наиболее часто встречающуюся категорию. Если тестируемая нами предсказательная модель не превосходит по своей эффективности нулевую, то результат моделирования можно трактовать как неудачу.
Насыщенная модель (англ. saturated model) – модель, включающая все имеющиеся предикторы. Если качество тестируемой модели по тому или иному критерию значимо превышает таковое у нулевой модели и приближается к качеству максимально насыщенной модели, то процесс отбора оптимальной модели можно считать завершенным.
Модель с одним наиболее информативным предиктором. Если в процессе отбора более сложные модели по своей эффективности не превосходят модель с одним наиболее информативным предиктором, то включение дополнительных предикторов вряд ли имеет смысл.

Процедуры оценки адекватности классической модели регрессии (обозначим ее как f) основаны на анализе остатков (англ. residuals), т.е. разностей между прогнозируемыми f(x[i,]) и наблюдаемыми y[i] значениями, где x - матрица независимых переменных (предикторов). Базовыми критериями являются сумма квадратов остатков RSS (англ. residual sum of squares), корень из среднеквадратичной ошибки RMSE (англ. root mean square error) и стандартное отклонение остатков RSE (англ. residual standard error): В рассматриваемом примере у нас есть лишь одна зависимая переменная (электропроводность) и один предиктор (содержание сока в мякоти фруктов). Для удобства обозначим их как y и x соответственно:

y <- fruitohms$juice
x <- fruitohms$ohms

#  Строим модель с одним предиктором:
n <- dim(as.matrix(x))[1]
m <- dim(as.matrix(x))[2] 
M_reg <- lm(y ~ x)
pred <- predict(M_reg)
RSS <- sum((y - pred) * (y - pred))
RMSE <- sqrt(RSS/n)
RSE <- sqrt(RSS/(n - m - 1)) 
c(RSS, RMSE, RSE)
[1] 12548.645657     9.901328     9.979600

Другими важными критериями качества аппроксимации данных являются коэффициент детерминации (Rsquared), дисперсионное отношение Фишера (F) и соответствующее ему р-значение:

Rsquared <- 1 - RSS/sum((mean(y) - y)^2)
F <- (sum((mean(y) - pred)^2)/m)/(RSS/(n - m - 1))
p <- pf(q = F, df1 = m, df2 = (n - m - 1), lower.tail = FALSE)
c(Rsquared, F, p)
[1] 6.387089e-01 2.227492e+02 1.234110e-29

Разумеется, все эти величины можно было бы получить и с помощью команды summary(M_reg), но нам показалось интересным показать всю "кухню" их расчетов. Отметим, что F-критерий интерпретируется как мера превышения качества построенной модели над качеством нулевой модели, которая не содержит объясняющих переменных.

summary(M_reg)
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 66.1364158  2.2556380   29.32   <2e-16 ***
x           -0.0071064  0.0004762  -14.93   <2e-16 ***
---
Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1  1 
Residual standard error: 9.98 on 126 degrees of freedom
Multiple R-squared: 0.6387,  Adjusted R-squared: 0.6358 
F-statistic: 222.7 on 1 and 126 DF,  p-value: < 2.2e-16

Чтобы выяснить, насколько объясненная нашей моделью сумма квадратов остатков превысила общую сумму квадратов остатков (т.е. сумму квадратов остатков у нулевой модели) можно выполнить дисперсионный анализ (см. также здесь):

anova(M_reg)
Analysis of Variance Table
Response: y
           Df Sum Sq Mean Sq F value    Pr(>F)    
x           1  22184 22184.1  222.75 < 2.2e-16 ***
Residuals 126  12549    99.6

В случае обобщенных линейных моделей (Generalized Linear Models, GLM) вместо минимизации суммы квадратов отклонений выполняется поиск экстремума логарифма функции максимального правдоподобия (англ. maximum likelihood), вид которой зависит от характера распределения данных. В общем случае значение функции правдоподобия численно равно вероятности того, что модель правильно предсказывает любое предъявленное ей наблюдение из заданной выборки. Логарифм функции правдоподобия (LL) всегда будет отрицательным, и поскольку $\log(1) = 0$, то для наиболее оптимальной модели, прогнозирующей значение отклика с наименьшей ошибкой, значение LL будет стремиться к 0. Часто для оценки расхождений между наблюдаемыми и прогнозируемыми значениями вместо LL используют девианс (англ. deviance; в русскоязычной литературе используются также термины "девиация" и "аномалия"), который определен как $D = -2(LL - S)$, где S – правдоподобие модели с минимально возможным уровнем ошибки (такую ошибку называют также неустранимой, или (в задачах классификации) байесовской). Обычно принимают S = 0, поскольку истинно оптимальная модель, как правило, выполняет верные предсказания.

Чем меньше остаточный девианс D некоторой модели, тем лучше эта модель. Аналогично можно рассчитать логарифм правдоподобия и девианс для нулевой модели (ниже обозначен как D.null). Тогда адекватность тестируемой модели можно описать при помощи псевдо-коэффициента детерминации (Prsquar), который рассчитывается с учетом того, во сколько раз D превышает D.null. Подобный анализ является обобщением техники дисперсионного анализа. При этом статистическую значимость разности девиансов можно оценить по критерию $\chi^2$:

M_glm <- glm(y ~ x)
lgLik <- logLik(M_glm)
D.null <- M_glm$null.deviance 
D <- M_glm$deviance
df <- with(M_glm, df.null - df.residual)
p <- pchisq(D.null - D, df, lower.tail = FALSE)
PRsquar = 1 - D/D.null
c(lgLik, D , D.null, p, PRsquar)
[1]  -475.08574  12548.64565 34732.77929 0.00000 0.63871

То же можно получить и с использованием функции summary():

summary(M_glm)
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 66.1364158  2.2556380   29.32   <2e-16 ***
x           -0.0071064  0.0004762  -14.93   <2e-16 ***
---
Signif. codes:  0 С***Т 0.001 С**Т 0.01 С*Т 0.05 С.Т 0.1 С Т 1 
(Dispersion parameter for gaussian family 99.59243)
    Null deviance: 34733  on 127  degrees of freedom
Residual deviance: 12549  on 126  degrees of freedom
AIC: 956.17

with(M_glm, null.deviance - deviance)
[1] 22184.13

anova(M_glm, glm(x ~ 1) , test = "Chisq")
Analysis of Deviance Table
     Df Deviance Resid. Df Resid. Dev P(>|Chi|)    
NULL                   127      34733              
x     1    22184       126      12549 < 2.2e-16 ***

Нетрудно заметить, что классическая (т.е. рассчитываемая по методу наименьших квадратов) и обобщенная линейные модели в случае нормального распределения дают идентичные результаты.

Оптимизация параметров моделей с использованием информационных критериев

Поиск оптимальной модели часто реализуется с использованием тех или иных критериев качества, налагающих "штраф" на добавление новых параметров. Если D - выборочный остаточный девианс, k – число свободных параметров модели, а n - объем обучающей выборки, то можно определить семейство следующих информационных критериев (Burnham, Anderson, 2002):

классический критерий Акаике: AIC = D + 2 * k
байесовский критерий Шварца: BIC = D + k * ln(n)
скорректированный критерий Акаике: AICс = AIC + 2 * k * (k+1)/(n-k-1)

При построении параметрических моделей с использованием R рассчитать эти критерии можно различными способами:

k <- extractAIC(M_glm)[1]
AIC <- extractAIC(M_glm)[2]
AIC <- AIC(M_glm)
AICc <- AIC(M_glm) + 2 * k * (k + 1)/(n - k - 1)
BIC <- BIC(M_glm)
BIC <- AIC(M_glm, k = log(n))
c(AIC, AICc, BIC)
[1] 956.1715 956.2675 964.7276

Оптимальной считается модель с минимальным значением соответствующего информационного критерия. Рассмотрим пример использования информационных критериев для выбора оптимального числа параметров полиномиальной регрессии:

# Построение моделей со степенью полинома от 1 до 7:
max.poly <- 7

# Создание пустой таблицы для хранения значений AIC и BIC,
# рассчитанных для всех моделей и ее заполнение:
AIC.BIC <- data.frame(criterion = c(rep("AIC", max.poly),
                                    rep("BIC", max.poly)),
                      value = numeric(max.poly * 2),
                      degree = rep(1:max.poly, times = 2))
for(i in 1:max.poly) {
     AIC.BIC[i, 2] <- AIC(lm(y ~ poly(x, i)))
     AIC.BIC[i + max.poly, 2] <- BIC(lm(y ~ poly(x, i)))
}

# График AIC и BIC для разных степеней полинома:
qplot(degree, value, data = AIC.BIC, 
      geom = "line", linetype = criterion) +
xlab("Степень полинома") + ylab("Значение критерия")

Минимальные значения информационных критериев приводят к выводу о том, что оптимальным является полином 5-й степени.

Использование алгоритмов ресэмплинга для тестирования и оптимизации параметров моделей

Минимизация ошибок на обучающем множестве, которое не бывает ни идеальным, ни бесконечно большим, неизбежно приводит к моделям, смещенным относительно истинной функции процесса, порождающего наблюдаемые данные. Преодолеть это смещение можно только с использованием "принципа внешнего дополнения", т.е. блоков "свежих" данных, не участвовавших в построении модели.

Перекрестная проверка или "кросс-проверка" (англ. cross-validation, CV) является общей процедурой эмпирического оценивания параметров моделей по выборочным данным. В ходе многократного случайного разбиения исходной выборки на обучающее и проверочное подмножества оценка параметров модели выполняется только на обучающей выборке, тогда как оценка погрешности прогноза выполняется для наблюдений из проверочной совокупности.

Одной из часто применяемых является методика r×k-кратной перекрестной проверки (англ. r×k-fold cross-validation), когда исходная выборка случайным образом r раз разбивается на k блоков (англ. folds) равной (или почти равной) длины. При реализации каждой повторности r (англ. replication) один из блоков по очереди становится проверочной последовательностью, а все остальные блоки – обучающей выборкой. При этом генерируется r×n модельных значений отклика (n – общее количество наблюдений), по которым оценивается ошибка перекрестной проверки.

Выполнить разбиение исходной выборки y на k блоков можно с использованием следующих функций из пакета caret (Kuhn, Johnson, 2013), о которых мы упоминали ранее при рассмотрении функции createDataPartition():

createFolds (y, k = 10)
createMultiFolds (y, k = 10, times = r)

Частным случаем является "скользящий контроль", или перекрестная проверка с последовательным исключением одного наблюдения (англ. leave-one-out CV), т.е. k = n. При этом строится n моделей по (n – 1) выборочным значениям, а исключенное наблюдение каждый раз используется для расчета ошибки прогноза. В. Н. Вапник (1984) теоретически обосновал применение скользящего контроля и показал, что если исходные выборки независимы, то средняя ошибка перекрестной проверки дает несмещенную оценку ошибки модели. Это выгодно отличает ее от средней ошибки на обучающей выборке, которая при переусложнении модели может оказаться оптимистично заниженной.

Выше мы рассмотрели пример выбора оптимального числа параметров полиномиальной регрессии с использованием информационных критериев. Попробуем достичь той же цели, выполнив перекрестную проверку моделей в режиме "leave-one-out" с использованием самостоятельно оформленной функции:

# Функция для выполнения скользящего контроля для модели y ~ poly(x, degree):
crossvalidate <- function(x, y, degree) {
   
   preds <- numeric(length(x))

   for (i in 1:length(x)) {
      x.in <- x[-i]
      x.out <- x[i]
      y.in <- y[-i]
      y.out <- x[i]
      m <- lm(y.in ~ poly(x.in, degree = degree) )
      new <- data.frame(x.in = seq(-3, 3, by = 0.1))
      preds[i] <- predict(m, newdata = data.frame(x.in = x.out))
   }

  # Тестовая статистика - сумма квадратов отклонений:
  return(sum((y - preds)^2))
}

# Заполнение таблицы результатами перекрестной проверки 
# и сохранение квадрата ошибки в таблице "a":
a <- data.frame(cross = numeric(max.poly))
for (i in 1:max.poly)
{
  a[i, 1] <- crossvalidate(x, y, degree = i)
}

# График суммы квадратов ошибки при перекрестной проверке:
qplot(1:max.poly, cross, data = a, geom = c("line")) +
      xlab("Степень полинома") + ylab("Квадратичная ошибка")

Из представленного выше графика следует, что в нашем конкретном случае оба метода в качестве оптимальной выбрали модель с полиномом 5-й степени. Проверим, имеется ли значимое превышение ее качества над моделью с одним предиктором:

M_poly <- glm(y ~ poly(x, 5))
summary(M_poly)
Coefficients:
(Intercept)  poly(x, 5)1  poly(x, 5)2  poly(x, 5)3  poly(x, 5)4  poly(x, 5)5  
      35.15      -148.94        -5.54        51.08       -13.90       -23.53  
Degrees of Freedom: 127 Total (i.e. Null);  122 Residual
Null Deviance:      34730 
Residual Deviance: 9162         AIC: 923.9 

anova(M_glm, M_poly, test = "Chisq")
Analysis of Deviance Table
Model 1: y ~ x
Model 2: y ~ poly(x, 5)
  Resid. Df Resid. Dev Df Deviance P(>|Chi|)    
1       126      12549                          
2       122       9162  4   3386.6 3.798e-09 ***

Как видим, имеет место статистически значимое снижение ошибки модели.

В R перекрестную проверку модели можно выполнить с использованием не только самостоятельно оформленных процедур, но и функций из нескольких общедоступных пакетов:

CVlm(df, form.lm, m = 3) из пакета DAAG, где df – исходная таблица с данными, form.lm – формула линейной модели, m – число блоков (сюда подставляется значение k);
cv.glm(df, glmfit, K = n) из пакета boot, где df – исходная таблица с данными, glmfit – объект обобщенной линейной модели класса glm, K – число блоков (по умолчанию выполняется скользящий контроль);
cvLm(lmfit, folds = cvFolds(n, K, R), cost) из пакета cvTools, где lmfit – объект обобщенной линейной модели, K – число блоков, R – число репликаций, cost – наименование одного из пяти разновидностей выводимых критериев перекрестной проверки;
train(form, df, method, trainControl, ...) из пакета caret, где df – исходная таблица с данными, form, method – формула и метод построения модели, trainControl – объект, в котором прописываются все необходимые параметры перекрестной проверки.

Использование функции train(...) будет предметом наших дальнейших сообщений.

Литературные источники

Вапник В. Н. (1984) Алгоритмы и программы восстановления зависимостей. М.: Наука, 816 с.
Ивахненко А. Г. (1982) Индуктивный метод самоорганизации моделей сложных систем. Киев: Наук. думка, 296 с.
Мастицкий С. Э. (2016) Визуализация данных с помощью ggplot2. М.: ДМК Пресс, 222 с.
Burnham K. P., Anderson D. R. (2002) Model selection and multimodel inference: a practical information-theoretic approach. N.Y.: Springer-Verlag, 496 p.
Kuhn M., Johnson K. (2013) Applied Predictive Modeling. Springer, 600 p.

Обзор R-пакетов за ноябрь 2016 г.

2016-12-10T12:49:00.001+00:00

Автор рисунка: Xiao Nan

По состоянию на 10 декабря 2016 г. на сайте CRAN было опубликовано около 12 тысяч пакетов для R. Такое разнообразие приложений, наряду с очень активным сообществом пользователей, является уникальным преимуществом R, по сравнению с другими языками и системами, предназначенными для анализа данных. В то же время становится все труднее находить пакеты, подходящие для решения тех или иных задач. Эта проблема особенно актуальна среди пользователей, недотостаточно владеющих английским языком. Чтобы как-то помочь таким пользователям, я буду периодически делать небольшие обзоры наиболее интересных пакетов, опубликованных или обновленных примерно за последний месяц (ниже представлен обзор для ноябрь - начало декабря 2016 г.). Безусловно, эти обзоры будут в значительной мере отражать мое собственное представление о том, что является интересным. Поэтому если я что-то упустил, сообщайте об этом в комментариях.

anytime - мощный инструмент для работы с любыми форматами дат.
Несколько утилит от Томаса Липера (Thomas Leeper) для работы с Amazon Web Services:

aws.signature: генерация подписей для различных AWS API;
aws.polly: R-клиент для работы с сервисом имитации речи Polly;
aws.sns: R-клиент для работы с сервисом отправки уведомений Simple Notification Service (SNS);
aws.sqs: R-клиент для работы с Simple Queue Service.

bayesplot - набор функций для визуализации апостериорных распределений параметров моделей и диагностики качества моделей, построенных с использованием принципов байесовской статистики. В частности, содержит функции, которые позволяют работать с моделями Stan.
bookdown - набор утилит для подготовки книг и технической документации с помощью R Markdown.
diverse - предназначен для вычисления целого ряда мер энтропии, используемых в социальных и других науках. Может оказаться особенно полезным для создания новых предикторов по имеющимся переменным при разработке предсказательных моделей.
fastcluster - эффективная (с точки зрения скорости вычислений) реализация иерархического кластерного анализа. Пакет одновременно содержит соответствующие функции как для R, так и для Python (см. руководство пользователя).
Radviz (от англ. radial coordinate visualisation) - позволяет выполнять проецирование многомерных данных на двухмерную плоскость, что, в частности, будет полезным при визуализации результатов кластерного анализа (см. руководство пользователя).
rex - набор функций для работы с регулярными выражениями. Стандартные возможности R для работы с регулярными выражениями трудно назвать "дружественными" для пользователя. Пакет rex решает эту проблему очень удачно (см. пример).
tsoutliers - функции для обнаружения выбросов при работе с временными рядами (по методу Chen & Liu, 1993).
Strategy - пакет для разработки и оценки эффективности торговых стратегий. Имеется возможность использовать как готовые стратегии, так и добавлять свои собственные.

Сергей Мастицкий

Новые книги по R на русском языке

2016-12-04T10:04:00.000+00:00

Осенью этого года в московском издательстве "ДМК Пресс" вышли две новые книги по R на русском языке. Краткое описание этих книг приводится ниже.

Мастицкий С. Э. (2016) Визуализация данных с помощью ggplot2. - М.: ДМК Пресс, 222 с.

Визуализация данных играет важную роль на всех этапах статистического анализа – от первичного ознакомления со свойствами данных до диагностики качества построенных моделей и представления полученных результатов. Из всего разнообразия статистических программ выделяется R - интенсивно развивающаяся и свободно распространяемая система статистических вычислений, в которой реализовано множество классических и современных методов анализа данных. Программные реализации алгоритмов, входящих в базовую версию R, проверены на практике не одним поколением пользователей и ученых. Кроме того, пользователи R постоянно разрабатывают многочисленные дополнения (т.н. «пакеты») для этой системы. Настоящая книга посвящена ggplot2 - одному из таких пакетов, который значительно расширяет и без того богатые базовые графические возможности R. В 2015 г. ggplot2 был установлен более миллиона раз. Такая популярность этого пакета обусловлена несколькими причинами, среди которых можно отметить эстетическую привлекательность и пригодное для публикации качество получаемых с его помощью графиков, возможность создавать пользовательские типы диаграмм, а также большой набор опций для тонкой настройки внешнего вида графиков. В этой книге описаны основы работы с ggplot2 и приведены многочисленные примеры кода, которые читатели легко могут модифицировать для собственных нужд. Книга окажется полезной для всех пользователей R, желающих освоить новый мощный инструмент анализа данных.

Храмов Д. А. (2017) Сбор данных в Интернете на языке R. - М.: ДМК Пресс, 283 с.

Всё, что регистрирует человек и созданные им машины, может считаться данными. Фиксируя новое и переводя архивы в цифровую форму, мы с каждым днём производим всё больше данных. Часть из них находится в специальных хранилищах, готовые к использованию. Но гораздо чаще случается так, что данные разбросаны по всемирной сети на многочисленных страницах онлайновых магазинов, заметках в социальных сетях, логах серверов и т. п. Прежде чем начать работать с такими данными, их необходимо собрать и сохранить в пригодном для анализа виде. Решению этих вопросов и посвящена данная книга. Основной материал книги разделён на две части. В первой части дано краткое введение в R - описание среды разработки, языка и основных пакетов-расширений. Вторая часть посвящена непосредственно сбору данных: работе с открытыми данными, извлечению данных из веб-страниц и из социальных сетей. Также рассмотрены необходимые технические вопросы: протокол HTTP, функции импорта данных различных форматов и регулярные выражения. Завершается рассказ созданием карт на основе собранных данных. Издание предназначено специалистам по анализу данных, а также программистам, интересующихся сбором данных в Интернете.

Распространение информации об этих книгах приветствуется!

Сергей Мастицкий

Перезагрузка

2016-11-19T05:39:00.000+00:00

В июне этого года я объявил о прекращении ведения блога "R: Анализ и визуализация данных" в связи с хронической нехваткой у меня свободного времени. Оказалось, однако, что многие читатели хотели бы продолжения этого проекта. Сегодня я рад сообщить о его "возвращении к жизни", хотя и в несколько измененном виде. У сайта появился новый редактор - доктор биологических наук, профессор Андрей Викторович Коросов, который любезно согласился взять на себя основной труд по подготовке новых публикаций. Я также буду участвовать в написании новых материалов, однако не так активно, как в прошлом. Вместе с Андреем Викторовичем мы постараемся привлечь и других авторов, желающих поделиться своим опытом работы с R. Если бы вы хотели бы стать одним из таких авторов, то свяжитесь с нами, пожалуйста, по адресу электронной почты, приведенному в разделе "О проекте".

Сергей Мастицкий