สวัสดีคุณผู้อ่านทุกท่านค้าบ 😆 ฮี่ๆ วันนี้พามาใช้ R ในการทำโมเดล Logistic regression กันฮะ โดยใช้ข้อมูลที่ Data Analyst น่าจะเคยผ่านหูผ่านตากันมาบ้างสักครั้งนึง เพราะ “ถ้าคุณโดด ผมโดด” 😭 และสุดท้ายหนุ่มต้องการจะสร้าง App เพื่อ Role play เป็นคนในเรือดูว่าจะรอด ไม่จม ไม่แข็งตายไปรึป่าว

กดไล่ดูตามหัวข้อได้เลยคับ

🤝 รู้จักข้อมูลกันก่อน
🥘 Prep ข้อมูลให้พร้อมใช้ในโมเดล
🚂 Train Model: พี่จะสอนน้องเอง
🧪 Test Model: ลองดูก่อน
🍌มา Role play แบบเรียลไทม์กัน

🤝 รู้จักข้อมูลกันก่อน

Data set นี้มาจาก library(titanic) ของ R มีเรคคอร์ดทั้งหมด 891 rows เป็น 891 คนที่ได้ขึ้นเรือไทนานิคและเก็บข้อมูลมาได้ (มั้งง?)

Columns	Description
PassengerId	หมายเลขประจำตัวผู้โดยสาร (Unique ID)
Survived	สถานะการรอดชีวิต (0 = ไม่รอด, 1 = รอด)
Pclass	ระดับชั้นโดยสาร (1 = ชั้น 1, 2 = ชั้น 2, 3 = ชั้น 3)
Name	ชื่อผู้โดยสาร
Sex	เพศ (male, female)
Age	อายุ
SibSp	จำนวนพี่น้อง/คู่สมรสที่เดินทางมาด้วย
Parch	จำนวนพ่อแม่/ลูกที่เดินทางมาด้วย
Ticket	หมายเลขตั๋วโดยสาร
Fare	ค่าโดยสาร
Cabin	หมายเลขห้องโดยสาร
Embarked	ท่าเรือที่ขึ้นเรือ (C = Cherbourg, Q = Queenstown, S = Southampton)

มาจับโค้ดดูข้อมูลกันซักนิดว่าเบื้องต้น เช่น ขนาดข้อมูลมีเท่าไหร่, ตัวแปรที่มีเป็นประเภทไหนเป็น Qualitative เชิงคุณภาพ หรือ Quantitative เชิงปริมาณ, มีช่องว่างหรือค่าที่หายไป (missing values) ในคอลัมน์มั้ย

# Import Library and Dataset Titanic
library(tidyverse)
library(titanic)

summary(titanic_train)
count(titanic_train)

ในตัวชุดข้อมูลมี หมายเลขประจำตัวผู้โดยสาร (PassengerId), สถานะการรอดชีวิต (Survived), ระดับชั้นโดยสาร (Pclass), ชื่อผู้โดยสาร (Name), เพศ (Sex), หมายเลขตั๋วโดยสาร (Ticket), หมายเลขห้องโดยสาร (Cabin) และ ท่าเรือที่ขึ้นเรือ (Embarked) ที่เป็นตัวแปรเชิงคุณภาพ (Qualitative variables) นอกนั้นเป็นตัวแปรเชิงปริมาณ (Quantitative Variables)

ตรวจสอบค่าว่าง กับค่า na กันซักหน่อย เพราะมันส่งผลต่อโมเดลเหมือนกัน

# ตรวจสอบค่า NA โดยใช้ sapply()
sapply(titanic_train, function(x) sum(is.na(x)))

# ตรวจสอบค่าว่าง "" โดยใช้ sapply
sapply(titanic_train, function(x) sum(x == ""))

