Category: R

สอนปลูกต้นไม้ในภาษา R (ภาค 2): วิธีสร้าง ประเมิน และปรับทูน random forest model ด้วย ranger package – ตัวอย่างการทำนายระดับการกินน้ำมันของรถใน mpg dataset

ในบทความนี้ เราจะไปทำความรู้จักกับ random forest รวมทั้งการสร้าง ประเมิน และปรับทูน random forest model ด้วย ranger package ในภาษา R

ถ้าพร้อมแล้ว ไปเริ่มกันเลย

🌲 Random Forest Model คืออะไร?

Random forest model เป็น tree-based model ซึ่งสุ่มสร้าง decision trees ขึ้นมาหลาย ๆ ต้น (forest) และใช้ผลลัพธ์ในภาพรวมเพื่อทำนายข้อมูลสุดท้าย:

Regression task: หาค่าเฉลี่ยของผลลัพธ์จากทุกต้น
Classification task: ดูผลลัพธ์ที่เป็นเสียงโหวตข้างมาก

Random forest เป็น model ที่ทรงพลัง เพราะใช้ผลรวมของหลาย ๆ decision trees แม้ว่า decision tree แต่ละต้นจะมีความสามารถในการทำนายนอยก็ตาม

💻 Random Forest Models ในภาษา R

ในภาษา R เรามี 2 packages ที่นิยมใช้สร้าง random forest model ได้แก่:

randomForest ซึ่งเป็น package ที่มีลูกเล่น แต่เก่ากว่า
ranger ซึ่งใหม่กว่า ประมวลผลได้เร็วกว่า และใช้งานง่ายกว่า

ในบทความก่อน เราดูวิธีการใช้ randomForest แล้ว

ในบทความนี้ เราจะไปดูวิธีใช้ ranger โดยใช้ mpg dataset เป็นตัวอย่างกัน

🚗 mpg Dataset

mpg dataset เป็น dataset จาก ggplots2 package และมีข้อมูลของรถ 38 รุ่น จากช่วงปี ช่วง ค.ศ. 1999 ถึง 2008 ทั้งหมด 11 columns ดังนี้:

No.	Column	Description
1	`manufacturer`	ผู้ผลิต
2	`model`	รุ่นรถ
3	`displ`	ขนาดถังน้ำมัน (ลิตร)
4	`year`	ปีที่ผลิต
5	`cyl`	จำนวนลูกสูบ
6	`trans`	ประเภทเกียร์
7	`drv`	ประเภทล้อขับเคลื่อน
8	`cty`	ระดับการกินน้ำมันเวลาวิ่งในเมือง
9	`hwy`	ระดับการกินน้ำมันเวลาวิ่งบน highway
10	`fl`	ประเภทน้ำมัน
11	`class`	ประเภทรถ

ในบทความนี้ เราจะลองใช้ ranger เพื่อทำนาย hwy กัน

เราสามารถเตรียม mpg เพื่อสร้าง random forest model ได้ดังนี้

โหลด dataset:

# Install ggplot2
install.packages("ggplot2")

# Load ggplot2
library(ggplot2)

# Load the dataset
data(mpg)

ดูตัวอย่าง dataset:

# Preview
head(mpg)

ผลลัพธ์:

# A tibble: 6 × 11
  manufacturer model displ  year   cyl trans      drv     cty   hwy fl    class  
  <chr>        <chr> <dbl> <int> <int> <chr>      <chr> <int> <int> <chr> <chr>  
1 audi         a4      1.8  1999     4 auto(l5)   f        18    29 p     compact
2 audi         a4      1.8  1999     4 manual(m5) f        21    29 p     compact
3 audi         a4      2    2008     4 manual(m6) f        20    31 p     compact
4 audi         a4      2    2008     4 auto(av)   f        21    30 p     compact
5 audi         a4      2.8  1999     6 auto(l5)   f        16    26 p     compact
6 audi         a4      2.8  1999     6 manual(m5) f        18    26 p     compact

สำรวจโครงสร้าง:

# View the structure
str(mpg)

ผลลัพธ์:

tibble [234 × 11] (S3: tbl_df/tbl/data.frame)
 $ manufacturer: chr [1:234] "audi" "audi" "audi" "audi" ...
 $ model       : chr [1:234] "a4" "a4" "a4" "a4" ...
 $ displ       : num [1:234] 1.8 1.8 2 2 2.8 2.8 3.1 1.8 1.8 2 ...
 $ year        : int [1:234] 1999 1999 2008 2008 1999 1999 2008 1999 1999 2008 ...
 $ cyl         : int [1:234] 4 4 4 4 6 6 6 4 4 4 ...
 $ trans       : chr [1:234] "auto(l5)" "manual(m5)" "manual(m6)" "auto(av)" ...
 $ drv         : chr [1:234] "f" "f" "f" "f" ...
 $ cty         : int [1:234] 18 21 20 21 16 18 18 18 16 20 ...
 $ hwy         : int [1:234] 29 29 31 30 26 26 27 26 25 28 ...
 $ fl          : chr [1:234] "p" "p" "p" "p" ...
 $ class       : chr [1:234] "compact" "compact" "compact" "compact" ...

จากผลลัพธ์จะเห็นว่า บาง columns (เช่น manufacturer, model) มีข้อมูลประเภท character ซึ่งเราควระเปลี่ยนเป็น factor เพื่อช่วยให้การสร้าง model มีประสิทธิภาพมากขึ้น:

# Convert character columns to factor

## Get character columns
chr_cols <- c("manufacturer", "model",
              "trans", "drv",
              "fl", "class")

## For-loop through the character columns
for (col in chr_cols) {
  mpg[[col]] <- as.factor(mpg[[col]])
}

## Check the results
str(mpg)

ผลลัพธ์:

tibble [234 × 11] (S3: tbl_df/tbl/data.frame)
 $ manufacturer: Factor w/ 15 levels "audi","chevrolet",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ model       : Factor w/ 38 levels "4runner 4wd",..: 2 2 2 2 2 2 2 3 3 3 ...
 $ displ       : num [1:234] 1.8 1.8 2 2 2.8 2.8 3.1 1.8 1.8 2 ...
 $ year        : int [1:234] 1999 1999 2008 2008 1999 1999 2008 1999 1999 2008 ...
 $ cyl         : int [1:234] 4 4 4 4 6 6 6 4 4 4 ...
 $ trans       : Factor w/ 10 levels "auto(av)","auto(l3)",..: 4 9 10 1 4 9 1 9 4 10 ...
 $ drv         : Factor w/ 3 levels "4","f","r": 2 2 2 2 2 2 2 1 1 1 ...
 $ cty         : int [1:234] 18 21 20 21 16 18 18 18 16 20 ...
 $ hwy         : int [1:234] 29 29 31 30 26 26 27 26 25 28 ...
 $ fl          : Factor w/ 5 levels "c","d","e","p",..: 4 4 4 4 4 4 4 4 4 4 ...
 $ class       : Factor w/ 7 levels "2seater","compact",..: 2 2 2 2 2 2 2 2 2 2 ...

ตอนนี้ เราพร้อมที่จะนำ dataset ไปใช้งานกับ ranger แล้ว

🐣 ranger Basics

การใช้งาน ranger มีอยู่ 4 ขั้นตอน:

ติดตั้งและโหลด ranger
สร้าง training และ test sets
สร้าง random forest model
ทดสอบความสามารถของ model

1️⃣ ติดตั้งและโหลด ranger

ในครั้งแรกสุด ให้เราติดตั้ง ranger ด้วยคำสั่ง install.packages():

# Install
install.packages("ranger")

และทุกครั้งที่เราต้องการใช้งาน ranger ให้เราเรียกใช้งานด้วย library():

# Load
library(ranger)

2️⃣ สร้าง Training และ Test Sets

ในขั้นที่ 2 เราจะแบ่ง dataset เป็น 2 ส่วน ได้แก่:

Training set สำหรับสร้าง model (70% ของ dataset)
Test set สำหรับทดสอบ model (30% ของ dataset)

# Split the data

## Set seed for reproducibility
set.seed(123)

## Get training rows
train_rows <- sample(nrow(mpg),
                     nrow(mpg) * 0.7)

## Create a training set
train <- mpg[train_rows, ]

## Create a test set
test <- mpg[-train_rows, ]

3️⃣ สร้าง Random Forest Model

ในขั้นที่ 3 เราจะสร้าง random forest ด้วย ranger() ซึ่งต้องการ input หลัก 2 อย่าง ดังนี้:

ranger(formula, data)

formula: ระบุตัวแปรต้นและตัวแปรตามที่ใช้ในการสร้าง model
data: dataset ที่ใช้สร้าง model (เราจะใช้ training set กัน)

เราจะเรียกใช้ ranger() ดังนี้:

# Initial random forest model

## Set seed for reproducibility
set.seed(123)

## Train the model
rf_model <- ranger(hwy ~ .,
                   data = train)

Note: เราใช้ set.seed() เพื่อให้เราสามารถสร้าง model ซ้ำได้ เพราะ random forest มีการสร้าง decision trees แบบสุ่ม

เมื่อได้ model มาแล้ว เราสามารถดูรายละเอียดของ model ได้แบบนี้:

# Print the model
rf_model

ผลลัพธ์:

Ranger result

Call:
 ranger(hwy ~ ., data = train) 

Type:                             Regression 
Number of trees:                  500 
Sample size:                      163 
Number of independent variables:  10 
Mtry:                             3 
Target node size:                 5 
Variable importance mode:         none 
Splitrule:                        variance 
OOB prediction error (MSE):       1.682456 
R squared (OOB):                  0.9584596

ในผลลัพธ์ เราจะเห็นลักษณะต่าง ๆ ของ model เช่น ประเภท model (type) และ จำนวน decision trees ที่ถูกสร้างขึ้นมา (sample size)

4️⃣ ทดสอบความสามารถของ Model

สุดท้าย เราจะทดสอบความสามารถของ model ในการทำนายข้อมูล โดยเริ่มจากใช้ model ทำนายข้อมูลใน test set:

# Make predictions
preds <- predict(rf_model,
                 data = test)$predictions

จากนั้น คำนวณตัวบ่งชี้ความสามารถ (metric) ซึ่งสำหรับ regression model มีอยู่ 3 ตัว ได้แก่:

MAE (mean absolute error): ค่าเฉลี่ยความคลาดเคลื่อนแบบสัมบูรณ์ (ยิ่งน้อยยิ่งดี)
RMSE (root mean square error): ค่าเฉลี่ยความคาดเคลื่อนแบบยกกำลังสอง (ยิ่งน้อยยิ่งดี)
R squared: สัดส่วนข้อมูลที่อธิบายได้ด้วย model (ยิ่งมากยิ่งดี)

# Get errors
errors <- test$hwy - preds

# Calculate MAE
mae <- mean(abs(errors))

# Calculate RMSE
rmse <- sqrt(mean(errors^2))

# Calculate R squared
r_sq <- 1 - (sum((errors)^2) / sum((test$hwy - mean(test$hwy))^2))

# Print the results
cat("Initial model MAE:", round(mae, 2), "\n")
cat("Initial model RMSE:", round(rmse, 2), "\n")
cat("Initial model R squared:", round(r_sq, 2), "\n")

ผลลัพธ์:

Initial model MAE: 0.79
Initial model RMSE: 1.07
Initial model R squared: 0.95

⏲️ Hyperparametre Tuning

ranger มี hyperparametre มากมายที่เราสามารถปรับแต่งเพื่อเพิ่มประสิทธิภาพของ random forest model ได้ เช่น:

num.trees: จำนวน decision trees ที่จะสร้าง
mtry: จำนวนตัวแปรต้นที่จะถูกสุ่มไปใช้ในแต่ละ node
min.node.size: จำนวนข้อมูลขั้นต่ำที่แต่ละ node จะต้องมี

เราสามารถใช้ for loop เพื่อปรับหาค่า hyperparametre ที่ดีที่สุดได้ดังนี้:

# Define hyperparametres
ntree_vals <- c(300, 500, 700)
mtry_vals <- 2:5
min_node_vals <- c(1, 5, 10)

# Create a hyperparametre grid
grid <- expand.grid(num.trees = ntree_vals,
                    mtry = mtry_vals,
                    min.node.size = min_node_vals)

# Instantiate an empty data frame
hpt_results <- data.frame()

# For-loop through the hyperparametre grid
for (i in 1:nrow(grid)) {
  
  ## Get the combination
  params <- grid[i, ]
  
  ## Set seed for reproducibility
  set.seed(123)
  
  ## Fit the model
  model <- ranger(hwy ~ .,
                  data = train,
                  num.trees = params$num.trees,
                  mtry = params$mtry,
                  min.node.size = params$min.node.size)
  
  ## Make predictions
  preds <- predict(model,
                   data = test)$predictions
  
  ## Get errors
  errors <- test$hwy - preds
  
  ## Calculate MAE
  mae <- mean(abs(errors))
  
  ## Calculate RMSE
  rmse <- sqrt(mean(errors^2))
  
  ## Store the results
  hpt_results <- rbind(hpt_results,
                       cbind(params,
                             MAE = mae,
                             RMSE = rmse))
}

# View the results
hpt_results

ผลลัพธ์:

   num.trees mtry min.node.size       MAE      RMSE
1        300    2             1 0.8101026 1.0971836
2        500    2             1 0.8012484 1.0973957
3        700    2             1 0.8039271 1.1001252
4        300    3             1 0.7434543 1.0051344
5        500    3             1 0.7417985 1.0069989
6        700    3             1 0.7421666 1.0028184
7        300    4             1 0.6989314 0.9074216
8        500    4             1 0.7130704 0.9314843
9        700    4             1 0.7141147 0.9292718
10       300    5             1 0.7157657 0.9370918
11       500    5             1 0.7131899 0.9266787
12       700    5             1 0.7091556 0.9238312
13       300    2             5 0.8570125 1.1673637
14       500    2             5 0.8515116 1.1736009
15       700    2             5 0.8522571 1.1756648
16       300    3             5 0.7885005 1.0654548
17       500    3             5 0.7872713 1.0664734
18       700    3             5 0.7859149 1.0581331
19       300    4             5 0.7561500 0.9790160
20       500    4             5 0.7623437 0.9869463
21       700    4             5 0.7611660 0.9813048
22       300    5             5 0.7615190 0.9777769
23       500    5             5 0.7615861 0.9804616
24       700    5             5 0.7613151 0.9788333
25       300    2            10 0.9257704 1.2391377
26       500    2            10 0.9292344 1.2611164
27       700    2            10 0.9258555 1.2635794
28       300    3            10 0.8790601 1.1635695
29       500    3            10 0.8704461 1.1594165
30       700    3            10 0.8704562 1.1507016
31       300    4            10 0.8609516 1.0887466
32       500    4            10 0.8672105 1.0962367
33       700    4            10 0.8624934 1.0875710
34       300    5            10 0.8558867 1.0811168
35       500    5            10 0.8567463 1.0783473
36       700    5            10 0.8536824 1.0751511

จะเห็นว่า เราจะได้ MSE และ RMSE ของส่วนผสมระหว่างแต่ละ hyperparametre มา

เราสามารถใช้ ggplot() เพื่อช่วยเลือก hyperparametres ที่ดีที่สุดได้ดังนี้:

# Visualise the results
ggplot(hpt_results,
       aes(x = mtry,
           y = RMSE,
           color = factor(num.trees))) +
  
  ## Use scatter plot
  geom_point(aes(size = min.node.size)) +
  
  ## Set theme to minimal
  theme_minimal() +
  
  ## Add title, labels, and legends
  labs(title = "Hyperparametre Tuning Results",
       x = "mtry",
       y = "RMSE",
       color = "num.trees",
       size = "min.node.size")

ผลลัพธ์:

จากกราฟ จะเห็นได้ว่า hyperparametres ที่ดีที่สุด (มี RMSE น้อยที่สุด) คือ:

num.trees = 300
mtry = 4
min.node.size = 2.5

เมื่อได้ค่า hyperparametres แล้ว เราสามารถใส่ค่าเหล่านี้กลับเข้าไปใน model และทดสอบความสามารถได้เลย

สร้าง model:

# Define the best hyperparametres
best_num.tree <- 300
best_mtry <- 4
best_min.node.size <- 2.5

# Fit the model
rf_model_new <- ranger(hwy ~ .,
                       data = train,
                       num.tree = best_num.tree,
                       mtry = best_mtry,
                       min.node.size = best_min.node.size)

ทดสอบความสามารถ:

# Evaluate the model

## Make predictions
preds_new <- predict(rf_model_new,
                     data = test)$predictions

## Get errors
errors_new <- test$hwy - preds_new

## Calculate MAE
mae_new <- mean(abs(errors_new))

## Calculate RMSE
rmse_new <- sqrt(mean(errors_new^2))

## Calculate R squared
r_sq_new <- 1 - (sum((errors_new)^2) / sum((test$hwy - mean(test$hwy))^2))

## Print the results
cat("Final model MAE:", round(mae_new, 2), "\n")
cat("Final model RMSE:", round(rmse_new, 2), "\n")
cat("Final model R squared:", round(r_sq_new, 2), "\n")

ผลลัพธ์:

Final model MAE: 0.71
Final model RMSE: 0.93
Final model R squared: 0.96

เราสามารถเปรียบเทียบความสามารถของ model ล่าสุด (final model) กับ model ก่อนหน้านี้ (initial model) ได้:

# Compare the two models
model_comp <- data.frame(Model = c("Initial", "Final"),
                         MAE = c(round(mae, 2), round(mae_new, 2)),
                         RMSE = c(round(rmse, 2), round(rmse_new, 2)),
                         R_Squared = c(round(r_sq, 2), round(r_sq_new, 2)))

# Print
model_comp

ผลลัพธ์:

    Model  MAE RMSE R_Squared
1 Initial 0.85 1.08      0.95
2   Final 0.71 0.93      0.96

ซึ่งจะเห็นว่า model ใหม่สามารถทำนายข้อมูลได้ดีขึ้น เพราะมี MAE และ RMSE ที่ลดลง รวมทั้ง R squared ที่เพิ่มขึ้น

🍩 Bonus: Variable Importance

ส่งท้าย ในกรณีที่เราต้องการดูว่า ตัวแปรต้นไหนมีความสำคัญต่อการทำนายมากที่สุด เราสามารถใช้ importance argument ใน ranger() คู่กับ vip() จาก vip package ได้แบบนี้:

# Fit the model with importance
rf_model_new <- ranger(hwy ~ .,
                       data = train,
                       num.tree = best_num.tree,
                       mtry = best_mtry,
                       min.node.size = best_min.node.size,
                       importance = "permutation") # Add importance

# Install vip package
install.packages("vip")

# Load vip package
library(vip)

# Get variabe importance
vip(rf_model_new)  +
  
  ## Add title and labels
  labs(title = "Variable Importance - Final Random Forest Model",
       x = "Variables",
       y = "Importance") +
  
  ## Set theme to minimal
  theme_minimal()

ผลลัพธ์:

จากกราฟ จะเห็นได้ว่า ตัวแปรต้นที่สำคัญที่สุด 3 ตัว ได้แก่:

cty: ระดับการกินน้ำมันเวลาวิ่งในเมือง
displ: ขนาดถังน้ำมัน (ลิตร)
cyl: จำนวนลูกสูบ

😎 Summary

ในบทความนี้ เราได้ดูวิธีการใช้ ranger package เพื่อ:

สร้าง random forest model
ปรับทูน model

พร้อมวิธีการประเมิน model ด้วย predict() และการคำนวณ MAE, RMSE, และ R squared รวมทั้งดูความสำคัญของตัวแปรต้นด้วย vip package

😺 GitHub

ดู code ทั้งหมดในบทความนี้ได้ที่ GitHub

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

🙋 ใครที่กำลังเรียนจิตวิทยาหรือทำงานสายจิตวิทยา และเบื่อที่ต้องใช้ software ราคาแพงอย่าง SPSS และ Excel เพื่อทำข้อมูล

💪 ผมขอแนะนำ R Book for Psychologists หนังสือสอนใช้ภาษา R เพื่อการวิเคราะห์ข้อมูลทางจิตวิทยา ที่เขียนมาเพื่อนักจิตวิทยาที่ไม่เคยมีประสบการณ์เขียน code มาก่อน

ในหนังสือ เราจะปูพื้นฐานภาษา R และพาไปดูวิธีวิเคราะห์สถิติที่ใช้บ่อยกัน เช่น:

Correlation
t-tests
ANOVA
Reliability
Factor analysis

🚀 เมื่ออ่านและทำตามตัวอย่างใน R Book for Psychologists ทุกคนจะไม่ต้องพึง SPSS และ Excel ในการทำงานอีกต่อไป และสามารถวิเคราะห์ข้อมูลด้วยตัวเองได้ด้วยความมั่นใจ

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists

2025-06-05

วิธีสร้าง linear regression ด้วย lm() ในภาษา R — ตัวอย่างการทำนายราคาเพชรใน diamonds dataset

Linear regression เป็นวิธีการทำนายข้อมูลด้วยสมการเส้นตรง:

y = a + bx

y = ตัวแปรตาม หรือข้อมูลที่ต้องการทำนาย
a = จุดตัดระหว่าง x และ y (intercept)
b = ค่าความชัด (slope)
x = ตัวแปรต้น

เนื่องจากเป็นเทคนิคที่ใช้งานและทำความเข้าใจได้ง่าย linear regression จึงเป็นวิธีที่นิยมใช้ในการทำนายข้อมูลในบริบทต่าง ๆ เช่น:

ทำนาย	จาก
กำไร	ค่าโฆษณา
ความสามารถของนักกีฬา	ชั่วโมงฝึกซ้อม
ความดันเลือด	ปริมาณยา + อายุ
ผลลิตทางการเกษตร	ปริมาณน้ำ + ปุ๋ย

ในบทความนี้ เราจะมาดูวิธีใช้ linear regression ในภาษา R กัน

ถ้าพร้อมแล้ว ไปเริ่มกันเลย

💎 Example Dataset: diamonds

ในบทความนี้ เราจะใช้ diamonds dataset เป็นตัวอย่างในการใช้ linear regression กัน

diamonds dataset เป็น built-in dataset จาก ggplot2 package ซึ่งมีข้อมูลเพชรมากกว่า 50,000 ตัวอย่าง และประกอบด้วย 10 columns ดังนี้:

