รู้จักการเติบโตแบบ exponential และ logistic จาก COVID-19

บล็อกที่แล้ว เราเล่าเกี่ยวกับแนวโน้มของจำนวนผู้ติดเชื้อ COVID-19 ว่ายังคงเพิ่มขึ้นต่อเนื่องทั่วโลก และยังพูดถึงอัตราการเสียชีวิตและการแพร่กระจายของโรคเมื่อเทียบกับโรคอื่นๆ ไปแล้ว ในบล็อกนี้เราจะมารู้จักแบบจำลองทางคณิตศาสตร์ (โมเดล) อย่างง่ายสำหรับอธิบายการระบาดของโรค โดยใช้ exponential และ logistic function

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

Dataset & source code

ในโปรเจ็กต์นี้นำข้อมูลจำนวนผู้ติดเชื้อแยกตามประเทศระหว่างวันที่ 22 มกราคม 2020 จนถึงปัจจุบัน จาก github นี้ ซึ่งดึงข้อมูลผู้ติดเชื้อจาก Johns Hopkins University Center for Systems Science and Engineering มาจัดรูปแบบให้พร้อมใช้งาน

ในบล็อกนี้จะนำเสนอแนวคิดของแบบจำลองแต่ละแบบไม่ได้มาสอนเขียนโค้ดทีละขั้นตอน แต่โค้ดที่เขียนในบล็อกนี้ทั้งหมดรวมอยู่ใน github ซึ่งสามารถเข้าไปโหลดมารันเองได้เลยครับ

หาเส้นที่ฟิตที่สุด

และนี่คือรายงานจำนวนผู้ติดเชื้อรายใหม่ ของแต่ละวัน (ไม่ใช่จำนวนผู้ติดเชื้อสะสมนะ) ในประเทศไทย นับตั้งแต่วันที่ 22 มกราคม 2020 จนถึงวันที่ 21 มีนาคม 2020

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

แล้วเราจะลากเส้นอะไรดี?

ถ้าดูแนวโน้มจากรูปด้านบนเราคงบอกได้ว่าเส้นตรงสีน้ำเงิน คงไม่ใช่เส้นที่ดีนัก เพราะบอกแนวโน้มได้ใกล้เคียงกับความเป็นจริงแค่จุดเริ่มต้นกับจุดสุดท้าย เส้นสีเขียวก็ยังไม่ค่อยดี เพราะฟิตกับข้อมูลแค่ช่วงแรก แต่ช่วงท้ายก็แยกทางกันไป เส้นที่ดีน่าจะเป็นเส้นโค้งสีส้ม ที่ลากผ่านจุดสีแดงได้ใกล้เคียงกว่า หากเราพิจารณาธรรมชาติของโรคติดต่อดูแล้ว เราอาจเริ่มต้นจากสมมติฐานง่ายๆ ว่า ตอนแรกเรามีผู้ติดเชื้อคนเดียว จากนั้นก็แพร่ให้คนรอบตัว 2 คน เพราะฉะนั้นเมื่อเวลาผ่านไปผู้ติดเชื้อหนึ่งคน จะเพิ่มเป็น 2 จาก 2 เพิ่มเป็น 4 และจาก 4 เพิ่มเป็น 8 เราอาจกล่าวได้ว่า เมื่อมีการแพร่เชื้อ จะมีผู้ป่วยรายใหม่เพิ่มแบบทวีคูณ (exponential)

exponential growth

ดังนั้นเส้นที่เราลากขึ้นมาก็น่าจะเป็นฟังก์ชัน exponential ซึ่งฟังก์ชันที่เลือกใช้คือ y = eˣ หน้าตาเป็นแบบในรูปด้านล่าง ในช่วงแรกๆ การเติบโตจะเกิดขึ้นช้าๆ และจะเติบโตเร็วขึ้นเรื่อยๆ เมื่อเวลาผ่านไป (เราสามารถเลือกฐานของฟังก์ชันเป็นเลขอะไรก็ได้ที่มากกว่า 1 แต่ขอเลือกเป็นค่า e เพื่อความสะดวกในการดิฟหรืออินทิเกรต)

เราจะใส่พารามิเตอร์ a, b, c เข้าไปเพื่อปรับแต่งฟังก์ชันนี้ให้ฟิตกับข้อมูลของเราที่สุดดังนี้

เป้าหมายของเราคือการหาค่า a, b, c ที่สอดคล้องกับจำนวนผู้ติดเชื้อของเราที่สุด การปรับค่า a, b, c จะส่งผลต่อรูปร่างของเส้นโค้งดังนี้

