R語言從入門到精通：Day11

• 2019 年 10 月 6 日
• 筆記

在上一次推文中，我們已經介紹了兩組獨立樣本的t檢驗，今天我們來介紹用於常見實驗設計的方差分析大全。

在上一次推文中，我們已經介紹了兩組獨立樣本的t檢驗，那麼多組（大於兩組）樣本的比較就要用到這次推文中介紹的方差分析（ANOVA）了。方差分析在各種實驗和准實驗設計的分析中都有廣泛應用，接下來將要介紹用於常見研究設計分析的 R函數。

在開始閱讀之前，我建議大家能回頭去複習一下方差分析的基礎知識，特別是一些基本術語，對接下來的學習會有很大幫助。

溫馨提示

1、本節內容重點內容較多，

務必緊跟紅色標記。

2、測試數據及程式碼

見文末客服小姐姐二維碼。

1、基礎方差分析

ANOVA和回歸方法（下一次推文的主要內容）雖然都是獨立發展而來，但是從函數形式上看，它們都是廣義線性模型的特例。用學習繪圖時用到的函數lm()也能分析ANOVA模型。不過，這裡我們基本都使用函數aov()。兩個函數的結果是等同的，但函數aov()函數展示結果的格式更容易理解。為保證完整性，最後會提供一個使用函數lm()的例子供大家參考。

函數aov()的用法和之前接觸過的函數略有不同，語法為aov(formula, data=dataframe)，參數formula代表了方差分析實驗設計中將要研究的變數之間的關係。表1和表2分別列舉了表達式中常見的符號和常見的表達式。

表1：表達式中常見符號

表2:常見表達式及對應實驗設計。

重點：表達式右邊的式子並不符合「交換律」，即式子右邊的變數順序很重要。變數順序在以下兩種情況下會造成影響:（1）因子不止一個，並且是非平衡設計;(2) 存在協變數。例如，對於雙因素方差分析，若不同處理方式中的觀測數不同， 那麼模型y ~ A*B與模型y ~ B*A的結果不同。此外，以表達式y~A+B+A:B為例，R會默認這樣理解它（序貫型）：（1）A對y的影響；（2）控制A時，B對y的影響;（3）控制A和B的主效應時，A與B的交互效應。（對於任意表達式而言，有三種理解方式，分別是序貫型、分層型、邊界型。其中R默認調用序貫型，而有些統計軟體（如SPSS）默認調用邊界型。也不用擔心必須應用其它兩種理解方式的情形，R中提供了很多函數包來應對（如包car中的函數Anova()）。）

2、單因素方差分析

從最簡單的情況出發，單因素方差分析中，你感興趣的是比較分類因子定義的兩個或多個組別中的因變數均值。以 multcomp包中的cholesterol數據集為例（實驗設計為：50 個患者均接受降低膽固醇藥物治療(trt)五種療法中的一種療法。其中三種治療條件使用藥物 3 相同，分別是20mg一天一次(1time)、10mg一天兩次(2times)和5mg一天四次(4times)。剩下的兩種方式(drugD和drugE)代表候選藥物。哪種藥物療法降低膽固醇(響應變數)最多呢?）。利用單因素方差分析對不同療法做F檢驗，結果非常顯著（p<0.0001）的說明五種療法的效果不同。（利用包gplot中的函數plotmeans()繪製了組均值圖形，圖中也可以清楚的看到它們之間的差異。程式碼見文末找胖雨小姐姐。）

圖1:單因素方差分析，五種療法的均值差異

上面的分析之後，你肯定會對究竟哪種療法與其他療法不同感興趣。多重比較可以解決這個問題。函數TukeyHSD()可以回答你的疑問。利用函數plot()將多重比較結果繪製為圖形，圖形中置信區間包含0的療法說明差異不顯著(p>0.5)。

圖2:多重比較示例1

包multcomp中的函數glht()函數提供了多重均值比較更為全面的方法，下面是圖形化的分析結果。程式碼中，函數cld()中的參數level設置了使用的顯著水平(0.05，即本例中的95%的置信區間)。有相同字母的組(用箱線圖表示)說明均值差異不顯著。（做完上面這些工作並不意味著分析結束了，必須要檢查數據是否滿足統計檢驗的假設條件的程度，單因素方差分析中，我們假設因變數服從正態分布，各組方差相等。可以使用Q-Q圖（函數qqPlot()）來檢驗正態性假設，用函數bartlett.test()檢查方差齊性。後台程式碼已提供示例。）

圖3：多重比較示例2

3、單因素協方差分析

單因素協方差分析(ANCOVA)擴展了單因素方差分析(ANOVA)，包含一個或多個定量的協變數。下面的例子來自於包multcomp中的litter數據集。懷孕小鼠被分為四個小組，每個小組接受不同劑量(0、5、50或500)的藥物處理。產下幼崽的體重均值為因變數，懷孕時間為協變數。ANCOVA的F檢驗表明:(1)懷孕時間與幼崽出生體重相關;(2)控制懷孕時間，藥物劑量與出生體重相關。控制懷孕時間，確實發現每種藥物劑量下幼崽出生體重均值不同。（由於使用了協變數，你可能想要獲取調整的組均值，即去除協變數效應後的組均值。可使用 effects包中的effects()函數來計算調整的均值。）

圖4:單因素協方差分析示例

