一、簡介
在實際工作中，遇到數據中帶有缺失值是非常常見的現象，簡單粗暴的做法如直接刪除包含缺失值的記錄、刪除缺失值比例過大的變量、用0填充缺失值等，但這些做法會很大程度上影響原始數據的分布或者浪費來之不易的數據信息，因此怎樣妥當地處理缺失值是一個持續活躍的領域，貢獻出眾多巧妙的方法，在不浪費信息和不破壞原始數據分布上試圖尋得一個平衡點，在R中用於處理缺失值的包有很多，本文將對最為廣泛被使用的mice和VIM包中常用的功能進行介紹，以展現處理缺失值時的主要路徑；
二、相關函數介紹
2.1 缺失值預覽部分
在進行缺失值處理之前，首先應該對手頭數據進行一個基礎的預覽：
1、matrixplot
效果類似matplotlib中的matshow，VIM包中的matrixplot將數據框或矩陣中數據的缺失及數值分布以色彩的形式展現出來，下面是利用matrixplot對R中自帶的airquality數據集進行可視化的效果：
rm(list=ls()) library(mice) library(VIM) #以airquality數據為例 data(airquality) data <- airquality #利用矩陣熱力圖查看缺失情況，紅色代表缺失 matrixplot(data)
紅色部分即代表數據缺失值所在位置，通過這個方法，可以在最開始對數據整體的缺失情況有一個初步認識，如通過上圖可以一眼看出變量Ozone缺失情況較為嚴重；
2、marginplot與marginmatrix
缺失值是否符合完全隨機缺失是在對數據進行插補前要着重考慮的事情，VIM中的marginplot包可以同時分析兩個變量交互的缺失關系，依然以airquality數據為例：
marginplot(data[c(1,2)])
如上圖所示，通過marginplot傳入二維數據框，這里選擇airquality中包含缺失值的前兩列變量，其中左側對應變量Solar.R的紅色箱線圖代表與Ozone缺失值對應的Solar.R未缺失數據的分布情況，藍色箱線圖代表與Ozone未缺失值對應的Solar.R未缺失數據的分布情況，下側箱線圖同理，當同一側紅藍箱線圖較為接近時可認為其對應考察的另一側變量缺失情況比較貼近完全隨機缺失，這種情況下可以放心大膽地進行之后的插補，否則就不能冒然進行插補；
與marginplot功能相似，marginmatrix在marginplot只能展現兩個變量的基礎上推廣到多個變量兩兩之間，效果類似相關性矩陣圖：
marginmatrix(data)
3、自編函數計算各個變量缺失比例
為了計算出每一列變量具體的缺失值比例，可以自編一個簡單的函數來實現該功能：
> #查看數據集中每一列的缺失比例 > miss.prop <- function(x){sum(is.na(x))/length(x)} > apply(data,2,miss.prop) Ozone Solar.R Wind Temp Month Day 0.24183007 0.04575163 0.00000000 0.00000000 0.00000000 0.00000000
通過自編的函數miss.prop，可以對每個變量中缺失值所占比例有個具體的了解；
2.2 mice函數
mice包中最核心的函數是mice()，其主要參數解釋如下：
data: 傳入待插補的數據框或矩陣，其中缺失值應表示為NA
m: 生成插補矩陣的個數，mice最開始基於gibbs采樣從原始數據出發為每個缺失值生成初始值以供之后迭代使用，而m則控制具體要生成的完整初始數據框個數，在整個插補過程最后需要利用這m個矩陣融合出最終的插補結果，若m=1，則唯一的矩陣就是插補的結果；
method: 這個參數控制了傳入數據框中每一個變量對應的插補方式，無缺失值的變量對應的為空字符串，帶有缺失值的變量默認方法為"pmm"，即均值插補
predictorMatrix: 因為mice中絕大部分方法是用擬合的方式以含缺失值變量之外的其他變量為自變量，缺失值為因變量構建回歸或分類模型，以達到預測插補的目的，而參數predictorMatrix則用於控制在對每一個含缺失值變量的插補過程中作為自變量的有哪些其他變量，具體用法下文示例中會詳細說明
maxit: 整數，用於控制每個數據框迭代插補的迭代次數，默認為5
seed: 隨機數種子，控制隨機數水平
在對缺失值插補過程中，非常重要的是為不同的變量選擇對應的方法，即method中對應的輸入，下表是每種算法對應的參數代號、適用數據類型和算法名稱：
方法代號 |
適用數值類型 | 對應的具體算法名稱 |
pmm |
any | Predictive mean matching |
midastouch |
any | Weighted predictive mean matching |
sample |
any | Random sample from observed values |
cart |
any | Classification and regression trees |
rf |
any | Random forest imputations |
mean |
numeric | Unconditional mean imputation |
norm |
numeric | Bayesian linear regression |
norm.nob |
numeric | Linear regression ignoring model error |
norm.boot |
numeric | Linear regression using bootstrap |
norm.predict |
numeric | Linear regression, predicted values |
quadratic |
numeric | Imputation of quadratic terms |
ri |
numeric | Random indicator for nonignorable data |
logreg |
binary | Logistic regression |
logreg.boot |
binary | Logistic regression with bootstrap |
polr |
ordered | Proportional odds model |
polyreg |
unordered | Polytomous logistic regression |
lda |
unordered | Linear discriminant analysis |
2l.norm |
numeric | Level-1 normal heteroscedastic |
2l.lmer |
numeric | Level-1 normal homoscedastic, lmer |
2l.pan |
numeric | Level-1 normal homoscedastic, pan |
2l.bin |
binary | Level-1 logistic, glmer |
2lonly.mean |
numeric | Level-2 class mean |
2lonly.norm |
numeric | Level-2 class normal |
在面對數據集具體情況時，對插補方法進行微調是很必要的步驟，在上面鋪墊了這么多之后，下面在具體示例上進行演示，並引入其他的輔助函數；
2.3 利用mice進行缺失值插補——以airquality數據為例
因為前面對缺失值預覽部分已經利用airquality進行演示，這里就不再贅述，直接進入正式插補部分，首先，我們將data傳入mice函數，注意這里設置maxit為0以取得未開始迭代的初始模型參數：
#初始化插補模型,這里最大迭代次數選0是為了取得未開始插補的朴素模型參數 init <- mice(data, maxit = 0)
下面我們來看看取得的需要進行調整的重要參數的初始情況：
> #取得對於每一個變量的初始插補方法 > methods <- init$method > methods Ozone Solar.R Wind Temp Month Day "pmm" "pmm" "" "" "" ""
可以看到對應缺失變量Ozone和Solar.R的插補擬合方法為pmm，下面我們把它們改成CART決策樹回歸：
#將變量Ozone的插補方法從pmm修改為norm methods[c("Ozone")] <- 'cart' #將變量Solar.R的插補方法從pmm修改為norm methods[c("Solar.R")] <- 'cart'
接着我們來查看predictorMatrix參數：
> #取得對每一個變量進行擬合用到的變量矩陣，0代表不用到，1代表用到
> predM <- init$predictorMatrix
> predM
Ozone Solar.R Wind Temp Month Day
Ozone 0 1 1 1 1 1
Solar.R 1 0 1 1 1 1
Wind 1 1 0 1 1 1
Temp 1 1 1 0 1 1
Month 1 1 1 1 0 1
Day 1 1 1 1 1 0
這里我們認為變量Month和Day是日期，與缺失變量無相關關系，因此將其在矩陣中對應位置修改為0使它們不參與擬合過程：
#調整參與擬合的變量 #這里認為日期對與其他變量無相關關系，因此令變量Month與變量Day不參與對其他變量的擬合插補過程 predM[, c("Month")] <- 0 predM[c("Month"),] <- 0 predM[, c("Day")] <- 0 predM[c("Day"),] <- 0
這樣我們就完成了兩個重要參數的初始化，下面我們進行正式的擬合插補：
#利用修改后的參數組合來進行擬合插補
imputed <- mice(data, method = methods, predictorMatrix = predM)
隨着程序運行完，我們需要的結果便呼之欲出，但在取得最終插補結果前，為了嚴謹起見，需要對模型的統計學意義進行分析，下面以Ozone為例：
1、查看模型中Ozone對應的擬合公式：
> #查看Ozone主導的擬合公式 > imputed$formulas['Ozone'] $Ozone Ozone ~ 0 + Solar.R + Wind + Temp <environment: 0x000000002424d410>
可以看到，Ozone對應的公式與前面predictorMatrix參數中經過修改的保持一致；
2、基於上述公式為合成出的m=5個數據框分別進行擬合：
> #把上面的公式填入下面的lm()內 > fit <- with(imputed, lm(Ozone ~ Solar.R + Wind + Temp)) > > #查看fit中對應每一個插補數據框的回歸顯著性結果 > fit call : with.mids(data = imputed, expr = lm(Ozone ~ Solar.R + Wind + Temp)) call1 : mice(data = data, method = methods, predictorMatrix = predM) nmis : Ozone Solar.R Wind Temp Month Day 37 7 0 0 0 0 analyses : [[1]] Call: lm(formula = Ozone ~ Solar.R + Wind + Temp) Coefficients: (Intercept) Solar.R Wind Temp -53.36269 0.05791 -3.37688 1.51550 [[2]] Call: lm(formula = Ozone ~ Solar.R + Wind + Temp) Coefficients: (Intercept) Solar.R Wind Temp -55.70273 0.06189 -3.25456 1.50816 [[3]] Call: lm(formula = Ozone ~ Solar.R + Wind + Temp) Coefficients: (Intercept) Solar.R Wind Temp -55.45601 0.05659 -2.90911 1.47244 [[4]] Call: lm(formula = Ozone ~ Solar.R + Wind + Temp) Coefficients: (Intercept) Solar.R Wind Temp -70.00005 0.07008 -2.56784 1.59856 [[5]] Call: lm(formula = Ozone ~ Solar.R + Wind + Temp) Coefficients: (Intercept) Solar.R Wind Temp -72.08381 0.05252 -2.71741 1.67060
3、將上述5個擬合模型融合並查看融合后的顯著性水平
> #計算fit中模型的顯著性 > pooled <- pool(fit) > summary(pooled) estimate std.error statistic df p.value (Intercept) -61.32105668 21.84851057 -2.806647 53.60143 6.964478e-03 Solar.R 0.05979673 0.02144304 2.788632 90.71532 6.449107e-03 Wind -2.96516062 0.67509037 -4.392243 29.06277 1.361726e-04 Temp 1.55305117 0.23531131 6.599985 78.14524 4.468489e-09
可以看到從截距項，到每一個變量的p值都遠遠小於0.05，至少在0.05顯著性水平下每個參數都具有統計學意義；
4、對5個合成出的數據框在缺失值位置進行融合，這里需要用到新的函數complete，其主要有下面三個參數：
data: 前面mice函數輸出的結果
action: 當只希望從合成出的m個數據框中取得某個單獨的數據框時，可以設置action參數，如action=3便代表取得m個數據框中的第3個
mild: 邏輯型變量，當為TRUE時，會輸出包含全部m個合成數據框的列表
獲悉上列參數意義后，若只想抽取某個數據框如第3個：
result <- complete(imputed, action = 3) matrixplot(result)
可以看到，取回第3個數據框，每個缺失值都已被補全，若希望得到5個合成數據框的融合結果，則需要自編函數：
#取得所有合成數據框組成的列表 complete(imputed, mild = T) all.imputed <- complete(imputed, action = 'all') #對得到的m個合成數據框缺失值部分求平均 avgDF <- function(data,all.imputed){ baseDF <- all.imputed[[1]][is.na(data)] for(child in 2:length(all.imputed)){ baseDF <- baseDF + all.imputed[[child]][is.na(data)] } data[is.na(data)] <- baseDF/length(all.imputed) return(data) } #得到最終平均數據框 result <- avgDF(data,all.imputed) matrixplot(result)
以上就是本文的全部內容，如有錯誤之處望斧正。
參考資料：
https://stefvanbuuren.name/Winnipeg/Lectures/Winnipeg.pdf
https://www.rdocumentation.org/packages/mice/versions/3.5.0/topics/mice