การสร้างโมเดล SIR เพื่อทำนายการติดเชื้อ COVID-19 (Part 1)
เนื่องจากช่วงนี้สถานการณ์ COVID-19 มาแรง รวมถึงมีการทำงานจากที่บ้าน ทำให้หลายๆคนได้เปลี่ยนแปลงพฤติกรรมการใช้ชีวิตประจำวัน ผมเองก็เป็นหนึ่งในนั้น จากปกติที่ต้องฝ่าวงรถติดกว่าจะถึงบ้านก็หมดเวลาแล้ว เลยไม่เคยได้มีเวลามานั่งเขียน Medium
วันนี้ผมจะเขียนเกี่ยวกับโมเดลของการจำลองการกระจายตัวของโรคระบาดด้วย SIR model (Susceptible-Infected-Recovered model) เนื่องจากผมเคยทำโปรเจคเกี่ยวกับโมเดลแนวนี้นิดหน่อยในห้องเรียน แต่นั่นนานมากมาแล้ว และไม่เคยคิดจะจับอีกเลย วันนี้ถือเป็นการดีได้มาปัดฝุ่นวิชาความรู้ที่เคยเรียนมา ประกอบกับผมชอบติดตาม Channel 3Blue1Brown บน YouTube อยู่แล้ว แล้วสัปดาห์ที่ผ่านมาเค้าได้ทำ Animation ที่ได้ลองหลายๆ scenario เลยคิดอยากลองเองบ้าง (ดูวิีดีโอข้างล่างโพสนี้ครับ)
ก่อนที่จะออกนอกเรื่องไปไกล เรามาเริ่มกันเลยดีกว่าครับ ก่อนอื่นต้องเริ่มจากไอเดียว่าเราจะแบ่งประชากรเป็น 3 กลุ่ม คือ a) กลุ่มที่มียังไม่ติดโรค (Susceptible), b) กลุ่มที่สามารถแพร่โรคต่อ (Infectious) และ c) กลุ่มที่หายจากโรคแล้ว (Recovered) ซึ่งตัวแปรที่สำคัญของโมเดลมี 2 ตัวคือ อัตราของการติดเชื้อจากกลุ่มที่ติดโรคสู่กลุ่มที่ยังไม่เป็น (Infection rate, β) กับ อัตราของผู้ติดเชื้อที่จะหาย (Recovery rate, γ)
ในที่นี้ ผมขอจำลองสถานการณ์ COVID-19 ด้วยสมมติฐานดังต่อไปนี้
- คนที่หายจากโรคนี้จะไม่กลับมาเป็นซ้ำ
- เนื่องจากช่วงเวลาที่เกิดโรคระบาดนี้โดยทั่วไปจะพิจารณาช่วงสั้นๆ เลยทำให้เราสามารถลดตัวแปรของอัตราการเกิดและอัตราการตายลงได้
เพื่อความง่ายของการอธิบายของบทความนี้ ผมจะใช้ solution ของ ODE ไปเลยซึ่งจะได้ว่า (ศึกษาเพิ่มเติมได้จากลิ้งค์ข้างล่าง)
โดยที่ S = # Susceptible group, I = # Infected group, R = # Recovered group, N = S+ I + R.
สำหรับผู้อ่านที่เป็นสาย Hardcore สนใจรายละเอียดสามารถศึกษาเพิ่มเติมรายละเอียดของโมเดลนี้ได้จาก
หรือชิลๆหน่อยก็ Wikipedia ครับ
ถ้าท่านผู้อ่านเคยติดตามข่าว หรือเคยดูหนังเรื่อง Contagion อาจจะเคยได้ยินการพูดถึงสัมประสิทธิ์การแพร่กระจาย (R₀) ซึ่งเราสามารถคำนวณจาก β, γ เช่นกัน โดยการแปลความหมายนั้นถ้า R₀ > 1 นั่นหมายถึงจะเกิดการระบาด นั่นคือ อัตราการติดต่อเร็วกว่าอัตราการหายจากโรค ยิ่ง R₀ มาก ยิ่งแปลว่าโรคนี้ระบาดหนัก
ผมได้ลองใช้โมเดล SIR และ optimizer ใน scipy ช่วยแก้สมการหาค่า β, γ ที่เหมาะสมเพื่อที่จะประมาณค่าของ R₀ ของประเทศไทยว่าเป็นเท่าไร โดยใช้ข้อมูลจากเวปไซต์ HDX ครับ
ส่วนตัว loss function นั้นตอนนี้ลองเล่นอยู่ระหว่าง MAE กับ RMSE ซึ่งได้ผลต่างกันพอสมควรอยู่ครับ ใครอยากให้ลอง Loss function แบบไหนอีก บอกได้ใน comment เลยนะครับ
ได้ผลอย่างไรไปชมกันเลยครับ…
จากกราฟข้างต้น สีดำคือข้อมูลที่เราเก็บได้จริง ซึ่งวันที่รันโมเดล คือวันที่ 1 เมษายน 2563 เลยมีข้อมูลเพียงแค่ 31 มีนาคม 2563, สีฟ้าคือกลุ่มคนที่มีโอกาสติดได้ ตอนนี้เซ็ตไว้ที่ 15000 คนซึ่งผมคิดว่ายังต้องปรับ parameter อีกเยอะ, สีแดงคือกลุ่มที่ติด, สีเขียวคือกลุ่มที่หายจากการติดโรคเเล้ว ส่วนแกน X เป็นแกนของเวลา ซึ่งผมได้ทำการพยากรณ์จนถึงประมาณเดือน มิถุนายน
ผลลัพธ์ที่โมเดลทำออกมานี่เห็นเเล้วเสียวๆเลยใช่ไหมครับว่าจะไปไกล ผมก็ไม่แน่ใจเหมือนกันว่าจะไปขนาดนั้นไหม แต่ขอให้ไม่ไกลไปกว่านี้เลย T.T
ค่าของ R₀ ของแต่ละ Loss function ได้ต่างกันครับ คือ R₀ ของ MAE จะอยู่ที่ ≈16 ส่วนของ RMSE นี่ไปไกลถึง 5000+ เลยเเหละครับ
คราวนี้เอาไว้เท่านี้ก่อน เดี๋ยวคราวหน้าคิดว่าจะมาทำต่อโดยใช้ข้อมูลของคนที่หาย (Recovery group) มาเป็น condition ในสมการด้วยครับ และจะลองนำผลลัพธ์ของคราวนี้ไปใส่กับประเทศอื่นๆต่อไป ถ้าสนใจกด Clap และ Follow ให้ด้วยนะครับ จะได้มีกำลังใจปั่นโมเดลเพิ่มเติมครับ
Edit: ลิ้งค์นี้สำหรับ part 2 ครับ