Глубокие нейросети (часть II). Выбор предикаторов

Глубокие нейросети (часть II). Разработка

1. Разработка признаков

Разработка признаков — это наука (и искусство) извлечения дополнительной информации из имеющихся данных. Здесь наша задача — не добавить новые данные, а сделать более полезными и “говорящими” те, которые у нас уже есть. Новые возможности позволяют получить новые признаки выборки данных. Эти признаки позволяют более точно характеризовать и разделять обучающие данные, и за счет этого дают улучшенную точность модели.

Технологически процесс может быть разделен на два этапа:

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

 

1.1. Трансформация признаков

 

1.1.1. Преобразование

 

Что такое преобразование переменной?

В моделировании данных преобразование заключается в замене переменной функцией. Например, замена переменной x на квадратный / кубический корень или логарифм x — это трансформация (или преобразование). Иными словами, трансформация/преобразование — это процесс, который изменяет распределение или связь переменной с другими.

Перечислим ситуации, когда преобразование переменной полезно.

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

 

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

 

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

 

    • Когда нужно преобразовать непрерывную переменную в дискретную. Метод такого преобразования — дискретизация.

 

Каковы общие методы преобразования переменной?

 

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

    1. Логарифмирование: это общий метод преобразования, используемый для изменения формы распределения переменной на распределительном участке. Обычно его применяют для сокращения правой асимметрии. К нулю или отрицательным значениям логарифмирование не может быть применено.

 

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

 

  1. Дискретизация/Бининг: используется для категоризации переменных. Выполняется на исходных значениях, процентилях или частотах. Решение о методе категоризации основывается на том, какова природа данных. Мы можем проводить совместные сегментирования взаимозависимых переменных.

Любая трансформация данных приводит к изменению распределения переменных. Рассмотрим это на примерах двух методов преобразования.

Две проблемы нашего исходного набора данных — выбросы и правая асимметрия. Удаление выбросов мы уже рассмотрели. Сейчас сначала попробуем удалить/уменьшить асимметрию, а после этого уберем выбросы.

Метод 1.

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

evalq({x.ln <- apply(x, 2, function(x) log2(x + 1))
       sk.ln <- skewness(x.ln)},
      env)
 > env$sk.ln
               ftlm      stlm      rbci      pcci   v.fatl
Skewness -0.2715663 -2.660613 -4.484301 0.4267873 1.253008
          v.satl   v.rftl     v.rstl    v.ftlm     v.stlm
Skewness 1.83489 2.065224 -0.0343451 -15.62414 0.01529019
            v.pcci
Skewness 0.1811206

 

Три переменных — stlm, rbci и v.ftlm — получили сильную левую асимметрию, три переменных — v.fatl, v.satl и v.rftl — остались с сильным правым скосом, у остальных переменных асимметрия выровнялась. Уберем и импутируем выбросы из этого набора и посмотрим, какая асимметрия и распределение переменных у нас получится:

evalq({
  foreach(i = 1:ncol(x.ln), .combine = "cbind") %do% {
    remove_outliers(x.ln[ ,i])
  } -> x.ln.out
  colnames(x.ln.out) <- colnames(x.ln)
  },  
env)
evalq({
  foreach(i = 1:ncol(x.ln), .combine = "cbind") %do% {
    capping_outliers(x.ln[ ,i])
  } -> x.ln.cap
  colnames(x.ln.cap) <- colnames(x.ln)
},  
env)
evalq({
  sk.ln.out <- skewness(x.ln.out) 
  sk.ln.cap <- skewness(x.ln.cap)
}, 
env)
> env$sk.ln.out
              ftlm       stlm       rbci        pcci
Skewness -0.119055 -0.3549119 -0.1099921 -0.01476384
              v.fatl      v.satl      v.rftl     v.rstl
Skewness -0.02896319 -0.03634833 -0.06259749 -0.2120127
              v.ftlm      v.stlm      v.pcci
Skewness -0.05819699 -0.01661317 -0.05420077
> env$sk.ln.cap
               ftlm       stlm       rbci        pcci
Skewness -0.1814781 -0.4582045 -0.1658855 -0.02849945
              v.fatl      v.satl     v.rftl     v.rstl
Skewness -0.04336238 -0.04400781 -0.0692754 -0.2269408
              v.ftlm      v.stlm      v.pcci
Skewness -0.06184128 -0.02856397 -0.06258243

 

В обоих наборах (x.out и x.cap) данные почти симметричны. Распределение продемонстрировано на рисунках ниже.

par(mfrow = c(2,2))
boxplot(env$x.ln, 
              main = "x.ln with outliers",
              xlab = "")
boxplot(env$x.ln.out, 
              main = "x.ln.out without outliers",
              xlab = "")
boxplot(env$x.ln.cap, 
              main = "x.ln.cap with imputed outliers",
              xlab = "")
par(mfrow = c(1,1))

 

x_ln__1x_ln_out__1

Рис.1. Логтрансформированные данные с выбросами и без

x_ln_cap__1

Рис.2. Логтрансформированные данные с импутированными выбросами

Результаты похожи на предыдущее преобразование, за одним исключением — диапазон изменения переменных стал шире.

Преобразуем x.ln.cap датафрейм и посмотрим вариацию и ковариацию набора:

 evalq(x.ln.cap %>% tbl_df() %>% 
        cbind(Data = dataSetClean$Data, .,
              Class = dataSetClean$Class) -> 
        dataSetLnCap, 
      env)

Построим графики:

require(GGally)
evalq(ggpairs(dataSetLnCap, columns = 2:7, 
              mapping = aes(color = Class),
              title = "PredLnCap1"), 
      env)
evalq(ggpairs(dataSetLnCap, columns = 8:13, 
              mapping = aes(color = Class),
              title = "PredLnCap2"), 
      env)

 

LnCap1__1

Рис.3. Параметры логтрансформированных данных, часть 1

LnCap2__1

Рис. 4. Параметры логтрансформированных данных, часть 2

Метод 2.

Трансформируем данные с помощью функции sin(2*pi*x), уберем и импутируем выбросы и посмотрим на асимметрию, распределение выбросов и ковариацию трансформированных переменных на графиках.

evalq({x.sin <- apply(x, 2, function(x) sin(2*pi*x))
      sk.sin <- skewness(x.sin)
      },
env)
#----------
evalq({
  foreach(i = 1:ncol(x.sin), .combine = "cbind") %do% {
    remove_outliers(x.sin[ ,i])
  } -> x.sin.out
  colnames(x.sin.out) <- colnames(x.sin)
},  
env)
#-----------------
evalq({
  foreach(i = 1:ncol(x.sin), .combine = "cbind") %do% {
    capping_outliers(x.sin[ ,i])
  } -> x.sin.cap
  colnames(x.sin.cap) <- colnames(x.sin)
},  
env)
#-----------
evalq({
  sk.sin.out <- skewness(x.sin.out) 
  sk.sin.cap <- skewness(x.sin.cap)
}, 
env)

Какова асимметрия в этих трансформированных наборах?

env$sk.sin
                ftlm        stlm        rbci         pcci
Skewness -0.02536085 -0.04234074 -0.00587189 0.0009679463
             v.fatl    v.satl     v.rftl      v.rstl
Skewness 0.03280465 0.5217757 0.05611136 -0.02825112
             v.ftlm     v.stlm     v.pcci
Skewness 0.04923953 -0.2123434 0.01738377
> env$sk.sin.out
                ftlm        stlm        rbci       pcci
Skewness -0.02536085 -0.04234074 -0.00587189 0.03532892
             v.fatl      v.satl      v.rftl      v.rstl
Skewness 0.00360966 -0.02380975 -0.05336561 -0.02825112
               v.ftlm     v.stlm       v.pcci
Skewness 0.0009366441 0.01835948 0.0008843329
> env$sk.sin.cap
                ftlm        stlm        rbci       pcci
Skewness -0.02536085 -0.04234074 -0.00587189 0.03283132
              v.fatl      v.satl      v.rftl      v.rstl
Skewness 0.007588308 -0.02424707 -0.04106469 -0.02825112
              v.ftlm      v.stlm      v.pcci
Skewness 0.007003051 0.009237835 0.002101687

Как видим, при этой трансформации все наборы симметричны. Посмотрим теперь, как эти наборы выглядят:

