Tag: RStan

Stan กับ ปัญหา linear regression

จากที่เคยไปโม้ไว้เยอะว่า Stan มันเจ๋งยังไงให้กับเพื่อนร่วมงานฟัง เมื่อวันศุกร์ที่แล้วผมก็เลยนัดโชว์แบบคร่าวๆพร้อมกับสาธิตการใช้งานกับปัญหาง่ายๆอย่าง linear regression ให้เค้าดูกันครับ ก็เลยคิดว่าจะเอาที่ไปโม้ไว้มาใส่ไว้ที่นี่ด้วยเผื่อว่าใครสนใจใช้งาน Stan แต่ไม่รู้จะเริ่มตรงไหน หวังว่าคงมีประโยชน์กันบ้าง Stan จัดว่าเป็น Probabilistic Programming Language อันหนึ่งที่ช่วยให้เราทำโมเดลที่เกี่ยวข้องกับความน่าจะเป็นในแบบที่ต้องอาศัยทฤษฎีเบย์ได้ง่ายขึ้น และปัญหาที่เราสามารถเอา Stan มาใช้ประโยชน์ก็เช่นด้าน optimization ครับ ซึ่งโปรแกรมด้านนี้มันก็มีหลายตัวครับเช่น WinBUGS/OpenBUGS, JAGS และอื่นๆอีกมาก ลองดูเพิ่มเติมที่ http://probabilistic-programming.org/wiki/Home ครับ ตัวอย่างที่เอามาโชว์นี้ก็มาจากที่ผมเคยพูดไว้แล้วที่  http://www.sakngoi.com/2017/03/22/probabilistic-programming-language/ เรื่องที่เราจะหาค่า จากข้อมูลที่เรามีโดยที่โมเดลของเราเป็นสมการเส้นตรง อันนี้ตัวอย่างcode R สำหรับสร้างข้อมูลที่เอามาfitกับ Stan ครับ set.seed(12345) fakedata<-function(alpha, beta, t){ err<-rnorm(length(t),mean = 0, sd = 5) fake<-data.frame(time=t,y=alpha + beta*t +err) return(fake) } ob<-fakedata(alpha = 1.5, beta = 0.5,

เหตุผลที่ทำไมกินยาartesunateทุก12ชม.ถึงไม่ช่วยเรื่องการรักษามาลาเรีย

