Comparing Means: ANOVA
Ekaphan Kraichak
บทเรียนรู้เราจะเรียนรู้การเปรียบเทียบค่าเฉลี่ยของมากกว่า 2 กลุ่มว่าต่างกันไหม โดยใช้การวิเคราะห์ที่เรียกว่า ANOVA หรือ F-test ซึ่งมีขั้นตอนต่าง ๆ ดังต่อไปนี้
1. เรียก ‘tidyverse’ package
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors2. นำเข้าข้อมูล
rice <- read_csv(file = "rice_yield.csv")## Rows: 71 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): province, nitrogen
## dbl (4): rep, yield, soil_moisture, plant_height
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.rice## # A tibble: 71 × 6
## province nitrogen rep yield soil_moisture plant_height
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Lopburi High 1 429 68 92
## 2 Lopburi High 3 442 71 94
## 3 Lopburi High 4 552 75 103
## 4 Lopburi High 5 517 73 98
## 5 Lopburi High 6 446 69 95
## 6 Lopburi High 7 442 70 93
## 7 Lopburi High 8 544 74 101
## 8 Lopburi High 9 576 77 105
## 9 Lopburi High 10 561 76 104
## 10 Lopburi High 11 482 72 97
## # ℹ 61 more rowsข้อมูลนี้เป็นข้อมูลผลผลิตของข้าวจาก 3 จังหวัด ที่ 2 ระดับความเข้มข้นของปุ๋ยไนโตรเจน โดยความหมายของแต่ละ column เป็นดังนี้
province= ชื่อจังหวัดที่มีแปลงปลูก จำนวน 3 จังหวัดnitrogen= ระดับปุ๋ยไนโตรเจนที่ใส่ (high และ low)rep= หมายเลข replicateyield= ปริมาณผลผลิต หน่วยเป็น กก. ต่อ ไร่soil moisture= ความชื้นในดิน หน่วยเป็น ร้อยละplant_height= ความสูงขอต้นข้าว หน่วยเป็น เซนติเมตร
3. ผลของจังหวัดต่อผลผลิต (3 จังหวัด - ANOVA)
1) Visualize data
ก่อนเริ่มวิเคราะห์ควรวาดกราฟก่อนเสมอ
1.1) boxplots in ggplot
rice %>%
ggplot(aes(x = province, y = yield)) +
geom_boxplot()
ภาพที่ 1: กราฟแท่งเทียนแสดงการกระจายของผลผลิต (yield) ข้าวที่ปลูกในพื้นที่สามจังหวัด
1.2) dynamite plots with ggpubr
library(ggpubr)
## add mean ± SE
ggbarplot(rice, x = "province", y = "yield", fill = "province", add = "mean_se") +
scale_fill_brewer(palette = "Accent")
ภาพที่ 1: กราฟแท่งแสดงผลผลิตเฉลี่ย (yield) ของข้าวที่ปลูกในการพื้นที่สามจังหวัด Error bars แสดงถึงค่าเฉลี่ย (mean) ± ส่วนเบี่ยงเบนมาตรฐาน
2) Testing
2.1) Summary Statistics
ก่อนเริ่มก็ควรทำสถิติสรุปก่อนเช่นเคย
rice %>%
group_by(province) %>%
summarize(mean = mean(yield),
sd = sd(yield))## # A tibble: 3 × 3
## province mean sd
## <chr> <dbl> <dbl>
## 1 Buriram 377 68.3
## 2 Lopburi 413. 102.
## 3 Surin 722. 73.52.1) ANOVA (Overall Test)
การวิเคราะห์ ANOVA เป็นการทดสอบโดยรวมว่ามีความแตกต่างระหว่างกลุ่มไหม แต่บอกไม่ได้ว่าใครต่างจากใคร การตีความผลให้ดูที่ค่า p-value เป็นสำคัญ
P ≤ 0.05 = ค่าเฉลี่ยแตกต่างกันระหว่างกลุ่มอย่างมีนัยสำคัญ
P > 0.05 = ไม่สามารถสรุปได้ว่าค่าเฉลี่ยแตกต่างกันระหว่างกลุ่มอย่างมีนัยสำคัญ
province.anova <- aov(yield ~ province, data = rice)
summary(province.anova)## Df Sum Sq Mean Sq F value Pr(>F)
## province 2 1722042 861021 127.8 <2e-16 ***
## Residuals 68 458302 6740
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.2) Post-Hoc Comparison
ถ้า P ≤ 0.05 แล้ว แสดงว่ามีอย่างน้อย 1 คู่ที่ต่างกัน เราสามารถตรวจสอบได้ว่ากลุ่มไหนต่างจากกลุ่มไหนจริง ๆ ด้วย post-hoc test ซึ่งมีหลายทางเลือก
1) Tukey Honest Significant Difference (HSD)
วิธีการนี้เป็นวิธีมาตรฐานในงานสายชีววิทยาและวิทยาศาสตร์สุขภาพ มีความ “โหด” ในการตัดสินใจว่าคู่ใดต่างกันจริง ๆ บ้าง
ใช้ function หลักใน R
TukeyHSD(province.anova)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = yield ~ province, data = rice)
##
## $province
## diff lwr upr p adj
## Lopburi-Buriram 35.86957 -21.52926 93.26839 0.2985767
## Surin-Buriram 345.33333 288.54842 402.11825 0.0000000
## Surin-Lopburi 309.46377 252.06494 366.86259 0.0000000function ที่ง่ายขึ้นใน package agricolae
library(agricolae)
province.HSD <- HSD.test(y = province.anova, trt = "province")
province.HSD## $statistics
## MSerror Df Mean CV
## 6739.734 68 505.3521 16.24528
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey province 3 3.388576 0.05
##
## $means
## yield std r se Min Max Q25 Q50 Q75
## Buriram 377.0000 68.30558 24 16.75775 278 493 321.75 362.5 437.25
## Lopburi 412.8696 101.51458 23 17.11818 266 576 323.00 381.0 499.50
## Surin 722.3333 73.50757 24 16.75775 566 859 683.50 742.5 761.00
##
## $comparison
## NULL
##
## $groups
## yield groups
## Surin 722.3333 a
## Lopburi 412.8696 b
## Buriram 377.0000 b
##
## attr(,"class")
## [1] "group"2) Least Significant Difference (LSD)
นิยมใช้ในสาขาการเกษตร เนื่องจากต้องการตรวจสอบความแตกต่างที่น้อยที่สุดที่ยังนับว่าต่าง
## function in agricolae
province.LSD <- LSD.test(y = province.anova, trt = "province")
province.LSD## $statistics
## MSerror Df Mean CV
## 6739.734 68 505.3521 16.24528
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none province 3 0.05
##
## $means
## yield std r se LCL UCL Min Max Q25 Q50
## Buriram 377.0000 68.30558 24 16.75775 343.5604 410.4396 278 493 321.75 362.5
## Lopburi 412.8696 101.51458 23 17.11818 378.7108 447.0284 266 576 323.00 381.0
## Surin 722.3333 73.50757 24 16.75775 688.8938 755.7729 566 859 683.50 742.5
## Q75
## Buriram 437.25
## Lopburi 499.50
## Surin 761.00
##
## $comparison
## NULL
##
## $groups
## yield groups
## Surin 722.3333 a
## Lopburi 412.8696 b
## Buriram 377.0000 b
##
## attr(,"class")
## [1] "group"3) Duncan Multiple Range Test (DMRT)
อันนี้ก็นิยมใช้ในทางการเกษตร เช่นกัน เหมาะสำหรับข้อมูลที่กระจายไม่ปกติ การคำนวณใช้ลำดับของข้อมูลมาตัดสินใจแทนก็ใช้การกระจาย
province.duncan <- duncan.test(y = province.anova, trt = "province")
province.duncan## $statistics
## MSerror Df Mean CV
## 6739.734 68 505.3521 16.24528
##
## $parameters
## test name.t ntr alpha
## Duncan province 3 0.05
##
## $duncan
## NULL
##
## $means
## yield std r se Min Max Q25 Q50 Q75
## Buriram 377.0000 68.30558 24 16.75775 278 493 321.75 362.5 437.25
## Lopburi 412.8696 101.51458 23 17.11818 266 576 323.00 381.0 499.50
## Surin 722.3333 73.50757 24 16.75775 566 859 683.50 742.5 761.00
##
## $comparison
## NULL
##
## $groups
## yield groups
## Surin 722.3333 a
## Lopburi 412.8696 b
## Buriram 377.0000 b
##
## attr(,"class")
## [1] "group"2.3) Visualize Posthoc
เราสามารถนำผล post-hoc comparison มาวาดกราฟได้ดังต่อไปนี้
1) Built-in agricolae way
plot(province.HSD, main = "Comparison between provinces",
ylim = c(200,1000),
ylab = "yield",
xlab = "province")
2) ggplot2 style
group.labels <- province.HSD$groups %>%
rownames_to_column("province")
rice %>%
ggplot(aes(x = province, y = yield)) +
geom_boxplot() +
geom_text(data = group.labels, aes(x = province, label = groups), y = max(rice$yield)*1.02)
3) ggpubr style
my.pairs <- list(c("Buriram","Lopburi"), c("Surin", "Buriram"))
ggbarplot(rice, x = "province", y = "yield", fill = "province", add = "mean_se") +
stat_compare_means(comparisons = my.pairs, method = "t.test") +
stat_compare_means(method = "anova", label.y = 1200)
วิธีการเขียนบรรยายผล ANOVA และ posthoc รวมกันได้ดังนี้
ผลผลิตข้าวจากการปลูกข้าวในพื้นที่ 3 จังหวัดมีความแตกต่างกันอย่างมีนัยสำคัญ (ANOVA, p < 0.001) โดยที่ผลผลิตข้าวในจังหวัดสุรินทร์มีค่าเฉลี่ย ( 722 ± 73 กก./ไร่) สูงกว่าผลผลิตจังหวัดลพบุรี และ จังหวัดบุรีรัมย์อย่างมีนัยสำคัญ (Tukey HSD test, P < 0.05)