par(mfrow = c(2, 2))
boxplot(env$x.sin, main = "x.sin with outlier")
abline(h = 0, col = 2)
boxplot(env$x.sin.out, main = "x.sin.out without outlier")
abline(h = 0, col = 2)
boxplot(env$x.sin.cap, main = "x.sin.cap with capping outlier")
abline(h = 0, col = 2)
par(mfrow = c(1, 1))

xSin__1

Рис.5. Набор данных, трансформированный функцией sin()

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

Посмотрим, как NA распределены в переменных после удаления выбросов.

require(VIM)
evalq(a <- aggr(x.sin.out), env)

SinMissAggr__1

Рис.6. Распределение NA в наборе

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

> print(env$a)

 Missings in variables:
 Variable Count
     pcci   256
   v.fatl   317
   v.satl   289
   v.rftl   406
   v.ftlm   215
   v.stlm   194
   v.pcci   201

Как распределены NA в переменных?

 par(mfrow = c(3, 4))
evalq(
  foreach(i = 1:ncol(x.sin.out)) %do% {
    barMiss(x.sin.out, pos = i, only.miss = TRUE, 
            main = "x.sin.out without outlier")
  }, env
)
par(mfrow = c(1, 1))

SinMissBar__1

Рис.7. Распределение NA в переменных

Синим цветом показаны наблюдаемые значения переменной, красным — количество NA других переменных в различных диапазонах значений текущей. Бар справа — вклад NA текущей переменной в общее количество NA всех переменных.

В заключение посмотрим вариацию и ковариацию трансформированного набора с импутированными выбросами.

#---------------
evalq(x.sin.cap %>% tbl_df() %>% 
        cbind(Data = dataSetClean$Data, .,
              Class = dataSetClean$Class) -> 
        dataSetSinCap, 
      env) 
require(GGally)
evalq(ggpairs(dataSetSinCap, columns = 2:7, 
              mapping = aes(color = Class),
              title = "dataSetSinCap1 with capping outlier "), 
      env)
evalq(ggpairs(dataSetSinCap, columns = 8:13, 
              mapping = aes(color = Class),
              title = "dataSetSinCap2 with capping outlier"), 
      env)
#---------------------------

 

SinCap1__1

Рис.8. Параметры трансформированных sin() данных, часть 1

SinCap2__3

Рис.9. Параметры трансформированных sin() данных, часть 2

 

1.1.2. Нормализация

 

Поскольку мы готовим данные для нейросети, нужно привести переменные в диапазон { -1..+1 }. Для этого используем функцию preProcess()::caret и method = “spatialSign”. Как вариант, перед нормализацией можно центрировать и шкалировать данные. Поскольку это довольно тривиальный процесс, мы не будем останавливаться на нем подробно.

Единственное, что нужно подчеркнуть: параметры нормализации получаем на тренировочном наборе, тестовый и валидационный наборы обрабатываем с ними же.

Для использования в дальнейших экспериментах разделим наш набор, полученный в предыдущих вычислениях (dataSet без удаления высокоррелированных ) на train/test/val и приведем в диапазон (-1,+1) без стандартизации.

При нормировании со стандартизацией нужно помнить, что при определении параметров нормализации (mean/median, sd/mad) мы определяем и параметры импутации выбросов. В дальнейшем они будут использованы на train/val/test . Ранее в статье мы уже написали 2 функции: prep.outlier() и treatOutlier(), которые предназначены именно для этих целей.

Последовательность выполнения операций:

  1. Определение параметров выбросов в train
  2. Удаление выбросов в train
  3. Определение параметров стандартизации в train
  4. Импутируем выбросы в train/val/test
  5. Нормализация train/val/test.

Этот вариант в статье мы рассматривать не будем, вы можете ознакомиться с ним самостоятельно.

Разделяем набор на train/val/test:

 evalq(
{
  train = 1:2000
  val = 2001:3000
  test = 3001:4000
  DT <- list()
  list(clean = data.frame(dataSet) %>% na.omit(), 
       train = clean[train, ], 
       val = clean[val, ], 
       test = clean[test, ]) -> DT
}, env)

Определим параметры нормализации для набора train и нормализуем наборы train/test/val:

 require(foreach)
evalq(
{
 preProcess(DT$train, method = "spatialSign") -> preproc 
 list(train = predict(preproc, DT$train), 
        val = predict(preproc, DT$val),
        test = predict(preproc, DT$test)
       ) -> DTn
}, 
env)

Посмотрим суммарную статистику набора train:

 > table.Stats(env$DTn$train %>% tk_xts())
Using column `Data` for date_var.
                     ftlm      stlm      rbci      pcci
Observations    2000.0000 2000.0000 2000.0000 2000.0000
NAs                0.0000    0.0000    0.0000    0.0000
Minimum           -0.5909   -0.7624   -0.6114   -0.8086
Quartile 1        -0.2054   -0.2157   -0.2203   -0.2110
Median             0.0145    0.0246    0.0147    0.0068
Arithmetic Mean    0.0070    0.0190    0.0085    0.0028
Geometric Mean    -0.0316   -0.0396   -0.0332   -0.0438
Quartile 3         0.2139    0.2462    0.2341    0.2277
Maximum            0.6314    0.8047    0.7573    0.7539
SE Mean            0.0060    0.0073    0.0063    0.0065
LCL Mean (0.95)   -0.0047    0.0047   -0.0037   -0.0100
UCL Mean (0.95)    0.0188    0.0333    0.0208    0.0155
Variance           0.0719    0.1058    0.0784    0.0848
Stdev              0.2682    0.3252    0.2800    0.2912
Skewness          -0.0762   -0.0221   -0.0169   -0.0272
Kurtosis          -0.8759   -0.6688   -0.8782   -0.7090
                   v.fatl    v.satl    v.rftl    v.rstl
Observations    2000.0000 2000.0000 2000.0000 2000.0000
NAs                0.0000    0.0000    0.0000    0.0000
Minimum           -0.5160   -0.5943   -0.6037   -0.7591
Quartile 1        -0.2134   -0.2195   -0.1988   -0.2321
Median             0.0015    0.0301    0.0230    0.0277
Arithmetic Mean    0.0032    0.0151    0.0118    0.0177
Geometric Mean    -0.0323   -0.0267   -0.0289   -0.0429
Quartile 3         0.2210    0.2467    0.2233    0.2657
Maximum            0.5093    0.5893    0.6714    0.7346
SE Mean            0.0058    0.0063    0.0062    0.0074
LCL Mean (0.95)   -0.0082    0.0028   -0.0003    0.0033
UCL Mean (0.95)    0.0146    0.0274    0.0238    0.0321
Variance           0.0675    0.0783    0.0757    0.1083
Stdev              0.2599    0.2798    0.2751    0.3291
Skewness          -0.0119   -0.0956   -0.0648   -0.0562
Kurtosis          -1.0788   -1.0359   -0.7957   -0.7275
                   v.ftlm    v.stlm    v.rbci    v.pcci
Observations    2000.0000 2000.0000 2000.0000 2000.0000
NAs                0.0000    0.0000    0.0000    0.0000
Minimum           -0.5627   -0.6279   -0.5925   -0.7860
Quartile 1        -0.2215   -0.2363   -0.2245   -0.2256
Median            -0.0018    0.0092   -0.0015   -0.0054
Arithmetic Mean   -0.0037    0.0036   -0.0037    0.0013
Geometric Mean    -0.0426   -0.0411   -0.0433   -0.0537
Quartile 3         0.2165    0.2372    0.2180    0.2276
Maximum            0.5577    0.6322    0.5740    0.9051
SE Mean            0.0061    0.0065    0.0061    0.0070
LCL Mean (0.95)   -0.0155   -0.0091   -0.0157   -0.0124
UCL Mean (0.95)    0.0082    0.0163    0.0082    0.0150
Variance           0.0732    0.0836    0.0742    0.0975
Stdev              0.2706    0.2892    0.2724    0.3123
Skewness           0.0106   -0.0004   -0.0014    0.0232
Kurtosis          -1.0040   -1.0083   -1.0043   -0.4159

Из таблицы мы видим, что переменные симметричны и имеют очень близкие параметры.

Посмотрим на распределение переменных в наборах train/val/test:

 boxplot(env$DTn$train %>% 
          dplyr::select(-c(Data, Class)),
        horizontal = T, main = "Train")
abline(v = 0, col = 2)
boxplot(env$DTn$test %>% 
          dplyr::select(-c(Data, Class)),
        horizontal = T, main = "Test")
