R語言：缺失值處理
前言
剛接觸缺失數據研究的讀者可能會被各式各樣的方法和言論弄得眼花繚亂。該領域經典的讀本是Little和Rubin的Statistical Analysis with Missing Data, Second Edition（2002）一書。其他比較優秀的專著還有Allison的Missing Data（2001）、Schafer和Graham的"Missing Data: Our View of the State of the Art"（2002），以及Schlomer、Bauman和Card的"Best Practices for Missing Data Management in Counseling Psychology"（2010）。
一個完整的處理方法通常包含以下幾個步驟：
(1) 識別缺失數據；
(2) 檢查導致數據缺失的原因；
(3) 刪除包含缺失值的實例或用合理的數值代替（插補）缺失值。
但遺憾的是，僅有識別缺失數據是最清晰明確的步驟。知道數據為何缺失依賴於你對數據生成過程的理解，而決定如何處理缺失值則需要判斷哪種方法的結果最為可靠和精確。
統計學家通常將缺失數據分為三類。它們都用概率術語進行描述，但思想都非常直觀。我們將用sleep研究中對做夢時長的測量（有12個動物有缺失值）來依次闡述三種類型。
(1) 完全隨機缺失 若某變量的缺失數據與其他任何觀測或未觀測變量都不相關，則數據為完全隨機缺失（MCAR）。若12個動物的做夢時長值缺失不是由於系統原因，那么可認為數據是MCAR。注意，如果每個有缺失值的變量都是MCAR，那么可以將數據完整的實例看做是對更大數據集的一個簡單隨機抽樣。
(2) 隨機缺失 若某變量上的缺失數據與其他觀測變量相關，與它自己的未觀測值不相關，則數據為隨機缺失（MAR）。例如，體重較小的動物更可能有做夢時長的缺失值（可能因為較小的動物較難觀察），“缺失”與動物的做夢時長無關，那么該數據就可以認為是MAR。此時，一旦你控制了體重變量，做夢時長數據的缺失與出現將是隨機的。
(3) 非隨機缺失 若缺失數據不屬於MCAR或MAR，則數據為非隨機缺失（NMAR）。例如，做夢時長越短的動物也更可能有做夢數據的缺失（可能由於難以測量時長較短的事件），那么數據可認為是NMAR。大部分處理缺失數據的方法都假定數據是MCAR或MAR。此時，你可以忽略缺失數據的生成機制，並且（在替換或刪除缺失數據后）可以直接對感興趣的關系進行建模。當數據是NMAR時，想對它進行恰當地分析比較困難，你既要對感興趣的關系進行建模，還要對缺失值的生成機制進行建模。
處理缺失數據的方法有很多，但不能保證都生成一樣的結果。下圖列出了一系列可用來處理不完整數據的方法，以及相應的R包。
目錄
1. 識別缺失數據
2. 探索缺失值模式
3. 理解缺失值數據的來由和影響
4. 處理缺失值完整實例分析
5. 處理缺失值的其他方法
主要程序包
install.packages(c("VIM","mice"))
library(VIM)
library(mice）
1. 識別缺失數據
首先，我們回顧一下前節的內容並地一步拓展。R使用NA（不可得）代表缺失值，NaN（不是一個數）代表不可能的值。另外，符號Inf和-Inf分別代表正無窮和負無窮。函數is.na()、is.nan()和is.infinite()可分別用來識別缺失值、不可能值和無窮值。每個返回結果都是TRUE或FALSE。表15-1給出了一些示例。
這些函數返回的對象與其自身參數的個數相同。若每個元素的類型檢驗通過，則由TRUE替換，否則用FALSE替換。例如，令y <- c(1, 2, 3, NA)，則is.na(y)返回向量c(FALSE, FALSE, FALSE,TRUE)。
函數complete.cases()可用來識別矩陣或數據框中沒有缺失值的行。若每行都包含完整的實例，則返回TRUE的邏輯向量；若每行有一個或多個缺失值，則返回FALSE。
對於識別缺失值，有兩點需要牢記。第一點，complete.cases()函數僅將NA和NaN識別為缺失值，無窮值（Inf和-Inf）被當做有效值。第二點，必須使用與本章中類似的缺失值函數來識別R數據對象中的缺失值。像myvar == NA這樣的邏輯比較無法實現。
NA：代表缺失值；
NaN：代表不可能的值；
Inf：代表正無窮；
-Inf：代表負無窮。
is.na（）：識別缺失值；
is.nan（）：識別不可能值；
is.infinite（）：無窮值。
is.na（）、is.nan（）和is.infinte（）函數的返回值示例
2. 探索缺失值模式
在決定如何處理缺失數據前，了解哪些變量有缺失值、數目有多少、是什么組合形式等信息非常有用。本節中，我們將介紹探索缺失值模式的圖表及相關方法。最后，如果知道了數據為何缺失，這將為后續深入研究提供許多啟示。
2.1 列表顯示缺失值
你已經學習了一些識別缺失值的基本方法。比如使用complete.cases()函數列出完整的實例，或者相反，列出含一個或多個缺失值的實例。但隨着數據集的增大，該方法就逐漸喪失了吸引力。此時你可以轉向其他R函數。
mice包中的md.pattern()函數可生成一個以矩陣或數據框形式展示缺失值模式的表格。將函數應用到sleep數據集，可得到：
> library(mice)
> data(sleep,package="VIM")
> md.pattern(sleep)
BodyWgt BrainWgt Pred Exp Danger Sleep Span Gest Dream NonD
42 1 1 1 1 1 1 1 1 1 1 0
2 1 1 1 1 1 1 0 1 1 1 1
3 1 1 1 1 1 1 1 0 1 1 1
9 1 1 1 1 1 1 1 1 0 0 2
2 1 1 1 1 1 0 1 1 1 0 2
1 1 1 1 1 1 1 0 0 1 1 2
2 1 1 1 1 1 0 1 1 0 0 3
1 1 1 1 1 1 1 0 1 0 0 3
0 0 0 0 0 4 4 4 12 14 38
表中1和0顯示了缺失值模式，0表示變量的列中有缺失值，1則表示沒有缺失值。第一行表述了“無缺失值”的模式（所有元素都為1）。第二行表述了“除了Span之外無缺失值”的模式。第一列表示各缺失值模式的實例個數，最后一列表示各模式中有缺失值的變量的個數。此處可以看到，有42個實例沒有缺失值，僅2個實例缺失了Span。9個實例同時缺失了NonD和Dream的值。數據集包含了總共(42 × 0) + (2 × 1) + … + (1 × 3) = 38個缺失值。最后一行給出了每個變量中缺失值的數目。
2.2 圖形探究缺失數據
雖然md.pattern()函數的表格輸出非常簡潔，但我通常覺得用圖形展示模式更為清晰。VIM包提供了大量能可視化數據集中缺失值模式的函數，本節我們將學習其中幾個：aggr()、matrixplot()和scattMiss()。aggr()函數不僅繪制每個變量的缺失值數，還繪制每個變量組合的缺失值數。
2.2.1 例1：使用函數aggr()函數繪圖。
library("VIM")
aggr(sleep,prop=FALSE,numbers=TRUE)
aggr(sleep, prop = TRUE, numbers = TRUE)
上述代碼的結果見下圖。
（VIM包將會打開GUI界面，你可以關閉它；本章我們使用代碼完成所有的工作。）可以看到，變量NonD有最大的缺失值數（14），有2個哺乳動物缺失了NonD、Dream和Sleep的評分。42個動物沒有缺失值。
代碼aggr(sleep, prop = TRUE, numbers = TRUE)將生成相同的圖形，但用比例代替了計數。選項numbers = FALSE（默認）刪去數值型標簽。
2.2.2 例2：使用函數matrixplot()函數繪圖。
matrixplot(sleep)
matrixplot()函數可生成展示每個實例數據的圖形。matrixplot(sleep)的圖形見下圖。
此處，數值型數據被重新轉換到[0, 1]區間，並用灰度來表示大小：淺色表示值小，深色表示值大。默認缺失值為紅色。
該圖形可以進行交互，單擊一列將會按其對應的變量重排矩陣。圖中的行便按BodyWgt降序排列。通過矩陣圖，你可以看出某些變量的缺失值模式是否與其他變量的真實值有關聯。此圖中可以看到，無缺失值的睡眠變量（Dream、NonD和Sleep）對應着較小的體重（BodyWgt）或腦重（BrainWgt）。
2.2.3 例3：使用函數marginplot()函數繪圖。
marginplot()函數可生成一幅散點圖，在圖形邊界展示兩個變量的缺失值信息。以做夢時長與哺乳動物妊娠期時長的關系為例，來看下列代碼：它的生成圖形見下圖。參數pch和col為可選項，控制着繪圖符號和使用的顏色。
library("VIM")
marginplot(sleep[c("Gest","Dream")],pch=c(20),col=c("darkgray","red","blue"))
圖形的主體是Gest和Dream（兩變量數據都完整）的散點圖。左邊界的箱線圖展示的是包含（深灰色）與不包含（紅色）Gest值的Dream變量分布。注意，在灰度圖上紅色是更深的陰影。四個紅色的點代表着缺失了Gest得分的Dream值。在底部邊界上，Gest和Dream間的關系反過來了。可以看到，妊娠期和做夢時長呈負相關，缺失妊娠期數據時動物的做夢時長一般更長。兩個變量均有缺失值的觀測個數在兩邊界交叉處（左下角）用藍色輸出。VIM包有許多圖形可以幫助你理解缺失數據在數據集中的模式，包括用散點圖、箱線圖、直方圖、散點圖矩陣、平行坐標圖、軸須圖和氣泡圖來展示缺失值的信息，因此這個包很值得探索。
2.3 用相關性探索缺失值
影子矩陣：用指示變量替代數據集中的數據（1表示缺失，0表示存在），這樣生成的矩陣有時稱作影子矩陣。
求這些指示變量間和它們與初始（可觀測）變量間的相關性，有且於觀察哪些變量常一起缺失，以及分析變量“缺失”與其他變量間的關系。
head(sleep)
str(sleep)
x<-as.data.frame(abs(is.na(sleep)))
head(sleep,n=5)
head(x,n=5)
y<-x[which(sd(x)>0)]
cor(y)
cor(sleep,y,use="pairwise.complete.obs")
3. 理解缺失值數據的來由和影響
識別缺失數據的數目、分布和模式有兩個目的：
（1）分析生成缺失數據的潛在機制；
（2）評價缺失數據對回答實質性問題的影響。
即：
（1）缺失數據的比例有多大？
（2）缺失數據是否集中在少數幾個變量上，抑或廣泛存在？
（3）缺失是隨機產生的嗎？
（4）缺失數據間的相關性或與可觀測數據間的相關性，是否可以表明產生缺失值的機制呢？
若缺失數據集中在幾個相對不太重要的變量上，則可以刪除這些變量，然后再進行正常的數據分析；
若有一小部分數據隨機分布在整個數據集中（MCAR），則可以分析數據完整的實例，這樣仍可得到可靠有效的結果；
若以假定數據是MCAR或MAR，則可以應用多重插補法來獲得有鏟的結論。
若數據是NMAR，則需要借助專門的方法，收集新數據，或加入一個相對更容易、更有收益的行業。
4. 處理缺失值完整實例分析
4.1 行刪除
函數complete.cases（）、na.omit（）可用來存儲沒有缺失值的數據框或矩陣形式的實例（行）：
# code1
newdata<-mydata[complete.cases(mydata),]
newdata<-na.omit(mydata)
# code2
options(digits=1)
cor(na.omit(sleep))
cor(sleep,use="complete.obs")
fit<-lm(Dream~Span+Gest,data=na.omit(sleep))
summary(fit)
4.2 多重插補
多重插補（MI）是一種基於重復模擬的處理缺失值的方法。
MI從一個包含缺失值的數據集中生成一組完整的數據集。每個模擬數據集中，缺失數據將使用蒙特卡洛方法來填補。
此時，標准的統計方法便可應用到每個模擬的數據集上，通過組合輸出結果給出估計的結果，以及引入缺失值時的置信敬意。
可用到的包Amelia、mice和mi包
mice() 函數首先從一個包含缺失數據的數據框開始，然后返回一個包含多個完整數據集的對象。每個完整數據集都是通過對原始數據框中的缺失數據進行插而生成的。
with() 函數可依次對每個完整數據集應用統計模型
pool() 函數將這些單獨的分析結果整合為一組結果。
最終模型的標准誤和p值都將准確地反映出由於缺失值和多重插補而產生的不確定性。
基於mice包的分析通常符合以下分析過程：
# 示例說明code
library(mice)
imp<-mice(mydata,m)
fit<-with(imp,analysis)
pooled<-pool(fit)
summary(pooled)
# mydata是一個飲食缺失值的矩陣或數據框；
# imp是一個包含m個插補數據集的列表對象，同時還含有完成插補過程的信息，默認的m=5
# analysis是一個表達式對象，用來設定應用於m個插補的統計分析方法。方法包括做線回歸模型的lm()函數、做廣義線性模型的glm()函數、做廣義可加模型的gam()、及做負二項模型的nbrm()函數。
# fit是一個包含m個單獨統計分析結果的列表對象；
# pooled是一個包含這m個統計分析平均結果的列表對象。
具體例子實現過程
> library(mice)
> data(sleep,package="VIM")
> imp <- mice(sleep,seed=1234)
iter imp variable
1 1 NonD Dream Sleep Span Gest
1 2 NonD Dream Sleep Span Gest
1 3 NonD Dream Sleep Span Gest
1 4 NonD Dream Sleep Span Gest
1 5 NonD Dream Sleep Span Gest
2 1 NonD Dream Sleep Span Gest
2 2 NonD Dream Sleep Span Gest
2 3 NonD Dream Sleep Span Gest
2 4 NonD Dream Sleep Span Gest
2 5 NonD Dream Sleep Span Gest
3 1 NonD Dream Sleep Span Gest
3 2 NonD Dream Sleep Span Gest
3 3 NonD Dream Sleep Span Gest
3 4 NonD Dream Sleep Span Gest
3 5 NonD Dream Sleep Span Gest
4 1 NonD Dream Sleep Span Gest
4 2 NonD Dream Sleep Span Gest
4 3 NonD Dream Sleep Span Gest
4 4 NonD Dream Sleep Span Gest
4 5 NonD Dream Sleep Span Gest
5 1 NonD Dream Sleep Span Gest
5 2 NonD Dream Sleep Span Gest
5 3 NonD Dream Sleep Span Gest
5 4 NonD Dream Sleep Span Gest
5 5 NonD Dream Sleep Span Gest
> fit <- with(imp,lm(Dream~Span+Gest))
> pooled <- pool(fit)
> summary(pooled)
est se t df Pr(>|t|) lo 95 hi 95 nmis fmi lambda
(Intercept) 2.546 0.255 10.0 52 1e-13 2.035 3.057 NA 0.09 0.05
Span -0.005 0.012 -0.4 52 7e-01 -0.029 0.020 4 0.09 0.05
Gest -0.004 0.001 -2.7 56 1e-02 -0.007 -0.001 4 0.05 0.02
> imp
Multiply imputed data set
Call:
mice(data = sleep, seed = 1234)
Number of multiple imputations: 5
Missing cells per column:
BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger
0 0 14 12 4 4 4 0 0 0
Imputation methods:
BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger
"" "" "pmm" "pmm" "pmm" "pmm" "pmm" "" "" ""
VisitSequence:
NonD Dream Sleep Span Gest
3 4 5 6 7
PredictorMatrix:
BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger
BodyWgt 0 0 0 0 0 0 0 0 0 0
BrainWgt 0 0 0 0 0 0 0 0 0 0
NonD 1 1 0 1 1 1 1 1 1 1
Dream 1 1 1 0 1 1 1 1 1 1
Sleep 1 1 1 1 0 1 1 1 1 1
Span 1 1 1 1 1 0 1 1 1 1
Gest 1 1 1 1 1 1 0 1 1 1
Pred 0 0 0 0 0 0 0 0 0 0
Exp 0 0 0 0 0 0 0 0 0 0
Danger 0 0 0 0 0 0 0 0 0 0
Random generator seed value: 1234
> imp$imp$Dream
1 2 3 4 5
1 1.0 0.5 0.5 0.5 0.3
3 2.6 2.1 1.5 1.8 1.3
4 3.4 3.1 3.4 1.2 3.4
14 0.3 0.5 0.5 0.3 1.2
24 1.8 1.3 3.6 0.9 5.6
26 2.3 3.1 2.0 2.6 2.1
30 1.2 0.3 3.4 2.6 2.3
31 3.4 0.5 0.6 1.0 0.5
47 0.5 1.5 1.5 2.2 3.4
53 0.3 0.5 0.5 0.5 0.6
55 0.5 0.9 2.6 2.7 2.4
62 1.0 2.1 0.5 3.9 3.6
>
>
> # 利用complete（）函數可觀察m個插補數據集中的任意一個，格式為：complete(imp,action=#)
> # eg:
>
> dataset3<-complete(imp,action=3)
> dataset3
BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger
1 7e+03 6e+03 2 0.5 3 39 645 3 5 3
2 1e+00 7e+00 6 2.0 8 4 42 3 1 3
3 3e+00 4e+01 11 1.5 12 14 60 1 1 1
4 9e-01 6e+00 13 3.4 16 2 25 5 2 3
5 3e+03 5e+03 2 1.8 4 69 624 3 5 4
6 1e+01 2e+02 9 0.7 10 27 180 4 4 4
7 2e-02 3e-01 16 3.9 20 19 35 1 1 1
8 2e+02 2e+02 5 1.0 6 30 392 4 5 4
9 3e+00 3e+01 11 3.6 14 28 63 1 2 1
10 5e+01 4e+02 8 1.4 10 50 230 1 1 1
11 4e-01 6e+00 11 1.5 12 7 112 5 4 4
12 5e+02 4e+02 3 0.7 4 30 281 5 5 5
13 6e-01 2e+00 8 2.7 10 18 46 2 1 2
14 2e+02 4e+02 3 0.5 3 40 365 5 5 5
15 7e-02 1e+00 6 2.1 8 4 42 1 1 1
16 3e+00 2e+01 9 0.0 9 50 28 2 2 2
17 8e-01 4e+00 7 4.1 11 6 42 2 2 2
18 2e-01 5e+00 10 1.2 11 10 120 2 2 2
19 1e+00 2e+01 5 1.3 6 34 28 1 2 1
20 6e+01 8e+01 12 6.1 18 7 21 1 1 1
21 5e+02 7e+02 11 0.3 11 28 400 5 5 5
22 3e+01 1e+02 3 0.5 4 20 148 5 5 5
23 1e-01 1e+00 11 3.4 14 4 16 3 1 2
24 2e+02 4e+02 8 3.6 12 39 252 1 4 1
25 8e+01 3e+02 5 1.5 6 41 310 1 3 1
26 4e+01 1e+02 11 2.0 13 16 63 1 1 1
27 1e-01 4e+00 10 3.4 14 9 28 5 1 3
28 1e+00 6e+00 7 0.8 8 8 68 5 3 4
29 5e+02 7e+02 2 0.8 3 46 336 5 5 5
30 1e+02 2e+02 7 3.4 11 22 100 1 1 1
31 4e+01 6e+01 3 0.6 4 16 33 3 5 4
32 5e-03 1e-01 8 1.4 9 3 22 5 2 4
33 1e-02 2e-01 18 2.0 20 24 50 1 1 1
34 6e+01 1e+03 6 1.9 8 100 267 1 1 1
35 1e-01 3e+00 8 2.4 11 13 30 2 1 1
36 1e+00 8e+00 8 2.8 11 4 45 3 1 3
37 2e-02 4e-01 12 1.3 13 3 19 4 1 3
38 5e-02 3e-01 11 2.0 13 2 30 4 1 3
39 2e+00 6e+00 14 5.6 19 5 12 2 1 1
40 4e+00 1e+01 14 3.1 17 6 120 2 1 1
41 2e+02 5e+02 8 1.0 8 24 440 5 5 5
42 5e-01 2e+01 15 1.8 17 12 140 2 2 2
43 1e+01 1e+02 10 0.9 11 20 170 4 4 4
44 2e+00 1e+01 12 1.8 14 13 17 2 1 2
45 2e+02 2e+02 6 1.9 8 27 115 4 4 4
46 2e+00 1e+01 8 0.9 8 18 31 5 5 5
47 4e+00 4e+01 11 1.5 12 14 63 2 2 2
48 3e-01 2e+00 11 2.6 13 5 21 3 1 3
49 4e+00 5e+01 7 2.4 10 10 52 1 1 1
50 7e+00 2e+02 8 1.2 10 29 164 2 3 2
51 8e-01 1e+01 6 0.9 7 7 225 2 2 2
52 4e+00 2e+01 5 0.5 5 6 225 3 2 3
53 1e+01 1e+02 2 0.5 3 17 150 5 5 5
54 6e+01 2e+02 3 0.6 4 20 151 5 5 5
55 1e+00 1e+01 8 2.6 11 13 90 2 2 2
56 6e-02 1e+00 8 2.2 10 4 100 3 1 2
57 9e-01 3e+00 11 2.3 13 4 60 2 1 2
58 2e+00 1e+01 5 0.5 5 8 200 3 1 3
59 1e-01 2e+00 13 2.6 16 2 46 3 2 2
60 4e+00 6e+01 10 0.6 10 24 210 4 3 4
61 4e+00 4e+00 13 6.6 19 3 14 2 1 1
62 4e+00 2e+01 18 0.5 19 13 38 3 1 1
5. 處理缺失值的其他方法
| 軟件包 | 描 述 |
| Hmisc | 包含多種函數，支持簡單插補、多重插補和典型變量插補 |
| mvnmle | 對多元正態頒數據中缺失值的最大似然估計 |
| cat | 對數線性模型中多元類別型變量的多重插補 |
| arrayImpute\arraryMissPattern、SeqKnn | 處理微陣列缺失值數據的實用函數 |
| longitudinalData | 相關的函數列表，比如對時間序列缺失值進行插補的一系列函數 |
| kmi | 處理生存分析缺失值的Kaplan-Meier多重插補 |
| mix | 一般位置模型中混合類別型和連續型數據的多重插補 |
| pan | 多元面板數據或聚類的多重插補 |