จะเห็นว่าแค่ตรวจ na อย่างเดียวไม่ได้ต้องตรวจค่าว่าง “” ด้วย มีตัวแปร Age, Cabin และ Embarked ที่มีค่าว่างอยู่ และในมุมมองของหนุ่ม หมายเลขประจำตัวผู้โดยสาร (PassengerId), ชื่อผู้โดยสาร (Name), หมายเลขตั๋วโดยสาร (Ticket), หมายเลขห้องโดยสาร (Cabin) และ ท่าเรือที่ขึ้นเรือ (Embarked) ไม่น่าจะส่งผลต่อการรอดชีวิตของการรอดชีวิตในครั้งนี้ ก็เลยจะ ✂️ ตัดตัวแปรนี้ออกไปนะครับ

ตัวแปรที่ต้องจัดการต่อก็จะเหลือเพียงอายุ (Age) เท่านั้น จึงพิจารณาตัดออกเช่นเดียวกันเพราะ Sample size ของเรา (หลังจากตัด age ว่างออกแล้ว) ก็ค่อนข้างจะใหญ่อยู่ โดยใช้ na.omit() ~เหลืออยู่ 714 rows

# copy original dataset to my object
titanic_df <- titanic_train

# Feature Selection:: Drop Not use Columns
titanic_df$PassengerId <- NULL
titanic_df$Name <- NULL
titanic_df$Ticket <- NULL
titanic_df$Cabin <- NULL
titanic_df$Embarked <- NULL

# Feature Selection:: Drop Age.na
titanic_df <- na.omit(titanic_df)

summary(titanic_df)
count(titanic_df) #714 rows

🥘 Prep ข้อมูลให้พร้อมใช้ในโมเดล

ขอแวะจัดข้อมูลสักแป๊ปนะครับ อยากให้หัวคอลัมน์เป็นพิมพ์เล็กทั้งหมด และก็แปลงตัวแปรเชิงคุณภาพที่เหลืออยู่อย่างเพศที่ยังเป็น male, female อยู่ให้เป็นตัวแปรเชิงปริมาณ (numerical variables) โดยใช้ dummy variables และเพิ่มตัวแปร surv แปลงตัวเลขเป็นข้อความ ใช้ตอน explorer ดูแนวโน้มข้อมูล

# Data Transformation
titanic_df <- titanic_df %>%
  select(pclass = Pclass, sex = Sex, age = Age, sibsp = SibSp, parch = Parch, fare = Fare, survived = Survived) %>%
  mutate( sex = if_else(sex=='male', 0, 1),
				  surv = if_else(survived==0, 'Died', 'Survived'))
  
head(titanic_df)

1. ส่อง Age: ค่า median ของอายุอยู่ที่เกือบถึง 30 ปี (~ประมาณ 28-30 ปี) ไม่ว่าจะ Died หรือ Survived โดยอายุมากสุดอยู่ที่ประมาณ 80 ปี+

# boxplot:: Age
titanic_df %>% 
  ggplot(aes(age)) +
  geom_boxplot()

# boxplot:: Age+surv
titanic_df %>% 
  ggplot(aes(x=age, y=surv)) +
  geom_boxplot()

2. ส่อง Fare: ในส่วนของค่าโดยสารมี median อยู่ 15USD มากสุดจ่ายอยู่ที่ 500USD++ หากดูจากกราฟดูเหมือนว่าคนที่จ่ายค่าโดยสารสูงจะมีโอกาสรอดมากกว่านะเนี่ย 🤔

# boxplot:: Fare
titanic_df %>% 
  ggplot(aes(fare)) +
  geom_boxplot()

# boxplot:: Fare+surv
titanic_df %>% 
  ggplot(aes(x=fare, y=surv)) +
  geom_boxplot()

3. ส่อง Pclass: คร่าวๆ แล้วในระดับชั้นโดยสาร ผู้โดยสารชั้น 3 มีจะการเสียชีวิตมากว่าชั้นอื่นเลยครับ ถ้าเรากลับไปดูโครงสร้างเรือไททานิคแล้ว ชั้น 3 จัดอยู่ส่วนลึกสุดของเรืออาจทำให้หนีเอาชีวิตรอดไม่ทัน 😭

