R教學 變異數分析R指令與範例 羅琪老師

變異數分析介紹 變異數分析(ANOVA)能用來分析得自觀察型研究的資料，以檢定三個或三個以上的母體平均數是否相等 2019/8/17 變異數分析介紹 變異數分析(ANOVA)能用來分析得自觀察型研究的資料，以檢定三個或三個以上的母體平均數是否相等 H0: 1=2=3=. . . = k H1: 不是所有母體平均都相等 如果拒絕 H0，我們不能說所有的母體平均數都不相等 拒絕 H0 是指至少有兩個母體平均數不相等

變異數分析的對資料的假設 1.(常態性)每個母體之反應變數均呈常態分配 2.(均質性)所有母體反應變數的變異數σ2均相等 2019/8/17 變異數分析的對資料的假設 1.(常態性)每個母體之反應變數均呈常態分配 2.(均質性)所有母體反應變數的變異數σ2均相等 3.(獨立性)由每個母體抽取之樣本必須互為獨立

2019/8/17

2019/8/17

變異數分析介紹 變異數分析是用來 檢定多組平均數是否相等的問題， 不是在檢定變異數相等的問題 檢定超過兩個以上的平均數的考驗。 2019/8/17 變異數分析介紹 變異數分析是用來 檢定多組平均數是否相等的問題， 不是在檢定變異數相等的問題 檢定超過兩個以上的平均數的考驗。 運用F檢定來檢定平均數間的變異量是否顯著的高於隨機變異量，又稱為變異數分析。 平均數間的變異數(組間變異)除以隨機變異(組內變異)得到的比值(F值)，來取代平均數差異與隨機差異的比值(t或Z值)

資料與敘述統計 總變異=組間變異+組內變異(殘差) n = n1 + n2 +. . . + nk 是總樣本數 1 2 k x11 x21 2019/8/17 資料與敘述統計 n = n1 + n2 +. . . + nk 是總樣本數 總變異=組間變異+組內變異(殘差) 1 2 k x11 x21 xk1 x12 x22 xk2 x1n1 x2n2 xknk 樣本數 n1 n2 nk n 總樣本數 樣本平均數 總平均 樣本變異數

變異的來源 ANOVA將總變異(total variation)劃分成兩部份： 2019/8/17 變異的來源 ANOVA將總變異(total variation)劃分成兩部份： 實驗變數的變異(treatment variation)：衡量在不同的實驗變數(即解釋變數)水準下樣本結果的差異，亦稱「組間變異」(between-groups variation)。 誤差(error)：衡量在個別樣本組內觀察值的變化，亦稱「組內變異」(within-groups variation)。

(Sum of squares between groups) 2019/8/17 SSB代表組間變異平方和 (Sum of squares between groups) SSW代表組內變異平方和 (Sum of squares within groups)

2019/8/17 母體變異數之處理間(組間)估計值 2 的處理間估計值，稱為組間均方 (mean square between groups)，記作 MSB，計算 MSB 的公式如下 k－ 1 為 SSB 的自由度

2019/8/17 母體變異數之處理內(組內)估計值 2 的處理間估計值，稱為組內均方 (mean square within groups)，記作 MSW，或稱為均方誤，記作 MSE，計算的公式如下 n - k 為 SSW 的自由度

2019/8/17 比較變異數之估計值：F 檢定 若虛無假設為真且 ANOVA 之假設均成立， MSB/MSE 的抽樣分配將會服從分子自由度為 k－1，分母自由度為n－k 的 F 分配 不論H0是真是假，則因MSE都不會高估也不會低估σ2 若虛無假設H0為假，則因MSB高估σ2，MSB/MSE 的值將提高

臨界值法：若 F > Fα(k-1, n-k)，則拒絕 H0 2019/8/17 檢定假設 H0: 1=2=3=. . . = k H1: 不是所有母體平均都相等 檢定統計量 F=MSB/MSE 決策法則 p 值法：若 p 值 < α，則拒絕 H0 臨界值法：若 F > Fα(k-1, n-k)，則拒絕 H0

2019/8/17

