tutorial - Johannes Karreth
Contents
1. iH THE betai 0 47904271 0 5804246 1 541857 beta2 1 89980599 0 2904349 1 577197 deviance 0 05770752 0 4806812 1 169611 Window From Start 0 10000000 0 4542000 0 719760 Window From Stop 0 50000000 0 0827600 0 226630 The Gelman Rubin diagnostic Potential scale reduction factors Point est Upper alpha 1 1 03 betai 15201 1 04 beta2 1 00 1 00 deviance 1 00 1 00 Multivariate psrf 1 01 The Heidelberger Welch diagnostic Chain 1 epsilon 0 1 alpha 0 05 Stationarity start p value test iteration alpha passed T 0 558 betai passed 1 0 441 beta2 passed 4 0 491 deviance passed 0 477 Halfwidth Mean Halfwidth test alpha passed 1 87 0 03087 betai passed 1 13 0 00553 beta2 passed e3 09 0 013850 deviance passed 271 76 0 18102 Chain 2 epsilon 0 062 alpha 0 1 Stationarity start p value test iteration alpha passed au 0 178 betai passed 1 0 300 beta2 passed f 0 907 deviance passed 1 0 752 Halfwidth Mean Halfwidth test alpha passed 1 86 0 03183 betai passed 113 0 00595 beta2 passed 3 09 0 01363 deviance passed 271 84 0 21074 Chain 3 epsilon 0 111 alpha 0 05 Stationarity start p value test iteration alpha passed a 0 219 betai passed f YNO beta2 passed 1 0 122 deviance passed 1 0 145 H2. Coda files loaded successfully Calculating summary statistics Calculating the Gelman Rubin statistic for 3 variables Finished running the simulation JAGS files were saved to the runjagsfiles 2 folder in your current working directory Note a few things e The model is read into R within quotation marks e setting keep jags files TRUE we get run jags to create a folder runjagsfiles in your WD This folder will contain the same files you obtain from JAGS in the next step below This can be useful for further processing e The option method allows for parallel chains on separate cores of your computer and other higher end computing options See the help file for run jags for more information 5 1 Output and diagnostics runjags also allows you to use the same set of commands for diagnostics and accessing the fitted objects as you explored above e Summarize the model object gt print bayes mod fit Inference for Bugs model at var folders 8j p9ngyxk5295bnlk0n50mm5ch0000gg T RtmpgkE7yZ model470543d32d1 txt fit using 3 chains each with 9000 iterations first 1000 discarded n thin 8 n sims 3000 iterations saved dt mu vect sd vect 25 50 75 alpha 1 854 0 289 1291 1 660 1 863 2 054 betal 1 134 0 052 1 036 1 097 1 134 1 168 beta2 3 092 210 34296 3 239 3 087 2 943 deviance 271 766 2 871 268 240 269 647 271 079 273 257 RH 97 5 Rhat n ef
3. Fraction in 1st window Fraction in 2nd window d HH alpha betai beta2 deviance 0 3987 0 1764 0 8337 1 0954 di d 3 di Fraction in 1st window Fraction in 2nd window d alpha betai beta2 deviance di 1 477 1 651 1 003 182 omi 0 5 027 0 5 0 1 0 5 gt geweke plot bayes mod fit mcmc gt raftery diag bayes mod fit mcmc 1 Quantile q 0 025 Accuracy r 0 005 Probability s 0 95 You need a sample size of at least 3746 with these values of q r and s 2 Quantile q 0 025 Accuracy r 0 005 Probability s 0 95 You need a sample size of at least 3746 with these values of q r and s 3 Quantile q 0 025 Accuracy r 0 005 Probability s 0 95 You need a sample size of at least 3746 with these values of q r and s gt heidel diag bayes mod fit mcmc 1 dH Stationarity start p value test iteration alpha passed ii 0 185 betal passed 1 0 138 10 beta2 passed 1 0 589 deviance passed 1 0 827 Halfwidth Mean Halfwidth test alpha passed 02971 betal passed 1 13 0 00527 beta2 passed 3 09 0 01308 deviance passed 271 74 0 18831 21 b Stationarity start p value test iteration alpha passed 1 0 157 betal passed 1 0 393 beta2 passed 1 0 982 deviance
4. e And we create and summarize a frequentist linear model fit on these data gt gt Ht Ht freq mod lt 1 x1 x2 data sim dat summary freq mod lm formula y 1 x2 data sim dat Residuals Min 1Q Median 3Q Max 21 3432 O0 6797 0 1112 0 5367 3 2204 Coefficients Estimate Std Error t value gt Intercept 1 84949 0 28810 6 42 5 04e 09 1 1 13511 0 05158 22 00 lt 2e 16 x2 3 09361 0 20650 14 98 lt 2e 16 Signif codes OJ O DO ee um cO dote 0b Residual standard error 0 9367 on 97 degrees of Multiple R squared 0 8772 Adjusted R squared F statistic 346 5 on 2 and 97 DF p value lt 2 freedom 0 8747 2e 16 Now we write a model for JAGS and save it as the text file bayes mod in your working directory Do not paste this model straight into R yet You can set your working directory via gt setwd Documents Uni 9 ICPSR 2015 Applied Bayes Tutorials 2 JAGS and R The model looks just like the JAGS models shown in throughout this course model for i in 1 N ylil mu i alpha 7 betai 7 beta2 dnorm mu i tau alpha betai x1 i beta2 x2 i dnorm 0 01 dunif 100 100 dunif 100 100 tau dgamma 01 01 Instead of saving the model in your WD you can also enter it in your R script gt bayes
5. 32 deviance 0 15 0 10 0 05 0 00 1 T 275 280 value e Or you can use the ggmcmc command to create a PDF file containing a variety of diagnostic plots gt ggmcmc bayes mod fit gg oF file Documents Uni 9 ICPSR 2015 Applied Bayes Tutorials 2 JAGS and R bayes fit ggmcmc pdf Plotting histograms Plotting density plots Plotting traceplots Plotting running means Plotting comparison of partial and full chain Plotting autocorrelation plots Plotting crosscorrelation plot Plotting Potential Scale Reduction Factors Plotting Geweke Diagnostic Plotting caterpillar plot pdf 2 5 Using runjags Instead of R2jags you can also use the runjags package to command JAGS via R Here is an abbre viated version of the workflow for runjags The workflow mimics what you saw for R2jags above gt library runjags For an example dataset we simulate our own data in R For this tutorial we aim to fit a linear model so we create a continuous outcome variable y as a function of two predictors x and x2 and a disturbance term e We simulate 100 observations e First we create the predictors gt n sim lt 100 set seed 123 gt 1 lt rnorm n sim mean 5 sd 2 gt x2 lt rbinom n sim size 1 prob 0 3 gt e lt rnorm n sim mean 0 sd 1 e Next we create the outcome y based on coefficients b and b for the respective predictors and an intercept
6. or Notepad for Windows http notepad plus plus org RStudio is a very neat integrated environment for using R but a separate text editor will be useful from time to time to inspect JAGS model files and other files 9 Note for users of Mac OS X 10 5 Leopard Due to a particular behavior of the JAGS installer on Leopard the JAGS files that rjags requires to run are not located where rjags is looking for them See http martynplummer wordpress com 2011 11 04 rjags 3 for mac cs x If you would like to use R2jags or rjags on Mac OS X 10 5 you need to manually relocate these files from usr to usr local Ask me if you would like help with this 10 You should now be able to run the following code in R taken directly from the help file for the jags function library R2jags model file lt system file package R2jags model schools txt J lt 8 0 Vi lt cC 29 27 9 2 9 60 8 0 0 11079 0561820 12 2 sd lt c 14 9 10 2 16 3 11 0 9 4 1 B dein 1 4 10 4 17 6 jags data lt list y sd J jags params lt c mu sigma theta jags inits lt function list mu rnorm 1 sigma runif 1 theta rnorm J gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt 1 f ae gt jagsfit lt jags data list y sd J inits jags inits jags params n iter 10 model file model file Compiling model graph Resolving unde
7. With an MCMC object you can use a variety of commands for diagnostics and presentation using the CODA package e Plot gt library lattice gt xyplot bayes mod fit mcmc e You can customize the plot layout you can use other Lattice options here as well gt xyplot bayes mod fit mcmc layout c 2 2 aspect fill e Density plot gt densityplot bayes mod fit mcmc gt densityplot bayes mod fit mcmc layout c 2 2 aspect fill beta2 deviance w 54 54 2g 84 s oJ 8 4 e S gt T T T T a 40 35 30 25 265 270 275 280 285 290 a alpha betat a or 4 o a4 84 T 10 15 20 25 3 0 09 10 11 12 13 e Trace and density in one plot print directly to your working directory gt pdf Documents Uni 9 ICPSR 2015 Applied Bayes Tutorials 2 JAGS and R bayes fit mcmc plot pdf gt plot bayes mod fit mcmc gt dev off e Autocorrelation plot print directly to your working directory gt pdf Documents Uni 9 ICPSR 2015 Applied Bayes Tutorials 2 JAGS and R bayes_fit_mcmc_autocorr pdf gt autocorr plot bayes mod fit mcmc gt dev off Other diagnostics using CODA gt gelman plot bayes mod fit mcmc gt geweke diag bayes mod fit mcmc 11 di Fraction in 1st window Fraction in 2nd window d alpha betai beta2 deviance 1 65474 1 71754 0 03899 0 30234 di d 2 di
8. a bi lt 1 2 b2 lt 3 1 a lt 1 5 gt gt gt gt y lt a table 10 26 15 e Now we combine the variables into one dataframe for processing later gt sim dat lt data frame y x1 x2 e And we summarize a frequentist linear model fit on these data gt freg mod lt lm y xi x2 data sim dat gt summary freq mod RH Call lm formula y x1 x2 data sim dat RH Residuals Min 1Q Median 39 Max 1 3432 0 6797 0 1112 0 5367 3 2304 RH Coefficients dt Estimate Std Error t value Pr gt t Intercept 1 84949 0 28810 6 42 5 04e 09 x1 1 19511 0 05158 22 00 lt 2e 16 x2 3 09361 0 20650 14 98 lt 2e 16 W Signif codes HH Taek 07001 ee OLOL 0 0b 7 0 1 t3 RH Residual standard error 0 9367 on 97 degrees of freedom Multiple R sguared 0 8772 Adjusted R squared 0 8747 F statistic 346 5 on 2 and 97 DF p value lt 2 2e 16 Create the model object gt bayes mod lt model for i in 1 N y i dnorm mu i tau mu i lt alpha betai x1 i beta2 x2 i alpha dnorm 0 01 betai dunif 100 100 beta2 dunif 100 100 tau dgamma 01 01 pu Convert the data frame to a list gt sim list lt as list sim dat and add the number of observations gt sim list N lt nrow sim dat Convert the data for runjags gt sim dat runjags dump f
9. course on my web site at http www jkarreth net bayes icpsr html The most recent version of this file is posted at http www jkarreth net files Lab2_JAGS Many JAGS related questions are answerd at http stackoverflow com questions tagged and the discussion board at http sourceforge net projects mcmc jags forums forum 610037 10 Other software solutions for Bayesian estimation Other options to fit Bayesian models include e WinBUGS OpenBUGS These programs have a GUI but can also be accessed from R with similar functionality A convenient way to fit Bayesian models using WinBUGS or OpenBUGS is to use R packages that function as frontends for WinBUGS or OpenBUGS such as R2WinBUGS or BRugs 22006 R2WinBUGS works very similar to R2jags so you can adapt the code in this tutorial easily Mac users can find setup instructions on my course website e Stan http mc stan org is anew and fast program that can be accessed via R and other statistical packages The documentation on the Stan website is very easy to follow and offers tutorials on fitting some example models in Stan We will show you how to use Stan in Week 4 of this course 23 References Curtis S McKay 2012 mcmcplots Create Plots from MCMC Output R package version 0 4 1 URL hittp CRAN R project org package mcmcplots Denwood Matthew J Under review runjags An R package providing interface utilities paral lel computing methods an
10. each line of the script does model clear data clear load dic Remove previous data and models if applicable load the dic module so you can monitor the model deviance later model in bayes mod Use the model bayes mod which is saved in your WD and looks like a regular JAGS model Make sure you use the exact and full name of the model file as it is in your working directory otherwise JAGS will not find it Look out for hidden file name extensions You wrote this model earlier and saved it in your WD as bayes mod model for i in 1 N y i 7 dnorm mu i tau mu i lt alpha betai x1 i beta2 x2 i alpha dnorm 0 01 betai dunif 100 100 beta2 dunif 100 100 tau dgamma 01 01 Y data in sim dat Use the data sim dat These data can be saved as vectors in one file that you name sim dat This file will look something like this y c 6 94259729574987 4 61661626624104 x1 c 3 87904870689558 4 53964502103344 x2 lt 0 1 N 100 You can create this data file by hand inconvenient or you can work with some R commands to do this automatically If sim dat is a data frame object in R you can do the following 3On the Mac you can set the Finder to display all file extensions by going to Finder gt Preferences gt check Show all filename extensions 19 sim dat list lt as list sim dat sim dat list N lt nrow sim dat dump sim dat
11. list file sim dat dump bugs2jags sim dat dump sim dat The first command converts the data frame into list the second command writes out the data in BUGS format as a file to your working directory the third command part of the coda package translates it into JAGS format Both angell dump and angell dat are written to your working directory so be sure that you have specified that in the beginning Also be sure to visually inspect whether your data file was created correctly compile nchains 2 Compile the models and run two Markov chains inits in bayes mod initsi chain 1 inits in bayes mod inits2 chain 2 inits in bayes mod inits3 chain 3 Use the starting values you provide in bayes mod inits1 bayes mod inits2 and bayes mod inits3 each of which could look like this alpha lt 0 betai lt 0 beta2 lt c 0 This way you can specify different starting values for each chain initialize Initialize and run the model update 2500 by 100 Specify 2500 updates before you start monitoring your parameters of interest monitor alpha thin 2 monitor betai thin 2 monitor beta2 thin 2 monitor deviance thin 2 Monitor the values for these parameters in this case the three regression coefficients and the model deviance thin 2 specifies that only each second draw from the chains will be used for your model output update 2500 by 100 c
12. passed 1 0 604 Halfwidth Mean Halfwidth test alpha passed 1 86 0 0296 betal passed 1 13 0 005B beta2 passed o 097070129 deviance passed 271 80 0 2042 3 Stationarity start p value test iteration alpha passed 1 0 4468 betal passed 1 0 3556 beta2 passed 1 0 0964 deviance passed 101 0 0839 Halfwidth Mean Halfwidth test alpha passed 1 83 0 03066 betal passed 1 14 0 00568 beta2 passed 3 10 0 01360 deviance passed 271 81 0 20054 Quick convergence diagnostics superdiag A very convenient function to analyze numerical representations of diagnostics in one sweep is the superdiag package Tsai Gill and Rapkin 2012 e First install the package gt install packages superdiag dependencies TRUE repos http cran us r project org e Then load it gt library superdiag e Next all you need to do is gt superdiag bayes mod fit mcmc burnin 100 Number of chains 3 Number of iterations 1000 per chain before discarding the burn in period The burn in period 100 per chain Sample size in total 2700 The Geweke diagnostic Z scores chain1 chain 2 chain 3 alpha 0 32702322 0 8681148 1 804617 11
13. the 1 command in R and produces the same MCMC output you analyzed above One potential downside of MCMCpack is that it does not offer as much options for customization as writing your full model JAGS For more information on MCMCpack see http mcmcpack Because MCMCpack uses an R like formula it is straightforward to fit a Bayesian linear model gt library MCMCpack For an example dataset we again simulate our own data in R We create a continuous outcome variable y as a function of two predictors x and x2 and a disturbance term e We simulate 100 Observations e First we create the predictors n sim 100 set seed 123 xi lt rnorm n sim mean 5 sd 2 x2 lt rbinom n sim size 1 prob 0 3 e lt rnorm n sim mean 0 sd 1 NONE NN e Next we create the outcome y based on coefficients b and b for the respective predictors and an intercept a bl lt 1 2 53011 lt 1 5 y lt a k pil ri 185206 2 tre MOM N e Now we combine the variables into one dataframe for processing later gt sim dat lt data frame y x1 x2 e And we summarize a frequentist linear model fit on these data freq mod lt lm y x2 data sim dat gt summary freq mod Call lm formula 1 x2 data sim dat dH Residuals Min 1Q Median 3Q Max 1 3432 0 6797 0 1112 0 5367 3 2304 Coefficients d Estimate Std Erro
14. 3 dmg library from the CRAN tools directory http cran r project org bin macosx tools 3 Install JAGS version 3 4 0 from Martyn Plummer s repository http sourceforge net projects mcmc jags files JAGS 3 x 4 Start the Terminal Mac or Console Windows and type jags If JAGS is installed you will receive the following message Welcome to JAGS 3 4 0 on Tue Dec 23 20 11 23 2014 JAGS is free software and comes with ABSOLUTELY NO WARRANTY Loading module basemod ok Loading module bugs ok You can guit the Terminal Console again 5 Install the packages R2jags coda R2WinBUGS lattice and lastly rjags from within R or RStudio via the Package Installer or by using gt install packages R2jags dependencies TRUE repos http cran us r project org The dependencies TRUE option will install the aforementioned packages automatically 6 Install the package runjags from within R or RStudio via the Package Installer or by using gt install packages runjags dependencies TRUE repos http cran us r project org 7 Install the package MCMCpack from within R or RStudio via the Package Installer or by using gt install packages MCMCpack dependencies TRUE repos http cran us r project org 8 Download a scientific text editor for writing R and JAGS code I recommend TextWrangler for Mac http www barebones com products textwrangler Sublime Text http www sublimetext com
15. 810 beta2 2 692 1 001 3000 deviance 278 737 1 001 3000 iH For each parameter n eff is a crude measure of effective sample size and Rhat is the potential scale reduction factor at convergence Rhat 1 DIC info using the rule pD var deviance 2 pD 4 1 and DIC 275 9 DIC is an estimate of expected predictive error lower deviance is better gt plot bayes mod fit gt traceplot bayes mod fit j pOngyxk5295bnlk0n50mm5ch0000ga T RtmpgkE7yZ model470543d32d1 txt fit using jags 3 chains each with 9000 ite 80 interval for each chain R hat medians and 80 intervals 4 2 0 2 4 1 15 2 alpha Ro n betat H DRRA dm v 28 4 2 0 2 4 1 15 2 2 alpha 1 5 1 1 254 12 betal 1 155 A 1 054 28 3 beta2 82 34 280 275 deviance 270 265 alpha a oO o T T T T T T T T T 1 100 200 300 400 500 600 700 800 900 1000 iteration beta1 o 100 200 300 400 500 600 700 800 900 1000 iteration beta2 2 5 3 0 beta2 3 5 T T T T T T T T T 100 200 300 400 500 600 700 800 900 1000 iteration deviance deviance T T T T T T T T T 100 200 300 400 500 600 700 800 900 1000 iteration If you want to print and save the plot you can use the following set of commands gt pdf Documents Uni 9 ICPSR 2015 Applied Bayes Tutorials 2 JAGS and R bayes
16. Using JAGS via R Johannes Karreth University at Albany SUNY jkarreth albany edu ICPSR Summer Program 2015 All code used in this tutorial can also be found on my course website at http www jkarreth Updated code will be posted there If you find any errors in this document or have suggestions for improvement please don t hesitate to email me Contents 1 Using as frontend 2 2 What are JAGS R2jags 2 3 Installing JAGS and R2jags 2 4 Fitting Bayesian models using R2jags 3 4 1 Preparing 4 4 2 Fitting t emodel 6 4 3 1 5 a YES sqb sees YSOL A UD EM s 6 5 Using runjags 15 5 1 Output and diagnostics 2 FF I I I I Y 17 6 Using JAGS via the command line 18 7 Using the coda package 20 8 MCMCpack 22 8 1 23 9 Following this course using JAGS 23 10 Other software solutions for Bayesian estimation 23 This tutorial focuses on using JAGS for fitting Bayesian models via R There are other options for fitting Bayesian models that we will briefly discuss during the workshop You can find more information about them at the end of this tutorial and on my workshop website 1 Using R as frontend A convenient way to fit Bayesian models using JAGS or WinBUGS or OpenBUGS is to use R packages that function a
17. _trace pdf defines that the plot will be saved as a PDF file with the name bayes_trace pdf in your working directory gt traceplot bayes mod fit creates the plot in the background you will not see it gt dev off finishes the printing process and creates the PDF file of the plot If successful R will display the message null device 1 More diagnostics are available when you convert your model output into an MCMC object You can generate an MCMC object for analysis with this command gt bayes mod fit mcmc lt as mcmc bayes mod fit gt summary bayes mod fit mcmc 2IATEX cannot process file names with periods so if you use and try to include the graphics file angell trace IATEX will not compile your document Iterations 1001 8993 Thinning interval 8 Number of chains 3 Sample size per chain 1000 1 Empirical mean and standard deviation for each variable plus standard error of the mean Mean SD Naive SE Time series SE alpha 1 854 0 28863 0 0052697 0 008836 betal 1 134 0 05185 0 0009466 0 001616 beta2 3 092 0 20974 0 0038294 0 003890 deviance 271 766 2 87084 0 0524142 0 057127 2 Quantiles for each variable 25 25 50 75h 97 5 alpha 1 281 1 660 1 863 2 051 27396 betal 1 036 1 097 1 134 1 168 15237 beta2 3 496 3 239 3 087 2 943 2 692 deviance 268 240 269 647 271 079 273 257 278 737
18. a 1250 valid values Abstracting betal 1250 valid values Abstracting beta2 1250 valid values Abstracting deviance 1250 valid values gt bayes chains lt as mcmc list list chain1 chain2 chain3 Now you can analyze using some of the commands listed in gt help package coda for instance gt summary bayes chains Iterations 2501 4999 Thinning interval 2 Number of chains 3 Sample size per chain 1250 1 Empirical mean and standard deviation for each variable plus standard error of the mean SD Naive SE Time series SE alpha 1 803 0 28359 0 0046310 0 014504 betal 1 143 0 05142 0 0008397 0 002786 beta2 3 088 0 20954 0 0034218 0 004091 deviance 271 768 2 93772 0 0479728 0 071076 2 Quantiles for each variable 2 25 50 75h 97 5 alpha 1 238 1 614 1 808 1 992 22355 betal 1 041 1 109 1 143 TS 1 244 beta2 3 496 3 228 3 089 2 949 2 678 deviance 268 138 269 671 271 082 273 063 279 708 gt traceplot bayes chains You should be able to play around with graphical parameters and different printing devices this way Convenient the superdiag mcmcplots and ggmcmc packages also work like described above once you have declared your JAGS output an MCMC object 21 8 MCMCpack MCMCpack Martin Quinn and Park 2011 is an R package that we will also use in this course It is as easy to use as
19. alfwidth Mean Halfwidth test alpha passed 1283 00328 betai passed 1 14 0 0060 beta2 passed 3 10 0 0137 deviance passed 271 81 0 2005 The Raftery Lewis diagnostic Chain 1 converge eps 0 001 Quantile q 0 025 Accuracy r 0 005 Probability s 0 95 12 You need a sample size of at least 3746 with these values of q r ands Chain 2 converge eps be 04 Quantile q 0 001 Accuracy r 0 005 Probability s 0 95 Burn in Total Lower bound Dependence M N Nmin factor I alpha 2 173 154 1512 betal 2 173 154 2 beta2 2 173 154 1 deviance 2 173 154 1512 Chain 3 converge eps 0 0025 Quantile q 0 25 Accuracy r 0 001 Probability s 0 99 You need a sample size of at least 1244044 with these values of q r ands Note that the R2jags object by default retains 1000 iterations through thinning hence the burn in period you provide for superdiag must be less than 1000 Quick diagnostic plots mcmcplots A convenient way to obtain graphical diagnostics and results is using the mcmcplots package 2012 e First install the package gt install packages mcmcplots dependencies TRUE repos http cran us r project org e Then load it gt library mcmcplots e Some commands for plots gt denplot bayes mod fit mcmc gt denplot bayes mod fit mcm
20. c parms c alpha betai beta2 gt traplot bayes mod fit mcmc parms c alpha betai beta2 alpha beta1 beta2 13 alpha beta1 0 0 2000 4000 6000 8000 beta2 0 0 2000 4000 6000 8000 As always check the help files for options to customize these plots e Or for quick diagnostics you can produce html files with trace density and autocorrelation plots all on one page The files will be displayed in your default internet browser gt mcmcplot bayes mod fit mcmc e If you want to produce a coefficient dot plot with credible intervals use caterplot gt caterplot bayes mod fit mcmc gt caterplot bayes mod fit mcmc parms c alpha betai beta2 labels c alpha betai beta2 alpha UU beta1 beta2 LE 1 3 2 1 0 1 2 More plots ggmcmc Yet another option for plotting output is the ggemcmc package Fernandez i Mar n 2013 e First install the package gt install packages ggmcmc dependencies TRUE repos http cran us r project org e Then load it gt library ggmcmc e To use this package you need to first convert your MCMC object into a ggs object using the ggs function This object can then be passed on to the various ggmcmc plotting commands 14 gt bayes mod fit gg lt ggs bayes mod fit mcmc gt ggs density bayes mod fit gg alpha beta2 154 8 1 04 0 5 0 0
21. clared variables Allocating nodes HH Graph Size 41 Initializing model This should take less than a second and you should see the output above in your R console 4 Fitting Bayesian models using R2jags Just like R2WinBUGS the purpose of R2jags is to allow fitting JAGS models from within R and to analyze convergence and perform other diagnostics right within R A typical sequence of using R2jags could look like this See Appendix C in Gelman and Hill 2007 or their online appendix http www stat columbia edu gelman bugsR runningbugs html for more info on how to use R2WinBUGS and 41 Preparing the data and model For an example dataset we simulate our own data in R For this tutorial we aim to fit a linear model so we create a continuous outcome variable y as a function of two predictors x and x2 and a disturbance term e We simulate a dataset with 100 observations e First we create the predictors gt gt gt gt n sim lt 100 set seed 123 xi lt rnorm n sim mean 5 sd 2 x2 lt rbinom n sim size 1 prob 0 3 e lt rnorm n sim mean 0 sd 1 e Next we create the outcome y based on coefficients b and b for the respective predictors and an intercept a gt gt gt gt E 152 b2 lt 3 1 lt 5 y lt a bl xl b2 x2 e e Now we combine the variables into one dataframe for processing later 2 sim dat lt data frame y x1 x2
22. d Jonathan Rapkin 2012 superdiag R Code for Testing Markov Chain Nonconvergence R package version 1 1 URL hittp CRAN R project org package superdiag 24
23. d additional distributions for MCMC models in JAGS Manuscript University of Copenhagen Fernandez i Marin Xavier 2013 ggmcmc Graphical tools for analyzing Markov Chain Monte Carlo simulations from Bayesian inference R package version 0 5 1 URL http xavier fim net packages ggmcmc Gelman Andrew and Jennifer Hill 2007 Data Analysis Using Regression and Multi level Hierarchical Models New York NY Cambridge University Press Jackman Simon 2009 Bayesian Analysis for the Social Sciences Chichester Wiley Kruschke John 2011 Doing Bayesian Data Analysis A Tutorial Introduction with R Oxford Academic Press Elsevier Martin Andrew D Kevin M Quinn and Jong Hee Park 2011 MCMCpack Markov Chain Monte Carlo in R Journal of Statistical Software 42 9 22 Plummer Martyn 2011 JAGS Version 3 1 0 User Manual Plummer Martyn 2013 rjags Bayesian graphical models using MCMC R package version 3 10 URL hittp CRAN R project org package rjags Sturtz Sibylle Uwe Ligges and Andrew Gelman 2005 R2WinBUGS A Package for Running WinBUGS from R Journal of Statistical Software 12 3 1 16 Su Yu Sung and Masanao Yajima 2012 R2jags A Package for Running jags from R R package version 0 03 08 URL hittp CRAN R project org package R2jags Thomas Andrew Bob O Hara Uwe Ligges and Sibylle Sturtz 2006 Making BUGS Open R News 6 1 12 17 Tsai Tsung han Jeff Gill an
24. f alpha 2 396 1 003 820 betal 1 237 1 003 810 beta2 2 692 1 001 3000 deviance 278 737 1 001 3000 17 THE Ht THE For each parameter n eff is crude measure of effective sample size and Rhat is the potential scale reduction factor at convergence Rhat 1 DIC info using the rule pD var deviance 2 pD 4 1 and DIC 275 9 DIC is an estimate of expected predictive error lower deviance is better e Convert to an MCMC object Here we use as mcmc list to keep the 3 chains separate in the resulting MCMC object gt bayes mod fit2 mcmc lt as mcmc list bayes mod fit2 gt summary bayes mod fit2 mcmc Ht Iterations 2001 7000 Thinning interval 1 Number of chains 3 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 alpha 1 829 0 2994 0 0024447 0 012116 betai 1 138 0 0536 0 0004377 0 002189 beta2 3 092 0 2075 0 0016942 0 002772 2 Quantiles for each variable 25b 25 50 97 5 alpha 1 242 1 635 1 830 2 022 2 420 betai 1 032 1 104 1 139 1 173 1 244 beta2 3 496 3 233 3 090 2 950 2 687 e Use mcmcplot and superdiag for diagnostics gt mcmcplot bayes mod fit2 mcmc gt superdiag bayes mod fit2 mcmc burnin 1000 6 Using JAGS via the command line JAGS can also be operated stra
25. ifies more than one chain each chain will start at a different random value for each parameter gt bayes mod inits lt function list alpha rnorm 1 betai rnorm 1 beta2 rnorm 1 gt Alternatively if you want to have control over which starting values are chosen you can provide specific separate starting values for each chain initsi lt list alpha 0 betai 0 beta2 0 inits2 lt list alpha 1 betai 1 beta2 1 inits3 lt list alpha 1 betai 1 beta2 1 bayes mod inits lt list initsi inits2 inits3 gt gt gt gt Before using R2jags the first time you need to load the package and you might need to set a random number seed To do this type gt library R2jags gt set seed 123 directly in R or in your R script You can choose any not too big number here Setting a random seed before fitting a model is also good practice for making your estimates replicable We will discuss replication in more detail in Weeks 3 4 4 2 Fitting the model Fit the model in JAGS gt bayes mod fit lt jags data sim dat jags inits bayes mod inits parameters to save bayes mod params n chains 3 n iter 9000 n burnin 1000 model file Documents Uni 9 ICPSR 2015 Applied Bayes Tutorials 2 JAGS and R bayes mod Compiling model graph Resolving undeclared variables Li Allocating nodes Graph Size 511 Initializ
26. ight from the command line on Windows and Unix systems alike This can be useful if you don t want to have R busy fitting models that take a longer time The most feasible way to do this is to write a script file with the following parts and save it for instance as bayes jags Before running this script you will need to create some files see below model clear data clear load dic model in data in compile inits n inits n inits in bayes mod sim dat nchains 3 bayes mod initsi chain 1 bayes mod inits2 chain 2 bayes mod inits2 chain 3 initialize update 2500 by 100 monitor alpha thin 2 monitor betal thin 2 monitor beta2 thin 2 monitor deviance thin 2 update 2500 by 100 coda The following instructions are for the Terminal on the Mac but you can reproduce the same steps on Windows or Linux You run this script file by opening a Terminal window changing to the working directory in which all the above files are located 18 cd Documents Uni 9 ICPSR 2015 Applied Bayes Tutorials 2 JAGS and and then simply telling JAGS to run the script jags bayes jags In your Terminal you will see something like this eoe angell demo bash 92x18 w ngell demo B Loading modu Loading module Reading data file Compili r file angelll ini Reading parameter file angell2 inits Initializing model Updating 10000 Updating 2500 In more detail this is what
27. ing model Note If you use as your model file the function you gave directly to R above then remove the quotation marks gt bayes mod fit lt jags data sim dat jags inits bayes mod inits parameters to save bayes mod params n chains 3 n iter 9000 n burnin 1000 model file bayes mod Compiling model graph Resolving undeclared variables Li Allocating nodes Graph Size 511 Initializing model Update your model if necessary e g if there is no little convergence gt bayes mod fit upd lt update bayes mod fit n iter 1000 gt bayes mod fit upd lt autojags bayes mod fit This function will auto update until convergence 4 3 Diagnostics One of the neat things about running JAGS or BUGS from R is that this offers seamless and therefore quick access to convergence diagnostics after fitting a model See for yourself gt print bayes mod fit Inference for Bugs model at var folders 8j p9ngyxk5295bnlk0n50mm5ch0000gg T RtmpgkE7yZ model470543d32d1 txt fit using jags iit 3 chains each with 9000 iterations first 1000 discarded n thin 8 n sims 3000 iterations saved mu vect sd vect 2 5 25 50 75 alpha 1 854 0 289 1 281 1 660 1 863 205 betal 1 134 0 052 1 036 1 097 1 134 1 168 beta2 3 092 0 210 3 496 322390 83007 2 943 deviance 271 766 2 871 268 240 269 647 271 079 273 257 RH 97 5 Rhat n eff alpha 2 396 1 003 820 betal 1 237 1 003
28. mod lt function in dN yli dnorm mu i tau mu i lt alpha betai xi i beta2 x2 i vy alpha dnorm 0 01 betai dunif 100 100 beta2 dunif 100 100 tau dgamma 01 01 Uu Now define the vectors of the data matrix for JAGS y lt sim dat y xi lt sim dat x1 x2 lt sim dat x2 gt gt gt gt N lt nrow sim dat Read in the data frame for JAGS gt sim datags lt Ist yl UN You could also do this more conveniently using the as list command on your data frame gt sim dat jags lt as list sim dat Note though that you will also need to specify any other variables not in the data like in this case N So here you would need to add gt sim dat jags N lt nrow sim dat Define the parameters whose posterior distributions you are interested in summarizing later gt bayes mod params lt c alpha betai beta2 Now we need to define the starting values for JAGS Per Gelman and 2007 370 you can use a function to do this This function creates a list that contains one element for each parameter Each parameter then gets assigned a random draw from a normal distributio as a starting value This random draw is created using the rnorm function The first argument of this function is the number of draws If your parameters are not indexed in the model code this argument will be 1 If your jags command below then spec
29. oda stem bayes out Tell JAGS to produce coda files that all begin with the stem bayes out and are put in your WD These can then be analyzed with the tools described in this tutorial 7 Using the coda package You can use the coda package in R to analyze your JAGS output These are the key steps to get your data into R to use CODA e Fit your model in JAGS and identify where your software saved the chains and index files most likely in the working directory where the other components of your model data model inits are e In R load the coda package gt library coda e Read in your JAGS output This requires that the chains and index files see above are in your working directory 20 gt setwd Documents Uni 9 ICPSR 2015 Applied Bayes Tutorials 2 JAGS and R gt chaini lt read coda output file CODAchaini txt index file CODAindex txt Abstracting alpha 1250 valid values Abstracting betal 1250 valid values Abstracting beta2 1250 valid values Abstracting deviance 1250 valid values gt chain2 lt read coda output file CODAchain2 txt index file CODAindex txt Abstracting alpha 1250 valid values Abstracting betal 1250 valid values Abstracting beta2 1250 valid values Abstracting deviance 1250 valid values gt chain3 lt read coda output file CODAchain3 txt index file CODAindex txt Abstracting alph
30. ormat sim list Specify initial values 1 betal 1 beta2 1 O betai O beta2 0 1 betal 1 beta2 1 gt initsi lt list alpha gt inits2 lt list alpha gt inits3 lt list alpha Fit the model in using run jags gt bayes mod fit2 lt run jags model bayes mod monitor c alpha betai beta2 an data sim dat runjags n chains 3 inits list initsi inits2 inits3 burnin 1000 sample 5000 keep jags files TRUE 16 Calling the simulation Welcome to JAGS 3 4 0 on Wed May 27 19 07 26 2015 JAGS is free software and comes with ABSOLUTELY NO WARRANTY Loading module basemod ok Loading module bugs ok it Reading data file data txt Compiling model graph Resolving undeclared variables Allocating nodes Graph Size 511 Reading parameter file initsi txt Reading parameter file inits2 txt Reading parameter file inits3 txt Initializing model Adapting 1000 UAU RRR OS REET SEIS RRS BRS SS 1000 HE 100 Adaptation successful Updating 1000 tM aad 1000 FER 100 Updating 5000 Ui 5000 FEE o aooaa ooo oo ooo 100 Updating 0 it Deleting model d Simulation complete Reading coda files
31. r t value Pr gt t Intercept 1 84949 0 28810 6 42 5 04e 09 1 151351 0 05158 22 00 lt 2e 16 x2 3 09361 0 20650 14 98 lt 2e 16 Signif codes HO eee SO TOTO To OB de Oal UD 5 Residual standard error 0 9367 on 97 degrees of freedom Multiple R squared 0 8772 Adjusted R squared 0 8747 F statistic 346 5 on 2 and 97 DF p value lt 2 2e 16 MCMCpack only requires us to specify one single model object For the linear model we use MCMCregress For other models available in MCMCpack see help package MCMCpack 22 gt bayes mod fit3 lt MCMCregress y xi x2 data sim dat burnin 1000 meme 5000 seed 193 beta start c 0 O 0 te 0 0 0 BO COS OUO gt summary bayes mod fit THE Length Class Mode model 8 jags list BUGSoutput 24 bugs list parameters to save 4 none character model file 1 none character n iter 1 none numeric DIC 1 none logical The resulting object is already an MCMC object 8 1 Output and diagnostics MCMCpack also allows you to use the same set of commands for diagnostics and accessing the fitted objects as you explored above For example you may use mcmcplot and superdiag for diagnostics gt mcmcplot bayes mod fit3 gt superdiag bayes mod fit3 burnin 1000 9 Following this course using JAGS You can find modified if necessary JAGS code for all models presented in this
32. s frontends for JAGS These packages make it easy to process the output of Bayesian models and present it in publication ready form In this tutorial I focus on the R2jags and runjags packages rjags is another package for the same purpose 2 What are JAGS R2jags JAGS 2011 is Just Another Gibbs Sampler that was mainly written by Martyn Plummer in order to provide a BUGS engine for Unix More information can be found in the excellent JAGS manual at http sourceforge net projects mcmc jags R2jags Su and Yajima 2012 is an R package that allows fitting JAGS models from within R Al most all examples in Gelman and Hill s Data Analysis Using Regression and Multilevel Hierarchical Models 2007 can be worked through equivalently in JAGS using R2jags rjags 2013 is another R package that allows fitting JAGS models from within R R2jags depends on it Simon Jackman s Bayesian Analysis for the Social Sciences 2009 provides many examples using rjags and so does John Kruschke s Doing Bayesian Data Analysis 2011 runjags allows some additional functionalities including parallel com puting In this tutorial I focus on the use of R2jags and runjags as well as using JAGS directly from the Terminal 3 Installing JAGS and R2jags 1 Install the most recent R version from the CRAN website http cran r project org bin macosx or http cran r project org bin windows 2 Mac users Install the GNU Fortran gfortran 4 2
Download Pdf Manuals
Related Search