與單因素方差分析類似，劑量的F檢驗雖然表明了不同的處理方式幼崽的體重均 值不同，但無法告知我們哪種處理方式與其他方式不同。用multcomp包來對所有均值進行成對比較或者對你感興趣的對照進行比較分析（程式碼中提供了未用藥與用藥條件之間的比較）。單因素協方差分析也需要正態性和方差齊性假設，檢查方法與單因素方差分析中一樣，此外，單因素協方差分析還需要假定回歸斜率相同，示例如下。

圖5：單因素協方差分析的斜率相同檢查

圖5中，可以看到交互效應不顯著，支援了斜率相等的假設。若假設不成立，可以嘗試變換協變數或 因變數，或使用能對每個斜率獨立解釋的模型，或使用不需要假設回歸斜率同質性的非參數 ANCOVA方法。包sm中的函數sm.ancova()為後者提供了一個例子。還可以用包HH中的函數ancova()對單因素協方差分析的結果進行可視化。從圖6中可以看出，用懷孕時間來預測出生體重的回歸線相互平行，只是截距項不同。隨著懷孕時間增加，幼崽出生體重也會增加。另外，還可以看到0劑量組截距項最大，5劑量組截距項最小。（不妨嘗試一下將程式碼中ancova(weight~gesttime+dose)改為ancova(weight~gesttime*dose)，圖片出現了什麼變化呢？）

圖6:單因素協方差分析可視化

4、雙因素方差分析

討論完單因素方差分析，我們來看一下更複雜的情形：雙因素方差分析和重複測量方差分析。在雙因素方差分析中，受試者被分配到兩因子的交叉類別組中。以基礎安裝中的 ToothGrowth數據集為例，隨機分配60隻豚鼠，分別採用兩種餵食方法(橙汁或維生素C)，各餵食方法中抗壞血酸含量有三種水平(0.5mg、1mg或2mg)，每種處理方式組合都被分配10隻豚鼠。其中，牙齒長度為因變數。分析結果可以採用函數interaction.plot()、包gplots中的函數plotmeans()或者包HH中的函數interaction2wt()來可視化，下面展示一下函數interaction2wt()的例子。（三種方法中，更推薦包HH中的函數interaction2wt()，因為它能展示任意複雜度設計(雙因素方差分析、三因素方差分析等)的主效應(箱線圖)和交互效應。程式碼中提供了三種方法的示例，大家可以自己選擇偏好的方式。）

圖7:函數interaction2wt()示例

5、重複測量方差分析

而所謂重複測量方差分析，即受試者被測量不止一次。示例來源於基礎安裝包中的CO2數據集，包含了北方和南方牧草類植物Echinochloa crus-galli的寒冷容忍度研究結果，在某濃度二氧化碳的環境中，對寒帶植物與非寒帶植物的光合作用率進行了比較。因變數是二氧化碳吸收量(uptake)，自變數是植物類型Type和七種水平的二氧化碳濃度(conc)。另外，Type是組間因子，conc是組內因子。Type已經被存儲為一個因子變數，還需要將 conc轉換為因子變數。方差分析表表明在0.01的水平下，主效應類型和濃度以及交叉效應類型×濃度都非常顯著，圖8中通過函數boxplot()展示了交互效應。程式碼找文末客服小姐姐。

圖8:函數boxplot()的結果可視化

6、多元方差分析

前面我們討論都是單個因變數的情形，當因變數(結果變數)不止一個時，可用多元方差分析(MANOVA)對它們同時進行分析。以MASS包中的UScereal數據集為例，研究穀物中的卡路里、脂肪和糖含量是否會因為儲存架位置的不同而發生變化。卡路里、脂肪和糖含量是因變數，貨架是三水平(1、2、3)的自變數。函數manova()能對組間差異進行多元檢驗。方差分析表中F值顯著，說明三個組的營養成分測量值不同。函數summary.aov()可以對每一個變數做單因素方差分析。（單因素多元方差分析有兩個前提假設，一個是多元正態性，一個是方差協方差矩陣同質性，可以用Q-Q圖來檢驗假設條件。還可以使用包mvoutlier中的函數ap.plot()來檢驗多元離群點。方差分析及假設檢驗程式碼見後台。）

每種方差分析後的評估假設檢驗都是很重要的，前面的例子都通過了檢查，可是現實中總會出現不符合假設條件的情況，那麼可以考慮用穩健或非參數版本的MANOVA檢驗。穩健單因素MANOVA可通過包rrcov中的 函數Wilks.test()實現。包vegan中的函數adonis()則提供了非參數MANOVA的等同形式。程式碼中提供函數Wilks.test()的例子。

最後，我們來兌現教程最開始的「諾言」，給大家介紹一個用回歸來做的方差分析。實際上，ANOVA和回歸都是廣義線性模型的特例，前面所有的設計都可以用函數lm()來分析。以單因素方差分析實例為例，即比較五種降低膽固醇藥物療法(trt)的影響。結果如下：

圖9:函數lm()做方差分析

小結

本次教程內容比較多，包括了單因素ANOVA、單因素ANCOVA、 雙因素ANOVA、重複測量ANOVA和多因素MANOVA。除了這些基本分析，還介紹了模型的假設檢驗，以及應用多重比較過程來進行綜合檢驗的方法，對各種結果可視化方法也進行了探索。下一次的重點將是回歸分析。還在堅持學習的各位，要加油哦～