（數據科學學習手札58）在R中處理有缺失值數據的高級方法

時間 2019-12-07

標籤數據科學學習手札處理缺失高級方法简体版

原文原文鏈接

1、簡介算法

　　在實際工做中，遇到數據中帶有缺失值是很是常見的現象，簡單粗暴的作法如直接刪除包含缺失值的記錄、刪除缺失值比例過大的變量、用0填充缺失值等，但這些作法會很大程度上影響原始數據的分佈或者浪費來之不易的數據信息，所以怎樣穩當地處理缺失值是一個持續活躍的領域，貢獻出衆多巧妙的方法，在不浪費信息和不破壞原始數據分佈上試圖尋得一個平衡點，在R中用於處理缺失值的包有不少，本文將對最爲普遍被使用的mice和VIM包中經常使用的功能進行介紹，以展示處理缺失值時的主要路徑；bootstrap

2、相關函數介紹數組

2.1 缺失值預覽部分app

　　在進行缺失值處理以前，首先應該對手頭數據進行一個基礎的預覽：dom

　　一、matrixplot函數

　　效果相似matplotlib中的matshow，VIM包中的matrixplot將數據框或矩陣中數據的缺失及數值分佈以色彩的形式展示出來，下面是利用matrixplot對R中自帶的airquality數據集進行可視化的效果：spa

rm(list=ls())
library(mice)
library(VIM)

#以airquality數據爲例
data(airquality)
data <- airquality
#利用矩陣熱力圖查看缺失狀況，紅色表明缺失
matrixplot(data)

　　紅色部分即表明數據缺失值所在位置，經過這個方法，能夠在最開始對數據總體的缺失狀況有一個初步認識，如經過上圖能夠一眼看出變量Ozone缺失狀況較爲嚴重；rest

　　二、marginplot與marginmatrixcode

　　缺失值是否符合徹底隨機缺失是在對數據進行插補前要着重考慮的事情，VIM中的marginplot包能夠同時分析兩個變量交互的缺失關係，依然以airquality數據爲例：orm

marginplot(data[c(1,2)])

　　如上圖所示，經過marginplot傳入二維數據框，這裏選擇airquality中包含缺失值的前兩列變量，其中左側對應變量Solar.R的紅色箱線圖表明與Ozone缺失值對應的Solar.R未缺失數據的分佈狀況，藍色箱線圖表明與Ozone未缺失值對應的Solar.R未缺失數據的分佈狀況，下側箱線圖同理，當同一側紅藍箱線圖較爲接近時可認爲其對應考察的另外一側變量缺失狀況比較貼近徹底隨機缺失，這種狀況下能夠放心大膽地進行以後的插補，不然就不能冒然進行插補；

　　與marginplot功能類似，marginmatrix在marginplot只能展示兩個變量的基礎上推廣到多個變量兩兩之間，效果相似相關性矩陣圖：

marginmatrix(data)

　　三、自編函數計算各個變量缺失比例

　　爲了計算出每一列變量具體的缺失值比例，能夠自編一個簡單的函數來實現該功能：

> #查看數據集中每一列的缺失比例
> miss.prop <- function(x){sum(is.na(x))/length(x)}
> apply(data,2,miss.prop)
     Ozone    Solar.R       Wind       Temp      Month        Day 
0.24183007 0.04575163 0.00000000 0.00000000 0.00000000 0.00000000

　　經過自編的函數miss.prop，能夠對每一個變量中缺失值所佔比例有個具體的瞭解；

2.2 mice函數

　　mice包中最核心的函數是mice()，其主要參數解釋以下：

data: 傳入待插補的數據框或矩陣，其中缺失值應表示爲NA

m: 生成插補矩陣的個數，mice最開始基於gibbs採樣從原始數據出發爲每一個缺失值生成初始值以供以後迭代使用，而m則控制具體要生成的完整初始數據框個數，在整個插補過程最後須要利用這m個矩陣融合出最終的插補結果，若m=1，則惟一的矩陣就是插補的結果；

method: 這個參數控制了傳入數據框中每個變量對應的插補方式，完好失值的變量對應的爲空字符串，帶有缺失值的變量默認方法爲"pmm"，即均值插補

