R語言中的缺失值處理

在處理一些真實(shí)數(shù)據(jù)時(shí)，樣本中往往會包含缺失值（Missing values)。我們需要對缺失值進(jìn)行適宜的處理，才能建立更為有效的模型，使得后續(xù)預(yù)測分析能有更小的偏差。本文將羅列不同的缺失值處理方法，并進(jìn)行具體應(yīng)用。

數(shù)據(jù)準(zhǔn)備和缺失模式設(shè)定

本文使用mlbench包中的BostonHousing數(shù)據(jù)集作為示例來演示不同的缺失值處理方法。由于原始的數(shù)據(jù)集并不包含缺失值，我們需要隨機(jī)刪除一些數(shù)據(jù)。通過這種方法，我們不僅可以評估由數(shù)據(jù)缺失帶來的精度損失，也可以比較不同處理方式的效果好壞。讓我們先加載這個(gè)數(shù)據(jù)集，并隨機(jī)刪除一些數(shù)據(jù)。

# 加載數(shù)據(jù)集
data ("BostonHousing", package="mlbench")
original <- BostonHousing  # backup original data

# 引入缺失值
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  

缺失值已經(jīng)生成好了，盡管我們已經(jīng)知道哪些位置的數(shù)據(jù)缺失，但還是用mice包中的md.pattern函數(shù)快速檢查下。

# 缺失值的模式
library(mice)
md.pattern(BostonHousing)  # 返回?cái)?shù)據(jù)的缺失值的模式

#>     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. 刪除觀測（記錄）

如果你的數(shù)據(jù)集擁有大量觀測，足以用來建立模型，那你可以把包含缺失值的觀測刪去（或者在建模時(shí)選擇不納入這些觀測，如設(shè)定na.action=na.omit）。在刪去相應(yīng)觀測后，請確保：

你有足夠的樣本點(diǎn)可以用來建模。

沒有引入偏差（譯者注：即認(rèn)為這些缺失值是隨機(jī)產(chǎn)生的，刪除對應(yīng)觀測后，樣本總體還是一個(gè)隨機(jī)樣本而非選擇樣本）。

# 例子
lm(medv ~ ptratio + rad, data=BostonHousing, na.action=na.omit)

2.刪除變量（字段）

如果某個(gè)變量包含大量的缺失值，我們可以直接刪除這個(gè)變量來保留更多的觀測，除非這個(gè)變量對于模型而言特別重要。應(yīng)用這個(gè)方法需要我們在變量的重要性和觀測的數(shù)量之間做權(quán)衡。

3.用均值、中位數(shù)或眾數(shù)插值

把缺失值用相應(yīng)變量的均值、中位數(shù)或眾數(shù)替換是一種比較粗糙的處理方法。其可行性也要取決于具體情境，如果變量的數(shù)值本身波動比較小或者對相應(yīng)變量的影響較小，使用這種粗略的插值法才可以得到使人滿意的結(jié)果。

library(Hmisc)
impute(BostonHousing$ptratio, mean)  # 均值替代
impute(BostonHousing$ptratio, median)  # 中位數(shù)替代
impute(BostonHousing$ptratio, 20)  # 用特殊值替代（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.用預(yù)測值插值

用預(yù)測值插值是一種比較前沿的方法，我們有很多模型可以實(shí)現(xiàn)這個(gè)過程，比如KNN插值，rpart還有mice。

4.1. KNN插值

DMwR包中的knnImputation函數(shù)會使用k近鄰方法來填補(bǔ)缺失值。具體流程如下：對于每個(gè)需要插值的觀測，先基于歐氏距離找到k個(gè)和它最近的觀測。再將這k個(gè)近鄰的數(shù)據(jù)利用距離逆加權(quán)得到插補(bǔ)值，最后用該值替代缺失值。

這種方式的優(yōu)勢在于你只要調(diào)用一次函數(shù)就能把所有缺失值插補(bǔ)好。該函數(shù)會把整個(gè)數(shù)據(jù)框作為參數(shù)，你不需要做其他設(shè)定。但在使用時(shí)請不要把響應(yīng)變量也一并輸入，因?yàn)樵谀銓y試集做處理時(shí)，你無法用未知的響應(yīng)變量來插值。

