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.
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:
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 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.