predictorMatrix: 由於mice中絕大部分方法是用擬合的方式以含缺失值變量以外的其餘變量爲自變量，缺失值爲因變量構建迴歸或分類模型，以達到預測插補的目的，而參數predictorMatrix則用於控制在對每個含缺失值變量的插補過程當中做爲自變量的有哪些其餘變量，具體用法下文示例中會詳細說明

maxit: 整數，用於控制每一個數據框迭代插補的迭代次數，默認爲5

seed: 隨機數種子，控制隨機數水平

在對缺失值插補過程當中，很是重要的是爲不一樣的變量選擇對應的方法，即method中對應的輸入，下表是每種算法對應的參數代號、適用數據類型和算法名稱：

`方法代號`	適用數值類型	對應的具體算法名稱
`pmm`	any	Predictive mean matching
`midastouch`	any	Weighted predictive mean matching
`sample`	any	Random sample from observed values
`cart`	any	Classification and regression trees
`rf`	any	Random forest imputations
`mean`	numeric	Unconditional mean imputation
`norm`	numeric	Bayesian linear regression
`norm.nob`	numeric	Linear regression ignoring model error
`norm.boot`	numeric	Linear regression using bootstrap
`norm.predict`	numeric	Linear regression, predicted values
`quadratic`	numeric	Imputation of quadratic terms
`ri`	numeric	Random indicator for nonignorable data
`logreg`	binary	Logistic regression
`logreg.boot`	binary	Logistic regression with bootstrap
`polr`	ordered	Proportional odds model
`polyreg`	unordered	Polytomous logistic regression
`lda`	unordered	Linear discriminant analysis
`2l.norm`	numeric	Level-1 normal heteroscedastic
`2l.lmer`	numeric	Level-1 normal homoscedastic, lmer
`2l.pan`	numeric	Level-1 normal homoscedastic, pan
`2l.bin`	binary	Level-1 logistic, glmer
`2lonly.mean`	numeric	Level-2 class mean
`2lonly.norm`	numeric	Level-2 class normal

　　在面對數據集具體狀況時，對插補方法進行微調是很必要的步驟，在上面鋪墊了這麼多以後，下面在具體示例上進行演示，並引入其餘的輔助函數；

2.3 利用mice進行缺失值插補——以airquality數據爲例

　　由於前面對缺失值預覽部分已經利用airquality進行演示，這裏就再也不贅述，直接進入正式插補部分，首先，咱們將data傳入mice函數，注意這裏設置maxit爲0以取得未開始迭代的初始模型參數：

#初始化插補模型,這裏最大迭代次數選0是爲了取得未開始插補的樸素模型參數
init <- mice(data, maxit = 0)

　　下面咱們來看看取得的須要進行調整的重要參數的初始狀況：

> #取得對於每個變量的初始插補方法
> methods <- init$method
> methods
  Ozone Solar.R    Wind    Temp   Month     Day 
  "pmm"   "pmm"      ""      ""      ""      ""

　　能夠看到對應缺失變量Ozone和Solar.R的插補擬合方法爲pmm，下面咱們把它們改爲CART決策樹迴歸：

#將變量Ozone的插補方法從pmm修改成norm
methods[c("Ozone")] <- 'cart'
#將變量Solar.R的插補方法從pmm修改成norm
methods[c("Solar.R")] <- 'cart'

　　接着咱們來查看predictorMatrix參數：

> #取得對每個變量進行擬合用到的變量矩陣，0表明不用到，1表明用到
> predM <- init$predictorMatrix
> predM
Ozone Solar.R Wind Temp Month Day
Ozone 0 1 1 1 1 1
Solar.R 1 0 1 1 1 1
Wind 1 1 0 1 1 1
Temp 1 1 1 0 1 1
Month 1 1 1 1 0 1
Day 1 1 1 1 1 0

　　這裏咱們認爲變量Month和Day是日期，與缺失變量無相關關係，所以將其在矩陣中對應位置修改成0使它們不參與擬合過程：

#調整參與擬合的變量
#這裏認爲日期對與其餘變量無相關關係，所以令變量Month與變量Day不參與對其餘變量的擬合插補過程
predM[, c("Month")] <- 0
predM[c("Month"),] <- 0
predM[, c("Day")] <- 0
predM[c("Day"),] <- 0

