หลังจากที่ได้ลองใช้ Python ทำรายงานมาแล้ว คราวนี้ก็มาลองใช้ R ในการทำรายงานดูบ้างครับ เพราะหนุ่มเชื่อว่าแค่การเรียนอย่างเดียว มันจำไม่ได้หรอก ต้องเรียนจากประสบการณ์ และการใช้งานกับสถานการณ์จริงด้วย 😎

ครั้งนี้มาพร้อมกับโจทย์ที่ว่า “จำนวนผู้ป่วยทั้งหมดที่รับไว้ในหอผู้ป่วยในและอัตราการครองเตียงเดือนนี้เป็นเท่าใด? เพิ่มขึ้นหรือลดลงแค่ไหนเมื่อเทียบกับเดือนที่แล้ว” สร้างกราฟ และตารางตัวเลขให้พี่ด้วย

หลังจากที่ไปดีลกับทีม DE รบกวนดึงข้อมูลให้หน่อยค้าบบ ได้ไฟล์มาแล้วก็มาเริ่มกันเลย ✌️

1. Import data
2. Data exploratory analysis
3. Data preparation
4. Data report & analysis

---

Import data

เกริ่นไว้ก่อนนะครับ ข้อมูลที่ใช้แสดงตัดมาใช้แค่สองเดือนนะครับเป็นข้อมูลผู้ป่วยใน (In-patient Department: IPD) ของโรงพยาบาลและเข้ารหัส (Encode) ไว้แล้วน้า .. โอเค เริ่ม

ใช้ Library readxl ในการอ่านไฟล์ Excel ที่ได้มาจากทีม DE ใช้ทั้งหมด 3 ชีท

library(readxl)
ipd <- read_excel('RefFiles/ipd-202410-202411.xlsx', sheet = 1)
dtype <- read_excel('RefFiles/ipd-202410-202411.xlsx', sheet = 2)
beds <- read_excel('RefFiles/ipd-202410-202411.xlsx', sheet = 3)

• ipd จะเป็นข้อมูลการรับเข้าเป็นผู้ป่วยใน (admit) เป็นราย individual
• dtype เป็น master table ที่เราจะใช้ join หาความหมายของคอลัมน์ ‘dchtype’
• beds เป็น master table ที่ใช้ join หาจำนวนเตียงของแต่ละ ward เพื่อคำนวนหาอัตราการครองเตียง

Data exploratory analysis

มาสำรวจดูข้อมูลกันก่อน 🔍 แต่ละชีทมีข้อมูลเป็นยังไงมั้ง

a. Raw Data

ipd มีประมาณ 6,500 แถว คร่าวๆ ก็ตกประมาณเดือนละ 3,200 รายที่รับเข้าเป็นผู้ป่วยใน มีรายละเอียดคอลัมน์คือ

• id หมายเลขรหัสผู้ป่วย
• admit_date วันที่รับเป็นผู้ป่วยใน (admit)
• dis_date วันที่จำหน่ายผู้ป่วยใน (discharge)
• ward เลขรหัสหอผู้ป่วย (จริงๆ ควรเป็นชื่อหอผู้ป่วยแหละ ขอปิดไว้ละกันนะครับ)
• dchtype ประเภทการจำหน่ายผู้ป่วย
• cost ค่าใช้จ่ายที่เกิดขึ้นในการรักษาพยาบาล

dtype มี 2 คอลัมน์ที่เป็นเลขรหัสและความหมาย

beds มี 2 คอลัมน์เหมือนกัน โดยมี ward เป็น primary key ใช้มองหาจำนวนเตียงในแต่ละหอผู้ป่วย

b. Data structure

เรามาดูโครงสร้างกันเล็กน้อยนะครับ

# Exploratory data analysis ####
# Data structure
str(ipd)
summary(ipd)

จากรูปมีคอลัมน์:

• id เป็น character
• admit_date กับ dis_date เป็น datetime
• ward กับ dchtype เป็น numeric ซึ่งเดี๋ยวเราจะเปลี่ยน type ให้เป็น character เพราะเราไม่ได้ใช้ในการคำนวณ
• cost เป็น numeric

จะเห็นว่า .. คอลัมน์ dchtype มีค่าว่างอยู่ 6 แถว และคอลัมน์ cost ค่าใช้จ่ายที่เป็น Max นี่ค่อนข้างสูงเลยนะเนี่ย (ประมาณ 1 ล้านบาท 😲)

c. Missing data

