การใช้ R ในการคำนวณ ANOVA

ปราโมทย์ คูวิจิตรจารุ
เพิ่มเติมล่าสุด 7 มกราคม 2552

Note สำหรับการคำนวณ ANOVA ตามแผนการทดลองแบบต่างๆ ด้วย R โดยการยกตัวอย่างให้เห็น ผู้สนใจอาจจะนำไปปรับใช้กับงานของตนเองได้ครับ ผู้ที่สนใจอ่านพื้นฐานการใช้ลองดูที่หน้านี้

หัวข้อ

One-factor CRD

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

ในการวิเคราะห์ ANOVA ด้วย R นั้นเราจะป้อนข้อมูล (อยู่ในรูปที่เรียกว่าเว็คเตอร์) 2 เว็คเตอร์ คือ treatment เป็นเว็คเตอร์ของชื่อสิ่งทดลองของเรา และ count เป็นจำนวนจุลินทรีย์ในแต่ละสิ่งทดลอง ข้อมูล treatment จะต้องถูกแปลงให้เป็น factor นั่นคือมีการแบ่งเป็น level

> treatment <- c("comm","comm","comm","vacc","vacc","vacc","mixed","mixed","mixed","co2","co2","co2")
> treatment <-factor(treatment)
> treatment
 [1] comm  comm  comm  vacc  vacc  vacc  mixed mixed mixed co2   co2   co2  
Levels: co2 comm mixed vacc
> count <- c(7.66,6.98,7.8,5.26,5.44,5.8,7.41,7.33,7.04,3.51,2.91,3.66)

จากนั้นจัดข้อมูลของทั้งสองเว็คเตอร์ให้อยู่ในรูปของ dataframe

> crd <- data.frame(treatment,count)
> crd
   treatment count
1       comm  7.66
2       comm  6.98
3       comm  7.80
4       vacc  5.26
5       vacc  5.44
6       vacc  5.80
7      mixed  7.41
8      mixed  7.33
9      mixed  7.04
10       co2  3.51
11       co2  2.91
12       co2  3.66

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

> m1=aov(count~treatment,data=crd)
> summary(m1)
            Df Sum Sq Mean Sq F value    Pr(>F)    
treatment    3 32.873  10.958  94.584 1.376e-06 ***
Residuals    8  0.927   0.116                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

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


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

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

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

จากค่า p เราสามารถดูได้ว่าคู่ของสิ่งทดลองคู่ใดมีความแตกต่างกันบ้าง

Reference


Back to home page