之前介紹了多個(gè)樣本均數(shù)的多重比較,，今天說說kruskal-Wallis H檢驗(yàn)后的多重比較,，Friedman M檢驗(yàn)后的多重比較,。
也是和課本對照著來,，孫振球,，徐勇勇《醫(yī)學(xué)統(tǒng)計(jì)學(xué)》第四版。本書電子版已上傳到qq群中,，大家加群即可免費(fèi)獲取,。
非參數(shù)檢驗(yàn)后的多重比較，我們也是用這個(gè)寶藏R包：PMCMRplus,。
kruskal-Wallis H檢驗(yàn)及多重比較
使用課本例8-5的數(shù)據(jù),。
rm(list = ls())
death_rate <- c(32.5,35.5,40.5,46,49,16,20.5,22.5,29,36,6.5,9.0,12.5,18,24)
drug <- rep(c("Drug_A","drug_B","drug_C"),each=5)
mydata <- data.frame(death_rate,drug)
# 分類變量因子化
mydata$drug <- factor(mydata$drug)
str(mydata)
## 'data.frame': 15 obs. of 2 variables:
## $ death_rate: num 32.5 35.5 40.5 46 49 16 20.5 22.5 29 36 ...
## $ drug : Factor w/ 3 levels "Drug_A","drug_B",..: 1 1 1 1 1 2 2 2 2 2 ...
進(jìn)行kruskal-Wallis H 檢驗(yàn)：
fit <- kruskal.test(death_rate ~ drug, data = mydata)
fit
##
## Kruskal-Wallis rank sum test
##
## data: death_rate by drug
## Kruskal-Wallis chi-squared = 9.74, df = 2, p-value = 0.007673
多重檢驗(yàn)，課本上用的是Nemenyi檢驗(yàn),，我們通過多重比較的全能R包PMCMRplus實(shí)現(xiàn),。
library(PMCMRplus)
也是提供兩種輸入方式，直接提供kruskal-Wallis H檢驗(yàn)的結(jié)果,，或者formula形式,，都可以。
res <- kwAllPairsNemenyiTest(fit)
res <- kwAllPairsNemenyiTest(death_rate ~ drug, data = mydata)
summary(res)
##
## Pairwise comparisons using Tukey-Kramer-Nemenyi all-pairs test with Tukey-Dist approximation
## data: death_rate by drug
## alternative hypothesis: two.sided
## P value adjustment method: single-step
## H0
## q value Pr(>|q|)
## drug_B - Drug_A == 0 2.5 0.1805089
## drug_C - Drug_A == 0 4.4 0.0052932 **
## drug_C - drug_B == 0 1.9 0.3710425
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
除此之外,，還提供了kwAllPairsConoverTest()和kwAllPairsDunnTest()方法,。
Friedman M檢驗(yàn)及多重比較
使用課本本例8-9的數(shù)據(jù)，這個(gè)方式適用于隨機(jī)區(qū)組設(shè)計(jì)資料的多樣本比較,。
進(jìn)行Friedman M檢驗(yàn)需要矩陣形式的數(shù)據(jù)(這個(gè)是R語言里為數(shù)不多的不支持formula形式的統(tǒng)計(jì)檢驗(yàn)函數(shù)之一),，可以自己輸入，也可以直接讀取spss格式數(shù)據(jù),，然后變成矩陣即可,。
df <- matrix(
c(8.4, 11.6, 9.4, 9.8, 8.3, 8.6, 8.9, 7.8,
9.6, 12.7, 9.1, 8.7, 8, 9.8, 9, 8.2,
9.8, 11.8, 10.4, 9.9, 8.6, 9.6, 10.6, 8.5,
11.7, 12, 9.8, 12, 8.6, 10.6, 11.4, 10.8
),
byrow = F, nrow = 8,
dimnames = list(1:8,LETTERS[1:4])
)
print(df)
## A B C D
## 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
## 7 8.9 9.0 10.6 11.4
## 8 7.8 8.2 8.5 10.8
進(jìn)行Friedman M檢驗(yàn)：
fit <- friedman.test(df)
fit
##
## Friedman rank sum test
##
## data: df
## Friedman chi-squared = 15.152, df = 3, p-value = 0.001691
使用q檢驗(yàn)（quade test）進(jìn)行多重比較:
res <- quadeAllPairsTest(df)
summary(res)
##
## Pairwise comparisons using Quade's test with TDist approximation
## data: df
## P value adjustment method: holm
## H0
## t value Pr(>|t|)
## B - A == 0 1.126 0.27307148
## C - A == 0 3.236 0.01583866 *
## D - A == 0 5.093 0.00028885 ***
## C - B == 0 2.110 0.14099045
## D - B == 0 3.967 0.00351141 **
## D - C == 0 1.857 0.15475978
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
可以看到非常簡單，函數(shù)名稱清晰易懂,，結(jié)果也是非常直觀,，直接給出了兩兩比較的P值和統(tǒng)計(jì)量。
除此之外,，還提供了非常多其他方法:
大家選