ค่า a กำหนดว่าตัวคูณเป็นกี่เท่าของ a=1, b กำหนดอัตราการเติบโตของ y (b ยิ่งสูงยิ่งเติบโตเร็ว) ในขณะที่การเปลี่ยนค่า c จะสร้างเส้นโค้งที่หน้าตาเหมือนกันทุกประการ แต่เลื่อนไปทางซ้ายหรือขวาตามแกน x

การหาค่า a,b,c ที่ฟิตกับข้อมูลเราที่สุด อาจใช้วิธีที่เป็นที่นิยมกันอยู่แล้ว อย่าง least square estimation (LSE) ซึ่งทำได้โดยวัด error ระหว่างจุดจริงๆ กับค่าที่ทำนายได้ แล้วพยายามปรับ a,b,c ให้ได้ error น้อยที่สุด

แต่เนื่องจากค่า y (จำนวนผู้ติดเชื้อ) เป็นจำนวนนับ ในงานวิจัยเกี่ยวกับโรคระบาดส่วนใหญ่ จะหาค่าพารามิเตอร์ด้วยวิธี maximum likelihood estimation (MLE) บนการแจกแจงแบบ Poisson distribution หรือ negative binomial distribution มากกว่า

ในบล็อกนี้เราจะไม่ได้อฺธิบายหลักการของ MLE มีหลายบทความที่อธิบายเรื่องนี้ไว้แล้วเช่น What is Maximum Likelihood Estimation และ Maximum Likelihood Estimation For Regression ถ้ามีโอกาสจะเขียนเล่าเรื่องนี้เป็นบล็อกแยกอีกที หากใครมีข้อสงสัย หลังไมค์มาคุยกันได้ครับ

การทำ curve fitting ด้วยวิธี maximum likelihood estimation จะพิจารณาว่าค่า y ที่เราต้องการทราบมีการแจกแจงอย่างไร และพยายามหาพารามิเตอร์ที่สอดคล้องกับการแจกแจงนั้นที่สุด (source)

กลับมาที่ข้อมูลของเรากันต่อสำหรับข้อมูล 69 วันแรก เราสามารถหาเส้นโค้ง exponential ด้วยวิธี maximum likelihood estimation โดยถือว่าจำนวนผู้ติดเชื้อรายใหม่มีการแจกแจงแบบ Poisson เส้นโค้งที่ฟิตกับข้อมูลเราที่สุดออกมาเป็นดังนี้

จำนวนผู้ติดเชื้อรายวันค่าจริง (จุดเขียว) ค่าที่ทำนาย (เส้นฟ้า) พารามิเตอร์ a=0.025, b = 0.122, c= -3.467

กราฟนี้บอกกับเราว่าจำนวนผู้ติดเชื้อรายวันจะเพิ่มขึ้นเรื่อยๆ และถ้าลากเส้นนี้ต่อไปเรื่อยๆ จะพบว่าหากจำนวนผู้ติดเชื้อเติบโตตามเส้นโค้งนี้ จะมีผู้ติดเชื้อใหม่ต่อวันมากถึง 300 ราย ประมาณวันที่ 4 เมษายน และพุ่งทะยานสูงขึ้นเรื่อยๆ

นั่นทำให้กราฟจำนวนผู้ติดเชื้อสะสมหน้าตาเป็นแบบกราฟด้านล่างนี้ กราฟนี้บอกกับเราว่าจะประเทศไทยจะมีผู้ติดเชื้อสะสมสูงถึง 2,500 คน ในวันที่ 3 เมษายน สังเกตว่าเมื่อผู้ติดเชื้อรายวัน เพิ่มขึ้นแบบ exponential จะทำให้จำนวนผู้ติดเชื้อสะสมเพิ่มแบบ exponential ตามไปด้วย

จำนวนผู้ติดเชื้อสะสมในไทย ค่าจริง (จุดเทา) ค่าที่ทำนายได้ (เส้นฟ้า)

และถ้าเราลากเส้นจำนวน ผู้ติดเชื้อสะสม ต่อไปถึงประมาณกลางเดือนสิงหาคม จะมีผู้ติดเชื้อสะสมทั้งหมดประมาณครึ่งหนึ่งของประชากรโลก!

ผลการทำนายจำนวนผู้ติดเชื้อสะสมในไทยเมื่อยึดฟังก์ชัน exponential เป็นสรณะ