# bar:: pclass
titanic_df %>% 
  ggplot(mapping = aes(x=pclass)) +
  geom_bar(fill = 'salmon', alpha=0.6)

# bar:: pclass+surv
titanic_df %>% 
  ggplot(mapping = aes(x=pclass)) +
  geom_bar(fill = 'salmon', alpha=0.6) +
  facet_wrap(~surv, ncol=2)

4. ส่อง Sex: จากที่เราได้ mutate sex ให้ male เป็น 0 และ female เป็น 1 แล้ว จะเห็นว่า male มีจำนวนผู้เสียชีวิตมากกว่า female ประมาณ 4 เท่า

# bar:: sex
titanic_df %>% 
  ggplot(mapping = aes(x=sex)) +
  geom_bar(fill = 'salmon', alpha=0.6)

# bar:: sex+surv
titanic_df %>% 
  ggplot(mapping = aes(x=sex)) +
  geom_bar(fill = 'salmon', alpha=0.6) +
  facet_wrap(~surv, ncol=2)

5. ส่อง Sibsp และ Parch: จำนวนพี่น้อง/คู่สมรสที่เดินทางมาด้วย และ จำนวนพ่อแม่/ลูกที่เดินทางมาด้วย มีความถี่ที่มากสุดคือ มาคนเดียวนะครับ ที่เป็นแบบนี้อาจเป็นเพราะข้อมูลนี้นับรวม cabin crew ด้วยที่ต้องขึ้นมาทำงานคนเดียว

# bar:: sibsp
titanic_df %>% 
  ggplot(mapping = aes(x=sibsp)) +
  geom_bar(fill = 'salmon', alpha=0.6)

# bar:: sibsp+surv
titanic_df %>% 
  ggplot(mapping = aes(x=sibsp)) +
  geom_bar(fill = 'salmon', alpha=0.6) +
  facet_wrap(~surv, ncol=2)

# bar:: parch
titanic_df %>% 
  ggplot(mapping = aes(x=parch)) +
  geom_bar(fill = 'salmon', alpha=0.6)

# bar:: parch+surv
titanic_df %>% 
  ggplot(mapping = aes(x=parch)) +
  geom_bar(fill = 'salmon', alpha=0.6) +
  facet_wrap(~surv, ncol=2)

6. ส่องตัวสุดท้ายเป็นเรื่อง สถานะการรอดชีวิต (0 = ไม่รอด, 1 = รอด) จะเห็นว่าอัตราการเสียชีวิตสูงกว่าบุคคลที่รอดชีวิต

# bar:: survived
titanic_df %>% 
  ggplot(mapping = aes(x=survived)) +
  geom_bar(fill = 'salmon', alpha=0.6)

เมื่อลองดูแต่ละตัวแปรแนวโน้มเบนไปทางไม่รอดชีวิตจะเหตุการณ์นี้ซะเยอะเลย แล้วตัวแปรไหนกันนะที่ส่งผลต่อการรอดชีวิตจริงๆ ก็เลยทดสอบทุกตัวแปรเข้าไปในโมเดล เพื่อดูว่าตัวแปรใดที่ Significant บ้าง (อย่าลืม Drop ตัวแปล surv ทิ้งไปก่อนนะครับ เดี๋ยวเพี้ยน) จะพบว่า ….

# Drop Label Variable titanic_df$surv
titanic_df$surv <- NULL

# Test Multivalidate
model_all <- glm(survived ~ ., data = titanic_df, family = 'binomial')
summary(model_all)

• มีตัวแปร 2 ตัวที่มีค่า p-value มากกว่า 0.05 เราจึงจะพิจารณา Drop ทิ้งไป คือ parch, fare
• ตัวแปรที่เราจะเอาไปใช้ในโมเดล คือ pclass, sex, age, sibsp

