單因素方差分析-aov()
單因素方差分析中，我們感興趣的是比較分類因子定義的兩個(gè)或多個(gè)組別中的因變量均值
以multcomp包中的cholesterol數(shù)據(jù)集為例，50個(gè)患者均接受降低膽固醇藥物治療(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
結(jié)果顯示五種療法效果極顯著不同（F<0.001）
gplots()包中的plotmeans()函數(shù)可以用來(lái)繪制帶有置信區(qū)間的組均值圖形
> library(gplots)
載入程輯包：‘gplots’
The following object is masked from ‘package:stats’:
lowess
> plotmeans(response~trt,data = cholesterol)
image.png
如圖可以看清楚他們之間的差異
多重比較
雖然單因素方差分析對(duì)各療法的F檢驗(yàn)表明五種藥物療法效果不同，但是并沒有告訴你哪種療法與其他療法不同，因此我們需要進(jìn)行多重比較
TukeyHSD的成對(duì)組間比較
> 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)
上述結(jié)果也可以用圖形·表示
>par(las=2)
> par(mar=c(5,8,4,2))
> plot(TukeyHSD(fit))
image.png
圖中置信區(qū)間包含0的療法差異不顯著(p>0.5)
glht()函數(shù)
mulicomp包中的glht()函數(shù)提供了更為全面的方法
> library(multcomp)
> par(mar=c(5,4,6,2))
> tuk<-glht(fit,linfct=mcp(trt="Tukey"))
image.png
如圖所示，含有相同字母的說(shuō)明均值差異不明顯
