秩轉(zhuǎn)換的非參數(shù)檢驗(yàn)

知識(shí)清單

• 配對(duì)樣本比較的Wilcoxon符號(hào)秩檢驗(yàn)（Wilcoxon signed-rank test）
• 兩個(gè)獨(dú)立樣本比較的Wilcoxon之和檢驗(yàn)（坑：也稱為成組設(shè)計(jì)兩樣本）
• 完全隨機(jī)設(shè)計(jì)多個(gè)樣本的Kruskal-Wallis H檢驗(yàn)
• 隨機(jī)區(qū)組設(shè)計(jì)多個(gè)樣本的Friedman M檢驗(yàn)

1. 配對(duì)樣本比較的Wilcoxon符號(hào)秩檢驗(yàn)

1.1 導(dǎo)入數(shù)據(jù)

12份血清分別使用原方法和新方法測(cè)定谷-丙轉(zhuǎn)氨酶句占，問(wèn)結(jié)果有無(wú)差別丈氓。
數(shù)據(jù)文件：例08-01.sav

# 導(dǎo)入spss數(shù)據(jù)
paired_data <- haven::read_sav(
  "E:/醫(yī)學(xué)統(tǒng)計(jì)學(xué)（第4版）/各章例題SPSS數(shù)據(jù)文件/例08-01.sav")
colnames(paired_data) <- c("原法","新法")
# 查看數(shù)據(jù)
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

# 計(jì)算差值創(chuàng)建新的一列d
paired_data$d <- paired_data$新法-paired_data$原法

1.2 對(duì)差值d進(jìn)行正態(tài)性檢驗(yàn)

直觀感受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圖就知道這幾個(gè)點(diǎn)明顯不在一條直線上，不過(guò)還是需要進(jìn)行統(tǒng)計(jì)推斷，用qq圖和pp圖對(duì)正態(tài)性進(jìn)行主觀上的推斷只有在較為顯著的情況下才可以很明顯肯定或否定其正態(tài)性。

偏度和峰度的假設(shè)檢驗(yàn)

• Shapiro-Wilk單值檢驗(yàn)法

shapiro.test(paired_data$d)

## 
##  Shapiro-Wilk normality test
## 
## data:  paired_data$d
## W = 0.87507, p-value = 0.07582

結(jié)果p值為0.07582，按檢驗(yàn)水準(zhǔn)為0.10蝶俱，拒接H0，有統(tǒng)計(jì)學(xué)意義饥漫，可以認(rèn)為這個(gè)樣本總體不服從正態(tài)分布榨呆。

• 其他檢驗(yàn)方法

# 使用normtest包分別檢驗(yàn)偏度和峰度
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頁(yè)的算法檢驗(yàn)偏度和峰度
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("偏度系數(shù)p值：", pvals, "\n", "峰度系數(shù)P值：", pvalk, sep="")

## 偏度系數(shù)p值：0.1899681
## 峰度系數(shù)P值：0.9664809

結(jié)果峰度和偏度p值為都大于0.10，按檢驗(yàn)水準(zhǔn)為0.10趾浅，可能這些算法不太準(zhǔn)確愕提，還是以最經(jīng)典的Shapiro-Wilk單值檢驗(yàn)法為準(zhǔn)。

1.3 進(jìn)行配對(duì)樣本的秩和檢驗(yàn)

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

結(jié)果出現(xiàn)兩條Warning：
第一條提示有秩相同的情況不能算出準(zhǔn)確的p值皿哨。
第二條提示有差值為0的情況（配對(duì)的兩個(gè)數(shù)值相等）浅侨，不能算出準(zhǔn)確的p值。
一般在秩相同的個(gè)數(shù)和差值為0的情況不是很多或者大樣本數(shù)據(jù)時(shí)证膨，可以直接忽略提示或者使用下面的語(yǔ)句：

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

結(jié)果計(jì)算p值為0.06175如输，按雙側(cè)檢驗(yàn)alpha=0.05水準(zhǔn)，不拒絕H0央勒，還不能認(rèn)為兩法有區(qū)別

參考：
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時(shí)不见，可近似u檢驗(yàn)

2. 兩個(gè)獨(dú)立樣本比較的Wilcoxon之和檢驗(yàn)（坑：也稱為成組設(shè)計(jì)兩樣本）

• 原始數(shù)據(jù)

與配對(duì)兩樣本的Wilcoxon實(shí)現(xiàn)相同，就是把paired參數(shù)設(shè)置為FALSE即可崔步，這和t.test函數(shù)也是一樣的稳吮。

n1>10或n2-n1>10時(shí)，可近似u檢驗(yàn)

• 等級(jí)資料

確定各等級(jí)的合計(jì)人數(shù)井濒，秩范圍和平均秩灶似，再計(jì)算秩和列林，其余同兩獨(dú)立樣本W(wǎng)ilcoxon檢驗(yàn)

3. 完全隨機(jī)設(shè)計(jì)多個(gè)樣本的Kruskal-Wallis H檢驗(yàn)

• 原始數(shù)據(jù)

3種藥物殺滅釘螺，其死亡率有無(wú)差別酪惭。
數(shù)據(jù)文件：例08-05.sav

# 導(dǎo)入spss數(shù)據(jù)
multi_samples_data <- haven::read_sav(
  "E:/醫(yī)學(xué)統(tǒng)計(jì)學(xué)（第4版）/各章例題SPSS數(shù)據(jù)文件/例08-05.sav")
colnames(multi_samples_data ) <- c("group","ratio")
# 查看數(shù)據(jù)
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

因?yàn)槭菙?shù)據(jù)是百分率數(shù)據(jù)希痴，所以肯定不服從正態(tài)分布，不能使用方差分析春感，直接使用H檢驗(yàn)

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

結(jié)果中的Kruskal-Wallis chi-squared也就是我們書上說(shuō)的H值砌创，p值為0.007673，按檢驗(yàn)水準(zhǔn)為0.05鲫懒，拒絕H0嫩实，接受H1，可認(rèn)為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

• 等級(jí)資料

確定各等級(jí)的合計(jì)人數(shù)舶赔，秩范圍和平均秩，再計(jì)算秩和谦秧，其余同原始數(shù)據(jù)的H檢驗(yàn)

4. 隨機(jī)區(qū)組設(shè)計(jì)多個(gè)樣本的Friedman M檢驗(yàn)

8名受試者竟纳，對(duì)4種聲音的反應(yīng)有無(wú)差別
數(shù)據(jù)文件：例08-05.sav

4.1 導(dǎo)入數(shù)據(jù)

# 導(dǎo)入spss數(shù)據(jù)
multi_samples_block_data <- haven::read_sav(
  "E:/醫(yī)學(xué)統(tǒng)計(jì)學(xué)（第4版）/各章例題SPSS數(shù)據(jù)文件/例08-09.sav")
colnames(multi_samples_block_data ) <- c("sound1", "sound2", "sound3", "sound4")
# 查看數(shù)據(jù)
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

# 清理數(shù)據(jù)
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 進(jìn)行M檢驗(yàn)

# 方法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