ANOVA表 變異數分析表(其中k為總組數，n為總樣本數) 單因子變異數分析表 變異 來源 平方和 自由度 2019/8/17 ANOVA表 變異數分析表(其中k為總組數，n為總樣本數) 單因子變異數分析表 變異 來源 平方和 自由度 (Degree of Freedom) 均方 (Mean Square) F值 組間 SSB k-1 MSB=SSB/k-1 F=MSB/MSE 殘差 (組內) SSW n-k MSE=SSW/n-k 總和 SST n-1

若將所有觀察值視為同一組樣本，則總平方和 SST之計算公式為 2019/8/17 若將所有觀察值視為同一組樣本，則總平方和 SST之計算公式為

2019/8/17 多重比較程序 整體檢定(overall test): 當變異數分析F檢定值達顯著水準，即推翻了平均數全相等的虛無假設，亦即表示至少有兩組平均數之間有顯著差異存在 當整體檢定顯著後必須檢驗哪幾個平均數之間顯著有所不同，即進行多重比較(multiple comparison)來檢驗 多重比較在進行F檢定之前進行，稱為事前比較(priori comparisons)，在獲得顯著的F值之後所進行的多重比較，稱為事後比較(posteriori comparisons)

2019/8/17 事後比較 常用的所有成對之比較檢定有： Fisher’s LSD費雪最小顯著差異 (least significant difference, LSD) Scheffe‘s多重比較法 Tukey’s多重比較法 Bonferroni’s多重比較法 可用以決定哪些母體平均數間存在差異

LSD= 𝑡 𝛼 2 (𝑛−𝑘) 𝑀𝑆𝑊( 1 𝑛 𝑖 + 1 𝑛 𝑗 ) 2019/8/17 費雪 LSD 最小顯著差異 檢定假設 H0 : μi = μj H1 : μi ≠ μj 檢定統計量 𝑥 𝑖 − 𝑥 𝑗 決策法則： 若 𝑥 𝑖 − 𝑥 𝑗 > LSD，拒絕 H0 最小顯著差異(least significant difference) LSD= 𝑡 𝛼 2 (𝑛−𝑘) 𝑀𝑆𝑊( 1 𝑛 𝑖 + 1 𝑛 𝑗 ) 缺點:若同時檢定多個配對 型I誤差的機率大於𝜶

LSD= 𝑡 𝛼 2𝑚 (𝑛−𝑘) 𝑀𝑆𝑊( 1 𝑛 𝑖 + 1 𝑛 𝑗 ) , m=k(k-1)/2所有可能配對數 2019/8/17 Bonferroni’s 多重比較法 檢定假設 H0 : μi = μj H1 : μi ≠ μj 檢定統計量 𝑥 𝑖 − 𝑥 𝑗 決策法則： 若 𝑥 𝑖 − 𝑥 𝑗 > LSD，拒絕 H0 最小顯著差異 LSD= 𝑡 𝛼 2𝑚 (𝑛−𝑘) 𝑀𝑆𝑊( 1 𝑛 𝑖 + 1 𝑛 𝑗 ) , m=k(k-1)/2所有可能配對數 優點:若同時檢定m個配對 型I誤差的機率等於𝜶 缺點:t分配表沒那麼精確 T值需用電腦得到

HSD= 𝑞 𝛼 𝑘, 𝑛−𝑘 𝑀𝑆𝑊 𝑛 𝑖 ∗ , 𝑛 𝑖 ∗ =min( 𝑛 𝑖 , 𝑛 𝑗 ) 2019/8/17 Tukey’s 多重比較法 檢定假設 H0 : μi = μj H1 : μi ≠ μj 檢定統計量 𝑥 𝑖 − 𝑥 𝑗 決策法則： 若 𝑥 𝑖 − 𝑥 𝑗 > HSD，拒絕 H0 誠實顯著差異(honest significant difference) HSD= 𝑞 𝛼 𝑘, 𝑛−𝑘 𝑀𝑆𝑊 𝑛 𝑖 ∗ , 𝑛 𝑖 ∗ =min( 𝑛 𝑖 , 𝑛 𝑗 ) 優點:若同時檢定多個配對 型I誤差的機率等於𝜶 缺點:需用專門的表得到q值 http://davidmlane.com/hyperstat/sr_table.html

