Loading the necessary libraries

library(rjags)
library(mcmcplots)

Binomial Proportion

modelString = 
  "model { 
  
    x~dbin(theta,n)        #  Likelihood function

    theta ~ dbeta(1,1)     #  Prior density for theta
    
  }"

#Write the model to a file
writeLines(modelString,con="model.txt")

#Compiling the model together with the data
jagsModel = jags.model("model.txt",data=list(x=6,n=20))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 4
## 
## Initializing model
#Burn-in
update(jagsModel,n.iter=5000)

# The parameter(s) to be monitored.
parameters = c( "theta")  

#Sampling from the posterior distribution:
output = coda.samples(jagsModel,variable.names=parameters,n.iter=10000)

# Examining the posterior sample
# Density plots
denplot(output, parms=c("theta"))

# History plot(s)
traplot(output, parms=c("theta"))

# Summary statistics 
summary(output) 
## 
## Iterations = 6001:16000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##      0.3183555      0.0964720      0.0009647      0.0012358 
## 
## 2. Quantiles for each variable:
## 
##   2.5%    25%    50%    75%  97.5% 
## 0.1435 0.2483 0.3141 0.3835 0.5166

Binomial Proportion Difference

modelString = 
  "model { x  ~ dbin(theta1, n1)   #  Likelihood for group 1
         theta1 ~ dbeta(1,1)         #  Prior for theta1

         y  ~ dbin(theta2, n2)   #  Likelihood for group 2
         theta2 ~ dbeta(1,1)         #  Prior for theta2

      propdiff <- theta1-theta2        #  Calculate difference for binomial parameters
      rr <- theta1/theta2              #  Calculate relative risk
                                       #  Calculate odds ratio
      or<- theta1*(1-theta2)/((1-theta1)*theta2)
  }"

#Write the model to a file
writeLines(modelString,con="model.txt")

#Compiling the model together with the data
jagsModel = jags.model("model.txt",data=list(x  = 6,
                                             n1 = 20,
                                             y  = 20,
                                             n2 = 25))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 2
##    Total graph size: 14
## 
## Initializing model
#Burn-in
update( jagsModel,n.iter=5000)

# The parameter(s) to be monitored.
parameters = c( "theta1", "theta2", "propdiff", "or", "rr")  

#Sampling from the posterior distribution:
output = coda.samples(jagsModel,variable.names=parameters,n.iter=10000)

# Examining the posterior sample
# History plots
traplot(output, parms=c("propdiff"))

# Density plots
denplot(output, parms=c("propdiff"))