R語言第九章——方差分析

01-28

9.1 術語速成

交互效應、因素方差分析設計、雙因素方差分析、三因素方差分析、組內和組間因子、混合模型方差分析、混淆因素、干擾變數、協變數、協方差分析（ANCOVA）、多元方差分析（MANOVA），多元協方差分析（MANCOVA）。

9.2 ANOVA 模型擬合

9.2.1 aov() 函數

aov() 函數的語法為 aov(formula, data=dataframe)

其中，y 是因變數，字母 A 、 B 、 C 代表因子。

其中，小寫字母表示定量變數，大寫字母表示組別因子， Subject 是對被試者獨有的標識變數。

9.2.2 表達式中各項的順序

一般來說，越基礎性的效應越需要放在表達式前面。具體來講，首先是協變數，然後是主效應，接著是雙因素的交互項，再接著是三因素的交互項，以此類推。

9.3 單因素方差分析

單因素方差分析中，比較分類因子定義的兩個或多個組別中的因變數均值。

此處以multcomp 包中的 cholesterol 數據集為例（取自Westfall、Tobia、Rom、Hochberg，1999），50個患者均接受降低膽固醇藥物治療（ trt ）五種療法中的一種療法。其中三種治療條件使用藥物相同，分別是20mg一天一次（1time）、10mg一天兩次（2times）和5mg一天四次（4times）。剩下的兩種方式（drugD和drugE）代表候選藥物。哪種藥物療法降低膽固醇（響應變數）最多呢？

cholesterol 數據集：

分析代碼：

圖形顯示drugE降低膽固醇最多，而1time降低膽固醇最少。且ANOVA對治療方式（ trt ）的F檢驗非常顯著（p<0.0001），說明五種療法的效果不同。

9.3.1 多重比較

多重比較TukeyHSD() 函數提供了對各組均值差異的成對檢驗。

1time和2times的均值差異不顯著（p=0.138），而1time和4times間的差異非常顯（p<0.001）。

且圖形中置信區間包含0的療法說明差異不顯著（p>0.5），結果跟使用TukeyHSD() 函數得出的結論是一樣的，只是圖形更直觀一些。

multcomp 包中的 glht() 函數提供了多重均值比較更為全面的方法，既適用於線性模型，也適用於廣義線性模型。

有相同字母的組（用箱線圖表示）說明均值差異不顯著。可以看到，1time和2times差異不顯著（有相同的字母a），2times和4times差異也不顯著（有相同的字母b），而1time和4times差異顯著（它們沒有共同的字母）。

9.3.2 評估檢驗的假設條件

可以使用Q-Q圖來檢驗正態性假設：

qqPlot() 要求用 lm() 擬合。如圖所示。數據落在95%的置信區間範圍內，說明滿足正態性假設。

Bartlett檢驗方差齊性：

Bartlett檢驗表明五組的方差並沒有顯著不同（p=0.97）。

用 car 包中的 outlierTest() 函數來檢測離群點：

說明膽固醇數據中不含有離群點（當p>1時將產生 NA ）。因此根據Q-Q圖、Bartlett檢驗和離群點檢驗，該數據似乎可以用ANOVA模型擬合得很好。

9.4 單因素協方差分析

單因素協方差分析（ANCOVA）擴展了單因素方差分析（ANOVA），包含一個或多個定量的協變數。

以 multcomp 包中的 litter 數據集（見Westfall et al.，1999）為例。懷孕小鼠被分為四個小組，每個小組接受不同劑量（0、5、50或500）的藥物處理。產下幼崽的體重均值為因變數，懷孕時間為協變數。

litter 數據集：

代碼分析：

ANCOVA的F檢驗表明：(a)懷孕時間與幼崽出生體重相關；(b)控制懷孕時間，藥物劑量與出生體重相關。控制懷孕時間，確實發現每種藥物劑量下幼崽出生體重均值不同。

使用 multcomp 包來對所有均值進行成對比較。另外， multcomp 包還可以用來檢驗用戶自定義的均值假設。

假如好奇對沒有用藥條件與其他三種用藥條件影響是否不同：

library(multcomp)ncontrast <- rbind("no drug vs. drug" = c(3, -1, -1, -1))nsummary(glht(fit, linfct=mcp(dose=contrast)))n

9.4.1 評估檢驗的假設條件

ANCOVA與ANOVA相同，都需要正態性和同方差性假設，另外，ANCOVA還假定回歸斜率相同。

檢驗回歸斜率的同質性：

9.4.2 結果可視化

HH 包中的 ancova() 函數可以繪製因變數、協變數和因子之間的關係圖。

從圖中可以看到，用懷孕時間來預測出生體重的回歸線相互平行，只是截距項不同。隨著懷孕時間增加，幼崽出生體重也會增加。另外，還可以看到0劑量組截距項最大，5劑量組截距項最小。

更改代碼：

生成的圖形斜率和截距項依據組別而發生變化，不再是互相平行。

9.5 雙因素方差分析

在雙因素方差分析中，受試者被分配到兩因子的交叉類別組中。

以ToothGrowth 數據集為例，隨機分配60隻豚鼠，分別採用兩種餵食方法（橙汁或維生素C），各餵食方法中抗壞血酸含量有三種水平（0.5mg、1mg或2mg），每種處理方式組合都被分配10隻豚鼠。其中，牙齒長度為因變數。

雙因素ANOVA：

數據集：

分析代碼：

還可以用 gplots 包中的 plotmeans() 函數來展示交互效應：

還能用 HH 包中的 interaction2wt() 函數來可視化結果，圖形對任意順序的因子設計的主效應和交互效應都會進行展示。

以上三幅圖形都表明隨著橙汁和維生素C中的抗壞血酸劑量的增加，牙齒長度變長。對於0.5mg和1mg劑量，橙汁比維生素C更能促進牙齒生長；對於2mg劑量的抗壞血酸，兩種餵食方法下牙齒長度增長相同。

9.6 重複測量方差分析

所謂重複測量方差分析，即受試者被測量不止一次。

以CO2 數據集為例，CO2 數據集包含了北方和南方牧草類植物Echinochloa crus-galli （Potvin、Lechowicz、Tardif，1990）的寒冷容忍度研究結果，在某濃度二氧化碳的環境中，對寒帶植物與非寒帶植物的光合作用率進行了比較。研究所用植物一半來自於加拿大的魁北克省（Quebec），另一半來自美國的密西西比州（Mississippi）。

首先，我們關注寒帶植物。因變數是二氧化碳吸收量（ uptake ），單位為ml/L，自變數是植物類型 Type （魁北克VS.密西西比）和七種水平（95~1000 umol/m^2 sec）的二氧化碳濃度（ conc ）。另外， Type 是組間因子， conc 是組內因子。 Type 已經被存儲為一個因子變數，但你還需要先將conc 轉換為因子變數。

含一個組間因子和一個組內因子的重複測量方差分析：