拨开数据迷云雾，得见缺失值真面目

碎碎念专栏之缺失值

大家好，我是风。这里是风风的碎碎念专栏。这个专栏主要用来分享一下杂七杂八不好归类，但是又非常普遍需要我们掌握的内容，所以将其命名为“碎碎念”专栏。在我们开启单细胞数据分析专栏之前，大概会先更新3-4期碎碎念专栏，跟大家介绍一些后面可能会用到的小知识，本期主要内容是“缺失值的识别与处理”。

什么是缺失值？

在上课时候，也许你会发现，老师们用来讲解和演示的数据大部分是完整的数据，也就是没有缺失值，这样的数据适合讲解，也能够直接应用各种统计方法，比如Cox回归，因为大部分统计方法都是假定处理的是完整的矩阵或者数据框。但是我们不得不面对一个问题，就是当你拿到一份自己的数据的时候，往往数据存在很多缺失值。

以我们最熟悉的临床数据为例：患者可能忘记回答了某个问题、没做某个检查或者拒绝回答某个敏感问题，也可能是问卷太长没有耐心，或者是数据记录错误等等，这些原因都可能导致缺失值，并且有些缺失值是有意的，有些是由于一些不可控因素。如果我们不对缺失值进行合理的处理，则可能会丢失某些重要信息，或者导致研究结果在一定程度上无效。

缺失值的分类

在统计学中，一般把缺失数据分为三类：

完全随机缺失数据（MCAR）：即某变量的缺失数据与其他任何观测或未观测的数据都不相关，则可以将其视为完全随机缺失。需要注意的是，如果所有有缺失值的变量都是完全随机缺失，那可以直接将完整的数据看作是一个更大数据集的简单随机抽样。
随机缺失数据（MAR）：若某个变量的缺失数据和其他观测相关，与它自己的观测值不相关，则数据为随机缺失数据，比如婴儿的睡眠质量程度的缺失值更多，这可能跟父母的观察程度有关，跟婴儿本身的睡眠质量无关，这种数据就是随机缺失数据。
非随机缺失（NMAR）：若缺失数据不属于MCAR和MAR两种，则数据为非随机缺失数据。例如既往是否有传染病史，这样的缺失数据可能是多种原因导致，即为非缺失数据。

其实缺失数据属于哪种类型并不是特别重要，因为大部分处理的方法都是假定缺失数据为MCAR或者MAR，这种情况下在处理完缺失值后可以直接对感兴趣的关系进行分析。但是，如果数据是非随机缺失数据，此时需要对感兴趣的关系进行建模，也需要对缺失值的生成机制进行建模。关于NMAR的处理比较复杂，有模型混合法和模型选择法，我们本期不进行讨论，只讨论关于MCAR和MAR类型缺失数据的处理，后文所说的缺失值也仅代表MCAR和MAR两种。

缺失值的处理方法

缺失值的处理可以分为两个步骤：1. 识别缺失值；2. 处理缺失值

。

识别缺失值在R语言中有几个简单的方法，分别是：

is.na()函数，!complete.cases()函数，VIM包

，前面两个都比较简单，VIM包的用法在后文介绍。

处理缺失值的方法也很多，在这里简单给大家列举一下：

直接删除法：对于无效的样本，或者有缺失值的样本不多的情况下，可以直接删除具有缺失值的样本，使用omit.na()，那有缺失值的样本要占有总样本的多少才算“不多”呢？看你数据的具体情况，如果数据只有30多个样本，那么每一个样本都很珍贵，需要谨慎处理；如果总数据有3000多个样本，数据大样本的情况下，具有缺失值的样本不超过总样本的20%-30%，都可以进行直接删除，将剩下的完整数据视为总数居的简单随机抽样子集处理。
最大似然估计法：这个方法可以使用mvmle包对缺失数据进行填补。
简单插补法：也就是单个插补，可以使用Hmisc包进行填补。
多重插补法：同样是插补缺失值的常用方法，可以使用的R包很多，经典的R包有：mice包、mi包、mitools包、amelia包等等。缺失值的处理是一个非常复杂的问题，完整的做法还要考虑到实际数据的用途，以及数据是规则数据还是不规则数据等等内容，想要全部介绍，那至少得再开一门精品课，我们这里这对最经典的几种方法进行简单介绍。

