【通俗向】方差分析--幾種常見的方差分析

Ryan_Yang_發表於2017-03-29

上一篇文章說了方差和t檢驗的差異,這篇說說幾種實用的方差分析方法和R語言實現。

一般情況下,基本的方差分析模型包含以下三類,三類下面會根據具體情況再進行細分,主要的三類為一元方差分析,協方差分析,多元方差分析。

1、一元方差分析
一元方差分為單因素、多因素兩類(協方差單獨分類),既然方差是檢驗各組差異的,那麼從一個最簡單的例子入手,探尋各類方差分析的適用條件和特點。

OK,正題開始,鑑於自己也算是酷愛籃球,就舉個籃球運動員的例子吧,以下資料純屬瞎編,如有雷同純屬巧合。

有一天,詹姆斯和科比遇見了,鑑於科比已經退役,詹姆斯說老科,我們比投籃吧,看看是你厲害還是我厲害。科比一聽這話頓時來了精神,就在我家球場,放學別走。

詹姆斯說,我的投籃命中率應該比你高,但是為了避免你蒙進去的多,我們比試5組,每一組投10次,看誰進的多,科比:Deal!

下面是相關程式碼

name<-c(rep('Kobe',5),rep('James',5))
good<-c(5,4,6,5,7,7,5,6,4,8)
a<-data.frame(name,good)
fit1<-aov(good~name,a)
summary(fit1)

其中建立了一個資料框:a,也就是兩個人的投籃比較:

 name good
1   Kobe    5
2   Kobe    4
3   Kobe    6
4   Kobe    5
5   Kobe    7
6  James    7
7  James    5
8  James    6
9  James    4
10 James    8

結果如下:

  Df Sum Sq Mean Sq F value Pr(>F)
name         1    0.9     0.9   0.474  0.511
Residuals    8   15.2     1.9               

可以看到,p=0.511,也就是說,二者的投籃沒有顯著差異(同樣可以使用t檢驗)。

這時候,喬丹老流氓來了,看到他們在比投籃,說,要不我也參與下,好久沒活動了。

隨著老流氓的連續空心,科比和詹姆斯馬上就要報警了,最後看看三人的命中數:

b<-rbind(a,data.frame(name=rep('jordan',5),good=c(10,9,10,10,10)))

name good
1    Kobe    5
2    Kobe    4
3    Kobe    6
4    Kobe    5
5    Kobe    7
6   James    7
7   James    5
8   James    6
9   James    4
10  James    8
11 jordan   10
12 jordan    9
13 jordan   10
14 jordan   10
15 jordan   10

接下來就做一個t檢驗做不了的事情:兩組以上的方差分析:

fit2<-aov(good~name,b)
summary(fit2)

結果是:

 Df Sum Sq Mean Sq F value   Pr(>F)    
name         2  56.93  28.467   21.35 0.000111 ***
Residuals   12  16.00   1.333                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我們看到P值是0.000111,也就是超級顯著了,但是我們知道三者有顯著差別,但是誰和誰有顯著差別暫時還不知道(雖然能猜到,但還是裝不知道),這時候需要藉助Tukey或Duncan方法進行檢驗下,這裡先使用Tukey方法檢測:

TukeyHSD(fit2)

結果為

Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = good ~ name, data = b)

$name
             diff       lwr      upr     p adj
Kobe-James   -0.6 -2.548332 1.348332 0.6973146
jordan-James  3.8  1.851668 5.748332 0.0005973
jordan-Kobe   4.4  2.451668 6.348332 0.0001637

可以看到,Kobe和James沒有差異,喬老爺子和兩個人都有差異,但是和Kobe差異最大,這個用圖形表示就是

plot(TukeyHSD(fit2))

這裡寫圖片描述

圖形中縱軸是均值的差異,包含0說明差異不明顯,不包含0說明差異明顯,可以看到喬丹的均值比詹姆斯和科比都大很多。

那麼這個結果說明了喬丹比科比和詹姆斯投籃準,最後還需要進行正態性假設,因為我們需要知道這5組投籃是隨機從他們職業生涯的千千萬萬個投籃中抽取的,這裡需要檢測下殘差是否符合正態分佈。

到這裡說的有點遠,因為殘差需要進行線性擬合,但是有三個因素性變數,但是我不能把名字作為橫軸吧,那麼喬丹=?科比=?,一般在這種情況下處理就是兩兩比較,比如把詹姆斯的投籃命中作為因變數,喬丹作為自變數1,科比作為自變數2,那麼其實捨棄了科比和喬丹的比較,比如載入car包中,看具體的因子編碼:

library(car)
contrasts(b$name)

結果:

  Kobe jordan
James     0      0
Kobe      1      0
jordan    0      1

