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

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

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

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