处理数据时，难免会遇到缺失值。如何处理缺失值呢？统计学中有个概念叫imputation，看到一些资料将其翻译成“填补”或“插补”，但我觉得翻译成“设算”似乎更合适些。
本文采用R自带的iris数据进行演示。首先修改iris数据集，将其中一些值更改为NA。
注意：不能直接修改iris数据集，要先将其赋予一个变量。
# Changes some values to NA ----
data <- iris
data[1,3] <- NA
data[c(1,4),c(1,2)] <- NA
data[c(6,9),2] <- NA
data[10,1] <- NA
data[c(8:9),5] <- NA
head(data,10)
> head(data,10)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 NA NA NA 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 NA NA 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 NA 1.7 0.4 setosa
7 4.6 3.4 1.4 0.3 setosa
8 5.0 3.4 1.5 0.2 <NA>
9 4.4 NA 1.4 0.2 <NA>
10 NA 3.1 1.5 0.1 setosa
1. 缺失值可视化
使用visdat对缺失值进行可视化。
# Visualizing missing values ----
install.packages("visdat")
install.packages("dplyr")
library(visdat)
library(dplyr)
data %>% vis_miss()
图1.缺失值可视化
2. 评估含缺失值的变量之间的关系
使用naniar检查含有缺失值的变量之间是否存在关系。
install.packages("naniar")
library(naniar)
## Is it often that we have both Sepal.Length and Sepal.Width missing at the same time?
data %>% gg_miss_upset()
图2.缺失值的数量统计
3. 设算填补
使用simputation对缺失值进行填补。常用的方法有：线性模型、平均值、均值等。
install.packages("simputation")
library(simputation)
3.1 线性模型填补（linear model imputation）
使用线性模型对Sepal.Length的缺失值进行填补。注意到第1行和第3行依然存在缺失值，因为它们的Sepal.Width（在该模型中作为predictor variable）也缺失，因此该线性模型无法基于Sepal.Width对Sepal.Length作出预测。但它仍可以返回部分填补后的结果。
# Impute Sepal.Length using a linear model use the `impute_lm` function
data_lm <- impute_lm(data, Sepal.Length~Sepal.Width + Species)
> head(data_lm, 10)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 NA NA NA 0.2 setosa
2 4.900000 3.0 1.4 0.2 setosa
3 4.700000 3.2 1.3 0.2 setosa
4 NA NA 1.5 0.2 setosa
5 5.000000 3.6 1.4 0.2 setosa
6 5.400000 NA 1.7 0.4 setosa
7 4.600000 3.4 1.4 0.3 setosa
8 5.000000 3.4 1.5 0.2 <NA>
9 4.400000 NA 1.4 0.2 <NA>
10 4.746835 3.1 1.5 0.1 setosa
3.2 中位数填补（median imputation）
按照Species分别计算中位数，并进行填补。
# Impute Sepal.Length using the group median
data_median <- impute_median(data_lm, Sepal.Length ~ Species)
> head(data_median, 10)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.000000 NA NA 0.2 setosa
2 4.900000 3.0 1.4 0.2 setosa
3 4.700000 3.2 1.3 0.2 setosa
4 5.000000 NA 1.5 0.2 setosa
5 5.000000 3.6 1.4 0.2 setosa
6 5.400000 NA 1.7 0.4 setosa
7 4.600000 3.4 1.4 0.3 setosa
8 5.000000 3.4 1.5 0.2 <NA>
9 4.400000 NA 1.4 0.2 <NA>
10 4.746835 3.1 1.5 0.1 setosa
3.3 决定树模型填补（decision tree model imputation）
这里，我们采用决定树模型对Species的缺失值进行填补。除了Species外的所有变量都被用作预测值（predictor）。
data_cart <- impute_cart(data_median, Species ~ .)
> head(data_cart, 10)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.000000 NA NA 0.2 setosa
2 4.900000 3.0 1.4 0.2 setosa
3 4.700000 3.2 1.3 0.2 setosa
4 5.000000 NA 1.5 0.2 setosa
5 5.000000 3.6 1.4 0.2 setosa
6 5.400000 NA 1.7 0.4 setosa
7 4.600000 3.4 1.4 0.3 setosa
8 5.000000 3.4 1.5 0.2 setosa
9 4.400000 NA 1.4 0.2 setosa
10 4.746835 3.1 1.5 0.1 setosa
3.4 链式填补（chaining imputation）
可以使用|>运算符将以上操作打包。
# Chaining imputation methods ----
data_chain <- data |>
impute_lm(Sepal.Length ~ Sepal.Width + Species) |>
impute_median(Sepal.Length ~ Species) |>
impute_cart(Species ~ .)
3.5 多元变量模型预测填补（imputation model for multiple variables）
我们可以根据不同物种的Petal.Length分别对它们的Sepal.Length和Sepal.Width进行预测。
data_rlm <- impute_rlm(data, Sepal.Length + Sepal.Width ~ Petal.Length + Species)
> head(data_rlm, 10)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 NA NA NA 0.2 setosa
2 4.900000 3.000000 1.4 0.2 setosa
3 4.700000 3.200000 1.3 0.2 setosa
4 5.042242 3.427884 1.5 0.2 setosa
5 5.000000 3.600000 1.4 0.2 setosa
6 5.400000 3.489194 1.7 0.4 setosa
7 4.600000 3.400000 1.4 0.3 setosa
8 5.000000 3.400000 1.5 0.2 <NA>
9 4.400000 NA 1.4 0.2 <NA>
10 5.042242 3.100000 1.5 0.1 setosa
3.6 分组数据的填补（grouping data for imputation）
在formula中使用|来设定分组。
注意：这里的分组数据不能含有缺失值，因此我们采用data_cart。
data_group <- impute_lm(data_cart, Sepal.Length ~ Petal.Width | Species)
该代码等同于：
data_group |> group_by(Species) |>
impute_lm(Sepal.Length ~ Petal.Width)
注意：使用group_by前要先加载dplyr包。
> head(data_group, 10)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.000000 NA NA 0.2 setosa
2 4.900000 3.0 1.4 0.2 setosa
3 4.700000 3.2 1.3 0.2 setosa
4 5.000000 NA 1.5 0.2 setosa
5 5.000000 3.6 1.4 0.2 setosa
6 5.400000 NA 1.7 0.4 setosa
7 4.600000 3.4 1.4 0.3 setosa
8 5.000000 3.4 1.5 0.2 setosa
9 4.400000 NA 1.4 0.2 setosa
10 4.746835 3.1 1.5 0.1 setosa
3.7 自定义填补方法（specify your own method）
simputation提供impute_proxy功能，可以自定义填补方法。
data_own <- impute_proxy(data, Sepal.Length ~ median(Sepal.Length,na.rm=TRUE)/median(Sepal.Width, na.rm=TRUE) * Sepal.Width | Species)
> head(data_own, 10)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 NA NA NA 0.2 setosa
2 4.900000 3.0 1.4 0.2 setosa
3 4.700000 3.2 1.3 0.2 setosa
4 NA NA 1.5 0.2 setosa
5 5.000000 3.6 1.4 0.2 setosa
6 5.400000 NA 1.7 0.4 setosa
7 4.600000 3.4 1.4 0.3 setosa
8 5.000000 3.4 1.5 0.2 <NA>
9 4.400000 NA 1.4 0.2 <NA>
10 4.558824 3.1 1.5 0.1 setosa
3.8 基于其它数据训练得到的模型对缺失值进行预测（inputing a dataset with models trained on another dataset）
比如我们已经有完整的iris数据，我们可以基于这个完整的数据建模：
model <- lm(Sepal.Length ~ Sepal.Width + Species, data=iris)
利用该模型预测缺失值：
data_model <- impute(data, Sepal.Length ~ model)
> head(data_model, 10)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 NA NA NA 0.2 setosa
2 4.900000 3.0 1.4 0.2 setosa
3 4.700000 3.2 1.3 0.2 setosa
4 NA NA 1.5 0.2 setosa
5 5.000000 3.6 1.4 0.2 setosa
6 5.400000 NA 1.7 0.4 setosa
7 4.600000 3.4 1.4 0.3 setosa
8 5.000000 3.4 1.5 0.2 <NA>
9 4.400000 NA 1.4 0.2 <NA>
10 4.742432 3.1 1.5 0.1 setosa
4. 参考资料
[1] https://www.business-science.io/code-tools/2021/01/19/handle-missing-data.html
[2] https://cran.r-project.org/web/packages/simputation/vignettes/intro.html
5. 软件及R包信息
> sessionInfo()R version 4.2.2 (2022-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C
[3] LC_TIME=C.UTF-8 LC_COLLATE=C.UTF-8
[5] LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[1] simputation_0.2.8 naniar_0.6.1 dplyr_1.0.10
[4] visdat_0.5.3
loaded via a namespace (and not attached):
[1] fansi_1.0.3 utf8_1.2.2 MASS_7.3-58.1
[4] grid_4.2.2 R6_2.5.1 gtable_0.3.1
[7] lifecycle_1.0.3 magrittr_2.0.3 scales_1.2.1
[10] ggplot2_3.4.0 pillar_1.8.1 rlang_1.0.6
[13] cli_3.4.1 rpart_4.1.19 vctrs_0.5.0
[16] generics_0.1.3 tools_4.2.2 gower_1.0.0
[19] glue_1.6.2 munsell_0.5.0 compiler_4.2.2
[22] colorspace_2.0-3 pkgconfig_2.0.3 tidyselect_1.2.0
[25] tibble_3.1.8
