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/ เรื่องที่เราจะหาค่า \alpha, \beta จากข้อมูลที่เรามีโดยที่โมเดลของเราเป็นสมการเส้นตรง y = \alpha+\beta x

อันนี้ตัวอย่าง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, 1:100)
plot(ob)

# using lm to fit the fake data
lmfit<-lm(ob$y ~ ob$time)
cat(paste0("alpha = ",format(round(lmfit$coefficients[1],4),nsmall = 4), " beta = ", format(round(lmfit$coefficients[2],4),nsmall = 4)))
abline(lmfit, col='red')

ซึ่งจากข้อมูลที่สร้างมานี้เรามารถใช้คำสั่ง lm หาได้ว่าค่า \alpha, \beta เป็นเท่าไหร่โดยคำตอบที่ได้นี้มาจากเทคนิคที่เรียกว่า QR decomposition ซึ่งจะต่างจากเทคนิค Markov chain Monte Carlo  ที่ใช้ใน Stan ที่คำตอบจะออกมาเป็นในลักษณะของdistributionที่เราประมาณได้ว่าค่าที่เป็นไปได้อาจอยู่แถวค่า median ของdistributionที่ได้นั้น

อันนี้เป็นตัวอย่างcode ที่ผมใช้fitกับข้อมูลด้านบนครับ

library(rstan)
#เป็นfunctionที่สร้างสำหรับให้ค่าเริ่มต้นของตัวแปร
initf<-function(){
 list(alpha=0.1,beta=0.1,sig=0.1)
}
#สร้างตัวแปรdataใหม่ในรูปแบบของlist
dat <- list(y=ob$y,x=ob$time,np=length(ob$y))

#codeของstan
mymodel <- "
functions{
 real fn(real x, real alpha, real beta){
 return alpha + beta*x;
 }
}
data{
 int np;
 vector[np] y;
 vector[np] x;
}
parameters{
 real alpha; 
 real beta;
 real sig;
}
model{
 vector[np] md;

alpha ~ uniform(-5,5);
 beta ~ uniform(-5,5);
 sig ~ uniform(0,30);

for(i in 1:np)
 md[i] <- fn(x[i],alpha,beta);

y ~ normal(md,sig);
}
generated quantities{
 vector[np] y_pred;
 vector[np] md;

for(i in 1:np)
 md[i] <- fn(x[i],alpha,beta);

for(i in 1:np)
 y_pred[i] <- normal_rng(md[i], sig);

}
"
stanout <- stan(model_code = mymodel, data = dat, chains = 1, iter = 10000, init = initf)
summary(stanout,pars=c("alpha","beta"),probs=c(0.025,0.5,0.975))$summary

ตัวอย่างของoutputของค่าที่ได้ครับ

บางคนพอเห็นcodeของ Stanแล้วบอกว่าจะมาเสียเวลาเขียนทำไมยาวๆใช้ lm สิบันทัดเดียวจบเลย ก็จริงครับสำหรับปัญหาง่ายๆแบบนี้ แต่ถ้าปัญหาที่มันซับซ้อนกว่านี้ล่ะ ผมหมายถึงโมเดลที่มันมีหลายตัวแปรที่ซับซ้อนในการคิดและแน่นอนไม่ได้เป็นแบบ linear ล่ะ Stan มันถูกออกแบบมาให้ใช้กับปัญหาที่หลากหลายครับ แต่มันก็ต้องเปลี่ยนวิธีมองปัญหาให้เป็นแบบเบย์ (Bayesian) ครับ อย่างเช่นในปัญหานี้ ….

เดี๋ยวมาต่อ

 

เหตุผลที่ทำไมกินยา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

##lorenz.stan
functions{
  vector lorenz(real t,vector y,vector parms){
    vector[3] ydot;
    // parms[1]=beta ,[2]=sigma,[3]=rho
    
    ydot[1]=parms[2]*(y[2]-y[1]);
    ydot[2]=y[1]*(parms[3]-y[3])-y[2];
    ydot[3]=y[1]*y[2]-parms[1]*y[3];

    return ydot;
  }
}
model{
}

จากนั้นก็ใช้ expose_stan_functions เรียกจากใน R ได้เลยตามนี้ครับ

library(rstan)

expose_stan_functions("lorenz.stan")

# parms[1]=beta ,[2]=sigma,[3]=rho
parameters1<-c(2.66666,10,28)

# x=-10, y=-12, z=30.05
init.state<-c(-10,-12,30.05)

nderivs<-function(t,state,parameters){
  tmp<-lorenz(t,state,parameters)
  return(list(tmp))
}

run1<-ode(y=init.state,times=seq(0,15,0.01), 
   func=nderivs, parms=parameters1)

plot(run1)

 

 

Stan’s Program Blocks

stanblocks

หลักการของ Stan คือจะใช้ตัวแปรหรือ function อะไรก็ต้อง define ก่อน ไม่เหมือนกับ Bugs

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

MCMC ใน R

RStudio-R-JAGS

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