Topic 17. 临床预测模型之缺失值识别及可视化

桓峰基因公众号推出基于R语言临床预测模型构建方法教程并配有视频在线教程，目前整理出来的教程目录如下：

Topic 1. _临床_标志物生信分析常规思路
Topic 2. 生存分析之 Kaplan-Meier
Topic 3. SCI文章第一张表格–基线表格
Topic 4. _临床_预测模型构建 Logistic 回归
Topic 5. 样本量确定及分割
Topic 6 计数变量泊松回归
Topic 7. _临床_预测模型–Cox回归
Topic 8. _临床_预测模型-Lasso回归
Topic 9. SCI 文章第二张表—单因素回归分析表
Topic 10. 单因素 Logistic 回归分析—单因素分析表格
Topic 11. SCI中多元变量筛选—单/多因素表
Topic 12 _临床_预测模型—列线表 (Nomogram)
Topic 13. _临床_预测模型—一致性指数 (C-index)
Topic 14. _临床_预测模型之校准曲线 (Calibration curve)
Topic 15. _临床_预测模型之决策曲线 (DCA)
Topic 16. _临床_预测模型之接收者操作特征曲线 (ROC)
Topic 17. 临床预测模型之如何处理缺失值——识别并可视化

前言

缺失值是一个在任何研究中几乎都存在的问题，你进行任何的调查、测量，总是不可能保证所有数据齐全。目前关于缺失值的研究已经发展出一个专门的领域，处理缺失值的方法也有很多,这里主要介绍几种实际中常见的方式。在统计软件中，如果你忽略缺失值，那么它会自动把缺失的观测删除，不管这个观测是有1个变量缺失还是有10个变量缺失。所以，如果你想假装看不见，那是不行的，必须采取一定的处理措施。

概念

不过在介绍缺失值处理技术之前，我们需要先了解几个与缺失有关的概念。

1．完全随机缺失(Missing Completely at Random，MCAR）

完全随机缺失的意思是：缺失的数据与自身和其他任何变量都没有关系。例如，某研究调查了收入、教育程度等变量，如果收入有缺失，而且这种缺失与收入本身无关(不管是收入高的人还是收入低的人，都有同样的缺失率)，与其他变量也无关（如不管教育程度是高还是低，都有同样的缺失率)，那么这种情况称为完全随机缺失，也就是说，任何人都有相同的机会产生缺失。在这种情况下，把缺失数据直接删除不会影响结果估计的准确性，只会影响精确性。例如，在不缺失的情况下系数估计值为0.6，在删除缺失值后系数估计值仍为0.6, 只是由于例数变少，标准误会增大，从而使得置信区间变宽。

