18 thoughts on “MCMC chain analysis and convergence diagnostics with coda in R

  1. Pingback: A simple Metropolis-Hastings MCMC in R « theoretical ecology

  2. Pingback: Monday, April 30 | Julie I. Santos

  3. Pingback: A simple Approximate Bayesian Computation MCMC (ABC-MCMC) in R « theoretical ecology

  4. Pingback: Explaining the ABC-Rejection Algorithm in R | theoretical ecology

  5. Pingback: Explaining the ABC-Rejection Algorithm in R ← Patient 2 Earn

  6. Pingback: Making the Most of MCMC | Sweet Tea, Science

    • Thanks for broken link, I fixed this. In coda and most other packages, burn-in is not done automatically, so if you don’t specify it the whole chain will be used. It depends on what you want to do with the chain if you discard burn-in (for sure, you want to discard if you want to estimate posterior summaries). In this example, I didn’t discard it so that you can see the full trace / Gelman plots.

  7. Hi Florian,

    I have a question about the section of Marginal densities hide correlations.
    Typically, people talk about the correlation of input variables, not the parameters. I am puzzled that how do we know the parameters are correlated and what’s the high correlation of parameters supposed to mean.

    By running the code, i got a narrower bandwidth for parameter 1 and a wider bandwidth for parameter 2. This is counter-intuitive by reading your post because we are supposed to get a wider bandwidth for both parameters, correct? Can you please shed more lights on it?

    Thank you very much.
    Lizy

    • Hi Lizy,

      when I say that there is a correlation between parameters, what I mean is that we can maintain a good fit / likelihood / posterior by increasing parameter X1 when we decrease parameter X2 at the same time. You could also call this a trade-off between the two parameters.

      Such trade-offs can appear through a number of reasons. ONE important reason is correlations between the predictors (collinearity). Other possible reasons are a) the model structure (as in the example with the intercept), or b) simply that the problem is underdetermined (not enough data to constrain the parameters).

      You may be right that people talk more about collinearity (correlation in predictor space), but identifying trade-offs in parameter space is crucial for a correct interpretation of the uncertainty and for correct forecasting. This is independent of whether you are working Bayesian or not. Also lme4, for example, returns the covariance matrix of the likelihood density (i.e. the correlation of the parameters) as a standard output with their summary() function.

      About the width – I just reran it and yes, the marginal uncertainty for the slope shrinks. I have to think a bit more about why this happens, I find it somewhat counterintuitive. In any case, what we should keep in mind is that the true uncertainty of the combination of slope / intercept is smaller than in appears marginally, because particular combinations of the two are not compatible with the data.

  8. Hi Florian,
    first of all thanks for the great articles on McMC sampling in R. They have been of great help so far also for an R-newbie like me.

    I would like to ask you in it’s possible to use CODA or other packages in order to run parallel independent chains using a multicore system and then combine the chains. Thanks a lot.

    • If you want to use the code from this post, just use the standard parallel functions of R to parallelize a for loop / apply for the chains, collect the results, and put the chains together.

      For the more serious MCMC packages for Bayesian analysis (e.g. Jags, STAN) there are options in the R interface to do this.

  9. Hi Florian,
    Thanks for the post, very helpful. I would like to use your version (betterPairs) of the mcmc pairs plot but I can’t figure out how to add in the code to select a subset of parameters from a stanreg object generated using stan_glmer in the rstanarm package. Any idea how I can add that to the arguments? For example, in pairs() I can achieve this using regex_pars.
    Thanks in advance for any help you can give

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