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

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

碎碎念专栏之缺失值

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

什么是缺失值?

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


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

缺失值的分类

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

  1. 完全随机缺失数据(MCAR:即某变量的缺失数据与其他任何观测或未观测的数据都不相关,则可以将其视为完全随机缺失。需要注意的是,如果所有有缺失值的变量都是完全随机缺失,那可以直接将完整的数据看作是一个更大数据集的简单随机抽样。
  2. 随机缺失数据(MAR:若某个变量的缺失数据和其他观测相关,与它自己的观测值不相关,则数据为随机缺失数据,比如婴儿的睡眠质量程度的缺失值更多,这可能跟父母的观察程度有关,跟婴儿本身的睡眠质量无关,这种数据就是随机缺失数据。
  3. 非随机缺失(NMAR:若缺失数据不属于MCAR和MAR两种,则数据为非随机缺失数据。例如既往是否有传染病史,这样的缺失数据可能是多种原因导致,即为非缺失数据。
其实缺失数据属于哪种类型并不是特别重要,因为大部分处理的方法都是假定缺失数据为MCAR或者MAR,这种情况下在处理完缺失值后可以直接对感兴趣的关系进行分析。但是,如果数据是非随机缺失数据,此时需要对感兴趣的关系进行建模,也需要对缺失值的生成机制进行建模。关于NMAR的处理比较复杂,有模型混合法和模型选择法,我们本期不进行讨论,只讨论关于MCAR和MAR类型缺失数据的处理,后文所说的缺失值也仅代表MCAR和MAR两种。

缺失值的处理方法

缺失值的处理可以分为两个步骤:1. 识别缺失值;2. 处理缺失值
识别缺失值在R语言中有几个简单的方法,分别是:is.na()函数,!complete.cases()函数,VIM包,前面两个都比较简单,VIM包的用法在后文介绍。
处理缺失值的方法也很多,在这里简单给大家列举一下:

  1. 直接删除法:对于无效的样本,或者有缺失值的样本不多的情况下,可以直接删除具有缺失值的样本,使用omit.na(),那有缺失值的样本要占有总样本的多少才算“不多”呢?看你数据的具体情况,如果数据只有30多个样本,那么每一个样本都很珍贵,需要谨慎处理;如果总数据有3000多个样本,数据大样本的情况下,具有缺失值的样本不超过总样本的20%-30%,都可以进行直接删除,将剩下的完整数据视为总数居的简单随机抽样子集处理。
  2. 最大似然估计法:这个方法可以使用mvmle包对缺失数据进行填补。
  3. 简单插补法:也就是单个插补,可以使用Hmisc包进行填补。
  4. 多重插补法:同样是插补缺失值的常用方法,可以使用的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)
拨开数据迷云雾,得见缺失值真面目
拨开数据迷云雾,得见缺失值真面目

缺失值的插补

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

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)
拨开数据迷云雾,得见缺失值真面目


2

多重插补法

多重插补法是一种基于重复模拟的处理缺失值的方法,也是复杂矩阵最常用的一种方法,主要是利用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            0Splice_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”获得本次代码和数据,那我们下次再见吧!*^_^*

风之美图系列传送门

手把手教你复现一篇CNS级别美图!附代码,建议收藏!

CNS级别的美图复现起来难不难?小白也能快速上手(附代码)(二)

小白课堂 | 学个R包,复现一张CNS级别美图(附代码)

学一个R包给文章加一张CNS级别的美图(四)

学一个R包给文章加一张CNS级别的美图(五)

学一个R包给文章加一张CNS级别的美图(六)

拆和叠?原来你和高阶论文的图片就差这点

可繁可简,可盐可甜,这才是upset的真面目吧!

高分文章不讲xx,竟然用这种图片

荡气回肠一柄剑,看T细胞舍生取义

Cancer Cell的一张美图送给你,新年也发个CC

嘀嘀嘀,美图系列末班车已到站,来收拾最后的图片下车

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

END
撰文丨风   风
排版丨四金兄
值班 | 菠小萝
主编丨小雪球

欢迎大家关注解螺旋生信频道-挑圈联靠公号~


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

2021年国自然申请指南已公布

2021申请国自然还很迷茫?
国自然申请指南找不到重点?
酸谈为你解读最新国自然申请指南!

本周四酸谈老谈直播主题
---2021年国自然申请指南解读
---

下方视频号点击预约直播
  酸菜老谈助你冲刺国自然!

👇👇👇

周四酸谈直播基本信息

主题:《2021年国自然申请指南解读》
日期:1月21日 18:00—20:00
(👆点击“预约”就完事了!!!)

生物医学科研方法

医院里,最忌讳的就是埋头苦干的勤奋

2021-1-18 6:35:42

生物医学科研方法

中科院研究生被杀案,凶手已被执行死刑

2021-1-18 6:36:31