abline(v = 0, col = 2)
boxplot(env$DTn$val %>% 
          dplyr::select(-c(Data, Class)),
        horizontal = T, main = "Val")
abline(v = 0, col = 2)

norm1__1

Рис.10. Распределение переменных в наборах train/val/test после нормализации

Распределение переменных во всех наборах практически одинаково. Дополнительно нужно рассмотреть корреляцию и ковариацию переменных в наборе train:

require(GGally)
evalq(ggpairs(DTn$train, columns = 2:7, 
              mapping = aes(color = Class),
              title = "DTn$train1 "), 
      env)
evalq(ggpairs(DTn$train, columns = 8:14, 
              mapping = aes(color = Class),
              title = "DTn$train2"), 
      env)

DTnktrain1__1

Рис.11. Вариация и ковариация набора 1 train

DTnktrain2__1

Рис.12. Вариация и ковариация набора 2 train

Высококоррелированных данных у нас нет, распределение компактное и без выбросов, данные хорошо разделяемые. На первый взгляд предварительно видны две проблемных переменных — stlm и v.rstl. В дальнейшем, когда мы будем оценивать релевантность предикторов, мы уточним это предварительное мнение. Сейчас мы можем посмотреть корреляцию этих предикторов к целевой:

> funModeling::correlation_table(env$DTn$train %>% tbl_df %>%
+                    select(-Data), str_target = 'Class')
   Variable Class
1     Class  1.00
2    v.fatl  0.38
3      ftlm  0.34
4      rbci  0.28
5    v.rbci  0.28
6    v.satl  0.27
7      pcci  0.24
8    v.ftlm  0.22
9    v.stlm  0.22
10   v.rftl  0.18
11   v.pcci  0.08
12     stlm  0.03
13   v.rstl -0.01

Видим, что указанные переменные находятся внизу таблицы с очень низкими показателями. Под вопросом также релевантность переменной v.pcci. Посмотрим ближе переменную v.fatl во всех наборах train/val/test:

require(ggvis)
evalq(
  DTn$train %>% ggvis(~v.fatl, fill = ~Class) %>% 
  group_by(Class) %>%  layer_densities() %>% 
  add_legend("fill", title = "DTn$train$v.fatl"),
  env)
evalq(
  DTn$val %>% ggvis(~v.fatl, fill = ~Class) %>% 
    group_by(Class) %>%  layer_densities() %>% 
    add_legend("fill", title = "DTn$val$v.fatl"),
  env)
evalq(
  DTn$test %>% ggvis(~v.fatl, fill = ~Class) %>% 
    group_by(Class) %>%  layer_densities() %>% 
    add_legend("fill", title = "DTn$test$v.fatl"),
  env)

Trainfvfatl__1

Рис.13. Распределение переменной v.fatl в наборе train после нормализации

Validqvfatl__1

Рис.14. Распределение переменной v.fatl в наборе valid после нормализации

Test5vfatl__1

Рис.15. Распределение переменной v.fatl в наборе test после нормализации

Проведенный анализ показывает, что нормализация часто дает хорошее распределение предикторов без выбросов и высококоррелированных данных. Во многом это зависит от характера сырых данных.

 

1.1.3. Дискретизация

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

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

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

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

Discretization. Этот пакет представляет собой набор алгоритмов дискретизации с учителем. Он также может быть сгруппирован в терминах «сверху вниз» или «снизу вверх», реализующих алгоритмы дискретизации. Рассмотрим некоторые из них на примере нашего набора dataSet.

Проведем очистку набора (без удаления высококоррелированных), разделим на train/val/test в соотношении 2000/1000/1000.

require(discretization)
require(caret)
require(pipeR)
evalq(
  {
    dataSet %>%
    preProcess(.,
               method = c("zv", "nzv", "conditionalX")) %>%
    predict(., dataSet) %>%
    na.omit -> dataSetClean
    train = 1:2000
    val = 2001:3000
    test = 3001:4000
    DT <- list()
    list(train = dataSetClean[train, ], 
         val = dataSetClean[val, ], 
         test = dataSetClean[test, ]) -> DT
  }, 
  env)

Будем использовать функцию mdlp()::discretization, которая описывает дискретизацию с использованием принципа минимальной длины описания. Эта функция дискретизирует непрерывные атрибуты матрицы данных с использованием критерия энтропии с минимальной длиной описания в качестве правила остановки.

evalq(
  pipeline({
    DT$train
    select(-Data)
    as.data.frame()
    mdlp()}) -> mdlp.train, envir = env)

Функция возвращает лист с двумя слотами: это cutp — датафрейм с точками разделения для каждой переменной и Disc.data — датафрейм с размеченными переменными.

> env$mdlp.train%>%str()
List of 2
 $ cutp     :List of 12
  ..$ : num [1:2] -0.0534 0.0278
  ..$ : chr "All"
  ..$ : num -0.0166
  ..$ : num [1:2] -0.0205 0.0493
  ..$ : num [1:3] -0.0519 -0.0055 0.019
  ..$ : num 0.000865
  ..$ : num -0.00909
  ..$ : chr "All"
  ..$ : num 0.0176
  ..$ : num [1:2] -0.011 0.0257
  ..$ : num [1:3] -0.03612 0.00385 0.03988
  ..$ : chr "All"
 $ Disc.data:'data.frame':      2000 obs. of  13 variables:
  ..$ ftlm  : int [1:2000] 3 3 3 3 3 2 1 1 1 1 ...
  ..$ stlm  : int [1:2000] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ rbci  : int [1:2000] 2 2 2 2 2 2 1 1 1 1 ...
  ..$ pcci  : int [1:2000] 2 2 1 2 2 1 1 2 3 2 ...
  ..$ v.fatl: int [1:2000] 4 4 3 4 3 1 1 2 3 2 ...
  ..$ v.satl: int [1:2000] 1 1 1 2 2 1 1 1 1 1 ...
  ..$ v.rftl: int [1:2000] 1 2 2 2 2 2 2 2 1 1 ...
  ..$ v.rstl: int [1:2000] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ v.ftlm: int [1:2000] 2 2 1 1 1 1 1 1 2 1 ...
  ..$ v.stlm: int [1:2000] 1 1 1 2 2 1 1 1 1 1 ...
  ..$ v.rbci: int [1:2000] 4 4 3 3 2 1 1 2 3 2 ...
  ..$ v.pcci: int [1:2000] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ Class : Factor w/ 2 levels "-1","1": 2 2 2 2 2 1 1 1 1 1 ...

Что мы можем увидеть в первом слоте?

У нас есть три неразмеченных переменных, значения которых никак не связаны с целевой. Это 2,8 и 12 (stlm, v.rstl, v.pcci). Их можно удалить без потери качества набора. Кстати, эти же переменные мы определили как нерелевантные и в предыдущем разделе.

4 переменных разделены на два класса, 3 переменных — на 3 класса и две переменных — на 4 класса.

Разметим наши наборы val/test на сегменты, используя точки разделения, полученные на train. Для этого уберем из набора train переменные, которые остались неразмеченными, и сохраним набор в датафрейм train.d. Затем, используя функцию findInterval(), разметим набор test/val точками разделения, полученными ранее.

evalq(
  {
    mdlp.train$cutp %>% 
    lapply(., function(x) is.numeric(x)) %>%
    unlist -> idx   # bool
    #----train-----------------
    mdlp.train$Disc.data[ ,idx] -> train.d
    #---test------------
    DT$test %>% 
      select(-c(Data, Class)) %>%
      as.data.frame() -> test.d
  
    foreach(i = 1:length(idx), .combine = 'cbind') %do% {
      if (idx[i]) {findInterval(test.d[ ,i], 
                   vec = mdlp.train$cutp,
                   rightmost.closed = FALSE, 
                   all.inside = F,
                   left.open = F)}
    } %>% as.data.frame() %>% add(1) %>%
      cbind(., DT$test$Class) -> test.d
    colnames(test.d) <- colnames(train.d)
    #-----val-----------------
    DT$val %>% 
      select(-c(Data, Class)) %>%
      as.data.frame() -> val.d
    foreach(i = 1:length(idx), .combine = 'cbind') %do% {
      if (idx[i]) {findInterval(val.d[ ,i], 
                                vec = mdlp.train$cutp,
                                rightmost.closed = FALSE, 
                                all.inside = F,
                                left.open = F)}
    } %>% as.data.frame() %>% add(1) %>%
      cbind(., DT$val$Class) -> val.d
    colnames(val.d) <- colnames(train.d)
  },env
)

