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

ต่อจากของเดิมนะครับ ใครที่ยังไม่ได้อ่าน part 1 เชิญอ่านที่นี่ได้เลยครับ

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

ก่อนที่จะเริ่ม part2 ผมขออธิบาย part1 เพิ่มเติมสำหรับบางท่านที่อาจจะมีคำถาม

• จุดแรกที่อยากพูดถึงคือ โมเดล SIR นี้จะพยายามประมาณค่า (parameter estimate) β, γ โดย fit กับข้อมูลจำนวนผู้ติดเชื้อให้มากที่สุด โดย loss function ต่างๆ ซึ่งตอนสุดท้ายจะสามารถประมาณค่า R₀ ได้ โดยมอง β, γ ได้เหมือนกับเป็นอัตราการเปลี่ยนจากกลุ่ม S ไป I และ I ไป R ตามลำดับ ซึ่ง R₀ จะเหมือนเป็น index เพื่อบอกว่าโรคนี้จะเป็นโรคระบาดหรือไม่ โดยมี threshold อยู่ที่ R₀ > 1
• จุดที่สอง คือตัวโมเดลนี้ได้จำลองเหตุการณ์ว่าถ้าปล่อยให้เป็นไปโดยธรรมชาติ ไม่ได้มีการแทรกแซงด้วยมนุษย์หรือมาตรการจากภาครัฐหรือเอกชน เหตุการณ์จะเป็นอย่างไร
• จุดที่สามคือ ข้อมูลของประเทศไทยนั้น ส่วนตัวคิดว่าไม่มีทางที่ผู้ติดเชื้อจะเจริญเติบโตเป็น exponential เพราะว่าด้วยความสามารถของการตรวจผู้ติดเชื้อต่อวันยังจำกัดอยู่เพียงแค่ linear (อย่างมาก ~700เคส ต่อวัน) ซึ่งจุดนี้ต้องบอกอีกว่ากว่าจะได้ตรวจก็ต้องดูอาการนานมาก ซึ่งเหมือนเป็น Bias sample แต่ว่าทางผู้เขียนไม่ได้นำมาคิดในสมการ

มาคราวนี้ตามที่สัญญาว่าจะลองรันโมเดลกับประเทศอื่นๆด้วย เนื่องจากข้อมูลที่มี มีมากกว่า 180 ประเทศ แต่เพื่อความง่ายต่อการดูกราฟ ผมเลยเลือกมาเพียง 12 ประเทศ คือ Thailand, South Korea, Japan, China, Taiwan, Malaysia, US, Switzerland, Canada, Spain, Iran, Italy ข้อมูลที่ใช้ทำคราวนี้จะเป็นเหมือนคราวที่เเล้วนะครับ คือ As of 31 March 2020 สำหรับผู้อ่านท่านใดอยากให้ลองทำประเทศที่ไม่อยู่ในลิสนี้ก็บอกในคอมเม้นได้นะครับ จะลองรันให้ ผลลัพธ์จะเป็นอย่างไร ไปดูกันข้างล่างนี่เลยครับ

ครั้งที่แล้วผมได้สรุปว่า loss function = MAE (Mean absolute error) ให้ผลลัพธ์ของประเทศไทย ตัวนี้ให้ผลดีกว่า RMSE (Root mean square error) ซึ่งถ้าดูประเทศอื่นด้วยแล้ว จะเห็นว่าจำนวนกลุ่มผู้ติดเชื้อจะหักหัวลงโดยส่วนใหญ่ ซึ่งเป็นเพราะตอนที่ทางผู้เขียนทำโมเดลนั้น ได้ใส่ bound ของการ estimate parameter เอาไว้เพียงแค่ 0.4 ทำให้การประมาณค่าตัวแปรตกไปอยู่ที่ขอบของลิมิตที่ให้ไว้

Parameter estimation of 12 selected countries w.r.t. MAE

ซึ่งถ้าเราเปรียบเทียบกับ RMSE จะเห็นว่ามีหลายประเทศที่ infectious group ติดเส้นบนสูงสุดที่กำหนดไว้เป็น upper bound ของแต่ละประเทศ และไม่มีทีท่าว่า recovery group จะเพิ่มจำนวนขึ้นเลย นั่นเป็นเพราะว่าเราได้ใช้เฉพาะข้อมูลของจำนวนผู้ติดเชื้อเพียงอย่างเดียวมาใช้ในการประมาณค่าตัวแปร

Parameter estimation of 12 selected countries w.r.t. RMSE

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

ตัวแรกนี้เกิดจากการใช้ loss function = MAE

Parameter estimation using MAE

ส่วนตัวนี้เกิดจากการใช้ loss function = RMSE

Parameter estimation using RMSE

