读书笔记_第十八章

#处理缺失数据的高级方法

#解决缺失数据的两种方法：
#1：删除含有缺失数据的实例
#2：用合理的替代之替换缺失值
#不管是哪种方法，最后的结果都是没有缺失值的数据集

#缺失数据的分类：
#1 完全随机缺失 MCAR
#若某变量的缺失数据与其他任何观测或未观测变量都不相关
#2 随机缺失 MAR
#若某变量上的缺失数据与其他观测变量相关，而与自己的未观测值不相关
#3 非随机缺失 NMAR
#若缺失数据不属于MCAR,或MAR

#处理缺失值的步骤：
#第一步：识别缺失数据
#第二步：检查导致数据缺失的原因
#第三步：删除包含缺失值的实例或用合理的数值代替(插补)缺失值

install.packages(c("VIM","mice"))

#识别缺失值

#NA,不可得，代表缺失值；NaN,不是一个数，代表不可能值；Inf 正无穷；-Inf 负无穷
x1 <- NA  #返回x1=NA
is.na(x1) #返回TRUE
is.nan(x1) #返回FALSE
is.infinite(x1) #返回FALSE

x2 <- 0/0 #返回x2=NaN
is.na(x2) #返回TRUE
is.nan(x2) #返回TRUE
is.infinite(x2) #返回FALSE

x3 <- 1/0 #返回x3=Inf
is.na(x3) #返回FALSE
is.nan(x3) #返回FALSE
is.infinite(x3) #返回TRUE

#范例 

#将VIM包中的sleep数据集，加载进内存，放入全局变量
data(sleep,package="VIM")
str(sleep)# 返回62行，10列

#complete.cases()可以用来识别矩阵或数据框中没有缺失值的行
#NA，NaN识别为缺失值
#若每行包含完整的实例，则返回TRUE的逻辑向量，等价于1
#若每行有一个或多个缺失值，则返回false，等价于0

#返回所有不包含缺失值的行
sleep[complete.cases(sleep),]

#返回所有包含缺失值的行
sleep[!complete.cases(sleep),]

#返回12，Dream列包含缺失值的行数为12个
#is.na，转换成TRUE,FALSE的向量，统计TRUE,对应缺失值
sum(is.na(sleep$Dream)) 
mean(is.na(sleep$Dream))
#返回0.1935484，等价于12/64=0.1935483，在Dream列上的缺失比例

nrow(sleep[!complete.cases(sleep),]) #返回 20，包含缺失值的行数
mean(!complete.cases(sleep))#返回0.3225806，等价于20/62=0.3225806

#探索缺失值模式

#范例一  列表显示缺失值
library(lattice)
library(mice)
data(sleep,package="VIM")
md.pattern(sleep)#输出的表格更直观
#结果解读：
#0,相当于is.na返回true,表示变量的列中有缺失值
#1，相当于is.na返回False,表示变量的列中没有缺失值
#第一列，表示各缺失值模式，对应的观测值总个数
#最后一列，表示当前缺失值模式下，一个实例存在缺失值列的个数
#最后一行，表示每个变量中缺失值的总个数

#数据集中总共包含的缺失值数：
#第一列*最后一列，求和：
#0 ,2,1,1,3,2,2,3
#42,9,3,2,1,1,2,2
#0*42+2*9+1*3+1*2+3*1+2*1+2*2+3*2=38

#sleep列下存在两个零，每一个零对应模式的观测值总数，2+2=4
#Dream列下存在三个零，第一零，对应9个观测值，第二个零，对应1个观测值，第三个零对应2个观测值 9+1+2=12

#范例二  aggr函数，图形探究缺失数据
library("colorspace")
library("grid")
library("data.table")
library("VIM")

aggr(sleep,prop=FALSE,numbers=TRUE)
#结果解读：
#左图，直方图
#横坐标对应各变量名
#纵坐标对应各变量名下，缺失值总个数，等价于最后一行的值

