Saturday, June 13, 2020

การสร้างแบบจำลองโรคติดเชื้อ (Infectious Disease Modelling) : กรณีตัวอย่าง COVID-19 ในประเทศไทย

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


เนื่องจากการเพิ่มขึ้นของการติดเชื้อมีลักษณะไม่ต่างจากการเพิ่มขึ้นของประชากรของสิ่งมีชีวิต ซึ่งมีลักษณะเป็น Exponential ไม่ใช่เส้นตรง การสร้างแบบจำลองจึงอาศัยวิธี Ordinary Differential Equation เพื่อการสร้างแบบจำลอง โดยสร้างสมการ Differential ที่แสดงถึงอัตราการเพิ่มขึ้นต่อหน่วยเวลา ซึ่งในที่นี้ใช้หน่วยเป็นหนึ่งวัน และ solve หา solution หรือสมการของจำนวนประชากรหรือการติดเชื้อ ณ เวลาต่างๆ ในลำดับต่อไป โชคดีที่มีซอฟต์แวร์ที่สามารถช่วยในการนี้ ซึ่งในกรณีนี้ใช้ซอฟต์แวร์ ODEINT เพื่อใช้ในการ solve Ordinary Differential Equation ซอฟต์แวร์ ODEINT เป็น C++ library ซึ่งถูกนำไปใช้ในหลายภาษา ในการศึกษาครั้งนี้เลือกใช้ภาษา Python ในการใช้งาน


แบบจำลองโรคติดเชื้อในกรณีนี้ได้ปรับปรุงมาจากผลงานของ Henri Froese (Froese, 2020) ซึ่งพัฒนาต่อมาจาก basic SIR model โดยแบบจำลองแบ่งประชากรออกเป็น 6 compartments ประกอบด้วย Susceptible, Exposed, Infected, Recovered, Critical และ Dead โดยมีความสัมพันธ์กันดังรูปข้างล่าง


โดย

Susceptible (S) คือกลุ่มประชากรที่ยังไม่ติดเชื้อและสามารถติดเชื้อได้

Exposed (E) คือกลุ่มประชากรที่ได้รับเชื้อแล้วแต่ยังไม่แพร่เชื้อ

Infected (I) คือกลุ่มผู้ติดเชื้อที่สามารถแพร่เชื้อได้

Recovered (R) คือกลุ่มที่หายจากการติดเชื้อและไม่แพร่เชื้อ

Critical (C) คือกลุ่มผู้ติดเชื้อที่มีอาการหนักและต้องการบริการ intensive care

Dead (D) คือกลุ่มผู้ติดเชื้อที่เสียชีวิต


เมื่อรวมทุกกลุ่มเข้าด้วยกันจะเท่ากับจำนวนประชากรทั้งหมด (N) เมื่อเริ่มการระบาดจะต้องมีอย่างน้อยสองกลุ่มคือ S และ I โดยกลุ่ม I จะเริ่มแพร่เชื้อไปยังกลุ่ม S และจะเปลี่ยนแปลงต่อไปดังแผนภูมิข้างต้นในอัตราที่ต่างๆ กัน ซึ่งอาจจะเป็นไปในทิศทางเพิ่มการระบาดคือจำนวนผู้ติดเชื้อเพิ่มขึ้น หรือลดการระบาดคือจำนวนผู้ติดเชื้อลดลงก็ได้


การหาอัตราการเพิ่มของประชากรในกลุ่มต่างๆ


เนื่องจากการจะสร้างแบบจำลองที่เป็น Ordinary Differential Equation การเขียนสมการที่แสดงอัตราการเปลี่ยนแปลงต่อวันของประชากรในกลุ่มนั้นๆ เป็นเรื่องที่สำคัญมากที่สุด โดยหลัก สมการจะประกอบไปด้วย 3 ส่วน คือ หนึ่ง อัตราหรือ rate ของ การเปลี่ยนแปลงต่อประชากรหนึ่งคนต่อหนึ่งวัน สอง โอกาสความน่าจะเป็นหรือ probability และ สาม จำนวนประชากรหรือ population ซึ่งสามารถเขียนเป็นสมการได้คือ


rate * probability * population


กรณีที่หนึ่ง กลุ่ม Susceptible (S)

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


