R語言中的t-test和ANOVA

方差分析(analysis ofvariance，簡寫為ANOVA)是生產和科學研究中分析試驗數據的一種有效的統計方法。引起觀測值不同(波動)的原因主要有兩類:

一類是試驗過程中隨機因素的幹擾或觀測誤差所引起不可控制的的波動；

另一類則是由於試驗中處理方式不同或試驗條件不同引起的可以控制的波動。

方差分析的主要工作就是將觀測數據的總變異(波動)按照變異的原因的不同分解為因子效應與試驗誤差，並對其作出數量分析，比較各種原因在總變異中所佔的重要程度，以此作為進一步統計推斷的依據。

1.aov()的調用格式

aov(formula, data=NULL, projections=FALSE,qr=TRUE,contrasts=NULL, ...)

2. 單因子方差分析

以澱粉為原料生產葡萄的過程中, 殘留許多糖蜜, 可作為生產醬色的原料.在生產醬色的過程之前應儘可能徹徹底底除雜, 以保證醬色質量.為此對除雜方法進行選擇. 在實驗中選用5種不同的除雜方法,每種方法做4次試驗, 即重複4次, 結果見表1

表1 不同除雜方法的除雜量

除雜方法Ai      除雜量Xij              均量Xi

A1            25.6 22.2 28.0 29.8      26.4

A2            24.4 30.0 29.0 27.5      27.7

A3            25.0 27.7 23.0 32.2      27.0

A4            28.8 28.0 31.5 25.9      28.6

A5            20.6 21.2 22.0 21.2      21.3

> X<-c(25.6, 22.2, 28.0,29.8, 24.4, 30.0, 29.0, 27.5, 25.0, 27.7,23.0, 32.2, 28.8, 28.0,31.5, 25.9, 20.6, 21.2, 22.0, 21.2)  #數據

> A<-factor(rep(1:5,each=4))   #分組輸出 A=(1 1 1 1 2 2 2 2 3 33 3 ……5 5 5 5）

>miscellany<-data.frame(X, A)  #拼接

> aov.mis<-aov(X~A,data=miscellany)  #進行anova

>summary(aov.mis)

           Df    Sum Sq  Mean Sq    F value   Pr(>F) 

A          4     132.0    32.99    4.306      0.0162 *

Residuals   15    114.9     7.66              

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*'0.05 '.' 0.1 ' ' 1

上述結果中, 

Df表示自由度; 

sum Sq表示平方和; 

Mean Sq表示均方和;

F value表示F檢驗統計量的值, 即F比; 

Pr(>F)表示檢驗的p值; A就是因素A;

Residuals為殘差.

可以看出, F=4.3061>F0.05(5-1,20-5)=3.06, 或者p=0.01618<0.05,

說明有理由拒絕原假設,即認為五種除雜方法有顯著差異. 

據上述結果可以填寫下面的方差分析表:

方差來源   自由度   平方和    均方和   F比    p值

因素A      4      131.957   32.989   4.3061 0.01618

誤差       15      114.915    7.661

總和       19     246.872

再通過函數plot( )繪圖可直觀描述5種不同除雜方法之間的差異, R中運行命令

> plot(miscellany$X~miscellany$A)

從圖形上也可以看出，5種除雜方法產生的除雜量有顯著差異，特別第5種與前面的4種，而方法1與3，方法2與4的差異不明顯。

3. 雙因子方差分析

3.1 無交互作用的方差分析

在R軟體中, 方差分析函數aov( )既適合於單因素方差分析, 也同樣適用於雙因素方差分析,其中方差模型公式為x  ~A+B,加號表示兩個因素具有可加的. 

原來檢驗果汁中含鉛量有三種方法A1、A2、A3, 現研究出另一種快速檢驗法A4,能否用A4代替前三種方法, 需要通過實驗考察. 觀察的對象是果汁, 不同的果汁當做不同的水平: B1為蘋果, B2為葡萄汁,B3為西紅柿汁, B4為蘋果飲料汁, B5桔子汁, B6菠蘿檸檬汁. 現進行雙因素交錯搭配試驗,即用四種方法同時檢驗每一種果汁,其檢驗結果如表2所示. 問因素A(檢驗方法)和B(果汁品種) 對果汁的含鉛量是否有顯著影響?

因素A    因素B

A       B1     B2     B3     B4    B5    B6    Xi

A1      0.05    0.46  0.12    0.16  0.84   1.30 2.93

A2      0.08    0.38  0.40    0.10  0.92   1.57 3.45

A3      0.11    0.43  0.05    0.10  0.94   1.10 2.73

A4      0.11    0.44  0.08    0.03  0.93   1.15 2.74