Как выглядят эти наборы?

> env$train.d %>% head()
  ftlm rbci pcci v.fatl v.satl v.rftl v.ftlm v.stlm v.rbci Class
1    3    2    2      4      1      1      2      1      4     1
2    3    2    2      4      1      2      2      1      4     1
3    3    2    1      3      1      2      1      1      3     1
4    3    2    2      4      2      2      1      2      3     1
5    3    2    2      3      2      2      1      2      2     1
6    2    2    1      1      1      2      1      1      1    -1
> env$test.d %>% head()
  ftlm rbci pcci v.fatl v.satl v.rftl v.ftlm v.stlm v.rbci Class
1    1    1    1      2      1      1      1      1      2    -1
2    1    1    3      3      1      1      2      2      3    -1
3    1    1    2      2      1      1      1      2      2    -1
4    2    1    2      3      1      1      2      2      3     1
5    2    2    2      3      1      1      1      2      3     1
6    2    2    2      4      1      1      2      2      3     1
> env$val.d %>% head()
  ftlm rbci pcci v.fatl v.satl v.rftl v.ftlm v.stlm v.rbci Class
1    2    2    2      2      2      2      1      2      2     1
2    2    2    2      2      2      2      1      2      2     1
3    2    2    2      3      2      2      1      2      2     1
4    2    2    2      4      2      2      2      2      3     1
5    2    2    2      3      2      2      1      2      2     1
6    2    2    2      3      2      2      2      2      2     1

> env$train.d$v.fatl %>% table()
.
  1   2   3   4 
211 693 519 577 
> env$test.d$v.fatl %>% table()
.
  1   2   3   4 
 49 376 313 262 
> env$val.d$v.fatl %>% table()
.
  1   2   3   4 
 68 379 295 258

Дальнейшее использование наборов с дискретными данными зависит от используемой модели. Если это нейросеть, то нужно будет превратить предикторы в dummy-переменные. Насколько хорошо разделены классы этими переменными? Насколько хорошо они коррелируют с целевой? Визуализируем эти зависимости с помощью cross-plot()::funModeling. Cross_plot показывает, как входная переменная коррелирует с целевой переменной, получая коэффициенты правдоподобия для каждого диапазона каждого входа.

Рассмотрим три переменных — v.fatl, ftlm и v.satl, разделенных на 4, 3 и 2 диапазона соответственно. Построим графики:

evalq(
  cross_plot(data = train.d, 
             str_input = c("v.fatl", "ftlm", "v.satl"), 
             str_target = "Class", 
             auto_binning = F,
             plot_type = "both"), #'quantity' 'percentual'
  env
  )

Discret1__1

Рис.16. Cross-plot переменной v.fatl/Class

Discret2__1

Рис.17. Cross-plot переменной ftlm/Class

Discret3__1

Рис.18. Cross-plot переменной v.satl/Class

Как видим, предикторы хорошо коррелированы с уровнями целевой, имеют хорошо выраженные пороги, разделяющие уровни целевой Class.

Можно просто (неоптимально) разделить предикторы на равные участки и посмотреть, как в таком случае они будут коррелировать с целевой. Возьмем три предыдущих переменных и две плохих (stlm, v.rstl) из набора train (недискретизованного). Разделим их на 10 равночастотных участков и посмотрим на cross-plot их с целевой:

evalq(
  cross_plot(
      DT$train  %>% select(-Data) %>%
      select(c(v.satl, ftlm, v.fatl, stlm, v.rstl, Class)) %>%
      as.data.frame(), 
      str_input = Cs(v.satl, ftlm, v.fatl, stlm, v.rstl), 
      str_target = "Class", 
      auto_binning = T,
      plot_type = "both"), #'quantity' 'percentual'
  env
)

Отрисуем пять графиков этих переменных:

Discret4__1

Рис.19. Cross-plot переменной v.satl (10 участков) vs Class

Discret5__1

Рис.20. Cross-plot переменной ftlml (10 участков) vs Class

Discret6__1

Рис.21. Cross-plot переменной v.fatl (10 участков) vs Class

Discret8__1

Рис.22. Cross-plot переменной stlm (10 участков) vs Class

Discret9__1

Рис.23. Cross-plot переменной v.rstl (10 участков) vs Class

Как видно из этих рисунков, и при разделении на 10 дискретных равночастотных участков переменные v.fatl, ftlm и v.satl имеют хорошо выраженный порог разделения уровней целевой. Становится понятно, почему две других переменных (stlm, v.rstl) нерелевантны. Кстати, это эффективный способ определения важности предикторов. Мы еще поговорим об этом ниже.

В заключение посмотрим, как входная переменная коррелирует с целевой переменной, сравнивая байесовским способом posterior conversion rate к целевой переменной. Полезно сравнивать категориальные ценности, которые не имеют внутреннего порядка. Для этого будем использовать функцию bayes_plot::funModeling. Возьмем переменные v.fatl, ftlm и v.satl с наборов train.d, val.d и test.d.

#------BayesTrain-------------------
evalq(
  {
    bayesian_plot(train.d, input = "v.fatl", 
                  target = "Class", 
                  title = "Bayesian comparison train$v.fatl/Class",
                  plot_all = F, extra_above = 5, 
                  extra_under = 5)
  },env
)
evalq(
  {
    bayesian_plot(train.d, input = "ftlm", 
                  target = "Class", 
                  title = "Bayesian comparison train$ftlm/Class",
                  plot_all = F, extra_above = 5, 
                  extra_under = 5)
  },env
)
evalq(
  {
    bayesian_plot(train.d, input = "v.satl", 
                  target = "Class", 
                  title = "Bayesian comparison train$v.satl/Class",
                  plot_all = F, extra_above = 5, 
                  extra_under = 5)
  },env
)
#------------BayesTest------------------------
evalq(
  {
    bayesian_plot(test.d, input = "v.fatl", 
                  target = "Class", 
                  title = "Bayesian comparison test$v.fatl/Class",
                  plot_all = F, extra_above = 5, 
                  extra_under = 5)
  },env
)
evalq(
  {
    bayesian_plot(test.d, input = "ftlm", 
                  target = "Class", 
                  title = "Bayesian comparison test$ftlm/Class",
                  plot_all = F, extra_above = 5, 
                  extra_under = 5)
  },env
)
evalq(
  {
    bayesian_plot(test.d, input = "v.satl", 
                  target = "Class", 
                  title = "Bayesian comparison test$v.satl/Class",
                  plot_all = F, extra_above = 5, 
                  extra_under = 5)
  },env
)
#-------------BayesVal---------------------------------
evalq(
  {
    bayesian_plot(val.d, input = "v.fatl", 
                  target = "Class", 
                  title = "Bayesian comparison val$v.fatl/Class",
                  plot_all = F, extra_above = 5, 
                  extra_under = 5)
  },env
)
evalq(
  {
    bayesian_plot(val.d, input = "ftlm", 
                  target = "Class", 
                  title = "Bayesian comparison val$ftlm/Class",
                  plot_all = F, extra_above = 5, 
                  extra_under = 5)
  },env
)
evalq(
  {
    bayesian_plot(val.d, input = "v.satl", 
                  target = "Class", 
                  title = "Bayesian comparison val$v.satl/Class",
                  plot_all = F, extra_above = 5, 
                  extra_under = 5)
  },env
)
#------------------------------------------

BayesCompTrain__1

Рис.24. Байесовское сравнение переменных с целевой в наборе train

BayesCompValid__1

Рис.25. Байесовское сравнение переменных с целевой в наборе val

BayesCompTest__1

Рис.26. Байесовское сравнение переменных с целевой в наборе test

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

Эту же задачу — разбиение переменной на отрезки, оптимально соответствующие уровням целевой — можно решить и другим путем, с помощью отличного пакета smbinning. Можете проверить это самостоятельно. Напомню и о том, что в предыдущей статье рассмотрен еще один интересный метод дискретизации — с помощью пакета “RoughSets“.

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

 

1.2. Создание новых признаков

Создание переменной — это процесс создания новых переменных на основе уже существующих. Например, у нас есть дата (ДД-мм-гг) в качестве входной переменной в наборе данных. Мы можем создавать новые переменные — день, месяц, год, неделя, день недели, которые будут лучше связаны с целевой переменной. Этот шаг используется, чтобы выделить скрытые взаимосвязи в переменной.

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

