方差分析主要通過(guò)F檢驗(yàn)來(lái)進(jìn)行效果評(píng)測(cè)蔬顾，若治療方案的F檢驗(yàn)顯著隅居，則說(shuō)明檢驗(yàn)樣本組間均值不同峡碉。
ANOVA模型擬合
從函數(shù)形式上看须揣，ANOVA和回歸方法都是廣義線性模型的特例盐股。因此回歸分析章節(jié)中提到的lm()函數(shù)也能分析ANOVA模型。不過(guò)耻卡，在這個(gè)章節(jié)中疯汁，我們基本使用aov()函數(shù)。最后卵酪，會(huì)提供了個(gè)lm()函數(shù)的例子幌蚊。
aov()函數(shù)
aov()函數(shù)的語(yǔ)法為aov(formula, data=dataframe)。下表列舉了表達(dá)式可以使用的特殊符號(hào)溃卡。
| 符號(hào) | 用途 |
|---|---|
| ~ | 分隔符號(hào)溢豆，左邊為響應(yīng)變量（因變量），右邊為解釋變量（自變量） |
| : | 表示預(yù)測(cè)變量的交互項(xiàng) |
| * | 表示所有可能交互項(xiàng)的簡(jiǎn)潔方式 |
| ^ | 表示交互項(xiàng)達(dá)到某個(gè)次數(shù) |
| . | 表示包含除因變量外的所有變量 |
下面是常見(jiàn)研究設(shè)計(jì)的表達(dá)式
| 設(shè)計(jì) | 表達(dá)式 |
|---|---|
| 單因素ANOVA | y ~ A |
| 含單個(gè)協(xié)變量的單因素ANOVA | y ~ x + A |
| 雙因素ANOVA | y ~ A * B |
| 含兩個(gè)協(xié)變量的雙因素ANOVA | y ~ x1 + x2 + A * B |
| 隨機(jī)化區(qū)組 | y ~ B + A (B是區(qū)組因子) |
| 單因素組內(nèi)ANOVA | y ~ A + Error(subject/A) |
| 含單個(gè)組內(nèi)因子(W)和單個(gè)組間因子的重復(fù)測(cè)量ANOVA | y ~ B * W + Error(Subject/W) |
表達(dá)式中各項(xiàng)的順序
當(dāng)
? 因子不止一個(gè)瘸羡，并且是非平衡設(shè)計(jì)漩仙；
? 存在協(xié)變量
兩者之一時(shí)，等式右邊的變量都與其他變量相關(guān)犹赖。此時(shí)队他，我們無(wú)法清晰地劃分它們對(duì)因變量的影響。
例如峻村，對(duì)于雙因素方差分析麸折，若不同處理方式中的觀測(cè)數(shù)不同，那么模型y ~ A*B與模型y ~ B*A的結(jié)果不同粘昨。
R默認(rèn)類型I（序貫型）方法計(jì)算ANOVA效應(yīng)（類型II和III分別為分層和邊界型垢啼，詳見(jiàn)R實(shí)戰(zhàn)（第2版）202頁(yè)）。R中的ANOVA表的結(jié)果將評(píng)價(jià)：
- A對(duì)y的影響
- 控制A時(shí)雾棺，B對(duì)y的影響
- 控制A和B的主效應(yīng)時(shí)膊夹，A與B的交互影響。
一般來(lái)說(shuō)捌浩，越基礎(chǔ)性的效應(yīng)需要放在表達(dá)式前面放刨。
car包的Anova()函數(shù)提供了三種類型方法，若想與其他軟件（如SAS SPSS）提供的結(jié)果保持一致尸饺，可以使用它进统，細(xì)節(jié)可參考helo(Anova, package="car")助币。
單因素方差分析
單因素方法分析中，你感興趣的是比較分類因子定義的兩個(gè)或多個(gè)組別中的因變量均值螟碎。以multcomp包中cholesterol數(shù)據(jù)集為例（包含50個(gè)患者接收5種降低膽固醇療法的一種眉菱，前三種是同樣的藥物不同的用法，后二者是候選藥物）掉分。哪種藥物療法降低膽固醇最多呢俭缓？
> library(mvtnorm)
> library(survival)
> library(TH.data)
> library(MASS)
> library(multcomp)
> attach(cholesterol)
> table(trt)
trt
1time 2times 4times drugD drugE
10 10 10 10 10
> aggregate(response, by=list(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
> aggregate(response, by=list(trt),FUN=sd)
Group.1 x
1 1time 2.878113
2 2times 3.483054
3 4times 2.923119
4 drugD 3.454636
5 drugE 3.345003
> fit <- aov(response ~ trt)
> 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
> library(gplots)
載入程輯包：‘gplots’
The following object is masked from ‘package:stats’:
lowess
> plotmeans(response ~ trt, xlab = "Treatment", ylab = "Response", main = "Mean Plot \nwith 95% CI")
> detach(cholesterol)
從結(jié)果可以看到，均值顯示drugE降低膽固醇最多酥郭，各組標(biāo)準(zhǔn)差相對(duì)恒定华坦。ANOVA對(duì)治療方式的F檢驗(yàn)非常顯著，說(shuō)明五種療法的效果不同不从。
多重比較
雖然ANOVA對(duì)各種療法的F檢驗(yàn)表明五種藥物的治療效果不同惜姐，但是沒(méi)有告訴你哪種療法與其他療法不同。多重比較可以解決這個(gè)問(wèn)題椿息。例如歹袁，TukeyHSD()函數(shù)提供了對(duì)各組均值差異的成對(duì)檢驗(yàn)。
> 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))
第一個(gè)par語(yǔ)句用來(lái)旋轉(zhuǎn)軸標(biāo)簽寝优，第二個(gè)用來(lái)增大左邊界的面積条舔，可使標(biāo)簽擺放更美觀。
成對(duì)比較圖形如下圖所示倡勇。
multcomp包中的glht()函數(shù)提供了多重均值比較更為全面的方法逞刷，既適用于線性模型，也適用于廣義線性模型妻熊。下面代碼重現(xiàn)了上述檢驗(yàn)結(jié)果，并用不同的圖形進(jìn)行展示仑最。
> library(multcomp)
> par(mar=c(5,4,6,2))
> tur <- glht(fit, linfct=mcp(trt="Tukey"))
> plot(cld(tur, level=0.05), col="lightgrey")
par語(yǔ)句增大了頂部邊界面積扔役，cld()函數(shù)中的level選項(xiàng)設(shè)置了使用的顯著水平。
有相同的字母的組說(shuō)明均值差異不顯著警医。
評(píng)估檢驗(yàn)的假設(shè)條件
可以使用Q-Q圖來(lái)檢驗(yàn)正態(tài)性
> library(car)
> qqPlot(lm(response ~ trt, data = cholesterol), simulate=T, main="Q-Q Plot", labels=FALSE)
注意qqPlot需要lm()擬合亿胸。
方差齊次性檢驗(yàn)：
例如，可以通過(guò)如下代碼做Bartlett檢驗(yàn)
> 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
結(jié)果表明五組方差沒(méi)有顯著地不同预皇。
注意侈玄，方差齊性分析對(duì)離群點(diǎn)非常敏感∫魑拢可以利用car包outlierTest()檢驗(yàn)序仙。
單因素協(xié)方差分析
ANCOVA擴(kuò)展了ANOVA，包含一個(gè)或多個(gè)定量的協(xié)變量鲁豪。
下面的例子來(lái)自multcomp包中的litter數(shù)據(jù)集潘悼。懷孕的小鼠被分為四個(gè)小組律秃，每組接受不同劑量的藥物處理。產(chǎn)下幼崽的體重均值為因變量治唤，懷孕時(shí)間為協(xié)變量棒动。
> data(litter, package = 'multcomp')
> attach(litter)
> table(dose)
dose
0 5 50 500
20 19 18 17
> aggregate(weight, by=list(dose), FUN=mean)
Group.1 x
1 0 32.30850
2 5 29.30842
3 50 29.86611
4 500 29.64647
> fit <- aov(weight ~ gesttime + dose)
> summary(fit)
Df Sum Sq Mean Sq F value Pr(>F)
gesttime 1 134.3 134.30 8.049 0.00597 **
dose 3 137.1 45.71 2.739 0.04988 *
Residuals 69 1151.3 16.69
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
由于使用了協(xié)變量，我們可以通過(guò)effect包中的effects()函數(shù)計(jì)算調(diào)整的組均值宾添。
> library(effects)
載入程輯包：‘effects’
The following object is masked from ‘package:car’:
Prestige
> effect("dose", fit)
dose effect
dose
0 5 50 500
32.35367 28.87672 30.56614 29.33460
同樣船惨，我們可以使用multcomp包對(duì)所有均值進(jìn)行成對(duì)比較。另外缕陕，該包還可以用來(lái)檢驗(yàn)用戶自定義的均值假設(shè)粱锐。
下面代碼清單可以用來(lái)檢驗(yàn)未用藥和其他三種藥條件影響是否不同。
> library(multcomp)
> contrast <- rbind("no drug vs. drug" = c(3, -1, -1, -1))
> summary(glht(fit, linfct=mcp(dose=contrast)))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: User-defined Contrasts
Fit: aov(formula = weight ~ gesttime + dose)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
no drug vs. drug == 0 8.284 3.209 2.581 0.012 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
對(duì)照c(3, -1, -1, -1)設(shè)定第一組和其他三組的均值進(jìn)行比較榄檬。假設(shè)檢驗(yàn)量t在p<0.05下顯著卜范。由此可以得出結(jié)論。詳見(jiàn)help(glht)鹿榜。
評(píng)估檢驗(yàn)的假設(shè)條件
ANCOVA與ANOVA相同海雪，都需要正態(tài)性和同方差假設(shè)，檢驗(yàn)可以參考上一節(jié)舱殿。另外ANCOVA還假定回歸斜率相同奥裸。
本例中，假定四個(gè)處理組通過(guò)懷孕時(shí)間來(lái)預(yù)測(cè)出生體重的回歸斜率都相同沪袭。ANCOVA模型包含懷孕時(shí)間X劑量的交互項(xiàng)時(shí)湾宙，可以對(duì)回歸斜率的同質(zhì)性進(jìn)行檢驗(yàn)。交互效果若顯著冈绊，則意味著時(shí)間和幼崽出生體重的關(guān)系依賴于藥物劑量的水平侠鳄。
fit2 <- aov(weight ~ gesttime*dose, data=litter)
summary(fit2)
HH包中的ancova()函數(shù)可以繪制因變量、協(xié)變量和因子之間的關(guān)系圖死宣。例如代碼：
library(HH)
ancova(weight ~ gesttime + dose, data=litter)
用回歸來(lái)做ANOVA
同樣是之前比較五種降低膽固醇藥物療法的影響的例子伟恶，我們分別用兩種不同的方法來(lái)做（aov()和lm()）。
> library(multcomp)
載入需要的程輯包：mvtnorm
載入需要的程輯包：survival
載入需要的程輯包：TH.data
載入需要的程輯包：MASS
載入程輯包：‘TH.data’
The following object is masked from ‘package:MASS’:
geyser
> levels(cholesterol)
NULL
> levels(cholesterol$trt)
[1] "1time" "2times" "4times" "drugD" "drugE"
> fit.aov <- aov(response ~ trt, data=cholesterol)
> summary(fit.aov)
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
> fit.lm <- lm(response ~ trt, data = cholesterol)
> summary(fit.lm) # 因子的第一個(gè)水平變成了參考組毅该，隨后的變量都以它為標(biāo)準(zhǔn)
Call:
lm(formula = response ~ trt, data = cholesterol)
Residuals:
Min 1Q Median 3Q Max
-6.5418 -1.9672 -0.0016 1.8901 6.6008
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.782 1.021 5.665 9.78e-07 ***
trt2times 3.443 1.443 2.385 0.0213 *
trt4times 6.593 1.443 4.568 3.82e-05 ***
trtdrugD 9.579 1.443 6.637 3.53e-08 ***
trtdrugE 15.166 1.443 10.507 1.08e-13 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.227 on 45 degrees of freedom
Multiple R-squared: 0.7425, Adjusted R-squared: 0.7196
F-statistic: 32.43 on 4 and 45 DF, p-value: 9.819e-13
