mean comparison ด้วยแพกเกจ agricolae

      ก่อนหน้านี้ผมได้โพสเกี่ยวกับการวิเคราะห์ ANOVA และการเปรียบเทียบทรีตเมนต์เมื่อพบอิทธิพลของทรีนเมนต์ที่ศึกษามีความสำคัญอย่างมีนัยสำคัญ โดยใช้แพกเกจ multcomp แต่ดูเหมือนจะยังไม่ตอบโจทย์ความต้องการในการวิเคราะห์เมื่อเทียบกับการใช้ซอฟท์แวร์สถิติอื่นๆ โดยเฉพาะการให้สัญลักษณ์ว่าทรีตเมนต์ใดต่างทรีตเมนต์ใดไม่ต่างกัน ในครั้งนี้ผมจะใช้แพกเกจ agricolae ซึ่งช่วยให้การวิเคราะห์ ANOVA และการเปรียบเทียบทรีตเมนต์ทำได้อย่างสะดวกมากทีเดียว

     เริ่มจากที่ข้อมูลที่ใช้เป็นตัวอย่างในไฟล์ Excel ดังรูป

     
     บน RStudio ใช้ File > Import Dataset > From Excel  เข้ามาในชื่อ packing

     ทำการติดตั้งแพกเกจ agricolae และเรียกใช้

> install.packages("agricolae")

> library(agricolae)

  สั่ง attach ชุดข้อมูลที่ต้องการใช้

> attach(packing)  

  
    วิเคราะห์ ANOVA ด้วย aov() ตามปรกติ

> aov.out<-aov(count~treatment)

> summary(aov.out)
            Df Sum Sq Mean Sq F value   Pr(>F)    
treatment    3  32.87  10.958   94.58 1.38e-06 ***
Residuals    8   0.93   0.116                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

  พบว่าทรีตเมนต์มีอิทธิพลอย่างมีนัยสำคัญ จึงทำการเปรียบเทียบค่าเฉลี่ยของแต่ละทรีตเมนต์ เช่น การใช้ LSD จะใช้ฟังก์ชั่น LSD.test()

> LSD.test(aov.out,"treatment",console=TRUE)

Study: aov.out ~ "treatment"

LSD t Test for count 

Mean Square Error:  0.11585 

treatment,  means and individual ( 95 %) CI

           count       std r      LCL      UCL  Min  Max
CO2         3.36 0.3968627 3 2.906844 3.813156 2.91 3.66
commercial  7.48 0.4386342 3 7.026844 7.933156 6.98 7.80
mixed gas   7.26 0.1946792 3 6.806844 7.713156 7.04 7.41
vacuum      5.50 0.2749545 3 5.046844 5.953156 5.26 5.80

alpha: 0.05 ; Df Error: 8
Critical Value of t: 2.306004 

t-Student: 2.306004
Alpha    : 0.05
Least Significant Difference 0.640859
Means with the same letter are not significantly different


Groups, Treatments and means
a   commercial   7.48 
a   mixed gas   7.26 
b   vacuum   5.5 
c   CO2   3.36 

จะเห็นได้ว่าฟังก์ชั่น LSD.test แสดงผลทั้งค่าเฉลี่ยของแต่ละทรีตเมนต์และผลการเปรียบเทียบด้วย LSD ในรูปของตารางในด้านล่างซึ่งช่วยให้ผู้วิเคราะห์สามารถอ่านผลได้ง่ายเช่นเดียวกับซอฟท์แวร์อื่นๆ

เช่นเดียวกัน เราสามารถใช้เทคนิค Duncan ในการเปรียบเทียบทรีตเมนต์ได้ด้วยฟังก์ชั่น duncan.test() ซึ่งในกรณีนี้ให้ผลไม่แตกต่างกัน

> duncan.test(aov.out,"treatment",console=TRUE)

Study: aov.out ~ "treatment"

Duncan's new multiple range test
for count 

Mean Square Error:  0.11585 

treatment,  means

           count       std r  Min  Max
CO2         3.36 0.3968627 3 2.91 3.66
commercial  7.48 0.4386342 3 6.98 7.80
mixed gas   7.26 0.1946792 3 7.04 7.41
vacuum      5.50 0.2749545 3 5.26 5.80

alpha: 0.05 ; Df Error: 8 

Critical Range
        2         3         4 
0.6408590 0.6678356 0.6829140 

Means with the same letter are not significantly different.

Groups, Treatments and means
a   commercial   7.48 
a   mixed gas    7.26 
b   vacuum       5.5 
c   CO2          3.36  