No.	Column	Description
1	`price`	ราคา (ดอลล่าร์สหรัฐฯ)
2	`caret`	น้ำหนัก
3	`cut`	คุณภาพ
4	`color`	สี
5	`clarity`	ความใสของเพชร
6	`x`	ความยาว
7	`y`	ความกว้าง
8	`z`	ความลึก
9	`depth`	สัดส่วนความลึก
10	`table`	สัดส่วนความกว้างของยอดเพชรต่อส่วนที่กว้างที่สุด

เป้าหมายของเรา คือ ทำนายราคาเพชร (price)

⬇️ Load diamonds

ในการใช้งาน diamonds เราสามารถเรียกใช้งาน dataset ได้ดังนี้:

ขั้นที่ 1. ติดตั้งและโหลด ggplot2:

# Install
install.packages("ggplot2")

# Load
library(ggplot2)

ขั้นที่ 2. โหลด diamonds dataset:

# Load dataset
data(diamonds)

ขั้นที่ 3. ดูตัวอย่างข้อมูล 10 rows แรกใน dataset:

# Preview the dataset
head(diamonds, 10)

ผลลัพธ์:

# A tibble: 10 × 10
   carat cut       color clarity depth table price     x     y     z
   <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl>
 1  0.23 Ideal     E     SI2      61.5    55   326  3.95  3.98  2.43
 2  0.21 Premium   E     SI1      59.8    61   326  3.89  3.84  2.31
 3  0.23 Good      E     VS1      56.9    65   327  4.05  4.07  2.31
 4  0.29 Premium   I     VS2      62.4    58   334  4.2   4.23  2.63
 5  0.31 Good      J     SI2      63.3    58   335  4.34  4.35  2.75
 6  0.24 Very Good J     VVS2     62.8    57   336  3.94  3.96  2.48
 7  0.24 Very Good I     VVS1     62.3    57   336  3.95  3.98  2.47
 8  0.26 Very Good H     SI1      61.9    55   337  4.07  4.11  2.53
 9  0.22 Fair      E     VS2      65.1    61   337  3.87  3.78  2.49
10  0.23 Very Good H     VS1      59.4    61   338  4     4.05  2.39

🍳 Prepare the Dataset

ก่อนจะทำนายราคาเพชรด้วย linear regression เราจะเตรียม diamonds dataset ใน 3 ขั้นตอนก่อน ได้แก่:

One-hot encoding
Log transformation
Split data

🪆 Step 1. One-Hot Encoding

ในกรณีที่ตัวแปรต้นที่เป็น categorical เราจะต้องแปลงตัวแปรเหล่านี้ให้เป็น numeric ก่อน ซึ่งเราสามารถทำได้ด้วย one-hot encoding ดังตัวอย่าง:

ก่อน one-hot encoding:

Data	Cut
1	Ideal
2	Good
3	Fair

หลัง one-hot encoding:

Data	Cut_Ideal	Cut_Good	Cut_Fair
1	1	0	0
2	0	1	0
3	0	0	1

ในภาษา R เราสามารถทำ one-hot encoding ได้ด้วย model.matrix() ดังนี้:

# Set option for one-hot encoding
options(contrasts = c("contr.treatment",
                      "contr.treatment"))

# One-hot encode
cat_dum <- model.matrix(~ cut + color + clarity - 1,
                        data = diamonds)

จากนั้น เราจะนำผลลัพธ์ที่ได้ไปรวมกับตัวแปรตามและตัวแปรต้นที่เป็น numeric:

# Combine one-hot-encoded categorical and numeric variables
dm <- cbind(diamonds[, c("carat",
                         "depth",
                         "table",
                         "x",
                         "y",
                         "z")],
            cat_dum,
            price = diamonds$price)

เราสามารถเช็กผลลัพธ์ของ one-hot encoding ได้ด้วย str():

# Check the results
str(dm)

ผลลัพธ์:

'data.frame':	53940 obs. of  25 variables:
 $ carat       : num  0.23 0.21 0.23 0.29 0.31 0.24 0.24 0.26 0.22 0.23 ...
 $ depth       : num  61.5 59.8 56.9 62.4 63.3 62.8 62.3 61.9 65.1 59.4 ...
 $ table       : num  55 61 65 58 58 57 57 55 61 61 ...
 $ x           : num  3.95 3.89 4.05 4.2 4.34 3.94 3.95 4.07 3.87 4 ...
 $ y           : num  3.98 3.84 4.07 4.23 4.35 3.96 3.98 4.11 3.78 4.05 ...
 $ z           : num  2.43 2.31 2.31 2.63 2.75 2.48 2.47 2.53 2.49 2.39 ...
 $ cutFair     : num  0 0 0 0 0 0 0 0 1 0 ...
 $ cutGood     : num  0 0 1 0 1 0 0 0 0 0 ...
 $ cutVery Good: num  0 0 0 0 0 1 1 1 0 1 ...
 $ cutPremium  : num  0 1 0 1 0 0 0 0 0 0 ...
 $ cutIdeal    : num  1 0 0 0 0 0 0 0 0 0 ...
 $ colorE      : num  1 1 1 0 0 0 0 0 1 0 ...
 $ colorF      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ colorG      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ colorH      : num  0 0 0 0 0 0 0 1 0 1 ...
 $ colorI      : num  0 0 0 1 0 0 1 0 0 0 ...
 $ colorJ      : num  0 0 0 0 1 1 0 0 0 0 ...
 $ claritySI2  : num  1 0 0 0 1 0 0 0 0 0 ...
 $ claritySI1  : num  0 1 0 0 0 0 0 1 0 0 ...
 $ clarityVS2  : num  0 0 0 1 0 0 0 0 1 0 ...
 $ clarityVS1  : num  0 0 1 0 0 0 0 0 0 1 ...
 $ clarityVVS2 : num  0 0 0 0 0 1 0 0 0 0 ...
 $ clarityVVS1 : num  0 0 0 0 0 0 1 0 0 0 ...
 $ clarityIF   : num  0 0 0 0 0 0 0 0 0 0 ...
 $ price_log   : num  5.79 5.79 5.79 5.81 5.81 ...

ตอนนี้ ตัวแปรต้นที่เป็น categorical ถูกแปลงเป็น numeric ทั้งหมดแล้ว

📈 Step 2. Log Transformation

ในกรณีที่ตัวแปรตามมีการกระจายตัว (distribution) ไม่ปกติ linear regression ทำนายข้อมูลได้ไม่เต็มประสิทธิภาพนัก

เราสามารถตรวจสอบการกระจายตัวของตัวแปรตามได้ด้วย ggplot():

# Check the distribution of `price`
ggplot(dm,
       aes(x = price)) +
  
  ## Instantiate a histogram
  geom_histogram(binwidth = 100,
                 fill = "skyblue3") +
  
  ## Add text elements
  labs(title = "Distribution of Price",
       x = "Price",
       y = "Count") +
  
  ## Set theme to minimal
  theme_minimal()

ผลลัพธ์:

จากกราฟ เราจะเห็นได้ว่า ตัวแปรตามมีการกระจายตัวแบบเบ้ขวา (right-skewed)

ดังนั้น ก่อนจะใช้ linear regression เราจะต้องแปรตัวแปรตามให้มีการกระจายตัวแบบปกติ (normal distribution) ก่อน ซึ่งเราสามารถทำได้ด้วย log transformation ดังนี้:

# Log-transform `price`
dm$price_log <- log(dm$price)

# Drop `price`
dm$price <- NULL

หลัง log transformation เราสามารถเช็กการกระจายตัวด้วย ggplot() อีกครั้ง:

# Check the distribution of logged `price`
ggplot(dm,
       aes(x = price_log)) +
  
  ## Instantiate a histogram
  geom_histogram(fill = "skyblue3") +
  
  ## Add text elements
  labs(title = "Distribution of Price After Log Transformation",
       x = "Price (Logged)",
       y = "Count") +
  
  ## Set theme to minimal
  theme_minimal()

ผลลัพธ์:

จะเห็นได้ว่า การกระจายตัวของตัวแปรตามใกล้เคียงกับการกระจายตัวแบบปกติมากขึ้นแล้ว

🚄 Step 3. Split the Data

ในขั้นสุดท้ายก่อนใช้ linear regression เราจะแบ่งข้อมูลออกเป็น 2 ชุด:

Training set สำหรับสร้าง linear regression model
Test set สำหรับประเมินความสามารถของ linear regression model

ในบทความนี้ เราจะแบ่ง 80% ของ dataset เป็น training set และ 20% เป็น test set:

# Split the data

## Set seed for reproducibility
set.seed(181)

## Training index
train_index <- sample(nrow(dm),
                      0.8 * nrow(dm))

## Create training set
train_set <- dm[train_index, ]

## Create test set
test_set <- dm[-train_index, ]

ตอนนี้ เราพร้อมที่จะสร้าง linear regression model กันแล้ว

🏷️ Linear Regression Modelling

การสร้าง linear regression model มีอยู่ 3 ขั้นตอน ได้แก่:

Fit the model
Make predictions
Evaluate the model performance

💪 Step 1. Fit the Model

ในขั้นแรก เราจะสร้าง model ด้วย lm() ซึ่งต้องการ input 2 อย่าง:

lm(formula, data)

formula = สูตรการทำนาย โดยเราต้องกำหนดตัวแปรต้นและตัวแปรตาม
data = ชุดข้อมูลที่ใช้สร้าง model

ในการทำนายราคาเพชร เราจะใช้ lm() แบบนี้:

# Fit the model
linear_reg <- lm(price_log ~ .,
                 data = train_set)

อธิบาย code:

price_log ~ . หมายถึง ทำนายราคา (price_log) ด้วยตัวแปรต้นทั้งหมด (.)
data = train_set หมายถึง เรากำหนดชุดข้อมูลที่ใช้เป็น training set

เราสามารถดูข้อมูลของ model ได้ด้วย summary():

# View the model
summary(linear_reg)

ผลลัพธ์:

Call:
lm(formula = price_log ~ ., data = train_set)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2093 -0.0930  0.0019  0.0916  9.8935 

Coefficients: (1 not defined because of singularities)
                 Estimate Std. Error  t value Pr(>|t|)    
(Intercept)    -2.7959573  0.0705854  -39.611  < 2e-16 ***
carat          -0.5270039  0.0086582  -60.867  < 2e-16 ***
depth           0.0512357  0.0008077   63.437  < 2e-16 ***
table           0.0090154  0.0005249   17.175  < 2e-16 ***
x               1.1374016  0.0055578  204.651  < 2e-16 ***
y               0.0290584  0.0031345    9.271  < 2e-16 ***
z               0.0340298  0.0054896    6.199 5.73e-10 ***
cutFair        -0.1528658  0.0060005  -25.476  < 2e-16 ***
cutGood        -0.0639105  0.0036547  -17.487  < 2e-16 ***
`cutVery Good` -0.0313800  0.0025724  -12.199  < 2e-16 ***
cutPremium     -0.0451760  0.0026362  -17.137  < 2e-16 ***
cutIdeal               NA         NA       NA       NA    
colorE         -0.0573940  0.0032281  -17.779  < 2e-16 ***
colorF         -0.0892633  0.0032654  -27.336  < 2e-16 ***
colorG         -0.1573861  0.0032031  -49.136  < 2e-16 ***
colorH         -0.2592763  0.0034037  -76.175  < 2e-16 ***
colorI         -0.3864526  0.0038360 -100.742  < 2e-16 ***
colorJ         -0.5258789  0.0047183 -111.455  < 2e-16 ***
claritySI2      0.4431577  0.0079170   55.976  < 2e-16 ***
claritySI1      0.6087513  0.0078819   77.234  < 2e-16 ***
clarityVS2      0.7523161  0.0079211   94.976  < 2e-16 ***
clarityVS1      0.8200656  0.0080463  101.918  < 2e-16 ***
clarityVVS2     0.9381319  0.0082836  113.252  < 2e-16 ***
clarityVVS1     1.0033931  0.0085098  117.910  < 2e-16 ***
clarityIF       1.0898015  0.0092139  118.277  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1825 on 43128 degrees of freedom
Multiple R-squared:  0.9677,	Adjusted R-squared:  0.9676 
F-statistic: 5.611e+04 on 23 and 43128 DF,  p-value: < 2.2e-16

Note: ดูวิธีการอ่านผลลัพธ์ได้ที่ Explaining the lm() Summary in R และ Understanding Linear Regression Output in R

🔮 Step 2. Make Predictions

ในขั้นที่สอง เราจะใช้ model เพื่อทำนายราคาด้วย predict():

# Predict in the outcome space
pred <- exp(pred_log)

# Preview predictions
head(pred_log)

ผลลัพธ์:

       2        5        9       16       19       22 
5.828071 5.816460 6.111859 5.777434 5.865820 6.088356

จะเห็นว่า ราคาที่ทำนายยังอยู่ในรูป log ซึ่งเราต้องแปลงกลับเป็นราคาปกติด้วย exp():

# Predict in the outcome space
pred <- exp(pred_log)

# Preview predictions
head(pred)

ผลลัพธ์:

       2        5        9       16       19       22 
339.7028 335.7812 451.1766 322.9295 352.7713 440.6961

เราสามารถเปรียบเทียบราคาจริงกับราคาที่ทำนาย พร้อมความคลาดเคลื่อน ได้ดังนี้:

# Compare predictions to actual
results <- data.frame(actual = round(exp(test_set$price_log), 2),
                      predicted = round(pred, 2),
                      diff = round(exp(test_set$price_log) - pred, 2))

# Print results
head(results)

ผลลัพธ์:

   actual predicted    diff
2     326    339.70  -13.70
5     335    335.78   -0.78
9     337    451.18 -114.18
16    345    322.93   22.07
19    351    352.77   -1.77
22    352    440.70  -88.70

🎯 Step 3. Evaluate the Model Performance

ในขั้นสุดท้าย เราจะประเมิน model โดยใช้ 2 ตัวชี้วัด ได้แก่:

Mean absolute error (MAE): ค่าเฉลี่ยความคลาดเคลื่อนโดยสัมบูรณ์
Root mean squared error (RMSE): ค่าเฉลี่ยความคลาดเคลื่อนแบบยกกำลังสอง

ทั้งสองตัวคำนวณความแตกต่างระหว่างสิ่งที่ทำนายและข้อมูลจริง ยิ่ง MAE และ RMSE สูง ก็หมายความว่า การทำนายมีความคาดเคลื่อนมาก แสดงว่า model ทำงานได้ไม่ดีนัก

ในทางกลับกัน ถ้า MAE และ RMSE น้อย ก็แสดงว่า การทำนายใกล้เคียงกับข้อมูลจริง และ model มีความแม่นยำสูง

(Note: เรียนรู้ความแตกต่างระหว่าง MAE และ RMSE ได้ที่ Loss Functions in Machine Learning Explained)

เราสามารถคำนวณ MAE และ RMSE ได้ดังนี้:

# Calculate MAE
mae <- mean(abs(results$diff))

# Calculate RMSE
rmse <- sqrt(mean((results$diff)^2))

# Print the results
cat("MAE:", round(mae, 2), "\n")
cat("RMSE:", round(rmse, 2))

ผลลัพธ์:

MAE: 491.71
RMSE: 1123.68

จากผลลัพธ์ เราจะเห็นว่า โดยเฉลี่ย model ทำนายราคาคลาดเคลื่อนไปประมาณ 492 ดอลล่าร์ (MAE)

😎 Summary

ในบทความนี้ เราได้ดูวิธีการทำ linear regression ในภาษา R กัน

เราดูวิธีการเตรียมข้อมูลสำหรับ linear regression:

One-hot encoding ด้วย model.matrix()
Log transformation ด้วย log()
Split data ด้วย sample()

สร้าง linear regression model ด้วย lm() พร้อมประเมิน model ด้วย predict() และการคำนวณค่า MAE และ RMSE

😺 GitHub

ดู code ทั้งหมดในบทความนี้ได้ที่ GitHub

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

Correlation
t-tests
ANOVA
Reliability
Factor analysis

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists

2025-05-29

วิธีวิเคราะห์และแปลผล principal component analysis (PCA) ในภาษา R — ตัวอย่างการใช้ prcomp() เพื่อลดตัวแปรใน wine dataset
ในบทความนี้ เราจะมาทำความรู้จักกับ principal component analysis (PCA) ในภาษา R กัน
🧐 PCA คืออะไร?

PCA เป็น machine learning algorithm ประเภท unsupervised learning สำหรับลดตัวแปรในข้อมูล (dimensionality reduction) ในขณะที่ยังเก็บข้อมูลที่มีความสำคัญเอาไว้

ยกตัวอย่างเช่น เรามีข้อมูลลูกค้าที่มี 50 ตัวแปร เช่น อายุ เงินเดือน ประวัติการใช้จ่าย ประวัติการเป็นสมาชิก … เราสามารถใช้ PCA เพื่อลดจำนวนตัวแปรลง เช่น ลดเหลือ 5 ตัวแปร ที่ยังให้ข้อมูลเกี่ยวกับลูกค้าได้เทียบเท่ากับ 50 ตัวแปร

ในการทำงาน เราสามารถใช้ PCA เพื่อ:
- ลดเวลาและ resource ในการประมวลผล เพราะมีตัวแปรที่ต้องประมวลผลน้อยลง
- ช่วยทำความเข้าใจข้อมูล เพราะมีตัวแปรที่ต้องพิจารณาน้อยลง
- ใช้เตรียมข้อมูลสำหรับสร้าง machine learning models อื่น ๆ เช่น regression model, hierarchical regression
💻 PCA ในภาษา R

.

🧑‍💻 prcomp()

ในภาษา R เราสามารถทำ PCA ได้ด้วย prcomp() ซึ่งเป็น function ใน base R (ซึ่งหมายความว่า เราไม่ต้องติดตั้ง package เพิ่มเติม)

prcomp() ต้องการ 3 arguments ได้แก่:
```
prcomp(x, center, scale.)
```
- x คือ dataset ที่ต้องการใช้
- center คือ ตัวเลือกว่า เราจะลบ mean ออกจากข้อมูลดิบ เพื่อให้ dataset มีค่า mean เป็น 0 ไหม (recommend ให้เป็น TRUE)
- scale. คือ ตัวเลือกว่า เราจะหารข้อมูลดิบด้วย variance เพื่อให้ทุก column อยู่ในช่วงข้อมูลเดียวกันไหม (recommend ให้เป็น TRUE)
.

🍷 wine dataset

เราลองมาดูตัวอย่าง PCA กับ wine dataset จาก rattle package กัน

wine เป็นชุดข้อมูลที่มีลักษณะต่าง ๆ ของไวน์ เช่น ระดับแอลกอฮอล์ สี และความเข้มข้น

(Note: ดูข้อมูลเพิ่มเติมเกี่ยวกับ wine dataset ได้ที่ wine: The wine dataset from the UCI Machine Learning Repository.)

เราสามารถเรียกใช้ wine dataset ได้จาก rattle package ดังนี้:
```
# Install and load the package

## Install
install.packages("rattle")

## Load
library(rattle)

# -----------------------------------

# Load the dataset

## Load
data(wine)
```
เรียกดูตัวอย่างข้อมูลใน wine dataset:
```
# Preview
head(wine)
```
ผลลัพธ์:
```
  Type Alcohol Malic  Ash Alcalinity Magnesium Phenols Flavanoids Nonflavanoids Proanthocyanins Color  Hue Dilution Proline
1    1   14.23  1.71 2.43       15.6       127    2.80       3.06          0.28            2.29  5.64 1.04     3.92    1065
2    1   13.20  1.78 2.14       11.2       100    2.65       2.76          0.26            1.28  4.38 1.05     3.40    1050
3    1   13.16  2.36 2.67       18.6       101    2.80       3.24          0.30            2.81  5.68 1.03     3.17    1185
4    1   14.37  1.95 2.50       16.8       113    3.85       3.49          0.24            2.18  7.80 0.86     3.45    1480
5    1   13.24  2.59 2.87       21.0       118    2.80       2.69          0.39            1.82  4.32 1.04     2.93     735
6    1   14.20  1.76 2.45       15.2       112    3.27       3.39          0.34            1.97  6.75 1.05     2.85    1450
```
ดูโครงสร้างของ wine dataset:
```
# View the structure
str(wine)
```
ผลลัพธ์:
```
'data.frame':	178 obs. of  14 variables:
 $ Type           : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
 $ Alcohol        : num  14.2 13.2 13.2 14.4 13.2 ...
 $ Malic          : num  1.71 1.78 2.36 1.95 2.59 1.76 1.87 2.15 1.64 1.35 ...
 $ Ash            : num  2.43 2.14 2.67 2.5 2.87 2.45 2.45 2.61 2.17 2.27 ...
 $ Alcalinity     : num  15.6 11.2 18.6 16.8 21 15.2 14.6 17.6 14 16 ...
 $ Magnesium      : int  127 100 101 113 118 112 96 121 97 98 ...
 $ Phenols        : num  2.8 2.65 2.8 3.85 2.8 3.27 2.5 2.6 2.8 2.98 ...
 $ Flavanoids     : num  3.06 2.76 3.24 3.49 2.69 3.39 2.52 2.51 2.98 3.15 ...
 $ Nonflavanoids  : num  0.28 0.26 0.3 0.24 0.39 0.34 0.3 0.31 0.29 0.22 ...
 $ Proanthocyanins: num  2.29 1.28 2.81 2.18 1.82 1.97 1.98 1.25 1.98 1.85 ...
 $ Color          : num  5.64 4.38 5.68 7.8 4.32 6.75 5.25 5.05 5.2 7.22 ...
 $ Hue            : num  1.04 1.05 1.03 0.86 1.04 1.05 1.02 1.06 1.08 1.01 ...
 $ Dilution       : num  3.92 3.4 3.17 3.45 2.93 2.85 3.58 3.58 2.85 3.55 ...
 $ Proline        : int  1065 1050 1185 1480 735 1450 1290 1295 1045 1045 ...
```
Note: จะเห็นว่า ทุก column เป็นข้อมูลประเภท numeric ยกเว้น Type ที่ะเป็น factor

.

⚒️ PCA

