#******************************************************************************* # 12.3 MCMC using STAN #******************************************************************************* # # Korner-Nievergelt, F., T. Roth, S. von Felten, J. Guelat, B. Almasi and # P. Korner-Nievergelt. 2015. Bayesian Data Analysis in Ecology using Linear # Models with R, BUGS and Stan. Elsevier. # # Last modification: 23 April 2014 #------------------------------------------------------------------------------- # Settings #------------------------------------------------------------------------------- library(rstan) #------------------------------------------------------------------------------- # 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 #------------------------------------------------------------------------------- # Bundle data into a list #------------------------------------------------------------------------------- datax <- list(n=length(y), y=y, x=x) #------------------------------------------------------------------------------- # Run STAN #------------------------------------------------------------------------------- fit <- stan(file = "STAN/linreg.stan.txt", data=datax, chains=10, iter=1000) #------------------------------------------------------------------------------- # Test convergence and make inference #------------------------------------------------------------------------------- print(fit, c("beta", "sigma")) traceplot(fit, "beta") traceplot(fit, "sigma") modsims <- extract(fit) str(modsims) # Make predictions for covariate values between 10 and 30 newdat <- data.frame(x=seq(10, 30, length=100)) Xmat <- model.matrix(~x, data=newdat) b <- apply(modsims$beta, 2, mean) newdat$fit <- Xmat%*%b nsim <- nrow(modsims$beta) fitmat <- matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim){ fitmat[,i] <- Xmat%*%modsims$beta[i,] } newdat$fitlwr <- apply(fitmat, 1, quantile, prob=0.025) newdat$fitupr <- apply(fitmat, 1, quantile, prob=0.975)