统计学第八章 秩转换的非参数检验
秩转换的非参数检验
知识清单
- 配对样本比较的Wilcoxon符号秩检验（Wilcoxon signed-rank test）
- 两个独立样本比较的Wilcoxon之和检验（坑：也称为成组设计两样本）
- 完全随机设计多个样本的Kruskal-Wallis H检验
- 随机区组设计多个样本的Friedman M检验
1. 配对样本比较的Wilcoxon符号秩检验
1.1 导入数据
12份血清分别使用原方法和新方法测定谷-丙转氨酶，问结果有无差别。
数据文件：例08-01.sav
# 导入spss数据
paired_data <- haven::read_sav(
"E:/医学统计学（第4版）/各章例题SPSS数据文件/例08-01.sav")
colnames(paired_data) <- c("原法","新法")
# 查看数据
head(paired_data)
## # A tibble: 6 x 2
## 原法 新法
## <dbl> <dbl>
## 1 60 76
## 2 142 152
## 3 195 243
## 4 80 82
## 5 242 240
## 6 220 220
# 计算差值创建新的一列d
paired_data$d <- paired_data$新法-paired_data$原法
1.2 对差值d进行正态性检验
直观感受qq图和pp图
# qq图
a <-
qqnorm(paired_data$d, pch=16)
# qqline(paired_data$d, col="red", lwd=2)
abline(mean(paired_data$d), sd(paired_data$d),
col=2, lwd=2)
# pp图
plot((rank(paired_data$d)-0.5)/length(paired_data$d),
pnorm(mean=mean(paired_data$d), sd=sd(paired_data$d), paired_data$d),
main="Normal P-P plot", xlab="Theoretical probability",
ylab="Sample probability", pch=20)
abline(0, 1, col="red", lwd=2)
看qq图和pp图就知道这几个点明显不在一条直线上，不过还是需要进行统计推断，用qq图和pp图对正态性进行主观上的推断只有在较为显著的情况下才可以很明显肯定或否定其正态性。
偏度和峰度的假设检验
- Shapiro-Wilk单值检验法
shapiro.test(paired_data$d)
##
## Shapiro-Wilk normality test
##
## data: paired_data$d
## W = 0.87507, p-value = 0.07582
结果p值为0.07582，按检验水准为0.10，拒接H0，有统计学意义，可以认为这个样本总体不服从正态分布。
- 其他检验方法
# 使用normtest包分别检验偏度和峰度
normtest::skewness.norm.test(paired_data$d)
##
## Skewness test for normality
##
## data: paired_data$d
## T = -0.4838, p-value = 0.364
normtest::kurtosis.norm.test(paired_data$d)
##
## Kurtosis test for normality
##
## data: paired_data$d
## T = 4.1142, p-value = 0.2225
# 按教材P40页的算法检验偏度和峰度
kurt <- agricolae::kurtosis(paired_data$d)
skew <- agricolae::skewness(paired_data$d)
length_test <- length(paired_data$d)
ses <- sqrt(6*length_test*(length_test-1)/
((length_test-1)*(length_test+1)*(length_test+3)))
sek <- sqrt(24*length_test*(length_test-1)^2/
((length_test-3)*(length_test-2)*(length_test+3)*(length_test+5)))
totests <- skew/ses
totestk <- kurt/sek
pvals <- pt(totests, (length(paired_data$d)-1))
pvalk <- pt(totestk, (length(paired_data$d)-1))
cat("偏度系数p值：", pvals, "\n", "峰度系数P值：", pvalk, sep="")
## 偏度系数p值：0.1899681
## 峰度系数P值：0.9664809
结果峰度和偏度p值为都大于0.10，按检验水准为0.10，可能这些算法不太准确，还是以最经典的Shapiro-Wilk单值检验法为准。
1.3 进行配对样本的秩和检验
wilcox.test(paired_data$原法, paired_data$新法, paired=TRUE)
## Warning in wilcox.test.default(paired_data$原法, paired_data$新法, paired =
## TRUE): cannot compute exact p-value with ties
## Warning in wilcox.test.default(paired_data$原法, paired_data$新法, paired =
## TRUE): cannot compute exact p-value with zeroes
##
## Wilcoxon signed rank test with continuity correction
##
## data: paired_data$原法 and paired_data$新法
## V = 11.5, p-value = 0.06175
## alternative hypothesis: true location shift is not equal to 0
结果出现两条Warning：
第一条提示有秩相同的情况不能算出准确的p值。
第二条提示有差值为0的情况（配对的两个数值相等），不能算出准确的p值。
一般在秩相同的个数和差值为0的情况不是很多或者大样本数据时，可以直接忽略提示或者使用下面的语句：
wilcox.test(paired_data$原法, paired_data$新法, paired=TRUE,
exact=F)
##
## Wilcoxon signed rank test with continuity correction
##
## data: paired_data$原法 and paired_data$新法
## V = 11.5, p-value = 0.06175
## alternative hypothesis: true location shift is not equal to 0
结果计算p值为0.06175，按双侧检验alpha=0.05水准，不拒绝H0，还不能认为两法有区别
参考：
https://stat.ethz.ch/pipermail/r-help/2009-December/415767.html
http://www.r-tutor.com/elementary-statistics/non-parametric-methods/wilcoxon-signed-rank-test
n>50时，可近似u检验
2. 两个独立样本比较的Wilcoxon之和检验（坑：也称为成组设计两样本）
- 原始数据
与配对两样本的Wilcoxon实现相同，就是把paired参数设置为FALSE即可，这和t.test函数也是一样的。
n1>10或n2-n1>10时，可近似u检验
- 等级资料
确定各等级的合计人数，秩范围和平均秩，再计算秩和，其余同两独立样本Wilcoxon检验
3. 完全随机设计多个样本的Kruskal-Wallis H检验
- 原始数据
3种药物杀灭钉螺，其死亡率有无差别。
数据文件：例08-05.sav
# 导入spss数据
multi_samples_data <- haven::read_sav(
"E:/医学统计学（第4版）/各章例题SPSS数据文件/例08-05.sav")
colnames(multi_samples_data ) <- c("group","ratio")
# 查看数据
head(multi_samples_data)
## # A tibble: 6 x 2
## group ratio
## <dbl+lbl> <dbl>
## 1 1 32.5
## 2 1 35.5
## 3 1 40.5
## 4 1 46.0
## 5 1 49.0
## 6 2 16.0
因为是数据是百分率数据，所以肯定不服从正态分布，不能使用方差分析，直接使用H检验
kruskal.test(ratio~group, data=multi_samples_data)
##
## Kruskal-Wallis rank sum test
##
## data: ratio by group
## Kruskal-Wallis chi-squared = 9.74, df = 2, p-value = 0.007673
结果中的Kruskal-Wallis chi-squared也就是我们书上说的H值，p值为0.007673，按检验水准为0.05，拒绝H0，接受H1，可认为3种药物效果不同。
参考：
http://www.r-tutor.com/elementary-statistics/non-parametric-methods/kruskal-wallis-test
https://stat.ethz.ch/R-manual/R-devel/library/stats/html/kruskal.test.html
- 等级资料
确定各等级的合计人数，秩范围和平均秩，再计算秩和，其余同原始数据的H检验
4. 随机区组设计多个样本的Friedman M检验
8名受试者，对4种声音的反应有无差别
数据文件：例08-05.sav
4.1 导入数据
# 导入spss数据
multi_samples_block_data <- haven::read_sav(
"E:/医学统计学（第4版）/各章例题SPSS数据文件/例08-09.sav")
colnames(multi_samples_block_data ) <- c("sound1", "sound2", "sound3", "sound4")
# 查看数据
head(multi_samples_block_data)
## # A tibble: 6 x 4
## sound1 sound2 sound3 sound4
## <dbl> <dbl> <dbl> <dbl>
## 1 8.4 9.6 9.8 11.7
## 2 11.6 12.7 11.8 12.0
## 3 9.4 9.1 10.4 9.8
## 4 9.8 8.7 9.9 12.0
## 5 8.3 8.0 8.6 8.6
## 6 8.6 9.8 9.6 10.6
# 清理数据
multi_samples_block_data2 <- data.frame(blocks=factor(rep(seq(1, 8, by=1), 4)),
groups=factor(rep(c("A", "B", "C", "D"), each=8)),
value=unlist(multi_samples_block_data))
head(multi_samples_block_data2)
## blocks groups value
## sound11 1 A 8.4
## sound12 2 A 11.6
## sound13 3 A 9.4
## sound14 4 A 9.8
## sound15 5 A 8.3
## sound16 6 A 8.6
4.2 进行M检验
# 方法1：
friedman.test(as.matrix(multi_samples_block_data))
##
## Friedman rank sum test
##
## data: as.matrix(multi_samples_block_data)
## Friedman chi-squared = 15.152, df = 3, p-value = 0.001691
# 方法2：
friedman.test(value~groups|blocks, data=multi_samples_block_data2)
##
## Friedman rank sum test
##
## data: value and groups and blocks
## Friedman chi-squared = 15.152, df = 3, p-value = 0.001691
参考：
http://tutorial.math.trinity.edu/content/friedman-test-r
https://www.statmethods.net/stats/nonparametric.html