พักก่อน! การเติบโตแบบ exponential จะเกิดขึ้นแค่ช่วงแรกๆ ของการระบาดเท่านั้น และจะค่อยๆ ช้าลง เนื่องจากประชากรมีจำกัด ไม่ได้ติดต่อหรือมีปฏิสัมพันธ์กันทุกคน และมีมาตรการควบคุมโรค ทำให้การระบาดที่มาแรงในช่วงแรก จะแผ่วลงเมื่อเวลาผ่านไป สำหรับประเทศไทยเอง หากเราดูกราฟต่อจากวันที่ 69 มาอีกประมาณหนึ่งสัปดาห์ จะพบว่าจำนวนผู้ติดเชื้อรายวัน ไม่ได้เติบโตแบบ exponential อีกต่อไป

หากเรายังฝืนใช้ฟังก์ชัน exponential กับประเทศที่จำนวนผู้ติดเชื้อรายใหม่เริ่มลดลงแล้วเช่น จีน จะเห็นได้ชัดว่าค่าที่ทำนายผิดไปจากความจริงมาก เพราะจำนวนผู้ติดเชื้อรายใหม่ จะเพิ่มถึงจุดสูงสุดค่าหนึ่ง แล้วค่อยลดลงเรื่อยๆ เหมือนโค้งรูประฆังคว่ำ

จำนวนผู้ติดเชื้อใหม่รายวันในจีน (จุดเขียว (training set) และแดง (test set)) ค่าที่โมเดลทำนาย (เส้นฟ้า) พารามิเตอร์ a=104.02 , b = 0.024, c= -116.601

เมื่อเราพิจารณากราฟจำนวนผู้ติดเชื้อสะสม ของจีน เราจะพบว่ากราฟเริ่มคงที่และไม่ได้เติบโตแบบ exponential อีกต่อไป แต่เติบโตในลักษณะ s-curve

การเติบโตของจำนวนผู้ติดเชื้อสะสมในจีน ตั้งแต่วันที่ 22 มกราคม — 13 เมษายน 2020 จำนวนผู้ติดเชื้อสะสมในจีนไม่ได้เติบโตแบบ exponential อีกต่อไป แต่เติบโตในลักษณะ s-curve

เราเรียกการเติบโตแบบรูปตัว S นี้ว่า logistic growth

ฟังก์ชัน logistic

แล้วถ้าเราลองฟังก์ชันที่หน้าตาคล้ายๆ s-curve มาใช้ดูล่ะ? มีฟังก์ชันหนึ่งที่ดูจะเข้ากับอัตราการเติบโตที่เร็วในช่วงแรกและเริ่มช้าลงเรื่อยๆ เมื่อเวลาผ่านไป นั่นคือฟังก์ชัน logistic หรืออีกชื่อหนึ่งว่าฟังก์ชัน sigmoid

ฟังก์ชัน sigmoid เป็นฟังก์ชันที่ให้ค่าอยู่ระหว่าง 0 และ 1 มีการเติบโตแบบตัว S ช่วงที่ x < 0 จะเติบโตในอัตราเร่ง แต่พอ x>0 จะเติบโตช้าลง เรื่อยๆและลู่เข้าค่า 1

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

เราสามารถพิสูจน์ทางคณิตศาสตร์ได้ว่า ถ้าจำนวนผู้ติดเชื้อสะสมมีรูปร่างแบบเส้นโค้ง sigmoid จำนวนผู้ติดเชื้อรายใหม่ในแต่ละวัน จะเพิ่มขึ้นแบบระฆังคว่ำ (พิสูจน์ได้โดยเอาฟังก์ชัน sigmoid มาดิฟหนึ่งทีเพื่อดูอัตราการเปลี่ยนแปลงที่เวลาใดๆ) เราเรียกฟังก์ชันระฆังคว่ำนี้ว่า Logistic distribution function

ผลจากการดิฟ Logistic function หนึ่งครั้ง เราจะเห็นอัตราการเปลี่ยนแปลงในแต่ละช่วงเวลา Logistic distribution function จะฟอร์มตัวเป็นรูประฆังคว่ำ

กราฟนี้บอกกับเราว่า หากจำนวนผู้ติดเชื้อใหม่รายวันหน้าตาเพิ่มขึ้นด้วยแพตเทิร์นนี้ (แรกๆ เพิ่มช้า และเร็วขึ้นเรื่อยๆ จนถึงจุดสูงสุดค่าหนึ่ง แล้วค่อยๆ กลับมาช้าลง) จำนวนผู้ป่วยสะสมจะหน้าตาออกมาเป็น s-curve แบบฟังก์ชัน sigmoid ดังจะเห็นได้จากตัวอย่างที่จำลองขึ้นมานี้