# Feature Selection:: Drop Not use Cols
titanic$parch <- NULL
titanic$fare <- NULL

🚂 Train Model: พี่จะสอนน้องเอง

เข้าสู่ขั้นตอนการเทรนแล้ว โดยเราต้องแบ่ง Train Dataset กับ Test Dataset โดยใช้โค้ดประมาณนี้นะครับ

# Split data 80:20
set.seed(137)
n <- nrow(titanic)
id <- sample(1:n, size = n*0.8)
train_df <- titanic[id, ]
test_df <- titanic[-id, ]

ซึ่งเราใช้ฟังก์ชัน sample ในการแบ่ง dataset ออกเป็น 80% เป็น Train dataset และอีก 20% เป็น Test dataset

# train model
logit_model <- glm(survived ~ pclass + sex + age + sibsp, data = train_df, 
										family = "binomial")
p_train <- predict(logit_model, type = "response") # probability %

train_df$p_train <- p_train
train_df$pred <- if_else(p_train >= 0.5, 1, 0)
train_df$pred_name <- if_else(p_train >= 0.5, "Survived", "Died")
tail(train_df, 10)

• ในสมการของเราจะได้ survived เป็นฟังก์ชันของ pclass + sex + age + sibsp
• family = ‘binomial’ ใช้สำหรับ logistic regression ที่ต้องการผลลัพธ์เป็น 2 ทาง (Survived, Died)
• ใน predict เราระบุ type = “response” เพื่อต้องการผลลัพธ์เป็น % Probability และเก็บผลลัพธ์ที่ได้เข้าใน titanic_df$p_train
• เพื่อให้อ่านผล % Prob. ง่ายขึ้นจึงแปลงให้เป็น Label โดยเงื่อนไข if_else(p_train >= 0.5, 1, 0) และ if_else(p_train >= 0.5, “Survived”, “Died”)
• เมื่อลองแสดงผลลัพธ์ออกมา จะได้ … (ทายถูกเยอะอยู่นะเนี่ยย 😝)

มาลองทดสอบ accuracy ของโมเดลกัน

## confusion matrix table() cross
conf <- table(train_df$pred, train_df$survived, 
              dnn=c('predicted', 'actual'))
conf

## evaluation
accu <- (conf[1,1] + conf[2,2]) / sum(conf)
precision <- conf[2,2] /  (conf[2,1] + conf[2,2])
recall <- conf[2,2] / (conf[1,2] + conf[2,2])
f1score <- 2 *((precision * recall) / (precision + recall))

cat('Accuracy : ', accu,
    '\nPrecision : ', precision,
    '\nRecall : ', recall,
    '\nF1 Score : ', f1score)

• เราเอามาสร้างเป็นตาราง Confusion Matrix และคำนวณ ..
• Accuracy: ไม่ว่าผลลัพธ์จะเป็น Survived หรือ Died แล้วโมเดลทายถูก / จำนวนทั้งหมด
• Precision: เมื่อเรา Predicted ว่า Survived → เราทายถูกกี่ครั้ง / จำนวนที่เราทายว่า Survived ทั้งหมด
• Recall: เมื่อข้อมูล Actual เป็น Survived → เราทายถูกกี่ครั้ง / จำนวนที่ข้อมูล Actual เป็น Survived ทั้งหมด
• F1 Score: วัดค่าผลลัพธ์โดยรวมของโมเดล
• พี่ว่าพี่โอเคกับผลลัพธ์แล้วนะ ทีนี้ไปลองกับ Test Dataset กัน

🧪 Test Model: ลองดูก่อน

เราจะใช้ test_df ที่ได้แบ่งไว้ก่อนหน้านี้ ในการทดสอบโมเดล