เรื่องมันมีอยู่ว่าแบบจำลองเรื่องการดื้อยาartemisininที่เสนอโดย Saralamba et al. (ผมและทีมงาน) ทำนายว่าการเพิ่มขึ้นของระยะเวลาในการกำจัดเชื้อมาลาเรียจากร่างกายคนไข้เป็นผลมาจากการที่ระยะหนึ่งของเชื้อมาลาเรียที่เรียกว่า Rings หรือวงแหวนนั้นมีการตอบสนองกับยาลดลงพร้อมกับเสนอว่าถ้าเปลี่ยนการให้ยาคนไข้ทุก 24 ชม.มาเป็น ทุก12ชม.จะช่วยลดเวลาของการกำจัดเชื้อได้ดีกว่า แต่จากการศึกษาโดย Das et al. ที่ทดลองให้ยาทุก 12 ชม. และเพิ่มขนาดเป็นของยากับคนไข้มาลาเรียพบว่าระยะในการกำจัดเชื้อไม่เปลี่ยนแปลงมากนัก จากผลที่ได้นี้มันแสดงให้เห็นว่าแบบจำลองที่เสนอไปนั้นต้องมีอะไรผิดพลาดแน่นอนและก็เลยทำให้มีการศึกษากันเพิ่มขึ้นพร้อมกับมีการเสนอสมมุติฐานต่างๆเพิ่มขึ้น ในแบบจำลองอันเดิม เราสมมุติว่าหลังจากที่เชื้อโดนยาแล้วจะตายพร้อมกับถูกกำจัดออกจากร่างการเลย อย่างรวดเร็วโดยม้าม และเป็นเพราะด้วยเหตุผลนี้ที่ทำให้มีการเถียงกันเรื่องการใช้ระยะของการกำจัดเชื้อเป็นตัวบอกถึงการดื้อยาartemisinin ว่าอาจไม่เหมาะสม (https://www.ncbi.nlm.nih.gov/pubmed/26239987) เพราะนั่นหมายความว่าถ้าเชื้อมันยังตอบสนองกับยาได้ดีอยู่ล่ะ แต่ม้ามมันกำจัดเชื้อที่ตายแล้วออกจากร่างการได้ช้ามันก็ดูเหมือนว่าเป็นเชื้อที่ดื้อยาถ้ายังใช้ระยะเวลาในการกำจัดเชื้อเป็นตัวบอกถึงการดื้อยา  แต่มันก็มีหลายงานวิจัยที่แสดงให้เห็นว่าม้ามนั้นกำจัดเชื้อมาลาเรียได้เร็วมากในผู้ป่วยที่รักษาด้วยยาartemisinin (https://www.ncbi.nlm.nih.gov/pubmed/28231817,https://www.ncbi.nlm.nih.gov/pubmed/11992295) ผมและทีมงานก็เลยเปลี่ยนวิธีที่เชื้อมาลาเรียจะตอบสนองกับยาในแบบจำลองโดยแทนที่เชื้อจะตายเลยหลังจากที่โดนยาก็เปลี่ยนมาเป็นเจ็บแทนแล้วค่อยๆตาย และก็ให้ว่าถ้าเชื้อเจ็บอยู่ถึงแม้จะโดนยาอีกก็ไม่ได้มีผลเปลี่ยนแปลงอะไรเลย ผลที่ได้จากแบบจำลองก็ดูเหมือนจะสอดคล้องกับทุกผลที่ได้จาก Das et al. นั่นก็คือถึงแม้จะให้ยาทุก 12 ชม.ก็ไม่ได้ลดระยะในการกำจัดเชื้อสักเท่าไหร่ กราฟด้านล่างนี้เป็นตัวอย่างผลที่ได้จากแบบจำลอง สีชมพูแสดงการลดลงของเชื้อในคนไข้ที่กินยาทุก 24 ชม.และสีฟ้าสำหรับคนไข้ที่กินยาทุก12ชม. ซึ่งแสดงให้เห็นว่าไม่ได้แตกต่างกันเท่าไหร่และกราฟนี้ก็แสดงให้เห็นว่าถ้าเราทำการศึกษาแบบที่ Das et al. ทำ โอกาสที่เราจะได้ว่าระยะของการกำจัดเชื้อแบบกินยาทุก24ชม.กับทุก12ชม.จะไม่แตกต่างกันมีสูงถึงเกือบ 60% ในการศึกษานี้ผมใช้ Stan ในการสร้างแบบจำลองและเปรียบเทียบผลที่ได้กับข้อมูลจริงจากคนไข้ครับ

ใช้ประโยชน์จาก rstan

ถ้าใครได้ใช้ deSolve สำหรับแก้ปัญหาเชิงตัวเลขของระบบสมการเชิงอนุพันธ์(ode) แล้วจะเห็นว่ามันมีข้อจำกัดเรื่องความเร็วอยู่ เพราะถ้าเขียน function ของสมการด้วย R แล้วจะช้ามาก ยิ่งถ้ามีการเอาfunction ไปใช้ fitting ด้วยเทคนิคอย่าง Markov chain monte carlo แล้วนี่เลิกคิดได้เลยสำหรับระบบใหญ่ ๆ deSolve เองก็พยายามแก้ปัญหานี้ โดยผู้ใช้สามารถเรียกใช้ function ที่เขียนใน C/C++ ที่คอมไพล์แล้วมาใช้ได้ แต่ก็ยังดูยากพอสมควรสำหรับคนที่ไม่เคยเขียน C/C++ เลย แต่มันก็ยังพอมีวิธีที่ดีกว่าหน่อยที่อยากจะแนะนำสำหรับคนที่อยากจะเพิ่มความเร็วสำหรับใช้กับ deSolve  นั่นก็คือเขียน ในภาษาของ rstan แทนครับ โดยหลักการก็คือเขียนสมการใน rstan แล้วใช้คำสั่งของ rstan ที่เรียกว่า expose_stan_functions เพื่อที่จะคอมไพล์ function แล้วเราก็จะยังเรียกใช้ได้ใน .GlobalEnv นั่นหมายความว่าเราสามารถเรียกใช้ function ที่คอมไพล์แล้วนี้กับ deSolve ได้ครับ ตัวอย่างเช่น ผมเขียนสมการ lorenz ในไฟล์ชื่อ lorenz.stan

Stan’s Program Blocks

หลักการของ Stan คือจะใช้ตัวแปรหรือ function อะไรก็ต้อง define ก่อน ไม่เหมือนกับ Bugs บล็อคที่จำเป็นสำหรับ Stan คือ model{} ส่วนอันอื่นเป็นแค่ตัวเสริม แต่ถ้าบล็อคตัวเสริมถูกใส่เข้ามาในโมเดลก็จะต้องเรียงตามลำดับที่แสดงไว้ในรูปคือ เริ่มจาก functions, data, transformed data, parameters, transformed parameters, model และ generated quantities

MCMC ใน R

ถ้าจะทำ Markov Chain Monte Carlo Monte Carlo (MCMC) ใน R นั้นนอกจากจะเขียนแล้วก็ยังมีโปรแกรมช่วยอีกหลายตัวครับ โดยโปรแกรมที่เป็นที่นิยมกันก็ได้แก่ WinBUGS, OpenBUGS, Jags, และก็ Stan ครับ โดย R มี package ที่ช่วยให้เราส่งผ่านหรือรับข้อมูล/โมเดล ระหว่าง R กับโปรแกรม MCMC เหล่านี้ครับ เช่น rjags, RStan, R2WinBUGS, R2OpenBUGS, BRugs ส่วน package ที่ช่วยในการวิเคราะห์ก็อย่างเช่น coda