一、方差分析的基本概念
方差分析是在20世紀20年代發展起來的一種統計方法，它是由英國統計學家費希爾在進行實驗設計時為解釋實驗數據而首先引入的。
從形式上看，方差分析是比較多個總體的均值是否相等；但是其本質上是研究變量之間的相互關系。方差分析主要用於研究一個數值因變量與一個或多個分類自變量的關系。
根據方差分析的計算方法給方差分析下一個定義：
方差分析（analysis of variance ,ANOVA）就是通過檢驗各總體的均值是否相等來判斷分類型自變量對數值型因變量是否有顯著影響。
二、利用實例講解方差分析的應用
假設某私立學校具有小學、初中、高中三個水平的學生，每個水平階段的學生每個學期都會進行一次期末考試。另外假設小學共抽取10個班級，初中抽取9個班級，高中抽取8個班級。
每個班級期末考試的平均分分別用Ai(i=1,2,...,10)、Bj(j=1,2,...,9)和Ck(k=1,2,...,8)表示，為了比較各個水平階段的班級平均分是否有顯著性差異，就可以使用方差分析。
這里的階段就是因素，也稱之為因子，因子的三個取值：小學、初中、高中，稱為水平或處理。
|
|
階段 |
||
|
|
小學 |
初中 |
高中 |
| 1 |
90 |
87 |
80 |
| 2 |
90 |
88 |
79 |
| 3 |
89 |
90 |
80 |
| 4 |
88 |
78 |
89 |
| 5 |
89 |
89 |
87 |
| 6 |
91 |
90 |
86 |
| 7 |
96 |
80 |
89 |
| 8 |
88 |
81 |
83 |
| 9 |
80 |
82 |
|
| 10 |
90 |
|
|
為了更為直觀地觀察各水平平均值，繪制箱線圖如下：
由於以上數據只涉及到一個分類自變量，即階段，因此屬於單因素方差分析。
從箱線圖可以看出，各水平學生的平均分存在一定的差異，但是這種差異顯著不顯著，還需要進一步分析。
同時，各個水平的方差看起來也不盡相同。
二、方差分析的基本思想
由以上分析可以看出，雖然各個水平的學生平均分存在差異，但是其方差也有差別，方差分析的基本思想就是弄清楚影響因變量取值的誤差來源，以判斷是否是分類自變量對因變量產生影響。
在上述數據中，各組數據的誤差主要來源於以下幾個部分。
首先，即使是同一組的數據，其取值也具有差別，這是因為班級是隨機抽取的，因此他們之間的差異可以看作是隨機因素的影響造成的，或者說是由抽樣的隨機性造成的，這種來自水平內部的誤差稱之為組內誤差，顯然，組內誤差只含有隨機誤差。
其次，各組的取值不同。來自不同水平之間的誤差稱為組間誤差，這種差異可能來自於隨機誤差，也可能來自於因子本身的系統性誤差造成的系統誤差。因此，組內誤差包含有可能包含兩個方面，即隨機誤差和系統誤差。
最后，總誤差為組內誤差與組間誤差之和。
這樣，就把造成因變量的差異的誤差分解成組內誤差和組間誤差。
即
總誤差=組內誤差+組間誤差
如果組內誤差與組間誤差相差太大，說明組間誤差存在很大成分的系統誤差，這時候就可以認為各水平均值顯著不等。
將組間誤差與總誤差的比值定義為關系強度R2，即
R2=
將各平方和除以對應的自由度，則得到相應的均方，也稱為方差。
SST的自由度為n-1
SSA的自由度為k-1
SSE的自由度為n-k
三、方差分析的基本假設
（1）各總體的方差必須相等
（2）各總體必須服從正態分布
（3）各觀測值相互獨立
四、方差分析的類型
根據影響因變量的因素個數，可以把方差分析分為單因素方差分析和雙因素方差分析。
如果是雙因素方差分析，根據兩個因素的交互作為是否對因變量產生影響可分為無交互作用的雙因素方差分析和有交互作用的雙因素方差分析。
五、方差分析的R語言實現
（一）方差分析基本假設的檢驗
將數據在R語言中以列表形式存儲，
> ave_score
$primary
[1] 90 90 89 88 89 91 96 88 80 90
$junior
[1] 87 88 90 78 89 90 80 81 82
$senior
[1] 80 79 80 89 87 86 89 83
1、方差齊性檢驗
（1）Bartlett檢驗（Bartlett檢驗也可以接受一個數據框為輸入，結構與下面的Levene檢驗相同。此方法比較適合用於總體服從正態分布的檢驗）
> bartlett.test(ave_score)
Bartlett test of homogeneity of variances
data: ave_score
Bartlett's K-squared = 0.28233, df = 2, p-value = 0.8683
P值為0.8683，由於p大於常用的a=0.05，因此，無法拒絕原假設，即認為方差相等。
（2）Levene檢驗 (Levene檢驗函數接受的數據結構為數據框結構，並且一列是各水平的取值score，另一列是所屬的水平level，用score~level表示score為因變量，level為自變量。此方法適用於總體非正態的檢驗）
> scores<-data.frame(score=c(ave_score$primary,ave_score$junior,ave_score$senior),level=rep(c("primary","junior","senior"),c(10,9,8)))
> scores
score level
1 90 primary
2 90 primary
3 89 primary
4 88 primary
5 89 primary
6 91 primary
7 96 primary
8 88 primary
9 80 primary
10 90 primary
11 87 junior
12 88 junior
13 90 junior
14 78 junior
15 89 junior
16 90 junior
17 80 junior
18 81 junior
19 82 junior
20 80 senior
21 79 senior
22 80 senior
23 89 senior
24 87 senior
25 86 senior
26 89 senior
27 83 senior
> leveneTest(score~level,data = scores)
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 2 1.0461 0.3668
24
P值大於0.05，可認為等方差
（3）Fligner-Killeen檢驗
> fligner.test(ave_score)（fliger.test函數可以接受列表輸入）
Fligner-Killeen test of homogeneity of variances
data: ave_score
Fligner-Killeen:med chi-squared = 2.0447, df = 2, p-value = 0.3597
結論同上
2、正態性檢驗
shapiro檢驗（輸入數據為一個向量，檢驗該向量的數據是否服從正態分布）
> sapply(ave_score,shapiro.test)（
primary junior
statistic 0.8332102 0.8725879
p.value 0.03657401 0.1310428
method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
data.name "X[[i]]" "X[[i]]"
senior
statistic 0.8764152
p.value 0.1739798
method "Shapiro-Wilk normality test"
data.name "X[[i]]"
primary 水平的檢驗拒絕原假設，即非正態。這個數據只是為了說明方差分析的原理，並沒有進行嚴格的驗證，但是不影響后續工作的推進。
3、獨立性檢驗
可以通過控制抽樣過程來控制獨立性，無具體的檢驗方法。
（二）方差分析
1、單因素方差分析
方差分析所需的數據結構一般是一個數據框，就像上面的那樣。
進行方差分析可以使用lm()函數，也可以使用aov()函數，再利用summary()函數或者anova()函數輸出最終結果。
對於上述的單因素方差分析，分別用這兩者方法分析如下：
（1）用lm()函數
> score_lm<-lm(score~level,data = scores)
> summary(score_lm)
Call:
lm(formula = score ~ level, data = scores)
Residuals:
Min 1Q Median 3Q Max
-9.100 -3.500 0.900 2.938 6.900
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 85.000 1.424 59.706 <2e-16 ***
levelprimaty 4.100 1.962 2.089 0.0475 *
levelsenior -0.875 2.075 -0.422 0.6770
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.271 on 24 degrees of freedom
Multiple R-squared: 0.2309, Adjusted R-squared: 0.1668
F-statistic: 3.602 on 2 and 24 DF, p-value: 0.04285
> anova(score_lm)
Analysis of Variance Table
Response: score
Df Sum Sq Mean Sq F value Pr(>F)
level 2 131.41 65.705 3.6021 0.04285 *
Residuals 24 437.78 18.241
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Summary()函數會輸出殘差和模型，anova()只會輸出結果。
（2）利用aov()函數
> score_aov<-aov(score~level,data = scores)
> summary(score_aov)
Df Sum Sq Mean Sq F value Pr(>F)
level 2 131.4 65.71 3.602 0.0429 *
Residuals 24 437.8 18.24
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(score_aov)
Analysis of Variance Table
Response: score
Df Sum Sq Mean Sq F value Pr(>F)
level 2 131.41 65.705 3.6021 0.04285 *
Residuals 24 437.78 18.241
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
從上述結果可以看出，anova()輸出的是標准的方差分析表，而利用lm()函數計算的方差分析會輸出回歸模型，以及殘差——更多的是回歸的信息。
2、無交互作用的雙因素方差分析
有的時候，因變量可能受到來自一個以上的因素的影響，最典型的就是雙因素方差分析。假如因素A與因素B沒有聯合效應，則稱為無交互作用的雙因素方差分析。
考慮以下的例子，假如某經銷商想知道一款飲料的銷售額與銷售地點和飲料包裝風格的關系，探究這兩個因素是否都在影響銷售額或者只有一個因素影響銷售額。
數據如下：
| 地區 | |||||||
| A | B | C | D | E | 合計 | ||
| 包裝 | 甲 | 30 | 33 | 29 | 41 | 37 | 170 |
| 乙 | 32 | 38 | 39 | 32 | 40 | 181 | |
| 丙 | 29 | 33 | 35 | 39 | 40 | 176 | |
| 丁 | 40 | 41 | 45 | 39 | 40 | 205 | |
| 合計 | 131 | 145 | 148 | 151 | 157 | 732 | |
無交互作用的雙因素方差分析與單因素方差分析類似，只是在寫公式的時候變成【銷售額～地區+包裝】就行。
如果為了省事，可以使用gl()函數生成因子，但是為了與原數據對應，減小閱讀壓力，建議自己寫因子水平。
gl(n, k, length = n*k, labels = seq_len(n), ordered = FALSE)n:一個整數，表示水平的個數
an integer giving the number of levels.
k：表示每個水平重復幾遍an integer giving the number of replications.
length：如果前面兩個都給出了，就不用給這項參數了，否則需要給出所有數據的個數an integer giving the length of the result.
labels：結果的標簽向量an optional vector of labels for the resulting factor levels.
ordered：是否排序a logical indicating whether the result should be ordered or not.
這里，我選擇自己手動生成水平數據。
> dat<-read.table("clipboard",header = TRUE,stringsAsFactors = FALSE)
> dat
包裝 A B C D E
1 甲 30 33 29 41 37
2 乙 32 38 39 32 40
3 丙 29 33 35 39 40
4 丁 40 41 45 39 40
> sales<-c(dat$A,dat$B,dat$C,dat$D,dat$E)
> areas<-rep(c("A","B","C","D","E"),each=4)
> style<-rep(c('甲','乙','丙','丁'),5)
> sales
[1] 30 32 29 40 33 38 33 41 29 39 35 45 41 32 39 39 37 40 40 40
> areas
[1] "A" "A" "A" "A" "B" "B" "B" "B" "C" "C" "C" "C" "D" "D" "D" "D" "E" "E" "E" "E"
> style
[1] "甲" "乙" "丙" "丁" "甲" "乙" "丙" "丁" "甲" "乙" "丙" "丁" "甲" "乙" "丙" "丁"
[17] "甲" "乙" "丙" "丁"
> drink_sale<-data.frame(sales,areas,style)
> drink_sale
sales areas style
1 30 A 甲
2 32 A 乙
3 29 A 丙
4 40 A 丁
5 33 B 甲
6 38 B 乙
7 33 B 丙
8 41 B 丁
9 29 C 甲
10 39 C 乙
11 35 C 丙
12 45 C 丁
13 41 D 甲
14 32 D 乙
15 39 D 丙
16 39 D 丁
17 37 E 甲
18 40 E 乙
19 40 E 丙
20 40 E 丁
分析之前需要對地區和包裝風格做方差齊性檢驗。
> bartlett.test(sales~areas,data = drink_sale)
Bartlett test of homogeneity of variances
data: sales by areas
Bartlett's K-squared = 4.833, df = 4, p-value = 0.3049
> bartlett.test(sales~style,data = drink_sale)
Bartlett test of homogeneity of variances
data: sales by style
Bartlett's K-squared = 2.017, df = 3, p-value = 0.5689
檢驗結果都無法拒絕原假設，即可以認為方差是相等的。
接下來進行方差分析
> drink_aov<-aov(sales~areas+style,data = drink_sale)
> anova(drink_aov)
Analysis of Variance Table
Response: sales
Df Sum Sq Mean Sq F value Pr(>F)
areas 4 93.8 23.450 1.6572 0.22397
style 3 141.2 47.067 3.3263 0.05655 .
Residuals 12 169.8 14.150
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
兩個p值都是大於0.05，因此可以認為銷售地區與包裝風格對銷售額沒有顯著影響。
3、有交互作用的雙因素方差分析
因素之間的交互作用在現實中很常見，比如胖胖的人喜歡藍色的衣服，南方的人更喜歡喝雪花啤酒等，前者是體重和顏色的交互作用，后者是地區和啤酒品牌的交互作用。
因此，如果兩個因素聯合在一起對因變量有顯著的影響，則稱這樣的方差分析為有交互作用的方差分析。
下面的數據展示的是各個路段在高峰期與非高峰期的車流量（數據來自《數據分析：R語言實戰》）。對其進行雙因素方差分析的過程如下。
| 路段1 | 路段2 | 路段3 | |
| 高峰期 | 25 24 27 25 25 | 19 20 23 22 21 | 29 28 31 28 30 |
| 非高峰期 | 20 17 22 21 17 | 18 17 13 16 12 | 22 18 24 21 22 |
> cars<-read.table("clipboard",header = TRUE,stringsAsFactors = TRUE)
> summary(cars)
車流量 路段 時期
Min. :12.00 路段1:10 非高峰期:15
1st Qu.:18.25 路段2:10 高峰期 :15
Median :22.00 路段3:10
Mean :21.90
3rd Qu.:25.00
Max. :31.00
方差齊性檢驗
> bartlett.test(車流量~路段,data=cars)
Bartlett test of homogeneity of variances
data: 車流量 by 路段
Bartlett's K-squared = 0.57757, df = 2, p-value = 0.7492
> bartlett.test(車流量~時期,data=cars)
Bartlett test of homogeneity of variances
data: 車流量 by 時期
Bartlett's K-squared = 0.053302, df = 1, p-value = 0.8174
上述檢驗顯示滿足方差齊性條件
接下來進行方差分析
> cars_aov<-aov(車流量~路段*時期,data=cars)
> anova(cars_aov)
Analysis of Variance Table
Response: 車流量
Df Sum Sq Mean Sq F value Pr(>F)
路段 2 261.600 130.800 35.3514 7.018e-08 ***
時期 1 313.633 313.633 84.7658 2.407e-09 ***
路段:時期 2 6.667 3.333 0.9009 0.4195
Residuals 24 88.800 3.700
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
路段和時期的p值都顯著小於0.05，但是二者的聯合效應的p值為0.4195，大於0.05，因此可以認為二者無顯著的交互作用。
交互效應圖可以更加直觀地看出兩個因素是否具有交互效應，可以用interaction.plot()繪制
interaction.plot(x.factor, trace.factor, response, fun = mean,
type = c("l", "p", "b", "o", "c"), legend = TRUE,
trace.label = deparse(substitute(trace.factor)),
fixed = FALSE,
xlab = deparse(substitute(x.factor)),
ylab = ylabel,
ylim = range(cells, na.rm = TRUE),
lty = nc:1, col = 1, pch = c(1:9, 0, letters),
xpd = NULL, leg.bg = par("bg"), leg.bty = "n",
xtick = FALSE, xaxt = par("xaxt"), axes = TRUE,
...)
x.factor |
a factor whose levels will form the x axis. |
trace.factor |
another factor whose levels will form the traces. |
response |
a numeric variable giving the response |
> names(cars)<-c("flow","path","time")
> attach(cars)
> interaction.plot(path,time,flow,legend = F)
> interaction.plot(time,path,flow,legend = F)
兩個圖中的曲線均沒有相交，可以初步認為沒有交互作用。