可以看到james都是0,也就是因變數,然後Kobe是1,然後喬丹是1,也就是先拿科比作為1,然後保持科比不變,拿喬丹作為1,看和詹姆斯命中的影響,既然看殘差符合正態分佈,所以可以通過兩個方式看qq圖或正態性檢測,先看第一個qqplot:

fit2.1<-lm(good~name,b)
summary(fit2.1)
qqPlot(fit2.1)

先看看fit2.1的情況:

Call:
lm(formula = good ~ name, data = b)

Residuals:
   Min     1Q Median     3Q    Max 
  -2.0   -0.6    0.2    0.4    2.0 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   6.0000     0.5164  11.619 6.92e-08 ***
nameKobe     -0.6000     0.7303  -0.822 0.427336    
namejordan    3.8000     0.7303   5.203 0.000221 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.155 on 12 degrees of freedom
Multiple R-squared:  0.7806,    Adjusted R-squared:  0.7441 
F-statistic: 21.35 on 2 and 12 DF,  p-value: 0.0001115

可以看到,namejordan那一列的p值是0.000221,說明保持科比不變的前提下二者有顯著性差異,但是反而科比和詹姆斯沒啥差異(p=0.427336)。然後看看qq圖:
這裡寫圖片描述

所有殘差都在兩個紅虛線,95%置信區間內,如果不放心,我們再用第二種正態檢測下,並畫出圖看看。

shapiro.test(rstudent(fit2.1))

第一個就是shapiro正態性檢測,結果為

Shapiro-Wilk normality test

data:  rstudent(fit2.1)
W = 0.97576, p-value = 0.9323

看到,有0.9323的概率是正態分佈。然後看密度圖:

plot(density(rstudent(fit2.1)))

這裡寫圖片描述

比較符合正態分佈。

本來想只說說單因素方差分析的,結果說了這麼多。

下面繼續剛才的話題,剛才喬丹加入了,我們看到三個人的投籃命中數:

 name good
1    Kobe    5
2    Kobe    4
3    Kobe    6
4    Kobe    5
5    Kobe    7
6   James    7
7   James    5
8   James    6
9   James    4
10  James    8
11 jordan   10
12 jordan    9
13 jordan   10
14 jordan   10
15 jordan   10

到現在為止,都是單因素方差分析,下面看看雙因素分析;

詹姆斯和科比在看了這篇文章之後,覺得老流氓簡直太厲害了,但是詹姆斯不服科比,我應該比你投籃準,但是因為這是你家後院,你家後院籃球場場地不好太晃眼,我們三個在每個人家的籃球場比比怎麼樣?

老流氓自然無所謂,科比倒是也覺得這樣公平,那麼把三個場地的每個人的投籃命中個數都統計下:

b1<-c(5,4,3,7,5,8,8,9,7,8,9,9,9,7,8)
b2<-c(6,6,7,6,7,8,7,6,7,8,9,10,10,8,9)
c<-cbind(b,b1,b2)
colnames(c)<-c('name','count.J','count.K','count.j')
library(reshape2)
d<-melt(c,variable_name = 'count')

d資料框就是整形好的資料:

 name   count value
1    Kobe count.J     5
2    Kobe count.J     4
3    Kobe count.J     6
4    Kobe count.J     5
5    Kobe count.J     7
6   James count.J     7
7   James count.J     5
8   James count.J     6
9   James count.J     4
10  James count.J     8
11 jordan count.J    10
12 jordan count.J     9
13 jordan count.J    10
14 jordan count.J    10
15 jordan count.J    10
16   Kobe count.K     5
17   Kobe count.K     4
18   Kobe count.K     3
19   Kobe count.K     7
20   Kobe count.K     5
21  James count.K     8
22  James count.K     8
23  James count.K     9
24  James count.K     7
25  James count.K     8
26 jordan count.K     9
27 jordan count.K     9
28 jordan count.K     9
29 jordan count.K     7
30 jordan count.K     8
31   Kobe count.j     6
32   Kobe count.j     6
33   Kobe count.j     7
34   Kobe count.j     6
35   Kobe count.j     7
36  James count.j     8
37  James count.j     7
38  James count.j     6
39  James count.j     7
40  James count.j     8
41 jordan count.j     9
42 jordan count.j    10
43 jordan count.j    10
44 jordan count.j     8
45 jordan count.j     9

count為場地

那麼我們看到三個人在三塊自己的主場上作戰,在這裡需要針對實際情況進行分析,在現有的情況下,我們看的是每個人在每塊場地的投籃命中,也就是投籃命中和場地和人都有關係,所以應該看的是兩個因素的互動項,也就是變數乘積:

fit3<-aov(value~name*count,data=d)
summary(fit3)