>juice<-data.frame(X = c(0.05, 0.46, 0.12,0.16, 0.84, 1.30, 0.08, 0.38, 0.4,0.10, 0.92, 1.57, 0.11, 0.43,0.05, 0.10, 0.94, 1.10,0.11, 0.44, 0.08, 0.03, 0.93, 1.15),A =gl(4, 6),B = gl(6, 1, 24))

#注: 這裡函數gl( )用來給出因子水平,其調用格式為

#gl(n, k, length=n*k,labels=1:n, ordered=FALSE)

#說明: n是水平數, k是每一水平上的重複次數,length是總觀測值數,ordered指明各水平是否先排序.

> juice.aov<-aov(X~A+B, data=juice)

>summary(juice.aov)

           Df Sum Sq      Mean Sq      F value    Pr(>F)    

A           3 0.057       0.0190       1.629      0.225   

B           5 4.902       0.9804       83.976     2e-10***

Residuals   15  0.175 0.0117                

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*'0.05 '.' 0.1 ' ' 1

結論: p值說明果汁品種(因素B)對含鉛量有顯著影響,而沒有充分理由說明檢驗方法(因素A)對含鉛量有顯著影響.

3.2有交互作用的方差分析

R軟體中仍用函數aov( )進行有交互作用的方差分析, 但其中的方差模型格式為x~A+B+A : B.下面用一個例子來全面展示有交互作用方差分析過程.

例：有一個關於檢驗毒品強弱的試驗, 給48隻老鼠注射I、II、III三種毒藥(因素A), 同時有A、B、C、D4種治療方案(因素B), 這樣的試驗在每一種因素組合下都重複四次測試老鼠的存活時間, 數據如表3所示.試分析毒藥和治療方案以及它們的交互作用對老鼠存活時間有無顯著影響.

表3 老鼠存活時間(年)的實驗報告

      A         B           C           D

I   0.31 0.45    0.82 1.10    0.43 0.45     0.45 0.71

    0.46 0.43    0.88 0.72    0.63 0.76    0.66 0.62

II  0.36 0.29    0.92 0.61    0.44 0.35     0.56 1.02

    0.40 0.23    0.49 1.24    0.31 0.40    0.71 0.38

III 0.22 0.21     0.30 0.37    0.23 0.25    0.30 0.36

    0.18 0.23    0.38 0.29    0.24 0.22    0.31 0.33

#建立數據框

> rats<-data.frame(

Time=c(0.31, 0.45, 0.46, 0.43, 0.82, 1.10, 0.88, 0.72, 0.43,0.45,

0.63, 0.76, 0.45, 0.71, 0.66, 0.62, 0.38, 0.29, 0.40,0.23,

0.92, 0.61, 0.49, 1.24, 0.44, 0.35, 0.31, 0.40, 0.56,1.02,

0.71, 0.38, 0.22, 0.21, 0.18, 0.23, 0.30, 0.37, 0.38,0.29,

0.23, 0.25, 0.24, 0.22, 0.30, 0.36, 0.31, 0.33),

Toxicant=gl(3, 16, 48, labels = c("I", "II", "III")),

Cure=gl(4, 4, 48, labels = c("A", "B", "C", "D")))

#下面再用函數interaction.plot( )作出交互效應圖, 以考查因素之間交互作用是否存在, R程序為

> op<-par(mfrow=c(1, 2))

> plot(Time~Toxicant+Cure, data=rats)

> with(rats, interaction.plot(Toxicant, Cure, Time,trace.label="Cure"))

> with(rats,interaction.plot(Cure, Toxicant, Time,trace.label="Toxicant"))

輸出結果如圖a和圖b. 兩圖中的曲線並沒有明顯的相交情況出現,因此我們初步認為兩個因素沒有交互作用.

儘管如此, 由於實驗誤差的存在, 我們用方差分析函數aov()對此進行確認, 其中方差模型格式為x~A*B,或A+B+A:B, 表示不僅考慮因素A、B各自的效應, 還考慮兩者的交互效應.若僅考慮A與B的交互效應則方差模型格式為A:B.

> rats.aov<-aov(Time~Toxicant*Cure, data=rats)

> summary(rats.aov)

             Df   Sum Sq    Mean Sq    F value   Pr(>F)   

Toxicant      2    1.0356    0.5178     23.225   3.33e-07 ***

Cure          3   0.9146    0.3049     13.674   4.13e-06 ***

Toxicant:Cure  6   0.2478    0.0413     1.853     0.116   

Residuals     36   0.8026    0.0223                  

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*'0.05 '.' 0.1 ' ' 1

根據p值知, 因素Toxicant和Cure對Time的影響是高度顯著的,而交互作用對Time的影響卻是不顯著的.