ในแพกเกจ agricolae ยังมีฟังก์ชั่นในการเปรียบเทียบทรีตเมนต์ที่ใช้กันอื่นๆ ไม่ว่าจะเป็น Tukey, SNK, bonferroni,Sche ffe

mean comparison ด้วย R

    หลังจากการวิเคราะห์ ANOVA  เมื่อพบความแตกต่างอย่างมีนัยสำคัญของ treatment เราก็อยากทราบว่า Treatment ใดบ้างที่แตกต่างกัน โดยการวิเคราะห์ที่เรียกว่า post-hoc mean comparison ก่อนหน้านี้ผมได้เขียนเกี่ยวกับการใช้วิธี Tukey สำหรับการเปรียบเทียบค่าเฉลี่ยไปแล้ว ด้วยฟังก์ชั่น TukeyHSD

    ใน R มีแพกเกจหนึ่งชื่อ multcomp ซึ่งใช้ในการเปรียบเทียบค่าเฉลี่ยได้หลายรูปแบบ เช่นจากตัวอย่างเดิม

บน R console

> install.packages("multcomp")
> library(multcomp)

ส่วนนี้เป็นการวิเคราะห์ ANOVA

> file<-file.choose()
> data<-read.csv(file,header=TRUE)
> data
    treatment count
1  commercial  7.66
2  commercial  6.98
3  commercial  7.80
4      vacuum  5.26
5      vacuum  5.44
6      vacuum  5.80
7   mixed gas  7.41
8   mixed gas  7.33
9   mixed gas  7.04
10        CO2  3.51
11        CO2  2.91
12        CO2  3.66
> model <- aov(count~treatment,data=data)
> summary(model)
            Df Sum Sq Mean Sq F value   Pr(>F)   
treatment    3  32.87  10.958   94.58 1.38e-06 ***
Residuals    8   0.93   0.116                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

การเปรียบเทียบแบบ pairwise comparison ด้วยฟังก์ชั่น glht
#ใน  treatment="Tukey" นั้น treatment คือชื่อของส่วนที่เราจะเปรียบเทียบซึ่งในที่นี่ชื่อ treatment ตามข้อมูลด้านบน

> compare<-glht(model,linfct=mcp(treatment="Tukey"))
> summary(compare)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = count ~ treatment, data = data)

Linear Hypotheses:
                            Estimate Std. Error t value Pr(>|t|)   
commercial - CO2 == 0         4.1200     0.2779  14.825   <0.001 ***
mixed gas - CO2 == 0          3.9000     0.2779  14.033   <0.001 ***
vacuum - CO2 == 0             2.1400     0.2779   7.700   <0.001 ***
mixed gas - commercial == 0  -0.2200     0.2779  -0.792    0.856   
vacuum - commercial == 0     -1.9800     0.2779  -7.125   <0.001 ***
vacuum - mixed gas == 0      -1.7600     0.2779  -6.333   <0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

จะเห็นว่าได้ผลลักษณะเดียวกันกับที่ใช้ TukeyHSD 


        เราสามารถทดสอบโดยใช้ Fisher's LSD โดยใช้ summary และกำหนด test=univariate() จะเห็นได้ว่าค่า p ที่ได้จะมีค่าน้อยลงเนื่องจากค่า critical difference ใน LSD จะมีค่าน้อย

> summary(compare,test=univariate())

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = count ~ treatment, data = data)

Linear Hypotheses:
                            Estimate Std. Error t value Pr(>|t|)   
commercial - CO2 == 0         4.1200     0.2779  14.825 4.22e-07 ***
mixed gas - CO2 == 0          3.9000     0.2779  14.033 6.45e-07 ***
vacuum - CO2 == 0             2.1400     0.2779   7.700 5.74e-05 ***
mixed gas - commercial == 0  -0.2200     0.2779  -0.792 0.451410   
vacuum - commercial == 0     -1.9800     0.2779  -7.125 9.95e-05 ***
vacuum - mixed gas == 0      -1.7600     0.2779  -6.333 0.000225 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Univariate p values reported)








หากต้องการใช้ Bonferroni test

> summary(compare,test=adjusted(type="bonferroni"))

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = count ~ treatment, data = data)

Linear Hypotheses:
                            Estimate Std. Error t value Pr(>|t|)   
commercial - CO2 == 0         4.1200     0.2779  14.825 2.53e-06 ***
mixed gas - CO2 == 0          3.9000     0.2779  14.033 3.87e-06 ***
vacuum - CO2 == 0             2.1400     0.2779   7.700 0.000345 ***
mixed gas - commercial == 0  -0.2200     0.2779  -0.792 1.000000   
vacuum - commercial == 0     -1.9800     0.2779  -7.125 0.000597 ***
vacuum - mixed gas == 0      -1.7600     0.2779  -6.333 0.001348 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- bonferroni method)

