R functions to filter rjags results

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.

Advertisements

9 thoughts on “R functions to filter rjags results

  1. Nice work John, now I just need to extract the elements from a print summary… eg the median and upper and lower quantiles, into an R object for manipulating them.

    • Thanks Pete. The object that holds the results of the jagsresults() function is a matrix containing mean, sd, a bunch of percentiles (2.5, 25, 50, 75, and 97.5), Rhat, and n.eff. These can then be directly accessed/subsetted as per usual. I’ve added an example of this in the post.

  2. Hi John. Thanks for your great work on ”jagstools”. I have a question about “jagstools” package. I’m now running a Bayesian model by R and I have tried both package “rjags” and package “R2jags”. I found that these 2 packages use different functions to create a jags model (“jags.model()” in “rjags” while “jags()” in “R2jags”). When I used “jags()” to create a model, I can successfully summary the results with “jagsresults”. But when used “jags.model()” to create a model, the output of “jagsresults” is “null”. I wonder whether I can use “jagsresults” to summary the results of “jags.model()” in “rjags” package. Can you tell me about that?
    Thanks a lot!

    • Hi Tang. Thanks for your comment. You won’t be able to use it directly on a jags model object, but I’ve modified the function so that it should now also work with mcmc.list objects that are created by the coda.samples function in the rjags package. Workflow would be (1) jags.model(), (2) coda.samples(), (3) jagsresults(). Hope that helps.

      • Thanks John. I have tried the workflow you advised. But it still doesn’t work. I don’t know whether there is any miss in my operation. However, I can use “jagstools” to analyze the results of “R2jags” and I think I will choose that way.
        Thank you for your help!

      • Sorry to hear it’s still not working out for you. I personally prefer using R2jags, but I tested jagsresults on an mcmc.list object (using version 1.2 of jagstools – on Github) and it worked for me. If you like, send me your code or the mcmc.list object and I’ll take a look.

      • Hi John. Thank you so much for your kind help. For some reasons, I have stopped this work for nealy one month. Sorry for late response. I would like to send you my code and data, if you have interest in it. Would you tell me your mailaddress that I can send them to? Actually it is a complicated model with large amount of data, at least to me. I’m a new player in the world of R and jags 🙂 So many things to learn.

  3. Pingback: These are QAEco’s favourite R packages, what’s yours? #favRpackage | The Quantitative & Applied Ecology Group

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s