R: обработка пропущенных значений

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

Набор данных и подготовка

Будем использовать набор данных BostonHousing из пакета mlbench для иллюстрации разных подходов к обработке пропущенных значений. Хотя в исходных данных BostonHousing нет пропущенных значений, я внесу их случайным образом. Благодаря этому мы сможем сравнивать вычисленные пропущенные значения с фактическими, чтобы оценить эффективность подходов к восстановлению данных. Давайте начнем с импорта данных из пакета mlbench и случайным образом внесем пропущенные значения (NA).

# инициализация данных
data ("BostonHousing", package="mlbench")
original <- BostonHousing  # сохранение исходных данных
# внесение пропущенных значений
set.seed(100)
BostonHousing[sample(1:nrow(BostonHousing), 40), "rad"] <- NA
BostonHousing[sample(1:nrow(BostonHousing), 40), "ptratio"]

#>      crim zn indus chas   nox    rm  age    dis rad tax ptratio      b lstat medv
#> 1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98 24.0
#> 2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14 21.6
#> 3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03 34.7
#> 4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94 33.4
#> 5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33 36.2
#> 6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21 28.7

Были внесены пропущенные значения. И хотя мы знаем, где они, давайте сделаем небольшую проверку с помощью mice::md.pattern.

# набор пропущенных значений
library(mice)
md.pattern(BostonHousing)  # набор пропущенных значений в данных

#>     crim zn indus chas nox rm age dis tax b lstat medv rad ptratio   
#> 431    1  1     1    1   1  1   1   1   1 1     1    1   1       1  0
#>  35    1  1     1    1   1  1   1   1   1 1     1    1   0       1  1
#>  35    1  1     1    1   1  1   1   1   1 1     1    1   1       0  1
#>   5    1  1     1    1   1  1   1   1   1 1     1    1   0       0  2
#>        0  0     0    0   0  0   0   0   0 0     0    0  40      40 80

В принципе, есть четыре способа обработки пропущенных значений.

1. Удаление данных

Если в вашем наборе сравнительно большое количество данных, где все требуемые классы достаточно представлены в данных режима обучения, то попробуйте удалить данные (строки), содержащие пропущенные значения (или не учитывать пропущенные значения при создании модели, например, установив na.action=na.omit). Убедитесь, что после удаления данных у вас:

1. достаточно точек, чтобы модель не потеряла достоверность;
2. не появилась погрешность (т.е. непропорциональность или отсутствие какого-либо класса).

# Пример
lm(medv ~ ptratio + rad, data=BostonHousing, na.action=na.omit)

2. Удаление переменной

Если в какой-то конкретной переменной больше пропущенных значений, чем в остальных вместе взятых, и если удалив ее, тем самым можно сохранить много данных, я бы предложил удалить эту переменную. Конечно, если она не является действительно значимым фактором. Фактически, это принятие решения — потерять переменную или часть данных.

3. Оценка средним, медианой, модой

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

library(Hmisc)
impute(BostonHousing$ptratio, mean)  # заменить средним
impute(BostonHousing$ptratio, median)  # медианой
impute(BostonHousing$ptratio, 20)  # заменить заданным числом
# или если хотите заменить вручную
BostonHousing$ptratio[is.na(BostonHousing$ptratio)] <- mean(BostonHousing$ptratio, na.rm = T)

Давайте посчитаем точность в случае замены средним:

library(DMwR)
actuals <- original$ptratio[is.na(BostonHousing$ptratio)]
predicteds <- rep(mean(BostonHousing$ptratio, na.rm=T), length(actuals))
regr.eval(actuals, predicteds)

#>        mae        mse       rmse       mape 
#> 1.62324034 4.19306071 2.04769644 0.09545664

4. Прогнозирование

Прогнозирование — самый сложный метод замены пропущенных значений. Он включает следующие подходы: kNN-оценка, rpart и mice.

4.1. kNN-оценка

DMwR::knnImputation использует метод k ближайших соседей для замены пропущенных значений. Проще говоря, kNN-оценка делает следующее. Для каждого данного, требующего замены, определяется k ближайших точек на основании евклидового расстояния, и рассчитывается их взвешенное (по расстоянию) среднее.

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

library(DMwR)
knnOutput <- knnImputation(BostonHousing[, !names(BostonHousing) %in% "medv"])  # выполнить knn-оценку
anyNA(knnOutput)

