2

I am using "Multivariate PSRF" statistics from gelman.diag() function to analyze my MCMC chains. Now I analyzed convergence 471 variables (parameters for each site).

  1. How is it possible that "Multivariate PSRF" says 1.248096, when psrf for each of 471 sites ranges from 0.9996 to 1.012 max? Psrf CI is 1.044 max!
  2. It seems Multivariate PSRF grows with the number of variables, which is definitely not desired of convergence assessment.

I think I can safely say that the model converged, but the Multivariate PSRF suggests otherwise... so is it computed wrong or something?

Tomas
  • 5,735
  • 11
  • 52
  • 93

2 Answers2

5

This is just a multiple comparisons problem. The multivariate PSRF is derived from the linear combination of variables that has the maximum scalar PSRF. With 471 variables, there is a good chance to find an "optimal" linear combination that has high PSRF just due to random fluctuations in the data.

You should feel free to ignore it in this case

Martyn Plummer
  • 211
  • 1
  • 1
  • Thanks Martyn. Then there is a question if this multivariate psrf shouldn't perhaps be computed differently to be useful? – Tomas Feb 25 '15 at 14:53
2

As Martyn says, I would not be too worried about the MPSRF in this case. I don't think that it is being computed incorrectly, but it occasionally gives false positives and false negatives. This is true of all convergence criteria, which is why it is important to always visualise trace plots (at least) for all key parameters in your model (NOT just the ones that you are interested in). For example, consider the following:

  library(coda)
  set.seed(1)
  chain1 <- mcmc(rnorm(100, 1:100, 10))
  chain2 <- mcmc(rnorm(100, 100:1, 10))

  gelman.diag(mcmc.list(chain1, chain2), autoburnin=FALSE)
  lattice::xyplot(mcmc.list(chain1, chain2), autoburinin=FALSE)

enter image description here

This is obviously a pathological example, but it illustrates my point :)

Matt

Tomas
  • 5,735
  • 11
  • 52
  • 93
Matt Denwood
  • 906
  • 6
  • 5
  • Thanks Matt! But 1) I have hundreds of models, cannot afford to visualise all chains. I definitely need some kind of statistics to do that for me 2) as for the *"NOT just the ones that you are interested in"*, I don't agree, I do exactly this. I have hierarchical GLM model with correlated predictors, the coefficients of particular predictors sometimes do not converge. But I am interested in predicted response at given sites which converges perfectly. – Tomas Feb 26 '15 at 09:12
  • (1) If you don't have time to visualise everything, then there will always be a risk (however small) that you have missed false-convergence somewhere. You could mitigate this by also using e.g. Geweke's diagnostic as well, but you will get more false positives that way... (2) I disagree with this - I think that it is dangerous to assume that a parameter has converged when another (especially influential) parameter has not. I have had a few cases where one parameter converged on a local optimum for some time before eventually finding the global optimum AFTER another parameter had converged. – Matt Denwood Feb 26 '15 at 11:11
  • thanks Matt. (2) this probably depends case by case. I think in my GLM case I can safely ignore divergence of the particular coefficients of the linear predictor *if* the predicted response converges. – Tomas Feb 26 '15 at 13:26