-1

Extracting stuff from objects has always been one of the most confusing aspects of R to me. I've fitted a bayesian linear regression model using rjags and have the following mcmc object:

summary(m_csim)
Iterations = 1:150000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 150000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

            Mean        SD  Naive SE Time-series SE
BR2     0.995805 0.0007474 1.930e-06      3.527e-06
BR2adj  0.995680 0.0007697 1.987e-06      3.633e-06
b[1]   -5.890842 0.1654755 4.273e-04      1.289e-02
b[2]    1.941420 0.0390239 1.008e-04      1.991e-03
b[3]    1.056599 0.0555885 1.435e-04      5.599e-03
sig2    0.004678 0.0008333 2.152e-06      3.933e-06

2. Quantiles for each variable:

            2.5%       25%       50%       75%    97.5%
BR2     0.994108  0.995365  0.995888  0.996339  0.99702
BR2adj  0.993932  0.995227  0.995765  0.996229  0.99693
b[1]   -6.210425 -6.000299 -5.894810 -5.784082 -5.55138
b[2]    1.867453  1.914485  1.940372  1.967466  2.02041
b[3]    0.942107  1.020846  1.057720  1.094442  1.16385
sig2    0.003321  0.004082  0.004585  0.005168  0.00657

In order to extract the coefficients' means I did b = colMeans(mod_csim)[3:5]. I want to calculate credible intervals so I need to extract the 0.025 and 0.975 quantiles too. How do I do that programmatically ?

MrFlick
  • 195,160
  • 17
  • 277
  • 295
FranGoitia
  • 1,965
  • 3
  • 30
  • 49
  • 2
    Usually there are methods for this. Not sure about `rjags`, but `coef` and `fixef` are generally used, either on the model object or the summary of the model object. There are also good packages, like `tidybayes`. – Axeman Nov 05 '18 at 15:53
  • 1
    It's easier to help you if you include a simple [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input and desired output that can be used to test and verify possible solutions. – MrFlick Nov 05 '18 at 15:53
  • 1
    I would start with `str(m_csim)` to see what that is exactly. Then, remember that [summary functions](https://cran.r-project.org/doc/FAQ/R-FAQ.html#How-should-I-write-summary-methods_003f) many times do their own computations before printing and also run `str(summary(m_csim))`. – Rui Barradas Nov 05 '18 at 15:55

3 Answers3

1

You can probably extract the quantiles directly. As others have pointed out, you can call str(m_csim), and you can limit the output of the str call with str(m_csim, max.level=1), and keep adding one to the max.level= argument until you see something that looks like quantiles.

What I like to do is to convert the MCMC output to a data.frame so it's easier to work with. I use jagsUI rather than rjags, but I often do something like

mcmc_df <- as.data.frame(as.matrix(MY_MCMC_OBJECT$samples))

Note: it might be a little different with rjags, but I'm sure you can find it with a little digging.

The benefit: I can then access a single vector for each mcmc_df$PARAMETER, and create a matrix of quantiles with

mcmc_quants <- apply(mcmc_df, 2, quantile, p=c(0.025, 0.25, 0.5, 0.75, 0.975))

or whatever quantiles you want.

Matt Tyers
  • 2,125
  • 1
  • 14
  • 23
1

You probably are looking for

model_summary_object <- summary(m_csim)
model_summary_object$quantiles[,c('2.5%','97.5%')]
Pepacz
  • 881
  • 3
  • 9
  • 24
0

I hope I'm not overstepping my knowledge bounds, but I want to answer from the point of view of 'in general' rather than specifically for rjags. m_csim is an object and many methods can possibly be used on it. You've used the summary method to see something. As people have commented, probably there is a coef method. However as someone else has commented (while I was replying !), using str() to see what an object contains is the best way to see what information is in an object and how to address it. I'd be very surprised if using str() doesn't show how to find not only the coefficients but also enough information on the confidence intervals to allow you to find the desired CI.

meh
  • 218
  • 1
  • 9