rjags를 이용한 베이시언Bayesian 입문 6

Introduction to Bayesian inference with JAGS 강의 정리 3

channel9.msdn.com/Events/useR-international-R-User-conferences/useR-International-R-User-2017-Conference/Introduction-to-Bayesian-inference-with-JAGS-II

Introduction to Bayesian inference with JAGS II

병원수술 완치 예제

n <- c(47, 148, 119, 810, 211, 196, 148, 215, 207, 97, 256, 360)
r <- c(0, 18, 8, 46, 8, 13, 9, 31, 14, 8, 29, 24)
N <- 12

N은 전체 병원 숫자이고

n은 N 각각의 병원에 대한 수술 횟수이고

r은 n번 수술했을때 사망한 환자수이다.

 

이제 위 데이터를 가지고 어느 병원이 가장 높은 사망률을 가졌는지 알아보자.

library(rjags)

surgical = list(N = N, n = n, r = r)

modelstring = "
model {
  for(i in 1:N) 
  {
    r[i] ~ dbin(p[i], n[i])
    p[i] ~ dbeta(1,1)
  }
}
"

m = jags.model(textConnection(modelstring), data = surgical)

update(m, 1000)
s = coda.samples(m, "p", n.iter = 1000)

summary(s)
summary(s, quantiles = c(0.05, 0.95))[[2]]
HPDinterval(s, prob = 0.90)

rjags 라이브러리를 불러오고,

n, r, N을 surgical이라는 list를 만들어 넣고

modelstring에 JAGS 모델을 작성하여 문자열로 넣는다.

JAGS에서 dbin( )은 binomial distribution을 나타낸다. 주의할점을 R에서 사용할때와 조금 다르다.

 

R에서는 y <- dbinom(x, size, prob)로 사용한다. 어떤 사건을 실행했을때 확률이 prob일때, 이를 size번 실행했을때 x번나올 probability density function값을 돌려준다.

 

JAGS에서는 x ~ dbin(prob, size)로 사용한다. Probability density function값을 돌려주지 않는다. 각각의 실행확률이 prob이고 이를 size번 실행했을때 x번 성공할 binomial distribution이라는 것을 알려준다.

 

따라서 r ~ dbin(p, n)이라 하면 사망률이 p인 수술을 n번 했을때 r번 사망할 binomial distribution임을 모델로 작성한 것이다.

 

p는 베타분포로 주었다.  shape parameter가 a=1, b=1이면 uniform distribution과 같다.

 

jags.model로 모델을 만들고

초기값의 영향을 없애주기 위하여 update(m, 1000)로 burn-in 1000번을 생략하고

coda.samples로 posterior distribution으로부터 1000개의 샘플을 취한다.

90% credible interval과 90% highest posterior density interval을 구하였다.