หลังโหลด dataset แล้ว เราสามารถใช้ prcomp() เพื่อทำ PCA ได้ดังนี้:
```
# PCA
pca <- prcomp(wine[, -1],
              center = TRUE,
              scale. = TRUE)
```
อธิบาย code:
- เราใส่ wine[, -1] เพราะ PCA ใช้งานได้กับ column ที่เป็น numeric เท่านั้น เราเลย subset ข้อมูลเพื่อนำ Type ที่เป็น factor ออก
- เรากำหนด center = TRUE เพื่อให้ mean เท่ากับ 0
- กำหนดให้ scale. = TRUE เพื่อให้ข้อมูลอยู่ในช่วงข้อมูลเดียวกัน และป้องกันไม่ให้ข้อมูลที่อยู่ในมีช่วงข้อมูลกว้าง (เช่น ช่วง 1-100) มีผลต่อการวิเคราะห์มากกว่าข้อมูลที่มีช่วงแคบ (เช่น 1-10)
.

📈 ดูผลลัพธ์

เราสามารถดูผลลัพธ์ของ PCA ได้ 2 วิธี ได้แก่:
1. ดูค่าทางสถิติ
2. สร้างกราฟ
วิธีที่ 1. เราสามารถดูค่าทางสถิติได้ด้วย summary():
```
# Print the results
summary(pca)
```
ผลลัพธ์:
```
Importance of components:
                         PC1    PC2    PC3     PC4     PC5     PC6     PC7     PC8     PC9   PC10    PC11    PC12    PC13
Standard deviation     2.169 1.5802 1.2025 0.95863 0.92370 0.80103 0.74231 0.59034 0.53748 0.5009 0.47517 0.41082 0.32152
Proportion of Variance 0.362 0.1921 0.1112 0.07069 0.06563 0.04936 0.04239 0.02681 0.02222 0.0193 0.01737 0.01298 0.00795
Cumulative Proportion  0.362 0.5541 0.6653 0.73599 0.80162 0.85098 0.89337 0.92018 0.94240 0.9617 0.97907 0.99205 1.00000
```
ในผลลัพธ์ เราจะเห็นรายละเอียดดังนี้:
- จำนวนตัวแปรที่ PCA สร้างให้ (PC ซึ่งย่อมาจาก principal component) ซึ่งในตัวอย่างมีทั้งหมด 13 ตัวแปร
- Standard deviation (SD) ของแต่ละ PC
- สัดส่วนข้อมูลที่อธิบายได้ด้วย PC แต่ละตัว (proportion of variance)
- สัดส่วนข้อมูลสะสมที่อธิบายได้ เมื่อเพิ่ม PC แต่ละตัว (cumulative proportion) อย่างในตัวอย่าง จะเห็นว่า เมื่อเพิ่ม PC ตัวที่ 5 เราสามารถอธิบายข้อมูลได้ถึง 80% แล้ว แสดงว่า เราสามารถเก็บตัวแปรไว้ 5 จาก 13 ตัวแปรได้
(Note: เราควรใช้ PC เท่ากับจำนวนที่สามารถอธิบายข้อมูลได้ตั้งแต่ 80% ขึ้นไป)

วิธีที่ 2. เราสามารถดูผลลัพธ์ผ่านกราฟได้ เช่น scree plot ที่แสดงจำนวน PC และสัดส่วนข้อมูลสะสมที่อธิบายได้

ในการเริ่มสร้างกราฟ ให้เราคำนวณหาสัดส่วน variance และสัดส่วนสะสมที่อธิบายได้ก่อน:
```
# Extract variance explained
pca_var <- pca$sdev^2
pca_var_exp <- pca_var / sum(pca_var)

# Compute cumulative variance explained
cum_var_exp <- cumsum(pca_var_exp)
```
จากนั้น นำผลลัพธ์ไปสร้างกราฟ:
```
# Plot a scree plot for cumulative variance explained
plot(cum_var_exp, 
     type = "b", col = "blue", pch = 19, lwd = 2,
     main = "Cumulative Variance Explained",
     xlab = "Number of Principal Components",
     ylab = "Cumulative Variance Explained",
     ylim = c(0, 1))
```
ผลลัพธ์:

จากกราฟ เราจะได้ข้อสรุปเดียวกันกับค่าทางสถิติ นั่นคือ เมื่อเรามี PC 5 ตัว เราจะสามารถอธิบายข้อมูลได้ถึง 80% ของข้อมูลเก่า

😺 GitHub

ดู code ทั้งหมดในบทความนี้ได้ที่ GitHub

📃 References
✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

🙋 ใครที่กำลังเรียนจิตวิทยาหรือทำงานสายจิตวิทยา และเบื่อที่ต้องใช้ software ราคาแพงอย่าง SPSS และ Excel เพื่อทำข้อมูล

💪 ผมขอแนะนำ R Book for Psychologists หนังสือสอนใช้ภาษา R เพื่อการวิเคราะห์ข้อมูลทางจิตวิทยา ที่เขียนมาเพื่อนักจิตวิทยาที่ไม่เคยมีประสบการณ์เขียน code มาก่อน

ในหนังสือ เราจะปูพื้นฐานภาษา R และพาไปดูวิธีวิเคราะห์สถิติที่ใช้บ่อยกัน เช่น:
- Correlation
- t-tests
- ANOVA
- Reliability
- Factor analysis
🚀 เมื่ออ่านและทำตามตัวอย่างใน R Book for Psychologists ทุกคนจะไม่ต้องพึง SPSS และ Excel ในการทำงานอีกต่อไป และสามารถวิเคราะห์ข้อมูลด้วยตัวเองได้ด้วยความมั่นใจ

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists
Share this:
X
Facebook
Like Loading…
2025-05-15
วิธีทำ Hierarchical Clustering ในภาษา R — ตัวอย่างการใช้ hclust() จัดกลุ่มข้อมูลอาชญากรรมจาก USArrests dataset
ในบทความนี้ เราจะมาทำความรู้จักกับ hierarchical cluster analysis (HCA) หรือ hierarchical clustering ในภาษา R กัน
🧐 HCA คืออะไร?

HCA เป็น machine learning algorithm ประเภท unsupervised learning และใช้จัดกลุ่มข้อมูล (clustering) แบบลำดับชั้น (hierarchical)

HCA มี 2 ประเภท:
1. Agglomerative (bottom-up): จับ data point ที่อยู่ใกล้กันให้อยู่ cluster เดียวกัน แล้วจับ clusters ที่อยู่ใกล้กันให้อยูากลุ่มเดียวกัน ทำอย่างนี้ไปเรื่อย ๆ จนได้ 1 cluster ใหญ่
2. Divisive (top-down): เริ่มจาก 1 cluster และแยก cluster ย่อยออกมาเรื่อย ๆ
ตัวอย่างการใช้งาน HCA เช่น จัดกลุ่มผู้บริโภค (customer segmentation):

💻 HCA ในภาษา R: hclust()

ในภาษา R เราสามารถทำ HCA ได้ด้วย 2 functions:
1. hclust() สำหรับ agglomerative HCA
2. diana() จาก cluster package สำหรับ divisive HCA
ในบทความนี้ เราจะดูวิธีใช้ hclust() ซึ่งเป็นที่นิยมสำหรับ HCA กัน

🔢 Example Dataset: USArrests

Dataset ที่เราจะใช้เป็นตัวอย่างในบทความนี้ คือ USArrests

USArrests มีข้อมูลจำนวนประชากรจาก 50 รัฐในอเมริกา และข้อมูลการจับกุมใน 3 ประเภทการกระทำผิด ได้แก่:
1. ฆาตกรรม (murder)
2. ทำร้ายร่างกาย (assault)
3. ข่มขืน (rape)
ในการใช้งาน เราสามารถโหลด USArrests ได้ด้วย data():
```
# Load
data(USArrests)
```
จากนั้น ดูตัวอย่างข้อมูลด้วย head():
```
# Preview
head(USArrests)
```
ผลลัพธ์:
```
           Murder Assault UrbanPop Rape
Alabama      13.2     236       58 21.2
Alaska       10.0     263       48 44.5
Arizona       8.1     294       80 31.0
Arkansas      8.8     190       50 19.5
California    9.0     276       91 40.6
Colorado      7.9     204       78 38.7
```
📏 Normalise the Data

เนื่องจาก HCA จัดกลุ่มข้อมูลโดยใช้ระยะห่างระหว่างข้อมูล เราจะต้อง normalise ข้อมูล เพื่อให้ข้อมูลมีช่วงข้อมูลที่เท่า ๆ กัน และป้องกันไม่ให้ข้อมูลที่มีระยะห่างมาก (เช่น ข้อมูลที่อยู่ในช่วง 1 ถึง 1,000 เทียบกับ ข้อมูลที่มีช่วง 1 ถึง 10) มีผลต่อการจัดกลุ่ม

สำหรับ HCA เราจะ normalise ด้วย z-score standardisation ซึ่งมีสูตรคำนวณดังนี้:
```
Z = (X - M) / SD
```
- Z = ข้อมูลที่ normalised แล้ว
- X = ข้อมูลตั้งต้น
- M = ค่าเฉลี่ยของข้อมูลตั้งต้น
- SD = standard deviation (SD) ของข้อมูลตั้งต้น
ทั้งนี้ เราสามารถทำ z-score standardisation ได้ด้วย scale():
```
# Perform z-score standardisation
USArrests_scaled <- scale(USArrests)
```
จากนั้น เช็กผลลัพธ์ด้วย colMeans() และ apply() กับ sd():
```
# Check the results

## Mean
colMeans(USArrests_scaled)

## SD
apply(USArrests_scaled, 2, sd)
```
ผลลัพธ์:
```
> ## Mean
> round(colMeans(USArrests_scaled), 2)
  Murder  Assault UrbanPop     Rape 
       0        0        0        0 
> 
> ## SD
> apply(USArrests_scaled, 2, sd)
  Murder  Assault UrbanPop     Rape 
       1        1        1        1 
```
จะเห็นได้ว่า ทุก column มี mean เป็น 0 และ SD เป็น 1 แสดงว่า เรา normalise ข้อมูลได้สำเร็จ และพร้อมไปขั้นตอนถัดไป

🍿 HCA With hclust()

ตอนนี้ เราจะเริ่มใช้ HCA เพื่อจัดกลุ่มข้อมูลกัน

hclust() ต้องการ input 2 อย่าง:
```
hclust(d, method)
```
Input #1. d หมายถึง distance matrix หรือ matrix ที่เก็บค่าระยะห่างระหว่างข้อมูลแบบ pair-wise ไว้

เราสามารถคำนวณ distance matrix ได้ด้วย dist() function

Input #2. method หมายถึง linkage method หรือวิธีในการจับ clusters รวมกันเป็น cluster ที่ใหญ่ขึ้น

Linkage methods มี 5 ประเภท:
1. Ward’s linkage: จับกลุ่ม cluster โดยลด variance ภายใน cluster
2. Complete linkage: ใช้ระยะห่างสูงสุด
3. Single linkage: ใช้ระยะห่างต่ำที่สุด
4. Average linkage: ใช้ระยะห่างเฉลี่ยระหว่างข้อมูลจาก 2 clusters
5. Centroid linkage: ใช้ระยะห่างระหว่างจุดศูนย์กลางของ 2 clusters
Linkage ที่มักใช้กันได้ แก่ Ward’s, complete, และ average และในบทความนี้ เราจะใช้ Ward’s linkage ("ward.D2") กัน

เราเรียกใช้ hclust() เพื่อจัดกลุ่มข้อมูล:
```
# Create a distance matrix
dm <- dist(USArrests_scaled)

# HCA
hc <- hclust(dm,
             method = "ward.D2")
```
ตอนนี้ เราก็ได้ข้อมูลที่จัดกลุ่มแล้ว

📈 Print & Dendrogram

สุดท้าย เราสามารถดูผลลัพธ์ของ HCA ได้ด้วย print():
```
# Print HCA
print(hc)
```
ผลลัพธ์:
```
Call:
hclust(d = dm, method = "ward.D2")

Cluster method   : ward.D2 
Distance         : euclidean 
Number of objects: 50 
```
และการสร้าง dendrogram:
```
# Plot a dendrogram
plot(hc,
     hang = -1,
     cex = 0.8,
     main = "Dendrogram of USArrests Data")
```
ผลลัพธ์:

😺 GitHub

ดู code ทั้งหมดในบทความนี้ได้ที่ GitHub

📃 References
✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

🙋 ใครที่กำลังเรียนจิตวิทยาหรือทำงานสายจิตวิทยา และเบื่อที่ต้องใช้ software ราคาแพงอย่าง SPSS และ Excel เพื่อทำข้อมูล

💪 ผมขอแนะนำ R Book for Psychologists หนังสือสอนใช้ภาษา R เพื่อการวิเคราะห์ข้อมูลทางจิตวิทยา ที่เขียนมาเพื่อนักจิตวิทยาที่ไม่เคยมีประสบการณ์เขียน code มาก่อน

ในหนังสือ เราจะปูพื้นฐานภาษา R และพาไปดูวิธีวิเคราะห์สถิติที่ใช้บ่อยกัน เช่น:
- Correlation
- t-tests
- ANOVA
- Reliability
- Factor analysis
🚀 เมื่ออ่านและทำตามตัวอย่างใน R Book for Psychologists ทุกคนจะไม่ต้องพึง SPSS และ Excel ในการทำงานอีกต่อไป และสามารถวิเคราะห์ข้อมูลด้วยตัวเองได้ด้วยความมั่นใจ

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists
Share this:
X
Facebook
Like Loading…
2025-05-08

วิธีจัดกลุ่มข้อมูลด้วย k-means ผ่าน kmeans() function ในภาษา R — ตัวอย่างการจัดกลุ่มหินจาก rock dataset

ในบทความนี้ เราจะมาทำความรู้จักกับ k-means และวิธีใช้ kmeans() ในภาษา R กัน

👉 Introduction to k-Means

🤔 k-Means คืออะไร?

k-means เป็น machine learning algorithm ประเภท unsupervised learning และใช้จัดกลุ่ม data (clustering) ที่เราไม่รู้จำนวนกลุ่มล่วงหน้า

ตัวอย่างการใช้ k-means ในโลกจริง:

Customer segmentation: จัดกลุ่มลูกค้าที่เข้ามาซื้อสินค้า/บริการ
Anomaly detection: ตรวจจับความผิดปกติในคอมพิวเตอร์ (แบ่งกลุ่มกิจกรรมในคอมพิวเตอร์ เป็น “ปกติ” และ “ไม่ปกติ”)
Document clustering: จัดกลุ่มเอกสาร โดยอ้างอิงจากเนื้อหา

🪜 Steps การทำงานของ k-Means

k-means มีการทำงานอยู่ 5 ขั้นตอน ได้แก่:

กำหนดจำนวนกลุ่ม หรือ clusters (k)
สุ่มวาง centroid หรือจุดศูนย์กลางของ cluster ลงในข้อมูล
จัดกลุ่มข้อมูล โดยข้อมูลจะอยู่กลุ่มเดียวกับ centroid ที่ใกล้ที่สุด
คำนวณหา centroid ใหม่
ทำขั้นที่ 3 และ 4 ซ้ำ ตามจำนวนครั้งที่กำหนด หรือจนกว่าข้อมูลจะไม่เปลี่ยนกลุ่ม

ตัวอย่างเช่น เราต้องการจัดกลุ่มข้อมูล 100 ตัวอย่าง:

ขั้นที่ 1. เรากำหนด k เช่น ให้ k = 3

ขั้นที่ 2. สุ่มวาง centroid (ดอกจันสีแดง) ลงในข้อมูล:

ขั้นที่ 3. จัดข้อมูลให้อยู่กลุ่มเดียวกัน โดยอิงจาก centroid ที่อยู่ใกล้ที่สุด:

ขั้นที่ 4. คำนวณหา centroids ใหม่:

จะสังเกตเห็นว่า centroids ของเราเปลี่ยนไป

ขั้นที่ 5. ทำขั้นที่ 3 และ 4 ซ้ำ ๆ ไปเรื่อย ๆ เช่น ทำไป 10 ครั้ง:

เราก็จะได้การจัดกลุ่มข้อมูลมา

🏫 Learn More

เรียนรู้เพิ่มเติมเกี่ยวกับ k-means ได้จากคอร์ส The Nuts and Bolts of Machine Learning จาก Google Career Certificates (เริ่มต้นที่ 1:33:49)

💻 k-Means ในภาษา R: kmeans()

ในภาษา R เราสามารถใช้ k-means ได้ผ่าน kmeans() function ซึ่งต้องการ 3 arguments หลัก ดังนี้:

kmeans(x, centers, nstart)

x = dataset ที่ต้องการจัดกลุ่ม
centers = จำนวนกลุ่มข้อมูล หรือ k
nstart = จำนวนครั้งที่ k-means จะสุ่มวาง centroids ลงใน dataset เพื่อหาการจัดกลุ่มที่ดีที่สุด

(Note: ศึกษา arguments เพิ่มเติมของ kmeans() ได้ที่ kmeans: K-Means Clustering)

เราไปดูวิธีใช้งาน kmeans() กัน

🔢 Dataset

🪨 rock

ในบทความนี้ เราจะใช้ rock dataset ซึ่งเป็น built-in dataset ในภาษา R เป็นตัวอย่าง

rock มีข้อมูลหิน 48 ตัวอย่าง และมีข้อมูลลักษณะ 4 อย่าง:

ลำดับ	ลักษณะ	คำอธิบาย
1	`area`	พื้นที่ผิว
2	`peri`	เส้นผ่านศูนย์กลาง
3	`shape`	ขนาด (เส้นผ่านศูนย์กลาง หารด้วย พื้นที่ผิว)
4	`perm`	ระดับความสามารถที่ให้น้ำซึมผ่านได้

เป้าหมายของเรา คือ จัดกลุ่มหิน 48 ตัวอย่างตามลักษณะทั้งสี่

เราสามารถโหลด rock dataset ได้โดยใช้ data() function:

# Load
data(rock)

หลังจากโหลดแล้ว เราสามารถ preview ข้อมูลได้ด้วย head():

# Preview the dataset
head(rock)

ผลลัพธ์:

  area    peri     shape perm
1 4990 2791.90 0.0903296  6.3
2 7002 3892.60 0.1486220  6.3
3 7558 3930.66 0.1833120  6.3
4 7352 3869.32 0.1170630  6.3
5 7943 3948.54 0.1224170 17.1
6 7979 4010.15 0.1670450 17.1

📏 Data Normalisation

เนื่องจาก k-means ใช้ระยะห่างระหว่างข้อมูลในการจัดกลุ่ม ข้อมูลที่มีระยะห่างมาก (เช่น ระยะห่าง 500 เมตร กับ 1,000 เมตร) อาจมีผลต่อการจัดกลุ่มมากกว่าข้อมูลที่มีระยะห่างน้อย (เช่น ระยะห่าง 1 ซม. กับ 5 ซม.) และทำให้เกิด bias ในการจัดกลุ่มได้

เพื่อป้องกัน bias เราควร normalise ข้อมูล หรือการปรับให้ข้อมูลอยู่ใน scale เดียวกัน

สำหรับ k-means เรามักจะ normalise ข้อมูล ด้วย z-score standardisation ซึ่งมีสูตรการคำนวณดังนี้:

Z = (X - M(X)) / SD(X)

Z = ข้อมูลที่ scaled แล้ว
X = ข้อมูลดิบ
M(X) = mean ของ X
SD(X) = SD ของ X

ในภาษา R เราสามารถทำ z-score standardisation ได้ด้วย scale():

# Scale
rock_scaled <- scale(rock)

Note: เราใส่ rock ไปใน argument และเก็บผลลัพธ์ไว้ใน data frame ใหม่ ชื่อ rock_scaled

เราสามารถเช็กผลลัพธ์ได้ด้วย:

colMeans() เพื่อเช็ก mean
apply() และ sd() เพื่อเช็ก SD

# Check the results

## Check mean
round(colMeans(rock_scaled), 2)

# Check SD
apply(rock_scaled, 2, sd)

ผลลัพธ์:

> ### Check mean
> round(colMeans(rock_scaled), 2)
 area  peri shape  perm 
    0     0     0     0 
> 
> ### Check SD
> apply(rock_scaled, 2, sd)
 area  peri shape  perm 
    1     1     1     1

จากผลลัพธ์ เราจะเห็นได้ว่า ทุกลักษณะมี mean เท่ากับ 0 และ SD เท่ากับ 1 แสดงว่า เราทำ z-score standardisation ได้สำเร็จ และพร้อมไปขั้นตอนถัดไป

🔎 Finding the Optimal k

ในการจัดกลุ่มด้วย k-means เราต้องเริ่มด้วยการกำหนด k

แต่เราจะรู้ได้ยังไงว่า k ควรมีค่าเท่าไร?

เรามี 3 วิธีในการหาค่า k ที่ดีที่สุด (optimal k):

Elbow method
Silhouette analysis
Gap analysis

ในบทความนี้ เราจะมาใช้วิธีแรกกัน: elbow method

(Note: เรียนรู้เกี่ยวทั้งสามวิธีได้ที่ ML | Determine the optimal value of K in K-Means Clustering)

Elbow method หาค่า k ที่ดีที่สุด โดยสร้างกราฟระหว่างค่า k และ within-cluster sum of squares (WSS) หรือระยะห่างระหว่างข้อมูลในกลุ่ม ค่า k ที่ดีที่สุด คือ ค่า k ที่ WSS เริ่มไม่ลดลง

ในภาษา R เราสามารถเริ่มสร้างกราฟได้ โดยเริ่มจากใช้ for loop หา WSS สำหรับช่วงค่า k ที่เราต้องการ

ในตัวอย่าง rock dataset เราจะใช้ช่วงค่า k ระหว่าง 1 ถึง 15:

# Initialise a vector for within cluster sum of squares (wss)
wss <- numeric(15)

# For-loop through the wss
for (k in 1:15) {
  
  ## Try the k
  km <- kmeans(rock_scaled,
               centers = k,
               nstart = 20)
  
  ## Get WSS for the k
  wss[k] <- km$tot.withinss
}

จากนั้น ใช้ plot() สร้างกราฟความสัมพันธ์ระหว่างค่า k และ WSS:

# Plot the wss
plot(1:15,
     wss,
     type = "b",
     main = "The Number of Clusters vs WSS",
     xlab = "Number of Clusters",
     ylab = "WSS")

ผลลัพธ์:

จากกราฟ จะเห็นว่า WSS เริ่มชะลอตัว เมื่อค่า k อยู่ที่ 3 และ 4 ซึ่งเป็นจุดที่เป็นข้อศอก (elbow) ของกราฟ:

แสดงว่า optimal k มีค่า 3 หรือ 4

สำหรับบทความนี้ เราจะกำหนด optimal k = 4:

# Set optiomal k = 4
opt_k <- 4

⏹️ k-Means

🔥 Train the Model

หลังเตรียมข้อมูลและหา optimal k แล้ว เราก็พร้อมที่จะใช้ kmeans() ในการจัดกลุ่มข้อมูลแล้ว:

# Set see for reproducibility
set.seed(100)

# Train the model
km <- kmeans(rock_scaled,
             centers = opt_k,
             nstart = 20)

🤓 Get the Results

เราสามารถดูผลการจัดกลุ่มได้ 2 วิธี:

วิธีที่ 1. ดูค่าทางสถิติ:

# Print the model
print(km)

ผลลัพธ์:

K-means clustering with 4 clusters of sizes 10, 18, 6, 14

Cluster means:
        area       peri
1 -0.4406840 -0.9164442
2  0.3496197  0.8022501
3  1.5646450  1.3101981
4 -0.8052989 -0.9383748
       shape       perm
1  1.5369800  1.1775435
2 -0.7319607 -0.8017914
3  0.2358392 -0.7425075
4 -0.2578245  0.5079897

Clustering vector:
 [1] 2 2 2 2 2 2 2 2 2 2 2 2 3
[14] 2 2 2 3 3 3 2 2 3 3 2 1 4
[27] 4 4 4 1 1 1 4 1 4 1 4 1 4
[40] 4 1 4 1 1 4 4 4 4

Within cluster sum of squares by cluster:
[1] 18.917267  8.104718
[3]  1.471788 22.095636
 (between_SS / total_SS =  73.1 %)

Available components:

[1] "cluster"     
[2] "centers"     
[3] "totss"       
[4] "withinss"    
[5] "tot.withinss"
[6] "betweenss"   
[7] "size"        
[8] "iter"        
[9] "ifault"

จากผลลัพธ์ จะเห็นได้ว่า ข้อมูลถูกแบ่งเป็น 4 กลุ่ม และเราสามารถดูค่า mean และ WSS ของแต่ละกลุ่มได้

วิธีที่ 2. สร้างกราฟ:

# Create a plot
plot(rock_scaled[, c("shape", "perm")], 
     col = km$cluster,
     pch = 19,
     main = "K-Means Clustering (Rock Dataset)",
     xlab = "Shape",
     ylab = "Permeability")

# Add cluster centers
points(km$centers[, c("shape", "perm")], 
       col = 1:5,
       pch = 4,
       cex = 2,
       lwd = 2)

ผลลัพธ์:

ตอนนี้ เราก็ได้ข้อมูลที่จัดกลุ่มด้วย k-means เรียบร้อยแล้ว 👍

😺 GitHub

ดู code ทั้งหมดในบทความนี้ได้ที่ GitHub

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

Correlation
t-tests
ANOVA
Reliability
Factor analysis

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists

2025-05-01

สอนปลูกต้นไม้ในภาษา R (ภาค 1): วิธีสร้าง tree-based models ใน 3 ขั้นตอนด้วย rpart และ randomForest packages — ตัวอย่างการทำนายประเภทเกียร์รถใน mtcars dataset

ในบทความนี้ เราจะมาทำความรู้จักและสร้าง tree-based models ในภาษา R กัน:

Classification tree
Random forest

ถ้าพร้อมแล้ว ไปเริ่มกันเลย

🌳 Tree-Based Algorithms คืออะไร?

Tree-based algorithm เป็น machine learning algorithm ประเภท supervised learning ที่ทำนายข้อมูลโดยใช้ decision tree

Decision tree เป็นเครื่องมือช่วยตัดสินใจที่ดูคล้ายกับต้นไม้กลับหัว เพราะประกอบด้วยจุดตัดสินใจ (node) ที่แตกยอด (ทางเลือก) ออกไปเรื่อย ๆ เพื่อใช้ทำนายผลลัพธ์ที่ต้องการ

ยกตัวอย่างเช่น เราอยากรู้ว่า เราควรจะสมัครงานกับบริษัทแห่งหนึ่งไหม เราอาจจะสร้าง decision tree จาก 3 ปัจจัย คือ:

เงินเดือน
ได้ใช้ทักษะที่มี
การทำงานแบบ hybrid หรือ remote

ได้แบบนี้:

จาก decision tree ถ้าเราเห็นงานที่ได้เงินเดือนไม่ตรงใจ เราจะไม่สมัครงานนั้น (เส้นการตัดสินใจซ้ายสุด)

ในทางกลับกัน ถ้างานมีเงินเดือนที่น่าสนใจ และได้ใช้ทักษะที่มีอยู่ เราจะสมัครงานนั้น (เส้นการตัดสินใจขวาสุด) เป็นต้น

💻 Tree-Based Models ในภาษา R

ในภาษา R เราสามารถสร้าง tree-based model ได้ด้วย rpart() จาก rpart package:

# Install
install.packages("rpart")
# Load
library(rpart)

rpart() ต้องการ 4 arguments ดังนี้:

rpart(formula, data, method, control)

formula = สูตรในการวิเคราะห์ (ตัวแปรตาม ~ ตัวแปรต้น)
data = dataset ที่ใช้สร้าง model
method = ประเภท algorithm ("anova" สำหรับ regression และ "class" สำหรับ classification)
control (optional) = เงื่อนไขควบคุม “การเติบโต” ของ decision tree เช่น ระดับชั้นที่มีได้ เป็นต้น

(Note: ศึกษาการใช้งาน rpart() เพิ่มเติมได้ที่ rpart: Recursive Partitioning and Regression Trees)

เราไปดูตัวอย่างการใช้งาน rpart() กัน

🚗 Dataset: mtcars

ในบทความนี้ เราจะลองสร้าง tree-based model ประเภท classification หรือ classification tree เพื่อทำนายประเภทเกียร์รถใน mtcars dataset

mtcars เป็นชุดข้อมูลรถจาก ปี ค.ศ. 1974 ซึ่งประกอบไปด้วยข้อมูล เช่น รุ่นรถ น้ำหนัก ระดับการกินน้ำมัน แรงม้า เป็นต้น

เราสามารถโหลด mtcars มาใช้งานได้ด้วย data() และ preview ด้วย head():

# Load
data(mtcars)
# Preview
head(mtcars)

ตัวอย่างข้อมูล:

                   mpg cyl disp  hp drat    wt  qsec vs        am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0    manual    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0    manual    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1    manual    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1 automatic    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0 automatic    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1 automatic    3    1

🔧 Prepare the Data

ก่อนนำ mtcars ไปใช้สร้าง classification tree เราจะต้องทำ 2 อย่างก่อน:

อย่างที่ #1. ปรับ column am ให้เป็น factor เพราะสิ่งที่เราต้องการทำนายเป็น categorical data:

# Convert `am` to factor
mtcars$am <- factor(mtcars$am,
                    levels = c(0, 1),
                    labels = c("automatic", "manual"))
# Check the result
class(mtcars$am)

ผลลัพธ์:

[1] "factor"

อย่างที่ #2. Split ข้อมูลเป็น 2 ชุด:

Training set สำหรับสร้าง model
Test set สำหรับประเมิน model

# Set seed for reproducibility
set.seed(500)
# Get training index
train_index <- sample(nrow(mtcars),
                      nrow(mtcars) * 0.7)
# Split the data
train_set <- mtcars[train_index, ]
test_set <- mtcars[-train_index, ]

🪴 Train the Model

ตอนนี้ เราพร้อมที่จะสร้าง classification tree ด้วย rpart() แล้ว

สำหรับ classification tree ในบทความนี้ เราจะลองตั้งเงื่อนไขในการปลูกต้นไม้ (control) ดังนี้:

Argument	Explanation
`cp = 0`	คะแนนประสิทธิภาพขั้นต่ำ ก่อนจะแตกกิ่งใหม่ได้ = 0
`minsplit = 1`	จำนวนกิ่งย่อยขั้นต่ำที่ต้องมี ก่อนจะแตกกิ่งใหม่ได้ = 1
`maxdepth = 5`	จำนวนชั้นที่ decision tree มีได้สูงสุด = 5

# Classification tree
ct <- rpart(am ~ .,
            data = train_set,
            method = "class",
            control = rpart.control(cp = 0, minsplit = 1, maxdepth = 5))

เราสามารถดู classification tree ของเราได้ด้วย rpart.plot():

# Plot classification tree
rpart.plot(ct,
           type = 3,
           extra = 101,
           under = TRUE,
           digits = 3,
           tweak = 1.2)

ผลลัพธ์:

Note:

ก่อนใช้งาน rpart.plot() เราต้องจะติดตั้งและเรียกใช้งานด้วย install.packages() และ library() ตามลำดับ
ศึกษาการใช้งาน rpart.plot() เพิ่มเติมได้ที่ rpart.plot: Plot an rpart model. A simplified interface to the prp function.

📏 Evaluate the Model

เมื่อได้ classification tree มาแล้ว เราลองมาประเมินความสามารถของ model ด้วย accuracy หรือสัดส่วนการทำนายที่ถูกต้องกัน

เราเริ่มจากใช้ predict() เพื่อใช้ model ทำนายประเภทเกียร์:

# Predict the outcome
test_set$pred_ct <- predict(ct,
                            newdata = test_set,
                            type = "class")

จากนั้น สร้าง confusion matrix หรือ matrix เปรียบเทียบคำทำนายกับข้อมูลจริง:

# Create a confusion matrix
cm_ct <- table(Predicted = test_set$pred_ct,
               Actual = test_set$am)
# Print confusion matrix
print(cm_ct)

ผลลัพธ์:

           Actual
Predicted   automatic manual
  automatic         5      1
  manual            0      4

สุดท้าย เราหา accuracy ด้วยการนำจำนวนคำทำนายที่ถูกต้องมาหารด้วยจำนวนคำทำนายทั้งหมด:

# Get accuracy
acc_ct <- sum(diag(cm_ct)) / sum(cm_ct)
# Print accuracy
cat("Accuracy (classification tree):", acc_ct)

ผลลัพธ์:

Accuracy (classification tree): 0.9

จะเห็นได้ว่า model ของเรามีความแม่นยำถึง 90%

🍄 Random Forest

Random forest เป็น tree-based algorithm ที่ช่วยเพิ่มความแม่นยำในการทำนาย โดยสุ่มสร้าง decision trees ต้นเล็กขึ้นมาเป็นกลุ่ม (forest) แทนการปลูก decision tree ต้นเดียว

Decision tree แต่ละต้นใน random forest มีความสามารถในการทำนายแตกต่างกัน ซึ่งบางต้นอาจมีความสามารถที่น้อยมาก

แต่จุดแข็งของ random forest อยู่ที่จำนวน โดย random forest ทำนายผลลัพธ์โดยดูจากผลลัพธ์ในภาพรวม ดังนี้:

Task	Predict by
Regression	ค่าเฉลี่ยของผลลัพธ์การทำนายของทุกต้น
Classification	เสียงส่วนมาก (majority vote)

ดังนั้น แม้ว่า decision tree บางต้นอาจทำนายผิดพลาด แต่โดยรวมแล้ว random forest มีโอกาสที่จะทำนายได้ดีกว่า decision tree ต้นเดียว

ในภาษา R เราสามารถสร้าง random forest ได้ด้วย randomForest() จาก randomForest package ซึ่งต้องการ 3 arguments:

randomFrest(formula, data, ntree)

formula = สูตรในการวิเคราะห์ (ตัวแปรตาม ~ ตัวแปรต้น)
data = dataset ที่ใช้สร้าง model
ntree = จำนวน decision trees ที่ต้องการสร้าง

Note:

เราไม่ต้องกำหนดว่า จะทำ classification หรือ regression model เพราะ randomForest() จะเลือก model ให้อัตโนมัติตามข้อมูลที่เราใส่เข้าไป
ศึกษาการใช้งาน randomForest() เพิ่มเติมได้ที่ randomForest: Classification and Regression with Random Forest

ก่อนใช้ randomForest() เราต้องเตรียมข้อมูลแบบเดียวกันกับ rpart() ได้แก่:

เปลี่ยน am ให้เป็น factor
Split the data

สมมุติว่า เราเตรียมข้อมูลแล้ว เราสามารถเรียกใช้ randomForest() ได้เลย โดยเราจะลองสร้าง random forest ที่ประกอบด้วย decision trees 100 ต้น:

# Random forest
rf <- randomForest(am ~ .,
                   data = train_set,
                   ntree = 100)

แล้วลองประเมินความสามารถของ model ด้วย accuracy

เริ่มจากทำนายประเภทเกียร์:

# Predict the outcome
test_set$pred_rf <- predict(rf,
                            newdata = test_set,
                            type = "class")

สร้าง confusion matrix:

# Create a confusion matrix
cm_rf <- table(Predicted = test_set$pred_rf,
               Actual = test_set$am)
# Print confusion matrix
print(cm_rf)

ผลลัพธ์:

           Actual
Predicted   automatic manual
  automatic         5      0
  manual            0      5

และสุดท้าย คำนวณ accuracy:

# Get accuracy
acc_rf <- sum(diag(cm_rf)) / sum(cm_rf)
# Print accuracy
cat("Accuracy (random forest):", acc_rf)

ผลลัพธ์:

Accuracy (random forest): 1

จะเห็นว่า random forest (100%) มีความแม่นยำในการทำนายมากกว่า classification tree ต้นเดียว (90%)

🐱 GitHub

ดู code ทั้งหมดในบทความนี้ได้ที่ GitHub:

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

Correlation
t-tests
ANOVA
Reliability
Factor analysis

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists

2025-04-24

วิธีสร้างและประเมิน logistic regression model ใน 5 ขั้นตอน ด้วย glm(), roc(), และ auc() ในภาษา R — ตัวอย่างการทำนายประเภทเกียร์รถใน mtcars dataset

ในบทความนี้ เราจะมาทำความรู้จักกับ logistic regression model รวมทั้งการสร้างและประเมิน logistic regression ด้วย glm(), roc(), และ auc() ในภาษา R กัน

📈 Logistic Regression คืออะไร?

Logistic regression เป็น machine learning algorithm ประเภท classification สำหรับทำนายผลลัพธ์ที่มี 2 ค่า (binary) เช่น:

Spam filter: เป็น spam vs เป็นอีเมลปกติ
Churn prediction: ลูกค้าจะเลิกใช้บริการ vs ใช้บริการต่อ
Fraud detection: เป็นธุรกรรมทุจริต vs ไม่ทุจริต
Disease diagnosis: เป็นโรค vs ไม่เป็นโรค

โดยการทำงานของ logistic regression เป็นการต่อยอดจาก linear regression

ใน linear regression เราหาเส้นตรงที่ช่วยทำนายตัวแปรตาม (y) ด้วยตัวแปรต้น (x) ได้ดีที่สุด ตามสมการเส้นตรง:

y = mx + c

y = ตัวแปรตาม
m = ความชันของเส้นตรง
x ตัวแปรต้น
c = จุดตัดระหว่าง x และ y

Linear regression ใช้งานได้ดีกับข้อมูลที่เป็น continuous data และความสัมพันธ์ระหว่าง x และ y เป็นเส้นตรง:

แต่เส้นตรงไม่ใช่เส้นที่ดีที่สุด ในการทำนาย y ที่มีเพียง 2 ระดับ:

เพื่อแก้ปัญหานี้ เราต้องใช้ logistic regression ซึ่งใช้ sigmoid function ที่ทำนายความเป็นได้ (probability) ที่ y จะตกอยู่ในระดับใดระดับหนึ่ง แทนการทำนายค่า y โดยตรง:

y = sigmoid_func(mx + c)

y = ตัวแปรตาม
sigmoid_func = sigmoid function
m = ความชันของเส้นตรง
x ตัวแปรต้น
c = จุดตัดระหว่าง x และ y

เมื่อใช้ sigmoid function เราจะได้เส้นในการทำนายข้อมูลเรียกว่า sigmoid curve ที่ดูคล้ายตัว S ดังภาพ:

🔨 วิธีสร้าง Logistic Regression ในภาษา R

ในภาษา E เราสร้าง logistic regression model ได้ง่าย ๆ ด้วย glm() function โดยทำตาม 5 ขั้นตอนของ machine learning workflow:

Load the dataset
Prepare the data
Split the data
Train the model
Evaluate the model

1️⃣ Step 1. Load the Dataset

ในขั้นแรก ให้เราโหลด dataset ที่จะใช้งานเข้ามาใน R

โดยในบทความนี้ เราจะทำงานกับ mtcars dataset ซึ่งมข้อมูลรถจาก ปี ค.ศ. 1974 เช่น รุ่นรถ ระดับการกินน้ำมัน แรงม้า ประเภทเกียร์ และน้ำหนัก

เป้าหมายของเรา คือ สร้าง logistic regression model เพื่อทำนายประเภทเกียร์ ว่า เป็นรถเกียร์ auto (0) หรือ manual (1)

เนื่องจาก mtcars เป็น built-in dataset ในภาษา R เราสามารถโหลดข้อมูลได้ด้วย data() function:

# Load the dataset
data(mtcars)

หลังการโหลด เราสามารถ preview ข้อมูลได้ด้วย head():

# Preview the dataset
head(mtcars)

ผลลัพธ์:

                   mpg cyl disp  hp drat    wt  qsec vs        am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0    manual    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0    manual    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1    manual    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1 automatic    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0 automatic    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1 automatic    3    1

2️⃣ Step 2. Prepare the Data

สำหรับการสร้าง logistic regression model ตัวแปรตามของเราจะต้องมี data type เป็น factor

เราสามารถเช็ก data type ของ dataset ได้ด้วย str():

# Check the data typs of the columns
str(mtcars)

ผลลัพธ์:

'data.frame':	32 obs. of  11 variables:
 $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
 $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
 $ disp: num  160 160 108 258 360 ...
 $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
 $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
 $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
 $ qsec: num  16.5 17 18.6 19.4 17 ...
 $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
 $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
 $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
 $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

จากผลลัพธ์ เราจะเห็นได้ว่า am ซึ่งบอกประเภทเกียร์รถ มี data type เป็น num (numeric) แสดงว่า เราจะต้องทำ operation เพื่อเปลี่ยน data type ก่อนไปขั้นต่อไป

เราสามารถเปลี่ยน data type ได้ด้วย factor():

# Convert column `am` to factor
mtcars$am <- factor(mtcars$am,
                    levels = c(0, 1),
                    labels = c("automatic", "manual"))

จากนั้น เช็กผลลัพธ์ที่ได้ด้วย str() อีกครั้ง:

# Check the result
str(mtcars)

ผลลัพธ์:

'data.frame':	32 obs. of  11 variables:
 $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
 $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
 $ disp: num  160 160 108 258 360 ...
 $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
 $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
 $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
 $ qsec: num  16.5 17 18.6 19.4 17 ...
 $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
 $ am  : Factor w/ 2 levels "automatic","manual": 2 2 2 1 1 1 1 1 1 1 ...
 $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
 $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

จะเห็นว่า ตอนนี้ am เป็น factor แล้ว และเราสามารถไปขั้นตอนถัดไปได้

3️⃣ Step 3. Split the Data

ในขั้นที่สาม เราจะต้องแบ่ง dataset เป็น 2 ชุด ได้แก่:

Training set สำหรับสร้าง model
Test set สำหรับประเมินความสามารถของ model

เราแบ่ง dataset ได้โดยเริ่มจากสร้าง row index ที่จะอยู่ใน training set:

# Set seed for reproducibility
set.seed(300)

# Create a training index
train_index <- sample(1:nrow(mtcars),
                      nrow(mtcars) * 0.7)

Note: เรากำหนด row index ใน training set เป็น 70% แสดงว่า test set จะมีข้อมูล 30% จาก dataset

จากนั้น ใช้ row index แบ่ง dataset เป็น 2 ชุด:

# Create training and test sets
train_set <- mtcars[train_index, ]
test_set <- mtcars[-train_index, ]

ตอนนี้ เราพร้อมที่จะสร้าง model กันแล้ว

4️⃣ Step 4. Train the Model

ในขั้นที่สี่ เราจะสร้าง logistic regression model โดยใช้ glm() function กัน

glm() ต้องการ input 3 อย่าง:

glm(formula, data, family)

No.	Input	คำอธิบาย	ในบทความนี้
1	formula	ตัวแปรตามที่ต้องการทำนาย และตัวแปรต้นที่จะใช้ทำนาย	`am ~ .`
2	data	ชุดข้อมูลที่จะใช้สร้าง model	`train_set`
3	family	Regression algorithm ที่เราจะเรียกใช้	`"binomial"`

# Build a logistic regression model
log_reg <- glm(am ~ .,
               data = train_set,
               family = "binomial")

ตอนนี้ เราก็จะได้ logistic regression model ในการทำนายประเภทเกียร์มาเรียบร้อย

5️⃣ Step 5. Evaluate the Model

ในขั้นสุดท้าย เราจะประเมินความสามารถของ model ในการทำนายประเภทเกียร์กัน

เราจะใช้ 3 metrics ในการประเมิน:

Accuracy: สัดส่วนของผลลัพธ์ที่ทำนายถูก
ROC curve: กราฟแสดงความสัมพันธ์ระหว่าง true positive rate (TPR) และ false positive rate (FPR) โดยยิ่งเส้นกราฟออกห่างเหนือเส้นตรงมากเท่าไรก็ยิ่งดี
AUC: พื้นที่ใต้กราฟของ ROC โดยมีค่าระหว่าง 0 กับ 1 โดยค่ายิ่งเข้าใกล้ 1 เท่าไรก็ยิ่งดี

⚖️ Threshold

ในการคำนวณทั้ง 3 metrics เราจะต้องกำหนด threshold เพื่อทำนายปรพเภทเกียร์ก่อน เพราะ glm() ให้ model ที่ส่งค่าความเป็นไปได้ (แทนประเภทเกียร์) กลับมา เช่น:

0.45
0.32
0.93

Threshold จะเป็นค่าที่ระบุประเภทเกียร์ให้เรา เช่น เรากำหนด threshold ที่ 0.30:

ความเป็นไปได้	เทียบกับ 0.30	ประเภทเกียร์ที่ทำนาย
0.45	> 0.30	Manual
0.32	> 0.30	Manual
0.93	> 0.30	Manual

