面對一份數據，我們會做數據的缺失值檢測和分析，根據數據的缺失程度，以知道數據的完整性和可用性。實際的數據，絕大部分會有缺失值現象。缺失值的產生與諸多因素有關聯，例如：數據採集不成功，數據採集成功了但是數據確實沒有值，數據的值受其它因素控制等。面對有缺失值的數據，我們要怎麼處理、分析和應用，是我們數據工作者要思考和實踐的命題。這篇文章，總結和分享R語言如何檢測數據中的缺失值，以及針對不同情形下處理缺失值的方法。
R語言可以使用is.na()函數來判斷數據集的取值是否為缺失值，並且利用一些匯總函數可以概述數據整體的缺失率、每個樣本的缺失率和每個變量的缺失率。R語言也提供可視化工具來直觀展示變量的缺失分布情況。
【代碼片段】
# R包library(pacman)p_load(tidyverse) # 數據科學套件p_load(mice) # 缺失值分析包p_load(VIM)
# 導入數據telecom_data <- read_csv('./data/telecom.csv')# 數據檢視telecom_data %>% glimpse()# 10行，5列# 4個字符型變量，1個雙精度型變量
# 缺失值檢驗is.na(telecom_data)# 數據整體缺失率sum(is.na(telecom_data))mean(is.na(telecom_data))
# 自定義缺失率計算函數var_na_ratio <- function(x){ return(mean(is.na(x)))}# 變量缺失率apply(telecom_data, 2, var_na_ratio)# 樣本缺失率apply(telecom_data, 1, var_na_ratio)【結果】
我們可以發現數據集為10行、5列，50個單元，整體缺失情況是5個缺失值，缺失率為10%。變量的缺失率情況，MonthlyCharges為30%、TotalCharges和PaymentMethod為10%，樣本的缺失情況，3個樣本有1個缺失值，1個樣本有2個缺失值。
【代碼片段】
md.pattern(telecom_data)aggr_plot <- aggr(telecom_data, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(data), cex.axis=.7, gap=3, ylab=c("Histogram of missing data","Pattern"))【結果】
橫向觀察，表示有6個樣本沒有缺失，2個樣本有個1個缺失值，1個樣本有2個缺失依次類推；縱向觀察，TotalChanges有1個缺失值，PaymentMethod有1個缺失，MonthlyCharges有3個缺失值，總共是有5個缺失。
左邊的圖按著變量缺失率做降序排列；右邊的圖從下往上看，第一行，表示60%的樣本無缺失值，第二行，表示20%的樣本在變量MonthlyCharges有1個缺失值，以此類推。
我們可以把缺失值化為兩種類型：
MCAR：完全隨機的缺失，這是缺失值的理想場景。
MNAR：非完全性隨機缺失，這是常見情形，這個時候進一步檢查數據收集過程並嘗試理解數據缺失的原因可能是明智而必要的。數據集含有太多缺失會帶來嚴重問題。在大數據集情況，通常我們會保留變量缺失率不超過5%或者樣本缺失率不超過5%的樣本，當然，這個最佳閾值，可以根據實際情況調整。
【代碼片段】
# 以自帶數據集airquality為例data <- airqualitydata[4:10,3] <- rep(NA,7)data[1:5,4] <- NA
data <- data[-c(5,6)]summary(data)
pMiss <- function(x){sum(is.na(x))/length(x)*100}apply(data,2,pMiss)apply(data,1,pMiss)【結果】
可以發現變量Ozone的缺失率有24.2%，而其它變量的缺失率是5%以下，我們可以刪除變量Ozone，這也是一種特徵選擇方法，基於數據的可用性來過濾特徵。因為，我們所觀察的數據集變量4個，只要一個有兩個或者以上變量缺失，就出現了樣本50%及以上的缺失率，因此這樣的樣本可以不予納入研究。
上面的這些做法，實際屬於一種常用R語言處理數據缺失值的方法，即根據變量的缺失率和樣本的缺失率，來過濾不符合閾值的變量和樣本集。若是在大數據和缺失率小比例的背景下，這種方法是有效的，也是簡單易行的。
有時候，我們為了保留變量和樣本數，需要採用別的方法來處理缺失值，常用的手段就是插補，不管是基於統計學裡面平均值、中位數或者眾數，還是基於某一種學習模型來完成缺失填充，這些都可以看作是插補法，只是方式不同而已。我們可以利用R語言mice包來做缺失值的插補法。
【代碼片段】
# 缺失值的插補tempData <- mice(data,m=5,maxit=50,meth='pmm',seed=500)completedData <- complete(tempData,1)sum(is.na(data))sum(is.na(completedData))缺失的值被替換為五個數據集中第一個數據集的輸入值。如果希望使用另一個數據集的結果，只需更改complete()函數中的第二個參數。
【結果】
mice包提供了很多進行缺失值插補的方法，具體方法查詢如下：
【代碼片段】
【結果】
關於這些方法的含義，可以查看mice包的幫助文檔。大家做一下，若是要用隨機森林做插補法，需要用哪個方法？
總結：本文分享缺失值的檢測方法和處理策略，以及缺失值兩種類型，還有R 語言處理缺失值的mice包。
【完整代碼片段】
library(pacman)p_load(tidyverse) p_load(mice) p_load(VIM)
telecom_data <- read_csv('./data/telecom.csv')telecom_data %>% glimpse()
is.na(telecom_data)sum(is.na(telecom_data))mean(is.na(telecom_data))
var_na_ratio <- function(x){ return(mean(is.na(x)))}apply(telecom_data, 2, var_na_ratio)apply(telecom_data, 1, var_na_ratio)
md.pattern(telecom_data)aggr_plot <- aggr(telecom_data, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(data), cex.axis=.7, gap=3, ylab=c("缺失值直方圖","缺失值模式"))
data <- airqualitydata[4:10,3] <- rep(NA,7)data[1:5,4] <- NA
data <- data[-c(5,6)]summary(data)
pMiss <- function(x){sum(is.na(x))/length(x)*100}apply(data,2,pMiss)apply(data,1,pMiss)
tempData <- mice(data,m=5,maxit=50,meth='pmm',seed=500)completedData <- complete(tempData,1)sum(is.na(data))sum(is.na(completedData))
methods(mice)
summary(tempData)dim(tempData$imp$Ozone)tempData$methxyplot(tempData,Ozone ~ Wind+Temp+Solar.R,pch=18,cex=1)densityplot(tempData)stripplot(tempData, pch = 20, cex = 1.2)tempData
modelFit1 <- with(tempData,lm(Temp~ Ozone+Solar.R+Wind))summary(pool(modelFit1))參考資料：
1https://datascienceplus.com/imputing-missing-data-with-r-mice-package/
2 https://raw.githubusercontent.com/dataoptimal/posts/master/data%20cleaning%20with%20R%20and%20the%20tidyverse/telecom.csv
好書推薦
1 R語言做商業智能，助你提高商業生產率
2 用RStudio做數據分析
3 用R、tidyverse和mlr做機器學習
4 推斷統計與數據科學，moderndive和tidyverse包
公眾號推薦