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)

วิเคราะห์ ANOVA สำหรับ RCBD ด้วย R

     ในโพสต์ก่อนหน้านี้ผมได้เขียนถึงการวิเคราะห์ ANOVA สำหรับแผนการทดลองแบบ CRD ไปแล้ว คราวนี้จะได้พูดถึงแผนการทดลองอีกแบบที่พบบ่อยคือ RCBD (Randomized Complete Block Design) โดยผมใช้ตัวอย่างจาก Kuehl (2001) เช่นกัน โดยการทดลองนี้มี 6 treatment และจัดเป็น 4  block

    บน R console
> file<-file.choose()   #เลือกไฟล์ csv ที่บันทึกข้อมูลไว้
> data <-read.csv(file,header=TRUE)    #อ่านไฟล์ที่เลือก
> data
       trt block nitrogen
1  control     1    34.98
2  control     2    41.22
3  control     3    36.94
4  control     4    39.97
5        2     1    40.89
6        2     2    46.69
7        2     3    46.65
8        2     4    41.90
9        3     1    42.07
10       3     2    49.42
11       3     3    52.68
12       3     4    42.91
13       4     1    37.18
14       4     2    45.85
15       4     3    40.23
16       4     4    39.20
17       5     1    37.99
18       5     2    41.99
19       5     3    37.61
20       5     4    40.45
21       6     1    34.89
22       6     2    50.15
23       6     3    44.57
24       6     4    43.29

ดูโครงสร้างข้อมูล

> str(data)
'data.frame':   24 obs. of  3 variables:
 $ trt     : Factor w/ 6 levels "2","3","4","5",..: 6 6 6 6 1 1 1 1 2 2 ...
 $ block   : int  1 2 3 4 1 2 3 4 1 2 ...
 $ nitrogen: num  35 41.2 36.9 40 40.9 ...


จะเห็นว่าข้อมูลที่อ่านเข้ามาอยู่ในรูปของ data.frame โดย trt (หมายถึง treatment ในการทดลองนี้) เป็นข้อมูลแบบจัดกลุ่ม (factor) ที่มี 6 ระดับอยู่แล้ว แต่ว่าข้อมูลของ block ยังไม่เป็น เราจะแปลงข้อมูลนี้เป็น factor

> data$block<-as.factor(data$block)
> str(data)
'data.frame':   24 obs. of  3 variables:
 $ trt     : Factor w/ 6 levels "2","3","4","5",..: 6 6 6 6 1 1 1 1 2 2 ...
 $ block   : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
 $ nitrogen: num  35 41.2 36.9 40 40.9 ...

จะเห็นว่าข้อมูล block กลายเป็น factor ที่มี 4 ระดับ

วิเคราะห์ ANOVA โดยใช้ฟังก์ชั่น aov คล้ายกับในกรณี CRD แต่เราเพิ่มเทอมของ block เข้าในโมเดลดังนี้

> model <-aov(nitrogen~trt+block,data=data)

แสดงผล
> summary(model)
            Df Sum Sq Mean Sq F value  Pr(>F)  
trt          5  201.3   40.26   5.592 0.00419 **
block        3  197.0   65.67   9.120 0.00112 **
Residuals   15  108.0    7.20                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
    จากผลการวิเคราะห์แสดงให้เห็นว่า treatment มีความแตกต่างกันอย่างมีนัยสำคัญ (p = 0.00419)


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

การใช้ 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.