วิธีของ Henri Froese จะคิดว่า ถ้าโดยเฉลี่ยผู้ติดเชื้อกลุ่ม I หนึ่งคนในหนึ่งวันจะแพร่เชื้อที่สามารถติดต่อได้ β ครั้งของการสัมผัสกับบุคคลอื่น (contacts) ดังนั้นถ้ามีกลุ่ม Infected จำนวนเท่ากับ I ดังนั้นในหนึ่งวันจะมี β*I contacts ที่สามารถแพร่เชื้อได้เกิดขึ้น หรือคือค่า population จะเป็นจำนวนของ contacts ทั้งหมด แต่การติดเชื้อจะเกิดขึ้นได้เฉพาะถ้าสัมผัสกับกลุ่ม S เท่านั้น ดังนั้นถ้ากลุ่มนี้มีจำนวนเท่ากับ S โอกาสที่จะสัมผัสกับกลุ่มนี้จึงมีค่า probability เท่ากับ S/N และเนื่องจากการติดเชื้อเกิดขึ้นโดยทันทีเมื่อเกิดการสัมผัส ค่า rate ซึ่งจะใช้ตัวย่อ k_S_E จึงเท่ากับ 1 ดังนั้นเราจึงสามารถเขียนสมการการเปลี่ยนแปลงในหนึ่งวันของกลุ่ม S (dSdt) ได้คือ


dSdt = -k_S_E * S/N * β*I


สมการเป็นค่าลบแสดงถึงการเปลี่ยนแปลงในทิศทางที่ลดลงของกลุ่ม S


อย่างไรก็ตาม มีข้อสงสัยว่าค่า probability ไม่ได้เท่ากันในแต่ละครั้งของ contact เช่นถ้า contact ครั้งที่หนึ่งและติดเชื้อ ค่า probability จะเท่ากับ S/N แต่ contact ครั้งที่สองค่า probability จะลดลงเหลือ (S-1)/N เนื่องจากจำนวน S ได้ลดลงไปแล้วหนึ่งคน และหากติดเชื้อครั้งที่สาม สี่ และเรื่อยไป ค่า probability ก็จะลดลงไปเป็นลำดับ ดังนั้นค่า probability จึงไม่ได้เท่ากันในทุกครั้ง


ในการศึกษาครั้งนี้จึงได้เสนอวิธีใหม่ โดยให้ population คือ S และคิดหาค่า probability ที่ S หนึ่งคนจะติดเชื้อได้ โดย S หนึ่งคนในหนึ่งวันจะมี contact โดยเฉลี่ยที่สามารถติดเชื้อได้เท่ากับ β ครั้ง และเขาจะติดเชื้อได้ถ้าสัมผัสประชากรในกลุ่ม I อย่างน้อยหนึ่งครั้งใน β ครั้ง ดังนั้นถ้าประชากรทั้งหมดมีผู้ติดเชื้อที่แพร่เชื้อได้เท่ากับ I จะสามารถคำนวณหาค่า probability นี้ได้เท่ากับ1-(1-I/N)จึงสามารถเขียนสมการการเปลี่ยนแปลงของ S ในหนึ่งวันได้ดังนี้


dSdt = -k_S_E * (1 - (1 - I/N)**β) * S


อย่างไรก็ตาม ค่า β ซึ่งคือค่าเฉลี่ยของจำนวน contacts ต่อคนต่อวันที่สามารถติดเชื้อได้ ค่านี้อาจเปลี่ยนแปลงได้ เช่นกรณีที่รัฐบาลออกมาตรการ lockdown หรือ social distancing จะทำให้ค่า β มีค่าลดลง ดังนั้นค่านี้จึงจะเปลี่ยนแปลง ณ ประมาณวันที่มีมาตรการออกมา ดังนั้นจึงได้กำหนดสมการการเปลี่ยนแปลงของ β หรือ beta(t) ไว้ดังนี้


beta(t) = (beta_start - beta_end) / (1 + np.exp(-k*(-t + x0))) + beta_end


โดย

beta_start คือ ค่า β ก่อนมาตรการ lockdown

beta_end คือ ค่า β หลังมาตรการ lockdown

k คือ อัตราการเปลี่ยนแปลง

x0 คือ ค่า t ณ วันที่เริ่มมาตรการ lockdown

np.exp คือ ค่า exponential ยกกำลัง


ดังนั้นสมการการเปลี่ยนแปลงของ S ในหนึ่งวัน คือ


dSdt = -k_S_E * (1 - (1 - I / N)**beta(t)) * S (1)


กรณีที่สอง กลุ่ม Exposed (E)