# Missing value
## missing ของแต่ละ cols
apply(ipd, 2, function(x) any(is.na(x)))

## missing ของแต่ละ col มีจำนวนเท่าไหร่
colSums(is.na(ipd))

เราสามารถตรวจสอบ missing รายคอลัมน์อีกครั้งนึง โดยใช้ is.na และผลลัพธ์ก็ยังคงเป็น dchtype มี 6 แถวที่เป็นค่าว่าง

Data preparation

a. Transform data type

# Transform data type
ipd$ward <- as.character(ipd$ward)
ipd$dchtype <- as.character(ipd$dchtype)

dtype$dchtype <- as.character(dtype$dchtype)
beds$ward <- as.character(beds$ward)

ward และ dchtype เป็นเพียง key ที่เราไม่ได้ใช้คำนวณ จึงเปลี่ยนให้เป็น character โดยใช้ as.character

b. Manage missing data

## dchtype ที่หายไปต้องการแทนที่ด้วย 5
ipd$dchtype[is.na(ipd$dchtype)] <- 5

## ตรวจอสบ missing data ของแต่ละ col อีกครั้ง
colSums(is.na(ipd))

assume ว่า dchtype ที่เป็นค่าว่างจะต้องถูกแทนที่ด้วย 5 นะครับและทำการตรวจสอบค่าว่างอีกครั้ง Make sure ว่าถูกแทนที่แล้วจริงๆ

c. Transform data

เรามาคำนวน LoS: Length of Stay (จำนวนวันนอน) กันก่อนนะครับ โดยใช้ difftime หาระยะห่างระหว่าง admit_date และ dis_date และหากผลลัพธ์ที่ได้เป็น 0 วัน ก็ให้ปรับเป็น 1 วัน (เพราะถือว่าผู้ป่วยได้ใช้ทรัพยากรของ รพ. แล้ว)

## Calc LoS: Length of Stay
ipd$los <- as.integer(difftime(ipd$dis_date, ipd$admit_date, units = "days"))
ipd$los[ipd$los == 0] <- 1 # replace 0 -> 1

ต่อมาก็ใช้ format.Date สร้าง month_id ซึ่งก็คือ ปี 4 digit และ เดือน 2 digit (’yyyy-mm’) เพราะโจทย์ของเราต้องการคำตอบเป็นรายเดือนครับ

# Make MonthID
ipd$month_id <- as.character(format.Date(ipd$dis_date, format = '%Y-%m'))

ขั้นถัดมาเป็นการ join ตารางเพื่อหาความหมายของ dchtype เพื่อดูว่าผู้ป่วยคนนี้ จำหน่ายออกอย่างไร โดยใช้ library dplyr และ left_join เข้าด้วยกัน ซึ่งคอมลัมน์ที่ใช้ join คือ dchtype เสร็จแล้วทำการเปลี่ยนชื่อคอลัมน์ name ให้เป็น dchtype_nm

# join data ipd + dtype
library(dplyr)
ipd <- left_join(ipd, dtype, by='dchtype')

# rename colName
names(ipd)[names(ipd) == "name"] <- "dchtype_nm"

มาดูตารางที่ได้กัน

ลอง plot bar ของ dchtype_nm สถิติประเภทการจำหน่ายผู้ป่วยกันสักนิ๊ดด

# สร้างตารางสรุปข้อมูลก่อน
df_count <- ipd %>%
  select(dchtype_nm, id) %>% 
  group_by(dchtype_nm) %>%  
  summarise(cnt_ = n_distinct(id))