2．随机缺失(Missing at Random,MAR）

完全随机缺失是一种最理想的情况，然而在实际中往往很难保证这一假定，更为常见的是随机缺失。随机缺失的意思是：缺失变量与自身无关，但与其他变量有关。如收入与教育程度，如果收入缺失与教育程度有关系（如教育程度高的人比教育程度低的人缺失得更多），但是与收入本身无关（高收入和低收入人群的缺失人数差不多）, 那么这种情况就是随机缺失。

3．非随机缺失(Not Missing at Random,NMAR）

非随机缺失是指缺失与自身变量有关，如收入的缺失，如果发现收入高的人更倾向于不填数据，而收入低的人一般不缺失，则说明收入的缺失是与自身变量有关的, 这种情况就是非随机缺失。

缺失值的处理函数

1. 常用缺失值处理函数

缺失数值在R中用符号“NA”来表示。当将一个含空单元格的Excel表格导入R控制台时,这些空的单元格将被NA替代。这不同于STATA用“.”替代空单元格。在R中数值变量及字符变量使用相同的缺失值符号。R提供了一些函数来处理缺失值(表1)。

is.na()是用来判断元素是不是NA类型最常用的方法。它返回一个长度和传入参数长度相同且所有数据都是逻辑值(FALSE 或TRUE)的对象。假设我们有6个患者,但只记录了5个乳酸值,有1个缺失了。

lactate <- c(0.2, 3.3, 4.5, NA, 6.1, 2.4)
is.na(lactate)
## [1] FALSE FALSE FALSE  TRUE FALSE FALSE
which(is.na(lactate))
## [1] 4

is.na()返回向量长度为6，第4位的值为TRUE，表示第4位患者乳酸值缺失。有人可能会使用逻辑检验(如lactateNA)去检测缺失数据,该方法不可能返回TRUE,因为缺失值是不能比较的,你必须用缺失值函数。“”操作符返回的值都是NA。通过使用which()函数,你可以找到哪个元素的向量包含 NA。在这个例子中which()函数返回4,表明第4位患者乳酸值丢失。接下来,你可能想要描述6个患者的乳酸水平。在统计描述中通常用到均值、方差和标准差。

mean(lactate)
## [1] NA
sum(lactate)
## [1] NA
var(lactate)
## [1] NA
sd(lactate)
## [1] NA

因为向量lactate包含缺失值，所以这些函数均返回NA，好在应用统计函数na.rm()能够移除NAs。

mean(lactate, na.rm = TRUE)
## [1] 3.3
sd(lactate, na.rm = TRUE)
## [1] 2.219234

上述的计算中我们加了“na.rm=TRUE”,表示将缺失数据移除后进行统计描述，结果5个患者的非缺失乳酸值用于计算平均值和标准偏差，而这些结果正是你想要的。

2. 含有缺失数据的数据框

在实际环境下，你会在数据框中遇到数据缺失。因此，本次将重点介绍如何处理数据框中的数据缺失。首先我们创建了一个包含3个变量和5个观察值的数据框。

ptid <- c(1, 2, 3, 4, 5)
sex <- c("m", "f", NA, "f", "m")
lactate <- c(0.2, 3.3, 4.5, NA, 6.1)
data <- data.frame(ptid, sex, lactate)
data
##   ptid  sex lactate
## 1    1    m     0.2
## 2    2    f     3.3
## 3    3 <NA>     4.5
## 4    4    f      NA
## 5    5    m     6.1

注意第3个患者性别缺失，第5个患者乳酸值缺失

na.fail(data)
#Error in na.fail.default(data) : missing values in object

上面语句中的na.fail()函数可以在数据集中检测出缺失数据。如果没有缺失数据,它会返回指定对象。如果有缺失数据，它返回一个错误信息，指出数据集中有一个或更多个缺失值。虽然回归模型中一些好的默认设置能有效忽略缺失数据，但创建一个新的排除了缺失数据的数据框也是有用的。

na.omit(data)
##   ptid sex lactate
## 1    1   m     0.2
## 2    2   f     3.3
## 5    5   m     6.1

上述的na.omit()返回的是移除了缺失值的新数据框。可以看到，含有缺失数据的第3和第4个患者在变量中被移除了。另外,通过使用以下的代码可以达到相同的目的。

complete.data <- data[complete.cases(data), ]
complete.data
##   ptid sex lactate
## 1    1   m     0.2
## 2    2   f     3.3
## 5    5   m     6.1

在以上代码中我们用complete.cases()函数获得了一个不包含NA的新数据框，当使用之前提到过的na.fail()函数进行检验时，返回了一个完整的complete.data数据框，而不再是之前显示的错误信息提示。

na.fail(complete.data)
##   ptid sex lactate
## 1    1   m     0.2
## 2    2   f     3.3
## 5    5   m     6.1

有时你可能想查找哪个患者的变量包含缺失数据，这时可以尝试前面介绍的函数is.na()。

is.na(data)
##       ptid   sex lactate
## [1,] FALSE FALSE   FALSE
## [2,] FALSE FALSE   FALSE
## [3,] FALSE  TRUE   FALSE
## [4,] FALSE FALSE    TRUE
## [5,] FALSE FALSE   FALSE
which(is.na(data))
## [1]  8 14

返回的8、14是缺失值的位置。R计算数据位置以列优先。当数据集不大时，这当然没有问题。如果这个数据集很大而且含有很多缺失值时，返回的一连串关于数据框位置的号码可能对你没有什么意义。下面的代码可以帮助你找出至少含有一个缺失数据的观察对象(行)。

unique(unlist(lapply(data, function(x) which(is.na(x)))))
## [1] 3 4

正如预期的那样,返回的结果表明第3和第4个患者至少含有一个缺失数据。函数lapply()返回了变量ptid、sex和llactate一系列缺失数据的位置。unlist()简化了list结构返回一个向量，lapply()函数返回的对象的属性为list。当一个观察对象(行)含有多个缺失值时，unique()函数只关心观测值中的一次数据缺失。有些情况下，你可能会对某个变量的缺失值数量(一列中缺失值的数量)感兴趣，从而将大于一定比例缺失值的变量排除。

unlist(lapply(data, function(x) sum(is.na(x))))/ncol(data)
##      ptid       sex   lactate
## 0.0000000 0.3333333 0.3333333
ptid
## [1] 1 2 3 4 5

上述语句提示，ptid没有缺失值，sex和lactate含有33%的缺失值。

3. 缺失值的回归模型

统计分析中建立回归模型是数据分析中较为常用的方法。然而，在回归模型拟合中缺失值的细节总是被忽略。少量数据的缺失，按照函数的默认处理方式是较为可靠的。但当出现在大量数据缺失时问题可能就会出现了，分析人员必须了解在建模中如何处理缺失数据。这里我将说明na.omit()和na.exclude()的差异。为了利于说明，我将用患者的id进行乳酸的回归分析。当然，这没有任何实际意义。

model.omit <- lm(ptid ~ lactate, data = data, na.action = na.omit)
model.exclude <- lm(ptid ~ lactate, data = data, na.action = na.exclude)
resid(model.omit)
##          1          2          3          5
##  0.3895652 -0.6052174 -0.3773913  0.5930435
resid(model.exclude)
##          1          2          3          4          5
##  0.3895652 -0.6052174 -0.3773913         NA  0.5930435
fitted(model.omit)
##         1         2         3         5
## 0.6104348 2.6052174 3.3773913 4.4069565
fitted(model.exclude)
##         1         2         3         4         5
## 0.6104348 2.6052174 3.3773913        NA 4.4069565

4. 一些关于缺失值的高级函数

有些情况下当你不想简单地删除缺失数据，或缺失数据可能具有特殊的临床意义时，我们需要使用一些关于缺失值的高级函数。在我们所举的例子中，缺失的乳酸值可能表示患者休克已经完全纠正(乳酸是组织灌注不足、缺氧的生物指标，对于血流动力学稳定的患者，临床医生通常不会去查乳酸值)。因此，乳酸NA显示病情稳定的患者。因为NA包含重要的信息，所以添加新的类别(称为“NA”的新的类别)来代替原有缺失值是明智的。你可以试试tree程序包中的na.tree.replace()函数来达到这种目的，虽然这个函数还有一定的局限性，因为它只能用于离散变量缺失值。我创建一个新的数据框进行说明。

if (!require(tree)) install.packages("tree")library(tree)
data.discrete <- data.frame(ptid = c(1, 2, 3, 4, 5), lactate = c("low", NA, "moderate",NA, "high"), death = c("y", "y", NA, "n", "y"))
data.discrete
na.fail(data.discrete)
# Error in na.fail.default(data) : missing values in object
newdata.discrete <- na.tree.replace(data.discrete)
# Error in na.tree.replace(data.discrete) : continuous variable lactate
# contained NAs
na.fail(newdata.discrete)
# Error in na.fail.default(data.discrete) : missing values in object

数据框data.discrete包含三个变量和最后两个为含有NA的离散变量。利用函数na.fail()返回的错误提示表明data.discrete数据框中含有缺失数据。利用函数na.tree.replace()，一个名为newdata.discrete的数据框被创建出。接着我们可以看到函数na.fail()返回了一个新的数据框。正如你所看到的,缺失值被字符串值“NA”所取代。na.tree.replace()的不足之处在于当任何一个连续变量包含NA，这个函数将停止运算。注意,因连续变量中包含NA，将返回一条错误信息。这个问题可以用gam程序包中的函数na.gam.replace()解决。注意安装gam程序包同时需安装foreach程序包。

if (!require(gam)) install.packages("gam")library(foreach)
library(gam)
newdata <- na.gam.replace(data)
newdata
##   ptid  sex lactate
## 1    1    m   0.200
## 2    2    f   3.300
## 3    3 <NA>   4.500
## 4    4    f   3.525
## 5    5    m   6.100
# na.fail(newdata) Error in na.fail.default(newdata) : missing values in object

注意：newdata通过函数na.fail()的检测,返回新的数据集。变量sex中的缺失值被字符串值“NA”所取代，变量 lactate中的缺失值被非缺失乳酸值的平均值所取代。

缺失数据的可视化方法

之前的章节对如何使用R语言来处理缺失数据进行了初步探讨，但有时候这些功能并不能满足复杂数据类型的处理，R语言为缺失数据的处理提供了许多高级的解决方案，本文就着重介绍一下这些方案，其中重点讲述了缺失数据的可视化方法"。本文首先建立了一个包含五个变量的数据框，其中三个变量含有缺失数据，且三个变量分别代表不同的数据缺失类型，接着我们讲述如何对这些数据进行探索研究。

1. 缺失数据分类

统计学家通常将缺失数据分为三类。完全随机缺失(MCAR)指一个变量缺失值的存在与其他任何观察到的和未观察到的变量均无关。换而言之，该缺失模式没有任何规律可循。随机缺失(MAR)指一个变量缺失值的存在与其他观察到的变量相关，而与未观测到的变量无关。非随机缺失(NMAR)指一个变量缺失值既不是MCAR，也不是MAR，即缺失值产生与观察到的变量也有关。例如，一个血流动力学稳定的患者一般会较少检测其乳酸值，因为乳酸值在这类患者中一般是正常的，这就导致了乳酸的缺失值与其自身值有关，而该自身值恰恰为未观察到的数据。

2. 模拟数据框

模拟建议一个包含200个观测值的数据框。该数据框仅仅用来说明R语言是如何工作的，并没有任何临床意义。以下是该数据框包含的五个变量：年龄(age)、性别(sex)、乳酸值(Iac)、白细胞数(wbc)和C反应蛋白(crp)。在每次模拟中，均设定一个seed，使读者可以复制结果。

set.seed(123456)
age <- round(abs(rnorm(200, mean = 67, sd = 19)))
set.seed(12345)
sex <- rbinom(200, 1, 0.45)
set.seed(12356)
sex.miss.tag <- rbinom(200, 1, 0.3)  #MCAR
sex[sex.miss.tag == 1] <- NA
sex[sex == 1] <- "male"
sex[sex == 0] <- "female"
set.seed(12456)
lac <- round(abs(rnorm(200, mean = 3, sd = 4)), 1)
set.seed(13456)
lac.miss.tag <- rbinom(200, 1, 0.3)
lac[lac <= 3 & lac.miss.tag == 1] <- NA# NMAR
set.seed(23456)
wbc <- round(abs(rnorm(200, mean = 12, sd = 4)), 1)
set.seed(123)
wbc.miss.tag <- rbinom(200, 1, 0.3)
wbc[wbc.miss.tag == 1] <- NA
set.seed(1234)
crp <- round(abs(rnorm(200, mean = 50, sd = 100)), 1)
set.seed(3456)
crp.miss.tag <- rbinom(200, 1, 0.4)
crp[wbc <= 128 & crp.miss.tag == 1] <- NA# MAR
data <- data.frame(age, sex, lac, wbc, crp)

所有观测值中年龄是完整的,不包含缺失数据。假设人口的平均年龄为67岁，而标准差为19。abs()函数用来避免负数出现。round()函数用于保留一定位数的有效数字。性别变量属于二分类变量，我们假设其分布遵从二项式分布。性别变量含有缺失数据，我们设定其缺失类型为MCAR。lac服从正态分布,平均值为3，标准差为4，其缺失类型设定为NMAR，我们假定其缺失值发生在Lac≤3时可能性更大。wbc呈正态分布,以及其缺失值属于MCAR。crp呈正态分布，其缺失值发生在wbc≤12时更多，属于MAR。crp缺失数据产生的基本原理是，在临床试验中，医师们首先会测定白细胞数值，对于白细胞数高的，他们才将会进一步测定crp，这样就会导致crp的缺失依赖白细胞计数。

3. 使用md.pattern()函数探索缺失数据模式

mice程序包中的md.pattern()函数可用于制作显示缺失模式的表格。

if (!require(mice)) install.packages("mice")library(mice)
md.pattern(data)

##    age lac crp wbc sex
## 42   1   1   1   1   1   0
## 39   1   1   1   1   0   1
## 32   1   1   1   0   1   1
## 16   1   1   1   0   0   2
## 36   1   1   0   1   1   1
## 12   1   1   0   1   0   2
## 7    1   0   1   1   1   1
## 5    1   0   1   1   0   2
## 4    1   0   1   0   1   2
## 3    1   0   1   0   0   3
## 3    1   0   0   1   1   2
## 1    1   0   0   1   0   3
##      0  23  52  55  76 206

在输出表格主体中，“1”代表非缺失值，“0”代表缺失值。第一纵列显示了独特缺失数据模式的数目。在我们的例子中,有58个观察对象不含缺失数据，42个观察对象仅性别变量有缺失数据。最右侧纵列显示了在特定缺失模式下缺失变量的个数。例如，在第一行中没有缺失值,则显示为“0”。最后一行统计了每个变量的缺失值数目。例如，年龄变量没有缺失值，显示“o”。而crp变量有33个缺失值。在研究中可能会剔除一些缺失值较多的变量，此时这个表格就可以提供有用的参考信息。

4. 缺失数据模式的视觉呈现

尽管以上表格简洁有效地展示了缺失模式，但人们仍希望以图形来显示它。VIM程序包内有功能强大的用于直观显示缺失数据模式的函数，这些函数对于研究缺失值或插补值有帮助。研究缺失数据模式对选择一个合适插补方式来估计缺失值是非常有必要的。因此可视化工具应在插补运算之前进行，通常在缺失数据插补后还会进行诊断以明确该插补值是否合理。能进行缺失数据可视化的函数有以下三个：matrixplot()， scattMiss()和aggr()。

if (!require(VIM)) install.packages("VIM")
library(VIM)
matrixplot(data)

matrixplot()函数提供了与用户的交互的功能。当“Click in a column to sortby the corresponding variable.”的消息弹出时，单击wbc纵列，结果显示如图。在该图形中，缺失值标记成红色。连续的变量重新调整，标记为灰色。浅色用于标记较小值，深色用于标记较大值。从图中可以看出crp的缺失值仅出现在低白细胞水平时，这与之前建立该数据框的规则相一致。因为选择了wbc变量来排序，所以首先显示的是缺失值，非缺失值根据wbc大小降序排列。

barMiss(data)  #similar function histMiss(data)

nrow(data[lac <= 1 & !is.na(lac), ])
## [1] 23
table(complete.cases(data[lac <= 1 & !is.na(lac), ]))
##
## FALSE  TRUE
##    20     3
table(complete.cases(data[is.na(lac), ][c("age", "sex", "wbc", "crp")]))
##
## FALSE  TRUE
##    16     7

barMiss()函数和histMiss()函数生成的图形是一致的，以下仅对barMiss()K数进行说明。默认情况下，barMiss()函数同样提供了交互功能，点击图形边缘(左右)可选择要显示的变量。接下来，将我们对lac变量进行展示。，上图图展示了用条形图来突出其他变量的缺失分布情况，其水平轴是lac值。每个条形图显示了除lac变量以外其他变量的缺失情况。此外，变量lac缺失值的信息在右边的条形图单独显示。接下来我们对该图进一步分析。当lac≤1时有23组观测值。其中，有18组存在其他变量的缺失，其余5组完整。右边的条形图显示了变量lac有23个缺失值。其中，有16组存在其他变量的缺失，其余7组完整。

aggr(data, numbers = TRUE, prop = FALSE)

aggr()函数产生的缺失数据模式显示在上图。左边的显示板展示了每个变量缺失值的个数。正如预期的那样，变量年龄没有缺失值，而lac的缺失值有20多个。右边的显示板表达了与md.pattern()函数生成的表格一样的信息。有58组完整的观测值,并没有数据缺失,有42组观测值仅存在性别变量缺失。

marginplot(data[c("wbc", "crp")], pch = c(20), col = c("green", "red", "blue"))

marginplot()函数的结果显示在上图。非缺失值标记为绿色，缺失值标记为红色。crp有33个缺失值，且相对CRP缺失值的wbc平均值为9左右，而crp值完整的wbc值平均值大约为13左右(对比横列的红、绿色箱线图)，显而易见，crp值缺失依赖于wbc值的大小。反过来，wbc值缺失并不受crp值影响(对比纵列的红、绿色箱线图)。

marginmatrix(data)

marginmatrix()函数是marginplot()函数的扩展如下图，只不过我们将所有的变量两两之间的关系用矩阵的形式表示，其中每幅图的解读和上图一样。

spineMiss(data)

spineMiss()函数生成的图形与barMiss()函数生成的相似。棘状图通过将每个单元格分成两部分来突出其他变量的缺失值上图。此外，关于指定变量(lac)的缺失值信息显示在右侧。该函数的纵轴以比例的形式来替代barMiss()函数条形图上的计数值。

scattmatrixMiss(data)

scattmatrixMiss()函数通过已知变量(age、 sex、 lac、wbc和lcrp)中被标记的缺失值来生成散点图矩阵上图。变量被标记的缺失值可以通过单击对角面板来添加或删除。这些对角面板展示了非标记观测值及标记观测值的密度图。红十字符号代表任一变量(年龄、性别、wbc、lac、crp)的缺失值。

5. 运用相关系数矩阵研究缺失数据模式

相关系数矩阵运用于研究当两个变量倾向于同时具有缺失值,，或某一变量缺失值的产生与另一变量观测值之间存在联系。为了完成该类对比，需要建立一个影子矩阵，所谓影子矩阵就是指用“1”代表缺失值，“o”代表非缺失值。如下一行语句用is.na()函数来判断对象数据框是否含缺失值,其中缺失值返回“1”，否则返回“o”。

shadow <- as.data.frame(abs(is.na(data)))

接着，创建一个新的数据框,该数据框仅保留含有缺失数据的变量。

miss.shadow <- shadow[, which(unlist(lapply(shadow, sum)) != 0)]
round(cor(miss.shadow), 3)
##        sex    lac    wbc    crp
## sex  1.000  0.008 -0.044 -0.159
## lac  0.008  1.000  0.024 -0.071
## wbc -0.044  0.024  1.000 -0.365
## crp -0.159 -0.071 -0.365  1.000

以上表格中这些变量没有很强的相关性,可以放心地得出以下结论：一个变量缺失值的产生与另一变量缺失值并无关联。接着，我们来研究某一变量缺失值产生和另一变量观察值之间是否存在关联。在运行cor()函数之前，cor()函数的第一个赋值对象仅仅保留连续变量，即将sex剔除，因为我们研究缺失数据是否依赖于其他变量的值，该命题要求其他变量为连续变量。再次使用round()函数使输出更加简洁。

round(cor(data[!names(data) %in% c("sex")], miss.shadow, use = "pairwise.complete.obs"),3)
##        sex    lac   wbc    crp
## age -0.031  0.055 0.082 -0.040
## lac -0.101     NA 0.163 -0.003
## wbc -0.073 -0.115    NA -0.088
## crp -0.019 -0.027 0.024     NA

通过以上表格,我们可以发现crp与wbc之间存在负相关(r=-0.413)，这意味着当wbc处于低水平值时，crp缺失值更容易出现。该命令有点复杂，“names(data)%in%c(“sex”)”代码返回一个逻辑向量，该逻辑向量对于与“sex”匹配的names(data)中的每个变量名行结果为TRUE，反之为FALSE。“!”符号可反转逻辑向量的值。然而，相关性的分析并不能取代采用专业知识来判断缺失值是否为NMAR。换而言之,运用学科知识来判断对于排除NMAR也是至关重要。

6. 总结

缺失值在临床大数据研究中无处不在，而有时缺失模式隐含的机制可能是复杂的。在这类情况下，一些高级技术在处理缺失值上是有帮助的。缺失数据按照机制来分可以分为以下三种：MCAR、MAR和NMAR。前两种类型允许使用其他协变量来估算插补值，而最后一种类型则需要外在的专业知识来判断。可运用表格的形式来探讨数据缺失模式。此外，VIM程序包具有多个函数，可用于缺失数据的图形显示。缺失值与其他变量之间的关系对于缺失值隐含的机制提供了进一步的线索，该机制可采用相关性分析来探索。

References:

van Buuren, S., Boshuizen, H.C., Knook, D.L. (1999) Multiple imputation of missing blood pressure covariates in survival analysis. Statistics in Medicine, 18, 681–694.
van Buuren, S., Brand, J.P.L., Groothuis-Oudshoorn C.G.M., Rubin, D.B. (2006) Fully conditional specification in multivariate imputation. Journal of Statistical Computation and Simulation, 76, 12, 1049–1064.
van Buuren, S., Groothuis-Oudshoorn, K. (2011). mice: Multivariate Imputation by Chained Equations in R. Journal of Statistical Software, 45(3), 1–67. doi: 10.18637/jss.v045.i03
Van Buuren, S. (2018). Flexible Imputation of Missing Data. Second Edition. Chapman & Hall/CRC. Boca Raton, FL.
Van Buuren, S., Groothuis-Oudshoorn, K. (2011). mice: Multivariate Imputation by Chained Equations in R. Journal of Statistical Software, 45(3), 1-67. doi: 10.18637/jss.v045.i03
Van Buuren, S. (2018). Flexible Imputation of Missing Data. Second Edition. Chapman & Hall/CRC. Boca Raton, FL.
Van Buuren, S., Brand, J.P.L., Groothuis-Oudshoorn C.G.M., Rubin, D.B. (2006) Fully conditional specification in multivariate imputation. Journal of Statistical Computation and Simulation, 76, 12, 1049–1064.
Van Buuren, S. (2007) Multiple imputation of discrete and continuous data by fully conditional specification. Statistical Methods in Medical Research, 16, 3, 219–242.
Van Buuren, S., Boshuizen, H.C., Knook, D.L. (1999) Multiple imputation of missing blood pressure covariates in survival analysis. Statistics in Medicine, 18, 681–694.
Brand, J.P.L. (1999) Development, implementation and evaluation of multiple imputation strategies for the statistical analysis of incomplete data sets. Dissertation. Rotterdam: Erasmus University.
参考《白话统计》
参考《聪明统计学》