Создание фиктивных переменных — еще один популярный способ работы с переменными. Чаще всего фиктивные переменные применяются при преобразовании категориальных переменных в числовые. Категориальная переменная может принимать значения 0 и 1. Мы можем создать фиктивные переменные для более чем двух классов категориальных переменных с N или N-1 переменных.

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

  1. Использовать значения даты и времени как переменные. Можно создавать новые переменные, учитывая различия в датах и времени.
  2. Создать новые соотношения и пропорции: вместо того, чтобы просто держать прошлые входы и выходы в наборе данных, мы можем занести в набор их соотношения, которые могут иметь большое значение.
  3. Применить стандартные преобразования. Глядя на колебания и участки переменной вместе с выходом, мы можем увидеть, улучшится ли корреляция после применения базовых преобразований. Наиболее часто используемые преобразования включают log, exp, квадратичную и тригонометрическую вариации.
  4. Проверить переменные на сезонность и создать модель за нужный период (неделя, месяц, сессия и т.п.).

Интуитивно понятно, что в понедельник поведение рынка не такое, как в среду или четверг. То есть, день недели — важный признак. Не менее важно, в каком периоде суток находится рынок (азиатская, европейская или американская сессии). Как мы можем определить эти признаки?

Используем возможности пакета timekit. Центральная функция пакета — tk_augment_timeseries_signature(), которая добавляет к временным меткам нашего начального набора pr целый ряд временных данных, которые могут быть полезны и как дополнительные признаки, и как параметры группировки. Какие это данные?

Index значение индекса, который был разложен
Index.num числовое значение индекса в секундах. База “1970-01-01 00:00:00”
diff разница в секундах от предыдущего числового значения индекса
Year год, компонента индекса
half половина составляющей индекса
quarter квартал, составляющая индекса
month месяц, составляющая индекса с основанием 1
month.xts месяц, составляющая индекс с базой 0, так как реализовано в xts
month.lbl метка месяца как упорядоченный фактор, начинающийся с января и заканчивающийся декабрем
day день, составляющая индекса
hour час, составляющая индекса
minute минута, составляющая индекса
second секунда, составляющая индекса
hour12 часовая компонента в 12-часовой шкале
am.pm утро (am) = 1, после полудня (pm) = 2
wday день недели с основанием 1. Воскресенье = 1, суббота = 7
wday.xts день недели с основанием 0, как реализовано в xts. Воскресенье = 0, суббота = 6
wday.lbl метка дня недели как упорядоченный фактор, начиная с воскресенья и заканчивая субботой
mday день месяца
qday день квартала
yday день года.
mweek неделя месяца
week номер недели в году (начало с воскресенья)
week.iso номер недели года по ISO (начало с понедельника)
week2 модуль для двухнедельной частоты
week3 модуль для трехнедельной частоты
week4 модуль для quad-недельной частоты

Возьмем начальный набор данных pr, усилим его функцией tk_augment_timeseries_signature(), добавим к начальному набору переменные mday, wday.lbl, hour, удалим неопределенные данные (NA) и сгруппируем данные по дням недели.

evalq(
  {
    tk_augment_timeseries_signature(pr) %>%
    select(c(mday, wday.lbl,  hour)) %>% 
    cbind(pr, .) -> pr.augm
    pr.compl <- pr.augm[complete.cases(pr.augm), ]
    pr.nest <- pr.compl %>% group_by(wday.lbl) %>% nest() 
  },
  env)
> str(env$pr.augm)
'data.frame':   8000 obs. of  33 variables:
 $ Data    : POSIXct, format: "2017-01-10 11:00:00" ...
 $ Open    : num  123 123 123 123 123 ...
 $ High    : num  123 123 123 123 123 ...
 $ Low     : num  123 123 123 123 123 ...
 $ Close   : num  123 123 123 123 123 ...
 $ Vol     : num  3830 3360 3220 3241 3071 ...
 ..................................................
 $ zigz    : num  123 123 123 123 123 ...
 $ dz      : num  NA -0.0162 -0.0162 -0.0162 -0.0162 ...
 $ sig     : num  NA -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
 $ mday    : int  10 10 10 10 10 10 10 10 10 10 ...
 $ wday.lbl: Ord.factor w/ 7 levels "Sunday"<"Monday"<..: 3 3 3 3 3 3 3 3 3 3 ...
 $ hour    : int  11 11 11 11 12 12 12 12 13 13 ...

Тот же результат можем получить с использованием библиотеки lubridate, удалив попутно данные за субботу.

require(lubridate)
evalq({pr %>% mutate(.,
                     wday = wday(Data), #label = TRUE, abbr = TRUE),
                     day = day(Data),
                     hour = hour(Data)) %>%
    filter(wday != "Sat") -> pr1
  pr1.nest <- pr1 %>% na.omit %>% 
    group_by(wday) %>% nest()}, 
  env
)
#-------
str(env$pr1)
'data.frame':   7924 obs. of  33 variables:
 $ Data  : POSIXct, format: "2017-01-10 11:00:00" ...
 $ Open  : num  123 123 123 123 123 ...
 $ High  : num  123 123 123 123 123 ...
 $ Low   : num  123 123 123 123 123 ...
 $ Close : num  123 123 123 123 123 ...
 $ Vol   : num  3830 3360 3220 3241 3071 ...
 ..........................................
 $ zigz  : num  123 123 123 123 123 ...
 $ dz    : num  NA -0.0162 -0.0162 -0.0162 -0.0162 ...
 $ sig   : num  NA -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
 $ wday  : int  3 3 3 3 3 3 3 3 3 3 ...
 $ day   : int  10 10 10 10 10 10 10 10 10 10 ...
 $ hour  : int  11 11 11 11 12 12 12 12 13 13 ...

Сгруппированные по дням недели данные выглядят так (воскресенье = 1, понедельник = 2 и т.д.):

> env$pr1.nest
# A tibble: 5 × 2
   wday                  data
  <int>                <list>
1     4 <tibble [1,593 Ч 32]>
2     5 <tibble [1,632 Ч 32]>
3     6 <tibble [1,624 Ч 32]>
4     2 <tibble [1,448 Ч 32]>
5     3 <tibble [1,536 Ч 32]>

Из данных в основном наборе pr можно дополнительно использовать переменные dL, dH на последних 3 барах.

2. Выбор предикторов

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

2.1. Визуальная оценка

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

Соберем в один набор log-трансформированные, sin-трансформированные, добавим tanh-преобразованные и нормализованные данные и оценим зависимость целевой и предикторов в этих наборах.

Последовательность обработки первичного набора (приведена на ниже показанном рисунке) следующая: сырые данные dataSet очищаем (без удаления высококоррелированных) и, разделив на train/val/test, получаем набор DT. Дальше по каждому виду трансформации выполним действия по схеме. Соберем все в один скрипт:

Preprocess__1

Рис.27. Структурная схема препроцессинга

Очистим набор, разделим на train/val/test и уберем ненужные данные:

#----Clean---------------------
require(caret)
require(pipeR)
evalq(
  {
    train = 1:2000
    val = 2001:3000
    test = 3001:4000
    DT <- list()
    dataSet %>%
      preProcess(., method = c("zv", "nzv", "conditionalX")) %>%
      predict(., dataSet) %>%
      na.omit -> dataSetClean
    list(train = dataSetClean[train, ], 
         val = dataSetClean[val, ], 
         test = dataSetClean[test, ]) -> DT
    rm(dataSetClean, train, val, test)
  }, 
  env)

Обработаем выбросы во всех наборах:

#------outlier-------------
evalq({
# определим новый лист для результата
  DTcap <- list()
# пройдемся по трем наборам
  foreach(i = 1:3) %do% {
    DT %>% 
# уберем колонки c(Data, Class)
      select(-c(Data, Class)) %>%
# преобразуем в  data.frame и сохраним во временную переменную x
      as.data.frame() -> x
    if (i == 1) {
# при первом входе определеим параметры выбросов
      foreach(i = 1:ncol(x), .combine = "cbind") %do% {
        prep.outlier(x[ ,i]) %>% unlist()
      } -> pre.outl
      colnames(pre.outl) <- colnames(x)
    } 
# заменим выбросы на 5/95 % и сохраним результат в x.cap
    foreach(i = 1:ncol(x), .combine = "cbind") %do% {
      stopifnot(exists("pre.outl", envir = env))
      lower = pre.outl['lower.25%', i] 
      upper = pre.outl['upper.75%', i]
      med = pre.outl['med', i]
      cap1 = pre.outl['cap1.5%', i] 
      cap2 = pre.outl['cap2.95%', i] 
      treatOutlier(x = x[ ,i], impute = T, fill = T, 
                   lower = lower, upper = upper, 
                   med = med, cap1 = cap1, cap2 = cap2) 
      } %>% as.data.frame() -> x.cap
    colnames(x.cap) <- colnames(x)
    return(x.cap)
  } -> Dtcap
#уберем ненужные переменные
  rm(lower, upper, med, cap1, cap2, x.cap, x)
}, env)

Трансформируем переменные во всех наборах Dtcap без выбросов функцией log(x+1). Получим DTLn лист с тремя наборами логарифмированных переменных.

#------logtrans------------
evalq({
  DTLn <- list()
  foreach(i = 1:3) %do% {
    DTcap %>% 
      apply(., 2, function(x) log2(x + 1)) %>%
      as.data.frame() %>%
      cbind(., Class = DT$Class)
  } -> DTLn
},
env)

Трансформируем переменные во всех наборах Dtcap без выбросов функцией sin(2*pi*x). Получим DTSin лист с тремя наборами sin-трансформированных переменных.

#------sintrans--------------
evalq({
  DTSin <- list()
  foreach(i = 1:3) %do% {
    DTcap %>% 
      apply(., 2, function(x) sin(2*pi*x)) %>%
      as.data.frame() %>%
      cbind(., Class = DT$Class)
  } -> DTSin
},
env)

Трансформируем переменные во всех наборах Dtcap без выбросов функцией tanh(x). Получим DTTanh лист с тремя наборами tanh-трансформированных переменных.

#------tanhTrans----------
evalq({
  DTTanh <- list()
  foreach(i = 1:3) %do% {
    DTcap %>% 
      apply(., 2, function(x) tanh(x)) %>%
      as.data.frame() %>%
      cbind(., Class = DT$Class)
  } -> DTTanh
},
env)

Нормализуем наборы DT, DTLn, DTSin, DTTanh.

#------normalize-----------
evalq(
  {
# определяем параметры нормализации
    preProcess(DT$train, method = "spatialSign") -> preproc 
    list(train = predict(preproc, DT$train), 
         val = predict(preproc, DT$val),
         test = predict(preproc, DT$test)
 ) -> DTn
  }, 
  env) 
#--ln---
evalq(
  {
    preProcess(DTLn, method = "spatialSign") -> preprocLn 
    list(train = predict(preprocLn, DTLn), 
         val = predict(preprocLn, DTLn),
         test = predict(preprocLn, DTLn)
    ) -> DTLn.n
  }, 
  env)
#---sin---
evalq(
  {
    preProcess(DTSin, method = "spatialSign") ->  preprocSin 
    list(train = predict(preprocSin, DTSin), 
         val = predict(preprocSin, DTSin),
         test = predict(preprocSin, DTSin)
    ) -> DTSin.n
  }, 
  env)
#-----tanh-----------------
evalq(
  {
    preProcess(DTTanh, method = "spatialSign") -> preprocTanh 
    list(train = predict(preprocTanh, DTTanh), 
         val = predict(preprocTanh, DTTanh),
         test = predict(preprocTanh, DTTanh)
    ) -> DTTanh.n
  }, 
  env)

Дискретизируем наш набор DT с помощью функции mdlp::discretization.

##------discretize----------
#--------preCut---------------------
# определяем параметры дискретизации (cutpoint)
require(pipeR)
require(discretization)
evalq(
  #require(pipeR) 
# занимает некоторое время !
  pipeline({
    DT$train
    select(-Data)
    as.data.frame()
    mdlp() 
  }) -> mdlp.train, 
  env)
#-------cut_opt----------
evalq(
  {
    DTd <- list()
    mdlp.train$cutp %>% 
# определяем колонки которые нужно дискретизировать
      lapply(., function(x) is.numeric(x)) %>%
      unlist -> idx   # bool
    #----train-----------------
    mdlp.train$Disc.data[ ,idx] -> DTd$train 
    #---test------------
    DT$test %>% 
      select(-c(Data, Class)) %>%
      as.data.frame() -> test.d
# переразметим данные в соответствии с расчитанными диапазонами    
    foreach(i = 1:length(idx), .combine = 'cbind') %do% {
      if (idx[i]) {
        findInterval(test.d[ ,i], 
        vec = mdlp.train$cutp,
        rightmost.closed = FALSE, 
        all.inside = F,
        left.open = F)
        }
    } %>% as.data.frame() %>% add(1) %>%
      cbind(., DT$test$Class) -> DTd$test
    colnames(DTd$test) <- colnames(DTd$train)
    #-----val-----------------
    DT$val %>% 
      select(-c(Data, Class)) %>%
      as.data.frame() -> val.d
# переразметим данные в соответствии с расчитанными диапазонами  
    foreach(i = 1:length(idx), .combine = 'cbind') %do% {
      if (idx[i]) {
        findInterval(val.d[ ,i], 
        vec = mdlp.train$cutp,
        rightmost.closed = FALSE, 
        all.inside = F,
        left.open = F)
        }
    } %>% as.data.frame() %>% add(1) %>%
      cbind(., DT$val$Class) -> DTd$val 
    colnames(DTd$val) <- colnames(DTd$train)
# подчистим за собой
    rm(test.d, val.d)
  }, 
  env
)

Напомним, какие переменные у нас в начальном наборе DT$train:

require(funModeling)
plot_num(env$DT$train %>% select(-Data), bins = 20)

FSelect_1__1

Рис.28. Распределение переменных набора DT$train

Используем возможности пакета smbinning для определения релевантных предикторов в поднаборах train всех нормализованных наборов, полученных нами ранее (Dtn, DTLn.n, DTSin.n и DTTanh.n). Поскольку целевая в этом пакете должна быть числовой со значениями (0, 1) напишем функцию для необходимых преобразований.

#--------------------------------
require(smbinning)
targ.int <- function(x){
  x %>% tbl_df() %>%
  mutate(Cl = (as.numeric(Class) - 1) %>%
           as.integer()) %>%
  select(-Class) %>% as.data.frame()
}

Кроме того, пакет не любит переменные с точкой в имени. Напишем функцию, которая переименует все переменные с точкой в переменные с нижней чертой.

renamepr <- function(X){
  X %<>% rename(v_fatl = v.fatl,
               v_satl = v.satl,
               v_rftl = v.rftl,
               v_rstl = v.rstl,
               v_ftlm = v.ftlm,
               v_stlm = v.stlm,
               v_rbci = v.rbci,
               v_pcci = v.pcci)
  return(X)
}

Вычислим и отрисуем графики релевантных предикторов.

par(mfrow = c(2,2))
#--Ln--------------
evalq({
  df <- renamepr(DTLn.n) %>% targ.int
  sumivt.ln.n = smbinning.sumiv(df = df, y = 'Cl')
  smbinning.sumiv.plot(sumivt.ln.n, cex = 0.7)
  rm(df)
}, 
env)
#---Sin-----------------
evalq({
  df <- renamepr(DTSin.n) %>% targ.int
  sumivt.sin.n = smbinning.sumiv(df = df, y = 'Cl')
  smbinning.sumiv.plot(sumivt.sin.n, cex = 0.7)
  rm(df)
  }, 
env)
#---norm-------------
evalq({
  df <- renamepr(DTn) %>% targ.int
  sumivt.n = smbinning.sumiv(df = df, y = 'Cl')
  smbinning.sumiv.plot(sumivt.n, cex = 0.7)
  rm(df)
  }, 
env)
#-----Tanh----------------
evalq({
  df <- renamepr(DTTanh.n) %>% targ.int
  sumivt.tanh.n = smbinning.sumiv(df = df, y = 'Cl')
  smbinning.sumiv.plot(sumivt.tanh.n, cex = 0.7)
  rm(df)
  }, 
env)
par(mfrow = c(1,1))

FSelect_2__1

Рис.29. Важность предикторов в поднаборе train нормализованных наборов

