R语言参数检验多重比较

本文首发于公众号：医学和生信笔记，完美观看体验请至公众号查看本文。

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。

对于多个样本均数的多重比较，比较常用的是LSD-t，SNK，Dunnett，Tukey等，这些方法在之前的推文中介绍过。

但是之前介绍的是用不同的R包完成的，整洁一致性不够，其实这些都是可以通过多重比较的全能R包：PMCMRplus完成的。

下面我们展示下~

还是使用课本例4-2的数据（孙振球，徐勇勇《医学统计学》第四版）。

trt<-c(rep("group1",30),rep("group2",30),rep("group3",30),rep("group4",30))

weight<-c(3.53,4.59,4.34,2.66,3.59,3.13,3.30,4.04,3.53,3.56,3.85,4.07,1.37,
          3.93,2.33,2.98,4.00,3.55,2.64,2.56,3.50,3.25,2.96,4.30,3.52,3.93,
          4.19,2.96,4.16,2.59,2.42,3.36,4.32,2.34,2.68,2.95,2.36,2.56,2.52,
          2.27,2.98,3.72,2.65,2.22,2.90,1.98,2.63,2.86,2.93,2.17,2.72,1.56,
          3.11,1.81,1.77,2.80,3.57,2.97,4.02,2.31,2.86,2.28,2.39,2.28,2.48,
          2.28,3.48,2.42,2.41,2.66,3.29,2.70,2.66,3.68,2.65,2.66,2.32,2.61,
          3.64,2.58,3.65,3.21,2.23,2.32,2.68,3.04,2.81,3.02,1.97,1.68,0.89,
          1.06,1.08,1.27,1.63,1.89,1.31,2.51,1.88,1.41,3.19,1.92,0.94,2.11,
          2.81,1.98,1.74,2.16,3.37,2.97,1.69,1.19,2.17,2.28,1.72,2.47,1.02,
          2.52,2.10,3.71)

data1<-data.frame(trt,weight)

# 分类变量因子化
data1$trt <- factor(data1$trt)

str(data1)
## 'data.frame':	120 obs. of  2 variables:
##  $ trt   : Factor w/ 4 levels "group1","group2",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ weight: num  3.53 4.59 4.34 2.66 3.59 3.13 3.3 4.04 3.53 3.56 ...

完全随机设计的多样本均数比较是用的one-way anova：

fit <- aov(weight ~ trt, data = data1)
summary(fit)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## trt           3  32.16  10.719   24.88 1.67e-12 ***
## Residuals   116  49.97   0.431                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

但是这样并不能判断到底是哪两个组之间有差别，所以需要进行两两比较（事后检验，多重比较）。

# 没安装的需要安装下这个包
library(PMCMRplus)

LSD

首先我们可以把方差分析的结果fit，直接作为输入：

res <- lsdTest(fit)
summary(res) # 结果非常直观
## 
## 	Pairwise comparisons using Least Significant Difference Test
## data: weight by trt
## alternative hypothesis: two.sided
## P value adjustment method: none
## H0
##                      t value   Pr(>|t|)    
## group2 - group1 == 0  -4.219 4.8872e-05 ***
## group3 - group1 == 0  -4.322 3.2889e-05 ***
## group4 - group1 == 0  -8.639 3.5772e-14 ***
## group3 - group2 == 0  -0.102    0.91871    
## group4 - group2 == 0  -4.420 2.2345e-05 ***
## group4 - group3 == 0  -4.318 3.3397e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

或者也可以使用formula形式：

res <- lsdTest(weight ~ trt, data = data1)
summary(res)
## 
## 	Pairwise comparisons using Least Significant Difference Test
## data: weight by trt
## alternative hypothesis: two.sided
## P value adjustment method: none
## H0
##                      t value   Pr(>|t|)    
## group2 - group1 == 0  -4.219 4.8872e-05 ***
## group3 - group1 == 0  -4.322 3.2889e-05 ***
## group4 - group1 == 0  -8.639 3.5772e-14 ***
## group3 - group2 == 0  -0.102    0.91871    
## group4 - group2 == 0  -4.420 2.2345e-05 ***
## group4 - group3 == 0  -4.318 3.3397e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

两种方法结果是一样的。

SNK

也是两种输入形式都可以。

res <- snkTest(fit)
res <- snkTest(weight ~ trt, data = data1)
summary(res)
## 
## 	Pairwise comparisons using SNK test
## data: weight by trt
## alternative hypothesis: two.sided
## P value adjustment method: step down
## H0
##                      q value   Pr(>|q|)    
## group2 - group1 == 0  -5.967 4.8872e-05 ***
## group3 - group1 == 0  -6.112 9.7010e-05 ***
## group4 - group1 == 0 -12.218 2.5524e-13 ***
## group3 - group2 == 0  -0.145    0.91871    
## group4 - group2 == 0  -6.251 6.6031e-05 ***
## group4 - group3 == 0  -6.106 3.3397e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tukey

也是两种输入形式都可以。

res <- tukeyTest(fit)
res <- tukeyTest(weight ~ trt, data = data1)
summary(res)
## 
## 	Pairwise comparisons using Tukey's test
## data: weight by trt
## alternative hypothesis: two.sided
## P value adjustment method: single-step
## H0
##                      q value   Pr(>|q|)    
## group2 - group1 == 0  -5.967 0.00028254 ***
## group3 - group1 == 0  -6.112 0.00019092 ***
## group4 - group1 == 0 -12.218 2.5524e-13 ***
## group3 - group2 == 0  -0.145 0.99961470    
## group4 - group2 == 0  -6.251 0.00013018 ***
## group4 - group3 == 0  -6.106 0.00019385 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dunnett

两种输入形式都可以。

res <- dunnettTest(fit)
res <- dunnettTest(weight ~ trt, data = data1)
summary(res)
## 
## 	Pairwise comparisons using Dunnett's-test for multiple	
##                          comparisons with one control
## data: weight by trt
## alternative hypothesis: two.sided
## P value adjustment method: single-step
## H0
##                      t value   Pr(>|t|)    
## group2 - group1 == 0  -4.219 0.00013734 ***
## group3 - group1 == 0  -4.322 8.9420e-05 ***
## group4 - group1 == 0  -8.639 3.3973e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

简单！

下次继续介绍非参数检验的多重比较，主要是kruskal-Wallis H检验后的多重比较，Friedman M检验后的多重比较。

本文首发于公众号：医学和生信笔记，完美观看体验请至公众号查看本文。

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。