กลุ่มนี้เปลี่ยนแปลงเข้าจากการติดเชื้อของกลุ่ม S และเปลี่ยนแปลงออกเป็นกลุ่ม I หลังผ่านระยะเวลาตั้งแต่ได้รับเชื้อจนสามารถแพร่เชื้อได้ ซึ่งอาจใกล้เคียงกับระยะฟักตัว (Incubation period) ซึ่งสมมติว่าเท่ากับ P วัน หรืออีกนัยหนึ่ง E หนึ่งคนจะเปลี่ยนเป็น I ในอัตรา 1/P คนต่อวัน ซึ่งจะใช้ตัวย่อว่า k_E_I ดังนั้นในหนึ่งวัน E จะเปลี่ยนเป็น I เท่ากับ k_E_I * 1 * E โดยค่า probability เท่ากับ 1 ดังนั้นสมการการเปลี่ยนแปลงในหนึ่งวันของกลุ่ม E (dEdt) เมื่อพิจารณาทั้งขาเข้าและออกจะเท่ากับ


dEdt = k_S_E * (1 - (1 - I / N)**beta(t)) * S - k_E_I * E (2)


จะสังเกตุได้ว่าสมการขาเข้าจะเป็นบวกเนื่องจากทำให้ E เพิ่มขึ้น ส่วนสมการขาออกจะเป็นลบเพราะทำให้ E ลดลง


กรณีที่สาม กลุ่ม Infected (I)

กลุ่มนี้เกิดจากการเปลี่ยนแปลงเข้าจากกลุ่ม E และเปลี่ยนแปลงออกเป็นกลุ่ม C หรือ R สมมติว่าผ่านระยะเวลาตั้งแต่แพร่เชื้อได้จนหายป่วย สมมติว่าเท่ากับ X วัน และกลายเป็นกลุ่ม R ดังนั้น rate ที่จะเป็นกลุ่ม R ต่อคนต่อวัน จะเท่ากับ 1/X ซึ่งจะใช้ตัวย่อ k_I_R ในขณะเดียวกันหากผ่านระยะเวลาตั้งแต่แพร่เชื้อได้แต่มีอาการหนัก สมมติว่าเท่ากับ Y วัน และกลายเป็นกลุ่ม C ดังนั้น rate ที่จะเป็นกลุ่ม C ต่อคนต่อวัน จะเท่ากับ 1/Y และใช้ตัวย่อ k_I_C


เนื่องจากกลุ่ม I สามารถเปลี่ยนแปลงเป็นกลุ่ม R หรือ C ดังนั้น ถ้าโอกาสที่จะเป็นกลุ่ม C เท่ากับ p_I_to_C โอกาสที่จะเป็นกลุ่ม R จึงเท่ากับ 1 - p_I_to_C สมการการเปลี่ยนแปลงในหนึ่งวันของกลุ่ม I (dIdt) จึงเท่ากับ


dIdt = k_E_I * E - k_I_C * p_I_to_C * I - k_I_R * (1 - p_I_to_C) * I (3)


กรณีที่สี่ กลุ่ม Critical (C)

กลุ่มนี้เปลี่ยนแปลงเข้าจากกลุ่ม I และเปลี่ยนแปลงออกเป็นสองกลุ่ม คือกลุ่ม D หากเสียชีวิต และกลุ่ม R หากหายป่วย โดยมีโอกาสที่จะเสียชีวิตหรือเป็นกลุ่ม D เท่ากับ p_C_to_D ดังนั้นโอกาสที่จะหายจึงเท่ากับ 1 - p_C_to_D และเช่นเดียวกับกลุ่มที่แล้ว หากรู้ระยะเวลาตั้งแต่ป่วยหนักจนเสียชีวิต จะสามารถหา rate ซึ่งจะใช้ตัวย่อ k_C_D และหากรู้ระยะเวลาตั้งแต่ป่วยหนักจนหาย ก็จะสามารถหา rate การหายได้เช่นกัน ซึ่งจะใช้ตัวย่อ k_C_R ดังนั้นสมการการเปลี่ยนแปลงในหนึ่งวันของกลุ่ม C (dCdt) จึงเท่ากับ


dCdt = k_I_C * p_I_to_C * I - k_C_D * p_C_to_D * C - k_C_R * (1 - p_C_to_D) * C


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


Beds(t) = beds_0 + s * beds_0 * t


โดย

beds_0 คือจำนวนเตียงเริ่มต้นที่มี และ s คืออัตราการเพิ่มของเตียง


ดังนั้นสมการการเปลี่ยนแปลงในหนึ่งวันของ C จึงเป็น


dCdt = k_I_C * p_I_to_C * I - k_C_D * p_C_to_D * min(Beds(t), C)

- max(0, C - Beds(t)) - k_C_R * (1 - p_C_to_D) * min(Beds(t), C) (4)


กรณีที่ห้า กลุ่ม Recovered (R)