Во всех наборах сильны пять предикторов — v_fatl, ftlm, v_satl, rbci, v_rbci, правда, в различном порядке. Средняя сила у четырех предикторов — pcci, v_ftlm, v_stlm, v_rftl. Слабые — v_pcci и stlm. Можно посмотреть по каждому набору распределение переменных в порядке их значимости:

env$sumivt.ln.n
     Char     IV               Process
5  v_fatl 0.6823    Numeric binning OK
1    ftlm 0.4926    Numeric binning OK
6  v_satl 0.3737    Numeric binning OK
3    rbci 0.3551    Numeric binning OK
11 v_rbci 0.3424    Numeric binning OK
10 v_stlm 0.2591    Numeric binning OK
4    pcci 0.2440    Numeric binning OK
9  v_ftlm 0.2023    Numeric binning OK
7  v_rftl 0.1442    Numeric binning OK
12 v_pcci 0.0222    Numeric binning OK
2    stlm     NA No significant splits
8  v_rstl     NA No significant splits

Последние три переменных можно спокойно отбросить. Т.е., останутся пять сильных и четыре средних. Определим имена лучших переменных (IV > 0.1).

evalq(sumivt.sin.n$Char[sumivt.sin.n$IV > 0.1] %>% 
        na.omit %>% as.character() -> best.sin.n, 
      env)
> env$best.sin.n
[1] "v_fatl" "ftlm"   "rbci"   "v_rbci" "v_satl" "pcci"  
[7] "v_ftlm" "v_stlm" "v_rftl"

Посмотрим более подробно переменные v_fatl и ftlm.

evalq({
    df <- renamepr(DTTanh.n) %>% targ.int
    x = 'v_fatl'
    y = 'Cl'
    res <- smbinning(df = df, 
                        y = y,
                        x = x) 
  #res$ivtable # Tabulation and Information Value
  #res$iv # Information value
  #res$bands # Bins or bands
  #res$ctree  # Decision tree from partykit
  par(mfrow = c(2,2))
  sub = paste0(x, "  vs  ", y) #rbci vs Cl"
  boxplot(df~df,
          horizontal = TRUE, 
          frame = FALSE, col = "lightblue",
          main = "Distribution")
  mtext(sub,3) #ftlm
  smbinning.plot(res, option = "dist",
                 sub = sub) #"pcci vs Cl")
  smbinning.plot(res, option = "goodrate", #"badrate"
                 sub = sub) #"pcci vs Cl")
  smbinning.plot(res, option = "WoE",
                 sub = sub) #"pcci vs Cl")
  par(mfrow = c(1, 1))
}, env)

FSelect_3__1

Рис.30. Связь диапазонов переменной v_fatl с целевой Cl

В объекте res, кроме прочей полезной информации, содержатся точки раздела переменной на диапазоны, оптимально связанные с целевой. В нашем конкретном случае диапазонов 4.

> env$res$cuts
[1] -0.3722 -0.0433  0.1482

Проведем аналогичный расчет для переменной ftlm и построим графики:

FSelect_4__1

Рис.31. Связь диапазонов переменной ftlm c целевой Cl

Точки разделения диапазонов:

> env$res$cuts
[1] -0.2084 -0.0150  0.2216

Имея точки разделения, мы можем дискретизировать переменные в наших наборах и посмотреть, насколько отличаются:

  • важные переменные, ранее определенные с помощью функции mdlp::discretization от переменных, определенных с помощью функции smbinning::smbinning;
  • разбиение переменных на диапазоны.

У нас уже есть один набор с данными, дискретизированными функцией mdlp::discretization DTd. Проделаем то же самое, но с помощью функции smbinning::smbinning и только train поднабором.

Определим параметры дискретизации:

evalq({
  res <- list()
  DT$train %>% renamepr() %>% targ.int() -> df
  x <- colnames(df)
  y <- "Cl"
  foreach(i = 1:(ncol(df) - 1)) %do% {
    smbinning(df, y = y, x = x[i])
  } -> res
  res %>% lapply(., function(x) x[1] %>% is.list) %>%
    unlist -> idx
}, env)

Дискретизируем поднабор DT$train:

evalq({
  DT1.d <- list()
  DT$train %>% renamepr() %>% 
    targ.int() %>% select(-Cl) -> train
  foreach(i = 1:length(idx), .combine = 'cbind') %do% {
    if (idx[i]) {
      findInterval(train[ ,i], 
                   vec = res$cuts,
                   rightmost.closed = FALSE, 
                   all.inside = F,
                   left.open = F)
    }
  } %>% as.data.frame() %>% add(1) %>%
    cbind(., DT$train$Class) -> DT1.d$train
  colnames(DT1.d$train) <- colnames(train)[idx] %>%
    c(., 'Class')
},
env)

Определим лучшие переменные с информационной значимостью больше 0.1, упорядоченные по возрастанию:

evalq({
  DT$train %>% renamepr() %>% targ.int() -> df
  sumivt.dt1.d = smbinning.sumiv(df = df, y = 'Cl')
  sumivt.dt1.d$Char[sumivt.dt1.d$IV > 0.1] %>% 
    na.omit %>% as.character() -> best.dt1.d
  rm(df)
}, 
env)

Отрисуем график разбиения переменных в наборе DTd$train:

require(funModeling)
plot_num(env$DTd$train)

FSelect_5__1

Рис.32. Переменные набора DT$ train, дискретизированные функцией mdlp

График набора DT1.d со всеми переменными и с лучшими приведен ниже.

plot_num(env$DT1.d$train)

FSelect_6__1

Рис.33. Переменные набора DT1 d$train, дискретизированные функцией smbinning (все)

plot_num(env$DT1.d$train[ ,env$best.dt1.d])

FSelect_7__1

Рис.34. Переменные набора DT1.d$train, дискретизированные функцией smbinning (лучшие упорядоченные по возрастанию информационной значимости)

Что мы видим из графиков? Переменные, определенные как важные — одинаковы в обеих случаях. А вот разбиение на диапазоны отличается. Нужно проверить, какой вариант дает лучший предсказательный результат с моделью.

2.2. Аналитическая оценка

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

Используем пакет varbvs. В одноименной функции реализованы быстрые алгоритмы для установки байесовских моделей выбора переменных и вычисления коэффициентов Байеса, в которых результат (или целевая) моделируется с использованием линейной регрессии или логистической регрессии. Алгоритмы основаны на вариационных приближениях, описанных в “Масштабируемом вариационном выводе для выбора байесовских переменных в регрессии и его точности в исследованиях генетических ассоциаций” (P. Carbonetto and M. Stephens, Bayesian Analysis 7, 2012, pages 73-108). Это программное обеспечение было применено к большим наборам данных с более чем миллионом переменных и с тысячами образцов.

Функция varbvs() принимает на вход матрицу и целевую — численный вектор (0, 1). Возьмем наш набор с нормализованными данными DTTanh.n$train и проверим, какие предикторы будут определяться как важные по этому методу.

require(varbvs)
evalq({
  train <- DTTanh.n$train %>% targ.int() %>%  as.matrix()
  fit <- varbvs(X = train[ ,-ncol(train)] , 
                Z = NULL,
                y = train[ ,ncol(train)] %>% as.vector(),
                "binomial", 
                logodds = seq(-2,-0.5,0.1),
                optimize.eta = T,
                initialize.params = T,
                verbose = T, nr = 100
                )
  print(summary(fit))
}, env)

Welcome to           --       *                              *               
VARBVS version 2.0.3 --       |              |               |               
large-scale Bayesian --       ||           | |    |          || |     |   |  
variable selection   -- |     || | |    |  | ||  ||        |||| ||    |   || 
****************************************************************************
Copyright (C) 2012-2017 Peter Carbonetto.
See http://www.gnu.org/licenses/gpl.html for the full license.
Fitting variational approximation for Bayesian variable selection model.
family:     binomial   num. hyperparameter settings: 16 
samples:    2000       convergence tolerance         1.0e-04
variables:  12         iid variable selection prior: yes 
covariates: 0          fit prior var. of coefs (sa): yes 
intercept:  yes        fit approx. factors (eta):    yes 
Finding best initialization for 16 combinations of hyperparameters.
-iteration-   variational    max.   incl variance params
outer inner   lower bound  change   vars   sigma      sa
 0016 00018 -1.204193e+03 6.1e-05 0003.3      NA 3.3e+00
