Two-way ANOVA and Analysis with Block
Two-way ANOVA
1) import library
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.6
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.1 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.2.0
## ── 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 errorslibrary(agricolae)2) import data
alfafa <- read_csv("alfalfa_rain.csv")## Rows: 36 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): cultivars, treatment
## dbl (6): AGB, BGB, PRL, LRL, MDA, Pro
##
## ℹ 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. alfafa## # A tibble: 36 × 8
## cultivars treatment AGB BGB PRL LRL MDA Pro
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 C1 D2 1.75 0.73 44.7 123 103. 778
## 2 C1 D2 1.65 0.61 40.9 114 111. 803
## 3 C1 D2 1.74 0.5 39.5 95 97.5 795
## 4 C1 D2 1.48 0.65 40.2 101 103. 884
## 5 C1 D2 1.61 0.87 45.7 108 116. 783
## 6 C1 D2 1.57 0.94 41.8 109 119. 845
## 7 C1 D4 2.1 1.34 43.5 111 100 1443
## 8 C1 D4 1.74 0.83 46.6 89 101 1219
## 9 C1 D4 1.73 0.63 44.2 92 99.9 1269
## 10 C1 D4 2.01 1.18 39.8 90 108. 1315
## # ℹ 26 more rows3) One-way ANOVA
3.1) ดูผลของแต่ละตัวแปร เสมือนว่าอีกตัวแปรไม่มีอยู่เลย (❌ เป็นวิธีที่ไม่ถูก)
อาจจะเป็นวิธีที่ไม่ถูกแต่ใช้ดูผลเบื้องต้น ได้
ผลของ cultivars: Graph
alfafa %>%
ggplot(aes(x = cultivars, y = Pro)) +
geom_boxplot()
ผลของ cultivars: t-test (2 กลุ่ม)
t.test(Pro ~ cultivars, data = alfafa)##
## Welch Two Sample t-test
##
## data: Pro by cultivars
## t = 0.39462, df = 33.967, p-value = 0.6956
## alternative hypothesis: true difference in means between group C1 and group C2 is not equal to 0
## 95 percent confidence interval:
## -677.8502 1004.5168
## sample estimates:
## mean in group C1 mean in group C2
## 1890.167 1726.833ผลของ treatment: ANOVA (3 กลุ่ม)
treatment_anova <- aov(Pro ~ treatment, data = alfafa)
summary(treatment_anova)## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 2 51872202 25936101 1082 <2e-16 ***
## Residuals 33 790699 23961
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1post-hoc test
treatment_hsd <- HSD.test(treatment_anova, trt = "treatment")
treatment_hsd## $statistics
## MSerror Df Mean CV MSD
## 23960.57 33 1808.5 8.559139 155.0641
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey treatment 3 3.470189 0.05
##
## $means
## Pro std r se Min Max Q25 Q50 Q75
## D2 825.4167 41.96851 12 44.68461 778 902 792.50 811.5 851.25
## D4 1101.5000 231.46038 12 44.68461 848 1443 886.25 1088.5 1280.50
## D8 3498.5833 128.63299 12 44.68461 3277 3640 3424.00 3546.5 3591.75
##
## $comparison
## NULL
##
## $groups
## Pro groups
## D8 3498.5833 a
## D4 1101.5000 b
## D2 825.4167 c
##
## attr(,"class")
## [1] "group"3.2) ดูผลของ combo สองตัวแปรพร้อมกัน (✅ วิธีนี้ทำได้)
ถ้าคำถามเราคือ “combo” ไหนของตัวแปรให้ผลดีสุด ให้ใช้วิธีนี้
ขั้นตอนแรก เริ่มจากการสร้างตัวแรปแบบรวมขึ้นมาใหม่ก่อน
alfafa2 <- alfafa %>%
unite("combo", cultivars:treatment)
alfafa2## # A tibble: 36 × 7
## combo AGB BGB PRL LRL MDA Pro
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 C1_D2 1.75 0.73 44.7 123 103. 778
## 2 C1_D2 1.65 0.61 40.9 114 111. 803
## 3 C1_D2 1.74 0.5 39.5 95 97.5 795
## 4 C1_D2 1.48 0.65 40.2 101 103. 884
## 5 C1_D2 1.61 0.87 45.7 108 116. 783
## 6 C1_D2 1.57 0.94 41.8 109 119. 845
## 7 C1_D4 2.1 1.34 43.5 111 100 1443
## 8 C1_D4 1.74 0.83 46.6 89 101 1219
## 9 C1_D4 1.73 0.63 44.2 92 99.9 1269
## 10 C1_D4 2.01 1.18 39.8 90 108. 1315
## # ℹ 26 more rowsขั้นที่ 2 วาดกราฟดูผลประกอบการ
alfafa2 %>%
ggplot(aes(x = combo, y = Pro)) +
geom_boxplot()
reorder graph by y-values
alfafa2 %>%
ggplot(aes(x = reorder(combo, Pro), y = Pro)) +
geom_boxplot() +
labs(x = "cultivar x treatment", y = "Proline level")
ขั้นที่ 3 ทำ ANOVA โดยใช้ตัวแปรใหม่
combo เป็นตัวแปรต้น
combo_anova <- aov(Pro ~ combo, data = alfafa2)
summary(combo_anova)## Df Sum Sq Mean Sq F value Pr(>F)
## combo 5 52436559 10487312 1390 <2e-16 ***
## Residuals 30 226342 7545
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ขั้นที่ 4 ทำ post-hoc analysis ถ้า p ≤ 0.05
combo_hsd <- HSD.test(combo_anova, trt = "combo")
combo_hsd## $statistics
## MSerror Df Mean CV MSD
## 7544.722 30 1808.5 4.802895 152.5325
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey combo 6 4.301464 0.05
##
## $means
## Pro std r se Min Max Q25 Q50 Q75
## C1_D2 814.6667 41.46645 6 35.46059 778 884 786.00 799.0 834.50
## C1_D4 1313.6667 91.84044 6 35.46059 1219 1443 1242.00 1292.0 1381.00
## C1_D8 3542.1667 87.78477 6 35.46059 3429 3640 3472.00 3555.5 3610.50
## C2_D2 836.1667 43.33782 6 35.46059 791 902 799.75 833.5 859.75
## C2_D4 889.3333 37.31845 6 35.46059 848 958 871.25 885.5 890.75
## C2_D8 3455.0000 155.36151 6 35.46059 3277 3631 3312.25 3480.5 3573.75
##
## $comparison
## NULL
##
## $groups
## Pro groups
## C1_D8 3542.1667 a
## C2_D8 3455.0000 a
## C1_D4 1313.6667 b
## C2_D4 889.3333 c
## C2_D2 836.1667 c
## C1_D2 814.6667 c
##
## attr(,"class")
## [1] "group"4) Two-way ANOVA (✅ วิธีนี้ถูกต้องเหมือนกัน)
เป็นกรณีที่เราสนใจว่า แต่ละปัจจัยมีผลอย่างไรกับตัวแปรตามของเรา โดยที่พิจารณาด้วยว่ามีอีกตัวแปรอยู่ด้วยกัน และสามารถบอกได้ด้วยว่า สองปัจจัยอยู่ด้วยกันแล้วมีผลทวีคูณขึ้นด้วยไหม ผลงานร่วมกันแบบนี้เรียกว่า “interaction”
4.1) แต่ละปัจจัยทำงานแยกกัน (ไม่มี interaction)
twoway_anova <- aov(Pro ~ cultivars + treatment, data = alfafa)
summary(twoway_anova)## Df Sum Sq Mean Sq F value Pr(>F)
## cultivars 1 240100 240100 13.95 0.000732 ***
## treatment 2 51872202 25936101 1507.37 < 2e-16 ***
## Residuals 32 550599 17206
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.2) สองปัจจัยทำงานร่วมกัน (มี interaction)
twoway_anova_interaction <- aov(Pro ~ cultivars + treatment + cultivars:treatment, data = alfafa)
summary(twoway_anova_interaction)## Df Sum Sq Mean Sq F value Pr(>F)
## cultivars 1 240100 240100 31.82 3.81e-06 ***
## treatment 2 51872202 25936101 3437.65 < 2e-16 ***
## cultivars:treatment 2 324257 162129 21.49 1.62e-06 ***
## Residuals 30 226342 7545
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.3) Post-hoc comparison จาก two-way ANOVA
ทุกปัจจัยที่ขึ้นว่า P ≤ 0.05 เราจะต้องมาทำ post-hoc comparison ต่อว่ามีผลอย่างไร เราจะไล่ทำไปทีละตัวแปร
1) เริ่มจาก cultivars
twoway_cultivars_hsd <- HSD.test(twoway_anova_interaction, trt = "cultivars")
twoway_cultivars_hsd## $statistics
## MSerror Df Mean CV MSD
## 7544.722 30 1808.5 4.802895 59.13084
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey cultivars 2 2.888209 0.05
##
## $means
## Pro std r se Min Max Q25 Q50 Q75
## C1 1890.167 1222.299 18 20.47318 778 3640 854.75 1292.0 3444.00
## C2 1726.833 1260.825 18 20.47318 791 3631 852.00 889.5 3279.25
##
## $comparison
## NULL
##
## $groups
## Pro groups
## C1 1890.167 a
## C2 1726.833 b
##
## attr(,"class")
## [1] "group"2) ตามด้วย treatment
twoway_treatment_hsd <- HSD.test(twoway_anova_interaction, trt = "treatment")
twoway_treatment_hsd## $statistics
## MSerror Df Mean CV MSD
## 7544.722 30 1808.5 4.802895 87.41998
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey treatment 3 3.48642 0.05
##
## $means
## Pro std r se Min Max Q25 Q50 Q75
## D2 825.4167 41.96851 12 25.07443 778 902 792.50 811.5 851.25
## D4 1101.5000 231.46038 12 25.07443 848 1443 886.25 1088.5 1280.50
## D8 3498.5833 128.63299 12 25.07443 3277 3640 3424.00 3546.5 3591.75
##
## $comparison
## NULL
##
## $groups
## Pro groups
## D8 3498.5833 a
## D4 1101.5000 b
## D2 825.4167 c
##
## attr(,"class")
## [1] "group"3) และดูผลของ interaction
twoway_inter_hsd <- HSD.test(twoway_anova_interaction, trt = c("cultivars", "treatment"))
twoway_inter_hsd## $statistics
## MSerror Df Mean CV MSD
## 7544.722 30 1808.5 4.802895 152.5325
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey cultivars:treatment 6 4.301464 0.05
##
## $means
## Pro std r se Min Max Q25 Q50 Q75
## C1:D2 814.6667 41.46645 6 35.46059 778 884 786.00 799.0 834.50
## C1:D4 1313.6667 91.84044 6 35.46059 1219 1443 1242.00 1292.0 1381.00
## C1:D8 3542.1667 87.78477 6 35.46059 3429 3640 3472.00 3555.5 3610.50
## C2:D2 836.1667 43.33782 6 35.46059 791 902 799.75 833.5 859.75
## C2:D4 889.3333 37.31845 6 35.46059 848 958 871.25 885.5 890.75
## C2:D8 3455.0000 155.36151 6 35.46059 3277 3631 3312.25 3480.5 3573.75
##
## $comparison
## NULL
##
## $groups
## Pro groups
## C1:D8 3542.1667 a
## C2:D8 3455.0000 a
## C1:D4 1313.6667 b
## C2:D4 889.3333 c
## C2:D2 836.1667 c
## C1:D2 814.6667 c
##
## attr(,"class")
## [1] "group"4.4) visualize two-way ANOVA
Boxplot
alfafa %>%
ggplot(aes(x = treatment, y = Pro, color = cultivars)) +
geom_boxplot()
ggbarplot
library(ggpubr)
ggboxplot(alfafa,
x = "cultivars",
y = "Pro",
facet.by = "treatment") +
stat_compare_means(method = "t")