การสร้างโมเดล 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, γ)

SIR Diagram (Source: Wikipedia)

ในที่นี้ ผมขอจำลองสถานการณ์ COVID-19 ด้วยสมมติฐานดังต่อไปนี้

1. คนที่หายจากโรคนี้จะไม่กลับมาเป็นซ้ำ
2. เนื่องจากช่วงเวลาที่เกิดโรคระบาดนี้โดยทั่วไปจะพิจารณาช่วงสั้นๆ เลยทำให้เราสามารถลดตัวแปรของอัตราการเกิดและอัตราการตายลงได้

เพื่อความง่ายของการอธิบายของบทความนี้ ผมจะใช้ solution ของ ODE ไปเลยซึ่งจะได้ว่า (ศึกษาเพิ่มเติมได้จากลิ้งค์ข้างล่าง)

สรุปอัตราการเปลี่ยนถ่ายประชากรในกลุ่มต่างๆ

โดยที่ S = # Susceptible group, I = # Infected group, R = # Recovered group, N = S+ I + R.

สำหรับผู้อ่านที่เป็นสาย Hardcore สนใจรายละเอียดสามารถศึกษาเพิ่มเติมรายละเอียดของโมเดลนี้ได้จาก

The SIR Model for Spread of Disease - The Differential Equation Model

หรือชิลๆหน่อยก็ Wikipedia ครับ

Compartmental models in epidemiology

ถ้าท่านผู้อ่านเคยติดตามข่าว หรือเคยดูหนังเรื่อง Contagion อาจจะเคยได้ยินการพูดถึงสัมประสิทธิ์การแพร่กระจาย (R₀) ซึ่งเราสามารถคำนวณจาก β, γ เช่นกัน โดยการแปลความหมายนั้นถ้า R₀ > 1 นั่นหมายถึงจะเกิดการระบาด นั่นคือ อัตราการติดต่อเร็วกว่าอัตราการหายจากโรค ยิ่ง R₀ มาก ยิ่งแปลว่าโรคนี้ระบาดหนัก

สมการสำหรับการหาค่า R₀

ผมได้ลองใช้โมเดล SIR และ optimizer ใน scipy ช่วยแก้สมการหาค่า β, γ ที่เหมาะสมเพื่อที่จะประมาณค่าของ R₀ ของประเทศไทยว่าเป็นเท่าไร โดยใช้ข้อมูลจากเวปไซต์ HDX ครับ

ส่วนตัว loss function นั้นตอนนี้ลองเล่นอยู่ระหว่าง MAE กับ RMSE ซึ่งได้ผลต่างกันพอสมควรอยู่ครับ ใครอยากให้ลอง Loss function แบบไหนอีก บอกได้ใน comment เลยนะครับ

Loss function

ได้ผลอย่างไรไปชมกันเลยครับ…

ผลการ predict ของประเทศไทย โดยใช้ Loss function เป็น MAE

ผลการ predict ของประเทศไทย โดยใช้ Loss function เป็น RMSE…น่ากลัวมาก

จากกราฟข้างต้น สีดำคือข้อมูลที่เราเก็บได้จริง ซึ่งวันที่รันโมเดล คือวันที่ 1 เมษายน 2563 เลยมีข้อมูลเพียงแค่ 31 มีนาคม 2563, สีฟ้าคือกลุ่มคนที่มีโอกาสติดได้ ตอนนี้เซ็ตไว้ที่ 15000 คนซึ่งผมคิดว่ายังต้องปรับ parameter อีกเยอะ, สีแดงคือกลุ่มที่ติด,​ สีเขียวคือกลุ่มที่หายจากการติดโรคเเล้ว ส่วนแกน X เป็นแกนของเวลา ซึ่งผมได้ทำการพยากรณ์จนถึงประมาณเดือน มิถุนายน

ผลลัพธ์ที่โมเดลทำออกมานี่เห็นเเล้วเสียวๆเลยใช่ไหมครับว่าจะไปไกล ผมก็ไม่แน่ใจเหมือนกันว่าจะไปขนาดนั้นไหม แต่ขอให้ไม่ไกลไปกว่านี้เลย T.T

ค่าของ R₀ ของแต่ละ Loss function ได้ต่างกันครับ คือ R₀ ของ MAE จะอยู่ที่ ≈16 ส่วนของ RMSE นี่ไปไกลถึง 5000+ เลยเเหละครับ

คราวนี้เอาไว้เท่านี้ก่อน เดี๋ยวคราวหน้าคิดว่าจะมาทำต่อโดยใช้ข้อมูลของคนที่หาย (Recovery group) มาเป็น condition ในสมการด้วยครับ และจะลองนำผลลัพธ์ของคราวนี้ไปใส่กับประเทศอื่นๆต่อไป ถ้าสนใจกด Clap และ Follow ให้ด้วยนะครับ จะได้มีกำลังใจปั่นโมเดลเพิ่มเติมครับ

อันนี้เป็น Clip จาก 3Blue1Brown ที่พูดถึงครับ

Edit: ลิ้งค์นี้สำหรับ part 2 ครับ

การสร้างโมเดล SIR เพื่อทำนายการติดเชื้อ COVID-19 (Part 2)