การใช้ R วิเคราะห์ ANOVA (How to do ANOVA in R)

     แผนการทดลองแบบง่ายที่นักวิจัยใช้กันบ่อยคือแผนการทดลองแบบ CRD (completely Randomized Design) โดยเมื่อเรามีหน่วยทดลองที่มีความเหมือนกันอยู่และสามารถสุ่มจัดสิ่งทดลอง ให้กับทุกหน่วยได้โดยไม่มีข้อจำกัดใดๆ เช่น จากตัวอย่างใน Kuehl (2001) ซึ่ง เป็นการทดลองการบรรจุ 4 แบบ คือบรรจุถุงปรกติ (Commercial) ถุึงสุญญากาศ (Vacuum) บรรจุโดยใช้แก๊สผสม (Mixed gas) และใช้คาร์บอนไดออกไซด์ (CO2) วางแผนการทดลองโดยใช้ชิ้นเนื้อสเต็ก 12 ชิ้น สุ่มบรรจุ 4 แบบ (3 ซ้ำ) ดังแสดงในรูป วัดจำนวนจุลินทรีย์หลังจากเก็บไว้ที่ 4 องศาเซลเซียสเป็นเวลา 9 วัน

โดยไฟล์ csv ที่ได้เตรียมไว้มีหน้าตาดังรูป

เลือกไฟล์แบบ Interactive ด้วย file.choose()

> data <- file.choose()
> test <- read.csv(data)
> test
    treatment count
1  commercial  7.66
2  commercial  6.98
3  commercial  7.80
4      vacuum  5.26
5      vacuum  5.44
6      vacuum  5.80
7   mixed gas  7.41
8   mixed gas  7.33
9   mixed gas  7.04
10        CO2  3.51
11        CO2  2.91
12        CO2  3.66
> str(test)
'data.frame':   12 obs. of  2 variables:
 $ treatment: Factor w/ 4 levels "CO2","commercial",..: 2 2 2 4 4 4 3 3 3 1 ...
 $ count    : num  7.66 6.98 7.8 5.26 5.44 5.8 7.41 7.33 7.04 3.51 ...


      เราจะได้ข้อมูล test ในรูปของ data.frame ซึ่งประกอบด้วยข้อมูล treatment ซึ่งเป็นข้อมูลแบบ factor 4 ระดับ (ตามชื่อ treatment)  และ count เป็นจำนวนจุลินทรีย์ในแต่ละ treatment จำนวน treatment ละ 3 ซ้ำ

     จากนั้นทำการวิเคราะห์ ANOVA ด้วยฟังก์ชั่น aov() และใช้ฟังก์ชั่น summary() เพื่อสรุปข้อมูลให้อยู่ในรูปตาราง ANOVA แบบที่พบโดยทั่วไป

> m <-aov(count~treatment,data=test)
> summary(m)
            Df Sum Sq Mean Sq F value   Pr(>F)   
treatment    3  32.87  10.958   94.58 1.38e-06 ***
Residuals    8   0.93   0.116                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

     จากตาราง ANOVA จะเห็นได้ว่าค่า Pr (>F) มีค่าต่ำมาก (1.38e-6) แสดงว่า treatment มีความแตกต่างกันอย่างมีนัยสำคัญ

      จากนั้นเราสามารถทำการวิเคราะห์ post-anova แบบต่างๆ ได้ เช่นการทำ multiple comparison โดยฟังก์ชันที่มีมากับ R คือฟังก์ชันสำหรับวิเคราะห์โดยใช้วิธี Tukey ด้วยฟังก์ชั่น TukeyHSD()

> TukeyHSD(m,"treatment")
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = count ~ treatment, data = test)

$treatment
                      diff       lwr       upr     p adj
commercial-CO2        4.12  3.230038  5.009962 0.0000020
mixed gas-CO2         3.90  3.010038  4.789962 0.0000031
vacuum-CO2            2.14  1.250038  3.029962 0.0002639
mixed gas-commercial -0.22 -1.109962  0.669962 0.8563618
vacuum-commercial    -1.98 -2.869962 -1.090038 0.0004549
vacuum-mixed gas     -1.76 -2.649962 -0.870038 0.0010160

   จากค่า p เราสามารถดูได้ว่าคู่ของสิ่งทดลองคู่ใดมีความแตกต่างกันบ้าง โดยในตัวอย่างนี้จะเห็นว่ามีเพียงคู่ mixed gas - commercial เท่านั้นที่ non-significant 


เอกสารอ้างอิง
Kuehl RO. 2000. Design of experiment: statistical principles of research design and analysis. Pacific Grove: Duxbury Press. 666 p.