###################################################################################### ###################################################################################### # Albert I. and Parent E.. - TD1: Initiation Stat Bayésienne # 2012 # R program for the basic Beta-Binomial example # applied to listeria in milk # Chapter 7 ###################################################################################### ###################################################################################### rm(list=ls()) library("R2WinBUGS") library(MASS) bugs.directory="c:/Program Files/WinBUGS14/" bugs.workingdir="C:/BookBiobayes" setwd(bugs.workingdir) ######################################################################### # Run basic Beta-Binomial model ######################################################################### ############################################################### # Define model model <- "model_Postulat1.txt" sink(model) cat(" model { for( i in 1 : N ) { r[i] ~ dbin(p,n[i])} p ~ dbeta(1,1) } ",fill=TRUE) sink() # Bundle new data win.data <- list(N=91, n=c(600,415,276,220,150,142,120,100,100,85,84,71,69,55,50,50,40,32, 20,1227053,635,4046,340,169,964,317,445,134,256,252,939,187,59,426,100, 100,100,98,190,560,540,1720,290,236,300,2511,123,177,361,350,561,80,150, 100,50,50,292,337,48,137,589,220,100,2009,256,961,113,315,69,100,1409,80, 200,81,50,445,30,100,124,121,640,176,97,77,16,51,17,14,16,95,21), r=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,278,2,14,2,1,9,4,6,2,4,4,15,3,1, 8,2,2,2,2,4,14,14,47,8,7,9,79,4,6,13,13,21,3,6,4,2,2,12,14,2,6,29,11,5,102, 13,50,6,17,4,6,85,5,14,6,4,38,3,12,15,15,90,27,15,14,3,10,5,6,7,43,17) ) # Inits function inits <- function(){ list(p = runif(1)) } # # Set the quantities we are interested in # Parameters to estimate params <- c("p") # MCMC settings ni <- 35000 # TOTAL Number of iterations including Burn in nb <- 10000 #Number of iterations for Burn in nt <- 1 nc <- 2 # Load/check model, load data and inits # Start Gibbs sampler with debug=TRUE to see WB working out1 <- bugs(data = win.data, inits = inits, parameters.to.save = params, model.file = model, n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni, debug = TRUE, bugs.directory= bugs.directory, working.directory=bugs.workingdir) # Inferences # After Burn in period let's work inference assuming ergodic regime is reached print(out2$summary,dig=4) # ------------------------------------------------------------------------- # Results # ------------------------------------------------------------------------- # Have a look on MCMC chains ################################################### read.bugs(c("CODA1.txt", "CODA2.txt"), start=19901, end=20000)->zz100 read.bugs(c("CODA1.txt", "CODA2.txt"))->zz varnames(zz) #zzz<-zz[,1:4,drop=F] #toutes les chaines variables 1 à 4 #Z <- ZZZ[[1]] #Select first chain zz<-zz[,2:2,drop=F] #toutes les chaines variables "p", zz100<-zz100[,2:2,drop=F] # variables "p", par(mfrow=c(2,1)) traceplot(zz100) traceplot(zz) par(mfrow=c(2,1)) acfplot(zz,lag.max=15,auto.layout=F) autocorr.plot(zz) densityplot(zz100) densityplot(zz) # Have a look on posterior statistics # -------------------------------------------------------- # Draw marginal posterior density for p # -------------------------------------------------------- # windows() res <- 6 name_figure <- "figTD1Beta1.png" png(filename = name_figure, height = 500*res, width = 500*res, res=72*res) cise.text <- 1.5 par(mfrow = c(1,2), oma = c(1,1,1,1)) pi <- out1$sims.list$p hist(pi, freq=F,main="", ylim = c(0,max(density(pi)$y)*1.2), xlim = range(density(pi)$x), ylab = "", xlab = "", bty="n",col="darkgrey") lines(density(pi), xlim = c(0,1), col = "black", lty = 1, lwd = 2) abline(a=1,b=0, lty=2) mtext(side=1, text = expression(p), bty="n", line=3, cex=cise.text, font = 1) mtext(side=2, text = "Tirages MCMC", bty="n", line=3, cex=cise.text, font = 1) xticks=seq(min(density(pi)$x),max(density(pi)$x),l=1000) a=1329 b=1255727 plot(xticks, y=dbeta(xticks,a,b),main="", ylim = c(0,max(density(pi)$y)*1.2), xlim = range(density(pi)$x), ylab = "", xlab = "", typ="l",col="black",bty="n") abline(a=1,b=0, lty=2) mtext(side=1, text = expression(p), bty="n", line=3, cex=cise.text, font = 1) mtext(side=2, text = "Loi beta théorique", bty="n", line=3, cex=cise.text, font = 1) dev.off() ###############sans la donnée 21 # Bundle data win.data <- list(N=90, n=c(600,415,276,220,150,142,120,100,100,85,84,71,69,55,50,50,40,32, 20,635,4046,340,169,964,317,445,134,256,252,939,187,59,426,100, 100,100,98,190,560,540,1720,290,236,300,2511,123,177,361,350,561,80,150, 100,50,50,292,337,48,137,589,220,100,2009,256,961,113,315,69,100,1409,80, 200,81,50,445,30,100,124,121,640,176,97,77,16,51,17,14,16,95,21), r=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,14,2,1,9,4,6,2,4,4,15,3,1, 8,2,2,2,2,4,14,14,47,8,7,9,79,4,6,13,13,21,3,6,4,2,2,12,14,2,6,29,11,5,102, 13,50,6,17,4,6,85,5,14,6,4,38,3,12,15,15,90,27,15,14,3,10,5,6,7,43,17) ) # Start Gibbs sampler out2 <- bugs(data = win.data, inits = inits, parameters.to.save = params, model.file = model, n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni, debug = FALSE, bugs.directory= bugs.directory, working.directory=bugs.workingdir) # Inferences # After Burn in period let's work inference assuming ergodic regime is reached print(out2$summary,dig=4) # ------------------------------------------------------------------------- # Results # ------------------------------------------------------------------------- # Have a look on MCMC chains ################################################### read.bugs(c("CODA1.txt", "CODA2.txt"), start=19901, end=20000)->zz100 read.bugs(c("CODA1.txt", "CODA2.txt"))->zz varnames(zz) #zzz<-zz[,1:4,drop=F] #toutes les chaines variables 1 à 4 #Z <- ZZZ[[1]] #Select first chain zz<-zz[,2:2,drop=F] #toutes les chaines variables "p", zz100<-zz100[,2:2,drop=F] # variables "p", par(mfrow=c(2,1)) traceplot(zz100) traceplot(zz) par(mfrow=c(2,1)) acfplot(zz,lag.max=15,auto.layout=F) autocorr.plot(zz) densityplot(zz100) densityplot(zz) # Draw marginal posterior density for p # -------------------------------------------------------- # windows() res <- 6 name_figure <- "figTD1Beta1bis.png" png(filename = name_figure, height = 500*res, width = 500*res, res=72*res) cise.text <- 1.5 par(mfrow = c(1,2), oma = c(1,1,1,1)) pi <- out2$sims.list$p hist(pi, freq=F,main="", ylim = c(0,max(density(pi)$y)*1.2), xlim = range(density(pi)$x), ylab = "", xlab = "", bty="n",col="darkgrey") lines(density(pi), xlim = c(0,1), col = "black", lty = 1, lwd = 2) abline(a=1,b=0, lty=2) mtext(side=1, text = expression(p), bty="n", line=3, cex=cise.text, font = 1) mtext(side=2, text = "Tirages MCMC", bty="n", line=3, cex=cise.text, font = 1) xticks=seq(min(density(pi)$x),max(density(pi)$x),l=1000) a=1051 b=28952 plot(xticks, y=dbeta(xticks,a,b),main="", ylim = c(0,max(density(pi)$y)*1.2), xlim = range(density(pi)$x), ylab = "", xlab = "", typ="l",col="black",bty="n") abline(a=1,b=0, lty=2) mtext(side=1, text = expression(p), bty="n", line=3, cex=cise.text, font = 1) mtext(side=2, text = "Loi beta théorique", bty="n", line=3, cex=cise.text, font = 1) dev.off() ############################################################### #Modelisation hierarchique ############################################################### # Define model model <- "model_Postulat2.txt" sink(model) cat(" model { for( i in 1 : N ) { p[i] ~ dbeta(alpha,beta) r[i] ~ dbin(p[i],n[i]) } alpha<-p.mean*taille beta<-(1-p.mean)*taille p.mean ~ dunif(0,1) taille ~ dexp(0.001) ppred ~ dbeta(alpha,beta) } ",fill=TRUE) sink() # Bundle new data win.data <- list(N=91, n=c(600,415,276,220,150,142,120,100,100,85,84,71,69,55,50,50,40,32, 20,1227053,635,4046,340,169,964,317,445,134,256,252,939,187,59,426,100, 100,100,98,190,560,540,1720,290,236,300,2511,123,177,361,350,561,80,150, 100,50,50,292,337,48,137,589,220,100,2009,256,961,113,315,69,100,1409,80, 200,81,50,445,30,100,124,121,640,176,97,77,16,51,17,14,16,95,21), r=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,278,2,14,2,1,9,4,6,2,4,4,15,3,1, 8,2,2,2,2,4,14,14,47,8,7,9,79,4,6,13,13,21,3,6,4,2,2,12,14,2,6,29,11,5,102, 13,50,6,17,4,6,85,5,14,6,4,38,3,12,15,15,90,27,15,14,3,10,5,6,7,43,17) ) # Inits function inits <- function(){ list(p = runif(91),ppred=runif(1),p.mean=runif(1),taille=rexp(1)) } # # Set the quantities we are interested in # Parameters to estimate params <- c("p","ppred","alpha","beta","p.mean","taille") # MCMC settings ni <- 60000 # TOTAL Number of iterations including Burn in nb <- 10000 #Number of iterations for Burn in nt <- 1 nc <- 2 # Load/check model, load data and inits # Start Gibbs sampler with debug=TRUE to see WB working out3 <- bugs(data = win.data, inits = inits, parameters.to.save = params, model.file = model, n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni, debug = FALSE, bugs.directory= bugs.directory, working.directory=bugs.workingdir) # Inferences # After Burn in period let's work inference assuming ergodic regime is reached print(out3$summary,dig=4) names(out3$sims.list) dim(out3$sims.list$p) p= out3$sims.list$p boxplot(log(p[1:100000,]),outline=F) res <- 6 name_figure <- "figTD1Boxplotp.png" png(filename = name_figure, height = 500*res, width = 500*res, res=72*res) cise.text <- 1.5 par(mfrow = c(1,1), oma = c(1,1,1,1)) xticks=seq(1,nc*(ni-nb),length.out=1000) boxplot(log(p[xticks,]/(1-p[xticks,])),outline=F) mtext(side=1, text = expression(p), bty="n", line=3, cex=cise.text, font = 1) mtext(side=2, text = "Incertitude intra-référence", bty="n", line=3, cex=cise.text, font = 1) dev.off() name_figure <- "figTD1Beta2.png" png(filename = name_figure, height = 500*res, width = 500*res, res=72*res) cise.text <- 1.5 par(mfrow = c(3,2), oma = c(1,1,1,1)) pi <- out3$sims.list$alpha hist(pi, freq=F,main="", ylim = c(0,max(density(pi)$y)*1.2), xlim = range(density(pi)$x), ylab = "", xlab = "", bty="n",col="darkgrey") lines(density(pi), xlim = c(0,1), col = "black", lty = 1, lwd = 2) mtext(side=1, text = expression(alpha), bty="n", line=3, cex=cise.text, font = 1) mtext(side=2, text = "Tirages MCMC", bty="n", line=3, cex=cise.text/2, font = 1) #beta pi <- out3$sims.list$beta hist(pi, freq=F,main="", ylim = c(0,max(density(pi)$y)*1.2), xlim = range(density(pi)$x), ylab = "", xlab = "", bty="n",col="darkgrey") lines(density(pi), xlim = c(0,1), col = "black", lty = 1, lwd = 2) mtext(side=1, text = expression(beta), bty="n", line=3, cex=cise.text, font = 1) #taille pi <- out3$sims.list$taille hist(pi, freq=F,main="", ylim = c(0,max(density(pi)$y)*1.2), xlim = range(density(pi)$x), ylab = "", xlab = "", bty="n",col="darkgrey") lines(density(pi), xlim = c(0,1), col = "black", lty = 1, lwd = 2) mtext(side=1, text = expression(taille), bty="n", line=3, cex=cise.text, font = 1) mtext(side=2, text = "Tirages MCMC", bty="n", line=3, cex=cise.text/2, font = 1) # p.mean pi <- out3$sims.list$p.mean hist(pi, freq=F,main="", ylim = c(0,max(density(pi)$y)*1.2), xlim = range(density(pi)$x), ylab = "", xlab = "", bty="n",col="darkgrey") lines(density(pi), xlim = c(0,1), col = "black", lty = 1, lwd = 2) mtext(side=1, text = expression(p.mean), bty="n", line=3, cex=cise.text, font = 1) # ppred pi <- out3$sims.list$ppred hist(pi, freq=F,main="", ylim = c(0,max(density(pi)$y)*1.2), xlim = range(density(pi)$x), ylab = "", xlab = "", bty="n",col="darkgrey") lines(density(pi), xlim = c(0,1), col = "black", lty = 1, lwd = 2) mtext(side=1, text = expression(ppred), bty="n", line=3, cex=cise.text, font = 1) mtext(side=2, text = "Tirages MCMC", bty="n", line=3, cex=cise.text/2, font = 1) dev.off() ############################################################### #Modelisation hierarchique avec deux sous-populations ############################################################### # Define model model <- "model_Postulat3.txt" sink(model) cat(" model { P[1:2] ~ ddirch(A[1:2]) for( i in 1 : N ) { w[i] ~ dcat(P[1:2]) p[i] ~ dbeta(alpha[w[i]],beta[w[i]]) r[i] ~ dbin(p[i],n[i]) proba2[i]<-step(w[i]-2) } for( j in 1 : 2 ) { alpha[j]<-p.mean[j]*taille[j] beta[j]<-(1-p.mean[j])*taille[j] taille[j] ~ dexp(0.001) ppred[j] ~ dbeta(alpha[j],beta[j]) } p.mean[1]~dbeta(2,18) p.mean[2]~dbeta(1,1) } ",fill=TRUE) sink() # Bundle new data win.data <- list(N=91,A=c(1,1), n=c(600,415,276,220,150,142,120,100,100,85,84,71,69,55,50,50,40,32, 20,1227053,635,4046,340,169,964,317,445,134,256,252,939,187,59,426,100, 100,100,98,190,560,540,1720,290,236,300,2511,123,177,361,350,561,80,150, 100,50,50,292,337,48,137,589,220,100,2009,256,961,113,315,69,100,1409,80, 200,81,50,445,30,100,124,121,640,176,97,77,16,51,17,14,16,95,21), r=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,278,2,14,2,1,9,4,6,2,4,4,15,3,1, 8,2,2,2,2,4,14,14,47,8,7,9,79,4,6,13,13,21,3,6,4,2,2,12,14,2,6,29,11,5,102, 13,50,6,17,4,6,85,5,14,6,4,38,3,12,15,15,90,27,15,14,3,10,5,6,7,43,17), w=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,1, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,2) ) # Inits function inits <- function(){ list(p = runif(91),ppred=runif(2),p.mean=runif(2),taille=rexp(2)) } # Set the quantities we are interested in # Parameters to estimate params <- c("P","ppred","alpha","beta","p.mean","taille") # MCMC settings ni <- 150000 # TOTAL Number of iterations including Burn in nb <- 50000 #Number of iterations for Burn in nt <- 1 nc <- 2 # Load/check model, load data and inits # Start Gibbs sampler with debug=TRUE to see WB working out3 <- bugs(data = win.data, inits = inits, parameters.to.save = params, model.file = model, n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni, debug = TRUE, bugs.directory= bugs.directory, working.directory=bugs.workingdir) # Inferences # After Burn in period let's work inference assuming ergodic regime is reached print(out3$summary,dig=4) names(out3$sims.list) name_figure <- "figTD1Beta3.png" res<-6 png(filename = name_figure, height = 500*res, width = 500*res, res=72*res) cise.text <- 1.5 par(mfrow = c(2,1), oma = c(1,1,1,1)) # ppred[1] pi <- out3$sims.list$ppred[,1] hist(pi, freq=F,main="", ylim = c(0,max(density(pi)$y)*1.2), xlim = range(density(pi)$x), ylab = "", xlab = "", bty="n",col="darkgrey") lines(density(pi), xlim = c(0,1), col = "black", lty = 1, lwd = 2) mtext(side=1, text = "ppred[1]", bty="n", line=3, cex=cise.text, font = 1) mtext(side=2, text = "Tirages MCMC", bty="n", line=3, cex=cise.text/2, font = 1) # ppred[2] pi <- out3$sims.list$ppred[,2] hist(pi, freq=F,main="", ylim = c(0,max(density(pi)$y)*1.2), xlim = range(density(pi)$x), ylab = "", xlab = "", bty="n",col="darkgrey") lines(density(pi), xlim = c(0,1), col = "black", lty = 1, lwd = 2) mtext(side=1, text = "ppred[2]", bty="n", line=3, cex=cise.text, font = 1) mtext(side=2, text = "Tirages MCMC", bty="n", line=3, cex=cise.text/2, font = 1) dev.off()