อย่างที่ผมได้เกริ่นไว้ตอนแรกเเล้วว่าได้จำกัดค่าของการ estimate β, γ ไว้ที่ 0.4 ซึ่งหลายๆประเทศได้ชนเพดานของลิมิตที่ผมให้ไว้เเล้ว ผลลัพธ์ที่ได้ออกมานี้ผลได้ลองเพิ่ม limit ให้มากกว่า 0.4 แล้วเหมือนกัน แต่ไม่เป็นที่ถูกใจผมสักเท่าไหร่ (ตอนเเรกเซ็ตไว้เท่านี้เพราะคิดว่าคอมพิวเตอร์ตัวเองอาจจะรันไม่ไหว search space ใหญ่เกิน)

ต้องบอกว่าหลังจากที่ผมทำโมเดล SIR ที่ใช้ข้อมูลผู้ติดเชื้อเสร็จ ก็ได้เกิดความคิดว่าน่าจะประมาณค่าได้ดีกว่านี้ และเนื่องจากเว็ปไซต์ HDX มีข้อมูลของผู้หายจากการติดเชื้อด้วย (Recovery group) ผลเลยได้ลองใส่ข้อมูลนี้ลงไปเป็น exogenous variable และเพิ่ม regularization ใน loss function ในการคำนวณค่านี้ด้วย … ปัญหาที่ตามมาคือจะถ่วงน้ำหนักอย่างไรระหว่าง Infectious data กับ Recovery data ซึ่งผมนั้นก็ได้ลองเปลี่ยนทั้ง loss functions, weights แต่สุดท้ายมาจบอยู่ที่ RMSE โดย weight อยู่ที่ประมาณ​ 0.65 ถ่วงให้ I group และ 0.35 ถ่วงให้ R group เหตุผลของการเลือกตัวนี้เป็นเพราะว่าผมต้องการประมาณค่าจำนวนผู้ติดเชื้อ เลยให้ความสำคัญกับ RMSE(Infectious data) มากกว่าครับ

0.65 * RMSE(Infectious data) + 0.35 * RMSE(Recovery data)

​เราไปดูผลลัพธ์กันว่าหน้าตาจะเป็นอย่างไร

Parameter estimation with Recovery data using RMSE for the loss function. Weight = 0.65 I + 0.35 R.

ตัวนี้ดูเข้าใกล้ความจริงมากขึ้นเนื่องจากเราใส่ค่าของคนที่หายจากโรคนี้เข้าไปด้วยเหมือนกัน

เรามาดูกันว่า R₀ จะเป็นเท่าไร

นอกจากนี้ผมยังได้ลอง loss function ตัวอื่นด้วย เช่น MAE, MAPE, SMAPE, customized loss แต่ผลลัพธ์ที่ได้ไม่ดีเท่าที่ควร แต่ไหนๆก็ทำมาแล้ว เอามาให้ดูเป็นตัวอย่างเล่นๆ หน้าตาจะเป็นอย่างไร ไปดูกันครับ

Parameter estimation with Recovery data using SMAPE for the loss function. Weight = 0.65 I + 0.35 R.

ตัวข้างบนนี้คือผลลัพธ์ของ SMAPE (Symmetric Mean Absolute Percentage Error) โดยที่

SMAPE definition (Image: Wikipedia)

A_t = actual value of observation t, F_t = forecast value of observation

ซึ่งผลลัพธ์ของจำนวนผู้ติดเชื้อส่วนใหญ่จะประมาณให้มีค่าที่มากเกินไป (overestimated) เพราะว่า สมมติค่าจริงคือ x แต่เราประมาณค่าจากโมเดลผิดจากค่าจริงไป a ซึ่งเป็นไปได้สองแบบคือ x-a กับ x+a โดยที่ x+a จะได้ % error ที่ต่ำกว่า x-a นั่นเลยทำให้โมเดลเรา bias ไปในทาง overestimation มากกว่า underestimation

สุดท้ายนี้ทางผู้เขียนอยากบอกว่าที่ทำมาทั้งหมดเป็นเพียงโมเดลที่ใส่ตัวแปรน้อยมาก ทำให้เข้าใจง่ายแต่อีกแง่หนึ่งคือไม่สามารถนำมาอธิบายสถานการณ์จริงได้ดีอย่างที่ต้องการ แต่ไม่ว่าโมเดลจะบอกว่าสถานการณ์ประเทศเราหรือทั้งโลกจะเป็นอย่างไร แต่ถ้าเราทุกคนร่วมมือคนละนิด ช่วยกันอยู่บ้านเพื่อลดการติดเชื้อ และไม่เพิ่มภาระให้บุคคลากรทางการแพทย์ เท่านี้เราก็จะดึง I group ในโมเดลให้ติดดินไม่เพิ่มตัวเลขกันเเล้วนะครับ ขออ่านผู้อ่านทุกท่านสุขภาพดี ผ่านวิกฤตินี้ไปด้วยกัน