หลังจากการวิเคราะห์ 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)
ไม่มีความคิดเห็น:
แสดงความคิดเห็น