《R語言實戰》第三部分第九章-方差分析複習筆記
1 術語
以焦慮症治療為例，現有兩種治療方案：認知行為療法（CBT）和眼動脫敏再加工法（EMDR）。我們招募10位焦慮症患者作為志願者，隨機分配一半的人接受為期五周的CBT，另外一半接受為期五周的EMDR，設計方案下表所示。在治療結束時，要求每位患者都填寫狀態特質焦慮問卷（STAI），也就是一份焦慮度測量的自我評測報告。
在這個實驗設計中，治療方案是兩水平（CBT、EMDR）的組間因子。之所以稱其為組間因子，是因為每位患者都僅被分配到一個組別中，沒有患者同時接受CBT和EMDR。表中字母s代表受試者（患者）。STAI是因變數，治療方案是自變數。由於在每種治療方案下觀測數相等，因此這種設計也稱為均衡設計（balanced design）；若觀測數不同，則稱作非均衡設計（unbalanceddesign）。
單因素方差分析（單因素方差分析 - MBA智庫百科）：指對單因素試驗結果進行分析，檢驗因素對試驗結果有無顯著性影響的方法，是兩個樣本平均數比較的引伸，它是用來檢驗多個平均數之間的差異，從而確定因素對試驗結果有無顯著性影響的一種統計方法。
單因素組內方差分析及重複測量方差分析：
在上述的實例中，假如只對CBT的效果感興趣，則需將10個患者都放在CBT組中，然後在治療五周和六個月後分別評價療效，設計方案如下表所示。
協方差分析（Analysis of covariance）：
假設招募患者時使用抑鬱症的自我評測報告，比如白氏抑鬱症量表（BDI），記錄了他們的抑鬱水平，那麼你可以在評測療法類型的影響前，對任何抑鬱水平的組間差異進行統計性調整。本案例中，BDI為協變數，該設計為協方差分析（ANCOVA）。協方差分析是關於如何調節協變數對因變數的影響效應，從而更加有效地分析實驗處理效應的一種統計技術，也是對實驗進行統計控制的一種綜合方差分析和回歸分析的方法。協方差是用來度量兩個變數之間「協同變異」大小的總體參數，即二個變數相互影響大小的參數，協方差的絕對值越大，二個變數相互影響越大。
多元方差分析（Multivariate analysis of variance）：
當因變數不止一個時，以上述實驗為例，為增強研究的有效性，可以對焦慮症進行其他的測量（比如家庭評分、醫師評分，以及焦慮症對日常行為的影響評價），此時稱作多元方差分析。
多元協方差分析（Multivariate analysis of covariance）：
在多元方差分析的基礎上，如果引入協變數，則稱作多元協方差分析。
以下兩個圖表，引自R語言實戰筆記--第九章 方差分析
能夠比較好的說明各種術語，向原作者表示致意。
2 ANOVA模型擬合
ANOVA模型擬合採用aov()函數，語法如下：
aov(formula, data = dataframe)n
其中formula為表達式，表達式在之前的章節中遇到過，這裡總結複習一下，R表達式中的符號意義如下：
此外，當因子不止一個時，以上表達式中因子出現的順序也非常重要。
3 單因素方差分析
單因素方差分析主要比較分類因子定義的不同組別的因變數均值。以multcomp包中的cholesterol數據集為例，研究哪種藥物療法降低膽固醇最多。
>install.packages("multcomp")n> library(multcomp)n載入需要的程輯包：mvtnormn載入需要的程輯包：survivaln載入需要的程輯包：TH.datan載入需要的程輯包：MASSnn載入程輯包：『TH.data』nnThe following object is masked from 『package:MASS』:nn geysernn> attach(cholesterol)n> head(cholesterol)n trt responsen1 1time 3.8612n2 1time 10.3868n3 1time 5.9059n4 1time 3.0609n5 1time 7.7204n6 1time 2.7139n> table(trt)ntrtn 1time 2times 4times drugD drugE n 10 10 10 10 10 n> aggregate(response,by=list(trt),FUN = mean)n Group.1 xn1 1time 5.78197n2 2times 9.22497n3 4times 12.37478n4 drugD 15.36117n5 drugE 20.94752n> aggregate(response,by=list(trt),FUN = sd)n Group.1 xn1 1time 2.878113n2 2times 3.483054n3 4times 2.923119n4 drugD 3.454636n5 drugE 3.345003n> fit_91 <- aov(response ~ trt)n> summary(fit_91)n Df Sum Sq Mean Sq F value Pr(>F) ntrt 4 1351.4 337.8 32.43 9.82e-13 ***nResiduals 45 468.8 10.4 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n> library(gplots)nError in library(gplots) : 不存在叫『gplots』這個名字的程輯包n> install.packages("gplots")n> library(gplots)nn載入程輯包：『gplots』nnThe following object is masked _by_ 『.GlobalEnv』:nn residplotnnThe following object is masked from 『package:stats』:nn lowessnn> plotmeans(response ~ trt, xlab = "Treatment", ylab = "Responese", main = "Mean Plotnwith 95% CI")n
gplots包中的plotmeans()函數可以用來繪製帶有置信區間的組均值圖形，默認置信區間為95%。
ANOVA僅僅表明五種藥物療法效果不同，但是並沒有告知哪種療法與其他療法不同。多重比較TukeyHSD()函數提供了對各組均值差異的成對檢驗。
> TukeyHSD(fit_91)n Tukey multiple comparisons of meansn 95% family-wise confidence levelnnFit: aov(formula = response ~ trt)nn$trtn diff lwr upr p adjn2times-1time 3.44300 -0.6582817 7.544282 0.1380949n4times-1time 6.59281 2.4915283 10.694092 0.0003542ndrugD-1time 9.57920 5.4779183 13.680482 0.0000003ndrugE-1time 15.16555 11.0642683 19.266832 0.0000000n4times-2times 3.14981 -0.9514717 7.251092 0.2050382ndrugD-2times 6.13620 2.0349183 10.237482 0.0009611ndrugE-2times 11.72255 7.6212683 15.823832 0.0000000ndrugD-4times 2.98639 -1.1148917 7.087672 0.2512446ndrugE-4times 8.57274 4.4714583 12.674022 0.0000037ndrugE-drugD 5.58635 1.4850683 9.687632 0.0030633nn> par(las=2)n> par(mar=c(5,8,4,2))n> plot(TukeyHSD(fit_91))n
multcomp包中的glht()函數提供了多重均值更為全面的方法。
> par(mar=c(5,4,6,2))n> tuk <- glht(fit_91, linfct = mcp(trt="Tukey"))n> par(las = 1)n> plot(cld(tuk, level = .05),col = "red")n
注意，記得par(las = 1)函數將軸標籤旋轉回來。
同樣的，我們對於結果的信心依賴於做統計檢驗時數據滿足假設條件的程度。單因素方差分析中，我們假設因變數服從正態分布，各組方差相等。使用Q-Q圖來檢驗正態性假設：
> library(car)n> qqPlot(lm(response ~ trt, data=cholesterol),nsimulate=TRUE, main="Q-Q Plot", labels=FALSE)n
qqPlot()要求用lm()擬合，從結果上看，數據落在95%的置信區間內，滿足正態性假設。此外，R還提供了一些可用來做方差齊性檢驗的函數，比如bartlett.test()函數（Bartlett檢驗）、fligner.test()函數（Fligner-Killeen檢驗）以及HH包中的hov()函數（Brown-Forsythe test ）。
考慮到，方差齊性分析對離群點非常敏感，可利用car包中的outlierTest(0函數檢測離群點。
> outlierTest(fit_91)nnNo Studentized residuals with Bonferonni p < 0.05nLargest |rstudent|:n rstudent unadjusted p-value Bonferonni pn19 2.251149 0.029422 NAn> n
4 單因素協方差分析
單因素協方差分析（ANCOVA）擴展了單因素方差分析（ANOVA），包含一個或多個定量的協變數。以multcomp包中的litter數據集為例，四組懷孕小鼠接受不同劑量（0、5、50或500）的藥物處理，懷孕時間為協變數，產下幼崽的體重均值為因變數。
> data(litter,package = "multcomp")n> attach(litter)n> table(dose)ndosen 0 5 50 500 n 20 19 18 17 n> aggregate(weight, by = list(dose), FUN = mean)n Group.1 xn1 0 32.30850n2 5 29.30842n3 50 29.86611n4 500 29.64647n> fit_93 <- aov(weight ~ gesttime + dose)n> summary(fit_93)n Df Sum Sq Mean Sq F value Pr(>F) ngesttime 1 134.3 134.30 8.049 0.00597 **ndose 3 137.1 45.71 2.739 0.04988 * nResiduals 69 1151.3 16.69 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n> n
我們不難發現，協方差分析同樣使用aov()函數，但是formula中的變數順序是將協變數放到實際實際水平變數前面。另外，使用aggregate()函數求均值時，由於均值收到協變數的影響，需要用effects包中的effects()函數計算去除協變數效應後的組均值。
> library(effects)nn載入程輯包：『effects』nnThe following object is masked from 『package:car』:nn Prestigenn> effect("dose",fit_93)nn dose effectndosen 0 5 50 500 n32.35367 28.87672 30.56614 29.33460 n> n
結果表明：(a)懷孕時間與幼崽出生體重相關；(b)控制懷孕時間，藥物劑量與出生體重相關。控制懷孕時間，確實發現每種藥物劑量下幼崽出生體重均值不同。
多重比較使用multcomp包中的glht()函數，mcp參數由rbind(c(3,-1,-1,-1))給出，意思為用量為0的和其它三組用量不為0的作比較，也可以使得c(1,-1,0,0)得到用量0和用量5的做比較。代碼及結果如下：
> library(multcomp)n> comtrast <- rbind("no drrug vs. drug" = c(3, -1, -1, -1))n> summary(glht(fit_93,linfct = mcp(dose = comtrast)))nnt Simultaneous Tests for General Linear HypothesesnnMultiple Comparisons of Means: User-defined ContrastsnnnFit: aov(formula = weight ~ gesttime + dose)nnLinear Hypotheses:n Estimate Std. Error t value Pr(>|t|) nno drrug vs. drug == 0 8.284 3.209 2.581 0.012 *n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n(Adjusted p values reported -- single-step method)n
最後用HH保重的ancova()函數繪製因變數、協變數和因子之間的關係圖。代碼和結果如下：
>library(HH)n載入需要的程輯包：latticen載入需要的程輯包：gridn載入需要的程輯包：latticeExtran載入需要的程輯包：RColorBrewern載入需要的程輯包：gridExtrann載入程輯包：『HH』nnThe following object is masked _by_ 『.GlobalEnv』:nn residplotnnThe following objects are masked from 『package:car』:nn logit, vifnnThe following object is masked from 『package:gplots』:nn residplotnn> ancova(weight ~ gesttime + dose, data = litter)nAnalysis of Variance TablennResponse: weightn Df Sum Sq Mean Sq F value Pr(>F) ngesttime 1 134.30 134.304 8.0493 0.005971 **ndose 3 137.12 45.708 2.7394 0.049883 * nResiduals 69 1151.27 16.685 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n> n
結果表明，在不同的劑量下，重量隨著懷孕時間增加而增加。，還可以看到0劑量組截距項最大，5劑量組截距項最小。由於你上面的設置，直線會保持平行，若用ancova(weight ~ gesttime*dose)，生成的圖形將允許斜率和截距項依據組別而發生變化，這對可視化那些違背回歸斜率同質性的實例非常有用。
> fit_932 <- aov(weight ~ gesttime*dose, data = litter)n> summary(fit_932)n Df Sum Sq Mean Sq F value Pr(>F) ngesttime 1 134.3 134.30 8.289 0.00537 **ndose 3 137.1 45.71 2.821 0.04556 * ngesttime:dose 3 81.9 27.29 1.684 0.17889 nResiduals 66 1069.4 16.20 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n> n
5 雙因素方差分析
與單因素協方差分析類似，雙因素方差分析的區別就在於他們有交互效應，這個與協方差分析中後面的回歸斜率檢驗有些許類似，區別在於交互效應是否顯著。
以基礎安裝中的ToothGrowth數據集為例，隨機分配60隻豚鼠，分別採用兩種餵食方法（橙汁或維生素C），各餵食方法中抗壞血酸含量有三種水平（0.5mg、1mg或2mg），每種處理方式組合都被分配10隻豚鼠。其中牙齒長度為因變數。
> attach(ToothGrowth)nThe following object is masked from litter:nn dosenn> table(supp, dose)n dosensupp 0.5 1 2n OJ 10 10 10n VC 10 10 10n> aggregate(len, by = list(supp,dose), FUN = mean)n Group.1 Group.2 xn1 OJ 0.5 13.23n2 VC 0.5 7.98n3 OJ 1.0 22.70n4 VC 1.0 16.77n5 OJ 2.0 26.06n6 VC 2.0 26.14n> aggregate(len, by = list(supp,dose), FUN = sd)n Group.1 Group.2 xn1 OJ 0.5 4.459709n2 VC 0.5 2.746634n3 OJ 1.0 3.910953n4 VC 1.0 2.515309n5 OJ 2.0 2.655058n6 VC 2.0 4.797731n> dose <- factor(dose)n> fit_96 <- aov(len ~ supp*dose)n> summary(fit_96)n Df Sum Sq Mean Sq F value Pr(>F) nsupp 1 205.4 205.4 15.572 0.000231 ***ndose 2 2426.4 1213.2 92.000 < 2e-16 ***nsupp:dose 2 108.3 54.2 4.107 0.021860 * nResiduals 54 712.1 13.2 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n
結果表明這是一個均衡設計，aggregate()函數獲得各單元的均值和標準差。summary()函數的結果表明主效應和交互效應都非常顯著。
雙因素方差分析可視化有以下三種方法：
> interaction.plot(dose, supp, len, type="b",n col=c("red","blue"), pch=c(16, 18),n main = "Interaction between Dose and Supplement Type")n
> library(gplots)n> plotmeans(len ~ interaction(supp, dose, sep = " "),n+ connect = list(c(1,3,5),c(2,4,6)),n+ col = c("red", "darkgreen"),n+ main = "Interaction Plot with 95% CIs",n+ xlab = "Treatment and Dose Combination")n> n
> library(HH)n> interaction2wt(len~supp*dose)n
三種繪圖方法中，我更推薦HH包中的interaction2wt()函數，因為它能展示任意複雜度設計（雙因素方差分析、三因素方差分析等）的主效應（箱線圖）和交互效應。
6 重複測量方差分析
所謂重複測量方差分析，即測試對象會被重複測量其結果，打個比方就很容易明白了，比如，兩個人在同一年中各月的身高變化，這裡就包含了一個因變數（身高），一個組間因子（人，有兩個，所以是兩個水平），一個組內因子（時間，12個月，就是12個水平），此時，對這兩個人來說（測試對象），他們每個人要測試12次，即組內因子的水平數，所以叫重複測量方差分析。
書中實例是研究生命系統的生理和生化過程如何響應環境因素的變異，對美國北方和南方牧草類植物Echinochloa crus-galli的寒冷容忍度研究結果進行了分析，比較了在某濃度二氧化碳的環境中，對寒帶植物與非寒帶植物的光合作用率。
#將conc轉換成因子變數n> CO2$conc <- factor(CO2$conc)n#去除寒帶植物的數據n> w1b1 <- subset(CO2, Treatment == chilled)n#aov擬合，formula帶Errorn> fit_97 <- aov(uptake ~ conc*Type + Error(Plant/(conc)),w1b1)n#顯示擬合結果n> summary(fit_97)nnError: Plantn Df Sum Sq Mean Sq F value Pr(>F) nType 1 2667.2 2667.2 60.41 0.00148 **nResiduals 4 176.6 44.1 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nnError: Plant:concn Df Sum Sq Mean Sq F value Pr(>F) nconc 6 1472.4 245.40 52.52 1.26e-12 ***nconc:Type 6 428.8 71.47 15.30 3.75e-07 ***nResiduals 24 112.1 4.67 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n#設置坐標軸旋轉n> par(las=2)n#圖面布局n> par(mar=c(10,4,4,2))n#繪製交互效應n> with(w1b1,interaction.plot(conc,Type,uptake,n+ type="b",col = c("red","blue"),pch = c(16,18),n+ main= "Interaction Plot for Plant Type and Concentration"))n#繪製箱型交互圖n> boxplot(uptake ~ Type*conc, data = w1b1,col=(c("gold","green")),n+ main = "Chilled Quebec and Mississippi Plants",n+ ylab = "Carbon dioxide uptake rate (umon/m^2 sec)")n
從以上任意一幅圖都可以看出，魁北克省的植物比密西西比州的植物二氧化碳吸收率高，而且隨著CO2濃度的升高，差異越來越明顯。
本例使用了傳統的重複測量方差分析。該方法假設任意組內因子的協方差矩陣為球形，並且任意組內因子兩水平間的方差之差都相等。但在現實中這種假設不可能滿足，於是衍生了一系列備選方法：
- 使用lme4包中的lmer()函數擬合線性混合模型（Bates，2005）
- 使用car包中的Anova()函數調整傳統檢驗統計量以彌補球形假設的不滿足（例如Geisser-Greenhouse校正）
- 使用nlme包中的gls()函數擬合給定方差-協方差結構的廣義最小二乘模型（UCLA，2009）
- 使用多元方差分析對重複測量數據進行建模（Hand，1987）。
7 多元方差分析
多元方差分析（MANOVA）用於應對因變數（結果變數）不止一個時的情形。以MASS包中的UScereal數據集為例，研究美國穀物中的卡路里、脂肪和糖含量是否會因為儲存架位置的不同而發生變化；其中1代表底層貨架，2代表中層貨架，3代表頂層貨架。卡路里、脂肪和糖含量是因變數，貨架是三水平（1、2、3）的自變數。
> library(MASS)n> attach(UScereal)n> shelf <- factor(shelf)n> y <- cbind(calories, fat, sugars)n> aggregate(y, by = list(shelf), FUN = mean)n Group.1 calories fat sugarsn1 1 119.4774 0.6621338 6.295493n2 2 129.8162 1.3413488 12.507670n3 3 180.1466 1.9449071 10.856821n> cov(y)n calories fat sugarsncalories 3895.24210 60.674383 180.380317nfat 60.67438 2.713399 3.995474nsugars 180.38032 3.995474 34.050018n> fit_98 <- manova(y ~ shelf)n> summary(fit_98)n Df Pillai approx F num Df den Df Pr(>F) nshelf 2 0.4021 5.1167 6 122 0.0001015 ***nResiduals 62 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n> n
需要注意的是，shelf需要因子化，manova()傳入的formula參數～號左邊的y必須用cbind()函數進行合成。從結果上看，F值顯著，三個組的營養成分測量值不同。在多元檢驗是顯著的前提下，可以使用summary.aov()函數對每一個變數再做單因素方差分析。
單因素單因素多元方差分析有兩個前提假設，一個是多元正態性，一個是方差-協方差矩陣同質性。第一個假設即指因變數組合成的向量服從一個多元正態分布。多元正態性可以由QQ圖檢驗，理論補充及代碼如下：
#若有一個p*1的多元正態隨機向量x，均值為μ，協方差矩陣為Σ，那麼x與μ的馬氏距離的平方服從自由度為p的卡方分布。Q-Q圖展示卡方分布n> center <- colMeans(y)n> n <- nrow(y)n> p <- ncol(t)n> cov <- cov(y)n> p <- ncol(y)n> d <- mahalanobis(y,center,cov)n> coord <- qqplot(qchisq(ppoints(n),df = p),n+ d, main = "Q-Q Plot Assessing Multivariate Normality",n+ ylab = "Mahalanobis D2")n> abline(a=0, b=1)n> identify(coord$x, coord$y, labels=row.names(UScereal))n警告: 已經找到了最近的點n警告: 已經找到了最近的點n警告: 已經找到了最近的點n警告: 沒有點在0.25英尺內n警告: 已經找到了最近的點n警告: 已經找到了最近的點n[1] 62 63 64 65n> n
方差-協方差矩陣同質性即指各組的協方差矩陣相同，通常可用Box』s M檢驗來評估該假設，不過遺憾的是目前R中並沒有Box『s M函數。
mvoutlier包中的ap.plot()函數可以用來檢驗多元離群點。代碼和圖形如下：
#需要提前按照mvoutlier包n> library(mvoutlier)n載入需要的程輯包：sgeostatnsROC 0.1-2 loadedn> outliers <- aq.plot(y)nProjection to the first and second robust principal components.nProportion of total variation (explained variance): 0.9789888n> outliersn$outliersn [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUEn[12] TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSEn[23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSEn[34] FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUEn[45] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSEn[56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSEnn> n
> library(rrcov)n載入需要的程輯包：robustbasenn載入程輯包：『robustbase』nnThe following object is masked from 『package:survival』:nn heartnnScalable Robust Estimators with High Breakdown Point (version 1.4-3)nn> Wilks.test(y,shelf,method = "mcd")nn#Windows運行超級慢n> Wilks.test(y, shelf, method = "mcd")nntRobust One-way MANOVA (Bartlett Chi2)nndata: xnWilks Lambda = 0.51073, Chi2-Value = 22.8480, DF = 4.8047,np-value = 0.0003013nsample estimates:n calories fat sugarsn1 119.8210 0.7010828 5.663143n2 128.0407 1.1849576 12.537533n3 160.8604 1.6524559 10.352646n> library(vegan)n載入需要的程輯包：permutenThis is vegan 2.4-2n> adonis(y~shelf)nnCall:nadonis(formula = y ~ shelf) nnPermutation: freenNumber of permutations: 999nnTerms added sequentially (first to last)nn Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) nshelf 2 0.34513 0.172567 6.7092 0.17792 0.002 **nResiduals 62 1.59470 0.025721 0.82208 nTotal 64 1.93983 1.00000 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n
Wilks.test()函數在Windows系統下運行超級慢，E31230 16G內存，大約1分28秒，CPU使用率大約15%，回頭我可以在Ubuntu系統下測試下做個對比。
從結果來看，這兩個方差分析的結果都是顯著的。再一次驗證了貨架的不同層對穀物的三種物質的含量是有影響的。
8 用回歸來做ANOVA
前面提到ANOVA和回歸都是廣義線性模型的特例，本章所有的設計都可以用lm()函數來進行分析，不過由於回歸要求預測變數是數值型，當lm()函數碰到因子時，它會用一系列與因子水平相對應的數值型對照變數來代替因子。如果因子有k個水平，將會創建k–1個對照變數。R提供了五種創建對照變數的內置方法。
> library(multcomp)n> levels(cholesterol$trt)n[1] "1time" "2times" "4times" "drugD" "drugE" n> fit.aov <- aov(response ~ trt, data = cholesterol)n> summary(fit.aov)n Df Sum Sq Mean Sq F value Pr(>F) ntrt 4 1351.4 337.8 32.43 9.82e-13 ***nResiduals 45 468.8 10.4 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1n> n> fit.lm <- lm(response ~ trt, data = cholesterol)n> summary(fit.lm)nnCall:nlm(formula = response ~ trt, data = cholesterol)nnResiduals:n Min 1Q Median 3Q Max n-6.5418 -1.9672 -0.0016 1.8901 6.6008 nnCoefficients:n Estimate Std. Error t value Pr(>|t|) n(Intercept) 5.782 1.021 5.665 9.78e-07 ***ntrt2times 3.443 1.443 2.385 0.0213 * ntrt4times 6.593 1.443 4.568 3.82e-05 ***ntrtdrugD 9.579 1.443 6.637 3.53e-08 ***ntrtdrugE 15.166 1.443 10.507 1.08e-13 ***n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nnResidual standard error: 3.227 on 45 degrees of freedomnMultiple R-squared: 0.7425,tAdjusted R-squared: 0.7196 nF-statistic: 32.43 on 4 and 45 DF, p-value: 9.819e-13n
不難發現，方差分析所得的F檢驗的P值與回歸得到的P值是一樣的，具體原因可以參考這篇文獻，Why ANOVA and Linear Regression are the Same Analysis。通常，解釋變數有連續的，只能用回歸分析。
在回歸分析時的變數對照，可以通過設置contrasts選項來修改lm()中默認的對照方法。
fit.lm <- lm(response ~ trt, data=cholesterol, contrasts="contr.helmert")n
上面的代碼將使用Helmert對照。
也可以通過options()函數修改R會話當中的默認對照方法，舉個例子：
options(contrasts = c("contr.SAS", "contr.helmert"))n
9 小結
我們回顧了基本實驗和准實驗設計的分析方法，包括ANOVA/ANCOVA/MANOVA。然後通過組內和組間設計的示例介紹了基本方法的使用，如單因素ANOVA、單因素ANCOVA、雙因素ANOVA、重複測量ANOVA和單因素MANOVA及各模型的假設檢驗。下一章，我們將介紹功效分析，功效分析可以幫助我們在給定置信度的情況下，判斷達到要求效果所需的樣本大小。
推薦閱讀：