หรือ threshold ที่ 0.80:

ความเป็นไปได้	เทียบกับ 0.80	ประเภทเกียร์ที่ทำนาย
0.45	< 0.80	Auto
0.32	< 0.80	Auto
0.93	> 0.80	Manual

จะเห็นว่า threshold ที่เรากำหนด มีผลต่อผลลัพธ์ในการทำนาย

Note: ทั้งนี้ threshold ที่ดีที่สุดขึ้นอยู่กับแต่ละ model

ในบทความนี้ เราจะลองกำหนด threshold ไว้ที่ 0.50

เมื่อกำหนด threshold แล้ว เราจะใช้ predict() เพื่อสร้าง column ที่เก็บค่าความเป็นไปได้ ใน test_set:

# Get predictive probability
test_set$pred_prob <- predict(log_reg,
                              newdata = test_set,
                              type = "response")

จากนั้น ใช้ ifelse() และ threshold เพื่อทำนายประเภทเกียร์:

# Predict the outcome
test_set$pred <- ifelse(test_set$pred_prob > 0.5,
                        1,
                        0)

# Set column `pred ` as factor
test_set$pred <- factor(test_set$pred,
                        levels = c(0, 1),
                        labels = c("automatic", "manual"))

จากนั้น เราดูผลลัพธ์ได้ด้วย head():

# Check the results
head(test_set[c("pred_prob", "pred", "am")])

ผลลัพธ์:

                      pred_prob      pred        am
Datsun 710         1.000000e+00    manual    manual
Hornet 4 Drive     2.220446e-16 automatic automatic
Merc 240D          9.999999e-01    manual automatic
Merc 280C          2.220446e-16 automatic automatic
Cadillac Fleetwood 2.220446e-16 automatic automatic
Chrysler Imperial  3.941802e-11 automatic automatic

ตอนนี้ เราได้ผลลัพธ์ในการทำนายมา และพร้อมที่จะคำนวณ metrics ในการประเมินแล้ว

🎯 Metric #1. Accuracy

เรามาคำนวณแต่ละ metric กัน โดยเริ่มจาก accuracy

สำหรับ accuracy เราจะสร้าง confusion matrix ที่เปรียบการทำนายและประเภทเกียร์จริง ๆ ในข้อมูล:

# Create a confusion matrix
cm <- table(Predicted = test_set$pred,
            Actual = test_set$am)

# Print cm
print(cm)

ผลลัพธ์:

           Actual
Predicted   automatic manual
  automatic         6      1
  manual            1      2

จากนั้น เรานำสัดส่วนผลลัพธ์ที่ทำนายถูกมาหารด้วยผลลัพธ์ทั้งหมด:

# Compute accuracy
accuracy <- sum(diag(cm)) / sum(cm)

# Print accuracy
cat("Accuracy:", round(accuracy, 2))

ผลลัพธ์:

Accuracy: 0.8

จะเห็นได้ว่า model มีความแม่นยำในการทำนายอยู่ที่ 80%

📈 Metric #2. ROC Curve

สำหรับการสร้าง ROC curve เราจะเรียกใช้ roc() function จาก pROC package

Note: ในการใช้งาน pROC package ให้รันคำสั่ง install.packages("pROC") เพื่อติดตั้ง และ library(pROC) เพื่อเรียกใช้งาน

roc() ต้องการ input 2 อย่าง:

roc(actual, predicted)

actual = ประเภทเกียร์ในข้อมูลจริง
predicted = ความเป็นไปได้ในการทำ

# Calculate ROC
ROC <- roc(test_set$am,
           pred_prob)

จากนั้น เราสร้างกราฟ ROC curve ด้วย plot():

# Plot ROC
plot(ROC,
     main = "ROC Curve",
     col = "blue",
     lwd = 2)

ผลลัพธ์:

จากกราฟ ดูเหมือนว่า model ของเราอาจจะไม่ได้ทำได้ดีเท่ากับ accuracy เพราะเส้นกราฟ (สีน้ำเงิน) ไม่ได้ออกห่างจากเส้นตรงมากนัก และมีบางส่วนที่อยู่ใต้เส้นตรง

📉 Metric #3. AUC

สุดท้าย เรามาคำนวณพื้นที่ใต้กราฟ ด้วย auc() จาก pROC() package เพื่อดูว่า การตีความ ROC curve เป็นอย่างอย่างที่เราคิดไหม:

# Calculate AUC
AUC <- auc(ROC)

# Print AUC
print(AUC)

ผลลัพธ์:

Area under the curve: 0.6429

จะเห็นได้ว่า model ของเรายังมีประสิทธิภาพไม่ดีมาก เพราะมีค่าน้อยกว่า 1 มาก ซึ่งเป็นไปตามการอ่าน ROC curve ก่อนหน้านี้

Note: จากผลลัพธ์ เราอาจจะอยากไปปรับปรุง model ของเรา ซึ่งเราสมารถทำได้ 2 วิธี:

เพิ่มข้อมูลใน training set
ปรับ threshold

😺 GitHub

ดู code ทั้งหมดในบทความนี้ได้ที่ GitHub

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

Correlation
t-tests
ANOVA
Reliability
Factor analysis

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists

2025-04-17

วิธีสร้าง Naïve Bayes classifier ใน 5 ขั้นตอนด้วย naivebayes package ในภาษา R – ตัวอย่างการทำนายประเภทดอกไม้ใน iris dataset
ในบทความนี้ เราจะมาทำความรู้กับ Naïve Bayes และดูวิธีสร้าง Naïve Bayes classifier ในภาษา R ด้วย naivebayes package กัน
👼 Naïve Bayes คืออะไร?

Naïve Bayes (NB) เป็น machine learning algorithm ประเภท supervised learning สำหรับทำนายกลุ่มของข้อมูล (classification)

NB ทำนายข้อมูล โดยวิเคราะห์ความน่าจะเป็นที่ 2 เหตุการณ์ (ตัวแปรต้นและตัวแปรตาม) จะเกิดขึ้นคู่กัน

ยกตัวอย่างเช่น การอนุมัติสินเชื่อ (อนุมัติ/ไม่อนุมัติ) ที่เกิดขึ้นอยู่กับ:
- รายได้: รายได้มาก -> อนุมัติ
- อายุ: อายุน้อย -> อนุมัติ
- ประวัติการขอสินเชื่อ: เคยขอแล้วผ่าน -> อนุมัติ
NB คำนวณความน่าจะเป็นของแต่ละปัจจัยต่อตัวแปรตาม (การอนุมัติ) โดยอนุมานว่า แต่ละปัจจัยเกิดแยกกัน (independent) แม้ว่าในความเป็นจริงหลายปัจจัยจะเกิดขึ้นร่วมกัน (dependent; เช่น อายุมาก -> รายได้มาก) ก็ตาม

เนื่องจาก NB มีการมองโลกอย่างง่าย จึงได้ชื่อว่า “naïve” หรือไร้เดียงนั่นเอง

แม้ว่า NB จะมีการมองโลกที่ไร้เดียงสา แต่ก็เป็น algorithm ที่:
1. เข้าใจได้ง่าย (เพราะไร้เดียงสา)
2. รวดเร็ว (เพราะใช้การวิเคราะห์อย่างง่าย)
3. น่าเชื่อถือ (มักจะทำนายได้แม่นยำ แม้จะมีมุมมองที่ไม่ตรงกับความเป็นจริง)
🧑‍💻 วิธีสร้าง NB Classifer ด้วย naivebayes Package

วิธีการสร้าง NB classifier ด้วย naivebayes package มี 5 ขั้นตอน ได้แก่:
1. Install and load the package
2. Load and preview the dataset
3. Split the dataset
4. Create a NB model
5. Evaluate the model
.

1️⃣ Step 1. Install & Load the Package

ในขั้นแรก เราต้องติดตั้งและเรียกใช้งาน naivebayes package:
```
# Install and load the package
## Install
install.packages("naivebayes")
## Load
library(naivebayes)
```
Note: เราทำการติดตั้งแค่ครั้งแรกครั้งเดียว แต่การเรียกใช้ต้องทำทุกครั้งทุกเริ่ม session ใหม่

.

2️⃣ Step 2. Load & Preview the Dataset

ในขั้นที่สอง ให้เราโหลด dataset ที่จะใช้งาน

โดยในบทความนี้ เราจะใช้ built-in dataset ในภาษา R ชื่อ iris dataset

iris dataset มีกลีบดอกของข้อมูลดอกไม้ 3 ชนิด ได้แก่:
- Iris setosa
- Iris virginica
- Iris versicolor
เป้าหมายของเรา คือ ทำนายประเภทของดอกไม้จากความกว้างและความยาวกลีบดอก

เนื่องจากเป็น built-in dataset เราสามารถโหลด iris dataset ด้วยคำสั่ง data():
```
# Load
data(iris)
```
จากนั้น เราสามารถใช้คำสั่ง head() เพื่อ preview ข้อมูล:
```
# Preview
head(iris)
```
ผลลัพธ์:
```
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa
```
Note: จะเห็นว่า iris dataset มีข้อมูลทั้งหมด 5 columns คือ ความยาวและความกว้างของกลีบดอกชั้นนอก (sepal) และความยาวและความกว้างของกลีบดอกชั้นใน (petal) รวมทั้งประเภทดอกไม้ (species)

เนื่องจาก NB classifier ต้องการตัวแปรตามที่เป็น factor เราต้องเช็กได้ว่า Species เป็น data type ไหนด้วยคำสั่ง str():
```
# Preview
str(iris)
```
ผลลัพธ์:
```
'data.frame':	150 obs. of  5 variables:
 $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
 $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
```
จะเห็นได้ว่า Species เป็น factor อยู่แล้ว และพร้อมที่ใช้สร้าง NB classifier

.

3️⃣ Step 3. Split the Dataset

ในขั้นที่สาม เราจะแบ่งข้อมูลเป็น 2 ชุด:
1. Training set เพื่อสร้าง NB classifier
2. Test set เพื่อทดสอบประสิทธิภาพของ NB classifier
ขั้นแรกในการแบ่งข้อมูล ให้เราสุ่มเลือก row index ที่จะอยู่ใน training set:
```
# Set seed for reproducibility
set.seed(2025)
# Create a training index
training_index <- sample(1:nrow(iris),
                         0.7 * nrow(iris))
```
จากนั้น ใช้ row index ที่ได้แบ่งข้อมูลเป็น 2 ชุด:
```
# Split the dataset
train_set <- iris[training_index, ]
test_set <- iris[-training_index, ]
```
.

4️⃣ Step 4. Create a NB Model

ในขั้นที่สี่ เราจะสร้าง NB classifier โดยใช้ naive_bayes() จาก naivebayes package

naive_bayes() ต้องการ input 2 อย่าง:
1. Formula: ตัวแปรต้นและตัวแปรตามที่ต้องการใช้สร้าง model โดยในกรณีนี้ คือ ทำนาย Species ด้วยข้อมูลที่เหลือทั้งหมด (.)
2. Data: dataset ที่จะใช้สร้าง model ซึ่งในกรณีนี้ คือ train_set
```
# Create a NB model
nb <- naive_bayes(Species ~ .,
                  data = train_set)
```
เมื่อเราได้ NB classifier มาแล้ว เรามาลองใช้ model ทำนายประเภทดอกไม้กัน:
```
# Predict the outcomes
pred <- predict(nb,
                newdata = test_set[, 1:4],
                type = "class")
print(pred)
```
ผลลัพธ์:
```
 [1] setosa     setosa     setosa     setosa     setosa    
 [6] setosa     setosa     setosa     setosa     setosa    
[11] setosa     setosa     setosa     setosa     setosa    
[16] virginica  versicolor versicolor versicolor versicolor
[21] virginica  versicolor versicolor versicolor versicolor
[26] versicolor virginica  versicolor versicolor versicolor
[31] versicolor versicolor versicolor versicolor virginica 
[36] virginica  virginica  virginica  virginica  virginica 
[41] virginica  virginica  virginica  virginica  virginica 
Levels: setosa versicolor virginica
```
.

5️⃣ Step 5. Evaluate the Model

ในขั้นสุดท้าย เราจะประเมินประสิทธิภาพของ model ด้วยการคำนวณ accuracy หรือ สัดส่วนของคำตอบที่ model ตอบถูก:

Accuracy = Correct predictions / Total predictions

เราสามารถหาค่า accuracy ได้โดยเริ่มสร้าง confusion matrix ก่อน:
```
# Create a confusion matrix
cm <- table(Predicted = pred, 
            Actual = test_set$Species)
print(cm)
```
ผลลัพธ์:
```
            Actual
Predicted    setosa versicolor virginica
  setosa         15          0         0
  versicolor      0         16         0
  virginica       0          3        11
```
จากนั้น นำจำนวนคำตอบที่ถูกต้องและจำนวนคำตอบทั้งหมดมาหาค่า accuracy:
```
# Calculate the accuracy
accuracy <- sum(diag(cm)) / sum(cm)
cat("Accuracy:", round(accuracy, 2))
```
ผลลัพธ์:
```
Accuracy: 0.93
```
จะเห็นได้ว่า NB classifier ที่เป็น model ไร้เดียงสาของเรา มีความแม่นยำในการทำนายสูงถึง 93%

😺 GitHub

ดู code ทั้งหมดในบทความนี้ได้ที่ GitHub

📃 References
✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

🙋 ใครที่กำลังเรียนจิตวิทยาหรือทำงานสายจิตวิทยา และเบื่อที่ต้องใช้ software ราคาแพงอย่าง SPSS และ Excel เพื่อทำข้อมูล

💪 ผมขอแนะนำ R Book for Psychologists หนังสือสอนใช้ภาษา R เพื่อการวิเคราะห์ข้อมูลทางจิตวิทยา ที่เขียนมาเพื่อนักจิตวิทยาที่ไม่เคยมีประสบการณ์เขียน code มาก่อน

ในหนังสือ เราจะปูพื้นฐานภาษา R และพาไปดูวิธีวิเคราะห์สถิติที่ใช้บ่อยกัน เช่น:
- Correlation
- t-tests
- ANOVA
- Reliability
- Factor analysis
🚀 เมื่ออ่านและทำตามตัวอย่างใน R Book for Psychologists ทุกคนจะไม่ต้องพึง SPSS และ Excel ในการทำงานอีกต่อไป และสามารถวิเคราะห์ข้อมูลด้วยตัวเองได้ด้วยความมั่นใจ

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists
Share this:
X
Facebook
Like Loading…
2025-04-10
วิธีสร้าง KNN model ใน 5 ขั้นตอนด้วย class package ในภาษา R – ตัวอย่างการทำนายประเภทดอกไม้ใน iris dataset
ในบทความนี้ เราจะมาทำความรู้กันกับ KNN model และวิธีสร้าง KNN model ด้วย class package ในภาษา R กัน
KNN หรือ K-Nearest Neighbors เป็น machine learning algorithm ประเภท supervised learning ซึ่งใช้ได้กับงานประเภท classification (จัดกลุ่ม) และ regression (ทำนายแนวโน้ม)

การทำงานของ KNN ตรงไปตรงมา คือ ทำนายข้อมูลโดยการดูจากข้อมูลใกล้เคียง จำนวน k ข้อมูล (เราสามารถกำหนดค่า k เองได้)

ยกตัวอย่างเช่น เราต้องการใช้ KNN ทำนายราคาหุ้น KNN จะดูราคาหุ้นที่อยู่ใกล้เคียงกับหุ้นที่เราต้องการทำนายราคา เช่น 5 หุ้นใกล้เคียง โดย:
- 2/5 หุ้น = ราคาลง
- 3/5 หุ้น = ราคาขึ้น
KNN จะทำนายว่า หุ้นในใจของเรา จะราคาขึ้น เพราะหุ้นใกล้เคียงส่วนใหญ่ปรับราคาขึ้น เป็นต้น

💻 KNN in R With class Package

เราสามารถสร้าง KNN model ในภาษา R ด้วย class package โดยมี 5 ขั้นตอนดังนี้:
1. Get class package
2. Load the dataset
3. Prepare the data
4. Train a KNN model
5. Evaluate the model
เราไปดูตัวอย่างการสร้าง KNN model เพื่อทำนายประเภทดอกไม้ ด้วย iris dataset กัน

1️⃣ Step 1. Get class Package

ในขั้นแรก ให้เราติดตั้งและเรียกใช้งาน class package

ติดตั้ง (ทำครั้งเดียว):
```
# Install
install.packages("class")
```
เรียกใช้งาน (ทำทุกครั้งที่เริ่ม session ใหม่):
```
# Load
library(class)
```
2️⃣ Step 2. Get the Dataset

ในขั้นที่สอง ให้เราโหลด dataset ที่จะใช้งาน

โดยในบทความนี้ เราจะใช้ iris ที่มีข้อมูลดอกไม้ 3 ชนิด ได้แก่:
- Iris setosa
- Iris virginica
- Iris versicolor
นอกจากประเภทดอกไม้ (species) iris dataset ยังประกอบด้วยความกว้าง (width) และความยาว (length) ของกลีบดอกชั้นนอก (sepal) และชั้นใน (petal)

เนื่องจาก iris เป็น built-in dataset ใน ภาษา R เราสามารถโหลดข้อมูลด้วยคำสั่ง data() ได้:
```
# Load
data(iris)
```
เมื่อโหลดแล้ว เราสามารถดูตัวอย่าง dataset ได้ด้วย head():
```
# Preview
head(iris)
```
ผลลัพธ์:
```
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa
```
3️⃣ Step 3. Prepare Data

ในขั้นที่สาม เราจะเตรียมข้อมูล เพื่อใช้สร้าง KNN model

โดยสำหรับ iris dataset เราจะทำ 3 อย่าง คือ:
1. Normalise data
2. Split data
3. Separate features from the label
.

📉 Normalise Data

เนื่องจาก KNN model อ้างอิงระยะห่างระหว่างข้อมูลในการทำนายผลลัพธ์ ข้อมูลที่มีช่วงกว้าง (เช่น 0-100) จะมีน้ำหนักกว่าข้อมูลที่มีช่วงแคบ (เช่น 0-10) และทำให้ model ของเรา bias ได้

เพื่อป้องกันความเสี่ยงนี้ เราต้อง normalise ข้อมูล หรือการทำให้ข้อมูลมีช่วงมีช่วงข้อมูลเดียวกัน

การทำ normalisation มีหลายวิธี แต่วิธีที่นิยมใช้กับ KNN คือ min-max normalisation ซึ่งเปลี่ยนช่วงข้อมูลให้อยู่ในช่วง 0-1 ผ่านสูตร:
```
(X - Min) / (Max - Min)
```
- X คือ ข้อมูล
- Min คือ ค่าที่น้อยที่สุด
- Max คือ ค่าที่มากที่สุด
ในภาษา R เราไม่มี function ที่ใช้ทำ min-max normalisation ได้โดยตรง เราเลยต้องเขียน function ขึ้นมาใช้เอง:
```
# Define a function for normalisation
normalise <- function(x) {
  return((x - min(x)) / (max(x) - min(x)))
}
```
จากนั้น ให้เราใช้ function นี้กับ columns ที่เป็นตัวเลข (numeric) เช่น Sepal.Length แล้วเก็บผลลัพธ์ไว้ใน data frame ใหม่ ชื่อ iris_normalises:
```
# Apply the function to the dataset
iris_normalised <- as.data.frame(lapply(iris[, 1:4],
                                        normalise))
```
เนื่องจากตอนที่เราใช้ function เราตัด column ที่ไม่ใช่ตัวเลขออก ทำให้ iris_normalised ไม่มี column Species

เราสามารถใส่ Species กลับเข้าไป เพื่อใช้ในขั้นตอนถัดไป แบบนี้:
```
# Add species column back into the data frame
iris_normalised$Species <- iris$Species
```
สุดท้าย เราเช็กผลลัพธ์ด้วยการเรียกดูสถิติของ iris_normalised:
```
# Check the results
summary(iris_normalised)
```
ผลลัพธ์:
```
  Sepal.Length     Sepal.Width    
 Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.2222   1st Qu.:0.3333  
 Median :0.4167   Median :0.4167  
 Mean   :0.4287   Mean   :0.4406  
 3rd Qu.:0.5833   3rd Qu.:0.5417  
 Max.   :1.0000   Max.   :1.0000  
  Petal.Length     Petal.Width     
 Min.   :0.0000   Min.   :0.00000  
 1st Qu.:0.1017   1st Qu.:0.08333  
 Median :0.5678   Median :0.50000  
 Mean   :0.4675   Mean   :0.45806  
 3rd Qu.:0.6949   3rd Qu.:0.70833  
 Max.   :1.0000   Max.   :1.00000  
       Species  
 setosa    :50  
 versicolor:50  
 virginica :50  
```
จะเห็นว่า:
1. Columns ที่เป็นตัวเลข มีช่วงอยู่ระหว่าง 0 และ 1
2. เรายังมี column Species อยู่
.

🪓 Split Data

ในการสร้าง KNN model เราควรแบ่ง dataset ที่มีเป็น 2 ส่วน คือ:
1. Training set: ใช้สำหรับสร้าง model
2. Test set: ใช้สำหรับประเมิน model
เราเริ่มแบ่งข้อมูลด้วยการสุ่ม row index ที่จะอยู่ใน training set:
```
# Set seed for reproducibility
set.seed(2025)

# Create a training index
train_index <- sample(1:nrow(iris_normalised),
                      0.7 * nrow(iris_normalised))
```
จากนั้น subset ข้อมูลด้วย row index ที่สุ่มไว้:
```
# Split the data
train_set <- iris_normalised[train_index, ]
test_set <- iris_normalised[-train_index, ]
```
.

🏷️ Separate Features From Label

ขั้นตอนสุดท้ายในการเตรียมข้อมูล คือ แยก features หรือ X (columns ที่จะใช้ทำนาย) ออกจาก label หรือ Y (สิ่งที่ต้องการทำนาย):
```
# Separate features from label

## Training set
train_X <- train_set[, 1:4]
train_Y <- train_set$Species

## Test set
test_X <- test_set[, 1:4]
test_Y <- test_set$Species
```
4️⃣ Step 4. Train a KNN Model

ขั้นที่สี่เป็นขั้นที่เราสร้าง KNN model ขึ้นมา โดยเรียกใช้ knn() จาก class package

