确定你会统计？大老粗别走，教你如何识别「离群值」和处理「缺失值」！

责编：科研猫 | 馋猫

无论是前瞻性数据收集还是回顾性数据收集，数据集中通常都会出现离群值或缺失值。对于统计学家来说，离群值和缺失值通常是一个棘手的问题，如果处理不当可能会导致错误。离群值可能会导致我们的结果偏离真实结果，而缺失值造成的信息损失可能会导致建模失败。因此，在执行数据分析之前，正确识别离群值并处理缺失值非常重要。本推文讨论的内容应该在建模之前执行。虽然本推文在整个统计模型系列中较为置后，却至关重要，望警醒。

01

离群值的识别

什么是离群值？简而言之就是，超越人类常识和不符合逻辑的变量的值即是离群值。例如，我们从一组患者中采集了空腹血糖，其中一名患者的空腹血糖超过50 mmol / L，这显然是一个异常值。再例如，我们调查了上海市徐汇区60岁以上老年人的高血压患病率。如果受试者的SBP超过1400 mmHg，则显然是异常值。可能是记录错误，实际SBP较可能是140.0 mmHg。

有时离群值是一个相对的概念，与我们的临床研究数据的收集环境有关。例如，如果我们的研究对象是10岁以下的孩子，那么这个年龄段的孩子不太可能是研究生，并且他们的身高不太可能超过170厘米，体重不太可能超过100公斤。当我们抽取的样本不好时，也可能存在产生异常值的情况。例如，从A区域中抽取了1,000个人，从B区域中抽取了100个人。如果该集合的值异常高于或异常低于区域A的值，B区域中的100个人很有可能是个孤独的集合。

当我们研究一项干预措施的效果时，如果只有部分患者有显著效果，这部分数据与其他疗效不太明显的患者相比是“离群值”，但这些异常值正是我们最关心的。因此，对于异常值的判断，要联系实际，不要武断，以免出现严重错误。当我们对数据不确定时，最好的解决方案是检查原始数据记录。

下面我将介绍几个常用的函数来识别数据集中的异常值。假设我们收集了1000个受试者的身高。首先，我们可以使用boxplot()函数绘制一个箱状图来描述数据。接下来使用range()函数帮助我们找到这些变量的最大值和最小值。

首先，我们模拟了1000名身高100-250厘米的受试者。使用range()查看这组患者的收缩压范围。

1set.seed(123) 
2height <- sample(100:250,1000,replace = TRUE) 
3boxplot(height)
4range(height)

使用min()和max()函数返回对象的最小值和最大值。

1min(height)
2max(height)

当处理含有缺失值的数据时，要设置参数na.rm = TRUE。

1height <- c(sample(100:250,998,replace = TRUE),NA,NA) 
2boxplot(height)
3range(height,na.rm = TRUE)
4mean(height,na.rm = TRUE)
5median(height,na.rm = TRUE)
6sd(height,na.rm = TRUE)
7quantile(height,na.rm = TRUE)
8fivenum(height)

