jags6.R

# I will now run the Bayesian analysis of the Marseilles trial using JAGS and MCMC.
# This time I use a different Bayesian model which I think is more intuitive.
# Moreover, I add simulated one new patient from each group and print their joint distribution.
#
# This time I created a more simple "slab and spike" prior distribution on the unit square.
# It's a 50-50 mixture of a uniform density on the unit square, and a uniform density
# on the main diagonal.

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

## Loading required package: coda

## Linked to JAGS 4.3.0

## Loaded modules: basemod,bugs

nA <- 2
nB <- 15
jags <- jags.model('raoult6.bug',
                    data = list('nA' = nA, 'nB' = nB), n.chains = 4, n.adapt = 100)

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

jags

## JAGS model:
## 
## model {
##  kA <- 16
##  kB <- 26
##  pAB ~ dunif(0, 1)
##  pA ~ dunif(0, 1)
##  pB ~ dunif(0, 1)
##  delta ~ dbern(0.5)
##  chi <- delta * pA + (1 - delta) * pAB
##  tau <- delta * pB + (1 - delta) * pAB
##  nA ~ dbin(chi, kA)
##  nB ~ dbin(tau, kB)
##  alpha ~ dbern(chi)
##  beta ~ dbern(tau)
##  twoAplusB <- 2 * alpha + beta
## }
## Fully observed variables:
##  nA nB

update(jags, 1000)
jags.samples(jags, c('alpha', 'beta', 'delta'), 1000)

## $alpha
## mcarray:
## [1] 0.18125
## 
## Marginalizing over: iteration(1000),chain(4) 
## 
## $beta
## mcarray:
## [1] 0.55725
## 
## Marginalizing over: iteration(1000),chain(4) 
## 
## $delta
## mcarray:
## [1] 0.9565
## 
## Marginalizing over: iteration(1000),chain(4)

samples <- coda.samples(jags, c('alpha', 'beta', 'delta'), 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
## alpha 0.1845 0.3879 0.006134       0.006570
## beta  0.5620 0.4962 0.007846       0.007833
## delta 0.9625 0.1900 0.003004       0.007861
## 
## 2. Quantiles for each variable:
## 
##       2.5% 25% 50% 75% 97.5%
## alpha    0   0   0   0     1
## beta     0   0   1   1     1
## delta    0   1   1   1     1

plot(samples)

jags.samples(jags, c('twoAplusB'), 1000)

## $twoAplusB
## mcarray:
## [1] 0.902
## 
## Marginalizing over: iteration(1000),chain(4)

samples <- coda.samples(jags, c('twoAplusB'), 1000)
summary(samples)

## 
## Iterations = 4101:5100
## 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 
##        0.92700        0.90187        0.01426        0.01426 
## 
## 2. Quantiles for each variable:
## 
##  2.5%   25%   50%   75% 97.5% 
##     0     0     1     1     3

plot(samples)