#******************************************************************************* # 12.2.2 Using BUGS from R #******************************************************************************* # # Fränzi Korner-Nievergelt, Tobias Roth, Stefanie von Felten, Jérôme Guélat, # Bettina Almasi, Pius Korner-Nievergelt, 2015. Bayesian data analysis in ecology # using linear models with R, BUGS and Stan. Elsevier. # # Last modification: 23 April 2014 #------------------------------------------------------------------------------- # Settings #------------------------------------------------------------------------------- library(R2OpenBUGS) bugslocation <- "C:/Program Files/OpenBUGS323/OpenBugs.exe" # location of OpenBUGS bugsworkingdir <- "C:/Users/name/BUGS" # BUGS working directory if different from R working directory #------------------------------------------------------------------------------- # Simulate fake data #------------------------------------------------------------------------------- n <- 50 # sample size sigma <- 5 # standard deviation of the residuals b0 <- 2 # intercept b1 <- 0.7 # slope x <- runif(n, 10, 30) # random numbers of the covariate simresid <- rnorm(n, 0, sd=sigma) # residuals y <- b0 + b1*x + simresid # calculate y, i.e. the data #------------------------------------------------------------------------------- # Function to generate initial values #------------------------------------------------------------------------------- inits <- function() { list(beta0=runif(1, -2, 2), beta1=runif(1, -2, 2), sigmaRes=runif(1, 0.1, 2)) } #------------------------------------------------------------------------------- # Bundle data into a list #------------------------------------------------------------------------------- datax <- list(n=length(y), y=y, x=x) #------------------------------------------------------------------------------- # Run OpenBUGS #------------------------------------------------------------------------------- parameters <- c("beta0", "beta1", "sigmaRes") fit <- bugs(datax, inits, parameters, model.file="normlinreg.txt", n.thin=1, n.chains=2, n.burnin=5000, n.iter=10000, debug=TRUE, OpenBUGS.pgm = bugslocation, working.directory=bugsworkingdir) # if OpenBUGS has been installed by the installer, the bugslocation # does not need to be indicated fit <- bugs(datax, inits, parameters, model.file="normlinreg.txt", n.thin=1, n.chains=2, n.burnin=5000, n.iter=10000, debug=TRUE, working.directory=bugsworkingdir) #------------------------------------------------------------------------------- # Test convergence and make inference #------------------------------------------------------------------------------- library(blmeco) # Make Figure 12.2 jpeg(filename="Figure12-2_bugs.jpg", width=6000, height=6000, units="px", res=1200) par(mfrow=c(3,1)) historyplot(fit, "beta0") historyplot(fit, "beta1") historyplot(fit, "sigmaRes") dev.off() # Parameter estimates print(fit$summary, 3) # Make predictions for covariate values between 10 and 30 newdat <- data.frame(x=seq(10, 30, length=100)) Xmat <- model.matrix(~x, data=newdat) predmat <- matrix(ncol=fit$n.sim, nrow=nrow(newdat)) for(i in 1:fit$n.sim) predmat[,i] <- Xmat%*%c(fit$sims.list$beta0[i], fit$sims.list$beta1[i]) newdat$lower.bugs <- apply(predmat, 1, quantile, prob=0.025) newdat$upper.bugs <- apply(predmat, 1, quantile, prob=0.975) # Make Figure 12.3 jpeg(filename="Figure12-3_bugs.jpg", width=6000, height=6000, units="px", res=1200) plot(y~x, pch=16, las=1, cex.lab=1.4, cex.axis=1.2, type="n", main="") polygon(c(newdat$x, rev(newdat$x)), c(newdat$lower.bugs, rev(newdat$upper.bugs)), col=grey(0.7), border=NA) abline(c(fit$mean$beta0, fit$mean$beta1), lwd=2) box() points(x,y) dev.off()