# plot chart ดู
library(ggplot2)
ggplot(df_count, aes(x = dchtype_nm, y = cnt_)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  labs(title = "จำนวนผู้ป่วยตามประเภทการจำหน่าย",
       x = "ประเภทการจำหน่าย", y = "จำนวนผู้ป่วย") +
  theme_minimal()

จากรูป จะเห็นว่าผู้ป่วยส่วนใหญ่ถูกจำหน่ายออก โดยแพทย์อนุญาต (with approval) เป็นส่วนใหญ่เลย และ 2 เดือนที่ผ่านมามีผู้ป่วยเสียชีวิต (Death) 179 ราย นอกจากนี้ยังมีการย้ายไปโรงพยาบาลอื่น (Transfer) จำนวน 151 ราย

## histogram
hist(ipd$los)
hist(ipd$cost)

## boxplot
boxplot(ipd$los)
boxplot(ipd$cost)

ลองพล็อต 📈 ดูการกระจายตัวของจำนวนวันนอน (LoS) และค่าใช้จ่าย (Cost) ที่เกิดขึ้นดูครับ ค่อนข้างจะมี outlier ค่อนข้างสูงนะครับ นั่นคือการจะหาค่ากลางที่เป็นตัวแทน อาจจะต้องเปลี่ยนไปใช้ค่า Median ที่ทนต่อค่า outlier ที่สูงแบบนี้

ถึงตอนนี้เราก็ได้ตารางที่เราเตรียมข้อมูลไว้ส่วนนึงแล้ว พร้อมต่อการนำไปคำนวณหารายงานตามโจทย์ที่ได้มาแล้วครับ ซึ่งกว่าจะได้รายงาน อาจจะต้องมีการ join ข้อมูลเพิ่ม, คำนวณคอลัมน์ใหม่ หรือสร้างกราฟสลับกันไปมานิดนึงนะครับ พร้อมแล้วไปต่อเลย✌️

Data report & analysis

จากโจทย์ของพี่พยาบาล “จำนวนผู้ป่วยทั้งหมดที่รับไว้ในหอผู้ป่วยในเดือนนี้เป็นเท่าใด?” มาลองทำดูกัน

ward_stat <- ipd %>%
  select(id, ward, month_id, dchtype_nm, los, cost) %>%
  group_by(ward, month_id) %>%
  summarise(
    cnt_pt = n_distinct(id),
    sum_los = sum(los),
    sum_cost = sum(cost),
    mean_los = mean(los),
    mean_cost = mean(cost),
    mead_los = median(los),
    mead_cost = median(cost))
ward_stat

จากโค้ดเป็นการใช้ pipe operator เราจะเลือก select เฉพาะคอลัมน์ที่จะใช้ขึ้นมา จากนั้นจัดกลุ่มตาม group_by ward และ month_id เพื่อหาค่า Summaries ต่างๆ คือ

• cnt_pt = n_distinct(id) จำนวนผู้ป่วยนับแบบไม่ซ้ำตาม id
• sum_los = sum(los) ผลรวมจำนวนวันนอน
• sum_cost = sum(cost) ผลรวมค่าใช้จ่ายในการรักษาพยาบาล
• mean_los = mean(los) ค่าเฉลี่ยจำนวนวันนอน
• mean_cost = mean(cost) ค่าเฉลี่ยค่าใช้จ่ายในการรักษาพยาบาล
• mead_los = median(los) ค่า Median ของจำนวนวันนอน
• mead_cost = median(cost) ค่า Median ของค่าใช้จ่ายในการรักษาพยาบาล

ได้ข้อมูล Summaries ออกมาประมาณนี้ครับ ที่นี้เรามา join กับตาราง beds เพื่อเอาจำนวนเตียงแต่ละหอผู้ป่วยไปคำนวณอัตราการครองเตียงต่อ

# join data ward_stat + beds
ward_stat <- left_join(ward_stat, beds, by='ward')

ต่อด้วยการหาจำนวนวันในแต่ละเดือน โดยเราจะประกาศ data frame mday ขึ้นมาก่อน จากนั้นค่อย join นะครับ

# df บอกว่าแต่ละเดือนมีกี่วัน
mday <- data.frame(
  month_id = c('2024-10', '2024-11'), 
  days = c(31, 30)
)

# join data ward_stat + mday
ward_stat <- left_join(ward_stat, mday, by='month_id')

ตามมาด้วยการคำนวณ “อัตราการครองเตียง” 🧮 ซึ่งเค้าใช้สูตรนี้

# calc อัตราการครองเตียง
# (ผลรวมจำนวนวันนอนของผู้ป่วยใน x 100 ) / (จำนวนเตียงของโรงพยาบาล x(จำนวนวันในเดือนนั้น)))
ward_stat <- ward_stat %>%
  mutate(bed_rate = round((sum_los*100)/(beds * days), digit=2)) 
ward_stat

😎 ได้มาแล้วอัตราการครองเตียง ซึ่งเราจะรู้การเพิ่มขึ้นหรือลดลงแค่ไหนเมื่อเทียบกับเดือนที่แล้ว โดยการพล็อตลง chart ซึ่งหนุ่มขอ assume แสดงเฉพาะหอผู้ป่วย “2301” เป็นตัวอย่างนะครับ

## assume ว่าต้องการแสดง ward = 2301
ward2301 <- 
  ward_stat %>%
  filter(ward == 2301)

## bar chart จำนวนผู้ป่วยใน
library(ggplot2)
ggplot(ward2301, aes(x = month_id, y = cnt_pt)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = cnt_pt), vjust = -0.5) +
  labs(title = "จำนวนผู้ป่วยที่รับไว้ในหอผู้ป่วยในแต่ละเดือน",
       x = "เดือน", y = "จำนวนผู้ป่วย") +
  theme_minimal()