#右图，方格图
#蓝色表示无缺失值，红色表示有缺失值
#横坐标对应各变量名，纵坐标对应各个模式下，存在缺失值的列的观测值个数，等价于第一列

aggr(sleep,prop=TRUE,numbers=TRUE)
#将生成相同图形，用比例代替了计数

aggr(sleep,prop=TRUE,numbers=FALSE)
#将生成相同图形，用比例代替了计数,并且删除了数值型标签

#范例三 

#可能没有纵横坐标，清空画图区，然后重新画，就解决了
matrixplot(sleep)
#结果解释：
#矩阵图，色块，颜色深浅
#纵轴，index,对应sleep数据集的行号
#横轴，变量名，横向一条记录，对应sleep数据集对应行号的记录
#红色表示当前观测值，当前变量对应为空值
#数红色窄条的个数，对应当前变量缺失值总数
#数值型数据被转换成[0,1]区间的灰度表示，灰度越浅表示数值小，灰度越深表示数值大

dim(sleep[c("Gest","Dream")]) #返回62行，2列
marginplot(sleep[c("Gest","Dream")],pch=c(20),
           col=c("darkgray","red","blue"))
#结果解读：
#图形的主体是两变量都完整的散点图
#Gest与Dream呈负相关

#左边界的箱线图展示了包含gest值(灰色)，以及不包含gest值(红色)的Dream变量分布
#缺失Gest数据时，动物的Dream时间比一般更长

#在底部边界上，Gest与Dream间的关系反过来了
#两个变量均有缺失值的观测个数在两边界交叉处用蓝色输出，左下角为0

#用相关性探索缺失值

#影子矩阵
#你可以用指示变量替代数据集中的数据，1表示缺失，0表示存在，这样生成的矩阵有时被称为因子矩阵
#求这些指示变量之间和它们与初始(可观测)变量之间的相关性，有助于观察哪些变量常一起缺失
#有助于观察哪些变量常一起缺失

#范例

#is.na(sleep)对缺失数据识别，缺失值对应TRUE,非缺失值对应FALSE
#取绝对值时，TRUE对应1，FALSE对应0
#影子矩阵 matrix 
class(abs(is.na(sleep))) #返回matrix
x <- as.data.frame(abs(is.na(sleep))) #影子矩阵强制转换成dataframe

#对照查看
head(sleep,n=5)
head(x,n=5)

#apply(x,2,sum)，针对列求和
#which(apply(x,2,sum)),返回和值>0的列名，即对应存在缺失值的列
#再根据需要的列名访问数据框x对应列，以及列值，提取存在缺失值的列，以及列值
y <- x[which(apply(x,2,sum)>0)]

#计算缺失列的相关系数
cor(y)
#结果解读：
#Dream，NonD，相关系数=0.90711474，说明常常一起缺失
#Sleep，NonD，相关系数=0.48626454，相比较而言一起缺失的可能性，小些
#Sleep，Dream，相关系数=0.20370138，相比较而言一起缺失的可能性，更小些

#sleep对应行变量，对应sleep全部10个预测变量
#y对应列变量，对应y全部含有缺失值列的变量
c <- cor(sleep,y,use="pairwise.complete.obs")
class(c) #返回matrix
c[order(c[,1]),] #基于第一列NonD,升序排列
c[order(-c[,1]),] #基于第一列NonD,降序排列
#结果解读：
#含有缺失值列，自相关系数为NA,即NonD与NonD的相关系数为NA
#针对第一列NonD,按照相关系数降序排序
#cor(EXP,NonD)=0.24546836
#cor(BodyWgt,NonD)=0.22682614
#cor(Gest,NonD)=0.20239201
#exp越大，BodyWgt越大，Gest越大，其Nond值更容易为空值
#表中的相关系数并不特别大，表明数据时MCAR的可能性比较小，更可能是MAR
#当缺乏强有力的外部证据时，我们通常假设数据时MCAR或者MAR