上述方法可以帮助我们识别最大值或最小值，但有时极限值并不是单独出现的，而是在聚类中，因此上述方法识别异常值是不够的。在实际的研究背景下，我们通常根据变量的均值和标准差，或中位数和四分位数（Tukey方法）来定义数据的异常值。例如，我们可以设置大于或小于mean±3sd均为异常值。当然，我们也可以对分类变量的某个值进行异常判断。例如，性别值为1=男性，2=女性。如果赋值为3，则为异常值。这里我们介绍一个自定义函数。该函数根据四分位Tukey方法判断异常值，有效地避免了极限值对均值和标准差的影响。函数如下：

 1outlierKD <- function(dt, var) {
 2  var_name <- eval(substitute(var),eval(dt))
 3  tot <- sum(!is.na(var_name))
 4  na1 <- sum(is.na(var_name))
 5  m1 <- mean(var_name, na.rm = T)
 6  par(mfrow=c(2, 2), oma=c(0,0,3,0))
 7  boxplot(var_name, main="With outliers")
 8  hist(var_name, main="With outliers", xlab=NA, ylab=NA)
 9  outlier <- boxplot.stats(var_name)$out
10  mo <- mean(outlier)
11  var_name <- ifelse(var_name %in% outlier, NA, var_name)
12  boxplot(var_name, main="Without outliers")
13  hist(var_name, main="Without outliers", xlab=NA, ylab=NA)
14  title("Outlier Check", outer=TRUE)
15  na2 <- sum(is.na(var_name))
16  cat("Outliers identified:", na2 - na1, "\n")
17  cat("Propotion (%) of outliers:", round((na2 - na1) / tot*100, 1), "\n")
18  cat("Mean of the outliers:", round(mo, 2), "\n")
19  m2 <- mean(var_name, na.rm = T)
20  cat("Mean without removing outliers:", round(m1, 2), "\n")
21  cat("Mean if we remove outliers:", round(m2, 2), "\n")
22  response <- readline(prompt="Do you want to remove outliers and to replace with NA? [yes/no]: ")
23  if(response == "y" | response == "yes"){
24    dt[as.character(substitute(var))] <- invisible(var_name)
25    assign(as.character(as.list(match.call())$dt), dt, envir = .GlobalEnv)
26    cat("Outliers successfully removed", "\n")
27    return(invisible(dt))
28  } else{
29    cat("Nothing changed", "\n")
30    return(invisible(var_name))
31  }
32}

自定义函数只有两个参数，第一个参数是数据集的名称，第二个参数是变量名；只要正确替换数据集和变量名，读取就可以直接运行代码。这里我们是以箱形图的外值为离群值，我们还可以根据专业知识重新设置离群值的定义，比如大于或小于mean±3sd。在函数结束时，还将设置用户输入的代码。用户可以通过键入“yes”或“no”来确定是否消除数据集中函数识别的异常值。

下面我们模拟一组数据来验证这个自定义异常值识别函数的功能。

1set.seed(123) 
2df <- data.frame(height = c(sample(100:250, 1000, replace = TRUE), NA, 380, 20)) 
3outlierKD(df, height)

02

缺失值的处理

数据丢失在临床研究中很常见。例如，护士在收集数据时，可能会因为工作繁忙而忘记记录某个时间点的尿量；当研究人员想研究乳酸变化对死亡率的影响时，患者可能只监测某个时间点的血乳酸值。缺乏数据的其他原因还包括编码错误、设备故障和调查研究中的应答者没有应答等。在统计软件包中，一些函数（如Logistic回归）可能会自动删除丢失的数据。如果只有少量的不完全观测，那么这种处理就不会有太大问题。

但是，当存在大量包含缺失值的观测值时，这些函数中的默认行删除可能会导致大量信息丢失。在这种情况下，分析人员应该仔细研究数据丢失可能导致的机制，并找到适当的处理方法。

如何处理缺失值是临床统计学家头疼的问题，所以我们也应该予以重视。数据的缺失或缺失程度直接影响到数据的质量，而数据的质量最终影响到我们的研究成果。如果对缺失数据的处理不当，很可能导致整个统计分析失败。本推文介绍了在R中如何处理丢失的数据，并介绍了处理丢失数据的一些基本技巧。

在R中，“NA”表示为一个缺失的值。当将带有空单元格的Excel表导入R控制台时，这些空单元格将被NA替换。这与STATA用“.”替换“空单元格”不同。R中的数值变量和字符变量使用相同的缺失值符号。R提供一些函数来处理缺失值。要确定向量是否包含缺少的值，可以使用is.na（）函数。“is.na（）”函数是用于确定元素是否为na类型的最常用方法。它返回与传入参数长度相同的对象，并且所有数据都是逻辑值（FALSE或TRUE）。假设我们有6个病人，但是只记录了4个值，而缺少了2个。

1x <- c(1.8,2.3,NA,4.1,NA,5.7) 
2is.na(x) 

03

缺失值的可视化

