R语言学习-单因素方差分析与多重比较

单因素方差分析-aov()

单因素方差分析中，我们感兴趣的是比较分类因子定义的两个或多个组别中的因变量均值
以multcomp包中的cholesterol数据集为例，50个患者均接受降低胆固醇药物治疗(trt)五种疗法中的一种疗法。其中三种药物相同，分别是20mg一天一次（1time），10mg一天两次（2times）和5mg一天4次（4times），剩下两种方式代表候选药物

> library(multcomp)
> table(cholesterol$trt)

 1time 2times 4times  drugD  drugE 
    10     10     10     10     10 

可以看出接受每种方式的患者为10人

> aggregate(cholesterol$response,by=list(cholesterol$trt),FUN=mean)
  Group.1        x
1   1time  5.78197
2  2times  9.22497
3  4times 12.37478
4   drugD 15.36117
5   drugE 20.94752

各组均值显示drugE的降低的胆固醇最多，而1time最少

> fit<-aov(response~trt,data = cholesterol) #aov(formula,data=dataframe)
> summary(fit)
            Df Sum Sq Mean Sq F value   Pr(>F)    
trt          4 1351.4   337.8   32.43 9.82e-13 ***
Residuals   45  468.8    10.4                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

结果显示五种疗法效果极显著不同（F<0.001）
gplots()包中的plotmeans()函数可以用来绘制带有置信区间的组均值图形

> library(gplots)

载入程辑包：‘gplots’

The following object is masked from ‘package:stats’:

    lowess

> plotmeans(response~trt,data = cholesterol)

image.png

如图可以看清楚他们之间的差异

多重比较

虽然单因素方差分析对各疗法的F检验表明五种药物疗法效果不同，但是并没有告诉你哪种疗法与其他疗法不同，因此我们需要进行多重比较

TukeyHSD的成对组间比较

> TukeyHSD(fit)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = response ~ trt, data = cholesterol)

$trt
                  diff        lwr       upr     p adj
2times-1time   3.44300 -0.6582817  7.544282 0.1380949
4times-1time   6.59281  2.4915283 10.694092 0.0003542
drugD-1time    9.57920  5.4779183 13.680482 0.0000003
drugE-1time   15.16555 11.0642683 19.266832 0.0000000
4times-2times  3.14981 -0.9514717  7.251092 0.2050382
drugD-2times   6.13620  2.0349183 10.237482 0.0009611
drugE-2times  11.72255  7.6212683 15.823832 0.0000000
drugD-4times   2.98639 -1.1148917  7.087672 0.2512446
drugE-4times   8.57274  4.4714583 12.674022 0.0000037
drugE-drugD    5.58635  1.4850683  9.687632 0.0030633

可以看出2time-1time差异不显著（p=0.138），而drugE-1times差异极显著(p<0.001)
上述结果也可以用图形·表示

>par(las=2)
> par(mar=c(5,8,4,2))
> plot(TukeyHSD(fit))

image.png

图中置信区间包含0的疗法差异不显著(p>0.5)

glht()函数

mulicomp包中的glht()函数提供了更为全面的方法

> library(multcomp)
> par(mar=c(5,4,6,2))
> tuk<-glht(fit,linfct=mcp(trt="Tukey"))

image.png

如图所示，含有相同字母的说明均值差异不明显