jags.R

# I will now run the Bayesian analysis of the Marseilles trial using JAGS and MCMC

setwd("/Users/richard")
library('rjags')

## Loading required package: coda

## Linked to JAGS 4.3.0

## Loaded modules: basemod,bugs

library('coda')
X <- 2
Y <- 15
jags <- jags.model('raoult.bug',
                    data = list('X' = X, 'Y' = Y), n.chains = 4, n.adapt = 100)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 3
##    Total graph size: 21
## 
## Initializing model

jags

## JAGS model:
## 
## model {
##  m <- 16
##  n <- 26
##  psi ~ dnorm(0, 1)
##  mu ~ dnorm(0, 1)
##  D ~ dbern(0.5)
##  Xodds <- exp(mu + D * psi / 2)
##  Yodds <- exp(mu - D * psi / 2)
##  Xprob <- Xodds / (1 + Xodds)
##  Yprob <- Yodds / (1 + Yodds)
##  X ~ dbin(Xprob, m)
##  Y ~ dbin(Yprob, n)
## }
## Fully observed variables:
##  X Y

update(jags, 1000)
jags.samples(jags, c('mu', 'psi', 'D'), 1000)

## $D
## mcarray:
## [1] 0.9225
## 
## Marginalizing over: iteration(1000),chain(4) 
## 
## $mu
## mcarray:
## [1] -0.5557471
## 
## Marginalizing over: iteration(1000),chain(4) 
## 
## $psi
## mcarray:
## [1] -1.31069
## 
## Marginalizing over: iteration(1000),chain(4)

samples <- coda.samples(jags, c('mu', 'psi', 'D'), 1000)
summary(samples)

## 
## Iterations = 2101:3100
## Thinning interval = 1 
## Number of chains = 4 
## Sample size per chain = 1000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##        Mean     SD Naive SE Time-series SE
## D    0.9230 0.2666 0.004216       0.007330
## mu  -0.5557 0.3414 0.005399       0.008094
## psi -1.3100 0.7490 0.011843       0.021629
## 
## 2. Quantiles for each variable:
## 
##       2.5%     25%     50%     75%   97.5%
## D    0.000  1.0000  1.0000  1.0000 1.00000
## mu  -1.261 -0.7809 -0.5409 -0.3266 0.08586
## psi -2.639 -1.7863 -1.3428 -0.9172 0.45581

plot(samples)

samples2 <- coda.samples(jags, c('Xprob', 'Yprob'), 1000)
summary(samples2)

## 
## Iterations = 3101:4100
## Thinning interval = 1 
## Number of chains = 4 
## Sample size per chain = 1000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean     SD Naive SE Time-series SE
## Xprob 0.2468 0.1006 0.001590       0.002840
## Yprob 0.5240 0.0942 0.001489       0.001844
## 
## 2. Quantiles for each variable:
## 
##          2.5%    25%    50%    75%  97.5%
## Xprob 0.08866 0.1750 0.2312 0.3053 0.4740
## Yprob 0.33974 0.4613 0.5238 0.5900 0.7072

plot(samples2)