#> FALSE

Давайте оценим точность:

actuals <- original$ptratio[is.na(BostonHousing$ptratio)]
predicteds <- knnOutput[is.na(BostonHousing$ptratio), "ptratio"]
regr.eval(actuals, predicteds)

#>        mae        mse       rmse       mape 
#> 1.00188715 1.97910183 1.40680554 0.05859526

Средняя абсолютная ошибка в процентах (mape) улучшилась примерно на 39% по сравнению с заменой средним. Неплохо.

4.2 rpart

Ограничение DMwR::knnImputation состоит в том, что иногда эту функцию не получится использовать, если пропущены значения факторной переменной. И rpart, и mice подходят для такого случая. Преимущество rpart в том, что достаточно хотя бы одной переменной, не содержащей NA.

Теперь воспользуемся rpart для замены пропущенных значений вместо kNN. Для того, чтобы обработать факторную переменную, нужно установить method=class при вызове rpart(). Для числовых значений будем использовать method=anova. В этом случае также надо убедиться, что в обучении rpart не используется выходная переменная (medv).

library(rpart)
class_mod <- rpart(rad ~ . - medv, data=BostonHousing[!is.na(BostonHousing$rad), ], method="class", na.action=na.omit)  # т.к. rad - факторная переменная
anova_mod <- rpart(ptratio ~ . - medv, data=BostonHousing[!is.na(BostonHousing$ptratio), ], method="anova", na.action=na.omit)  # т.к. ptratio - числовая переменная
rad_pred <- predict(class_mod, BostonHousing[is.na(BostonHousing$rad), ])
ptratio_pred <- predict(anova_mod, BostonHousing[is.na(BostonHousing$ptratio), ])

Посчитаем точность для ptratio:

actuals <- original$ptratio[is.na(BostonHousing$ptratio)]
predicteds <- ptratio_pred
regr.eval(actuals, predicteds)

#>        mae        mse       rmse       mape 
#> 0.71061673 0.99693845 0.99846805 0.04099908 

Средняя абсолютная ошибка в процентах (mape) улучшилась еще примерно на 30% по сравнению с kNN-оценкой. Очень хорошо.

Точность для rad:

actuals <- original$rad[is.na(BostonHousing$rad)]
predicteds <- as.numeric(colnames(rad_pred)[apply(rad_pred, 1, which.max)])
mean(actuals != predicteds)  # расчет ошибки неправильной классификации

#> 0.25

Ошибка неправильной классификации — 25%. Неплохо для факторной переменной!

4.3 mice

mice — сокращение от Multivariate Imputation by Chained Equations (многомерная оценка цепными уравнениями) — пакет R, предоставляющий сложные функции для работы с пропущенными значениями. Он использует немного необычный способ оценки в два шага: mice() для построения модели и complete() для генерации данных. Функция mice(df) создает несколько полных копий df, каждая со своей оценкой пропущенных данных. Функция complete() возвращает один или несколько наборов данных, набор по умолчанию будет первым. Давайте посмотрим, как заменить rad и ptratio:

library(mice)
miceMod <- mice(BostonHousing[, !names(BostonHousing) %in% "medv"], method="rf")  # оценка mice на основе случайных лесов
miceOutput <- complete(miceMod)  # сгенерировать полные данные
anyNA(miceOutput)

#> FALSE

Рассчитаем точность ptratio:

actuals <- original$ptratio[is.na(BostonHousing$ptratio)]
predicteds <- miceOutput[is.na(BostonHousing$ptratio), "ptratio"]
regr.eval(actuals, predicteds)

#>        mae        mse       rmse       mape 
#> 0.36500000 0.78100000 0.88374204 0.02121326

Средняя абсолютная ошибка в процентах (mape) улучшилась еще примерно на 48% по сравнению с rpart. Отлично!

Рассчитаем точность для rad:

actuals <- original$rad[is.na(BostonHousing$rad)]
predicteds <- miceOutput[is.na(BostonHousing$rad), "rad"]
mean(actuals != predicteds)  # расчет ошибки неправильной классификации

#> 0.15

Ошибка неправильной классификации сократилась до 15%, т.е. 6 из 40 наблюдений. Это значительное улучшение по сравнению с 25% для rpart.

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