Scheffe‘s 多重比較法 優點:若同時檢定所有配對或對比 型I誤差的機率等於𝜶 缺點:此法則太保守, 不夠敏感, 不易拒絕H0 2019/8/17 Scheffe‘s 多重比較法 檢定假設 H0 : μi = μj H1 : μi ≠ μj 檢定統計量 𝑥 𝑖 − 𝑥 𝑗 決策法則： 若 𝑥 𝑖 − 𝑥 𝑗 > 𝑘−1 𝐹 𝛼 𝑘−1, 𝑛−𝑘 𝑀𝑆𝑊( 1 𝑛 𝑖 + 1 𝑛 𝑗 ) 則拒絕 H0 優點:若同時檢定所有配對或對比 型I誤差的機率等於𝜶 缺點:此法則太保守, 不夠敏感, 不易拒絕H0

2019/8/17 變異數分析範例 檢定三個方法每星期之平均產量是否相等

在Chemitech公司的例子中，每個樣本數均為 5。使用表 13.1 的資料，我們可以得到下列結果 2019/8/17 在Chemitech公司的例子中，每個樣本數均為 5。使用表 13.1 的資料，我們可以得到下列結果

2019/8/17

若使用顯著水準 α ＝0.05來進行假設檢定，則檢定統計量的值 2019/8/17 若使用顯著水準 α ＝0.05來進行假設檢定，則檢定統計量的值 其分子自由度為 k－1＝3－1＝2，分母自由度為 n－k＝15－3＝12 p 值法：若 p 值 < α，則拒絕 H0 p 值為 F 分配在檢定統計量 F＝9.18 的右尾區域的面積

此檢定結果為三個母體平均數不全相等。換言之，Chemitech公司三種組裝方式所生產的平均產品數量不全相同 2019/8/17 由於 F＝9.18 大於 6.93，因此 F＝9.18的右尾區域會小於 0.01，亦即 p 值小於 0.01。因為 p 值 < α＝0.05，所以拒絕 H0。 此檢定結果為三個母體平均數不全相等。換言之，Chemitech公司三種組裝方式所生產的平均產品數量不全相同

2019/8/17 費雪 LSD 程序(實例) 就 Chemitech 公司之例子而言，LSD 之值為 當樣本大小均相同時，我們只需計算一個 LSD 值。在此情況下，我們僅需將兩樣本平均數之差異值與 LSD 值進行比較。

不拒絕H0，方法A與方法B之母體每週平均產量沒有顯著差異 2019/8/17 檢定假設 檢定統計量 不拒絕H0，方法A與方法B之母體每週平均產量沒有顯著差異

拒絕H0，方法A與方法C之母體每週平均產量有顯著差異 2019/8/17 檢定假設 檢定統計量 拒絕H0，方法A與方法C之母體每週平均產量有顯著差異

拒絕H0，方法B與方法C之母體每週平均產量有顯著差異 2019/8/17 檢定假設 檢定統計量 拒絕H0，方法B與方法C之母體每週平均產量有顯著差異

C A B 52 62 66 方法C與另外兩種方法的平均產量不同 B>C, A>C 2019/8/17 C A B 52 62 66 方法C與另外兩種方法的平均產量不同 B>C, A>C

2019/8/17 變異數分析範例 設陽明貨運公司想從甲、乙、丙、丁四家輪胎廠商中選一家廠商採購輪胎，各從四家廠商隨機抽樣10個輪胎做測試，試問此四家廠商輪胎平均壽命是否有顯著差異？(α= 0.05)

2019/8/17 檢定假設 H0: 1=2=3= 4 H1: 不是所有母體平均都相等

2019/8/17 組間平方和

2019/8/17 組內平方和 檢定統計量的值 其分子自由度為 k－1＝4－1＝3，分母自由度為 n－k＝40－4＝36

2019/8/17 ANOVA表

