Let’s fix NAs!

download.file("https://bitbucket.org/kailexx/fixnas/raw/ae65f7939974e709f10aa50c96c368120487a7f2/train.csv", destfile="train.csv", method= "wget") train <- read.csv("train.csv", na.strings = c(NA, "")) 

Посмотрим, что из себя представляет содержимое файла:

str(train) 'data.frame':	891 obs. of  12 variables:  $ PassengerId: int  1 2 3 4 5 6 7 8 9 10 ...  $ Survived   : int  0 1 1 1 0 0 0 0 1 1 ...  $ Pclass     : int  3 1 3 1 3 3 1 3 3 2 ...  $ Name       : Factor w/ 891 levels "Abbing, Mr. Anthony",..: 109 191 358 277 16 559 520 629 416 581 ...  $ Sex        : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...  $ Age        : num  22 38 26 35 35 NA 54 2 27 14 ...  $ SibSp      : int  1 1 0 1 0 0 0 3 0 1 ...  $ Parch      : int  0 0 0 0 0 0 0 1 2 0 ...  $ Ticket     : Factor w/ 681 levels "110152","110413",..: 525 596 662 50 473 276 86 396 345 133 ...  $ Fare       : num  7.25 71.28 7.92 53.1 8.05 ...  $ Cabin      : Factor w/ 147 levels "A10","A14","A16",..: NA 82 NA 56 NA NA 130 NA NA NA ...  $ Embarked   : Factor w/ 3 levels "C","Q","S": 3 1 3 3 3 2 3 3 3 1 ...  sum(is.na(train)) [1] 866 

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

Два простых решения

Первое решение — для ленивых: делать практически ничего не надо, т.к. многие функции в R имеют аргумент na.rm, выставив значение которого в TRUE, мы заставим R просто удалить NA перед тем, как что-то делать с данным. Другие же функции имеют параметр na.action, который может принимать такие значения:
na.fail — возвращать ошибку, если в данных присутствуют NA.
na.omit, na.exclude — удалить все переменные со значением NA.
na.pass — оставить данные как есть.
Собственно, если нам надо подсчитать медиану возраста пассажиров Титаника, то можно поступить так:

median(train$Age, na.rm=T) [1] 28 

Второе решение, в общем-то, вытекает из первого — заранее очистить данные от всех неопределенных значений и больше не думать о них:

train.nopain <- na.omit(train) nrow(train.nopain) [1] 183 sum(is.na(train.nopain)) [1] 0 median(train.nopain$Age) [1] 36 

Как видим, вместе с NA испарились и 708 строк. Не страшно, если нас не особенно интересует результат, но стоит рассмотреть и другие варианты.

Используем подручные средства

Приемлемый вариант — заменить NA на какое-то заранее выбранное значение; обычно это среднее или медиана. Для этого напишем простую функцию:

simpleFix <- function(x, imputeFn=mean){   return(ifelse(is.na(x), imputeFn(x, na.rm=TRUE), x)) }  train.median <- train nas.idx <- which(is.na(train.median$Age)) train.median$Age <- simpleFix(train.median$Age, median) 

head(train.median$Age[nas.idx]) [1] 28 28 28 28 28 28 

Все неопределенные значения Age теперь заменены на медиану. Иногда есть смысл вычислять среднее значение или медиану с каким-то условием. Например, если мы корректируем NA поля Age пассажира мужского пола, который занимает каюту первого класса, то и среднее надо вычислять с учетом этого обстоятельства:

fixAge <- function(tdf, imputeFn=mean) {   tdf$Age[is.na(tdf$Age)] <- sapply(which(is.na(tdf$Age)),                                      function(i)                                        imputeFn(tdf$Age[tdf$Pclass == tdf$Pclass[i] &                                                         tdf$Sex == tdf$Sex[i]],                                                 na.rm=T))   return(tdf) }  nas.idx <- which(is.na(train$Age)) train.cond <- fixAge(train, median) 

head(train.cond$Age[nas.idx]) [1] 25.0 30.0 21.5 25.0 21.5 25.0 

