数据异常值，肿么办，要不要删除呢？｜转载

dt=carData::Davis$weight
set.seet <- 123
x <- data.frame(id=seq(1,15,1),value=ceiling(c(rnorm(10,5,1),1,2,20,30,40)))

############### 任意分布 ##############
summary(dt) #先做简单描述
hist(dt,breaks=20) #设置20个条

#------一、图形法--------
#-----1 2D3D散点图-----
dt=carData::Davis
scatter.smooth(dt$weight,dt$height)
#3个变量 3d散点图
scatterplot3d::scatterplot3d(dt$weight,dt$height,dt$repwt)
library(rgl)
plot3d(dt$weight,dt$height,dt$repwt,col=rainbow(20),size=10) #注意rgl包不支持中文变量名

#-----2 boxplot-----
boxplot(dt$weight~dt$sex) #绘制根据性别分类的箱式图

#-----3 残差图-----
#标准化残差的绝对值大于2的点就是离群点
dt2=as.data.frame(datasets::state.x77) #调用state.x77数据集
#线性回归
mod=lm(Murder~Population+Illiteracy+Income+Frost,data = dt2) #构建回归方程
par(mfrow=c(2,2)) #把4幅图整合在一个画面中
plot(mod) #绘图
#----3.1 QQ图带区间-----
car::qqPlot(mod,labels=row.names(dt2),simulate=T,main="Q-Q图")
#-----4 强影响点-----
car::outlierTest(mod)
#-----5 强杠杆hat value-----
#该点的hat value大于其均值的2~3倍,就被认为是异常点
hat.plot <- function(mod){
  p <- length(coefficients(mod))
  n <- length(fitted(mod))
  plot(hatvalues(mod))
  abline(h = c(2, 3) * p/n, col = "red", lty = 2)
  identify(1:n, hatvalues(mod), names(hatvalues(mod)))
}
hat.plot(mod)
#-----6 CooK距离-----
#Cook's D距离的值大于4/(n-k-1)为异常
#需要手动再确认？
cooksd <- cooks.distance(mod)
plot(cooksd, pch="*", cex=2) #绘制Cook距离
abline(h = 4*mean(cooksd, na.rm=T), col="red")  #添加参考线
text(x=1:length(cooksd)+1, y=cooksd, labels=ifelse(cooksd>4*mean(cooksd, na.rm=T),names(cooksd),""), col="red") #添加标签
#-----7 汇总图-----
#influencePlo函数把离群点,高杠杆值点,强影响点这些异常值都整合在一个图
#纵坐标超过+2或小于-2的点可被认为是离群点,水平轴超过0.2或0.3的点是高杠杆点,圆圈大小与影响成比例,圆圈过大的点是对模型参数的估计造成的不成比例影响的强影响点
car::influencePlot(mod)

#-----8 KM聚类-----
dt=carData::Davis
dt=dt[2:5] #还是使用那个身高体重的数据,但是先去掉性别那一列
dt=na.omit(dt) #去掉缺失值
rt=kmeans(dt, centers=3) #进行聚类
ct=rt$centers[rt$cluster, ]
jl=sqrt(rowSums((dt-ct)^2)) #计算数据与簇中心的距离
outliers=order(jl, decreasing = T)[1:5] #挑选出前5个最大距离
print(dt[outliers,]) #输出异常值

#------二、分布法--------
#-----1 3σ标准差-----
mean(dt)-3*sd(dt) #计算下限
mean(dt)+3*sd(dt) #计算上限
#找出>上限和<下限的行号
dt[which(dt<mean(dt)-3*sd(dt) )]
dt[which(dt>mean(dt)+3*sd(dt) )]

#-----2 推荐！boxplot.stats 偏态分布-----
boxplot(dt) #箱式图
bout <- boxplot.stats(dt)[["out"]]#输出的结果$out异常值,3个
#异常值的标准是：小于P25-1.5*IQR,或大于P75+1.5*IQR,
dt[which(dt$weight %in% bout),]
#dataframe案例
set.seet <- 123
x <- data.frame(id=seq(1,15,1),value=ceiling(c(rnorm(10,5,1),1,2,20,30,40)))
bout <- boxplot(x$value)$out
bout
x[which(x$value %in% bout),]

#-----3 Hampel滤波-----
median(dt)-3*mad(dt) #中位数减去3倍的绝对中位差
median(dt)+3*mad(dt) #中位数加上3倍的绝对中位差
dt[which(dt>median(dt)+3*mad(dt))] #取出大于上限的值,5个
dt[which(dt<median(dt)-3*mad(dt))] #取出小于下限的值

#-----5 LOF,不好看-----
#局部异常因子分析（Local Outlier Factor, LOF）
#一个点的局部密度明显的比它周围的点的密度小异常
pacman::p_load(DMwR2)
ot <- lofactor(dt, k=10) #k是计算局部异常因子所需要判断异常点周围的点的个数,ot就是每个数据的比值
plot(density(ot)) #绘制比值的核密度图,可画可不画
wz <- order(ot, decreasing = T)[1:3] #挑出比值最高的前五名数据作为异常值
dt[wz] #输出异常值

#----三、经验法-------
#-----3.1 盖帽法P1 P99替换-----
#盖帽法
ul=quantile(dt,0.99)
ll=quantile(dt,0.01)
dt[dt<ll]<-ll #小于P1的数都被P1替换
dt[dt>ul]<-ul #大于P99的数都被P99替换
#---3.2 盖帽pcap函数！批量替换连续变量P1 P99----
x$id <- factor(x$id)
pcap <- function(x){
  for (i in which(sapply(x, is.numeric))) {
    quantiles <- quantile( x[,i], c(.01, .99 ), na.rm =TRUE)
    x[,i] = ifelse(x[,i] < quantiles[1] , quantiles[1], x[,i])
    x[,i] = ifelse(x[,i] > quantiles[2] , quantiles[2], x[,i])}
  x}
pcap(x)

#---- 四、检验法(正态) ----
#-----4.1 Grubbs正态异常检验-----
pacman::p_load(outliers)
#outliers：：Grubbs检验
grubbs.test(dt)#P<0.05为异常

#-----4.2 EnvStats包,不好看-----
pacman::p_load(EnvStats)
rosnerTest(dt) #函数默认找3个数判断是不是离群值,也可以设置参数k=5改成测五个
rosnerTest(x$value,3,warn=F)