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