## bar chart อัตราการครองเตียง
ggplot(ward2301, aes(x = month_id, y = bed_rate)) +
  geom_bar(stat = "identity", fill = "lightgreen") +
  geom_text(aes(label = bed_rate), vjust = -0.5) +
  labs(title = "อัตราการครองเตียงของหอผู้ป่วยในแต่ละเดือน",
       x = "เดือน", y = "อัตราการครองเตียง") +
  theme_minimal()

เรามีการประกาศตัวแปรใหม่ ward2301 ขึ้นมา แล้วใช้ library ggplot2 ขึ้นมาพล็อต bar chart แสดงจำนวนผู้ป่วยที่ admit และเทียบอัตราการครองเตียงได้ดังนี้

หรือจะพล็อตเป็น Line chart ก็ได้

## line chart
ggplot(ward2301, aes(x = month_id, y = cnt_pt, group = 1)) +
  geom_line(color = "skyblue", size = 1) +
  geom_point(color = "skyblue", size = 3) +
  geom_text(aes(label = cnt_pt), vjust = -0.5) +
  labs(title = "จำนวนผู้ป่วยที่รับไว้ในหอผู้ป่วยในแต่ละเดือน",
       x = "เดือน", y = "จำนวนผู้ป่วย") +
  scale_y_continuous(limits = c(50, NA)) +
  theme_minimal()

จริงๆ เราสามารถ export file .csv หรือ รูปกราฟต่างๆ ไปใช้งานได้เลย แต่หนุ่มขอแวะไปพล็อต scatter plot หาความสัมพันธ์ระหว่างค่าใช้จ่ายกับจำนวนวันนอนดูนิดนึง ใช้ค่า median มาพล็อตดู ได้โค้ดดังนี้

## อยากรู้การกระจายตัวของ จำนวนวันนอน + ค่าใช้จ่ายที่เกิดขึ้น
# make scatter plot 
ggplot(ward_stat, aes(mead_los, mead_cost)) +     
  geom_point() + 
  geom_smooth()

ผลที่ได้ทั้ง 2 ตัวแปรก็มีความสัมพันธ์กันอยู่ที่ยิ่งจำนวนวันนอนรักษาพยาบาลนานขึ้น → ค่าใช้จ่ายก็ยิ่งเพิ่มสูงขึ้น

ทีนี้หนุ่มสงสัยที่หอผู้ป่วยไหนน้า 🧐 ที่จำนวนวันนอนสูงเกิน 15 วันบ้าง ก็ทำการ filter และ join ดูนิดนึงว่าเป็นหอผู้ป่วยไหนและทำไมกัน

## เอ๊ะใครกันนะ
los_hg <- ward_stat %>% filter(mead_los > 15)

# df ward
ward_nm <- data.frame(
  ward = c('2418', '2455', '2458', '303041001'), 
  ward_nm = c('วิกฤตศัลยกรรมประสาท', 'วิกฤตศัลยกรรมหัวใจ', 'วิกฤตศัลยกรรมหัวใจ', 'โรคหลอดเลือดสมอง')
)

los_hg <- left_join(los_hg, ward_nm, by='ward')
los_hg

หลังจากไปสืบมาทั้ง 4 รายนี้ พักอยู่ที่หอวิกฤตศัลยกรรมและหอ Stroke ที่นอนเกิน 15 วันนั้น เกิดจากต้องต่อคิวเข้ารับการผ่าตัด และหลังผ่าตัดอาจจะมีการนอน observe ดูอาการด้วยทำให้จำนวนวันนอนและค่าใช้จ่ายก็เพิ่มสูงขึ้นไปด้วย

---

เย้ๆ ก็ได้ลองใช้ภาษา R ในการจัดการและทำรายงานประจำเดือนแล้วนะครับ ก็เป็นภาษานึงที่โค้ดและใช้งานได้ง่าย กราฟที่พล็อตออกมาก็สวยดีด้วย ถ้าได้มีโอกาสใช้งาน R ในการทำสถิติอีกครั้งจะเอาวิธีการมาฝากนะครับ 😁