Стоит отметить, что замена неопределенных значений средним или медианой приводит к уменьшению дисперсии данных. Данный метод можно несколько усовершенствовать (и усложнить), используя сингулярное разложение (SVD) и приближение матрицей меньшего ранга. При этом неопределенные значения сначала заменяются средними, потом исходная матрица аппроксимируется матрицей меньшего ранга, и неопределенные значения, которые были заменены на средние, заменяются на значения, взятые из разложения. Процедура аппроксимации повторяется несколько раз; для этого нам надо заранее задать ранг аппроксимации и количество шагов. Рассмотрим несколько сферический пример с матрицей случайных чисел размером 10х10 с 15 неопределенными значениями.

k <- 6 # Ранг аппроксимирующей матрицы n.iters <- 10 nrows <- 10 set.seed(100500) train.mat <- runif(nrows * nrows)  train.mat[sample(1:length(train.mat), 15)] <- NA train.mat <- matrix(train.mat, nrows) # Матрица случайных чисел с 15 NA значениями nas.idx <- which(is.na(train.mat)) train.svd <- train.mat train.svd <- apply(train.svd, 2, simpleFix) # Заменяем NA на среднее по столбцу for (i in 1:n.iters){   s <- svd(train.svd, k, k)   train.svd[nas.idx] <- (s$u %*% diag(s$d[1:k], nrow=k, ncol=k) %*%  t(s$v))[nas.idx] } 

head(train.svd[nas.idx]) [1] 0.3020229 0.4475467 0.3114711 0.7161445 0.4379184 0.6734933 

Что делать с нечисловыми значениями

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

fixSample <- function(x) {   x[is.na(x)] <- sample(x, sum(is.na(x)), replace = T)   return(x) }  set.seed(111) nas.idx <- which(is.na(train.cond$Embarked)) train.cond$Embarked <- fixSample(train.cond$Embarked) 

train.cond$Embarked[nas.idx] [1] S C Levels: C Q S  sum(is.na(train.cond$Embarked)) [1] 0 

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

Более универсальный подход

R обладает развитыми средствами для построения различных моделей — от простой линейной регрессии до техники «3B» (bagging, boosting, blending). Напишем функцию, которая использует RandomForest для вычисления неопределенных переменных; также уберем некоторые переменные PassengerId, Name, Ticket, Cabin (я, если честно, не нашел им достойного применения, а в поле Cabin столько неопределенных значений, что «довычисление» теряет смысл). Если пакет randomForest у вас не установлен, то командой install.packages("randomForest") R установит его из CRAN.

fixNA <- function(y, x) {   require(randomForest)   fixer <- randomForest(x[!is.na(y), ], y[!is.na(y)])   y[is.na(y)] <- predict(fixer, x[is.na(y), ])   return(y) }  set.seed(111) train.rf <- subset(train, select=-c(PassengerId, Name, Ticket, Cabin)) ageNA.idx <- which(is.na(train.rf$Age)) embNA.idx <- which(is.na(train.rf$Embarked)) sum(is.na(train.rf)) train.rf$Age <- fixNA(train.rf$Age, cbind(train.rf$Pclass, train.rf$Sex, train.rf$Parch)) train.rf$Embarked <- fixNA(train.rf$Embarked, cbind(train.rf$Pclass, train.rf$Sex, train.rf$Parch)) 

head(train.rf$Age[ageNA.idx]) [1] 29.65873 31.67546 26.64918 29.65873 26.64918 29.65873  head(train.rf$Embarked[embNA.idx]) [1] S S  sum(is.na(train.rf)) [1] 0 

В таком виде данные практически готовы для более детального анализа.

Заключение

Существует множество подходов и методов при работе с неполными данными — проблема далеко не тривиальная. В CRAN, в частности, есть и специализированные пакеты для обработки значений NA (например, Amelia, imputation), и пакеты, которые помимо прочего позволяют проводить манипуляции с NA (функции impute в пакете Hmisc, rfImpute в randomForest). Отдельного рассмотрения требует метод, основанный на алгоритме Еxpectation Maximization.

Ссылки и литература

1. Working with missing data
2. Data Imputation
3. Missing Data and Small-Area Estimation: Modern Analytical Equipment for the Survey Statistician. Nicholas T. Longford.