library(rjags)
s1=22;s2=10;n1=158;n2=78
c(s1/n1,s2/n2,s1/n1-s2/n2)
shape1=shape2=0.5
data=list(s1=s1,s2=s2,n1=n1,n2=n2,shape1=shape1,shape2=shape2)
modelString='
model{
# data
s1 ~ dbin(p1,n1)
s2 ~ dbin(p2,n2)
# model
p1 ~ dbeta(shape1,shape2)
p2 ~ dbeta(shape1,shape2)
delta = p1 - p2
# priors
pri.delta= pri.p1 - pri.p2
pri.p1 ~ dbeta(shape1,shape2)
pri.p2 ~ dbeta(shape1,shape2)
}
'
writeLines(modelString,'tmp.txt')
n.iter=1e4 ; thin=1
jagsModel=jags.model('tmp.txt',data=data,n.chains=4,n.adapt=n.iter/5)
update(jagsModel)
codaSamples=coda.samples(jagsModel,c('delta','pri.delta'),n.iter,thin)
plot(codaSamples)
js=as.data.frame(as.matrix(codaSamples))
BEST::plotPost(js$delta,compVal = 0)