(ซ้าย) หากจำนวนผู้ติดเชื้อใหม่เพิ่มขึ้นแบบ logistic distribution จำนวนผู้ติดเชื้อสะสม (ขวา)จะกลายเป็น sigmoid

เราสามารถใส่พารามิเตอร์ต่างๆ เข้าไปเพื่อปรับแต่งฟังก์ชันได้ดังนี้

ผลของการปรับค่าพารามิเตอร์ a, b, c ส่งผลต่อรูปร่างของเส้นโค้งดังนี้

หากเราตีความพารามิเตอร์ a, b, c ให้เข้ากับข้อมูลผู้ติดเชื้อรายวันของเราจะพบว่า
a กำหนดจำนวนผู้ติดเชื้อต่อวันสูงสุด (ยิ่ง a มาก ค่าสูงสุดยิ่งมากตาม)
b บอกว่าจำนวนผู้ติดเชื้อจะเพิ่มขึ้นเร็วแค่ไหน (ยิ่ง b มาก ยิ่งโตเร็ว)
c ระบุว่าการระบาดจะถึงจุดสูงสุด ณ วันที่ x=c (เมื่อแกน x แทนจำนวนวัน)

ลองกับประเทศไทยก่อน ให้ training set มี 69 วันเหมือนเดิม ได้ผลดังนี้

จำนวนผู้ติดเชื้อรายวันค่าจริง (จุดเขียว(training set) และแดง (test set)) ค่าที่ทำนาย (เส้นฟ้า) พารามิเตอร์ a=ุ582.772, b = 0.168, c= 71.069

กราฟที่ได้บอกกับเราว่า จำนวนผู้ติดเชื้อใหม่ต่อวัน จะถึงจุดสูงสุดที่เกือบ 150 คนต่อวัน ในวันที่ 1 เมษายน 2020 หลังจากนั้น จำนวนผู้ติดเชื้อรายวันจะลดลงเรื่อยๆ แม้จะเห็นได้ชัดว่าไม่ค่อยตรงกับค่าจริงๆ ใน test set (สีแดง)เท่าไหร่ แต่ก็ใกล้เคียงความจริงมากกว่าการเติบโตแบบ exponential

หากเราพิจารณาจำนวนผู้ติดเชื้อสะสมจะพบว่า จำนวนผู้ติดเชื้อสะสมในไทยจะลู่เข้าที่จำนวนประมาณ 3,500 คน

ถ้าดูจากแนวโน้มคร่าวๆ แล้ว จำนวนผู้ติดเชื้อสะสมน่าจะต่ำกว่าเส้นโค้งที่ทำนายได้เสียอีก

นี่เป็นนิมิตหมายอันดีว่าการระบาดจะยุติแล้วใช่มั้ย?

อาจจะใช่หรือไม่ใช่ก็ได้ เพราะรูปแบบการระบาดไม่จำเป็นต้องเติบโตแบบ ฟังก์ชัน sigmoid เสมอไป ดังเช่นกรณีของเกาหลีใต้

จำนวนผู้ติดเชื้อใหม่รายวันในเกาหลีใต้ (บน) ไม่เป็นตามเส้นโค้ง logistic distribution ซึ่งกำหนดด้วยพารามิเตอร์ a= 2381.637, b=0.269, c=41.98 ทำให้จำนวนผู้ติดเชื้อสะสม (ล่าง) ยังคงเพิ่มขึ้นอย่างต่อเนื่อง

ประเทศที่ดูจะมีการเติบโตแบบ sigmoid มากที่สุดคือประเทศจีน ซึ่งกว่าเราจะรู้ว่ามันเติบโตแบบ sigmoid จริงๆ มันก็ปรากฏรูปร่างแบบ sigmoid ให้เห็นแล้ว และไม่ได้การันตีด้วยว่าการระบาดจะไม่กลับมาเกิดซ้ำอีก

(บน) จำนวนผู้ติดเชื้อใหม่รายวันในจีน เส้นโค้งกำหนดด้วยพารามิเตอร์ a=16278.383,b=-0.193, c=19.1 (ล่าง) จำนวนผู้ติดเชื้อสะสมในจีน

การทำนายจะแม่นหรือไม่ยังขึ้นอยู่กับขนาดของ training set ด้วย หากเราทดลองแบ่งขนาด training set ออกเป็นช่วงไล่ตั้งแต่ 10,15,20,25,30 วัน ขอแสดงเฉพาะผลลัพธ์ของการทำนายจำนวนผู้ติดเชื้อสะสม เพื่อให้เห็นข้อผิดพลาดชัดเจน