若 F > Fα(k-1, n-k)= F0.05(3, 36)= 2.87，則拒絕 H0 2019/8/17 決策法則 p 值法：若 p 值 < α=0.05，則拒絕 H0 臨界值法： 若 F > Fα(k-1, n-k)= F0.05(3, 36)= 2.87，則拒絕 H0 由於 F＝3.49 大於 2.87，因此拒絕 H0。 此檢定結果為四家廠商輪胎平均壽命不全相等。換言之，此四家廠商輪胎平均壽命有顯著差異

2019/8/17 變異數分析R範例 > tire<-read.csv("c:/RData/tire.csv", header=T) > tire brand life 1 A 85 2 A 83 3 A 75 4 A 92 5 A 83 6 A 82 7 A 80 8 A 78 9 A 84 10 A 84 11 B 76 12 B 88 13 B 74 14 B 79 15 B 86 16 B 89 17 B 95 18 B 88 19 B 84 20 B 90 21 C 85 22 C 82 23 C 77 24 C 84 25 C 66 26 C 81 27 C 79 28 C 76 29 C 78 30 C 83 31 D 83 32 D 91 33 D 92 34 D 88 35 D 85 36 D 84 37 D 75 38 D 89 39 D 93 40 D 87

2019/8/17 變異數分析R範例 > attach(tire) The following objects are masked from tire (pos = 3): brand, life > str(tire) 'data.frame': 40 obs. of 2 variables: $ brand: Factor w/ 4 levels "A","B","C","D": 1 1 1 1 1 1 1 1 1 1 ... $ life : int 85 83 75 92 83 82 80 78 84 84 ... > detach(tire) Compactly Display the Structure of an Arbitrary R Object 簡潔地顯示任意R物件的結構

2019/8/17 變異數分析R範例 > tire1<-read.csv("c:/RData/tire_1.csv", header=T) > str(tire1) 'data.frame': 40 obs. of 2 variables: $ brand: int 1 1 1 1 1 1 1 1 1 1 ... $ life : int 85 83 75 92 83 82 80 78 84 84 ... > tire1$brand1<-as.factor(tire1$brand1) $ brand1: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ... 廠牌是數字 需將廠牌設定成因子 廠牌是因子

變異數分析R範例 > tire<-read.csv("c:/RData/tire.csv", header=T) > attach(tire) > table(brand) brand A B C D 10 10 10 10

2019/8/17 變異數分析R範例 Apply a function to each cell of a ragged array, that is to each (non-empty) group of values given by a unique combination of the levels of certain factors. > tapply(life, brand, summary) $A Min. 1st Qu. Median Mean 3rd Qu. Max. 75.0 80.5 83.0 82.6 84.0 92.0 $B 74.00 80.25 87.00 84.90 88.75 95.00 $C 66.00 77.25 80.00 79.10 82.75 85.00 $D 75.00 84.25 87.50 86.70 90.50 93.00

2019/8/17 變異數分析R範例 > tapply(life, brand, sd) A B C D 4.526465 6.657494 5.506562 5.313505

分組的盒子圖 boxplot(life~brand, main="tire life by brand", xlab=“brand”, ylab=“life”, col=terrain.colors(4))

分組的盒子圖

變異數分析R範例 SSB=322.475 SSW=1110.300 k-1=4-1=3 n-k=40-4=36 𝝈 =5.553527 2019/8/17 變異數分析R範例 > fit <- aov(life ~ brand, data=tire) > fit Call: aov(formula = life ~ brand, data = tire) Terms: brand Residuals Sum of Squares 322.475 1110.300 Deg. of Freedom 3 36 Residual standard error: 5.553527 Estimated effects may be unbalanced SSB=322.475 SSW=1110.300 k-1=4-1=3 n-k=40-4=36 𝝈 =5.553527

2019/8/17 變異數分析R範例 > anova(fit) Analysis of Variance Table Response: life Df Sum Sq Mean Sq F value Pr(>F) brand 3 322.48 107.492 3.4853 0.02553 * Residuals 36 1110.30 30.842 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 經ANOVA 檢定F值為3.485， p-value = 0.02553<0.05，已達顯著水準，故應拒絕H0 即四家廠商輪胎平均壽命不全相等

