Comparing Means: t-test

บทเรียนรู้เราจะเรียนรู้การเปรียบเทียบค่าเฉลี่ยของสองกลุ่มว่าต่างกันไหม โดยใช้การวิเคราะห์ที่เรียกว่า t-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 errors

2. อ่านข้อมูลจากไฟล์

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 = หมายเลข replicate

• yield = ปริมาณผลผลิต หน่วยเป็น กก. ต่อ ไร่

• soil moisture = ความชื้นในดิน หน่วยเป็น ร้อยละ

• plant_height = ความสูงขอต้นข้าว หน่วยเป็น เซนติเมตร

3. ผลของ Nitrogen (2 กลุ่ม เปรีบบเทียบด้วย t-test)

1) Visualize data

ก่อนที่จะเริ่มทดสอบสถิติ อยากให้วาดกราฟดูก่อน

1.1) Boxplots in ggplot

rice %>% 
  ggplot(aes(x = nitrogen, y = yield)) +
  geom_boxplot()

ตัวอย่างคำบรรยายใต้ภาพ (capture)

ภาพที่ 1: กราฟแท่งเทียนแสดงการกระจายของผลผลิต (yield) ข้าวที่ปลูกในการทดลองเพิ่มไนโตรเจน 2 ระดับ (High และ Low)

1.2) dynamite plots with ggpubr

library(ggpubr)

## add mean ±SD (standard Deviation)
ggbarplot(data = rice, 
          x = "nitrogen", 
          y = "yield", 
          fill = "lightblue", 
          add = "mean_sd") +
  labs(x = "Nitrogen Level", y = "Yield (kg/rai)") 

ตัวอย่างคำบรรยายใต้ภาพ (capture)

ภาพที่ 1: กราฟแท่งแสดงผลผลิตเฉลี่ย (yield) ของข้าวที่ปลูกในการทดลองเพิ่มไนโตรเจน 2 ระดับ (High และ Low) Error bars แสดงถึงค่าเฉลี่ย (mean) ± ส่วนเบี่ยงเบนมาตรฐาน

## add mean ± SE (standard error)
ggbarplot(data = rice, 
          x = "nitrogen", 
          y = "yield", 
          fill = "nitrogen", 
          add = "mean_se")

2) Testing

2.1) summary statistics

ก่อนเริ่มทำ t-test ควรคำนวณค่า means และ standard deviation ไว้ด้วย

rice %>% 
  group_by(nitrogen) %>% 
  summarize(mean = mean(yield), 
            sd = sd(yield))

## # A tibble: 2 × 3
##   nitrogen  mean    sd
##   <chr>    <dbl> <dbl>
## 1 High      550.  129.
## 2 Low       462.  205.

2.2) t-test

ผลการทดสอบให้ดูที่ค่า p-value เป็นสำคัญ

• P ≤ 0.05 = ค่าเฉลี่ยของสองกลุ่มแตกต่างกันอย่างมีนัยสำคัญ

• P > 0.05 = ไม่สามารถสรุปได้ว่าค่าเฉลี่ยแตกต่างกันอย่างมีนัยสำคัญ

t.test(yield ~ nitrogen, data = rice)

## 
##  Welch Two Sample t-test
## 
## data:  yield by nitrogen
## t = 2.1456, df = 59.199, p-value = 0.03602
## alternative hypothesis: true difference in means between group High and group Low is not equal to 0
## 95 percent confidence interval:
##    5.879608 168.428329
## sample estimates:
## mean in group High  mean in group Low 
##           549.5429           462.3889

ตัวอย่างการเขียนรายงานในเล่ม ควรมี descriptive statistics ควบคู่กันไปด้วย

ระดับไนโตรเจนสูงส่งผลให้ได้ผลผลิตข้าวเฉลี่ย 550 ± 129 กก./ไร่ (ค่าเฉลี่ย ± ส่วนเบี่ยงเบนมาตรฐาน) ซึ่งสูงกว่าระดับไนโตรเจนต่ำที่ให้ผลผลิตเฉลี่ย 462 ± 205 กก./ไร่ อย่างมีนัยสำคัญทางสถิติ (t-test, p = 0.04)

3.2) วาดกราฟพร้อมรายงานค่าสถิติ

This way you can get both statistics and graph on the same place.

ggbarplot(data = rice, 
          x = "nitrogen", 
          y = "yield", 
          fill = "pink", 
          add = "mean_sd") + 
  stat_compare_means(
    comparison = list(c("High","Low")), 
    method = "t.test")