mean comparison สำหรับ RCBD

     โพสต์ก่อนนี้ได้กล่าวถึงการทำ mean comparison สำหรับการทดลองแบบ CRD คราวนี้ลองมาดูกรณีของการทดลองแบบ RCBD บ้าง จากตัวอย่างก่อนหน้านี้ เราพบว่า treatment มีความแตกต่างกันอย่างนัยสำคัญเราจะใช้แพกเกจ  multcomp ในการเปรียบเทียบค่าเฉลี่ยของ treatment ดูว่า Treatment ใดแตกต่างกันบ้าง

   บน R console
เพื่อให้การใช้การเปรียบเทียบ treatment กับ control ทำได้สะดวกเราจะเปลี่ยนวิธีการกรอกข้อมูลเล็กน้อย โดยใช้ 0 แทน control

> data
   trt block nitrogen
1    0     1    34.98
2    0     2    41.22
3    0     3    36.94
4    0     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

 ทำหให้ทั้ง trt และ block เป็น factor
> data$block<-as.factor(data$block)
> data$trt<-as.factor(data$trt)


 > str(data)
'data.frame':   24 obs. of  3 variables:
 $ trt     : Factor w/ 6 levels "0","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
 $ 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 ...

> 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

> library(multcomp)
Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data

คำส่งจะคล้ายการใช้ Tukey เพียงแต่ใช้ "Dunnett"

> compare<-glht(model,linfct=mcp(trt="Dunnett"))
> summary(compare)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = nitrogen ~ trt + block, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)  
2 - 0 == 0    5.755      1.897   3.033  0.03331 *
3 - 0 == 0    8.492      1.897   4.476  0.00181 **
4 - 0 == 0    2.337      1.897   1.232  0.62267  
5 - 0 == 0    1.232      1.897   0.650  0.94550  
6 - 0 == 0    4.947      1.897   2.607  0.07415 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


ในกรณีที่ 0 ไม่ใช่ตัว control ที่เราอยากเปรียบเทียบด้วย เช่นเป็น treatment 4 ตามตามตัวอย่างใน Kuehl (2000) เราสามารถกำหนดให้เปรียบเทียบ treatment ต่างๆ กับ 4 ได้โดยการสร้าง contrast matrix ดังนี้

> contrast <-rbind("0-4"=c(1,0,0,-1,0,0),
                   "2-4"=c(0,1,0,-1,0,0),
                   "3-4"=c(0,0,1,-1,0,0),
                   "5-4"=c(0,0,0,-1,1,0),
                   "6-4"=c(0,0,0,-1,0,1))
> contrast
    [,1] [,2] [,3] [,4] [,5] [,6]
0-4    1    0    0   -1    0    0
2-4    0    1    0   -1    0    0
3-4    0    0    1   -1    0    0
5-4    0    0    0   -1    1    0
6-4    0    0    0   -1    0    1

> summary(glht(model,linfct=mcp(trt=contrast)))

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: aov(formula = nitrogen ~ trt + block, data = data)

Linear Hypotheses:
         Estimate Std. Error t value Pr(>|t|) 
0-4 == 0   -2.337      1.897  -1.232   0.6226 
2-4 == 0    3.417      1.897   1.801   0.2941 
3-4 == 0    6.155      1.897   3.244   0.0218 *
5-4 == 0   -1.105      1.897  -0.582   0.9643 
6-4 == 0    2.610      1.897   1.376   0.5284 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

ทำให้สรุปได้ว่ามีเพียง treatment 3 เท่านั้นที่ให้ค่าแตกต่างกับ treatment 4

วิเคราะห์ 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.