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

การเปิดไฟล์ csv ด้วย R แบบ interactive (How to select file interactively in R)

     โดยทั่วไปการเรียกเปิดไฟล์ csv (comma separated value)  ซึ่งมักจะใช้เก็บข้อมูลสำหรับการวิเคราะห์ (โดยใช้โปรแกรมอย่าง MS Excel) เรามักนิยมการสั่งคำสั่งด้วย read.csv() ซึ่งมักจะเจอปัญหาว่าเราต้องพิมพ์ชื่อไฟล์ให้ถูกต้อง วิธีการที่ง่ายกว่าคือการใช้ฟังก์ชั่น file.choose() เพื่อเลือกเปิดไฟล์แบบ interactive นั่นคือสามารถ browse หาไฟล์ที่ต้องการได้และใช้ mouse ในการคลิกเลือกไฟล์ที่ต้องการ เช่น สมมติว่าเรามีข้อมูลในไฟล์ Book1.csv ซึ่งสร้างบน Excel ดังรูป

บน R console  เราจะสั่งให้เปิดไฟล์ขึ้นมาเพื่อให้เก็บในชื่อ data

> data <- file.choose()


จะเห็นหน้าต่าง dialog box ขึ้นมาเพื่อให้เลือกไฟล์ได้โดยใช้ mouse

 
ลองดูข้อมูลของ data

> data
[1] "C:\\Users\\Pramote\\Documents\\Book1.csv"

อ่านไฟล์ csv ในชื่อ data เก็บในชื่อ testdata

> testdata <- read.csv(data)

ดูข้อมูลของ testdata ซึ่งอยู่ในรูปของ data.frame


> testdata
  x  y
1 1 10
2 2 20
3 3 30

> str(testdata)
'data.frame':   3 obs. of  2 variables:
 $ x: int  1 2 3
 $ y: int  10 20 30

ทดสอบ plot ข้อมูล

 > plot(testdata$x,testdata$y)