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)