看到fit3結果是

  Df Sum Sq Mean Sq F value   Pr(>F)    
name         2  97.91   48.96  47.891 7.18e-11 ***
count        2   2.84    1.42   1.391  0.26181    
name:count   4  18.76    4.69   4.587  0.00427 ** 
Residuals   36  36.80    1.02                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

不同姓名間的投籃命中確實有差別,同時,在不同場地的不同人之間也有顯著差別,那麼跟第一個分析一樣,差別在哪裡呢?

這裡用HH包的互動圖說明下:

library(HH)
interaction2wt(value~name*count,data=d)

圖形如下:
這裡寫圖片描述

先看左上圖,橫軸是三個人,縱軸是命中數量,顏色是三塊場地,從這個圖上看,交錯比較明顯,所以不同場地對所有人來說,差異不是很明顯;和summary中的count p值為0.26相符;

看左下,橫縱不變,顏色為三個人,可以明顯看出,不管在那個場地,Jordan命中都比其他兩人高;

看右上,橫軸為三個場地,場地和命中差異不明顯;

看右下,在任何的場地上,James表現比科比好,而在KOBE主場,這個差異最明顯。

下面,看另一種情況,在前面我們假設的是場地和投籃命中有互動關係,那麼如果三個人住在一起,在每個不同場地投籃的心態都一樣,那麼二者沒有互動,我們的方差分析只看兩個條件的差異:場地的差異和球員的差異:

fit4<-aov(value~name+count,data=d)
summary(fit4)

結果:

 Df Sum Sq Mean Sq F value   Pr(>F)    
name         2  97.91   48.96  35.248 1.49e-09 ***
count        2   2.84    1.42   1.024    0.368    
Residuals   40  55.56    1.39                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

看到就是姓名之間命中率有差異,而場地之間無差異;

這次利用duncan.test進行不同姓名間比對:

duncan.test(fit4,'name',console=T)

結果為

Study: fit4 ~ "name"

Duncan's new multiple range test
for value 

Mean Square Error:  1.388889 

name,  means

          value       std  r Min Max
James  7.066667 1.3345233 15   4   9
jordan 9.133333 0.9154754 15   7  10
Kobe   5.533333 1.2459458 15   3   7

alpha: 0.05 ; Df Error: 40 

Critical Range
        2         3 
0.8697324 0.9144844 

Means with the same letter are not significantly different.

Groups, Treatments and means
a    jordan      9.133 
b    James       7.067 
c    Kobe        5.533 

差異主要看後面最後三行,a/b/c,都是單個的值,則說明三個人都有差異,這次在我修改數值後,三個人都有了較大的差異,比如如果James和Kobe沒差異的話,應該是bc,cb這樣的表示。

2、多元方差分析

其實多元方差只是因變數多了一個,比如我在這個資料集後面加上籃板的引數(實在資料太多,隨機生成吧),

set.seed(1000)
e<-cbind(d,rebound=abs(round(rnorm(45,10,3))))
perform<-cbind(e$value,e$rebound)

最新的資料集是

 name   count value rebound
1    Kobe count.J     5       9
2    Kobe count.J     4       6
3    Kobe count.J     6      10
4    Kobe count.J     5      12
5    Kobe count.J     7       8
6   James count.J     7       9
7   James count.J     5       9
8   James count.J     6      12
9   James count.J     4      10
10  James count.J     8       6
11 jordan count.J    10       7
12 jordan count.J     9       8
13 jordan count.J    10      10
14 jordan count.J    10      10
15 jordan count.J    10       6
16   Kobe count.K     5      11
17   Kobe count.K     4      10
18   Kobe count.K     3      10
19   Kobe count.K     7       4
20   Kobe count.K     5      11
21  James count.K     8      18
22  James count.K     8       6
23  James count.K     9      13
24  James count.K     7      12
25  James count.K     8       8
26 jordan count.K     9      12
27 jordan count.K     9       5
28 jordan count.K     9      11
29 jordan count.K     7      12
30 jordan count.K     8      14
31   Kobe count.j     6       9
32   Kobe count.j     6      12
33   Kobe count.j     7       8
34   Kobe count.j     6       9
35   Kobe count.j     7       5
36  James count.j     8      11
37  James count.j     7       9
38  James count.j     6      13
39  James count.j     7       8
40  James count.j     8       6
41 jordan count.j     9       8
42 jordan count.j    10      14
43 jordan count.j    10      11
44 jordan count.j     8      10
45 jordan count.j     9       9

上資料集的value就是投籃good,很鬱悶reshape2包的value.name不好用了。。。

然後如果我們假設人和場地對籃球運動員的表現有顯著地影響,那麼:

fit5<-manova(perform~name*count,data=e)
summary(fit5)

得出

 Df  Pillai approx F num Df den Df   Pr(>F)    
name        2 0.76600  11.1735      4     72 4.17e-07 ***
count       2 0.14389   1.3954      4     72  0.24419    
name:count  4 0.40204   2.2644      8     72  0.03213 *  
Residuals  36                                            
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看出,還是姓名對錶現有很強的影響,以及人和場地的互動影響。

下面細分來看:

summary.aov(fit5)

結果:

Response 1 :
            Df Sum Sq Mean Sq F value    Pr(>F)    
name         2 97.911  48.956 47.8913 7.178e-11 ***
count        2  2.844   1.422  1.3913  0.261805    
name:count   4 18.756   4.689  4.5870  0.004266 ** 
Residuals   36 36.800   1.022                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response 2 :
            Df  Sum Sq Mean Sq F value Pr(>F)
name         2   9.644  4.8222  0.5865 0.5615
count        2  21.111 10.5556  1.2838 0.2894
name:count   4  14.222  3.5556  0.4324 0.7842
Residuals   36 296.000  8.2222               

這個就和一元方差類似,response1就是投籃,2是籃板,可以看出,投籃跟之前一樣有很大差異,但是籃板和人、場地以及互動都沒有差異。

3.協方差分析

如果考慮一下這種情況:詹姆斯對喬丹不服,說你進聯盟這麼長時間了,所以經驗豐富,科比經驗第二豐富,如果我們剔除掉訓練因素,看看究竟誰天賦更高?這就好像幾個學生考試,小明考了100分,小美考了90分,但是小明天天學習24小時,小美天天學習1小時,如果我們比較小明和小美的本身天賦(及去掉努力的因素),這時候就該用協方差分析。

簡單按照年齡大小排序,喬丹50,科比40,詹姆斯30,ok,建立資料框

f<-cbind(d,age=rep(c(rep(40,5),rep(30,5),rep(50,5)),3))

建好的資料框如下:

name   count value age
1    Kobe count.J     5  40
2    Kobe count.J     4  40
3    Kobe count.J     6  40
4    Kobe count.J     5  40
5    Kobe count.J     7  40
6   James count.J     7  30
7   James count.J     5  30
8   James count.J     6  30
9   James count.J     4  30
10  James count.J     8  30
11 jordan count.J    10  50
12 jordan count.J     9  50
13 jordan count.J    10  50
14 jordan count.J    10  50
15 jordan count.J    10  50
16   Kobe count.K     5  40
17   Kobe count.K     4  40
18   Kobe count.K     3  40
19   Kobe count.K     7  40
20   Kobe count.K     5  40
21  James count.K     8  30
22  James count.K     8  30
23  James count.K     9  30
24  James count.K     7  30
25  James count.K     8  30
26 jordan count.K     9  50
27 jordan count.K     9  50
28 jordan count.K     9  50
29 jordan count.K     7  50
30 jordan count.K     8  50
31   Kobe count.j     6  40
32   Kobe count.j     6  40
33   Kobe count.j     7  40
34   Kobe count.j     6  40
35   Kobe count.j     7  40
36  James count.j     8  30
37  James count.j     7  30
38  James count.j     6  30
39  James count.j     7  30
40  James count.j     8  30
41 jordan count.j     9  50
42 jordan count.j    10  50
43 jordan count.j    10  50
44 jordan count.j     8  50
45 jordan count.j     9  50

先看看年齡的影響

fit6<-aov(value~age+name+count,f)
summary(fit6)

結果為

  Df Sum Sq Mean Sq F value   Pr(>F)    
age          1  32.03   32.03  23.064 2.22e-05 ***
name         1  65.88   65.88  47.432 2.69e-08 ***
count        2   2.84    1.42   1.024    0.368    
Residuals   40  55.56    1.39                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看到count依然沒差異,但是age的差異比較大,先看看包含age的三個人的投籃命中均值:

aggregate(f$value,by=list(f$name),mean)

得出

Group.1        x
1   James 7.066667
2    Kobe 5.533333
3  jordan 9.133333

ok,暫時喬丹領先,那麼看看調整後的均值:

library(effects)
effect('name',fit6)

結果為:

name effect
name
   James     Kobe   jordan 
8.100000 5.533333 8.100000 

看到詹姆斯和喬丹的命中相當,也就是說假以時日,詹姆斯的投籃能趕上喬丹。

OK,這篇文章看起來好長,雖然作為勒布朗的球迷給他貼了金,但是三者的原理有些相似但略有不同,據說線性迴歸和方差分析都是廣義迴歸模型的特例,以後可能回頭再看看有什麼聯絡。

相關文章