ทั้งนี้ knn() ต้องการ input 3 อย่าง:
1. train: fatures จาก training set
2. test: feature จาก test set
3. cl: label จาก training set
4. k: จำนวนข้อมูลใกล้เคียงที่จะใช้ทำนายผลลัพธ์
```
# Train a KNN model
pred <- knn(train = train_X,
            test = test_X,
            cl = train_Y,
            k = 5)
```
ในตัวอย่าง เรากำหนด k = 5 เพื่อทำนายผลลัพธ์โดยดูจากข้อมูลที่ใกล้เคียง 5 อันดับแรก

5️⃣ Step 5. Evaluate the Model

หลังจากได้ model แล้ว เราประเมินประสิทธิภาพของ model ในการทำนายผลลัพธ์ ซึ่งเราทำได้ง่าย ๆ โดยคำนวณ accuracy หรือสัดส่วนของข้อมูลที่ model ตอบถูกต่อจำนวนข้อมูลทั้งหมด:
```
Accuracy = Correct predictions / Total predictions
```
ในภาษา R ให้เราสร้าง confusion matrix หรือ matrix แสดงจำนวนคำตอบที่ถูกและผิด ด้วย table() function:
```
# Create a confusion matrix
cm <- table(Predicted = predictions,
            Actual = test_Y)

# Print the matrix
print(cm)
```
ผลลัพธ์:
```
            Actual
Predicted    setosa versicolor virginica
  setosa         15          0         0
  versicolor      0         17         0
  virginica       0          2        11
```
Note: จะเห็นได้ว่า model เราตอบถูกเป็นส่วนใหญ่ (ราว 90%)

จากนั้น คำนวณ accuracy ด้วย sum() และ diag():
```
# Calculate accuracy
acc <- sum(diag(cm)) / sum(cm)

# Print accuracy
cat("Accuracy:", round(acc, 2))
```
ผลลัพธ์:
```
Accuracy: 0.96
```
จากผลลัพธ์ เราจะเห็นว่า model มีความแม่นยำสูงถึง 96%

🍩 Bonus: Fine-Tuning

ในบางครั้ง ค่า k ที่เราตั้งไว้ อาจไม่ได้ทำให้เราได้ KNN model ที่ดีที่สุด

แทนที่เราจะแทนค่า k ใหม่ไปเรื่อย ๆ เราสามารถใช้ for loop เพื่อหาค่า k ที่ทำให้เราได้ model ที่ดีที่สุดได้

ให้เราเริ่มจากสร้าง vector ที่มีค่า k ที่ต้องการ:
```
# Create a set of k values
k_values <- 1:20
```
และ vector สำหรับเก็บค่า accuracy ของค่า k แต่ละตัว:
```
# Createa a vector for accuracy results
accuracy_results <- numeric(length(k_values))
```
แล้วใช้ for loop เพื่อหาค่า accuracy ของค่า k:
```
# For-loop through the k values
for (i in seq_along(k_values)) {
  
  ## Set the k value
  k <- k_values[i]
  
  ## Create a KNN model
  predictions <- knn(train = train_X,
                     test = test_X,
                     cl = train_Y,
                     k = k)
  
  ## Create a confusion matrix
  cm <- table(Predicted = predictions,
              Actual = test_Y)
  
  ## Calculate accuracy
  accuracy_results[i] <- sum(diag(cm)) / sum(cm)
}

# Find the best k and the corresponding accuracy
best_k <- k_values[which.max(accuracy_results)]
best_accuracy <- max(accuracy_results)

# Print best k and accuracy
cat(paste("Best k:", best_k),
    paste("Accuracy:", round(best_accuracy, 2)),
    sep = "\n")
```
ผลลัพธ์:
```
Best k: 12
Accuracy: 0.98
```
แสดงว่า ค่า k ที่ดีที่สุด คือ 12 โดยมี accuracy เท่ากับ 98%

นอกจากนี้ เรายังสามารถสร้างกราฟ เพื่อช่วยทำความเข้าใจผลของค่า k ต่อ accuracy:
```
# Plot the results
plot(k_values,
     accuracy_results,
     type = "b",
     pch = 19,
     col = "blue",
     xlab = "Number of Neighbors (k)",
     ylab = "Accuracy",
     main = "KNN Model Accuracy for Different k Values")
grid()
```
ผลลัพธ์:

จะเห็นได้ว่า k = 12 ให้ accuracy ที่ดีที่สุด และ k = 20 ให้ accuracy ต่ำที่สุด ส่วนค่า k อื่น ๆ ให้ accuracy ในช่วง 93 ถึง 96%

😺 GitHub

ดู code ทั้งหมดในบทความนี้ได้ที่ GitHub

📃 References
✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

🙋 ใครที่กำลังเรียนจิตวิทยาหรือทำงานสายจิตวิทยา และเบื่อที่ต้องใช้ software ราคาแพงอย่าง SPSS และ Excel เพื่อทำข้อมูล

💪 ผมขอแนะนำ R Book for Psychologists หนังสือสอนใช้ภาษา R เพื่อการวิเคราะห์ข้อมูลทางจิตวิทยา ที่เขียนมาเพื่อนักจิตวิทยาที่ไม่เคยมีประสบการณ์เขียน code มาก่อน

ในหนังสือ เราจะปูพื้นฐานภาษา R และพาไปดูวิธีวิเคราะห์สถิติที่ใช้บ่อยกัน เช่น:
- Correlation
- t-tests
- ANOVA
- Reliability
- Factor analysis
🚀 เมื่ออ่านและทำตามตัวอย่างใน R Book for Psychologists ทุกคนจะไม่ต้องพึง SPSS และ Excel ในการทำงานอีกต่อไป และสามารถวิเคราะห์ข้อมูลด้วยตัวเองได้ด้วยความมั่นใจ

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists
Share this:
X
Facebook
Like Loading…
2025-04-03
แนะนำ dtplyr: package เชื่อม dplyr และ data.table เพื่อการทำ data wrangling กับข้อมูลขนาดใหญ่ที่ง่ายและรวดเร็ว — ตัวอย่างจาก mtcars dataset
ในบทความนี้ เรามาความรู้จักกับ package ในภาษา R ที่เรียกว่า dtplyr กัน
🏎️ dtplyr คืออะไร?

dtplyr เป็น package ที่เชื่อม 2 libraries ยอดนิยมสำหรับ data wrangling เข้าด้วยกัน ได้แก่:

.

Library #1. dplyr ซี่งมี verb-based syntax ที่ใช้งานง่าย

เช่น จาก mtcars dataset เลือกรถเกียร์ manual (am == 1) ที่มี miles/gallon (mpg) ตั้งแต่ 20 ขึ้นไป โดยเรียงจากมากไปน้อย:
```
# dplyr
mtcars |>
  filter(am == 1 & mpg >= 20) |>
  select(model, mpg) |>
  arrange(desc(mpg))
```
Note: mtcars ถูกปรับให้มี column ชื่อ model (ดู code ในการเพิ่ม column ได้ใน GitHub)

.

Library #2. data.table ซึ่งมี syntax ที่เข้าใจยากกว่า:
```
# data.table
mtcars_dt[am == 1 & mpg >= 20, .(model, mpg)][order(-mpg)]
```
Note: mtcars_dt เป็น dataset ที่ถูกเปลี่ยนจาก data.frame เป็น data.table object เพื่อใช้งานกับ data.table

แต่มีจุดเด่น คือ ประมวลผลได้เร็ว เหมาะกับการทำงานกับ data ขนาดใหญ่

.

dtplyr เชื่อม dplyr เข้ากับ data.table โดยการแปล dplyr syntax ให้เป็น data.table syntax ทำให้เราทำงานได้ง่ายและรวดเร็ว

🧑‍💻 วิธีใช้ dytplyr

เราสามารถใช้ dtplyr ได้ใน 3 ขั้นตอนง่าย ๆ:
1. Install and load dtplyr
2. Create a lazy data.table
3. Execute the syntax
.

1️⃣ ขั้นที่ 1. Install & Load dtplyr

ในการใช้งาน dtplyr เราต้องเรียกใช้ dplyr ด้วยเสมอ

ดังนั้น ในการติดตั้ง เราต้องติดตั้งทั้ง dtplyr และ dplyr:
```
# Install
install.packages("dplyr")
install.packages("dtplyr")
```
เมื่อติดตั้งแล้ว ให้เราเรียกใช้งาน dtplyr และ dplyr ทุกครั้งที่เริ่ม session ในการทำงาน:
```
# Load
library(dplyr)
library(dtplyr)
```
Note: คู่มือการใช้ dtplyr แนะนำให้โหลด data.table ด้วย เพื่อใช้ functions ของ data.table ในการทำงาน

.

2️⃣ ขั้นที่ 2. Create a Lazy data.table

หลังเรียกใช้งาน dtplyr ให้เราสร้าง dataset ที่เป็น lazy data.table object ขึ้นมา เพื่อช่วยให้ dtplyr เปลี่ยนคำสั่งของ dplyr เป็นคำสั่งของ data.table ได้:
```
# Convert mtcars dataset into lazy data.table object
mtcars_ldt <- lazy_dt(mtcars)
```
.

3️⃣ ขั้นที่ 3. Execute the Syntax

ในขั้นสุดท้าย เราสามารถเริ่มทำงานโดยใช้ syntax ของ dplyr ได้เลย:
```
# Execute the syntax
mtcars_ldt |>
  filter(am == 1 & mpg >= 20) |>
  select(model, mpg) |>
  arrange(desc(mpg)) |>
  as.data.table()
```
ผลลัพธ์:
```
             model   mpg
            <char> <num>
 1: Toyota Corolla  33.9
 2:       Fiat 128  32.4
 3:    Honda Civic  30.4
 4:   Lotus Europa  30.4
 5:      Fiat X1-9  27.3
 6:  Porsche 914-2  26.0
 7:     Datsun 710  22.8
 8:     Volvo 142E  21.4
 9:      Mazda RX4  21.0
10:  Mazda RX4 Wag  21.0
```
จะสังเกตว่า ใน code บรรทัดสุดท้าย เราเพิ่ม as.data.table() เข้ามาเพื่อบอกให้ dtplyr รู้ว่า เราเขียน syntax เสร็จแล้ว และพร้อมให้แปลงให้เป็น data.table syntax

ทั้งนี้ ถ้าเราเขียนโดยไม่มี as.data.table():
```
# Execute the syntax, without as.data.table()
mtcars_ldt |>
  filter(am == 1 & mpg >= 20) |>
  select(model, mpg) |>
  arrange(desc(mpg))
```
เราจะได้ผลลัพธ์ที่เป็นแค่ preview ที่ไม่สามารถนำไปใช้งานต่อได้:
```
  model            mpg
  <chr>          <dbl>
1 Toyota Corolla  33.9
2 Fiat 128        32.4
3 Honda Civic     30.4
4 Lotus Europa    30.4
5 Fiat X1-9       27.3
6 Porsche 914-2   26  
# ℹ 4 more rows
# ℹ Use `print(n = ...)` to see more rows
```
ดังนั้น ถ้าเราต้องการผลลัพธ์ทั้งหมด เราต้องใช้ as.data.table() ในการรันคำสั่งทุกครั้ง

Note: เราสามารถใช้คำสั่งอื่น ๆ แทน as.data.table() ได้ โดยจะให้ผลลัพธ์ที่แตกต่างกันไปดังนี้:

Function ผลลัพธ์ที่ส่งกลับมา
as.data.table() data.table
as.data.frame() data.frame
as.tibble() tibble
collect() tibble
pull() Single column

😺 GitHub

ดู code ทั้งหมดในบทความนี้ได้ที่ GitHub

📃 References
✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

🙋 ใครที่กำลังเรียนจิตวิทยาหรือทำงานสายจิตวิทยา และเบื่อที่ต้องใช้ software ราคาแพงอย่าง SPSS และ Excel เพื่อทำข้อมูล

💪 ผมขอแนะนำ R Book for Psychologists หนังสือสอนใช้ภาษา R เพื่อการวิเคราะห์ข้อมูลทางจิตวิทยา ที่เขียนมาเพื่อนักจิตวิทยาที่ไม่เคยมีประสบการณ์เขียน code มาก่อน

ในหนังสือ เราจะปูพื้นฐานภาษา R และพาไปดูวิธีวิเคราะห์สถิติที่ใช้บ่อยกัน เช่น:
- Correlation
- t-tests
- ANOVA
- Reliability
- Factor analysis
🚀 เมื่ออ่านและทำตามตัวอย่างใน R Book for Psychologists ทุกคนจะไม่ต้องพึง SPSS และ Excel ในการทำงานอีกต่อไป และสามารถวิเคราะห์ข้อมูลด้วยตัวเองได้ด้วยความมั่นใจ

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists
Share this:
X
Facebook
Like Loading…
2025-03-06

Function	ผลลัพธ์ที่ส่งกลับมา
`as.data.table()`	`data.table`
`as.data.frame()`	`data.frame`
`as.tibble()`	`tibble`
`collect()`	`tibble`
`pull()`	Single column

วิธีใช้ ggplot2 เพื่อสร้างกราฟอย่างมืออาชีพระดับโลก แบบ BBC และ Financial Times ในภาษา R — ตัวอย่างการสำรวจข้อมูลเพนกวินจาก palmerpenguins

ggplot2 เป็น package สำหรับ data visualisation ในภาษา R และเป็นเครื่องมือสร้างกราฟที่มืออาชีพนิยม ตั้งแต่นักวิจัยในการตีพิมพ์ผลงาน ไปจนถึงสำนักข่าวระดับโลกอย่าง BBC และ Financial Times

ggplot2 มีจุดเด่น 4 ข้อ:

ใช้งานง่าย
สร้างกราฟได้หลากหลาย
ปรับแต่งกราฟได้ดังใจ
ได้กราฟที่ดูดีและมีคุณภาพ

ในบทความนี้ เราจะมาดูวิธีใช้ ggplot2 เพื่อสร้างกราฟแบบมืออาชีพกัน:

ggplot2 syntax
Basic plotting
Plot customisations

ถ้าพร้อมแล้ว ไปเริ่มกันเลย

🔤 gg for “Grammar of Graphics”

gg ใน ggplot2 ย่อมาจาก “Grammar of Graphics”

หนังสือ Grammar of Graphics ของ Leland Wilkinson บน Amazon

Grammar of Graphics เป็นแนวคิดที่มองกราฟเป็นเหมือนภาษา คือ มีโครงสร้างและองค์ประกอบที่ตายตัว ซึ่งเมื่อเรานำมารวมกัน เราก็จะได้กราฟที่ต้องการขึ้นมา

โดยกราฟใน ggplot2 ประกอบด้วย 7 ส่วน หรือ layers ได้แก่:

No.	Layer	Description
1	Data	ชุดข้อมูลสำหรับสร้างกราฟ
2	Aesthetics	จับคู่ข้อมูลกับกราฟ (เช่น แกน x และ y)
3	Geometric objects	ประเภทกราฟ (เช่น กราฟเส้น กราฟแท่ง)
4	Facets	สร้างกราฟย่อย
5	Statistical transformations	วิเคราะห์ข้อมูล (เช่น หาค่าเฉลี่ย)
6	Coordinates	แกนในการสร้างกราฟ
7	Theme	หน้าตาของกราฟ (เช่น สีพื้นหลัง)

ในการทำงาน เรามักจะเรียกใช้งาน 3 layers แรกเป็นหลัก ได้แก่:

Data
Aesthetics
Geometric objects

🏁 Getting Started With ggplot2

ในการเริ่มต้นใช้งาน ggplot2 เราต้องทำ 3 อย่างก่อน คือ:

ข้อที่ 1. ติดตั้ง ggplot2 บน environment ของเรา:

install.packages("ggplot2")

Note: ถ้าใครเคยติดตั้งแล้ว สามารถข้ามไปขั้นถัดไปได้เลย

ข้อที่ 2. เรียกใช้งาน ggplot2:

library(ggplot2)

Note: เราต้องเรียกใช้งาน ggplot2 ทุกครั้งที่เริ่ม session ใหม่

ข้อที่ 3. โหลด dataset ที่เราจะใช้สร้างกราฟ

สำหรับบทความนี้ เราจะใช้ penguins dataset ที่มีข้อมูลของเพนกวิน 3 สายพันธุ์ (เช่น สปีชีส์ น้ำหนัก ความยาวปีก) กัน

เราสามารถโหลด dataset ได้ดังนี้:

install.packages("palmerpenguins")
library(palmerpenguins)

จากนั้น เราสามารถ preview ข้อมูลได้ด้วย head():

head(penguins)

ผลลัพธ์:

Note: เราสามารถอ่านคู่มือ penguins ได้ด้วยคำสั่ง ?penguins

เมื่อทำครบทั้ง 3 ขั้นตอน เราก็พร้อมที่จะสร้างกราฟใน ggplot2 แล้ว

✍️ Basic Syntax

ก่อนไปดูการสร้างกราฟ เรามาดู syntax ของ ggplot2 กันก่อน:

ggplot(data, aes(x, y, other)) +
	geom_*() +
	...

ggplot() คือ การเรียกใช้งาน ggplot2
data คือ ชุดข้อมูลในการสร้างกราฟ
aes() คือ ส่วนจับคู่ข้อมูลกับลักษณะของกราฟ
- x, y คือ ข้อมูลที่แสดงบนแกน x และ y
- other คือ ข้อมูลที่จะแสดงผ่านส่วนอื่น ๆ ของกราฟ เช่น สี ขนาด รูปทรง
geom_* คือ ประเภทกราฟ
… คือ function อื่น ๆ ในการตั้งค่ากราฟ (เช่น theme, facet)

📊 Basic Plotting: Data, Aesthetics, & Geom

ในการสร้างกราฟ เรามี 3 input เบื้องต้นที่เราต้องกำหนด ได้แก่:

No.	Input	Description
1	Data	Dataset ในการสร้างกราฟ
2	Aesthetics	ลักษณะของกราฟที่ใช้แสดงข้อมูล (เช่น แกน x และ y)
3	Geom	ประเภทกราฟ (เช่น กราฟเส้น กราฟแท่ง)

Note: ทั้งสามอย่างสะท้อนถึง 3 layers แรกของ ggplot2

ยกตัวอย่างเช่น

สร้าง scatter plot ที่แสดงความสัมพันธ์ระหว่างน้ำหนักตัว (body_mass_g) และความยาวปีก (flipper_length_mm) ของเพนกวิน:

ggplot(penguins, aes(x = body_mass_g,
                     y = flipper_length_mm)) +
  geom_point()

ผลลัพธ์:

🎨 Aesthetics

นอกจากแกน x และ y เราแสดงข้อมูลผ่านลักษณะอื่น ๆ ของกราฟได้ เช่น:

No.	Parametre	Description
1	`color`	สีขอบรูปทรง
2	`fill`	สีในรูปทรง
3	`shape`	รูปทรง
4	`size`	ขนาด
5	`alpha`	ความโปร่งใส

ยกตัวอย่างเช่น

เพิ่มสปีชีส์เข้าไปใน scatter plot:

ggplot(penguins, aes(x = body_mass_g,
                     y = flipper_length_mm,
                     color = species)) +
  geom_point()

ผลลัพธ์:

Note: เราใช้ color แทน fill เพราะ จุด default ของ scatter plot ไม่รองรับการเติมสี

⏹️ Geom

เราสามารถเปลี่ยนประเภทกราฟได้ โดยเปลี่ยน geom_*() เช่น:

No.	Geom	Graph
1	`geom_histogram()`	Histogram
2	`geom_boxplot()`	Box plot
3	`geom_line()`	Line plot
4	`geom_col()`	Bar plot
5	`geom_density()`	Density plot

Note: ggplot2 มีรูปแบบกราฟกว่า 40+ แบบให้เลือก เราสามารถดูรูปแบบกราฟต่าง ๆ ได้ที่ Function reference และดูตัวอย่างกราฟได้ที่ The R Gallery

นอกจากนี้ เรายังสามารถปรับแต่งหน้าตาของกราฟได้ด้วย parametre ต่าง ๆ ใน geom_*() เช่น:

No.	Parametre	Description
1	`color`	สีขอบรูปทรง
2	`fill`	สีในรูปทรง
3	`shape`	รูปทรง
4	`size`	ขนาด
5	`alpha`	ความโปร่งใส

Note:

จะสังเกตว่า parametre เหล่านี้ (เรียกว่า attribute) เหมือนกับ parametres ของ aes()
ความแตกต่างอยู่ที่ parametre ใน geom_*() จะคงที่ ในขณะที่ parametre ใน aes() จะเปลี่ยนแปลงตามข้อมูล

ยกตัวอย่างเช่น

เปลี่ยนจุดข้อมูลให้เป็นกล่องสี่เหลี่ยมใส:

ggplot(penguins, aes(x = body_mass_g,
                     y = flipper_length_mm,
                     color = species)) +
  geom_point(shape = 22)

ผลลัพธ์:

จะเห็นได้ว่า จุดข้อมูลเปลี่ยนเป็นสี่เหลี่ยมเหมือนกันหมด และไม่เปลี่ยนตามประเภทเพนกวิน เพราะเรากำหนด argument ใน geom_*()

Note: เราสามารถดู argument ของ shape และ parametres อื่น ๆ ได้ที่ Aesthetic specifications

🔧 More Customisations

จนถึงตอนนี้ เรารู้วิธีสร้างและปรับแต่งกราฟเบื้องต้นกันแล้ว

เรามาดู 3 วิธีในการปรับแต่งเพิ่มเติม เพื่อให้กราฟของเราดูเป็นมืออาชีพกัน:

Theme
Text
Facet

🖼️ Theme

Theme ใน ggplot2 หมายถึง หน้าตากราฟที่ไม่เกี่ยวข้องกับข้อมูล เช่น สีแกน x และ y และสีพื้นหลัง

เราสามารถปรับ theme ได้โดยการเรียกใช้ theme_*() เช่น:

ggplot(penguins, aes(x = body_mass_g,
                     y = flipper_length_mm,
                     color = species)) +
  geom_point() +
  
  # Apply classic theme
  theme_classic()

ผลลัพธ์:

จะเห็นว่า กราฟของเราดูสะอาดตามากขึ้น เมื่อไม่มี gridline

Note:

ggplot2 มี built-in themes ให้เลือกใช้ 5 แบบหลัก ได้แก่:

No.	Theme	Background	Gridline
1	`theme_gray()`	สีเทา	สีขาว
2	`theme_classic()`	สีขาว	None
3	`theme_bw()`	สีขาว	สีเทา
4	`theme_light()`	สีขาว	สีเทา
5	`theme_dark()`	สีเทาเข้ม	สีเทาอ่อน

theme_gray() เป็น default theme ของ ggplot2
ดู built-in themes อื่น ๆ ได้ที่ Complete themes
เราสามารถเลือกใช้ themes จาก packages อื่นใน R ได้ เช่น ggthemes และ bbcplot

🔤 Text

ข้อความเป็นส่วนสำคัญของกราฟ ซึ่งสามารถช่วยให้คนที่ดูกราฟของเราเข้าใจกราฟได้ง่ายขึ้น

ใน ggplot2 เราสามารถปรับแต่งข้อความได้ 3 วิธี:

No.	Syntax	Description
1	`theme()`	กำหนดขนาด และ typeface ของข้อความ
2	`labs()`	เพิ่มชื่อกราฟ ชื่อแกน x, y และ legend
3	`annotate()`	เพิ่มโน้ตในกราฟ

ตัวอย่างเช่น:

ggplot(penguins, aes(x = body_mass_g,
                       y = flipper_length_mm,
                       color = species)) +
  geom_point() +
  theme_classic() +
  
  # Adjust text size
  theme(plot.title = element_text(size = 16, face = "bold"),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14, face = "bold")) +
  
  # Add a title, labels, and a legend
  labs(title = "Penguin Body Mass vs. Flipper Length",
       x = "Body Mass (g)",
       y = "Flipper Length (mm)",
       color = "Penguin Species") +
  
  # Add annotation
  annotate("text",
           x = 3000,
           y = 225,
           label = "Larger penguins tend to \\n have longer flippers",
           size = 5,
           color = "gray",
           hjust = 0)

ผลลัพธ์:

จะเห็นได้ว่า กราฟของเรามีชื่อกราฟ คำประกอบแกน x, y และ legend รวมทั้งโน้ตที่ช่วยในการอ่านกราฟ “Larger penguins tend to have longer flippers” เพิ่มขึ้นมา

ในบางครั้ง เราอาจต้องการสร้างกราฟย่อย ซึ่งเราสามารถทำได้ด้วย 2 คำสั่งใน ggplot2:

No.	Function	Dividing Factor
1	`facet_wrap()`	1 categorical variable
2	`facet_grid()`	2 categorical variables

ตัวอย่างเช่น

ใช้ facet_wrap() แบ่งกราฟตามชนิดเพนกวิน:

ggplot(penguins, aes(x = body_mass_g,
                     y = flipper_length_mm,
                     color = species)) +
  geom_point() +
  theme_bw() +
  
  # Use facet_wrap()
  facet_wrap(~species)

ผลลัพธ์:

ใช้ facet_grid() แบ่งตามเพศและชนิดของเพนกวิน:

ggplot(penguins, aes(x = body_mass_g,
                     y = flipper_length_mm,
                     color = species)) +
  geom_point() +
  theme_bw() +
  
  # Use facet_grid()
  facet_grid(sex~species)

ผลลัพธ์:

Note: ในกราฟ เรามี row ที่ 3 (NA) เพราะบาง record ไม่มีข้อมูลเพศของเพนกวิน

🔥 Summary

ในบทความนี้ เราได้เรียนรู้วิธีใช้ ggplot2 เพื่อสร้างกราฟอย่างมืออาชีพ

Syntax ของ ggplot2
การตั้งค่า aesthetics
- x
- y
- Aesthetics อื่น ๆ
การตั้งค่า geom
- geom_*()
- paramatres ของ geom_*()
การปรับแต่งกราฟ
- Theme
- Text
- Facet

📚 Learn More About ggplot2

😺 GitHub

สำหรับผู้ที่สนใจ สามารถดู code ทั้งหมดในบทความนี้ได้ที่ GitHub

📖 Read More About ggplot2

📰 Cheat Sheets

📃 References

Articles:

Courses

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

Correlation
t-tests
ANOVA
Reliability
Factor analysis

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists

2025-02-20

data.table: แนะนำ 3 ส่วนวิธีใช้งาน package ทรงพลัง เพื่อทำงานกับข้อมูลขนาดใหญ่อย่างมีประสิทธิภาพและรวดเร็ว ในภาษา R แบบครบ — ตัวอย่างการทำงานกับข้อมูลเที่ยวบินจาก flights dataset

R ได้ชื่อว่าเป็น “statistical programming language” เพราะออกแบบมาเพื่อทำงานกับ data โดยเฉพาะ

ภาษา R มี packages มากมายที่รองรับการทำงานกับ data ในรูปแบบตาราง (tabular data) หรือ data frame อย่างตัวอย่างในภาพ:

หนึ่งใน packages ที่ได้รับความนิยมในการทำงานกับ data frame โดยเฉพาะงาน data science ได้แก่ data.table

data.table เป็น package สำหรับ data manipulation ที่มีจุดเด่น 3 อย่าง คือ:

ใช้งานง่าย (แม้จะมีการเขียนที่แตกต่างจาก functions ทั่วไป แต่มี syntax ที่ตายตัว)
ประมวลผลเร็วและประหยัด resource เพราะ data.table optimises code ที่เราเขียน
รองรับการทำงานกับ data ขนาดใหญ่ (เช่น data ที่ใช้ RAM ขนาด 10 ถึง 100 GB ในการประมวลผล)

ในบทความนี้ เราจะมาดูวิธีการใช้งาน data.table แบบครบจบใน 3 ส่วนกัน:

การใช้งาน data.table เบื้องต้น (syntax และ arguments)
การใช้งาน data.table ขั้นสูง (combining และ chaining)
Special symbols ใน data.table

ถ้าพร้อมแล้ว มาเริ่มกันเลย

🏁 Get Started With data.table

ในการเริ่มต้นใช้งาน data.table เราจะต้องทำ 2 อย่าง ได้แก่:

ติดตั้ง data.table (ทำครั้งแรกครั้งเดียว)
โหลด data.table (ทำทุกครั้งที่เริ่ม session ใหม่)

# Install data.table
install.packages("data.table")

# Load data.table
library(data.table)

เมื่อทำทั้ง 2 อย่างนี้แล้ว เราก็พร้อมที่จะใช้งาน data.table กันแล้ว

✈️ Example Dataset: flights

1️⃣ Intro to flights

ในบทความนี้ เราจะมาดูตัวอย่างการใช้ data.table กับ flights ซึ่งมีข้อมูลเที่ยวบินที่ออกจาก New York City ในปี ค.ศ. 2014 เช่น:

วันที่
สนามบินต้นทางและปลายทาง
ระยะเวลาที่เที่ยวบิน delay
ระยะเวลาบิน
ระยะทาง

Note: ดูรายละเอียดทั้งหมดของ dataset ได้ที่ GitHub

2️⃣ Get flights

เราสามารถเริ่มใช้งาน flights ใน 2 ขั้นตอน:

ขั้นที่ 1. Download ไฟล์จาก link โดยไฟล์ที่ได้จะเป็น CSV (comma-separated values):

ขั้นที่ 2. Import ข้อมูลเข้าใน R ด้วย fread() (”fast read”) ซึ่งเป็น function สำหรับโหลดข้อมูลของ data.table:

# Import the dataset
flights <- fread("flights14.csv")

เราสามารถดูตัวอย่างข้อมูลได้ด้วย head():

# Preview the dataset
head(flights)

ผลลัพธ์:

    year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
   <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
1:  2014     1     1        14        13      AA    JFK    LAX      359     2475     9
2:  2014     1     1        -3        13      AA    JFK    LAX      363     2475    11
3:  2014     1     1         2         9      AA    JFK    LAX      351     2475    19
4:  2014     1     1        -8       -26      AA    LGA    PBI      157     1035     7
5:  2014     1     1         2         1      AA    JFK    LAX      350     2475    13
6:  2014     1     1         4         0      AA    EWR    LAX      339     2454    18

✍️ Syntax: DT[i, j, by]

การใช้งาน data.table ประกอบด้วย 4 ส่วน ดังนี้:

DT[i, j, by]

DT คือ dataset ที่เราต้องการใช้งาน
i ใช้ทำงานกับ rows
j ใช้ทำงานกับ columns
by ใช้จับกลุ่มข้อมูล

เราไปดูรายละเอียดการใช้งาน i, j, และ by กัน

🍦 1st Argument: i

เราสามารถใช้งาน i ได้ 3 อย่าง:

Select rows: เลือก rows ที่ต้องการ
Filter: กรองข้อมูล
Sort: จัดลำดับข้อมูล

1️⃣ Select Rows

เราสามารถใช้ i เพื่อเลือก rows ได้ 3 แบบ:

เลือก 1 row
เลือกหลาย rows
คัด rows ที่ไม่ต้องการออก

แบบที่ #1. เลือก 1 row

เช่น เลือก row ที่ 5 จาก flights:

# Select a row
flights[5]

ผลลัพธ์:

    year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
   <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
1:  2014     1     1         2         1      AA    JFK    LAX      350     2475    13

แบบที่ #2 – เลือกหลาย rows

ซึ่งเราทำได้โดยการใช้ vector เช่น เลือก rows 1 ถึง 10:

# Select a range of rows
flights[1:10]

ผลลัพธ์:

     year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
    <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
 1:  2014     1     1        14        13      AA    JFK    LAX      359     2475     9
 2:  2014     1     1        -3        13      AA    JFK    LAX      363     2475    11
 3:  2014     1     1         2         9      AA    JFK    LAX      351     2475    19
 4:  2014     1     1        -8       -26      AA    LGA    PBI      157     1035     7
 5:  2014     1     1         2         1      AA    JFK    LAX      350     2475    13
 6:  2014     1     1         4         0      AA    EWR    LAX      339     2454    18
 7:  2014     1     1        -2       -18      AA    JFK    LAX      338     2475    21
 8:  2014     1     1        -3       -14      AA    JFK    LAX      356     2475    15
 9:  2014     1     1        -1       -17      AA    JFK    MIA      161     1089    15
10:  2014     1     1        -2       -14      AA    JFK    SEA      349     2422    18

หรือเลือก rows 1, 3, 5, 7, และ 9:

# Select multiple rows at interval
flights[c(1, 3, 5, 7, 9)]

ผลลัพธ์:

    year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
   <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
1:  2014     1     1        14        13      AA    JFK    LAX      359     2475     9
2:  2014     1     1         2         9      AA    JFK    LAX      351     2475    19
3:  2014     1     1         2         1      AA    JFK    LAX      350     2475    13
4:  2014     1     1        -2       -18      AA    JFK    LAX      338     2475    21
5:  2014     1     1        -1       -17      AA    JFK    MIA      161     1089    15

แบบที่ #3 – คัด rows ที่ไม่ต้องการออก

ซึ่งเราสามารถทำได้ 2 แบบ คือ:

ใช้ -
ใช้ !

เช่น คัด row ที่ 1 ออก โดยใช้ -:

# Deselect a row
flights[-1]

หรือใช้ !:

# Deselect a row
flights[!1]

โดยทั้งสองแบบจะให้ผลลัพธ์แบบเดียวกัน แบบนี้:

         year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
        <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
     1:  2014     1     1        -3        13      AA    JFK    LAX      363     2475    11
     2:  2014     1     1         2         9      AA    JFK    LAX      351     2475    19
     3:  2014     1     1        -8       -26      AA    LGA    PBI      157     1035     7
     4:  2014     1     1         2         1      AA    JFK    LAX      350     2475    13
     5:  2014     1     1         4         0      AA    EWR    LAX      339     2454    18
    ---                                                                                    
253311:  2014    10    31         1       -30      UA    LGA    IAH      201     1416    14
253312:  2014    10    31        -5       -14      UA    EWR    IAH      189     1400     8
253313:  2014    10    31        -8        16      MQ    LGA    RDU       83      431    11
253314:  2014    10    31        -4        15      MQ    LGA    DTW       75      502    11
253315:  2014    10    31        -5         1      MQ    LGA    SDF      110      659     8

2️⃣ Filter

นอกจากการเลือก rows เรายังสามารถใช้ i เพื่อกรอง data ได้ 3 แบบ:

กรองโดยใช้ 1 เงื่อนไข
กรองโดยใช้หลายเงื่อนไข
กรองโดยใช้ helpers

แบบที่ #1 – กรองข้อมูลโดยใช้ 1 เงื่อนไข

เช่น กรอง rows ที่มีระยะทางบิน 500 miles ขึ้นไป:

# Filter with 1 condition
flights[distance >= 500]

ผลลัพธ์:

         year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
        <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
     1:  2014     1     1        14        13      AA    JFK    LAX      359     2475     9
     2:  2014     1     1        -3        13      AA    JFK    LAX      363     2475    11
     3:  2014     1     1         2         9      AA    JFK    LAX      351     2475    19
     4:  2014     1     1        -8       -26      AA    LGA    PBI      157     1035     7
     5:  2014     1     1         2         1      AA    JFK    LAX      350     2475    13
    ---                                                                                    
198323:  2014    10    31        18       -14      UA    EWR    LAS      291     2227    16
198324:  2014    10    31         1       -30      UA    LGA    IAH      201     1416    14
198325:  2014    10    31        -5       -14      UA    EWR    IAH      189     1400     8
198326:  2014    10    31        -4        15      MQ    LGA    DTW       75      502    11
198327:  2014    10    31        -5         1      MQ    LGA    SDF      110      659     8

แบบที่ #2 – กรองข้อมูลด้วยหลายเงื่อนไข

เราสามารถเพิ่มเงื่อนไขการกรองได้ด้วย logical operators:

Operator	Meaning
`&`	AND
`\|`	OR
`!`	NOT

เช่น กรอง rows ที่:

มีระยะทางบิน 500 miles ขึ้นไป และ
ออกจากสนามบิน LaGuardia (LGA):

# Filter with multiple conditions
flights[distance >= 500 & origin == "LGA"]

ผลลัพธ์:

        year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
       <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
    1:  2014     1     1        -8       -26      AA    LGA    PBI      157     1035     7
    2:  2014     1     1        -7        -6      AA    LGA    ORD      142      733     5
    3:  2014     1     1        -7         0      AA    LGA    ORD      143      733     6
    4:  2014     1     1        -8       -17      AA    LGA    ORD      139      733     6
    5:  2014     1     1        -2        15      AA    LGA    ORD      145      733     7
   ---                                                                                    
63251:  2014    10    31        14       -17      UA    LGA    IAH      200     1416    17
63252:  2014    10    31        24        -5      UA    LGA    IAH      198     1416     6
63253:  2014    10    31         1       -30      UA    LGA    IAH      201     1416    14
63254:  2014    10    31        -4        15      MQ    LGA    DTW       75      502    11
63255:  2014    10    31        -5         1      MQ    LGA    SDF      110      659     8

แบบที่ #3 – กรองโดยใช้ helpers

เราสามารถกรองข้อมูลโดยใช้ helpers หรือ operators พิเศษ 3 อย่างนี้:

Helper	For	Syntax
`%between%`	กรองข้อมูลตาม range	`col %between% range`
`%like%`	กรองข้อมูลตาม text pattern	`col %like% pattern`
`%chin%`	กรองข้อมูลอยู่ใน set ที่กำหนด	`col %chin% set`

ตัวอย่าง:

ใช้ %between% เพื่อกรองข้อมูลที่มีระยะทางบินระหว่าง 500 ถึง 1,000 miles:

# Filter using %between%
flights[distance %between% c(500, 1000)]

ผลลัพธ์:

        year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
       <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
    1:  2014     1     1        18        69      AA    JFK    ORD      155      740    17
    2:  2014     1     1        -7        -6      AA    LGA    ORD      142      733     5
    3:  2014     1     1        -7         0      AA    LGA    ORD      143      733     6
    4:  2014     1     1        -8       -17      AA    LGA    ORD      139      733     6
    5:  2014     1     1        -2        15      AA    LGA    ORD      145      733     7
   ---                                                                                    
79754:  2014    10    31        10        -5      UA    EWR    ORD      110      719     6
79755:  2014    10    31         5         2      UA    EWR    ORD      132      719    10
79756:  2014    10    31       427       393      UA    EWR    ORD      100      719    21
79757:  2014    10    31        -4        15      MQ    LGA    DTW       75      502    11
79758:  2014    10    31        -5         1      MQ    LGA    SDF      110      659     8

Note: code นี้ให้ผลลัพธ์เดียวกับการเขียน flights[distance >= 500 & distance <= 1000] แต่การใช้ %between% ทำให้ code สั้นและอ่านง่ายกว่า

ใช้ %like% เพื่อกรองข้อมูลที่สนามบินปลายทางขึ้นต้นด้วย “A” เช่น “ABQ”, “ACK”, “AGS”:

# Filter using %like%
flights[dest %like% "^A"]

ผลลัพธ์:

        year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
       <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
    1:  2014     1     1        -1         1      AA    JFK    AUS      232     1521    17
    2:  2014     1     1        -5        16      B6    JFK    AUS      247     1521    20
    3:  2014     1     1        21        21      B6    JFK    AUS      237     1521     9
    4:  2014     1     1        10         4      B6    JFK    ABQ      280     1826    20
    5:  2014     1     1        10        10      DL    LGA    ATL      126      762    18
   ---                                                                                    
15630:  2014    10    31        50        43      UA    EWR    ATL      113      746    15
15631:  2014    10    31        -5       -38      UA    EWR    ATL      111      746     5
15632:  2014    10    31        -5        -2      UA    EWR    AUS      211     1504    15
15633:  2014    10    31        -9       -15      UA    EWR    ATL      119      746    11
15634:  2014    10    31        11       -10      UA    EWR    ATL      109      746     8

ใช้ %chin% เพื่อกรองข้อมูลที่สนามบินปลายทาง คือ ATL, LAX, หรือ ORD:

# Filter using %chin%
flights[dest %chin% c("ATL", "LAX", "ORD")]

ผลลัพธ์:

        year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
       <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
    1:  2014     1     1        14        13      AA    JFK    LAX      359     2475     9
    2:  2014     1     1        -3        13      AA    JFK    LAX      363     2475    11
    3:  2014     1     1         2         9      AA    JFK    LAX      351     2475    19
    4:  2014     1     1         2         1      AA    JFK    LAX      350     2475    13
    5:  2014     1     1         4         0      AA    EWR    LAX      339     2454    18
   ---                                                                                    
38827:  2014    10    31        10        -5      UA    EWR    ORD      110      719     6
38828:  2014    10    31         3       -32      UA    EWR    LAX      320     2454    20
38829:  2014    10    31         5         2      UA    EWR    ORD      132      719    10
38830:  2014    10    31       427       393      UA    EWR    ORD      100      719    21
38831:  2014    10    31        10       -27      UA    EWR    LAX      326     2454    10

3️⃣ Sort

สุดท้าย เราสามารถใช้ i เพื่อเรียงลำดับข้อมูลได้ 3 แบบ:

Sort ascending (A—Z)
Sort descending (Z—A)
Sort by multiple columns

แบบที่ #1 – Sort ascending

เช่น จัดเรียงตามชื่อสนามบินต้นทาง จาก A—Z:

# Sort ascending
flights[order(origin)]

ผลลัพธ์:

         year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
        <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
     1:  2014     1     1         4         0      AA    EWR    LAX      339     2454    18
     2:  2014     1     1        -5       -17      AA    EWR    MIA      161     1085    16
     3:  2014     1     1       191       185      AA    EWR    DFW      214     1372    16
     4:  2014     1     1        -1        -2      AA    EWR    DFW      214     1372    14
     5:  2014     1     1        -3       -10      AA    EWR    MIA      154     1085     6
    ---                                                                                    
253312:  2014    10    31        24        -5      UA    LGA    IAH      198     1416     6
253313:  2014    10    31         1       -30      UA    LGA    IAH      201     1416    14
253314:  2014    10    31        -8        16      MQ    LGA    RDU       83      431    11
253315:  2014    10    31        -4        15      MQ    LGA    DTW       75      502    11
253316:  2014    10    31        -5         1      MQ    LGA    SDF      110      659     8

แบบที่ #2 – Sorting descending

เราเรียงข้อมูลแบบ descending (Z—A) ได้ 2 วิธี:

ใช้ decreasing = TRUE
ใช้ -

เช่น จัดเรียงตามชื่อสนามบินต้นทาง จาก Z-A โดยใช้ decreasing = TRUE:

# Sort descending with decreasing = TRUE
flights[order(origin, decreasing = TRUE)]

หรือโดยใช้ -:

# Sort descending with -
flights[order(-origin)]

เราจะได้ผลลัพธ์แบบเดียวกัน:

         year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
        <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
     1:  2014     1     1        -8       -26      AA    LGA    PBI      157     1035     7
     2:  2014     1     1        -7        -6      AA    LGA    ORD      142      733     5
     3:  2014     1     1        -7         0      AA    LGA    ORD      143      733     6
     4:  2014     1     1        -8       -17      AA    LGA    ORD      139      733     6
     5:  2014     1     1        -2        15      AA    LGA    ORD      145      733     7
    ---                                                                                    
253312:  2014    10    31        41        19      UA    EWR    SFO      344     2565    12
253313:  2014    10    31       427       393      UA    EWR    ORD      100      719    21
253314:  2014    10    31        10       -27      UA    EWR    LAX      326     2454    10
253315:  2014    10    31        18       -14      UA    EWR    LAS      291     2227    16
253316:  2014    10    31        -5       -14      UA    EWR    IAH      189     1400     8

แบบที่ #3 – Sort by multiple columns

เช่น จัดเรียงตามชื่อสนามบินต้นทางและปลายทาง ตามลำดับ:

# Sort by multiple columns
flights[order(origin, dest)]

ผลลัพธ์:

         year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
        <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
     1:  2014     1     2        -2       -25      EV    EWR    ALB       30      143     7
     2:  2014     1     3        88        79      EV    EWR    ALB       29      143    23
     3:  2014     1     4       220       211      EV    EWR    ALB       32      143    15
     4:  2014     1     4        35        19      EV    EWR    ALB       32      143     7
     5:  2014     1     5        47        42      EV    EWR    ALB       26      143     8
    ---                                                                                    
