作者：科研猫 | 西红柿
责编：科研猫 | 馋猫
无论是前瞻性数据收集还是回顾性数据收集，数据集中通常都会出现离群值或缺失值。对于统计学家来说，离群值和缺失值通常是一个棘手的问题，如果处理不当可能会导致错误。离群值可能会导致我们的结果偏离真实结果，而缺失值造成的信息损失可能会导致建模失败。因此，在执行数据分析之前，正确识别离群值并处理缺失值非常重要。本推文讨论的内容应该在建模之前执行。虽然本推文在整个统计模型系列中较为置后，却至关重要，望警醒。
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()查看这组患者的收缩压范围。
当处理含有缺失值的数据时，要设置参数na.rm = TRUE。
上述方法可以帮助我们识别最大值或最小值，但有时极限值并不是单独出现的，而是在聚类中，因此上述方法识别异常值是不够的。在实际的研究背景下，我们通常根据变量的均值和标准差，或中位数和四分位数（Tukey方法）来定义数据的异常值。例如，我们可以设置大于或小于mean±3sd均为异常值。当然，我们也可以对分类变量的某个值进行异常判断。例如，性别值为1=男性，2=女性。如果赋值为3，则为异常值。这里我们介绍一个自定义函数。该函数根据四分位Tukey方法判断异常值，有效地避免了极限值对均值和标准差的影响。函数如下：
1outlierKD <-function(dt, var) {
2var_name <- eval(substitute(var),eval(dt))
3tot <- sum(!is.na(var_name))
4na1 <- sum(is.na(var_name))
5m1 <- mean(var_name, na.rm =T)
6par(mfrow=c(2,2), oma=c(0,0,3,0))
7boxplot(var_name, main="With outliers")
8hist(var_name, main="With outliers", xlab=NA, ylab=NA)
9outlier <- boxplot.stats(var_name)$out
10mo <- mean(outlier)
11var_name <- ifelse(var_name %in% outlier,NA, var_name)
12boxplot(var_name, main="Without outliers")
13hist(var_name, main="Without outliers", xlab=NA, ylab=NA)
14title("Outlier Check", outer=TRUE)
15na2 <- sum(is.na(var_name))
16cat("Outliers identified:", na2 - na1,"\n")
17cat("Propotion (%) of outliers:", round((na2 - na1) / tot*100,1),"\n")
18cat("Mean of the outliers:", round(mo,2),"\n")
19m2 <- mean(var_name, na.rm =T)
20cat("Mean without removing outliers:", round(m1,2),"\n")
21cat("Mean if we remove outliers:", round(m2,2),"\n")
22response <- readline(prompt="Do you want to remove outliers and to replace with NA? [yes/no]: ")
23if(response =="y"| response =="yes"){
24dt[as.character(substitute(var))] <- invisible(var_name)
25assign(as.character(as.list(match.call())$dt), dt, envir = .GlobalEnv)
26cat("Outliers successfully removed","\n")
27return(invisible(dt))
28}else{
29cat("Nothing changed","\n")
30return(invisible(var_name))
31}
32}
自定义函数只有两个参数，第一个参数是数据集的名称，第二个参数是变量名；只要正确替换数据集和变量名，读取就可以直接运行代码。这里我们是以箱形图的外值为离群值，我们还可以根据专业知识重新设置离群值的定义，比如大于或小于mean±3sd。在函数结束时，还将设置用户输入的代码。用户可以通过键入“yes”或“no”来确定是否消除数据集中函数识别的异常值。
下面我们模拟一组数据来验证这个自定义异常值识别函数的功能。
02
缺失值的处理
数据丢失在临床研究中很常见。例如，护士在收集数据时，可能会因为工作繁忙而忘记记录某个时间点的尿量；当研究人员想研究乳酸变化对死亡率的影响时，患者可能只监测某个时间点的血乳酸值。缺乏数据的其他原因还包括编码错误、设备故障和调查研究中的应答者没有应答等。在统计软件包中，一些函数（如Logistic回归）可能会自动删除丢失的数据。如果只有少量的不完全观测，那么这种处理就不会有太大问题。
但是，当存在大量包含缺失值的观测值时，这些函数中的默认行删除可能会导致大量信息丢失。在这种情况下，分析人员应该仔细研究数据丢失可能导致的机制，并找到适当的处理方法。
如何处理缺失值是临床统计学家头疼的问题，所以我们也应该予以重视。数据的缺失或缺失程度直接影响到数据的质量，而数据的质量最终影响到我们的研究成果。如果对缺失数据的处理不当，很可能导致整个统计分析失败。本推文介绍了在R中如何处理丢失的数据，并介绍了处理丢失数据的一些基本技巧。
在R中，“NA”表示为一个缺失的值。当将带有空单元格的Excel表导入R控制台时，这些空单元格将被NA替换。这与STATA用“.”替换“空单元格”不同。R中的数值变量和字符变量使用相同的缺失值符号。R提供一些函数来处理缺失值。要确定向量是否包含缺少的值，可以使用is.na（）函数。“is.na（）”函数是用于确定元素是否为na类型的最常用方法。它返回与传入参数长度相同的对象，并且所有数据都是逻辑值（FALSE或TRUE）。假设我们有6个病人，但是只记录了4个值，而缺少了2个。
03
缺失值的可视化
缺失值的可视化可以帮助我们更直观地观察数据集中的缺失值，这将有助于我们以后对缺失值进行插值。在本推文中，笔者将主要向读者介绍VIM包的使用。以下的演示数据集是R语言的内置数据集"airquality"。
"airquality"数据集包含了153个观测值和6个变量。从以上结果中，我们可以看到该数据集中有缺失值。在可视化之前，首先使用mice包中的md.pattern()函数探索缺失的数据模式。
在输出表格中，“1”表示非缺失值，“0”表示缺失值。第一列显示了唯一缺失数据模式的数目。在我们的例子中，111个观测值没有缺失数据，35个观测值仅在Ozone变量中有缺失数据，5个观测值仅在Solar. R变量中有缺失数据。最右边的一列显示了特定缺失模式中缺失变量的数目。例如，如果第一行中没有缺失值，则显示为“0”。最后一行计算每个变量缺失值的数量。例如，“Wind”变量没有缺失值，显示“0”，而Ozone变量有37个缺失值。在研究中，一些含有更多缺失值的变量可能会被剔除。显然，表格可以提供有用的参考信息。
下面我们调用VIM包来实现缺失值的可视化。研究缺失数据模式对于选择合适的插值方法来估计缺失值是必要的。因此，需要在插值操作之前执行可视化工具，并且通常应该在缺失数据插值之后进行诊断，以确定插值是否合理。可用于可视化丢失数据的函数如下:aggr()、matrixplot()、scattMiss()和marginplot()，以下我们主要演示了aggr()和marginplot()两个函数。
aggr()函数的作用是帮助我们可视化缺失的值。左图是缺失值比例直方图。从下图中可以看出Ozone和Solar. R有缺失值，其中Ozone的缺失值比率超过20%。右图反映了缺失值的模式，红色表示没有删除，蓝色表示删除。从图中可以看出，仅Ozone变量缺失值占了22.9%，仅Solar. R变量缺失值占了3.3%，两个变量都缺失的占了1.3%。数据完整的观测值占72.5%。
此外，marginplot()函数可以帮助我们可视化缺失值的分布。
在下图中，湖蓝色圆圈表示未缺失值，红色的实心点表示缺失值，而深紫色点表示两个变量都缺失。图左侧的红色方框图显示了在Ozone含有缺失值的情况下Solar.R的分布。蓝色方框图显示去除Ozone的缺失值后Sloar.R的分布。图表底部的方框图正好相反，反映了Solar.R含有缺失值和去除缺失值时Ozone的分布。
04
小结
还是那句话，“统计是一门严谨的科学”。只有做好了原始数据的清洗和处理，才能保证后续结果的准确性，否则“地基”没有打好，用再好的钢筋水泥也不能搭建起来稳固的房子。选好数据，处理好数据，选好方法，用对统计方法，只有这样，才是一个合格的“数据分析师”。好了，关于离群值和缺失值的处理我们今天先讲到这里，我们的《临床模型构建》系列文章也快要接近尾声了，不知道你的学习进度怎么样呢？
文末说个正事儿
由于微信平台开始试行乱序推送，为了让您第一时间看到我们的推送，强烈建议将“科研猫”设为星标。首先进入科研猫公众号，第一步点击右上方···，第二步点击设为星标，即可完成订阅啦。
本期干货
-离群值及缺失值处理代码 -
领取方法
关注“科研猫”公众号
公众号主页点击“更多信息”-“联系客服”领取干货
网友评论