กลุ่มนี้เปลี่ยนแปลงเข้าจาก I และ C และเป็นกลุ่มสิ้นสุดไม่มีเปลี่ยนแปลงออก การเปลี่ยนแปลงเข้าจะเป็นตามสมการเดิม จาก I->R และ C->R เพียงเปลี่ยนเครื่องหมายจากลบเป็นบวก ดังนั้นสมการการเปลี่ยนแปลงในหนึ่งวันของ R จึงเป็น


dRdt = k_I_R * (1 - p_I_to_C) * I + k_C_R * (1 - p_C_to_D) * min(Beds(t), C) (5)


กรณีที่หก กลุ่ม Dead (D)

กลุ่มนี้เปลี่ยนแปลงเข้าจาก C และเป็นกลุ่มสิ้นสุดไม่มีการเปลี่ยนแปลงออกเช่นกัน การเปลี่ยนแปลงเข้าจะเป็นตามสมการเดิม จาก C->D เพียงเปลี่ยนเครื่องหมายจากลบเป็นบวก ดังนั้นสมการการเปลี่ยนแปลงในหนึ่งวันของ D จึงเป็น


dDdt = k_C_D * p_C_to_D * min(Beds(t), C) + max(0, C-Beds(t)) (6)


การกำหนดค่าพารามิเตอร์และค่าเริ่มต้น


จากทั้ง 6 สมการ จะเห็นได้ว่ามีค่าพารามิเตอร์ที่ต้องกำหนด และในการศึกษานี้ได้กำหนดไว้ ดังนี้


ค่าพารามิเตอร์ในการคำนวณ rate การเปลี่ยนแปลงระหว่างกลุ่ม

ระยะเวลาเฉลี่ยตั้งแต่ได้รับเชื้อจนแพร่เชื้อได้ (E->I) เท่ากับ 3 วัน ดังนั้นค่า k_E_I = 1/3

ระยะเวลาเฉลี่ยตั้งแต่แพร่เชื้อได้จนหายป่วย (I->R) เท่ากับ 9 วัน ดังนั้นค่า k_I_R = 1/9

ระยะเวลาเฉลี่ยตั้งแต่แพร่เชื้อได้จนมีอาการหนัก (I->C) เท่ากับ 12 วัน ดังนั้นค่า k_I_C = 1/12

ระยะเวลาเฉลี่ยตั้งแต่มีอาการหนักจนเสียชีวิต (C->D) เท่ากับ 7.5 วัน ดังนั้นค่า k_C_D = 1/7.5

ระยะเวลาเฉลี่ยตั้งแต่มีอาการหนักจนหายป่วย (C->R) เท่ากับ 6.5 วัน ดังนั้นค่า k_C_R = 1/6.5

*ข้อมูลข้างต้นได้จากการรวบรวมและศึกษาการระบาด COVID-19 ของ Henri Froese


ค่า probability ในการเปลี่ยนแปลงระหว่างกลุ่ม

prob_I_to_C ค่า probability ที่ I จะเปลี่ยนเป็น C

prob_C_to_D ค่า probability ที่ C จะเปลี่ยนเป็น D



ค่าพารามิเตอร์ในการคำนวณจำนวนเตียง ICU

beds_0 คือจำนวนเตียงเริ่มต้นที่มี ซึ่งข้อมูลนี้ได้จากการประมวลข้อมูลจากรายงานข้อมูลทรัพยากรสาธารณสุข

ประจําปี 2561 ของกระทรวงสาธารณสุข




ที่มา: รายงานข้อมูลทรัพยากรสาธารณสุข ประจำปี 2561 กองยุทธศาสตร์และแผนงาน กระทรวงสาธารณสุข


ซึ่งมีจำนวนเตียง ICU ทั่วประเทศ 7,718 เตียง เมื่อคิดอัตราครองเตียง 80% จะเหลือเตียง ICU ว่าง 1,543.6 เตียง หรือ 2.212 ต่อประชากรแสนคน


s อัตราการเพิ่มของเตียง


ค่าพารามิเตอร์ในการคำนวณค่า β

beta_start ค่า β ก่อนมาตรการ lockdown

beta_end ค่า β หลังมาตรการ lockdown

k อัตราการเปลี่ยนแปลงของค่า β

x0 วันที่เริ่มใช้มาตรการ lockdown


กำหนดค่าเริ่มต้น

N จำนวนประชากร 66,600,000 คน

S = N-1, E = 1, I = 0, C = 0, R = 0, D = 0


การทำ Model Fitting


