This example will go through the basics of using JAGS (https://sourceforge.net/projects/mcmc-jags/files/JAGS/3.x/) by way of the rjags library, for estimation of simple linear and generalized linear models. You must install both JAGS and rjags for this to work.
We will use the BRFSS data for the state of Texas for our example, and use BMI as a continous outcome, and obesity status outcome (BMI >= 30) as a dichotomous outcome.
First we load our data and recode some variables:
library(rjags)## Loading required package: coda
## Loading required package: lattice
## Linked to JAGS 3.3.0
## Loaded modules: basemod,bugslibrary(dplyr)##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(car)
load("~/Google Drive/dem7903_App_Hier/data/brfss_11.Rdata")
nams<-names(brfss_11)
newnames<-gsub("_", "", nams)
names(brfss_11)<-tolower(newnames)
brfss_11$statefip<-sprintf("%02d", brfss_11$state )
brfss_11$cofip<-sprintf("%03d", brfss_11$cnty )
brfss_11$cofips<-paste(brfss_11$statefip, brfss_11$cofip, sep="")
brfss_11$obese<-ifelse(brfss_11$bmi5/100 >=30, 1,0)
brfss_11$black<-recode(brfss_11$racegr2, recodes="2=1; 9=NA; else=0", as.factor.result=F)
brfss_11$white<-recode(brfss_11$racegr2, recodes="1=1; 9=NA; else=0", as.factor.result=F)
brfss_11$other<-recode(brfss_11$racegr2, recodes="3:4=1; 9=NA; else=0", as.factor.result=F)
brfss_11$hispanic<-recode(brfss_11$racegr2, recodes="5=1; 9=NA; else=0", as.factor.result=F)
#education level
brfss_11$lths<-recode(brfss_11$educa, recodes="1:3=1;9=NA; else=0", as.factor.result=F)
brfss_11$coll<-recode(brfss_11$educa, recodes="5:6=1;9=NA; else=0", as.factor.result=F)
brfss_11$agez<-as.numeric(scale(brfss_11$age, center=T, scale=T))
brfss_11$lowinc<-recode(brfss_11$incomg, recodes = "1:3=1; 4:5=0; else=NA")Next, I use tselect the observations from Texas that are not missing on some key variables
tx<-brfss_11[brfss_11$statefip=="48"& is.na(brfss_11$obese)==F& is.na(brfss_11$black)==F& is.na(brfss_11$lths)==F& is.na(brfss_11$lowinc)==F ,]
nwncos<-table(tx$cofips)
nwncos #Number of people within counties##
## 48029 48133 48157 48201 48303 48329 48423 48439 48453
## 863 498 772 1185 603 421 452 483 859tx$conum<-rep(1:length(unique(tx$cofips)), nwncos[nwncos!=0])
length(unique(tx$conum)) #Number of counties## [1] 9Here is a linear mixed model for bmi. This model includes a random intercept and a few predictor variables There a loads of ways to do this, but I like doing it this way. It uses what is known as “nested indexing” in the BUGS language This is basically a way of explicitly nesting individuals within groups. I write my code as a big string, then feed it to jags.
model1<-"
model{
#Likelihood
for( i in 1:n)
{
bmi[i]~dnorm(mu[i], tau)
mu[i]<-b0+b[1]*black[i]+b[2]*hisp[i]+b[3]*other[i]+b[4]*lths[i]+b[5]*coll[i]+u[cofips[i]]
}
for (j in 1:ncos)
{
u[j]~dnorm(0, tau_u)
}
#priors
b0~dnorm(0, .01)
for(j in 1:5) { b[j]~dnorm(0, .01)}
tau<-pow(sd, -2)
sd~dunif(0,100)
tau_u<-pow(sd_u, -2)
sd_u~dunif(0,100)
}
"Next, we have to make a data list for jags, which contains anything we are reading into jags as data. I z-score the bmi variable prior to putting into the model
dat<-list(bmi=as.numeric(scale(tx$bmi5/100, center=T, scale=T)), obese=tx$obese, black=tx$black, hisp=tx$hispanic, other=tx$other, lths=tx$lths, coll=tx$coll, age=tx$agez, lowinc=tx$lowinc, n=length(tx$obese),cofips=tx$conum, ncos=length(unique(tx$cofips)))
#quick summary
lapply(dat, summary)## $bmi
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.510 -0.680 -0.182 0.000 0.508 9.030
##
## $obese
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.283 1.000 1.000
##
## $black
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.0942 0.0000 1.0000
##
## $hisp
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.169 0.000 1.000
##
## $other
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.0582 0.0000 1.0000
##
## $lths
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.0922 0.0000 1.0000
##
## $coll
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 0.691 1.000 1.000
##
## $age
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.5100 -0.6170 0.1620 0.0903 0.7740 2.5000
##
## $lowinc
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.376 1.000 1.000
##
## $n
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6140 6140 6140 6140 6140 6140
##
## $cofips
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 3.00 4.00 4.76 7.00 9.00
##
## $ncos
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 9 9 9 9 9 9To use jags, we have to create a jags.model object, which contains the text representation of our model, our data, and some other parameters for the MCMC run
init.rng1<-list(".RNG.seed" = 1234, ".RNG.name" = "base::Mersenne-Twister")
init.rng2<-list(".RNG.seed" = 5678, ".RNG.name" = "base::Mersenne-Twister")
mod<-jags.model(file=textConnection(model1), data=dat,inits =list(init.rng1, init.rng2) , n.chains=2)## Warning: Unused variable "obese" in data
## Warning: Unused variable "age" in data
## Warning: Unused variable "lowinc" in data## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 43094
##
## Initializing model#next, we update the model, this is the "burn in" period
update(mod, 1000)Next, we examine a few other elements of the model, including the posterior densities of the parameters, and First, we must collect some samples of each parameter using the coda.samples() function.
#burn in for 10,000 iterations
update(mod, 10000)
#sample 2000 samples from each chain (5,000/5 = 1000 * 2 chains = 2000)
samps<-coda.samples(mod, variable.names=c("b0", "b", "sd", "sd_u"), n.iter=5000, n.thin=5)
#Numerical summary of each parameter, here I also include the 90% credible interval:
summary(samps, quantiles = c(.025, .05, .95, .975))##
## Iterations = 12001:17000
## Thinning interval = 1
## Number of chains = 2
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## b[1] 0.4971 0.04417 4.42e-04 0.000552
## b[2] 0.1742 0.03720 3.72e-04 0.000571
## b[3] -0.1563 0.05453 5.45e-04 0.000661
## b[4] 0.1100 0.05096 5.10e-04 0.000900
## b[5] -0.0672 0.03138 3.14e-04 0.000887
## b0 -0.0277 0.03552 3.55e-04 0.001214
## sd 0.9836 0.00902 9.02e-05 0.000109
## sd_u 0.0563 0.02542 2.54e-04 0.000688
##
## 2. Quantiles for each variable:
##
## 2.5% 5% 95% 97.5%
## b[1] 0.4093 0.4243 0.5690 0.5845
## b[2] 0.1016 0.1128 0.2356 0.2470
## b[3] -0.2635 -0.2463 -0.0673 -0.0505
## b[4] 0.0105 0.0261 0.1936 0.2097
## b[5] -0.1271 -0.1182 -0.0148 -0.0048
## b0 -0.0975 -0.0863 0.0311 0.0424
## sd 0.9664 0.9689 0.9986 1.0014
## sd_u 0.0165 0.0217 0.1036 0.1173The “effective sample size” tells us how many independent samples we have, out of the 1000* 2 chains (maximum = 2000 here). If this number is low, then we have a lot of autocorrelation in the chains
effectiveSize(samps)## b[1] b[2] b[3] b[4] b[5] b0 sd sd_u
## 6415.1 4247.9 6832.1 3211.1 1252.8 855.3 6806.2 1397.1#Examine convergence of the Markov chains using the Gelman-Brooks-Rubin diagnostic
gelman.diag(samps)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## b[1] 1 1.00
## b[2] 1 1.01
## b[3] 1 1.00
## b[4] 1 1.00
## b[5] 1 1.02
## b0 1 1.01
## sd 1 1.00
## sd_u 1 1.00
##
## Multivariate psrf
##
## 1#here's a way to get p-values for each of the beta's from the model using the samples:
sampmat<-as.matrix(samps)
str(sampmat)## num [1:10000, 1:8] 0.479 0.544 0.449 0.363 0.534 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:8] "b[1]" "b[2]" "b[3]" "b[4]" ...head(sampmat)## b[1] b[2] b[3] b[4] b[5] b0 sd sd_u
## [1,] 0.4792 0.1900 -0.10163 0.07904 -0.07368 -0.021300 0.9842 0.03113
## [2,] 0.5441 0.1926 -0.19263 0.08153 -0.08015 -0.003315 0.9716 0.03409
## [3,] 0.4490 0.1223 -0.12613 0.10205 -0.07990 0.027161 0.9737 0.03913
## [4,] 0.3627 0.1395 -0.09562 0.07516 -0.08660 0.043738 0.9665 0.09570
## [5,] 0.5344 0.1467 -0.22094 0.09766 -0.10558 0.040411 0.9697 0.07939
## [6,] 0.5515 0.1588 -0.11270 0.06226 -0.10785 0.027470 0.9958 0.07199apply(sampmat[, 1:5], 2, function(x) mean(x > 0)) # Get p(beta > 0)## b[1] b[2] b[3] b[4] b[5]
## 1.0000 1.0000 0.0021 0.9851 0.0183apply(sampmat[, 1:5], 2, function(x) mean(x < 0)) # Get p(beta < 0)## b[1] b[2] b[3] b[4] b[5]
## 0.0000 0.0000 0.9979 0.0149 0.9817Now, we expand the model to include random slopes and intercepts. We could just sample the random coefficients from independent normal distributions, but we should sample them from a Multivariate Normal. This takes a little more coding in order to make the prior on the covariance matrix of the MN Normal, but also to transform it back into the variance scale (vs. precision), and calculate the correlation parameter.
model2<-"
model{
#Likelihood
for( i in 1:n)
{
bmi[i]~dnorm(mu[i], tau)
mu[i]<-b[1]*black[i]+b[2]*hisp[i]+b[3]*other[i]+b[4]*lths[i]+b[5]*coll[i]+u[cofips[i], 1]*lowinc[i]+u[cofips[i], 2]
}
for (j in 1:ncos)
{
u[j, 1:2] ~ dmnorm(meanu[], prec.Sigma[,])
}
#priors
#NOTICE I removed the separate intercept, b0
for(j in 1:5) { b[j]~dnorm(0, .01)}
meanu[1]~dnorm(0, .001)
meanu[2]~dnorm(0, .001)
prec.Sigma[1:2, 1:2] ~ dwish(Omega[,], 2) #Wishart is MV form of gamma on precision, it will be square, with nrow = ncol = # of correlated random effects, so if we have another random slope, it would be prec.Sigma[1:3, 1:3] ~ dwish(Omega[,], 3)
Sigma[1:2, 1:2]<-inverse(prec.Sigma[,]) #get the covariance matrix on the variance scale
rho12<-Sigma[1,2]/ sqrt(Sigma[1,1]* Sigma[2,2])
#Set some initial values for the covariance matrix
for (j in 1:2){ for (k in 1:2){ Omega[j,k] <-equals(j,k)*.1 } }
tau<-pow(sd, -2)
sd~dunif(0,100)
}
"
mod2<-jags.model(file=textConnection(model2), data=dat, n.chains=2,inits =list(init.rng1, init.rng2) )## Warning: Unused variable "obese" in data
## Warning: Unused variable "age" in data## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 49371
##
## Initializing modelupdate(mod2, 10000)
#collect 2000 samples of the parameters
samps2<-coda.samples(mod2, variable.names=c( "b", "Sigma", "rho12", "sd", "u"), n.iter=5000, n.thin=5)
effectiveSize(samps2)## Sigma[1,1] Sigma[2,1] Sigma[1,2] Sigma[2,2] b[1] b[2]
## 6387 6399 6399 7388 6136 4176
## b[3] b[4] b[5] rho12 sd u[1,1]
## 7031 3581 1001 6305 6574 6802
## u[2,1] u[3,1] u[4,1] u[5,1] u[6,1] u[7,1]
## 6847 8145 7963 7451 8079 7644
## u[8,1] u[9,1] u[1,2] u[2,2] u[3,2] u[4,2]
## 7699 6994 2137 4106 1792 1833
## u[5,2] u[6,2] u[7,2] u[8,2] u[9,2]
## 2587 3335 3275 2532 1800#Numerical summary of each parameter:
summary(samps2, quantiles = c(.025, .05, .95, .975))##
## Iterations = 11001:16000
## Thinning interval = 1
## Number of chains = 2
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## Sigma[1,1] 0.02671 0.01774 1.77e-04 0.000222
## Sigma[2,1] -0.00717 0.01149 1.15e-04 0.000144
## Sigma[1,2] -0.00717 0.01149 1.15e-04 0.000144
## Sigma[2,2] 0.02262 0.01421 1.42e-04 0.000165
## b[1] 0.47758 0.04510 4.51e-04 0.000577
## b[2] 0.15330 0.03791 3.79e-04 0.000587
## b[3] -0.15309 0.05513 5.51e-04 0.000658
## b[4] 0.09529 0.05100 5.10e-04 0.000857
## b[5] -0.04030 0.03242 3.24e-04 0.001024
## rho12 -0.26456 0.31746 3.17e-03 0.003997
## sd 0.98244 0.00879 8.79e-05 0.000108
## u[1,1] 0.10737 0.06319 6.32e-04 0.000773
## u[2,1] -0.08537 0.08016 8.02e-04 0.000969
## u[3,1] 0.06973 0.07953 7.95e-04 0.000881
## u[4,1] 0.10778 0.05619 5.62e-04 0.000630
## u[5,1] 0.09490 0.07219 7.22e-04 0.000838
## u[6,1] 0.04877 0.08627 8.63e-04 0.000960
## u[7,1] 0.14258 0.07959 7.96e-04 0.000910
## u[8,1] -0.00270 0.08057 8.06e-04 0.000918
## u[9,1] 0.15818 0.06879 6.88e-04 0.000824
## u[1,2] -0.02597 0.05075 5.07e-04 0.001101
## u[2,2] 0.04501 0.06715 6.72e-04 0.001050
## u[3,2] -0.13925 0.04847 4.85e-04 0.001149
## u[4,2] -0.08121 0.04542 4.54e-04 0.001062
## u[5,2] -0.03382 0.05494 5.49e-04 0.001080
## u[6,2] -0.02231 0.05898 5.90e-04 0.001021
## u[7,2] -0.14203 0.06242 6.24e-04 0.001093
## u[8,2] 0.00133 0.05933 5.93e-04 0.001186
## u[9,2] -0.16325 0.04857 4.86e-04 0.001147
##
## 2. Quantiles for each variable:
##
## 2.5% 5% 95% 97.5%
## Sigma[1,1] 0.00878 0.01002 0.05852 0.07237
## Sigma[2,1] -0.03473 -0.02679 0.00656 0.01001
## Sigma[1,2] -0.03473 -0.02679 0.00656 0.01001
## Sigma[2,2] 0.00790 0.00893 0.04843 0.05942
## b[1] 0.38936 0.40259 0.55108 0.56424
## b[2] 0.07910 0.09077 0.21501 0.22830
## b[3] -0.26073 -0.24253 -0.06050 -0.04654
## b[4] -0.00421 0.01364 0.18050 0.19578
## b[5] -0.10328 -0.09313 0.01272 0.02270
## rho12 -0.77977 -0.72819 0.30988 0.42641
## sd 0.96529 0.96819 0.99707 0.99993
## u[1,1] -0.01653 0.00400 0.21255 0.23246
## u[2,1] -0.24647 -0.21895 0.04394 0.06582
## u[3,1] -0.08645 -0.06175 0.19881 0.22783
## u[4,1] -0.00109 0.01571 0.19976 0.21920
## u[5,1] -0.04635 -0.02412 0.21328 0.23618
## u[6,1] -0.11938 -0.09208 0.19095 0.21932
## u[7,1] -0.01155 0.01238 0.27411 0.30150
## u[8,1] -0.16260 -0.13561 0.12774 0.15185
## u[9,1] 0.02428 0.04469 0.27291 0.29363
## u[1,2] -0.12764 -0.10931 0.05870 0.07350
## u[2,2] -0.08499 -0.06414 0.15744 0.17856
## u[3,2] -0.23329 -0.21833 -0.05993 -0.04389
## u[4,2] -0.17033 -0.15719 -0.00698 0.00755
## u[5,2] -0.14097 -0.12505 0.05665 0.07226
## u[6,2] -0.13598 -0.11827 0.07603 0.09674
## u[7,2] -0.26444 -0.24454 -0.04010 -0.02033
## u[8,2] -0.11405 -0.09544 0.09949 0.12065
## u[9,2] -0.25653 -0.24242 -0.08274 -0.06650#Examine convergence of the Markov chains using the Gelman-Brooks-Rubin diagnostic
gelman.diag(samps2, multivariate = F)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## Sigma[1,1] 1 1.01
## Sigma[2,1] 1 1.01
## Sigma[1,2] 1 1.01
## Sigma[2,2] 1 1.00
## b[1] 1 1.00
## b[2] 1 1.00
## b[3] 1 1.00
## b[4] 1 1.00
## b[5] 1 1.00
## rho12 1 1.00
## sd 1 1.00
## u[1,1] 1 1.00
## u[2,1] 1 1.00
## u[3,1] 1 1.00
## u[4,1] 1 1.00
## u[5,1] 1 1.00
## u[6,1] 1 1.00
## u[7,1] 1 1.00
## u[8,1] 1 1.00
## u[9,1] 1 1.00
## u[1,2] 1 1.00
## u[2,2] 1 1.00
## u[3,2] 1 1.00
## u[4,2] 1 1.00
## u[5,2] 1 1.00
## u[6,2] 1 1.00
## u[7,2] 1 1.00
## u[8,2] 1 1.00
## u[9,2] 1 1.00So, again, our model looks good after burning in for 10,000 iterations
The densities don’t look too out of wack
our traceplots reveal good mixing in the chains
Our effective sample sizes are all pretty large, (i.e. the chains are mixing well and providing independent samples at each iteration)
The Gelman-Brooks-Rubin diagnostics show that numerically, there is little to no variation between the chains
In Bayesian models, there isn’t a direct equivalent of a likelihood ratio test, and traditionally people use the Deviance Information Criteria (DIC) as a measure of relative model fit (Spiegelhalter et al, 2002). You can think of this as a Bayesian equivalent of the AIC, if you’re used to working with that. It’s a measure of model information (deviance), penalized for complexity (# of parameters).
dic1<-dic.samples(mod, n.iter = 1000, type = "pD")
dic2<-dic.samples(mod2, n.iter = 1000, type = "pD")
dic1## Mean deviance: 17208
## penalty 12.4
## Penalized deviance: 17221dic2## Mean deviance: 17195
## penalty 21.6
## Penalized deviance: 17217We see the penalized deviance in model 2 (the random slopes + intercepts model) is only lower than the model with only random intercepts. While there is no hard and fast rule of thumb for how big a difference there has to be, in their orginal paper, they suggest using differences in DIC greater than 7 to indicate that one model is preferred to another.
Here, we fit the hierarchical logistic regression model, here I just fit the random slopes and
model3<-"
model{
#Likelihood
for( i in 1:n)
{
obese[i]~dbern(p[i])
logit(p[i])<-b[1]*black[i]+b[2]*hisp[i]+b[3]*other[i]+b[4]*lths[i]+b[5]*coll[i]+u[cofips[i], 1]*lowinc[i]+u[cofips[i], 2]
}
for (j in 1:ncos)
{
u[j, 1:2] ~ dmnorm(meanu[], prec.Sigma[,])
}
#priors
#NOTICE I removed the separate intercept, b0
for(j in 1:5) { b[j]~dnorm(0, .01)}
meanu[1]~dnorm(0, .001)
meanu[2]~dnorm(0, .001)
prec.Sigma[1:2, 1:2] ~ dwish(Omega[,], 2) #Wishart is MV form of gamma on precision, it will be square, with nrow = ncol = # of correlated random effects, so if we have another random slope, it would be prec.Sigma[1:3, 1:3] ~ dwish(Omega[,], 3)
Sigma[1:2, 1:2]<-inverse(prec.Sigma[,]) #get the covariance matrix on the variance scale
rho12<-Sigma[1,2]/ sqrt(Sigma[1,1]* Sigma[2,2])
#Set some initial values for the covariance matrix
for (j in 1:2){ for (k in 1:2){ Omega[j,k] <-equals(j,k)*.1 } }
}
"And now we fit the model:
load.module("glm")## module glm loadedmod3<-jags.model(file=textConnection(model3), data=dat, n.chains=2, inits =list(init.rng1, init.rng2) )## Warning: Unused variable "bmi" in data
## Warning: Unused variable "age" in data## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 49560
##
## Initializing model#GLMM updating can take longer the LMM's, so make yourself a drink
update(mod3, 10000)
#collect 2000 samples of the parameters
samps3<-coda.samples(mod3, variable.names=c( "b", "Sigma", "rho12", "sd", "u"), n.iter=5000, n.thin=5)## Warning: Failed to set trace monitor for sd
## Node not foundeffectiveSize(samps3)## Sigma[1,1] Sigma[2,1] Sigma[1,2] Sigma[2,2] b[1] b[2]
## 3509 3268 3268 3840 6289 5837
## b[3] b[4] b[5] rho12 u[1,1] u[2,1]
## 4871 5209 2769 2827 4739 2897
## u[3,1] u[4,1] u[5,1] u[6,1] u[7,1] u[8,1]
## 4503 4435 4486 4438 3912 5087
## u[9,1] u[1,2] u[2,2] u[3,2] u[4,2] u[5,2]
## 3861 3200 2721 3317 3208 3209
## u[6,2] u[7,2] u[8,2] u[9,2]
## 3423 2810 3636 2889#Numerical summary of each parameter:
summary(samps3, quantiles = c(.025, .05, .95, .975))##
## Iterations = 10001:15000
## Thinning interval = 1
## Number of chains = 2
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## Sigma[1,1] 0.05092 0.0420 0.000420 0.000710
## Sigma[2,1] -0.02450 0.0318 0.000318 0.000556
## Sigma[1,2] -0.02450 0.0318 0.000318 0.000556
## Sigma[2,2] 0.04893 0.0361 0.000361 0.000586
## b[1] 0.78988 0.0950 0.000950 0.001197
## b[2] 0.25223 0.0829 0.000829 0.001085
## b[3] -0.24772 0.1384 0.001384 0.001983
## b[4] 0.15362 0.1090 0.001090 0.001515
## b[5] -0.12526 0.0735 0.000735 0.001418
## rho12 -0.42170 0.3314 0.003314 0.006272
## u[1,1] 0.29198 0.1251 0.001251 0.001816
## u[2,1] -0.00202 0.1714 0.001714 0.003193
## u[3,1] 0.23332 0.1525 0.001525 0.002273
## u[4,1] 0.22194 0.1131 0.001131 0.001712
## u[5,1] 0.17300 0.1409 0.001409 0.002104
## u[6,1] 0.26902 0.1594 0.001594 0.002399
## u[7,1] 0.37123 0.1633 0.001633 0.002613
## u[8,1] 0.24381 0.1492 0.001492 0.002091
## u[9,1] 0.36945 0.1439 0.001439 0.002317
## u[1,2] -0.95162 0.1111 0.001111 0.001982
## u[2,2] -0.81915 0.1485 0.001485 0.002902
## u[3,2] -1.18337 0.1098 0.001098 0.001906
## u[4,2] -1.07138 0.1030 0.001030 0.001829
## u[5,2] -1.01012 0.1202 0.001202 0.002133
## u[6,2] -1.02069 0.1265 0.001265 0.002201
## u[7,2] -1.20257 0.1410 0.001410 0.002662
## u[8,2] -1.06812 0.1262 0.001262 0.002094
## u[9,2] -1.28017 0.1148 0.001148 0.002135
##
## 2. Quantiles for each variable:
##
## 2.5% 5% 95% 97.5%
## Sigma[1,1] 0.011976 0.014230 0.12589 0.1582
## Sigma[2,1] -0.106460 -0.080552 0.00616 0.0111
## Sigma[1,2] -0.106460 -0.080552 0.00616 0.0111
## Sigma[2,2] 0.013283 0.015542 0.11453 0.1429
## b[1] 0.602134 0.632823 0.94534 0.9748
## b[2] 0.092640 0.117423 0.38913 0.4151
## b[3] -0.524596 -0.476269 -0.02158 0.0188
## b[4] -0.062017 -0.025978 0.33193 0.3666
## b[5] -0.268153 -0.244717 -0.00433 0.0202
## rho12 -0.884572 -0.851601 0.21651 0.3447
## u[1,1] 0.046049 0.086269 0.49796 0.5387
## u[2,1] -0.365504 -0.299176 0.26323 0.3091
## u[3,1] -0.075701 -0.020083 0.48049 0.5271
## u[4,1] -0.000197 0.035877 0.40814 0.4471
## u[5,1] -0.109276 -0.061293 0.40244 0.4464
## u[6,1] -0.041505 0.008495 0.53481 0.5878
## u[7,1] 0.071796 0.113991 0.65269 0.7171
## u[8,1] -0.048158 -0.000546 0.49018 0.5404
## u[9,1] 0.099448 0.140274 0.61445 0.6654
## u[1,2] -1.170413 -1.132977 -0.76920 -0.7337
## u[2,2] -1.098305 -1.056944 -0.56741 -0.5206
## u[3,2] -1.401253 -1.365394 -1.00308 -0.9722
## u[4,2] -1.275440 -1.244081 -0.90193 -0.8707
## u[5,2] -1.240625 -1.206127 -0.81338 -0.7729
## u[6,2] -1.273956 -1.229580 -0.81422 -0.7772
## u[7,2] -1.490931 -1.441643 -0.97886 -0.9395
## u[8,2] -1.316369 -1.276431 -0.86387 -0.8250
## u[9,2] -1.512763 -1.473872 -1.09469 -1.0640#Examine convergence of the Markov chains using the Gelman-Brooks-Rubin diagnostic
gelman.diag(samps3, multivariate = F)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## Sigma[1,1] 1 1.00
## Sigma[2,1] 1 1.00
## Sigma[1,2] 1 1.00
## Sigma[2,2] 1 1.00
## b[1] 1 1.01
## b[2] 1 1.00
## b[3] 1 1.00
## b[4] 1 1.00
## b[5] 1 1.00
## rho12 1 1.00
## u[1,1] 1 1.00
## u[2,1] 1 1.00
## u[3,1] 1 1.00
## u[4,1] 1 1.00
## u[5,1] 1 1.00
## u[6,1] 1 1.00
## u[7,1] 1 1.00
## u[8,1] 1 1.00
## u[9,1] 1 1.00
## u[1,2] 1 1.00
## u[2,2] 1 1.00
## u[3,2] 1 1.00
## u[4,2] 1 1.00
## u[5,2] 1 1.00
## u[6,2] 1 1.00
## u[7,2] 1 1.00
## u[8,2] 1 1.00
## u[9,2] 1 1.00The model looks converged, other diagnostics are all good. The dic is:
dic.samples(mod3, n.iter=1000)## Mean deviance: 7145
## penalty 17.8
## Penalized deviance: 7163