缺失值的可视化可以帮助我们更直观地观察数据集中的缺失值，这将有助于我们以后对缺失值进行插值。在本推文中，笔者将主要向读者介绍VIM包的使用。以下的演示数据集是R语言的内置数据集"airquality"。

1data(“airquality”) 
2str(airquality)
3## ‘data.frame’:    153 obs. of  6 variables: 
4##  $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ... 
5##  $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ... 
6##  $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ... 
7##  $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ... 
8##  $ Month  : int  5 5 5 5 5 5 5 5 5 5 ... 
9##  $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

"airquality"数据集包含了153个观测值和6个变量。从以上结果中，我们可以看到该数据集中有缺失值。在可视化之前，首先使用mice包中的md.pattern()函数探索缺失的数据模式。

1install.packages(“mice”) 
2library(mice)
3md.pattern(airquality)
4##     Wind Temp Month Day Solar.R Ozone    
5## 111    1    1     1   1       1     1  0 
6## 35     1    1     1   1       1     0  1 
7## 5      1    1     1   1       0     1  1 
8## 2      1    1     1   1       0     0  2 
9##        0    0     0   0       7    37 44

在输出表格中，“1”表示非缺失值，“0”表示缺失值。第一列显示了唯一缺失数据模式的数目。在我们的例子中，111个观测值没有缺失数据，35个观测值仅在Ozone变量中有缺失数据，5个观测值仅在Solar. R变量中有缺失数据。最右边的一列显示了特定缺失模式中缺失变量的数目。例如，如果第一行中没有缺失值，则显示为“0”。最后一行计算每个变量缺失值的数量。例如，“Wind”变量没有缺失值，显示“0”，而Ozone变量有37个缺失值。在研究中，一些含有更多缺失值的变量可能会被剔除。显然，表格可以提供有用的参考信息。

下面我们调用VIM包来实现缺失值的可视化。研究缺失数据模式对于选择合适的插值方法来估计缺失值是必要的。因此，需要在插值操作之前执行可视化工具，并且通常应该在缺失数据插值之后进行诊断，以确定插值是否合理。可用于可视化丢失数据的函数如下:aggr()、matrixplot()、scattMiss()和marginplot()，以下我们主要演示了aggr()和marginplot()两个函数。

 1install.packages(‘VIM’) 
 2library(VIM)
 3aggr_plot <- aggr(airquality , 
 4  col=c(‘red’,’blue’),
 5  numbers=TRUE, 
 6  sortVars=TRUE, 
 7  labels=names(airquality),
 8  cex.axis=.7, 
 9  gap=3, 
10  ylab=c(“Histogram of missing data”,”Pattern”))

aggr()函数的作用是帮助我们可视化缺失的值。左图是缺失值比例直方图。从下图中可以看出Ozone和Solar. R有缺失值，其中Ozone的缺失值比率超过20%。右图反映了缺失值的模式，红色表示没有删除，蓝色表示删除。从图中可以看出，仅Ozone变量缺失值占了22.9%，仅Solar. R变量缺失值占了3.3%，两个变量都缺失的占了1.3%。数据完整的观测值占72.5%。

此外，marginplot()函数可以帮助我们可视化缺失值的分布。

marginplot(airquality[1:2])

在下图中，湖蓝色圆圈表示未缺失值，红色的实心点表示缺失值，而深紫色点表示两个变量都缺失。图左侧的红色方框图显示了在Ozone含有缺失值的情况下Solar.R的分布。蓝色方框图显示去除Ozone的缺失值后Sloar.R的分布。图表底部的方框图正好相反，反映了Solar.R含有缺失值和去除缺失值时Ozone的分布。

04

小结

还是那句话，“统计是一门严谨的科学”。只有做好了原始数据的清洗和处理，才能保证后续结果的准确性，否则“地基”没有打好，用再好的钢筋水泥也不能搭建起来稳固的房子。选好数据，处理好数据，选好方法，用对统计方法，只有这样，才是一个合格的“数据分析师”。好了，关于离群值和缺失值的处理我们今天先讲到这里，我们的《临床模型构建》系列文章也快要接近尾声了，不知道你的学习进度怎么样呢？