图形数据检验工具R_SPSS实战笔记(二)
数据分析领域初期需要特别注意,目前大多数的数据分析软件都要求数据的存储形式为"宽格式",即每一列都应当是一个变量,而每一行则代表一个单独的观测值。且需要“长格式”数据的时候,可以通过宽格式数据轻易进行转换;存储格式,推荐使用.text或.csv
另外,任何形式的数据检验(异常值识别[缺失/错误值],正态性检验,方差齐性检验,共线性问题)都具有一定的局限性,大多数情况下,最行之有效的方式便是使用图形工具。
注:本文部分内容参考自公众号"生态文献分享 石亚飞"
1.异常/错误/缺失值识别(箱线图/Cleveland点阵图/散点图矩阵)
1.1 异常描述与R实战
【缺失值识别】对于缺失值,其实很容易识别到,当数据量不多的时候,可以使用Excel的筛选或排序功能[注意一定要先定义表,然后再筛选排序,以免产生错序情况;tips:建议使用NA()代替,不建议直接删除];当数据量较多的时候,可以使用R中内置函数进行识别,详细代码见下(此外,缺失值的处理也是一个很头疼的问题,这里推荐使用R语言中的mice包进行缺失值插补,这里仅展示了利用mice包可视化缺失值)
【图例】
#载入相关R包
library(mice)
library(VIM)
#定义数据集
data <- data.frame(
Pitcher_height = c(744, 700, 714, 667, 600, 777, 640, 440, 715, 573, 1500, 650, 480, 545, 845, 560, 450, 600, 607, 675, NA, NA, 5.1, 534, 655, 65.5),
Mouth_diameter = c(34.3, 34.4, 28.9, 32.4, 29.1, 33.4, 34.5, 29.4, 39.5, 33.0, 33.8, 36.3, 27.0, 30.3, 37.3, 42.1, 31.2, 34.6, 33.5, 31.4, NA, NA, 0.3, 30.2, 35.8, 3.52),
Tube_diameter = c(18.6, 20.9, 19.7, 19.5, 17.5, 21.1, 18.6, 18.4, 19.7, 15.8, 19.1, 20.2, 18.1, 17.3, 19.3, 14.6, 20.6, 17.1, 14.8, 16.3, NA, NA, 0.1, 16.5, 15.7, 1.77)
)
#缺失值识别
md.pattern(data)
aggr_plot <- aggr(data,
col=c('navyblue','red'),
numbers=TRUE,
sortVars=TRUE,
labels=names(data),
cex.axis=.7,
gap=3,
ylab=c("Histogram of missing data","Pattern"))
【异常/错误值识别】错误值通常是以异常值的形式出现,但是异常值不一定是错误值。一般情况下,如果某些数据超过了整体数据值分布的上10%分位数和下95%分位数,即可以认为这些数据是潜在的异常值,而如果数据超过了上5%分位数和下95%分位数,则可以认为数据大概率是异常值。首先可以使用"箱线图/茎叶图"快速识别异常值,并进一步通过两两变量之间的"散点图矩阵"对错误值进行识别。
#异常值识别
df1 <- data.frame(
id = rep(1:3, each = 4),
variable1 = rnorm(12),
variable2 = rnorm(12),
variable3 = rnorm(12)
)
stem(df[,2],scale = 2) #绘制第2列变量的茎叶图,茎部分显示2位数
# 使用 pivot_longer() 函数转换为长格式
library(tidyr)
long_df <- pivot_longer(df1, cols = -id, names_to = "variable", values_to = "value")#在这里,cols = -id 表示除了 id 列之外的所有列都应该被转换,names_to 是新列的名称,用于存储变量名,values_to 是新列的名称,用于存储变量值
# 绘制横向箱线图
ggplot(long_df, aes(y = value, x = variable)) +
geom_boxplot() +
labs(y = "Value", x = "Variable", title = "Horizontal Boxplot") +
theme_minimal()
#---环境建设---#
rm(list = ls());#清空变量空间
library(datasets)
library(export)
#---数据读入---#
dataname <- "iris"
rawdata <- datasets::iris
#散点图矩阵——Iris(GGally_ggpairs)
library(GGally)
ggpairs(rawdata, columns=1:5, aes(color=Species)) +
ggtitle("散点图矩阵——Iris(GGally_ggpairs)")+
theme_bw()
graph2png(file = paste("散点图矩阵(GGally_ggpairs) of", dataname))
2.正态性(直方图,QQ图)
2.1通过直方图与正态性检验了解数据的分布特征
计量数据的直方图,在SPSS中有多种方法可以获得,比如在图2-1-3中,点击Frequencies对话框中的【Charts】按钮,就可以选择输出直方图:
也可以在进行正态性检验时,一并输出直方图,操作如下:
点击菜单:Analyze => Descriptive Statistics => Explore
设置变量及分组因素(本例中为sex,将分别输出男性与女性的正态性检验结果),点击【Plots】按钮,在Expore: Plots对话框中,勾选直方图和正态性检验两个选项(如上图),点击【Continue】=>【OK】,就能输出结果,因为Expore的默认选项我们没有去除,所以输出的结果中内容非常多,截取我们需要的信息如下:
(1)正态性检验的结果
SPSS中提供了两种正态性检验的方法:
KS(Kolmogorov-Smirnov)检验和W检验(Shapiro-Wilk),一般认为样本量的范围在4~2000时,W检验的检验效能较高,而样本量超过2000时应采用KS检验结果。
本例中总样本量为199,因此选择W检验的结果:
对于Female,P = 0.043 < 0.05,因此拒绝原假设,认为女性的BMI不服从正态分布;同理,认为男性的BMI也不服从正态分布。
(2)正态Q-Q图(Quantile-Quantile Plot)
正态Q-Q图是直观地检查数据是否服从正态分布的方法,如下图:
如果数据呈正态分布,则Q-Q图中的点应位于对角线上。相反,图中的点与对角线的偏差越大,说明数据服从正态分布的可能性就越小。
当然,Q-Q图的方法是一种直观的目视法,通过正态Q-Q图判断数据是否服从正态分布有些主观,但是可以结合正态性检验的结果,对数据的正态性做一个综合判断。
通过正态性检验与正态Q-Q图,我们可以判断:无论哪种性别,BMI都是不服从正态分布的。事实上,从下面的直方图可以看出,男性和女性的BMI都有一点右偏,但是右偏并不严重,所以我们看到的Q-Q图中,大多数的点都在直线上或附近,只有少数离群值(与其它值相比,异常小或异常大的值)脱离对角线较远。
(3)直方图
本例输出的直方图:
直方图是观察数据分布最直观的方法,上图显示出BMI的分布:女性的BMI对称性稍差,男性的因为右侧的离群值因而显得比女性的更右偏一些。
(4)根据数据的分布特征,选择适当的描述性统计量
本例中BMI呈右偏态分布,宜选择中位数来描述其集中趋势,选择Q1和Q3来描述其离散趋势;相反的,如果数据服从正态分布,就可以选择算术均数(𝑥ˉ_x_ˉ)和标准差(S)来描述这个样本数据。
注意!
如果样本量很大(如下图),正态性检验的方法,几乎总是拒绝原假设,得到数据不服从正态分布的推断。
所以在大样本量的情况下,最好通过直方图或正态Q-Q图直观判断数据分布。
而在样本量较小的情况下,Shapiro-Wilk正态性检验的检验效能往往较低,如:
在R4.1中模拟生成一个beta分布的数据
dd=rbeta(10000,2,5)
数据分布:
从中随机抽取50个数据,进行Shapiro-Wilk正态性检验,只有大约50%的样本,能拒绝𝐻0_H_0,得到数据不服从正态分布的推论,也就是检验效能约为50%。如果样本量降低到20,针对此数据的检验效能只有不到20%。
所以对于小样本计量资料,除非有比较充分的证据,或者图形方法显示出明显的正态分布的钟型,“分布未知”可能是对其数据分布最好的判断。
2.2R实战
#方法一
# load library rcompanion
library(rcompanion)
# create data vector
x= c(rnorm(10000, 2000, 45))
# draw plot using plotNormalHistogram() function
# use colour and size arguments for formatting plot
plotNormalHistogram( x, prob = FALSE, col="black", border="green",
main = "Normal Distribution overlay on Histogram",
length = 10000, linecol="red", lwd=3 )
#绘制QQ图
qqnorm(x)
#方法二
#利用ggplot直接绘制
library(ggplot2)
# Make this example reproducible
set.seed(0)
# Define data
data <- data.frame(x = rnorm(1000))
# Create histogram and overlay normal curve
ggplot(data, aes(x)) +
geom_histogram(aes(y = ..density..), fill = 'lightgray', col = 'black') +
stat_function(fun = dnorm, args = list(mean = mean(data$x), sd = sd(data$x)))
3.方差齐性(箱线图)
3.1Levene检验
✔依次单击“分析-描述统计-探索”选项。
✔进入Levene分析界面,在设置中,将体重设为因变量、饲料类型纳入因子列表。
✔单击探索设置中的“图”设置,在其“含莱文检验的分布-水平图”中选择“未转换”选项,即可生成方差分析的结果。
✔结果解读
生成结果如下表,P=0.087,说明,尚不能认为两个总体方差有差异,即方差齐性。
4.共线性(VIF,散点图,相关性,PCA)
library(FactoMineR)
library(factoextra)
library(ggplot2)
library(ggpubr)
#----------------------------------------------------------------------------------------------------------------
#*****第一步:数据标准化。在FactoMineR的PCA()函数中,默认设置对数据执行标准化,即scale.unit=TRUE。*****
#加载数据
data("iris")
#查看数据结构
str(iris)
#PCA只对数值矩阵生效,我们不需要第五列(物种)
res.pca <- PCA(iris[,-5])
#查看结果
summary(res.pca)
#使用factoextra包的get_pca_var()函数来提取变量的分析结果
var <- get_pca_var(res.pca)
#------------------------------------
#{结果如下:这个结果是一个列表,其中:
# coord表示用于创建散点图的变量坐标。coord实际上就是成分载荷,指观测变量与主成分的相关系数
# cor表示相关系数
# cos2表示因子质量,其计算公式为:var.cos2 = var.coord * var.coord
# contrib表示包含变量对主成分的贡献(百分比),其计算公式为: (var.cos2 * 100) / (total cos2 of the component)。
# 此外,我们可以使用factoextra包别的函数来有针对性地提取别的结果,基本函数有这些:
# get_eigenvalue(res.pca):提取主成分的特征值/方差
# fviz_eig(res.pca):可视化特征值
# get_pca_ind(res.pca),get_pca_var(res.pca):分别提取个体和变量的结果。
# fviz_pca_ind(res.pca),fviz_pca_var(res.pca):分别可视化结果个体和变量。
# fviz_pca_biplot(res.pca):制作主成分分析散点图biplot图。
#----------------------------------------------------------------------------------------------------------------
#****第二步:选取主成分****
#绘制主成分碎石图,查看每一个主成分能在多大程度上代表原来的特征
fviz_screeplot(res.pca, addlabels = TRUE)
#提取特征值
get_eigenvalue(res.pca)
#绘制主成分碎石图(以特征值为纵轴)
fviz_screeplot(res.pca, choice = "eigenvalue",addlabels = TRUE)
#-------------------------------------------------------------
#断棍模型(一是选取特征根大于所对应的断棍长度的轴;第二种是选取特征根的总和大于所对应断棍长度总和前几轴)
#提取特征根
res <- res.pca$eig[,1]
#断棍模型
n <- length(res)
bsm <- data.frame(j=seq(1:n), p=0)
bsm$p[1] <- 1/n
for (i in 2:n) {
bsm$p[i] = bsm$p[i-1] + (1/(n + 1 - i))
}
bsm$p <- 100*bsm$p/n
bsm
#可视化
barplot(t(cbind(100*res/sum(res),bsm$p[n:1])), beside=TRUE,
main="% 变差", col=c("bisque",2), las=2)
legend("topright", c("% 特征根", "断棍模型"),
pch=15, col=c("bisque",2), bty="n")
#--------------------------------------------------------------------------------------------------------------
#函数PCA()会返回我们2张图,即****PCA得分图/观测值坐标图和载荷图****。但这是作为探索性分析默认的显示方式;
#****第三步:PCA结果得分图****
#观测值坐标图对应了factoextra包中的fviz_pca_ind()函数,而变量分布图对应fviz_pca_var()函数
#****得分图绘制/观测值坐标图(I型标尺)****
fviz_pca_ind(res.pca,geom.ind = "point", # 只显示点而不显示文本,默认都显示
col.ind = iris$Species,# 设定颜色分类种类
palette = c("#00AFBB", "#E7B800", "#FC4E07"), # 设定颜色
addEllipses = TRUE,# 添加置信椭圆
legend.title = "Groups" #添加标题
)
#绘制全部样本的边界,多边形
fviz_pca_ind(res.pca,
mean.point = F, # 去除分组的中心点,否则每个群中间会有一个比较大的点
label = "none", # 隐藏每一个样本的标签
habillage = iris$Species, # 根据样本类型来着色
palette = c("purple","orange","blue"), # 三个组设置三种颜色
addEllipses = TRUE, # 添加边界线
ellipse.type = "convex" # 设置边界线为多边
)
#调整样本点大小
fviz_pca_ind(res.pca, pointsize = "cos2",
pointshape = 21, fill = "#E7B800",
repel = TRUE # 避免文本相互重叠
)
#-------------------------------------------------------------------------------------------------------------
#****第四步:PCA载荷图****
#*最简单的图只需要一行代码
fviz_pca_var(res.pca)
#展示变量对主成分的贡献
fviz_pca_var(res.pca,
col.var = "contrib", #根据贡献度着色
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07") #设置色带
)
#*-------------------------------------------------------------------------
#更直观的方式是绘制相关热图
# 获取样本的主成分分析结果
var <- get_pca_var(res.pca)
#加载绘制热图的包
library("corrplot")
#绘制热图
corrplot(var$contrib, is.corr = FALSE)
#-------------------------------------------------------------------------------------------------------------
#****第五步:绘制双标图****
#绘制一个简单的双标图
fviz_pca_biplot(res.pca, repel = TRUE,
col.var = "#2E9FDF", # 变量颜色
col.ind = "#696969" # 样本颜色
)
# 只保留变量的标签
fviz_pca_biplot(res.pca, label ="var")
#增加分组和置信区间
fviz_pca_biplot(res.pca,
col.ind = iris$Species,
addEllipses = TRUE,
label = "var",
col.var = "black",
repel = TRUE,
legend.title = "Species")
#进一步美化图形
fviz_pca_biplot(res.pca,
col.ind = iris$Species,
addEllipses = TRUE,
label = "var",
col.var = "black",
repel = TRUE,
legend.title = "Species")+
ggpubr::fill_palette("jco")+ #样本填充颜色,使用jco期刊的风格
ggpubr::color_palette("npg") #变量颜色,使用npg期刊的风格
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 【.NET】调用本地 Deepseek 模型
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库