การทำ Model Fitting ทำเพื่อให้แบบจำลองที่สร้างขึ้น มีความใกล้เคียงกับความเป็นจริงให้มากที่สุด โดยการปรับแต่งค่าพารามิเตอร์บางค่า โชคดีอีกเช่นกันที่มีซอฟต์แวร์ LMFIT ช่วยในการนี้ ในการศึกษานี้แบ่งพารามิเตอร์ในหัวข้อที่แล้วเป็นสองส่วน ส่วนแรกมีค่าคงที่ (fixed) ซึ่งผู้ศึกษาจะกำหนดค่าเอง และส่วนที่สองแปรผันได้ โดยจะให้ซอฟต์แวร์ LMFIT ทำการ fitting และเลือกค่าที่ดีที่สุดมาให้ ส่วนที่สองนี้ประกอบด้วยค่าพารามิเตอร์ดังต่อไปนี้ beta_start, k, x0, beta_end, prob_I_to_C, prob_C_to_D และ s โดยจะใช้วิธี least square ในการทำ fitting


การทำ Model Fitting นี้ได้เลือกใช้ข้อมูลจริงของจำนวนการเสียชีวิตสะสมจาก COVID-19 ในแต่ละวันของไทย โดยกำหนดให้วันเริ่มต้นของการระบาดคือวันที่ 3 มกราคม 2563 ที่เริ่มเฝ้าระวังคัดกรองและรายงาน ซึ่งจากรายงานสถานการณ์การระบาดของกระทรวงสาธารณสุข พบผู้ป่วยยืนยันรายแรกในรายงานของวันที่ 16 มกราคม 2563 เหตุที่เลือกใช้จำนวนผู้เสียชีวิต เพราะเห็นว่าน่าจะมีความถูกต้องมากกว่าจำนวนผู้ติดเชื้อที่ได้รับการตรวจยืนยัน


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





[[Fit Statistics]]

# fitting method = leastsq

# function evals = 1181

# data points = 151

# variables = 6

chi-square = 128.079039

reduced chi-square = 0.88330371

Akaike info crit = -12.8594734

Bayesian info crit = 5.24420563

## Warning: uncertainties could not be estimated:

s: at initial value

[[Variables]]

beta_start: 0.36594311 (init = 0.355)

k: 4.15 (fixed)

x0: 85.5726555 (init = 85)

beta_end: 1.6864e-05 (init = 0.1)

prob_I_to_C: 0.01922841 (init = 0.01)

prob_C_to_D: 0.01673402 (init = 0.05)

s: 0.01000000 (init = 0.01)



อย่างไรก็ตามเพื่อให้ LMFIT ทำ fitting ได้ผลดี ผู้ศึกษาได้ปรับให้ ค่าพารามิเตอร์ k เป็นแบบ fixed และให้เท่ากับ 4.15 ซึ่งผลจากการทำ model fitting มีข้อสังเกตุดังต่อไปนี้ หนึ่ง ค่า β ก่อน lockdown มีค่าประมาณ 0.36 และหลัง lockdown มีค่าประมาณ 0.00002 ซึ่งเข้าใกล้ศูนย์ แสดงว่ามาตรการ lockdown ได้ผลดี และ สอง ค่า x0 หรือวันที่เริ่มใช้มาตรการ lockdown ได้ค่าประมาณวันที่ 85 หรือวันที่ 29 มีนาคม 2563 ในขณะที่ไทยเริ่มพรก.ฉุกเฉินในวันที่ 25 มีนาคม 2563 และ curfew ในวันที่ 3 เมษายน 2563 ซึ่งนับว่าใกล้เคียงความเป็นจริง


การทำการพยากรณ์


เมื่อได้ค่าพารามิเตอร์ทั้งแบบที่กำหนดเองและแบบที่ให้ LMFIT ช่วยหาให้แล้ว ก็จะทำการพยากรณ์ไปข้างหน้า โดยจะทำการพยากรณ์ไปข้างหน้า 500 วันนับจากวันที่กำหนดให้เป็นวันเริ่มการระบาด (3 มกราคม 2563) ได้ผลออกมาดังนี้


หนึ่ง จากจำนวนประชากรทั้งหมด 66.6 ล้านคน เมื่อหยุดการระบาดโดยดูจากกลุ่ม R พบว่าจะมีประชากรที่ติดเชื้อและหายดีทั้งหมดประมาณ 270,000 ถึง 280,000 คน นอกจากนี้จำนวนกลุ่ม E และ I มีแนวโน้มลดลงเกือบทันทีหลังเริ่มมาตรการ lockdown (ประมาณวันที่ 85) และลดลงเป็นลำดับ โดยค่า beta (β) ลดลงจากประมาณ 0.36 เป็นเข้าใกล้ศูนย์ ดังรูป







สอง จำนวนประชากรกลุ่มผู้ติดเชื้อและแพร่เชื้อได้ (I) จะหมดไป ประมาณวันที่ 250 หรือประมาณหลังเดือนสิงหาคม 2563 เป็นต้นไป และมีจำนวนผู้เสียชีวิตไม่เกิน 60 คน ดังรูป