2019/8/17 變異數分析R範例 Fitting Linear Models > model <- lm(life ~ brand, data=tire) > model Call: lm(formula = life ~ brand, data = tire) Coefficients: (Intercept) brandB brandC brandD 82.6 2.3 -3.5 4.1 𝐥𝐢𝐟𝐞 =𝟖𝟐.𝟔+𝟐.𝟑𝐛𝐫𝐚𝐧𝐝𝐁−𝟑.𝟓𝐛𝐫𝐚𝐧𝐝𝐂+𝟒.𝟏𝐛𝐫𝐚𝐧𝐝𝐃

2019/8/17 變異數分析R範例 > anova(model) Analysis of Variance Table Response: life Df Sum Sq Mean Sq F value Pr(>F) brand 3 322.48 107.492 3.4853 0.02553 * Residuals 36 1110.30 30.842 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 經ANOVA 檢定F值為3.485， p-value = 0.02553<0.05，已達顯著水準，故應拒絕H0 即四家廠商輪胎平均壽命不全相等

2019/8/17 費雪 LSD 程序(實例) > install.packages("agricolae") > library(agricolae) > fit<-aov(life ~ brand, data=tire) > out <- LSD.test(fit,"brand", p.adj="none") > out $statistics Mean CV MSerror LSD 83.325 6.664899 30.84167 5.037001 $parameters Df ntr t.value alpha test name.t 36 4 2.028094 0.05 Fisher-LSD brand LSD= 𝑡 𝛼 2 (𝑛−𝑘) 𝑀𝑆𝑊( 1 𝑛 𝑖 + 1 𝑛 𝑗 ) LSD <- qt(1-0.025, 40-4)*sqrt(30.84617*(1/10+1/10))=5.037001 t0.025,36=2.0208094

2019/8/17 費雪 LSD 程序(實例) $means life std r LCL UCL Min Max A 82.6 4.526465 10 79.0383 86.1617 75 92 B 84.9 6.657494 10 81.3383 88.4617 74 95 C 79.1 5.506562 10 75.5383 82.6617 66 85 D 86.7 5.313505 10 83.1383 90.2617 75 93 $comparison NULL $groups trt means M 1 D 86.7 a 2 B 84.9 a 3 A 82.6 ab 4 C 79.1 b

2019/8/17 費雪 LSD 程序(實例) D(丁) B(乙) A(甲) C(丙) 86.7 84.9 82.6 79.1 Fisher LSD方法之檢定結果認為廠商為B(乙)與C(丙)之輪胎平均壽命呈現顯著差異，廠商為D(丁)與C(丙) 之輪胎平均壽命呈現顯著差異 乙>丙, 丁>丙

2019/8/17 Scheffe‘s多重比較法 > install.packages("agricolae") > library(agricolae) > fit<-aov(life ~ brand, data=tire) > comparison <- scheffe.test(fit,"brand", group=TRUE, console=TRUE, main="Life of tire with different brand") $statistics Mean CV MSerror CriticalDifference 83.325 6.664899 30.84167 7.282873 $parameters Df ntr F Scheffe alpha test name.t 36 4 2.866266 2.93237 0.05 Scheffe brand diff <- sqrt((4-1)*qf(1-0.05, 4-1, 40-4)*30.84617*(1/10+1/10))=7.282873

2019/8/17 Scheffe‘s多重比較法 $means life std r Min Max A 82.6 4.526465 10 75 92 B 84.9 6.657494 10 74 95 C 79.1 5.506562 10 66 85 D 86.7 5.313505 10 75 93 $comparison NULL $groups trt means M 1 D 86.7 a 2 B 84.9 ab 3 A 82.6 ab 4 C 79.1 b

2019/8/17 Scheffe‘s多重比較法 D(丁) B(乙) A(甲) C(丙) 86.7 84.9 82.6 79.1 Scheffe方法之檢定結果認為廠商為C(丙)與D(丁)之輪胎平均壽命呈現顯著差異 丁>丙

付出最多的人，也是收穫最多的人 ~共勉之~