if (FALSE) {
library(jarbes)
# Example: Stemcells
data("stemcells")
stemcells$TE = stemcells$effect.size
stemcells$seTE = stemcells$se.effect
# Beta(0.5, 1)
a.0 = 0.5
a.1 = 1
# alpha.max
N = dim(stemcells)[1]
alpha.max = 1/5 *((N-1)*a.0 - a.1)/(a.0 + a.1)
alpha.max
# K.max
K.max = 1 + 5*alpha.max
K.max = round(K.max)
K.max
set.seed(20233)
bcmix.2.stemcell = bcmixmeta(stemcells,
mean.mu.0=0, sd.mu.0=100,
B.lower = -15,
B.upper = 15,
alpha.0 = 0.5,
alpha.1 = alpha.max,
a.0 = a.0,
a.1 = a.1,
K = K.max,
sort.priors = FALSE,
df.scale.between = 1,
scale.sigma.between = 0.5,
nr.chains = 4,
nr.iterations = 50000,
nr.adapt = 1000,
nr.burnin = 10000,
nr.thin = 4)
diagnostic(bcmix.2.stemcell, y.lim = c(-1, 15), title.plot = "Default priors")
bcmix.2.stemcell.mcmc <- as.mcmc(bcmix.1.stemcell$BUGSoutput$sims.matrix)
theta.names <- paste(paste("theta[",1:31, sep=""),"]", sep="")
theta.b.names <- paste(paste("theta.bias[",1:31, sep=""),"]", sep="")
theta.b.greek.names <- paste(paste("theta[",1:31, sep=""),"]^B", sep="")
theta.greek.names <- paste(paste("theta[",1:31, sep=""),"]", sep="")
caterplot(bcmix.2.stemcell.mcmc,
parms = theta.names, # theta
labels = theta.greek.names,
greek = T,
labels.loc="axis", cex =0.7,
col = "black",
style = "plain",
reorder = F,
val.lim =c(-6, 16),
quantiles = list(outer=c(0.05,0.95),inner=c(0.16,0.84)),
x.lab = "Effect: mean difference"
)
title( "95% posterior intervals of studies' effects")
caterplot(bcmix.2.stemcell.mcmc,
parms = theta.b.names, # theta.bias
labels = theta.greek.names,
greek = T,
labels.loc="no",
cex = 0.7,
col = "grey",
style = "plain", reorder = F,
val.lim =c(-6, 16),
quantiles = list(outer=c(0.025,0.975),inner=c(0.16,0.84)),
add = TRUE,
collapse=TRUE, cat.shift= -0.5,
)
attach.jags(bcmix.2.stemcell, overwrite = TRUE)
abline(v=mean(mu.0), lwd =2, lty =2)
legend(9, 20, legend = c("bias corrected", "biased"),
lty = c(1,1), lwd = c(2,2), col = c("black", "grey"))
}
Run the code above in your browser using DataLab