# 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 their joint distribution.
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('raoult3.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: 36
##
## Initializing model
jags
## JAGS model:
##
## model {
## kA <- 16
## kB <- 26
## lambdaAB ~ dnorm(0, 1)
## lambdaA ~ dnorm(0, 1)
## lambdaB ~ dnorm(0, 1)
## delta ~ dbern(0.5)
## pA <- exp(lambdaA) / (1 + exp(lambdaA))
## pB <- exp(lambdaB) / (1 + exp(lambdaB))
## pAB <- exp(lambdaAB) / (1 + exp(lambdaAB))
## 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)
## AplusBplus <- (alpha == 1) * (beta == 1)
## AplusBmin <- (alpha == 1) * (beta == 0)
## AminBplus <- (alpha == 0) * (beta == 1)
## AminBmin <- (alpha == 0) * (beta == 0)
## }
## Fully observed variables:
## nA nB
update(jags, 1000)
jags.samples(jags, c('alpha', 'beta', 'delta'), 1000)
## $alpha
## mcarray:
## [1] 0.22125
##
## Marginalizing over: iteration(1000),chain(4)
##
## $beta
## mcarray:
## [1] 0.571
##
## Marginalizing over: iteration(1000),chain(4)
##
## $delta
## mcarray:
## [1] 0.95025
##
## 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.2233 0.4165 0.006585 0.006984
## beta 0.5605 0.4964 0.007849 0.007997
## delta 0.9577 0.2012 0.003181 0.007934
##
## 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('AplusBplus', 'AplusBmin', 'AminBplus', 'AminBmin'), 1000)
## $AminBmin
## mcarray:
## [1] 0.351
##
## Marginalizing over: iteration(1000),chain(4)
##
## $AminBplus
## mcarray:
## [1] 0.42675
##
## Marginalizing over: iteration(1000),chain(4)
##
## $AplusBmin
## mcarray:
## [1] 0.10225
##
## Marginalizing over: iteration(1000),chain(4)
##
## $AplusBplus
## mcarray:
## [1] 0.12
##
## Marginalizing over: iteration(1000),chain(4)
samples <- coda.samples(jags, c('AplusBplus', 'AplusBmin', 'AminBplus', 'AminBmin'), 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
## AminBmin 0.3520 0.4777 0.007552 0.007691
## AminBplus 0.4255 0.4945 0.007818 0.007818
## AplusBmin 0.1022 0.3030 0.004791 0.004864
## AplusBplus 0.1202 0.3253 0.005143 0.005133
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## AminBmin 0 0 0 1 1
## AminBplus 0 0 0 1 1
## AplusBmin 0 0 0 0 1
## AplusBplus 0 0 0 0 1