Site icon R-bloggers

R functions to filter rjags results

[This article was first published on John Baumgartner's Research » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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:

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, para, 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, para, 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 about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.