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

Presentation on theme: "R教學 變異數分析R指令與範例 羅琪老師."— Presentation transcript:

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

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

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

4 2019/8/17

5 2019/8/17

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

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

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

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

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

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

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

13 臨界值法：若 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

14 2019/8/17

15 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

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

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

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

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

20 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 𝑛 𝑖 𝑛 𝑗 ) , m=k(k-1)/2所有可能配對數 優點:若同時檢定m個配對 型I誤差的機率等於𝜶 缺點:t分配表沒那麼精確 T值需用電腦得到

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

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

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

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

25 2019/8/17

26 若使用顯著水準 α ＝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 的右尾區域的面積

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

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

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

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

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

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

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

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

35 2019/8/17 組間平方和

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

37 2019/8/17 ANOVA表

38 若 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。 此檢定結果為四家廠商輪胎平均壽命不全相等。換言之，此四家廠商輪胎平均壽命有顯著差異

39 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 A 84 B 76 B 88 B 74 B 79 B 86 B 89 B 95 B 88 B 84 B 90 C 85 C 82 C 77 C 84 C 66 C 81 C 79 C 76 C 78 C 83 D 83 D 91 D 92 D 88 D 85 D 84 D 75 D 89 D 93 D 87

40 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": $ life : int > detach(tire) Compactly Display the Structure of an Arbitrary R Object 簡潔地顯示任意R物件的結構

41 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 $ life : int > tire1$brand1<-as.factor(tire1$brand1) $ brand1: Factor w/ 4 levels "1","2","3","4": 廠牌是數字 需將廠牌設定成因子 廠牌是因子

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

43 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 $B $C $D

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

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

46 分組的盒子圖

47 變異數分析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 Deg. of Freedom 3 36 Residual standard error: Estimated effects may be unbalanced SSB= SSW= k-1=4-1=3 n-k=40-4=36 𝝈 =

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

49 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 𝐥𝐢𝐟𝐞 =𝟖𝟐.𝟔+𝟐.𝟑𝐛𝐫𝐚𝐧𝐝𝐁−𝟑.𝟓𝐛𝐫𝐚𝐧𝐝𝐂+𝟒.𝟏𝐛𝐫𝐚𝐧𝐝𝐃

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

51 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 $parameters Df ntr t.value alpha test name.t Fisher-LSD brand LSD= 𝑡 𝛼 2 (𝑛−𝑘) 𝑀𝑆𝑊( 1 𝑛 𝑖 𝑛 𝑗 ) LSD <- qt( , 40-4)*sqrt( *(1/10+1/10))= t0.025,36=

52 2019/8/17 費雪 LSD 程序(實例) $means life std r LCL UCL Min Max A B C D $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

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

54 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 $parameters Df ntr F Scheffe alpha test name.t Scheffe brand diff <- sqrt((4-1)*qf(1-0.05, 4-1, 40-4)* *(1/10+1/10))=

55 2019/8/17 Scheffe‘s多重比較法 $means life std r Min Max A B C D $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

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

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