R語言數據分析與挖掘(第五章):方差分析(2)——多因素方差分析
- 2019 年 12 月 16 日
- 筆記
在實際應用中，更多出現的是包含多因素的試驗和處理。多因素試驗與雙因素試驗背後的基本思想是一致的。與單因素方差分析不同，在雙因素方差分析中因素間可能會有交互作用。假設有兩個因素A和B，因素A和B沒有交互作用指的是A的水平值不取決於B的水平值，反之亦然。對於有交互作用的因素，我們不可孤立地看待這些因素。對於雙因素的情形，一般從影像上看，沒有交互作用的因素水平圖表現為兩條不相交的線段，而有交互作用的因素水平圖為兩相交的線段。例如，下圖顯示的是在研究年齡和性別對身高是否有顯著作用過程中，因素年齡與性別之間的交互作用。從影像上看，兩曲線沒有明顯相交，據此可以推測二者間不存在相互作用。當然，要判定是否存在或者不存在交互作用，還需要根據相應的統計量來分析。
如果因素A和B間不存在交互作用，那麼可以對因素A和B各自進行獨立分析，在後續的分析中去除不顯著的因素。如果方差分析的結果顯示因素A和B間存在交互作用，這時候要對數據進行進一步的分析，具體包括：
在因素A的某個水平下，因素B對響應變數的作用。
在因素B的某個水平下，因素A對響應變數的作用。
在所有因素（A，B）的組合中，哪兩組的差異最大。
需要指出的是，多因素的情況與雙因素處理方法類似，只不過分析的因素會增多。
下面以不同劑量下老鼠妊娠重量的差異性分析為例進行介紹。
在R語言中，實現雙因素方差分析的函數與實現單因素方差分析的函數一致，可以實現aov()和anova()函數，不同之處在於模型公式的設定，雙因素方差分析的模型公式應設定為"X~A+B"或"X~A*B"的形式，後者表示考慮因素A與B的交互作用。
下面利用包multcomp中的litter數據集進行操作演練，該數據集是關於老鼠妊娠重量與劑量反應的研究數據，共74個行觀測值和4個變數，列變數包括，劑量等級，共有4個等級，0,5,50和500，妊娠時間，作為協變數共有4個水平，分別是21.5、22、22.5和23；妊娠動物的個數作為協變數；整個妊娠過程中出生的小鼠平均重量，數據的基本情況如下：
> library(multcomp) > data(litter) > dim(litter) [1] 74 4 > head(litter) dose weight gesttime number 1 0 28.05 22.5 15 2 0 33.33 22.5 14 3 0 36.37 22.0 14 4 0 35.52 22.0 13 5 0 36.77 21.5 15 6 0 29.60 23.0 5 > summary(litter) dose weight gesttime number 0 :20 Min. :19.22 Min. :21.50 Min. : 5.00 5 :19 1st Qu.:27.77 1st Qu.:21.50 1st Qu.:12.00 50 :18 Median :30.76 Median :22.00 Median :14.00 500:17 Mean :30.33 Mean :22.09 Mean :13.43 3rd Qu.:33.30 3rd Qu.:22.50 3rd Qu.:15.00 Max. :38.75 Max. :23.00 Max. :17.00 > table(litter$gesttime) 21.5 22 22.5 23 20 24 27 3 > aggregate(litter$weight,by=list(litter$dose),FUN=mean) Group.1 x 1 0 32.30850 2 5 29.30842 3 50 29.86611 4 500 29.64647 > aggregate(litter$weight,by=list(litter$dose),FUN=sd) Group.1 x 1 0 2.695119 2 5 5.092352 3 50 3.762529 4 500 5.404372 > par(mfrow=c(1,2)) > boxplot(weight~dose,data=litter,col=2:5,main="按變數dose分組") > boxplot(weight~gesttime,data=litter,col=2:5,main="按變數gesttime分組")
上述程式碼表示:利用函數summary()對數據集進行描述性統計，發現變數dose的4個不同水平的數量不相等，分別為20、19、18和17;需要注意的是，原始數據中變數gesttime中存儲的數據類型為數值型，故利用函數table()統計該變數4個水平的數據量，發現妊娠時間為21.5、22、22.5和23的老鼠數量分別為20、24、27和30;利用函數aggregate()對變數weight進行分組統計，並計算每-組的均值和方差， 分組依據為變數dose的4個水平，並根據兩個協變數的不同分組繪製變數weight的箱線圖，結果下圖所示。
下面對數據進行方差齊性檢驗，由於原始變數gesttime的數據類型為數值型，在進行方差齊性檢驗的時候需要將其轉化成因子，具體操作如下：
> bartlett.test(weight~as.factor(gesttime),data=litter) Bartlett test of homogeneity of variances data: weight by as.factor(gesttime) Bartlett's K-squared = 0.26778, df = 3, p-value = 0.966 > bartlett.test(weight~dose,data=litter) Bartlett test of homogeneity of variances data: weight by dose Bartlett's K-squared = 9.6497, df = 3, p-value = 0.02179
輸出結果顯示：變數gesttime的值為0.966，大於給定的顯著性水平0.05，故不能拒絕原假設，即認為不同水平下是等方差的，然而dose的p值為0.02179，沒有通過檢驗。接下來對數據進行方差分析，分別按照不考慮gesttime和dose的交互和考慮其交互進行分析。
> fit1<-aov(weight~gesttime+dose,data=litter) > summary(fit1) 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
在不考慮交互作用的情況下，協變數gesttime和dose的檢驗p值均小於給定的顯著水平，故拒絕原假設，即認為變數gesttime和dose對變數weight的影響均顯著的。
> fit2<-aov(weight~gesttime*dose,data=litter) > summary(fit2) Df Sum Sq Mean Sq F value Pr(>F) gesttime 1 134.3 134.30 8.289 0.00537 ** dose 3 137.1 45.71 2.821 0.04556 * gesttime:dose 3 81.9 27.29 1.684 0.17889 Residuals 66 1069.4 16.20 --- Signif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1
上述程式碼表示:在考慮交互作用的情況下進行方差分析，結果顯示，交互作用項的p值為0.17889，大於給定的顯著性水平0.05,故不認為交互項對交量weigh的影響是顯著的，該交互項可以刪除。
交互作用的效果還可以進行可視化展示圖，利用程輯包HH中的函數intrecion2wt()即可，其函數的用法與函數aov()類似，具體操作程式碼如下:
install.packages("HH") library(HH) interaction2wt(weight~gesttime*dose,data=litter)
結果如下圖所示。主要看圖形的第一、四象限的曲線是否存在明顯相交的情況，若存在，則說明兩因素間的交互作用顯著，否則認為不顯著，本圖中第一、四象限的曲線有一定的相交，說明在後續的方差分析中需要添加交互項，但是根據前面的分析結果，交互項並沒有通過檢驗，即交互項的作用並不明顯。
