例子：以R语言multcomp包中的cholesterol数据集为例，50名患者均接收降低胆固醇药物治疗（trt）五种疗法中的一种疗法。其中三种治疗条件使用药物相同，分别是20mg一天一次（1time）、10mg一天两次（2times）和5mg一天四次（4times）。剩下的两种方式（drugD和drugE）代表候选药物。哪种药物疗法降低胆固醇最多呢？

>library(multcomp)

>attach(cholesterol)

>ge zu yang ben da xiao

>table(trt)

>#各组样本大小

>aggregate(response,by=list(trt),FUN=mean)

>#各组均值

>aggregate(response,by=list(trt),FUN=sd)

>#各组标准差

>fit<-aov(response~trt)

>#使用上一章方差分析表函数“anova.tab.R ”

>source('anova.tab.R')

>anova.tab(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

Total 49 1820.1

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#绘制各组均值及其置信区间的图形

>library(gplots)

>plotmeans(response~trt,xlab = "Treatment",ylab = "Response",main="Mean Plot\nwith 95% CI")

>detach(cholesterol)

和《用R语言做数据分析——单因素方差分析》中方差分析不同的是，我们使用了gplots包中的plotmeans()可以用来绘制带有置信区间的组均值图形，上述图形中展示了带有95%的置信区间的各疗法均值，可以清楚看到它们之间的差异。

多重比较

虽然anova()函数对各疗法的F检验表明了五种药物疗法效果不同，但是并没有告诉我们哪种疗法与其他疗法不同。多重比较可以解决这个问题，TukeyHSD()函数提供了各组均值差异的成对检验

> TukeyHSD(fit)

Tukey multiple comparisons of means

95% family-wise confidence level

Fit: aov(formula = response ~ trt)

$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

> par(las=2)

> par(mar=c(5,8,4,2))

> plot(TukeyHSD(fit))

从中可以得出结论：1time和2times的均值差异不显著（p=0.138），而1time和4times间的差异非常显著（p<0.001）。

multcomp包中的glht()函数提供了多重均值比较更为全面的方法，既适用于线性模型，也适用于广义线性模型，下面的代码重现了Tukey HSD，并用箱线图对结果进行展示：

> tuk<-glht(fit,linfct=mcp(trt="Tukey"))

> plot(cld(tuk,level=0.05),col="lightgrey")

其中cld()函数中的level选项设置了使用的显著水平（0.05，即95%的置信区间）。有相同字母的组说明均值差异不显著。可以看到：1time和2times差异不明显（有相同的字母a），2times和4times差异也不显著（有相同的字母b），而1time和4times差异显著（没有共同的字母）。

最后可得出结论：使用降低胆固醇的药物时，一天四次5mg剂量比一天一次20mg剂量效果更佳，也优于候选药物drugD，但药物drugE比其他药物和疗法都更优。

评估检验的假设条件

我们对于结果的信心依赖于做统计检验时数据满足假设条件的程度，单因素方差分析中，我们假设因变量服从正态分布，各组方差想的呢个，可以使用Q-Q图来检验正态性假设：

> library(car)

> qqPlot(lm(response~trt,data=cholesterol),simulate = TRUE,main = "Q-Q Plot", labels = FALSE)

数据落在95%的置信区间范围内，说明满足正态性假设。

R语言提供了一些可用来作方差齐性检验的函数，例如，可以通过如下代码来做Bartlett检验：

> bartlett.test(response~trt,data = cholesterol)

Bartlett test of homogeneity of

variances

data: response by trt

Bartlett's K-squared = 0.57975, df = 4,

p-value = 0.9653

Bartle检验表明五组的方差并没有显著不同（p=0.97）。因为方差齐性分析对离群点非常敏感，可利用car包中的outlierTest()函数来检测离群点：

> outlierTest(fit)

No Studentized residuals with Bonferonni p < 0.05

Largest |rstudent|:

rstudent unadjusted p-value Bonferonni p

19 2.251149 0.029422 NA

从输出结果来看，并没有证据说明数据中含有离群点（当p>1时将产生NA）。因此，根据Q-Q图、Bartlett检验和离群点检验，该数据可以用anova模型拟合达到很好的效果。