缺失值的处理

按照我们前面的讲解，处理缺失值的第一步是识别缺失值，所以我们先来看看任何识别缺失值。在R语言中，缺失值有多种表达形式，所包含的意义也不一样，NaN代表不可能值，NA代表缺失值，Inf代表正无穷大，-Inf代表负无穷大，is.na()用来识别缺失值，is.nan()用来识别不可能值，is.infinite()用来识别无穷值，三个函数都会返回逻辑值TRUE or FALSE，TRUE代表有，FALSE代表没有。在之前的推文中有留言说想看看可变剪切数据的内容，从安德森肿瘤中心的网站中下载到的可变剪切数据，为了保证数据量，都会含有NA值，生信全书下篇第三段位我给大家展示的方法是直接用na.omit()，这里我们同样用一份可变剪切数据来讲解，看看其他处理缺失值的方法可以怎么操作。

识别缺失值

rm(list=ls()) #清除环境Sys.setenv(LANGUAGE = "en") #显示英文报错信息options(stringsAsFactors = FALSE) #禁止chr转成factorlibrary(tidyverse)Splice_ME <- read_tsv("ME.txt") %>% column_to_rownames("symbol") # 使用可变剪切数据进行演示Splice_ME[1:3,1:3] # 行名为样本名，列名为可变剪切事件## CAMTA1-506-ME MTFR1L-1211-ME COL16A1-1492-ME## TCGA_3B_A9HI 0.0165 0.7805 1## TCGA_3B_A9HJ 0.0048 0.4940 1## TCGA_3B_A9HL 0.0000 0.7234 1# 识别缺失值# is.na(Splice_ME) # 识别NA# complete.cases(Splice_ME) # 识别数据是否为完整数据# Splice_ME[complete.cases(Splice_ME),] # 识别没有缺失值的行# !complete.cases(Splice_ME) # 识别数据是否为不完整数据# Splice_ME[!complete.cases(Splice_ME),] # 识别有缺失值的行# 识别缺失值的有用信息sum(is.na(Splice_ME$`ABI1-11046-ME`)) ## [1] 23# [1] 23 表明ABI1-11046-ME事件有23个缺失值mean(is.na(Splice_ME$`ABI1-11046-ME`)) ## [1] 0.08812261# [1] 0.08812261 表明8.8%的样本在ABI1-11046-ME事件有缺失值mean(!complete.cases(Splice_ME)) ## [1] 0.954023# [1] 0.954023 表明Splice_ME数据中，95.4%的样本至少有一个缺失值

# 需要注意的是，complete.cases()函数只会认为NA和NaN是缺失值，对于无穷值Inf和-Inf都会被认为是有效值，那有没更直观的方法来识别缺失值呢？比如把识别结果变成一个另一个数据集？我们来试试：

# install.packages("mice") # 安装mice包library(mice)NA_table <- as.data.frame(md.pattern(Splice_ME))

NA_table[1:3,1:3]## CAMTA1-506-ME ALKBH3-15465-ME CDK4-22735-ME## X12 1 1 1## X1 1 1 1## X2 1 1 1

# 使用md.pattern识别数据中的缺失值，然后生成表格，0是缺失值，1是非缺失值，然后赋值为NA_table，这样我们就得到了一个0-1矩阵，可以直观知道有哪些缺失值了

# 如果想用图形来展示呢？也可以，用我们前面提到的VIM包# install.packages("VIM")library(VIM)# 使用aggr()函数绘制每个变量的缺失值和计算每个变量组合的缺失值数目aggr(Splice_ME, prop = TRUE, numbers = TRUE)

# matrixplot()可以将数值型数据重新转换到[0,1]之间，并且用灰度表示大小，红色为缺失值，浅色为数值较小，深色为数值较大

matrixplot(Splice_ME)

# 生成一个交互矩阵，鼠标可以点击相应的样本重新排列矩阵。# 如果想展示特定两个变量之间的关系，也可以使用matrixplotmatrixplot(Splice_ME[c("status","time")], pch = c(20), col = c("darkgray","red","skyblue")) # 红色代表缺失值