253312:  2014    10    29         0         9      MQ    LGA    XNA      174     1147     6
253313:  2014    10    29        -5       -16      MQ    LGA    XNA      162     1147    14
253314:  2014    10    30        -4       -23      MQ    LGA    XNA      154     1147     6
253315:  2014    10    30        -7       -11      MQ    LGA    XNA      157     1147    14
253316:  2014    10    31        -5       -11      MQ    LGA    XNA      165     1147     6

🧮 2nd Argument: j

เราสามารถใช้ j เพื่อทำงานได้ 3 อย่าง ได้แก่:

Select columns: เลือก columns ที่ต้องการ
Compute: วิเคราะห์ข้อมูล
Create columns: สร้าง columns ใหม่

1️⃣ Select Columns

เราสามารถใช้ j เพื่อเลือก columns ได้ 3 แบบ:

เลือก 1 column
เลือกหลาย columns
คัด columns ที่ไม่ต้องการออก

แบบที่ #1 – เลือก 1 column

เช่น เลือก column สนามบินต้นทาง:

# Select a column
flights[, "origin"]

ผลลัพธ์:

        origin
        <char>
     1:    JFK
     2:    JFK
     3:    JFK
     4:    LGA
     5:    JFK
    ---       
253312:    LGA
253313:    EWR
253314:    LGA
253315:    LGA
253316:    LGA

Note: เราสามารถใช้ตำแหน่ง (1, 2, 3, …) แทนชื่อ columns ("origin") ได้ แต่ไม่เป็นที่นิยม เพราะ columns อาจขยับตำแหน่งได้ และทำให้ผลลัพธ์เปลี่ยนไปได้

แบบที่ #2 – เลือกหลาย columns

เราใช้ j เพื่อเลือกหลาย columns ได้ 3 วิธี:

Vector
List หรือ .()
..

เช่น เลือก 3 columns คือ:

สนามบินต้นทาง
สนามบินปลายทาง
ระยะเวลาบิน

โดยใช้ vector:

# Select multiple columns with a vector
flights[, c("origin", "dest", "air_time")]

ผลลัพธ์:

        origin   dest air_time
        <char> <char>    <int>
     1:    JFK    LAX      359
     2:    JFK    LAX      363
     3:    JFK    LAX      351
     4:    LGA    PBI      157
     5:    JFK    LAX      350
    ---                       
253312:    LGA    IAH      201
253313:    EWR    IAH      189
253314:    LGA    RDU       83
253315:    LGA    DTW       75
253316:    LGA    SDF      110

โดยใช้ list:

# Select multiple columns with a list
flights[, list(origin, dest, air_time)]

หรือใช้ .() ซึ่งเป็น shorthand สำหรับ list:

# Select multiple columns using .()
flights[, .(origin, dest, air_time)]

ผลลัพธ์:

        origin   dest air_time
        <char> <char>    <int>
     1:    JFK    LAX      359
     2:    JFK    LAX      363
     3:    JFK    LAX      351
     4:    LGA    PBI      157
     5:    JFK    LAX      350
    ---                       
253312:    LGA    IAH      201
253313:    EWR    IAH      189
253314:    LGA    RDU       83
253315:    LGA    DTW       75
253316:    LGA    SDF      110

โดยใช้ .. ซึ่งเรามักใช้วิธีนี้เมื่อต้องการเลือก columns แบบ dynamic:

# Select multiple columns using ..

## Create a vector of column names
cols <- c("origin", "dest", "air_time")

## Select using ..
flights[, ..cols]

ผลลัพธ์:

        origin   dest air_time
        <char> <char>    <int>
     1:    JFK    LAX      359
     2:    JFK    LAX      363
     3:    JFK    LAX      351
     4:    LGA    PBI      157
     5:    JFK    LAX      350
    ---                       
253312:    LGA    IAH      201
253313:    EWR    IAH      189
253314:    LGA    RDU       83
253315:    LGA    DTW       75
253316:    LGA    SDF      110

Note: ความแตกต่างระหว่าง vector และ list หรือ .() และ .. คือ:

vector ให้ผลลัพธ์ที่เป็น vector
List หรือ .() และ .. ให้ผลลัพธ์เป็น data.table (data frame ของ data.table)

แบบที่ #3 – คัด columns ที่ไม่ต้องการออก

เราสามารถคัด columns ที่ไม่ต้องการออกได้ 2 วิธี:

ใช้ -
ใช้ !

เช่น เอาชื่อสายการบินออก โดยใช้ -:

# Deselect a column using -
flights[, -c("carrier")]

หรือใช้ !:

# Deselect a column using !
flights[, !c("carrier")]

ผลลัพธ์:

         year month   day dep_delay arr_delay origin   dest air_time distance  hour
        <int> <int> <int>     <int>     <int> <char> <char>    <int>    <int> <int>
     1:  2014     1     1        14        13    JFK    LAX      359     2475     9
     2:  2014     1     1        -3        13    JFK    LAX      363     2475    11
     3:  2014     1     1         2         9    JFK    LAX      351     2475    19
     4:  2014     1     1        -8       -26    LGA    PBI      157     1035     7
     5:  2014     1     1         2         1    JFK    LAX      350     2475    13
    ---                                                                            
253312:  2014    10    31         1       -30    LGA    IAH      201     1416    14
253313:  2014    10    31        -5       -14    EWR    IAH      189     1400     8
253314:  2014    10    31        -8        16    LGA    RDU       83      431    11
253315:  2014    10    31        -4        15    LGA    DTW       75      502    11
253316:  2014    10    31        -5         1    LGA    SDF      110      659     8

2️⃣ Compute

นอกจากการเลือก columns เรายังสามารถใช้ j เพื่อวิเคราะห์หรือข้อมูล (summarise, aggregate) ได้

เช่น หาค่าเฉลี่ยของระยะเวลาบิน:

# Calculate mean
flights[, mean(air_time)]

ผลลัพธ์:

[1] 156.7228

หรือ หาค่าเฉลี่ย และ standard deviation (SD) ของระยะเวลาบิน:

# Calculate mean and SD
flights[, .(avg_air_time = mean(air_time),
            sd_air_time = sd(air_time))]

ผลลัพธ์:

   avg_air_time sd_air_time      n
          <num>       <num>  <int>
1:     156.7228    96.12978 253316

Note: เราสามารถวิเคราะห์ข้อมูลโดยไม่ตั้งชื่อให้กับผลลัพธ์ได้ (avg_air_time, sd_air_time, n) แต่การตั้งชื่อจะช่วยให้เราอ่าน output ได้ง่ายขึ้น

3️⃣ Create Columns

สุดท้าย เราสามารถใช้ j สร้าง columns ใหม่ได้ร่วมกับ:

:= สำหรับสร้าง 1 column
`:=` สำหรับสร้างหลาย columns

เช่น สร้าง column ความเร็วในการบิน (ระยะทาง / เวลา) ด้วย :=:

# Creating 1 new column
flights[, speed := distance / (air_time / 60)]

ผลลัพธ์:

    year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour    speed
   <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>    <num>
1:  2014     1     1        14        13      AA    JFK    LAX      359     2475     9 413.6490
2:  2014     1     1        -3        13      AA    JFK    LAX      363     2475    11 409.0909
3:  2014     1     1         2         9      AA    JFK    LAX      351     2475    19 423.0769
4:  2014     1     1        -8       -26      AA    LGA    PBI      157     1035     7 395.5414
5:  2014     1     1         2         1      AA    JFK    LAX      350     2475    13 424.2857
6:  2014     1     1         4         0      AA    EWR    LAX      339     2454    18 434.3363

Note:

เราหาร air_time ด้วย 60 เพื่อแปลงหน่วยจากนาทีเป็นชั่วโมง และทำให้ผลลัพธ์ที่ได้เป็นหน่วย miles/hour
สังเกตว่า column ใหม่จะอยู่ท้ายสุดของ data.table

หรือ สร้าง 2 columns พร้อมกัน เช่น:

ความเร็วในการบิน (ระยะทาง / เวลา)
ระยะเวลาที่ delay โดยรวม (delay ขาออก + delay ขาเข้า)

ด้วย `:=`:

# Creating multiple new column
flights[, `:=`(speed = distance / (air_time / 60),
               total_delay = dep_delay + arr_delay)]

ผลลัพธ์:

    year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour    speed total_delay
   <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>    <num>       <int>
1:  2014     1     1        14        13      AA    JFK    LAX      359     2475     9 413.6490          27
2:  2014     1     1        -3        13      AA    JFK    LAX      363     2475    11 409.0909          10
3:  2014     1     1         2         9      AA    JFK    LAX      351     2475    19 423.0769          11
4:  2014     1     1        -8       -26      AA    LGA    PBI      157     1035     7 395.5414         -34
5:  2014     1     1         2         1      AA    JFK    LAX      350     2475    13 424.2857           3
6:  2014     1     1         4         0      AA    EWR    LAX      339     2454    18 434.3363           4

🤝 3rd Argument: by

เราใช้ by เพื่อจับกลุ่มข้อมูล ซึ่งมีประโยชน์มากเวลาที่เราต้องการวิเคราะห์ข้อมูลเป็น ๆ กลุ่ม

เราสามารถใช้ by ได้ 2 แบบ:

จับกลุ่มด้วย 1 column
จับกลุ่มด้วยหลาย columns

แบบที่ #1 – จับกลุ่มด้วย 1 column

เช่น หาค่าเฉลี่ยของ delay ขาออก ตามชื่อสนามบินต้นทาง:

# Group by 1 column
flights[, mean(dep_delay), by = origin]

ผลลัพธ์:

   origin       V1
   <char>    <num>
1:    JFK 11.44617
2:    LGA 10.60500
3:    EWR 15.21248

แบบที่ #2 – จับกลุ่มด้วยหลาย columns

เราจับกลุ่มด้วยหลาย columns ได้ด้วย 2 วิธี คือ:

Vector
List หรือ .()

เช่น หาค่าเฉลี่ยของ delay ขาออก โดยจับกลุ่มตามชื่อสนามบินต้นทางและปลายทาง ตามลำดับ

วิธีที่ 1. ใช้ vector:

# Group by with a vector
flights[, mean(dep_delay), by = c("origin", "dest")]

ผลลัพธ์:

     origin   dest        V1
     <char> <char>     <num>
  1:    JFK    LAX  8.359718
  2:    LGA    PBI 10.168617
  3:    EWR    LAX 15.882631
  4:    JFK    MIA 10.008364
  5:    JFK    SEA 10.858953
 ---                        
217:    LGA    AVL -6.500000
218:    LGA    GSP  6.000000
219:    LGA    SBN  5.000000
220:    EWR    SBN -1.500000
221:    LGA    DAL -6.266667

วิธีที่ 2. ใช้ list:

# Group by with a list
flights[, mean(dep_delay), by = list(origin, dest)]

หรือใช้ .():

# Group by with .()
flights[, mean(dep_delay), by = .(origin, dest)]

ผลลัพธ์:

     origin   dest        V1
     <char> <char>     <num>
  1:    JFK    LAX  8.359718
  2:    LGA    PBI 10.168617
  3:    EWR    LAX 15.882631
  4:    JFK    MIA 10.008364
  5:    JFK    SEA 10.858953
 ---                        
217:    LGA    AVL -6.500000
218:    LGA    GSP  6.000000
219:    LGA    SBN  5.000000
220:    EWR    SBN -1.500000
221:    LGA    DAL -6.266667

Note: เช่นเดียวกับการเลือก columns …

ถ้าเราใช้ vector เราจะได้ผลลัพธ์เป็น vector
ถ้าใช้ list หรือ .() เราจะได้ data.table

🔗 Combining & Chaining

เราสามารถปลอดล็อกพลังที่แท้จริงของ data.table ได้ด้วย 2 วิธี:

Combining: ใช้ i, j, by ร่วมกัน เพื่อตอบโจทย์ที่ซับซ้อน
Chaining: เชื่อมต่อ data.table เข้าด้วยกัน

1️⃣ Combining

ยกตัวอย่างเช่น:

หาค่าเฉลี่ยของความเร็ว
เฉพาะเที่ยวบินที่มีระยะทางตั้งแต่ 500 miles ขึ้นไป
โดยจับกลุ่มตามชื่อสนามบินต้นทาง:

# Combining
flights[distance >= 500, 
        .(avg_speed = mean(distance / (air_time / 60))), 
        by = origin]

อธิบาย code:

distance >= 500 เลือกเฉพาะข้อมูลที่มีระยะทางตั้งแต่ 500 miles ขึ้นไป
.(avg_speed = mean(distance / (air_time / 60))) หาค่าเฉลี่ยของความเร็ว
by = origin] จับกลุ่มข้อมูลด้วยสนามบินต้นทาง

ผลลัพธ์:

   origin avg_speed
   <char>     <num>
1:    JFK  436.6284
2:    LGA  407.5462
3:    EWR  417.1765

2️⃣ Chaining

Chaining คือ การนำ code มาต่อกัน เพื่อส่ง output ต่อกันเป็นทอด ๆ:

DT[...][...][...]

เช่น ค้นหาสนามบินปลายทาง 5 อันดับแรกที่มีค่าเฉลี่ย delay ขาเข้ามากที่สุด ในเดือนสิงหาคม:

# Chaining
flights[month == 8,
        .(avg_arr_delay = mean(arr_delay)),
        by = dest][order(-avg_arr_delay)][1:5]

อธิบาย code:

month == 8 เลือกข้อมูลจากเดือนสิงหาคม
.(avg_arr_delay = mean(arr_delay)) คำนวณค่าเฉลี่ย delay ขาเข้า
by = dest จับกลุ่มตามสนามบินปลายทาง
[order(-avg_arr_delay)] จัดลำดับตามค่าเฉลี่ย delay ขาเข้า แบบ descending
[1:5] เลือกเฉพาะ 5 rows แรกมาแสดง

ผลลัพธ์:

     dest avg_arr_delay
   <char>         <num>
1:    LIT      37.06452
2:    DSM      22.85714
3:    CAK      19.60976
4:    TYS      19.44681
5:    TVC      19.00000

🍩 Special Symbols

เพื่อช่วยให้เราทำงานได้ง่ายขึ้น data.table มี special symbols 3 ตัวที่เราสามารถเรียกใช้ในการทำงานได้:

.N
.SD
.SDcols

1️⃣ .N

.N เป็น special symbol ที่เก็บจำนวน rows ของ data.table เอาไว้ (เช่น 500)

เรามักใช้งาน .N ใน 2 กรณีหลัก ๆ ได้แก่:

เลือก rows
นับจำนวนข้อมูล

กรณีที่ 1. เลือก rows

เช่น เลือก row ที่ 500 จนถึง row สุดท้าย:

# Select rows with .N
flights[500:.N]

ผลลัพธ์:

         year month   day dep_delay arr_delay carrier origin   dest air_time distance  hour
        <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
     1:  2014     1     1        81        86      WN    EWR    HOU      222     1411    17
     2:  2014     1     1        -3        -4      WN    EWR    MDW      128      711     7
     3:  2014     1     1         0        22      WN    EWR    MDW      144      711    12
     4:  2014     1     1        88       190      WN    EWR    MDW      130      711    21
     5:  2014     1     1        45        63      WN    EWR    MDW      141      711    16
    ---                                                                                    
252813:  2014    10    31         1       -30      UA    LGA    IAH      201     1416    14
252814:  2014    10    31        -5       -14      UA    EWR    IAH      189     1400     8
252815:  2014    10    31        -8        16      MQ    LGA    RDU       83      431    11
252816:  2014    10    31        -4        15      MQ    LGA    DTW       75      502    11
252817:  2014    10    31        -5         1      MQ    LGA    SDF      110      659     8

กรณีที่ 2. นับจำนวนข้อมูล

เช่น นับจำนวนข้อมูลการบินตามสนามบินต้นทางแต่ละแห่ง:

# Compute with .N
flights[, .N, by = origin]

ผลลัพธ์:

   origin     N
   <char> <int>
1:    JFK 81483
2:    LGA 84433
3:    EWR 87400

2️⃣ .SD

.SD ย่อมาจาก “Subset of Data” ซึ่งหมายถึง ชุดข้อมูลย่อยที่เกิดจากการจับกลุ่มด้วย by

เรามักใช้ .SD ในการคำนวณคู่กับ lapply() function เพื่อวิเคราะห์ข้อมูลตามกลุ่ม

เช่น หาค่าสูงสุดของแต่ละ columns ในแต่ละเดือน:

# Compute with .SD
flights[,
        lapply(.SD, max, na.rm = TRUE),
        by = month]

ผลลัพธ์:

    month  year   day dep_delay arr_delay carrier origin   dest air_time distance  hour
    <int> <int> <int>     <int>     <int>  <char> <char> <char>    <int>    <int> <int>
 1:     1  2014    31       973       996      WN    LGA    XNA      688     4983    24
 2:     2  2014    28      1014      1007      WN    LGA    XNA      685     4983    24
 3:     3  2014    31       920       925      WN    LGA    XNA      706     4983    24
 4:     4  2014    30      1241      1223      WN    LGA    XNA      664     4983    24
 5:     5  2014    31       889       879      WN    LGA    XNA      650     4983    24
 6:     6  2014    30      1071      1073      WN    LGA    XNA      640     4983    24
 7:     7  2014    31      1087      1090      WN    LGA    XNA      638     4983    24
 8:     8  2014    31       978       964      WN    LGA    XNA      635     4983    24
 9:     9  2014    30      1056      1115      WN    LGA    XNA      635     4983    23
10:    10  2014    31      1498      1494      WN    LGA    XNA      662     4983    24

3️⃣ .SDcols

.SDcols เป็น special symbol ที่เก็บชื่อ columns ของ .SD เอาไว้

เรามักใช้ .SDcols คู่กับ .SD และ lapply() เพื่อวิเคราะห์เฉพาะ columns ที่ต้องการ

เช่น หาเวลา delay สูงสุดของขาเข้าและขาออกในแต่ละเดือน:

# Compute with .SDcols
flights[,
        lapply(.SD, max, na.rm = TRUE),
        by = month,
        .SDcols = c("arr_delay", "dep_delay")]

ผลลัพธ์:

    month arr_delay dep_delay
    <int>     <int>     <int>
 1:     1       996       973
 2:     2      1007      1014
 3:     3       925       920
 4:     4      1223      1241
 5:     5       879       889
 6:     6      1073      1071
 7:     7      1090      1087
 8:     8       964       978
 9:     9      1115      1056
10:    10      1494      1498

😎 Conclusion

ในบทความนี้ เราได้ไปดูวิธีใช้ data.table เพื่อทำทำงานกับข้อมูลขนาดใหญ่กัน

เราได้เห็นว่า data.table มีหลักการเขียนดังนี้:

DT[i, j, by]

และแต่ละ argument มีการใช้งานดังนี้:

Argument i:

Use Case	Example
Select rows	`flights[1:5]`
Filter	`flights[distance >= 500]`
Sort	`flights[order(origin)]`

Argument j:

Use Case	Example
Select columns	`flights[, .(origin, dest)]`
Compute	`flights[, mean(air_time)]`
Create columns	`flights[, speed := distance / (air_time / 60)]`

Argument by:

Use Case	Example
Group data	`flights[by = origin]`

Special symbols:

Symbol	Meaning
`.N`	จำนวน rows
`.SD`	Subset of Data
`.SDcols`	columns ใน Subset of Data

เราสามารถใช้ทุก arguments ร่วมกัน (combining) หรือเชื่อมต่อ data.table (chaining) เพื่อตอบโจทย์ที่ซับซ้อนได้

🎒 Learn More About data.table

😺 GitHub

ดู code ตัวอย่างทั้งหมดในบทความนี้ได้ที่ GitHub

📚 อ่านเพิ่มเติมเกี่ยวกับ data.table

📑 Cheat Sheets ในการใช้งาน data.table

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

📕 ขอฝากหนังสือเล่มแรกในชีวิตด้วยนะครับ 😆

Correlation
t-tests
ANOVA
Reliability
Factor analysis

แล้วทุกคนจะแปลกใจว่า ทำไมภาษา R ง่ายขนาดนี้ 🙂‍↕️

👉 สนใจดูรายละเอียดหนังสือได้ที่ meb:

ดูรายละเอียดหนังสือ R Book for Psychologists

2025-02-13

Category: R

🌲 Random Forest Model คืออะไร?

💻 Random Forest Models ในภาษา R

🚗 mpg Dataset

🐣 ranger Basics

1️⃣ ติดตั้งและโหลด ranger

2️⃣ สร้าง Training และ Test Sets

3️⃣ สร้าง Random Forest Model

4️⃣ ทดสอบความสามารถของ Model

⏲️ Hyperparametre Tuning

🍩 Bonus: Variable Importance

😎 Summary

😺 GitHub

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

Share this:

💎 Example Dataset: diamonds

⬇️ Load diamonds

🍳 Prepare the Dataset

🪆 Step 1. One-Hot Encoding

📈 Step 2. Log Transformation

🚄 Step 3. Split the Data

🏷️ Linear Regression Modelling

💪 Step 1. Fit the Model

🔮 Step 2. Make Predictions

🎯 Step 3. Evaluate the Model Performance

😎 Summary

😺 GitHub

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

Share this:

🧐 PCA คืออะไร?

💻 PCA ในภาษา R

🧑‍💻 prcomp()

🍷 wine dataset

⚒️ PCA

📈 ดูผลลัพธ์

😺 GitHub

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

Share this:

🧐 HCA คืออะไร?

💻 HCA ในภาษา R: hclust()

🔢 Example Dataset: USArrests

📏 Normalise the Data

🍿 HCA With hclust()

📈 Print & Dendrogram

😺 GitHub

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

Share this:

👉 Introduction to k-Means

🤔 k-Means คืออะไร?

🪜 Steps การทำงานของ k-Means

🏫 Learn More

💻 k-Means ในภาษา R: kmeans()

🔢 Dataset

🪨 rock

📏 Data Normalisation

🔎 Finding the Optimal k

⏹️ k-Means

🔥 Train the Model

🤓 Get the Results

😺 GitHub

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

Share this:

🌳 Tree-Based Algorithms คืออะไร?

💻 Tree-Based Models ในภาษา R

🚗 Dataset: mtcars

🔧 Prepare the Data

🪴 Train the Model

📏 Evaluate the Model

🍄 Random Forest

🐱 GitHub

📃 References

✅ R Book for Psychologists: หนังสือภาษา R สำหรับนักจิตวิทยา

Share this:

📈 Logistic Regression คืออะไร?

🔨 วิธีสร้าง Logistic Regression ในภาษา R