โพสต์ก่อนนี้ได้กล่าวถึงการทำ 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
ไม่มีความคิดเห็น:
แสดงความคิดเห็น