# test model
p_test <- predict(logit_model, newdata = test_df, type = "response") # probability %
test_df$p_test <- p_test
test_df$pred <- if_else(p_test >= 0.5, 1, 0)
test_df$pred_name <- if_else(p_test >= 0.5, "Survived", "Died")
tail(test_df, 10)

## confusion matrix table() cross
conf <- table(test_df$pred, test_df$survived, 
              dnn=c('predicted', 'actual'))
conf

## evaluation
accu <- (conf[1,1] + conf[2,2]) / sum(conf)
precision <- conf[2,2] /  (conf[2,1] + conf[2,2])
recall <- conf[2,2] / (conf[1,2] + conf[2,2])
f1score <- 2 *((precision * recall) / (precision + recall))

cat('Accuracy : ', accu,
    '\nPrecision : ', precision,
    '\nRecall : ', recall,
    '\nF1 Score : ', f1score)

• จะเห็นว่าค่า Accuracy เพิ่มขึ้นมาจาก Train dataset อาจเป็นเพราะขนาดของ Test dataset มีขนาดเล็กมาก ค่า Accuracy อาจมีความผันผวน หากจำนวนของ Test dataset เพิ่มมากขึ้นอาจจะทำให้โมเดลทายผิดเพิ่มขึ้นก็ได้ ไม่ก็ dataset ที่ใช้ทดสอบนั้นง่ายเกินไป ฮ่าๆๆ
• แต่เรามั่นใจได้ว่ามีการแบ่ง Train & Test โดยวิธีการ Random sampling แล้ว

🍌มา Role play แบบเรียลไทม์กัน

จากค่า coefficient ที่ได้จาก Train model จะสามารถเขียนสมการได้ ดังนี้

• logit(p) = 3.07026 – 1.29951 * pclass + 2.44141 * sex – 0.04553 * age – 0.37714 * sibsp
• ลองกำหนดค่าตัวแปรให้โมเดลลองทำนายดู เช่น ผู้โดยสารชั้น 2 (pclass_ = 2) เป็นผู้ชาย (sex_ = 0) อายุประมาณ 47 ปี (age_ = 2) และมีมากับคู่สมรส 1 คน (sibsp_ = 1)

coefs <- coef(logit_model)

pclass_ <- 2
sex_ <- 0
age_ <- 47
sibsp_ <- 1

# Calc Realtime
survived_pred <- coefs[[1]] + coefs[[2]]*pclass_ + coefs[[3]]*sex_ + coefs[[4]]*age_ + coefs[[4]]*sibsp_ 
survived_pred

# Calc Realtime -> Prop
survived_prop <- exp(survived_pred) / (1 + exp(survived_pred))
survived_prop

# Calc Realtime -> Prop -> Survived?
is_survived <- if_else(survived_prop >= 0.5, 1, 0)
is_survived

# Calc Realtime -> Prop -> SurvivedName?
is_survived_name <- if_else(survived_prop >= 0.5, "Survived", "Died")
is_survived_name

ผลลัพธ์คือ ชายท่านนี้ไม่รอดชีวิตจากเหตุการณ์ Titanic พอได้ลองเล่นข้อมูลแล้ว หนุ่มก็คิดว่าน่าจะลอง build application ขึ้นมาลองใช้งานดูโดยใช้ shiny app จะได้หน้าต่าง app เป็นแบบนี้เลย คุณผู้อ่านสามารถคลิกไปลองเล่นกันได้นะครับ ดีไม่ดียังไงคอมเม้นต์บอกกันได้นะครับ

ไปเล่นได้ที่ลิงก์นี้ค้าบบ: Data Modeling DSB Batch#11

---

จบแล้วครับ สำหรับการสร้างโมเดล Logistic Regression แบบง่ายๆ เพิ่มเสริมความเข้าใจหลังจากที่เข้าเรียนมาอย่างหนัก 😭 ครั้งหน้าจะเป็นโมเดลไหน ทำยังไง สามารถกดติดตามข่าวสารรอได้เลยนะครับ 🤩