背景与问题
在数据分析中最常用的的参数方法,都需要先进行正态性、方差齐性检验。然而R中Shapiro-Wilk(夏皮罗)W正态性检验【命令shapiro.test()
】、Kolmogorov-Smirnov正态性检验【ks.test()
】分析的对象都只是对一个数值型向量可用,没有相应的S3 method,在面临多分组单变量、多分组多变量时,我等新手往往束手无策。具体说来,使用正态性检验的命令,不像是方差齐性/异质性检验(Bartlett Test of Homogeneity of Variances)【bartlett.test( )
】命令有S3 method、S3 method for class ‘formula’,可以直接针对数据框多分组单变量进行分析。
实现方法与途径
实现多分组单变量、多分组多变量分析的关键问题在于数据的操作方法,而不是检验方法本身。其实在R里面有不少命令可以实现分组操作,如by()
、"doBy"
程序包等,有相关命令可实现。
与下面使用的代码不同的是,by()
计算分析结果以提示信息形式呈现(而不是表格形式列出),"doBy"
程序包不是R内置的包,需要安装后加载使用。内置程序的命令中,数据分组/分割/切片可使用split()
命令进行,分割后的数据为'list'
格式;分析结果表提取可利用sapply()
命令直接获得,该命令对'list'
进行计算后输出结果为'data.frame'
格式,如果需要获得'list'
格式的,使用lapply()
命令。结合split()
与sapply()
即可进行数据分割及批量结果输出。
单变量多样本组分析例1
# 单变量多样本组分析
## Shapiro-Wilk(夏皮罗)W正态性检验
sapply(
split(iris$Sepal.Length, # 响应变量
iris$Species), # 分组变量
shapiro.test) # 分析使用的命令
单变量多样本组分析例2
# 单变量多样本组分析
## Kolmogorov-Smirnov正态性检验
sapply(
split(iris$Sepal.Length, # 响应变量
iris$Species), # 分组变量
function(d.f) ks.test(d.f, 'pnorm', mean(d.f), sd(d.f)) # 分析使用的命令
)
多变量多样本组分析例1
# 多变量多样本组分析
## Shapiro-Wilk(夏皮罗)W正态性检验
iris.a = list() # 生成空列表用于存储分析结果
data.clip <- split(iris[1:4], # 需要检验的响应变量在数据集中的哪些列
iris$Species) # 分组变量所在列
for(i in names(data.clip)) { # 使用for循环进行多变量分析
iris.a[[i]] =
sapply(data.clip[[i]],
shapiro.test) # 分析使用的命令
}
iris.a # 查看计算结果
多变量多样本组分析例2
# 多变量多样本组分析
## Kolmogorov-Smirnov正态性检验
iris.a = list() # 生成空列表用于存储分析结果
data.clip <- split(iris[1:4], # 需要检验的响应变量在数据集中的哪些列
iris$Species) # 分组变量所在列
for(i in names(data.clip)) { # 使用for循环进行多变量分析
iris.a[[i]] =
sapply(data.clip[[i]],
function(d.f) ks.test(d.f, 'pnorm', mean(d.f), sd(d.f))) # 分析使用的命令
}
iris.a # 查看计算结果
这两个例子的代码,可以拓展、移植到很多分析命令上,只需要改动分析使用的命令即可(见拓展部分)。对于单变量多样本组的分析,这个命令效率应该是最高的了。而对于多变量多样本组分析,因为其中涉及到使用for
循环,在数据量大的时候可能速度慢(没有测试过,一般认为R中循环计算的效率低下),目前没有想到更高效的方法,欢迎指教。
拓展
比如在进行描述统计时,统计指标较多,基于向量的计算分析自编分析命令脚本,脚本的编写上要简单很多,如下:
自编描述统计脚本
# 描述统计脚本
## 计算众数函数定义
majority <- function(x){
result = tryCatch({
as.numeric(names(table(x))[table(x) == max(table(x))])
}, warning = function(w) {
names(table(x))[table(x) == max(table(x))]
})
return(result)
}
# 描述统计脚本主体函数
data.describe <- function(x){
N = length(x) # 样本容量(数)
x.b = mean(x) # 均值
V = var(x) # 方差
S = sd(x) # 标准差
Me = median(x) # 中位数
CV = 100 * S/x.b # 变异系数
Max = max(x) # 最大值
Min = min(x) # 最小值
R = max(x) - min(x) # 极差
SE = S/sqrt(N) # 标准误
skew = N/((N-1)*(N-2)*S^3) * sum((x - x.b)^3) # 偏度系数
kurt = N*(N+1)/((N-1)*(N-2)*(N-3)*S^4) *
sum((x - x.b)^4) -
3*(N-1)^2/((N-2)*(N-3)) # 峰度系数
maj = paste(majority(x), collapse="/") # 众数(需计算众数脚本),可能有多个,需组合
data.frame( # 建立数据框输出结果
N = N, Mean = x.b, Var = V, Sd = S,
Median = Me, CV = CV, SE = SE,
Max = Max, Min = Min, Range = R,
Majority = maj, Skew = skew, Kurt = kurt,
row.names = c(substitute(x))) # 提取分析对象名称用于结果标识
}
基于自编脚本进行多分组多变量的描述统计
组合向量计算分析的脚本(data.describe()
)与数据操作命令(split()
与sapply()
)使用,同样以iris数据集为例。
单变量描述统计
# 单变量描述统计
data.describe(iris$Sepal.Length)
多变量描述统计
# 多变量描述统计
sapply(iris[1:4], data.describe) # 对多列数据进行描述统计分析
多组样本单变量描述统计
# 多组样本单变量描述统计
sapply(
split(iris$Sepal.Length, # 对此数据进行分割(切片、分组)
iris$Species), # 分割(切片、分组)的变量
data.describe # sapply计算使用的命令
)
多组样本多变量描述统计
# 多组样本多变量描述统计
iris.a = list() # 生成空列表
data.clip <- split(iris[1:4], iris$Species) # 对多列数据进行分割
for(i in names(data.clip)) {
iris.a[[i]] =
sapply(data.clip[[i]],
data.describe) # 对多列数据进行描述统计分析
}
iris.a # 查看计算结果
瞎说两句
对于非科班、没有系统学习的人来说,学习任何一门技能、知识本身要走过很多弯路。然而走过弯路进入坦途的人,无论是何原因,很少有人再回来引导跌跌撞撞的初学者,进入坦途后,多数忙着有意无意的建壁垒拉差距、圈地盘坐享其成。
在课程改革探索中,有幸接了两个学期的《生物统计附试验设计》课程,倒逼自己系统学习基础内容,真正深刻体会教学相长,今天看到以前写的博文,发现原来想要解决的问题,原本有很多简单的组合、方法、途径可以实现的东西,比如要进行多分组数据进行正态分布检验。反思过来,实际上是自己对数据操作的学习不够系统造成。