#理解缺失数据的来由和影响
#目的：
#1：分析生成缺失数据的潜在机制
#2：评价缺失数据对回答实质性问题的影响

#需要弄清楚：
#1：缺失数据的比例多大？
#2：缺失数据是否集中在少数几个变量上，抑或广泛存在？
#3：缺失是随机产生的吗？
#4：缺失数据间的相关性或与可观测数据间的相关性，是否可以表明产生缺失值的机制

#处理方法
#如果缺失值数据集中在几个相对不太重要的变量上，那么你可以删除这些变量，
#再进行正常的数据分析

#如果一小部分数据(如小于10%)随机分布在整个数据集中MCAR,那么你可以
#分析数据完整的实例，这样仍可以得到可靠且有效的结果

#方法一：理性处理不完整数据
#推理方法会根据变量间的数学或者逻辑关系来填补或恢复缺失值
#推理法常常需要创造性和想法，同时还需要许多数据处理技巧，而且数据的恢复
#可能是准确的，或者是近似的

#范例   行删除，计算相关系数

options(digits=1) #全局变量设置，保留小数位后1

#将sleep中存在缺失值的行，全部删除后，剩余的无缺失值行，参与相关系数计算
#newdata <- mydata[complete.cases(mydata),]
#newdata <- na.omit(mydata)
#两条命令等价，对于存在缺失值的行进行完全删除
cor(na.omit(sleep))

#范例  行删除，应用线性模型
fit <- lm(Dream~Span+Gest,data=na.omit(sleep))
summary(fit)
#基于线性模型，返回残差自由度39,样本数n=39+3=42
#残差的自由度就是n-k啊。你的例子中n=42,k=3（截距项一个，x两个）
fit$df.residual 

#结果解读：
#Pr(>|t|) =0.041，落入拒绝域，显著不等于零
#-0.004394,Gest的回归系数，说明Gest每增加一单位， Dream将减少0.004394
#Residual standard error: 1 on 39 degrees of freedom

#lm使用有限的行删除法定义
#44=42+2
#42对应无任何缺失值的观测值个数
#2，对应变量Dream,Span,Gest全部无缺失值，但其他列有缺失值的观测值个数
fit <- lm(Dream~Span+Gest,data=sleep)
summary(fit)
#基于线性模型，返回残差自由度41,样本数n=41+3=44
#残差的自由度就是n-k啊。你的例子中n=44,k=3（截距项一个，x两个）
fit$df.residual 
#结果解读：
#Residual standard error: 1 on 41 degrees of freedom
#(18 observations deleted due to missingness)

#多重插补 MI
#多重插补，是一种基于重复模拟的处理缺失值的方法
#在面对复杂的缺失值问题时，MI是最常选用的方法

#范例 
library(lattice)
library(mice)
data(sleep,package="VIM") #加载VIM包中的sleep数据进入全局环境变量

#mice函数首先从一个包含缺失数据的数据框开始，然后返回一个包含多个(默认为5个)
#完整数据集的对象
#每个完整数据数据都是通过对原始数据框中的缺失数据进行插补而生成的
#因为每个插补有随机的成分，因此每个完整数据集都略有不同
imp <- mice(sleep,seed=1234)
class(imp)  #返回"mids",等价于Multiply imputed data set

#imp对象是一个包含m个插补数据集的列表对象，同时还含有完成插补过程的信息
#默认为5
imp
#结果解读：
#Number of multiple imputations:  5 

#Imputation methods:
#BodyWgt BrainWgt     NonD    Dream    Sleep     Span     Gest 
#""       ""         "pmm"    "pmm"    "pmm"    "pmm"    "pmm" 
#Pred      Exp   Danger 
#""       ""       "" 
#BodyWgt，BrainWgt，Pred，Exp，Danger的插补方法对应空，显示没有缺失值，无需插补
#NonD,Dream,Sleep,Span,Gest的插补方法对应pmm,存在缺失值，采用pmm方法
#pmm,预测均值匹配法来处理每个含缺失值数据的变量