Computing marginal likelihood for 16 combinations of hyperparameters.
-iteration-   variational    max.   incl variance params
outer inner   lower bound  change   vars   sigma      sa
 0016 00002 -1.204193e+03 3.2e-05 0003.3      NA 3.3e+00
Summary of fitted Bayesian variable selection model:
family:     binomial   num. hyperparameter settings: 16
samples:    2000       iid variable selection prior: yes
variables:  12         fit prior var. of coefs (sa): yes
covariates: 1          fit approx. factors (eta):    yes
maximum log-likelihood lower bound: -1204.1931
Hyperparameters: 
        estimate Pr>0.95             candidate values
sa          3.49 [3.25,3.6]          NA--NA
logodds    -0.75 [-1.30,-0.50]       (-2.00)--(-0.50)
Selected variables by probability cutoff:
>0.10 >0.25 >0.50 >0.75 >0.90 >0.95 
    3     3     3     3     3     3 
Top 5 variables by inclusion probability:
  index variable   prob PVE coef*  Pr(coef.>0.95)
1     1     ftlm 1.0000  NA 2.442 [+2.104,+2.900]
2     4     pcci 1.0000  NA 2.088 [+1.763,+2.391]
3     3     rbci 0.9558  NA 0.709 [+0.369,+1.051]
4    10   v.stlm 0.0356  NA 0.197 [-0.137,+0.529]
5     6   v.satl 0.0325  NA 0.185 [-0.136,+0.501]
*See help(varbvs) about interpreting coefficients in logistic regression.

Как видим, определены пять лучших предикторов (ftlm, pcci, rbci, v.stlm, v.satl). Они входят в девятку лучших, которую мы определили раньше, но в другом порядке и с другими весами значимости. Поскольку модель у нас уже есть, давайте проверим, какое качество результата мы получим на валидационном и тестовом наборе.

Валидационный набор:

#-----------------
evalq({
  val <- DTTanh.n$val %>% targ.int() %>%
    as.matrix()
  y = val[ ,ncol(val)] %>% as.vector()
  pr <- predict(fit, X = val[ ,-ncol(val)] , 
                Z = NULL)
  
}, env)
cm.val <- confusionMatrix(table(env$y, env$pr))
> cm.val
Confusion Matrix and Statistics

   
      0   1
  0 347 204
  1 137 312
                                          
               Accuracy : 0.659           
                 95% CI : (0.6287, 0.6884)
    No Information Rate : 0.516           
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.3202          
 Mcnemar's Test P-Value : 0.0003514       
                                          
            Sensitivity : 0.7169          
            Specificity : 0.6047          
         Pos Pred Value : 0.6298          
         Neg Pred Value : 0.6949          
             Prevalence : 0.4840          
         Detection Rate : 0.3470          
   Detection Prevalence : 0.5510          
      Balanced Accuracy : 0.6608          
                                          
       'Positive' Class : 0

Результат, прямо скажем, не ахти. Тестовый набор:

evalq({
  test <- DTTanh.n$test %>% targ.int() %>%
    as.matrix()
  y = test[ ,ncol(test)] %>% as.vector()
  pr <- predict(fit, X = test[ ,-ncol(test)] , 
                Z = NULL)
  
}, env)
cm.test <- confusionMatrix(table(env$y, env$pr))
> cm.test
Confusion Matrix and Statistics

   
      0   1
  0 270 140
  1 186 404
                                        
               Accuracy : 0.674         
                 95% CI : (0.644, 0.703)
    No Information Rate : 0.544         
    P-Value [Acc > NIR] : < 2e-16       
                                        
                  Kappa : 0.3375        
 Mcnemar's Test P-Value : 0.01269       
                                        
            Sensitivity : 0.5921        
            Specificity : 0.7426        
         Pos Pred Value : 0.6585        
         Neg Pred Value : 0.6847        
             Prevalence : 0.4560        
         Detection Rate : 0.2700        
   Detection Prevalence : 0.4100        
      Balanced Accuracy : 0.6674        
                                        
       'Positive' Class : 0

Результат практически такой же. Это значит, что модель не переучена и хорошо обобщает данные.

Итак: лучшие по версии varbvs — ftlm, pcci, rbci, v.stlm, v.satl.

2.3. Нейросеть

Поскольку мы исследуем область нейросетей, проверим, какие предикторы выберет в качестве наиболее важных сама нейросеть.

Используем пакет FCNN4R, который предоставляет интерфейс для программ ядра из библиотеки FCNN C ++. FCNN основан на совершенно новом представлении искусственной нейронной сети, которое предполагает эффективность, модульность и расширяемость. FCNN4R обеспечивает стандартное обучение (backpropagation, Rprop, имитация отжига, стохастический градиент) и алгоритмы обрезки (minimum magnitude, Optimal Brain Surgeon), но прежде всего я вижу в ней эффективный вычислительный движок.

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

Сети могут быть экспортированы в функции C, чтобы интегрировать их в практически любое программное решение.

Создаем полносвязную нейросеть с двумя скрытыми слоями. Количество нейронов в каждом слое: входной = 12 (количество предикторов), выходной = 1. Инициируем нейроны случайными весами в диапазоне +/- 0.17. Установим функции активации в каждом слое нейросети (кроме входного) = c(“tanh”, “tanh”, “sigmoid”). Подготовим наборы train/val/test.

Ниже представлен скрипт, выполняющий эту последовательность действий.

require(FCNN4R)
evalq({
mlp_net(layers = c(12, 8, 5, 1), name = "n.tanh") %>%
  mlp_rnd_weights(a = 0.17) %>% 
  mlp_set_activation(layer = c(2, 3, 4), 
  activation = c("tanh", "tanh", "sigmoid"), #"threshold", "sym_threshold",
                                            #"linear", "sigmoid", "sym_sigmoid",
                                            #"tanh", "sigmoid_approx",
                                            #"sym_sigmoid_approx"), 
                 slope = 0) -> Ntanh #show() 
#-------
train <- DTTanh.n$train %>% targ.int() %>% as.matrix()
test <- DTTanh.n$test %>% targ.int() %>%  as.matrix()
val <- DTTanh.n$val %>% targ.int() %>% as.matrix()
}, env)

Из предлагаемых методов обучения будем использовать rprop. Зададим константы: tol — уровень ошибки, при достижении которого нужно остановить обучение, max_ep — количество эпох, после выполнения которых нужно остановить обучение, l2reg — коэффициент регуляризации. Обучим модель с этими параметрами, посмотрим визуально, какую сеть и ошибку обучения мы получили.

evalq({
  tol <- 1e-1
  max_ep = 1000
  l2reg = 0.0001
net_rp <- mlp_teach_rprop(Ntanh, 
                          input = train[ ,-ncol(train)], 
                          output = train[ ,ncol(train)] %>% as.matrix(),
                          tol_level = tol, 
                          max_epochs = max_ep, 
                          l2reg = l2reg,
                          u = 1.2, d = 0.5, 
                          gmax = 50, gmin = 1e-06, 
                          report_freq = 100)
}, env)
plot(env$net_rp$mse, t = "l", 
     main = paste0("max_epochs =", env$max_ep, " l2reg = ", env$l2reg))

NN1__1

Рис.35. Ошибка обучения нейросети

evalq(mlp_plot(net_rp$net, FALSE), envir = env)

NN2__1

Рис.36. Структура нейросети

Обрезка (Prune)

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

evalq({
  tol <- 1e-1
  max_ep = 1000
  l2reg = 0.0001
  mlp_prune_mag(net_rp$net, 
                input = train[ ,-ncol(train)], 
                output = train[ ,ncol(train)] %>% as.matrix(),
                tol_level = tol,  
                max_reteach_epochs = max_ep, 
                report = FALSE,
                plots = TRUE) -> net_rp_prune
  
}, env)

NN3__1

Рис.37. Обрезанная нейросеть

Мы видим, что при конкретной структуре нейросети, начальной инициализации, функциях активации и ошибке обучения мы вполне можем обойтись сетью структуры (12, 2, 1, 1). Какие предикторы при этом выбрала нейросеть?

evalq(
  best <- train %>% tbl_df %>%  select(c(1,5,7,8,10,12)) %>% colnames(),
           env)
env$best
[1] "ftlm"   "v.fatl" "v.rftl" "v.rstl" "v.stlm"
[6] "v.pcci"

Двух переменных — v.rstl и v.pcci — нет в составе лучших 9, определенных ранее.

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

Экспериментируйте!

Comments ()