ผลลัพธ์ของการทำนายจำนวนผู้ติดเชื้อสะสมที่ขนาด training set ต่างๆ

จะเห็นว่ามองดูเส้นไหนก็ฟิตกับข้อมูลอย่างน้อยช่วงหนึ่งเสมอ ปัญหาคือเราไม่มีทางรู้อนาคตได้เลยว่าเส้นที่เราทำนายออกมาจะฟิต หรือไม่ฟิตกับสิ่งที่เกิดขึ้นจริง

เราอาจพอสรุปได้ว่าโมเดลอย่างง่ายที่เล่ามาทั้งหมดนี้ อาจไม่สามารถใช้คาดการณ์สิ่งที่เกิดขึ้นในอนาคตได้เลย!

แม้จะไม่เหมาะกับการทำนายอนาคต แต่งานวิจัยบางงานก็นำฟังก์ชัน exponential หรือ logistic (หรือโมเดลที่ซับซ้อนกว่านี้) มาฟิตกับเส้นโค้งเมื่อการระบาดผ่านมาแล้ว และมีข้อมูลทั้งหมดแล้ว เพื่อย้อนกลับมาดูลักษณะการระบาด และหาค่าอื่นๆที่เกี่ยวข้องกับโรคระบาด เช่นค่า R₀ ระยะฟักตัวของโรค หรือแม้แต่จำนวนผู้ติดเชื้อที่ไม่ได้รายงาน (unreported cases) เป็นเท่าไหร่ ถ้าใครอยากศึกษาต่อ แนะนำส่วน introduction ของ paper ฉบับนี้

เราอาจทดลองฟังก์ชันนอกเหนือจากที่เล่ามานี้ เช่น ฟังก์ชันพหุนาม (polynomial) ซึ่งอธิบายการเติบโตของผู้ติดเชื้อที่โตช้ากว่าทวีคูณ (sub-exponential) หรือฟังก์ชันอื่นๆ ที่ใช้คาดคะเนการเติบโตของประชากร เช่น Von Bertalanffy function และ Gompertz function อีกวิธีการหนึ่งคือใช้โมเดลสำหรับ time series เช่น ARIMA ไปเลย

แต่ไม่ว่าจะเลือกฟังก์ชันอะไรมาใช้ก็ตาม ข้อจำกัดหนึ่งของการทำนายการระบาดของโรคด้วยการเอาฟังก์ชันที่มีอยู่แล้วมาฟิตกับจำนวนผู้ติดเชื้อแบบทื่อๆ นี้ เราจะได้พารามิเตอร์ที่ควบคุมลักษณะของเส้นโค้งมาเช่น ลู่เข้าค่าไหน ตัวคูณเท่าไหร่ พุ่งขึ้นเร็วแค่ไหน หักลงเมื่อไหร่ แต่ค่า a,b,c ที่เราหามาได้เหล่านี้แทบไม่ได้ช่วยให้เรารู้จักธรรมชาติของโรคดีขึ้นเลย

แบบจำลองโรคระบาดที่ดี ควรสร้างขึ้นจากความเข้าใจในตัวโรคระบาด โดยมีพารามิเตอร์ที่สะท้อนสมบัติต่างๆ ของโรคและพฤติกรรมของประชากรเอง เช่น ปฏิสัมพันธ์ของผู้ป่วยกับคนที่ยังไม่ป่วย การเดินทางของประชากร อัตราการแพร่เชื้อ ระยะฟักตัว อัตราการหายป่วยจากโรค การเกิดและตายของคนในสังคม ฯลฯ เราจึงมีอีกวิธีการหนึ่งที่ใช้จำลองการระบาดของโรค โดยแบ่งประชากรออกเป็นกลุ่ม แล้วพิจารณาว่าพารามิเตอร์ที่ว่ามาส่งผลต่อการเปลี่ยนแปลงของประชากรแต่ละกลุ่มอย่างไร เราเรียกแบบจำลองประเภทนี้ว่า compartmental model ตัวอย่างของโมเดลประเภทนี้คือ SIR และ SEIR model

เล่าอย่างนี้อาจจะไม่เห็นภาพ เราจะมาต่อกันเรื่อง SIR model กันในบล็อกหน้า

โปรดติดตามตอนต่อไป และรักษาสุขภาพด้วยครับ…