前面说过，是MCAR还是MAR主要要看跟其他观测之间的相关性，那我们可不可以计算各个观测之间的相关性呢？当然可以，但是注意，由于原始矩阵含有缺失值，所以我们使用映射缺失值的 0-1 矩阵计算相关性，这种相关性结果也可以认为是原始观测数据之间的相关性，来看看怎么操作：

Splice_ME1 <- is.na(Splice_ME)x <- as.data.frame(abs(Splice_ME1))# head(Splice_ME, n = 4); head(x, n= 4) # 观察两份数据的差异,发现x是Splice_ME转换后的0-1矩阵# 在x中，1代表缺失值NA，0代表非缺失值y <- x[which(apply(x, 2, sum) >= 0)] # 提取含有缺失值或者全部是缺失值的变量# 计算变量间的相关性z1 <- cor(y) # 计算所有观测之间的相关性矩阵z2 <- cor(Splice_ME1,y,use = "pairwise.complete.obs") # 计算缺失值和其他观测之间的相关性矩阵# head(as.data.frame(z1), n = 4); head(as.data.frame(z2), n= 4)# 我们进行简单可视化看一下图片library(corrplot)corrplot(z1,tl.cex = 0.1);corrplot(z2,tl.cex = 0.1)

缺失值的插补

前面我们已经对缺失值进行了简单的探讨，具体的内容等以后有机会我们在训练营再一一展开，接下来我们看看如果对数据进行插补：

行删除法

# 我们先画个热图看一下还没填补的数据总体是什么样的情况library(pheatmap)pheatmap(Splice_ME, show_rownames = FALSE, show_colnames = FALSE)

# 发现缺失值影响了我们对数据的观测，因此我们需要对缺失值进行处理# 首先是行删除法，假设我们要删掉含有NA的可变剪切事件Splice_ME3 <- as.data.frame(t(Splice_ME)) # 转置行列Splice_ME4 <- Splice_ME3[complete.cases(Splice_ME3),] # 提取没有缺失值的行Splice_ME5 <- as.data.frame(t(na.omit(Splice_ME3))) # na.omit() 也可以直接得到相同结果pheatmap(log(Splice_ME4+0.01), show_rownames = FALSE, show_colnames = FALSE)

多重插补法

多重插补法是一种基于重复模拟的处理缺失值的方法，也是复杂矩阵最常用的一种方法，主要是利用Gibbs的原理，也就是利用其他变量的数据来预测缺失值的有效值，主要的方法是通过不断迭代直到所有缺失值都收敛，我们看看一般怎么操作：

library(mice) # 使用mice包进行多重插补imp <- mice(Splice_ME3, seed = 20210109)fit <- with(imp, glm(TCGA_3R_A8YX~.,data = Splice_ME3)) # 使用广义相加模型进行预测pooled <- pool(fit) # 得到相应的列表结果summary(pooled)imp$imp$TCGA_3B_A9HI# 展示具体变量的多次插补imp$predictorMatrix # 行表示插补变量，列表示为插补提供信息的变量，1和0代表使用和未使用## TCGA_3B_A9HI TCGA_3B_A9HJ TCGA_3B_A9HL TCGA_3B_A9HO## TCGA_3B_A9HI 0 0 0 0## TCGA_3B_A9HJ 0 0 0 0## TCGA_3B_A9HL 0 0 0 0## TCGA_3B_A9HO 0 0 0 0## TCGA_3B_A9HP 0 0 0 0## TCGA_3B_A9HQ 0 0 0 0

## TCGA_3B_A9HR                 0            0            0            0

Splice_ME6 <- complete(imp, action = 1) # 提取完整数据集：通过action的值展示m个完整数据集中的其中一个

# head(Splice_ME6, n= 10) pheatmap(log(Splice_ME6+0.01), show_rownames = FALSE, show_colnames = FALSE)

这样我们缺失值处理的内容就完成啦，在本章节，我们学习了什么是缺失值，以及如何识别缺失值，还有怎么处理缺失值，需要主要的是，缺失值的处理一定一定一定要结合实际情况来处理，我们这里也只展示了行删除和多重插补两种方法，实际上缺失值处理有n多种方法可以选择，如何结合数据选择合适的方法才是一门大学问！好了，今天就到这里，如果你有具体的问题，欢迎在留言区打出来一起探讨。

后台回复“feng33”获得本次代码和数据，那我们下次再见吧！*^_^*