背景与问题
在数据分析中最常用的的参数方法，都需要先进行正态性、方差齐性检验。然而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 # 查看计算结果
瞎说两句
对于非科班、没有系统学习的人来说，学习任何一门技能、知识本身要走过很多弯路。然而走过弯路进入坦途的人，无论是何原因，很少有人再回来引导跌跌撞撞的初学者，进入坦途后，多数忙着有意无意的建壁垒拉差距、圈地盘坐享其成。 在课程改革探索中，有幸接了两个学期的《生物统计附试验设计》课程，倒逼自己系统学习基础内容，真正深刻体会教学相长，今天看到以前写的博文，发现原来想要解决的问题，原本有很多简单的组合、方法、途径可以实现的东西，比如要进行多分组数据进行正态分布检验。反思过来，实际上是自己对数据操作的学习不够系统造成。
好文链接
您阅读本篇文章共花了：
发表评论