Linear Regression

หากเราต้องการทราบว่าตัวแปรต้น เป็นสาเหตุ ให้เกิดตัวแปรตามจริงหรือไม่เราจะต้องตรวจสอบด้วย linear regression หรือ ชื่อไทยคือสมการถดถอยกันนะครับ เรามาดูตัวอย่างการวิเคราะห์ด้วย code กันนะครับ

ข้อมูลวันนี้เป็นข้อมูลของต้นถั่วแขก (Phaseolus vulgaris) ที่มีขนาดต่าง ๆ กัน (pot_diameter หน่วยเป็น cm) และใส่ปุ๋ยปริมาณไม่เท่ากัน (dose หน่วยเป็น g/L) และวัดความสูง (height หน่วย cm) และจำนวนใบ (leaf_count หน่วย ใบ) ที่เวลาหลังปลูก 4 สัปดาห์

0. Libraries and Data Import

เรียก library

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

นำเข้าข้อมูล

plant <- read_csv("potplant.csv")

## Rows: 40 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (4): dose, pot_diameter, height, leaf_count
## 
## ℹ 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.

plant

## # A tibble: 40 × 4
##     dose pot_diameter height leaf_count
##    <dbl>        <dbl>  <dbl>      <dbl>
##  1     0         10.3   9.28          5
##  2     0         12.1  10.1           5
##  3     0         11.1  11.9           3
##  4     0         11.6  10.1           5
##  5     0          9.3  14.4           3
##  6    10          8.4  21.1           1
##  7    10          8.4  15.7           7
##  8    10         11.5  20.2           3
##  9    10          8.8  17.5           6
## 10    10         10.4  15.1           7
## # ℹ 30 more rows

1. Visualize data first

เราควรเริ่มด้วยการวาดกราฟก่อนเสมอเนอะครับ เราจะเริ่มตรวจสอบสมมติฐานว่า ปริมาณปุ๋ยที่ใส่เข้าไป (dose) มีผลต่อความสูงของต้นไม่หรือไม่

ข้อมูลเราเป็นตัวแปรต่อเนื่อง (continuous) ทั้งคู่ เราจึงต้องใช้ geom_point ในการวาดกราฟครับ

plant %>% 
  ggplot(aes(x = dose, y = height)) +
  geom_point(color = "red") +
  labs(x = "fertilizer dose (g/L)", y = "height (cm)") +
  theme_bw()

เราจะเห็นว่า ยิ่ง dose เพิ่ม เราก็ยิ่งมีความสูงมากขึ้นตามไปด้วยในกรณีนี้

Challenge 1

ลองวาดกราฟเพื่อตรวจสอบสมมติฐานที่ว่า “ขนาดกระถาง (pot_diameter) มีผลต่อความสูงของต้นไม้ (height)” อย่าลืมปรับชื่อแกนให้ถูกต้องด้วยครับ

2. Perform linear regression

จากด้านบน ถ้าเราต้องการทดสอบว่า dose มีผลต่อ height จริง ๆ จะต้องทำ linear regression ด้วยคำสั่ง lm() ซึ่งย่อมาจาก “linear model”

model1 <- lm(height ~ dose, data = plant)
summary(model1)

## 
## Call:
## lm(formula = height ~ dose, data = plant)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.1085  -2.3343  -0.2562   2.8323   8.4211 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.04834    1.10463   10.00  3.4e-12 ***
## dose         0.77402    0.02641   29.31  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.827 on 38 degrees of freedom
## Multiple R-squared:  0.9576, Adjusted R-squared:  0.9565 
## F-statistic: 859.2 on 1 and 38 DF,  p-value: < 2.2e-16

ตัวเลขหลักที่เราควรสนใจได้แก่

• จุดตัดแกน y (intercept) = 11.04834

• ความชัน (slope) ที่แสดงผลของ dose = 0.774402 ที่แสดงว่า ถ้าเพิ่ม dose 1 g/L แล้ว จะได้ความสูงเพิ่มขึ้น 0.77 cm

และพิจารณาเพิ่มเติมหน่อยว่าค่าเหล่านี้แตกต่างจาก 0 จริง ๆ ไหม โดยดูที่ p-value Pr(>|t|)

• p-value ของ (Intercept) = 3.4e-12 ซึ่ง < 0.05 แสดงว่าต่างจาก 0 อย่างมีนัยสำคัญ

• p-value ของ dose < 2e-16 ซึ่ง < 0.05 แสดงว่าต่างจาก 0 อย่างมีนัยสำคัญ

อีกหนึ่งสิ่งที่ควรตรวจสอบคือ ค่า \(R^2\) ที่บอกว่าโมเดลนี้ทำงานได้ดีหรือไม่ดีในการทำนายความสูง ค่า \(R^2\) ที่เข้าใกล้ 1 แปลว่า โมเดลนี้ทำงานในการทำนายได้ดี ในกรณีนี้ของเรา \(R^2\) = 0.9565 แสดงว่าทำงานได้ดี

ถ้าเราได้ข้อมูลว่าค่า intercept, slope ทำงานได้ดี โมเดลมีความสามารถในการทำนายค่อนข้างสูง เราสามารถนำข้อมูลมาเขียนเป็นสมการเพื่อทำนายค่าความสูงจากปริมาณปุ๋ยได้ดังนี้

\[ Height = slope\times Dose + intercept \\ Height = 0.77402\times Dose + 11.04834 \]

และสามารถเขียนรายงานได้ลักษณะนี้

ความสูงของต้นถั่วแขกเพิ่มขึ้นตามปริมาณปุ๋ยที่ใส่ไปอย่างมีนัยสำคัญ (linear regression, \(R^2\) = 0.96, P < 0.001)

Challenge 2

ลองทำ linear regression เพื่อตรวจสอบสมมติฐานว่า “ขนาดกระถาง (pot_diameter) มีผลต่อความสูงของต้นไม้ (height)”

เขียนรายงานผลต่อไปนี้ด้วย

• intercept = _________ ค่า p-value ของ intercept = _________

• slope = _________ ค่า p-value ของ slope = _________

• ค่า \(R^2\) = ________

• โดยรวมแล้วคิดว่า ขนาดกระถางมีผลต่อความสูงของต้นถั่วแขกไหมครับ

3. Add regression line

ถ้าเราต้องการแสดงเส้นสมการที่ได้จาก regression line เราสามารถทำได้โดยการเพิ่ม geom_smooth() เข้าไปในกราฟได้

plant %>% 
  ggplot(aes(x = dose, y = height)) +
  geom_point(color = "red") +
  geom_smooth(method = "lm", se = F, color = "black") +
  labs(x = "fertilizer dose (g/L)", y = "height (cm)") +
  theme_bw()

## `geom_smooth()` using formula = 'y ~ x'

ภาพที่ 2: กราฟจุดแสดงผลของปริมาณความเข้มข้นของปุ๋ยต่อความสูงของต้นถั่วแขก (Phaseolus vulgaris) หลังจากการปลูกได้ 4 สัปดาห์ เส้นถึบสีดำแสดงสมการถดถอย (slope = 0.77, \(R^2\) = 0.96)

Challenge 3

ลองตอบคำถามด้วยตัวเองว่า อะไรมีผลต่อจำนวนใบ (leaf_count) บ้าง โดยวาดกราฟและตรวจสอบสมมติฐานด้วยการวาดกราฟที่มีเส้น linear regression ด้วย