#PredictorMatrix:
#预测变量矩阵展示了进行插补过程的含有缺失数据的变量，它们利用了数据集中
#其他变量的信息。
#在矩阵中，行代表插补对象，列代表插补提供信息的变量，
#1和0分别表示使用和未使用

imp$imp$Dream
nrow(imp$imp$Dream) #返回12，Dream列共计12个缺失值
#结果解读：
#列编号，1~5，表示对应5次插补
#行编号，在Dream列对应缺失值的观测值上，行号分别对应12种动物

#complete函数，可以观察m个插补数据集中的任意一个
#将所有缺失值，在第三次的值取出来，填充到原始数据集的对应位置，所拼凑的完整数据集
#以Dream列为例
sleep_3 <- complete(imp,action=3)
Dream_na <- c(1,3,4,14,24,26,30,31,47,53,55,62) #保存sleep中Dream列值为空的行号
a <- sleep_3[Dream_na,"Dream"] #提取拼凑后的数据集，在Dream列空值的情况下，填充进去的数据
b <- imp$imp$Dream[,3] #提取imp对象，针对Dream列，第三次插补值
ifelse((sleep_3[Dream_na,"Dream"] == imp$imp$Dream[,3]),print("Same"),print("Not Same"))
identical(a,b) #判断两个向量是否完全相等
all(a==b) #判断两个向量是否完全相等

#with函数可依次对每个完整数据集应用统计模型
fit <- with(imp,
            lm(Dream~Span+Gest))

#-------异常处理代码-----------
#执行过程中，报错Error in mutate_impl(.data, dots) : 
#Evaluation error: `as_dictionary()` is defunct as of rlang 0.3.0.
#Please use `as_data_pronoun()` instead.
#解决方法：
#查看当前dplyr包的版本
packageVersion("dplyr") #返回‘0.7.4’

remove.packages("dplyr") #移除旧版本‘0.7.4’
install.packages("dplyr") #安装新版本
packageVersion("dplyr") #返回‘0.7.8’
#-----------------------------------

#正常代码段
#pool函数将这些单独的分析结果整合为一组结果
#最终模型的标准误和p值都将准确的反映出由于缺失值和多重插补而产生的不确定性
pooled <- pool(fit)
summary(pooled)
#结果解读：
#e，2.71828,约等于10，e-01约等于0,01
#span的P=8e-01,约等于0.8，未落入拒绝域，显著为0
#Gest p=8e-03,约等于0.008，落入拒绝域，不显著为0
#当控制寿命不变时，Gest与Dream成显著的，负相关关系

#成对删除
#处理缺失值数据集时，成对删除常作为行删除的备选方法使用
#对于成对删除，观测只是当它含缺失值变量涉及某个特定分析时才会被删除
#任何两个变量的相关系数都只利用了仅这两个变量的可用观测，忽略其他变量
cor(sleep,use="pairwise.complete.obs")
#结果解读：
cor(sleep$BodyWgt,sleep$BrainWgt) #不包含空值0.93

x <- sleep[!is.na(sleep$NonD),c("BodyWgt","NonD")]
dim(x) #返回48行，2列
cor(x$BodyWgt,x$NonD) #返回

x <- sleep[!is.na(sleep$NonD),c("BodyWgt","NonD")]
dim(x) #返回48行，2列
cor(x$BodyWgt,x$NonD) #返回

y <- na.omit(sleep[,c("Dream","NonD")])
dim(y) #返回48行，2列
cor(y$Dream,y$NonD) #返回0.5

#简单非随机插补
#即用某个值（如均值，中位数，众数）来替换变量中的缺失值
#简单插补的一个优点是，解决缺失值问题时不会减少分析过程中可用的样本量
#若缺失数值的数目非常大，那么简单插补很可能会低估标准差，曲解变量间的相关性
#并会生成不正确的统计检验的p值
#与成对删除一样，不建议使用