สาม จำนวนผู้ป่วยวิฤติ C ในแต่ละวัน สูงสุดประมาณไม่เกิน 1,000 รายและไม่เกินจำนวนเตียง ICU ที่รองรับได้ ดังรูป







รุป

การศึกษาครั้งนี้เป็นการใช้แบบจำลองเพื่อจำลองการติดเชื้อโดยใช้การติดเชื้อ COVID-19 ในประเทศไทยเป็นกรณีศึกษา การศึกษาใช้วิธี Ordinary Differential Equation และใช้ซอฟต์แวร์ ODEINT เพื่อการจำลอง และ LMFIT เพื่อการทำ Model Fitting ผลการศึกษาพบว่ามาตรการ lockdown มีผลให้การระบาดลดลงและคาดว่าผู้ติดเชื้อจะหมดไปหลังเดือนสิงหาคม 2563


อ้างอิง

Froese, H.(2020). Infectious Disease Modelling: Understanding the models that are used to model Coronavirus. Retrieved from https://towardsdatascience.com/infectious-disease-modelling-part-i-understanding-sir-28d60e29fdfc


Source code


Friday, December 6, 2019

SkinDx: Machine learning android app for pigmented skin lesion diagnosis with HAM10000 dataset

Summary
SkinDx is an android mobile app for pigmented skin lesion diagnosis. It is standalone machine learning app that uses models trained with HAM10000 dataset for diagnosis prediction without connecting to any server. The app can acheive 76.7% in accuracy.

Background
HAM10000 is a dermatoscopic image dataset for machine learning. It is created by Philipp Tschandl and others in 2018. Its aim is to solve the problem of small size and lack of diversity of available dataset of dermatoscopic images. There are 10,015  rows in the CSV file of this dataset. Each row describes a patient episode of skin disease, and each column describes an attribute eg. age, sex, localization, diagnosis and image file name of skin lesion. This dataset also contains 10,015 image files of related skin lesion.

The diagnosis is categorized into 7 groups: (1) actinic keratoses and intraepithelial carcinoma/bowen's disease (akiec), (2) basal cell carcinoma (bcc), (3) benign keratosis-like lesions (solar lentigines/seborrheic keratoses and lichen-planus like keratoses) (bkl), (4) dermatofibroma (df), (5) melanoma (mel), (6) melanocytic nevi (nv) and (7) vascular lesions (angiomas, angiokeratomas, pyogenic granulomas and hemorrhage) (vasc). The information from this dataset is used for developing android mobile app for diagnosis prediction of pigmented skin lesion.

Methods
1. HAM10000 dataset is checked for missing values and drop them out and then the remaining data is randomly divided into training and testing groups with ratio 80:20.
2.The trained model for classifying lesion image is firstly developed and tested by using transfer training from Mobilenet model.
3. The trained model from above is used to convert lesion image data into new 7 attribute columns in CSV file.
4. One-Hot encoding is used to convert categorical data.
5. The second trained model for diagnosis prediction is then developed and tested.
6. Both trained model are converted to Tensorflow lite model and are used for Android app development.

Results








The app is tested on android mobile device  with 1,992 unseen testing cases for 1st order prediction. The results are as follow:

Accuracy: 76.7% (1527/1992).


*Accuracy: correct prediction / all cases
*Sensitivity: TP / (TP + FN) *Precision: TP / (TP + FP)
*TP: True positive *FP: False positive *FN: False negative 

Discussion
1. Practical use of this app depends mainly on precision score for each diagnostic group. For example if the app predicts that the skin lesion is melanoma, there is a chance of 50.3% to be correct.
2. If using only lesion image without age, sex and localization data, only 69.7% in accuracy can be achieved.

Reference
Tschandl, P., Rosendahl, C. & Kittler, H. The HAM10000 dataset, a large collection of multi-source dermatoscopic images of common pigmented skin lesions. Sci. Data 5, 180161 doi:10.1038/sdata.2018.161 (2018).

Get it on Google Play

Monday, June 4, 2018

Project scheduling with Monte Carlo simulation on mobile

Background
Uncertainty is a common situation in project scheduling. It happens because we lack 
experience from doing some activities before or have some uncontrollable factors on 
the project. The classic technique we normally use is PERT (Program Evaluation and 
Review Technique). But there are some drawbacks with this technique especially when 
project has activities which are nearly equal in duration and also run in parallel. Nowadays 
the technique that should be used is Monte Carlo simulation. However using simulation 
requires at least desktop computer for constructing network model and experimenting 
with it. Microsoft Project can help us constructing network model. But we need  another 
special software for performing simulation with Microsoft Excel, e.g., RiskAMP, @Risk
or without Microsoft Excel, e.g., Full Monte. In this article we will show you the way easier 
by using AlmanacSoft SchemeSim app. Only one app for doing all these steps and also 
on mobile device.