library(DMwR)
knnOutput <- knnImputation(BostonHousing[, !names(BostonHousing) %in% "medv"])  # 使用KNN插值.
anyNA(knnOutput)
#> FALSE

檢驗(yàn)該方法的精度

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個(gè)百分點(diǎn)?？傮w還不錯(cuò)。

4.2 rpart

利用knn插值的局限在于它對于因子類變量的插補(bǔ)效果可能不盡如人意。這種情況下rpart和mice就提供了更靈活的解決方案。rpart的優(yōu)勢是你只需要一個(gè)未缺失值就可以插補(bǔ)整個(gè)樣本。

插值思路是利用rpart（決策樹）替代knn來預(yù)測缺失值。對于因子類變量而言，我們在調(diào)用rpart函數(shù)式可以把method設(shè)為class（譯者注：即用分類樹)，數(shù)值型變量就設(shè)定method=anova（回歸樹）。當(dāng)然，我們也要避免把響應(yīng)變量傳入函數(shù)。

library(rpart)
class_mod <- rpart(rad ~ . - medv, data=BostonHousing[!is.na(BostonHousing$rad), ], method="class", na.action=na.omit)  # 因?yàn)閞ad是因子
anova_mod <- rpart(ptratio ~ . - medv, data=BostonHousing[!is.na(BostonHousing$ptratio), ], method="anova", na.action=na.omit)  # ptratio是數(shù)值變量
rad_pred <- predict(class_mod, BostonHousing[is.na(BostonHousing$rad), ])
ptratio_pred <- predict(anova_mod, BostonHousing[is.na(BostonHousing$ptratio), ])

ptratio的插補(bǔ)精度

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

與knn相比，mape值又額外下降了30%，可喜可賀。

rad的插補(bǔ)精度

actuals <- original$rad[is.na(BostonHousing$rad)]
predicteds <- as.numeric(colnames(rad_pred)[apply(rad_pred, 1, which.max)])
mean(actuals != predicteds)  # 計(jì)算誤分類比率
#> 0.25 

僅有25%的缺失值被誤分類，這個(gè)結(jié)果也不壞。

4.3 mice

mice是鏈?zhǔn)椒匠潭嘣逯档暮唽懀∕ultivariate Imputation by Chained Equations）。R中有個(gè)同名包提供了多種先進(jìn)的缺失值處理方法。它使用一種頗不常見的方法來進(jìn)行兩步插值：先利用mice函數(shù)建模再用complete函數(shù)生成完整數(shù)據(jù)。mice(df)操作會返回df的多個(gè)完整副本，每個(gè)副本都對缺失的數(shù)據(jù)插補(bǔ)了不同的值。complete()函數(shù)則會返回這些數(shù)據(jù)集中的一個(gè)（默認(rèn)）或多個(gè)。讓我們看看如何對rad和ptratio兩個(gè)變量插值：

library(mice)
miceMod <- mice(BostonHousing[, !names(BostonHousing) %in% "medv"], method="rf")  # 基于隨機(jī)森林模型進(jìn)行mice插值
miceOutput <- complete(miceMod)  # 生成完整數(shù)據(jù)
anyNA(miceOutput)
#> FALSE

計(jì)算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值與rpart相比又提升了48個(gè)百分點(diǎn)，亦可賽艇。

再看看rad的插值效果：

actuals <- original$rad[is.na(BostonHousing$rad)]
predicteds <- miceOutput[is.na(BostonHousing$rad), "rad"]
mean(actuals != predicteds)  # compute misclass error.
#> 0.15

誤分類比率降低到了15%，也就是說40個(gè)缺失觀測里插補(bǔ)錯(cuò)誤的只有6個(gè)。相較于rpart的錯(cuò)誤率（25%），這是一個(gè)了不起的提升。

如果你想了解的更深入，這里是mice包的 手冊 和DataScience+上另一篇關(guān)于mice包的 文章 。

盡管通過本文你已經(jīng)對各類處理方法有了初步了解，可這些還不足以幫助你判斷每種方法的優(yōu)劣。但當(dāng)你下次處理缺失值的時(shí)候，逐一測試這些方法是值得一試的。數(shù)據(jù)分析師培訓(xùn)