　　這樣咱們就完成了兩個重要參數的初始化，下面咱們進行正式的擬合插補：

#利用修改後的參數組合來進行擬合插補
imputed <- mice(data, method = methods, predictorMatrix = predM)

　　隨着程序運行完，咱們須要的結果便呼之欲出，但在取得最終插補結果前，爲了嚴謹起見，須要對模型的統計學意義進行分析，下面以Ozone爲例：

　　一、查看模型中Ozone對應的擬合公式：

> #查看Ozone主導的擬合公式
> imputed$formulas['Ozone']
$Ozone
Ozone ~ 0 + Solar.R + Wind + Temp
<environment: 0x000000002424d410>

　　能夠看到，Ozone對應的公式與前面predictorMatrix參數中通過修改的保持一致；

　　二、基於上述公式爲合成出的m=5個數據框分別進行擬合：

> #把上面的公式填入下面的lm()內
> fit <- with(imputed, lm(Ozone ~ Solar.R + Wind + Temp))
> 
> #查看fit中對應每個插補數據框的迴歸顯著性結果
> fit
call :
with.mids(data = imputed, expr = lm(Ozone ~ Solar.R + Wind + 
    Temp))

call1 :
mice(data = data, method = methods, predictorMatrix = predM)

nmis :
  Ozone Solar.R    Wind    Temp   Month     Day 
     37       7       0       0       0       0 

analyses :
[[1]]

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)

Coefficients:
(Intercept)      Solar.R         Wind         Temp  
  -53.36269      0.05791     -3.37688      1.51550  


[[2]]

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)

Coefficients:
(Intercept)      Solar.R         Wind         Temp  
  -55.70273      0.06189     -3.25456      1.50816  


[[3]]

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)

Coefficients:
(Intercept)      Solar.R         Wind         Temp  
  -55.45601      0.05659     -2.90911      1.47244  


[[4]]

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)

Coefficients:
(Intercept)      Solar.R         Wind         Temp  
  -70.00005      0.07008     -2.56784      1.59856  


[[5]]

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)

Coefficients:
(Intercept)      Solar.R         Wind         Temp  
  -72.08381      0.05252     -2.71741      1.67060

　　三、將上述5個擬合模型融合並查看融合後的顯著性水平

> #計算fit中模型的顯著性
> pooled <- pool(fit)
> summary(pooled)
                estimate   std.error statistic       df      p.value
(Intercept) -61.32105668 21.84851057 -2.806647 53.60143 6.964478e-03
Solar.R       0.05979673  0.02144304  2.788632 90.71532 6.449107e-03
Wind         -2.96516062  0.67509037 -4.392243 29.06277 1.361726e-04
Temp          1.55305117  0.23531131  6.599985 78.14524 4.468489e-09

　　能夠看到從截距項，到每個變量的p值都遠遠小於0.05，至少在0.05顯著性水平下每一個參數都具備統計學意義；

　　四、對5個合成出的數據框在缺失值位置進行融合，這裏須要用到新的函數complete，其主要有下面三個參數：

data: 前面mice函數輸出的結果

action: 當只但願從合成出的m個數據框中取得某個單獨的數據框時，能夠設置action參數，如action=3便表明取得m個數據框中的第3個

mild: 邏輯型變量，當爲TRUE時，會輸出包含所有m個合成數據框的列表

　　獲悉上列參數意義後，若只想抽取某個數據框如第3個：

result <- complete(imputed, action = 3)
matrixplot(result)

　　能夠看到，取回第3個數據框，每一個缺失值都已被補全，若但願獲得5個合成數據框的融合結果，則須要自編函數：

#取得全部合成數據框組成的列表
complete(imputed, mild = T)
all.imputed <- complete(imputed, action = 'all')
#對獲得的m個合成數據框缺失值部分求平均
avgDF <- function(data,all.imputed){
  baseDF <- all.imputed[[1]][is.na(data)]
  for(child in 2:length(all.imputed)){
    baseDF <- baseDF + all.imputed[[child]][is.na(data)]
  }
  data[is.na(data)] <- baseDF/length(all.imputed)
  
  return(data)
}

#獲得最終平均數據框
result <- avgDF(data,all.imputed)
matrixplot(result)