What is AlmanacSoft SchemeSim (ASS)
ASS is android app that performs Monte Carlo simulation for project scheduling. It answers 
three important questions: firstly, when project will be completed, secondly, when each 
activity will start and end, and finally, the critical probability for each activity in the project. 
ASS automatically creates network model and do simulation within one step and produces 
scheduling report as an output. ASS supports two types of random generating for activity 
durations: Triangular and Pert distribution. It detects how many CPU core the mobile has 
and uses all of them for parallel native simulation. That is necessary because mobile has 
limited resource. To perform simulation efficiently needs high performance computing by 
creating native (binary) code and run it in concurrency. ASS also supports calendar time 
and working day and time. That means we can specify which day in a week we actually 
work, e.g., Monday to Friday,  and in which time interval, e.g., 8:30 - 16:30.

ASS overview and installation
ASS can be download and install from google play store. ASS offers in-app purchases. 
That means by default this app is in basic edition and can be upgraded to professional 
edition by paying subscription fee. In basic edition it can work for free with project not 
more than five activities. You also have a free trial period for 30 days after upgrading to 
professional edition. During this trial period, if you don’t satisfy with it you can cancel 
subscription with no charge.


When you open the app you will see this following screen. This is the Input page for 
type in your project data. 


There are other two pages, e.g., Report and About, you can navigate to them when you tap on three line menu. After simulation has finished, report will be generated in Report page.


 

 Project data preparation

You have to prepare project data as table above before inputing. In this case our project consists of 5 activities. We estimate activity duration by using three-point estimating (a = the best-case, 
b = the worst-case, m = the most likely) which a <= m <= b. We also analyze for predecessor of each activity.

Entering data

Type in number of activities (5) and click add button.

If you type wrong, you can check the box on the right upper side of each activity and click remove button, If you want to add more activities, you can do it again as many as you want. Type in activity name, distribution, a, b, m and predecessor. For distribution you can choose Triangular or Pert distribution. If you have no experience of this activity before, choose Triangular is better. But if you have it,  Pert distribution is better. If you have predecessor more than one activities, type activity name separated by comma with no space.

Now select project start date and time by click the right button. If you don't select it will use current date & time as a project starting point. In this case I choose 1 October this year and at 9 AM.
For confidence level, I choose 99% and use 100,000 iterations for simulation.


Simulation and report generation
You can save input file by click on three dot menu and select save input file for later use. Now it is time to perform simulation, click on the pink button and click run.




This project will take 20.258 days with confidence level 99% and will be finished on Sun, 21 October 2018 at 11:29. There are also detail of starting and ending date&time for each activity. Activity d is the most critical with probability 97.604% and the lowest critical is activity e which is 2.396%. You can see that we use just only 2.07 seconds for 100,000 iterations.

Tuesday, March 28, 2017

Asio ThreadPool Performance Test

Updated: 7 April 2017

Category: C++, Concurrency Programming, Thread pool, Asio

I have tested Asio ThreadPool performance by comparing processing time of three concurrency methods. The first method is Standard C++ Multithreading (mthread), the second is Asio ThreadPool (asio) and the last is Microsoft Parallel Patterns Library (ppl). All of them were used for calculating PI and Fibonacci numbers on the same machine for both Windows and Linux. They all used the same code except for PPL that is only available on Windows. I have also measured serial processing times for baseline comparison. The C++ compilers used are MSVC 19.10.25017 for Windows, GCC 6.3.1 and Clang 3.9.1 for Linux. The results are as follow:



Comparing time used between MSVC, GCC and Clang




GCC is the fastest in all scenarios (not including ppl). The second is MSVC except for fibonacci with asio that Clang comes the second.
 

Comparing each concurrency method with serial processing


Using MSVC
There is no big difference in time used for PI among mthread, asio and ppl. However for Fabonacci, mthread and ppl approximately process at the same speed and faster than asio which is the slowest.





Using GCC
There is no big difference between mthread and asio in time used for PI. But for Fabonacci, mthread has a little bit faster than asio.




Using Clang
There is no big difference between mthread and asio for both PI and Fibonacci.


Conclusions
  • The concurrency processing speed depends mainly on three factors: compiler, computational type and concurrency method respectively. 
  • Asio ThreadPool does not perform well with MSVC on some computational type (Fibonacci) when comparing with GCC and Clang.
The code is here.


