R functions to filter rjags results

June 7, 2012
By

(This article was first published on John Baumgartner's Research » R, and kindly contributed to R-bloggers)

A while back I was running a bunch of JAGS models through R, using the rjags (written by Martyn Plummer) and R2jags (by Yu-Sung Su) packages. These packages provide a great interface to the JAGS software, which allows analysis of Bayesian models (written in the BUGS language) through Markov chain Monte Carlo simulation.

Running a JAGS model using these tools returns an rjags object, which when printed to the screen, summarises the posterior distribution of each monitored node, giving its mean and standard deviation, a range of quantiles, and its Gelman-Rubin convergence diagnostic statistic (Rhat), which indicates the ratio of variance within chains to that among chains. The summary is great, but when monitoring a large number of nodes, printing these to the screen can cause R to hang, and can exceed the screen buffer (not to mention making it painful to find the nodes you’re immediately interested in).

To help deal with this I wrote a couple of simple R functions:

  • jagsresults: return a matrix containing summary results for just the nodes you are interested in (using regular expression pattern-matching, if desired).
  • rhats: sort the output by the nodes’ Rhat values, making it easy to show the n least converged nodes.

Examples

Consider an rjags object that has the following summary (this is the result of the example model given in ?jagsresults). Note that in the examples below I’ve rounded the output to 3 decimal places to prevent line-wrapping/scrollbars on my website (although I’ve excluded the call to round()).

> jagsfit
Inference for Bugs model at "C:/Users/John/AppData/Local/Temp/RtmpIDmQOb/model112833e66869.txt", fit using jags,
 3 chains, each with 10000 iterations (first 5000 discarded), n.thin = 5
 n.sims = 3000 iterations saved
           mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
a            0.041   0.052  -0.057   0.005   0.041   0.076   0.146 1.001  3000
beta.rain    0.997   0.055   0.888   0.960   0.999   1.034   1.105 1.001  3000
beta.temp    1.414   0.054   1.310   1.377   1.415   1.451   1.519 1.001  2400
beta.wind   -0.458   0.057  -0.570  -0.494  -0.458  -0.420  -0.342 1.001  3000
resid[1]    -0.090   0.028  -0.144  -0.109  -0.090  -0.071  -0.035 1.001  3000
resid[2]    -0.048   0.024  -0.094  -0.064  -0.048  -0.032   0.001 1.001  3000
resid[3]     0.076   0.031   0.015   0.055   0.075   0.097   0.138 1.001  3000
resid[4]     0.170   0.050   0.069   0.137   0.171   0.204   0.267 1.001  3000
resid[5]     0.110   0.028   0.054   0.091   0.109   0.129   0.165 1.001  3000
resid[6]     0.086   0.045   0.001   0.055   0.086   0.115   0.175 1.001  3000
resid[7]    -0.016   0.025  -0.066  -0.033  -0.017   0.001   0.033 1.001  3000
resid[8]    -0.226   0.027  -0.280  -0.243  -0.226  -0.207  -0.174 1.001  3000
resid[9]     0.250   0.023   0.204   0.235   0.250   0.266   0.294 1.001  3000
resid[10]   -0.017   0.022  -0.060  -0.032  -0.017  -0.002   0.027 1.001  3000
resid[11]   -0.069   0.025  -0.117  -0.085  -0.069  -0.052  -0.020 1.001  3000
resid[12]    0.145   0.022   0.102   0.131   0.145   0.159   0.187 1.001  3000
resid[13]   -0.184   0.030  -0.246  -0.204  -0.184  -0.163  -0.128 1.001  3000
resid[14]   -0.027   0.029  -0.085  -0.046  -0.027  -0.007   0.030 1.001  3000
resid[15]   -0.015   0.027  -0.067  -0.032  -0.015   0.003   0.038 1.001  3000
resid[16]   -0.176   0.028  -0.231  -0.195  -0.176  -0.158  -0.122 1.001  3000
resid[17]   -0.083   0.032  -0.146  -0.105  -0.083  -0.062  -0.021 1.001  3000
resid[18]   -0.010   0.030  -0.069  -0.030  -0.010   0.010   0.046 1.001  3000
resid[19]   -0.122   0.030  -0.180  -0.142  -0.122  -0.101  -0.063 1.001  3000
resid[20]    0.111   0.024   0.064   0.096   0.111   0.127   0.158 1.001  3000
resid[21]   -0.009   0.039  -0.085  -0.036  -0.010   0.018   0.070 1.001  3000
resid[22]   -0.171   0.026  -0.223  -0.189  -0.171  -0.154  -0.118 1.001  3000
resid[23]   -0.063   0.043  -0.147  -0.092  -0.064  -0.035   0.024 1.001  2800
resid[24]    0.069   0.037  -0.002   0.044   0.070   0.094   0.140 1.001  3000
resid[25]    0.061   0.029   0.003   0.042   0.062   0.081   0.119 1.001  3000
resid[26]    0.046   0.031  -0.017   0.025   0.045   0.066   0.106 1.001  3000
resid[27]    0.138   0.031   0.074   0.117   0.138   0.160   0.197 1.001  3000
resid[28]    0.116   0.020   0.076   0.102   0.116   0.129   0.155 1.001  3000
resid[29]   -0.069   0.035  -0.141  -0.093  -0.069  -0.045  -0.002 1.001  3000
resid[30]   -0.238   0.029  -0.294  -0.257  -0.238  -0.218  -0.183 1.001  3000
resid[31]    0.131   0.036   0.063   0.106   0.131   0.156   0.203 1.001  3000
resid[32]    0.024   0.036  -0.046   0.000   0.024   0.049   0.094 1.001  3000
resid[33]    0.073   0.036   0.002   0.049   0.073   0.097   0.140 1.001  3000
resid[34]   -0.053   0.036  -0.123  -0.077  -0.053  -0.029   0.019 1.001  2700
resid[35]    0.106   0.030   0.047   0.086   0.106   0.126   0.164 1.001  3000
resid[36]    0.154   0.028   0.102   0.136   0.154   0.173   0.211 1.001  3000
resid[37]    0.166   0.029   0.111   0.147   0.166   0.186   0.223 1.001  2100
resid[38]   -0.045   0.043  -0.129  -0.074  -0.045  -0.016   0.040 1.002  2700
resid[39]   -0.076   0.036  -0.149  -0.101  -0.076  -0.052  -0.005 1.001  3000
resid[40]    0.181   0.032   0.117   0.159   0.181   0.203   0.243 1.001  3000
resid[41]    0.135   0.051   0.033   0.101   0.136   0.169   0.233 1.001  3000
resid[42]    0.118   0.026   0.068   0.101   0.119   0.136   0.170 1.001  3000
resid[43]    0.026   0.031  -0.034   0.005   0.025   0.046   0.085 1.001  3000
resid[44]   -0.059   0.032  -0.120  -0.082  -0.059  -0.038   0.005 1.001  3000
resid[45]   -0.035   0.041  -0.117  -0.062  -0.035  -0.007   0.043 1.001  3000
resid[46]    0.058   0.035  -0.013   0.035   0.058   0.082   0.126 1.001  3000
resid[47]    0.035   0.033  -0.030   0.013   0.036   0.057   0.100 1.001  3000
resid[48]    0.087   0.024   0.041   0.070   0.086   0.103   0.135 1.001  3000
resid[49]    0.070   0.026   0.017   0.053   0.070   0.087   0.119 1.001  3000
resid[50]    0.085   0.042   0.001   0.058   0.086   0.113   0.168 1.001  3000
resid[51]    0.018   0.034  -0.047  -0.005   0.018   0.040   0.085 1.002  2000
resid[52]   -0.125   0.044  -0.210  -0.154  -0.125  -0.097  -0.040 1.001  2900
resid[53]    0.443   0.030   0.384   0.423   0.443   0.463   0.502 1.001  3000
resid[54]   -0.351   0.036  -0.421  -0.375  -0.351  -0.326  -0.281 1.001  3000
resid[55]   -0.204   0.021  -0.244  -0.218  -0.204  -0.190  -0.162 1.001  3000
resid[56]   -0.079   0.027  -0.131  -0.098  -0.079  -0.062  -0.027 1.001  3000
resid[57]    0.025   0.029  -0.034   0.005   0.025   0.044   0.081 1.000  3000
resid[58]   -0.149   0.043  -0.233  -0.177  -0.148  -0.120  -0.064 1.001  3000
resid[59]    0.058   0.041  -0.020   0.031   0.057   0.085   0.136 1.001  2200
resid[60]    0.050   0.019   0.012   0.037   0.050   0.063   0.088 1.001  3000
resid[61]    0.237   0.019   0.198   0.223   0.237   0.249   0.275 1.001  3000
resid[62]   -0.062   0.031  -0.121  -0.083  -0.062  -0.041  -0.004 1.001  3000
resid[63]   -0.196   0.030  -0.257  -0.216  -0.197  -0.175  -0.139 1.001  3000
resid[64]    0.009   0.032  -0.055  -0.012   0.010   0.031   0.073 1.001  3000
resid[65]    0.138   0.041   0.057   0.109   0.137   0.165   0.220 1.001  3000
resid[66]    0.124   0.030   0.066   0.103   0.124   0.144   0.183 1.002  1800
resid[67]    0.217   0.022   0.175   0.202   0.217   0.232   0.262 1.001  3000
resid[68]    0.043   0.036  -0.028   0.018   0.043   0.067   0.114 1.001  3000
resid[69]   -0.224   0.034  -0.290  -0.246  -0.224  -0.200  -0.158 1.001  3000
resid[70]    0.072   0.028   0.019   0.053   0.072   0.090   0.128 1.001  3000
resid[71]   -0.184   0.028  -0.238  -0.203  -0.184  -0.165  -0.129 1.001  3000
resid[72]    0.183   0.033   0.117   0.161   0.183   0.205   0.250 1.003  1600
resid[73]    0.041   0.028  -0.015   0.023   0.041   0.061   0.094 1.001  3000
resid[74]    0.131   0.025   0.079   0.114   0.131   0.147   0.180 1.001  3000
resid[75]   -0.138   0.022  -0.182  -0.153  -0.138  -0.123  -0.095 1.001  3000
resid[76]    0.341   0.023   0.295   0.326   0.341   0.357   0.387 1.001  3000
resid[77]   -0.394   0.045  -0.481  -0.424  -0.394  -0.365  -0.306 1.001  2700
resid[78]   -0.306   0.027  -0.359  -0.324  -0.306  -0.287  -0.252 1.001  3000
resid[79]    0.041   0.035  -0.029   0.017   0.041   0.065   0.110 1.001  3000
resid[80]   -0.295   0.030  -0.356  -0.315  -0.294  -0.274  -0.239 1.001  3000
resid[81]    0.159   0.036   0.087   0.135   0.160   0.183   0.226 1.002  3000
resid[82]    0.111   0.033   0.044   0.089   0.110   0.133   0.179 1.002  2000
resid[83]   -0.075   0.025  -0.125  -0.092  -0.075  -0.058  -0.028 1.001  3000
resid[84]   -0.336   0.031  -0.397  -0.356  -0.335  -0.313  -0.275 1.001  3000
resid[85]    0.203   0.019   0.167   0.191   0.203   0.216   0.241 1.001  3000
resid[86]    0.171   0.020   0.132   0.157   0.171   0.185   0.211 1.001  3000
resid[87]   -0.021   0.032  -0.083  -0.042  -0.021  -0.001   0.041 1.001  3000
resid[88]   -0.112   0.031  -0.173  -0.133  -0.112  -0.091  -0.052 1.001  3000
resid[89]    0.012   0.019  -0.026  -0.001   0.012   0.025   0.048 1.001  3000
resid[90]   -0.025   0.031  -0.086  -0.046  -0.025  -0.004   0.034 1.001  3000
resid[91]    0.194   0.037   0.119   0.170   0.194   0.218   0.264 1.001  3000
resid[92]   -0.060   0.036  -0.133  -0.084  -0.060  -0.035   0.009 1.001  3000
resid[93]   -0.073   0.022  -0.117  -0.088  -0.072  -0.058  -0.028 1.001  3000
resid[94]    0.180   0.020   0.140   0.167   0.180   0.193   0.220 1.001  2900
resid[95]    0.050   0.031  -0.012   0.029   0.050   0.071   0.110 1.001  3000
resid[96]   -0.246   0.023  -0.292  -0.261  -0.246  -0.230  -0.202 1.001  3000
resid[97]    0.020   0.037  -0.053  -0.004   0.021   0.045   0.093 1.001  3000
resid[98]   -0.298   0.035  -0.368  -0.322  -0.298  -0.274  -0.229 1.001  3000
resid[99]   -0.038   0.025  -0.089  -0.055  -0.038  -0.021   0.011 1.001  3000
resid[100]  -0.209   0.033  -0.273  -0.231  -0.209  -0.187  -0.145 1.001  3000
sd           0.159   0.012   0.138   0.151   0.158   0.166   0.183 1.001  2900
deviance   -85.830   3.245 -90.105 -88.214 -86.519 -84.102 -78.114 1.001  2500

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 = 5.3 and DIC = -80.6
DIC is an estimate of expected predictive error (lower deviance is better).

The jagsresults() function can be used to return a matrix containing just a and beta.wind:

> jagsresults(x=jagsfit, params=c('a', 'beta.wind'))
            mean    sd   2.5%    25%    50%    75%  97.5%  Rhat n.eff
a          0.041 0.052 -0.057  0.005  0.041  0.076  0.146 1.001  3000
beta.wind -0.458 0.057 -0.570 -0.494 -0.458 -0.420 -0.342 1.001  3000

Or, by specifying invert=TRUE, we can return a matrix containing all parameters except for resid:

> jagsresults(x=jagsfit, params='resid', invert=TRUE)
             mean    sd    2.5%     25%     50%     75%   97.5%  Rhat n.eff
a           0.041 0.052  -0.057   0.005   0.041   0.076   0.146 1.001  3000
beta.rain   0.997 0.055   0.888   0.960   0.999   1.034   1.105 1.001  3000
beta.temp   1.414 0.054   1.310   1.377   1.415   1.451   1.519 1.001  2400
beta.wind  -0.458 0.057  -0.570  -0.494  -0.458  -0.420  -0.342 1.001  3000
deviance  -85.830 3.245 -90.105 -88.214 -86.519 -84.102 -78.114 1.001  2500
sd          0.159 0.012   0.138   0.151   0.158   0.166   0.183 1.001  2900

The argument exact can be set to FALSE to return a matrix containing results for all parameters whose names contain (rather than match exactly) the strings given in params:

> jagsresults(x=jagsfit, params='beta', exact=FALSE)
            mean    sd   2.5%    25%    50%    75%  97.5%  Rhat n.eff
beta.rain  0.997 0.055  0.888  0.960  0.999  1.034  1.105 1.001  3000
beta.temp  1.414 0.054  1.310  1.377  1.415  1.451  1.519 1.001  2400
beta.wind -0.458 0.057 -0.570 -0.494 -0.458 -0.420 -0.342 1.001  3000

Finally, regex=TRUE can be specified to use regular expression pattern-matching, where param is now the pattern to be matched. In this case, additional arguments accepted by the grep function can be supplied to jagsresults, such as perl=TRUE, which is required for the syntax of some patterns:

> jagsresults(x=jagsfit, param='^(?!res)', regex=TRUE, perl=TRUE)
             mean    sd    2.5%     25%     50%     75%   97.5%  Rhat n.eff
a           0.041 0.052  -0.057   0.005   0.041   0.076   0.146 1.001  3000
beta.rain   0.997 0.055   0.888   0.960   0.999   1.034   1.105 1.001  3000
beta.temp   1.414 0.054   1.310   1.377   1.415   1.451   1.519 1.001  2400
beta.wind  -0.458 0.057  -0.570  -0.494  -0.458  -0.420  -0.342 1.001  3000
deviance  -85.830 3.245 -90.105 -88.214 -86.519 -84.102 -78.114 1.001  2500
sd          0.159 0.012   0.138   0.151   0.158   0.166   0.183 1.001  2900

Note, though, that the above (which returns all parameters except for resid) can be more simply produced by the earlier example demonstrating the use of the invert argument.

Standard matrix subsetting is possible for the resulting object, e.g.:

> jags.resids <- jagsresults(x=jagsfit, param='resid')
> jags.resids[c(1, 3, 10:15), c('mean', '2.5%', '97.5%')]
                 mean        2.5%       97.5%
resid[1]  -0.08970919 -0.14412395 -0.03469326
resid[3]   0.07606610  0.01462996  0.13840041
resid[10] -0.01695619 -0.06031680  0.02725152
resid[11] -0.06851323 -0.11697671 -0.02024998
resid[12]  0.14516502  0.10249481  0.18689894
resid[13] -0.18398696 -0.24553446 -0.12795644
resid[14] -0.02688139 -0.08512743  0.02994298
resid[15] -0.01459422 -0.06701636  0.03812518

Installation

The simplest way to install the jagstools package is with install_github in the devtools package, which will grab the source directly from Github:

install.packages('devtools')
library(devtools)
install_github(repo='jagstools', username='johnbaums')

The package can also be downloaded from github.

Filed under: R Tagged: Bayesian, code, JAGS, MCMC, model, R, R2jags, rjags, rstats

To leave a comment for the author, please follow the link and comment on their blog: John Baumgartner's Research » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers


Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)