一、缺失值
1. 缺失值的识别和处理
缺失值存在的影响：无法进行求均值、求和等运算。（一般这些函数中都会有na.rm参数，设置成TRUE时可以自动去除缺失值）
x <- c(1,2,3,NA,NA,4)
x
# [1] 1 2 3 NA NA 4
sum(x)
# [1] NA
sum(x,na.rm = T)
# [1] 10
1.1 is.na函数的简单使用：
判断一个数据框中有没有缺失值，有多少缺失值
#判断有没有缺失值
is.na(x)
# [1] FALSE FALSE FALSE TRUE TRUE FALSE
#判断有多少缺失值（TRUE=1, FALSE=0）
sum(is.na(x))
# [1] 2
#感叹号起到取非的作用
x[!is.na(x)]
# [1] 1 2 3 4
#如果没有!就只返回NA
x[is.na(x)]
# [1] NA NA
以iris数据集为基础生成一个新的含缺失值的数据集iris_na
#在iris数据集的1-4列（数值型变量），1到最后一行中随机抽取五个值，用缺失值NA替换（i in 1:4是在1到4列中对i值做一个遍历）
iris_na <- iris
for(i in 1:4){
iris_na[sample(1:nrow(iris),5),i]=NA
}
1.2 使用sapply函数➕which来找到这5个缺失值
# which就是一个指针函数，意思是在哪里
sapply(iris_na[,1:4],function(x)which(is.na(x)))
# Sepal.Length Sepal.Width Petal.Length Petal.Width
# [1,] 28 12 10 19
# [2,] 54 39 20 24
# [3,] 68 76 45 26
# [4,] 132 117 119 69
# [5,] 141 128 127 143
# 每一列的缺失值的行数被显示出来
1.3 使用psych包中的describe函数(describe函数可以极大程度的返回了数据集的基本统计值)⚠️
describe(iris_na)
# vars n mean sd median trimmed mad min max range skew kurtosis se
# Sepal.Length 1 145 5.83 0.82 5.8 5.80 1.04 4.3 7.7 3.4 0.26 -0.68 0.07
# Sepal.Width 2 145 3.06 0.44 3.0 3.04 0.44 2.0 4.4 2.4 0.31 0.06 0.04
# Petal.Length 3 145 3.77 1.75 4.4 3.79 1.63 1.0 6.7 5.7 -0.32 -1.40 0.15
# Petal.Width 4 145 1.21 0.76 1.3 1.20 1.04 0.1 2.5 2.4 -0.12 -1.34 0.06
# Species* 5 150 2.00 0.82 2.0 2.00 1.48 1.0 3.0 2.0 0.00 -1.52 0.07
#前四列都只返回了145个值，说明存在5个缺失值
1.4 计算缺失值的比例
sapply(iris_na[,1:4],function(x)sum(is.na(x)/nrow(iris_na)))
# Sepal.Length Sepal.Width Petal.Length Petal.Width
# 0.03333333 0.03333333 0.03333333 0.03333333
# 缺失值的和/总行数
1.5 存在缺失值的回归分析
lm(Sepal.Length~Sepal.Width,data = iris_na)
# Call:
# lm(formula = Sepal.Length ~ Sepal.Width, data = iris_na)
# Coefficients:
# (Intercept) Sepal.Width
# 6.4824 -0.2105
lm(Sepal.Length~Sepal.Width,data = iris_na,na.action = na.omit)
# Call:
# lm(formula = Sepal.Length ~ Sepal.Width, data = iris_na, na.action = na.omit)
# Coefficients:
# (Intercept) Sepal.Width
# 6.4824 -0.2105
na.action默认是na.omit，因此在使用lm求回归时写不写这个参数都可以。
2. 填补缺失值
❗️在进行缺失值的填补的时候一定要紧密结合数据背景和专业背景，否则易出现矫枉过正，适得其反。
2.1 简单操作系列：使用均值填补缺失值
# 求均值
mean_value <- sapply(iris_na[,1:4],mean,na.rm=TRUE)
mean_value
# Sepal.Length Sepal.Width Petal.Length
# 5.844138 3.064138 3.785517
# Petal.Width
# 1.201379
for(i in 1:4){
iris_na[is.na(iris_na[,i]),i]=mean_value[i]
}
summary(iris_na)
describe(iris_na)
is.na（iris_na[,i]）是在iris_na的1到4列对进行遍历，判断数据框中有没有缺失值，返回的是TRUE或FALSE。iris_na[is.na（iris_na[,I],i]就相当于根据is.na（iris_na[,i]）, i的结果从iris_na中进行提取
如下方式是使用均值对所有缺失值进行填补，在缺失值存在于多个不同变量中，如肝癌，肺癌，胰腺癌等患者的术后存活时间不适合用所有数值的均值来填补，需要分类处理。
cancer <- data.frame(id=1:1000,sur_days=sample(100:1000,1000,replace = TRUE),
+ type=sample(c('colon','liver','lung'),1000,replace = TRUE),
+ treatment=sample(c('chemo','surg'),1000,replace = TRUE))
head(cancer)
# id sur_days type treatment
# 1 1 717 lung chemo
# 2 2 699 liver chemo
# 3 3 762 liver surg
# 4 4 400 liver chemo
# 5 5 599 lung chemo
# 6 6 933 lung surg
#生成缺失值
cancer[sample(1:1000,90),2] <- NA
mean_value <- tapply(cancer$sur_days,list(cancer$type,cancer$treatment),mean,na.rm=TRUE)
mean_value
# chemo surg
# colon 529.7751 550.5411
# liver 533.8811 563.5965
# lung 568.9627 517.8435
#缺失值填补
for(i in 1:3){
for(j in 1:2){
cancer$sur_days[is.na(cancer$sur_days)&cancer$type==rownames(mean_value)[i]
&cancer$treatment==colnames(mean_value)[j]]=mean_value[i,j]
}
}
summary(cancer)
describe(cancer)
2.2 填补缺失值的高级操作
演示数据集准备
mlbench包中的BostonHousing包
library(mlbench)
data("BostonHousing")
head(BostonHousing)
# crim zn indus chas nox rm age dis rad tax ptratio b lstat medv
# 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90 4.98 24.0
# 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14 21.6
# 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03 34.7
# 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94 33.4
# 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33 36.2
# 6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12 5.21 28.7
original_data <- BostonHousing
#随机生成随机数
set.seed(2021)
BostonHousing[sample(1:nrow(BostonHousing),80),'rad'] <- NA
BostonHousing[sample(1:nrow(BostonHousing),80),'ptratio'] <- NA
对缺失值进行操作
- 2.2.1
mice包⚠️
mice包中的md.pattern函数，可以反映数据集中缺失值情况
library(mice)
md.pattern(BostonHousing)
# crim zn indus chas nox rm age dis tax b lstat medv rad ptratio
# 359 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
# 67 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1
# 67 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1
# 13 1 1 1 1 1 1 1 1 1 1 1 1 0 0 2
# 0 0 0 0 0 0 0 0 0 0 0 0 80 80 160
0表示缺失，1表示不缺失。最下面一行是每一个变量中缺失值的总数。最下面一行显示rad和ptratio各有80个缺失，共有160个缺失。左边一列359表示359个观测都没有缺失，有67个观测rad变量缺失，ptratio没有缺失。有67个变量ptratio缺失，rad没有缺失。有13个观测red和ptratio都缺失。
对缺失值进行插补
mice_mod <- mice(BostonHousing[,!names(BostonHousing)%in% 'medv'],method = 'rf')
# rf是随机森林算法
# BostonHousing[,!names(BostonHousing)%in% 'medv']是为了禁止medv这一列进入随机森林模型
# iter imp variable
# 1 1 rad ptratio
# 1 2 rad ptratio
# 1 3 rad ptratio
# 1 4 rad ptratio
# 1 5 rad ptratio
# 2 1 rad ptratio
# 2 2 rad ptratio
# 2 3 rad ptratio
# 2 4 rad ptratio
# 2 5 rad ptratio
# 3 1 rad ptratio
# 3 2 rad ptratio
# 3 3 rad ptratio
# 3 4 rad ptratio
# 3 5 rad ptratio
# 4 1 rad ptratio
# 4 2 rad ptratio
# 4 3 rad ptratio
# 4 4 rad ptratio
# 4 5 rad ptratio
# 5 1 rad ptratio
# 5 2 rad ptratio
# 5 3 rad ptratio
# 5 4 rad ptratio
# 5 5 rad ptratio
判断是否还存在缺失值
mice包中的complete()函数用于接受前面生成的模型，并生成一个完整数据。
mice_output <- complete(mice_mod)
anyNA(mice_output)
[1] FALSE
# 已经不存在缺失值了
head(mice_output)
# crim zn indus chas nox rm age dis rad tax ptratio b lstat
# 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 18.3 396.90 4.98
# 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14
# 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03
# 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94
# 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33
# 6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12 5.21
检验替换后的值和原本数据框中的值的是否一样，结果发现，预测的值80%都和未被NA替换的值一样，精度很高（高于均数、众数、中位数等等）。
predict <- mice_output[is.na(BostonHousing$rad),'rad']
actuals <- original_data$rad[is.na(BostonHousing$rad)]
mean(actuals!=predict)
# 0.2
- 2.2.2
Hmisc包
Hmisc包的impute()函数（较粗略，不推荐）
# 使用均值进行插补
im_mean <- impute(BostonHousing$ptratio,mean) #除了mean也可以使用median等，或直接使用特定值
head(im_mean)
# 1 2 3 4 5 6
# 18.49765* 17.80000 17.80000 18.70000 18.70000 18.70000
# 带星号表示是插补的值
#缺失值替换
BostonHousing$ptratio <- NULL
BostonHousing$im_mean <- im_mean
- 2.2.3
VIM包⚠️
使用本身就含有缺失值的airquality数据集进行操作
data('airquality')
head(airquality)
# Ozone Solar.R Wind Temp Month Day
# 1 41 190 7.4 67 5 1
# 2 36 118 8.0 72 5 2
# 3 12 149 12.6 74 5 3
# 4 18 313 11.5 62 5 4
# 5 NA NA 14.3 56 5 5
# 6 28 NA 14.9 66 5 6
# 缺失值主要集中在1-2列
采用md_pattern进行查看
md.pattern(airquality)
# Wind Temp Month Day Solar.R Ozone
# 111 1 1 1 1 1 1 0
# 35 1 1 1 1 1 0 1
# 5 1 1 1 1 0 1 1
# 2 1 1 1 1 0 0 2
# 0 0 0 0 7 37 44
#一共有44个缺失值
VIM的过人之处是可以使用aggr()函数对缺失值进行可视化
library(VIM)
aggr_plot <- aggr(airquality,col=c('red','green'),numbers=TRUE,sortVars=TRUE,
labels=names(airquality),cex.axis=0.7,gap=3)
# color定义缺失值和非缺失值的颜色
# numbers=TRUE表示显示确实值的比例
# sortVars=TRUE表示对变量按照缺失值的多少进行排序
# labels=names(airquality)定义图的标题
# cex.axis=0.7,gap=3定义坐标轴的大小和图之间的间距
# Variables sorted by number of missings:
# Variable Count
# Ozone 0.24183007
# Solar.R 0.04575163
# Wind 0.00000000
# Temp 0.00000000
# Month 0.00000000
# Day 0.00000000
# col是设置非缺失值和缺失值颜的参数。
# numbers=TRUE是显示缺失值和非缺失值的比例。
# sortVars=TRUE是根据缺失值的多少进行排序。
# labels是对生成的图进行标签化。
# cex.axis是定义坐标轴的大小。
绿色代表缺失，左图中左边两个变量有缺失，右边三个没有缺失。柱子的长短显示了缺失值的多少。右边的图反应缺失值的模式。Ozone单独缺失的比例是0.229，Solar.R单独缺失的比例是0.033，两个一起缺失的比例是0.013。
marginplot()是另一个对缺失值进行可视化的函数
marginplot(airquality[1:2]) #对存在缺失值的1-2列进行操作
蓝色空心点表示没有缺失，红色实心点表示缺失。紫色（37个）表示两者都缺失。图左边的盒形图，红色表示Ozone缺失时Solar.R的分布，蓝色表示不缺失时的分布，下面的盒形图红色表示Solar.R缺失时Ozone的分布，蓝色表示Solar.R不缺失时Ozone的分布。
VIM 除了进行可视化以外，还可以采取高级的方法如线性回归、KNN等来对缺失值进行插补。
如：使用线性回归法进行插补
# 使用R语言内置sleep数据集进行演示
sleepIm <- regressionImp(Sleep + Gest + Span + Dream + NonD ~ BodyWgt + BrainWgt,data = sleep)
# 使用BodyWgt和BrainWgt这两个不存在缺省的变量来对Sleep、Gest、Span、Dream和NonD这五个存在缺省的参数进行预测（线性回归的方式）和插补
插补之后
head(sleepIm)
# BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger Sleep_imp Gest_imp Span_imp Dream_imp NonD_imp
# 1 6654.000 5712.0 -11.732867 -0.6897314 3.3 38.60000 645 3 5 3 FALSE FALSE FALSE TRUE TRUE
# 2 1.000 6.6 6.300000 2.0000000 8.3 4.50000 42 3 1 3 FALSE FALSE FALSE FALSE FALSE
# 3 3.385 44.5 8.987353 2.0132372 12.5 14.00000 60 1 1 1 FALSE FALSE FALSE TRUE TRUE
# 4 0.920 5.7 9.017324 2.0148478 16.5 15.50179 25 5 2 3 FALSE FALSE TRUE TRUE TRUE
# 5 2547.000 4603.0 2.100000 1.8000000 3.9 69.00000 624 3 5 4 FALSE FALSE FALSE FALSE FALSE
# 6 10.550 179.5 9.100000 0.7000000 9.8 27.00000 180 4 4 4 FALSE FALSE FALSE FALSE FALSE
#上面一部分显示原有缺失值已被插补，下面一部分是判断原有的值是否为缺省值，FALSE表示不是缺省值，TRUE表示是缺省值，但已经被填补。
regressionImp函数中还有个参数family，选'auto'。因为因变量包括连续型和离散型等，他们采取的回归方式是不一样的。auto模式可以自动采取合适的方法。
二、异常值
生成一个数据集
set.seed(2017)
mmhg <- sample(60:250,1000,replace = TRUE)
range(mmhg)
# [1] 60 250
min(mmhg)
# [1] 60
max(mmhg)
# [1] 250
quantile(mmhg)
# 0% 25% 50% 75% 100%
# 60.0 106.0 153.5 198.0 250.0
fivenum(mmhg)
# [1] 60.0 106.0 153.5 198.0 250.0
处理异常值的自定义函数
在使用时，下述代码唯一需要修改的就是var_name>230这里。在该代码中大于230被定义为极端值，可根据实际情况进行修改。
outlierHH <- function(dt,var){
var_name <- eval(substitute(var),eval(dt))
tot <- sum(!is.na(var_name))
na1 <- sum(is.na(var_name))
m1 <- mean(var_name,na.rm=T)
par(mfrow=c(2,2),oma=c(0,0,3,0))
boxplot(var_name,main='With outliers')
hist(var_name,main='With outliers',xlab=NA,ylab=NA)
outlier <- var_name[var_name>230]
mo <- mean(outlier)
var_name <- ifelse(var_name %in% outlier,NA,var_name)
boxplot(var_name,main='Without outliers')
hist(var_name,main='Without outliers',xlab=NA,ylab=NA)
title('Outlier Check',outer = TRUE)
na2 <- sum(is.na(var_name))
cat('Outliers identified:',na2-na1,'\n')
cat('Proportion (%) of outliers:',round((na2-na1)/tot*100,1),'\n')
cat('Mean of the outliers:',round(mo,2),'\n')
m2 <- mean(var_name,na.rm=T)
cat('Mean without removing outliers:',round(m1,2),'\n')
cat('Mean if we remove outliers:',round(m2,2),'\n')
response <- readline(prompt = 'Do you want to remove outliers
to replace with NA?[yes/no]:')
if(response=='y'|response=='yes'){
dt[as.character(substitute(var))] <- invisible(var_name)
assign(as.character(as.list(match.call())$dt),dt,envir = .GlobalEnv)
cat('Outliers successfully removed','\n')
return(invisible(dt))
}else{
cat('Nothing changed','\n')
return(invisible(var_name))
}
}
生成一个数据框
运行上述代码
df <- data.frame(bp=c(sample(80:250,1000,replace = T),NA,390,100))
outlierHH(df,bp)
返回了图并询问是否将大于230的值用NA替换
上：未去除异常值 下：去除异常值
三、重复值
1. 针对向量
- 1.1unique函数（常用）
x <- c(1,2,3,4,5,1,2,3)
unique(x)
# [1] 1 2 3 4 5
- 1.2 duplicated函数（常用）
duplicated(x)
# [1] FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE
x[duplicated(x)]
# [1] 1 2 3
x[!duplicated(x)]
# [1] 1 2 3 4 5
2. 针对数据框
paste函数
使用paste函数，把用来判断是否存在重复值的变量贴在一起。对paste之后生成的test数据框进行判断，看它是否存在重复值就可以了。
网友评论