Saturday, March 4, 2017

Using Asio C++ library based ThreadPool class

Category: C++, Concurrency Programming, Thread pool, Asio

Prerequisites: C++, Asio C++ programming concept, Multithreading

Requirement: Asio C++ library

What is it?
It is only a C++ header file that defines ThreadPool class based on Asio C++ Library.

Why is it created?
It is created to serve the following purposes:
  1. To use C++ multithreading with thread pool.
  2. To use task based concurrency programming.
  3. Can use strictly sequential invocation of handlers.
  4. Can use C++ Exception handling.

How can you use it?
  1. Create the client class that will use thread pool.
  2. Create any tasks that you want to execute in sequential or parallel order. In our case we use sequential for simplicity.
  3. If you want to handle exception. Modify previous code, add final task to start and stop MainIoService as shown in the following code.
  4. In the main code, create ThreadPool instance. Then create client instance by using shared_ptr. Make a call from the client and try to throw an exception in the client code to test exception handling.
  5. If you use _threadPool.strand() instead of _threadPool.enqueue(), any tasks called can not be executed concurrently.
How to dowload the ThreadPool Class?
The ThreadPool Class is in ThreadPool.h header file on github.


Thursday, September 12, 2013

Introduction to AlmanacSoft Payer 1.0

Category: Windows Store app, PayPal, REST API, eCommerce

What is AlmanacSoft Payer 1.0?

AlmanacSoft Payer 1.0 is Windows Store app running on Windows 8.1 Preview or later. It is payment app that uses new PayPal's REST APIs with standards-based technologies such as OAuth and JSON for paying money on the Internet. It supports both Direct Credit Card Payment and PayPal Account Payment. The app is native code written in C++/CX and developed by Poom Malakul Na Ayudhya.



How to Install
  1. Download AlmanacSoft Payer 1.0 and you will get the file named "Payer_1.0.0.0_Win32_Test.zip".
  2. Extract it and right click on "Add-AppDevPackage.ps1" to Run with PowerShell. 
  3.  You may be asked to acquire the developer license if you don't have it yet. Use your Microsoft  account to log in and get the license for free. It will be expired in one month and you can renew it.
  4. You may be asked for Execution Policy Change, you have to answer "[Y] Yes"
  5. You may be asked for installing the signing certificate, you have to answer "[Y] Yes". 
  6. AlmanacSoft Payer will be installed successfully.

How to use

For Merchants
  1. You have to register at PayPal for PayPal merchant account.
  2. Log in at https://developer.paypal.com/ with your merchant account. Then create an application to get merchant's Client Id and Secret.
  3. Run AlmanacSoft Payer and select credentials page. Use your Client ID and Secret to apply and export your encrypted merchant data file (MDF).
  4. Send encrypted MDF to your customers by e-mail or let them download from your website.
  5. Tell your customers to use AlmanacSoft Payer to import or download your MDF and use it to make payment for you.
For Customers
  1. Run AlmanacSoft Payer. If you are behind proxy server, set your proxy credentials first.
  2. Using AlmanacSoft Payer to import or download encrypted MDF provided by the merchant you want to pay.
  3. Select appropriate payment method for merchant in the countries supported by PayPal.

Wednesday, March 27, 2013

C++ AMP: How fast is it?

Category: Windows Store app, C++ AMP, GPU Programming
Prerequisites: C++, C++ AMP

Full Text in Thai (PDF 1.13 MB)

This study measures time used in millisecond for calculating square matrix multiplication at different dimension sizes starting from 256x256 to 2048x2048. The C++ AMP tested engines are two GPUs (Intel HD Graphics 4000 and NVIDIA Geforce GT-650M) and one software engine (Microsoft Basic Render Driver).  Two C++ AMP methods (simple and tiling) are used.  The study also measures time used by normal sequential code for using as a baseline comparison. The testing software is C++ Windows Store app running on Intel i7 RAM 8 MB.


The figure above shows Windows Store app used in this study. It also has 3D rotating cube in background for testing with DirectX.



This is the result table of time used measured in milliseconds. The table also shows computed ratio (in red) comparing between MS Basic Render Driver and Sequential code and also between both GPUs and MS Basic Render Driver. The size means the matrix dimension starting from 256x256 to 2048x2048.


From the above figure, MS Render Driver speed is around five to ten times comparing with sequential code when using simple method and five to twenty times when using tile method.


When using simple method, NVIDIA's speed is ten to thirty times comparing with MS Render Driver while Intel's speed is around five times comparing with MS Render Driver.


When using tile method, NVIDIA's speed is from fifteen to twenty five times comparing with MS Render Driver while Intel's speed is around five times comparing with MS Render Driver.