Never blindly trust Rhat for convergence in a Bayesian analysis !

[Kery Marc]

Dear all,

introducing a new "Cautionary Picture" feature: here is a particularly striking example of how Rhat can sometimes tell you an MCMC algorithm has converged ..... but it hasn't ! OK, this may be a pathological example, but nevertheless it's *really* good practice to always give a quick visual check to the chains for at least some key parameters in a model fit with MCMC.

Best regards  --- Marc

[kkellner]

Hi Marc,

Split-chain R-hat should detect this issue. It was introduced somewhere in Gelman's Bayesian Data Analysis book. Basically you take each MCMC chain and split it in half by iteration #, and treat each half as a separate chain. In your example you'd then probably have a large contrast between the 3 chains that came from the beginning of the original chains, and 3 chains from the end, which would result in a very large between-chain variance and a large Rhat.

This is the approach now used by rstan, but you can actually apply the relevant rstan function (rstan::monitor) to any MCMC output after a bit of formatting. For example, with jagsUI output:

library(jagsUI)

library(rstan)

data(longley)

head(longley)

gnp <- longley$GNP

employed <- longley$Employed

n <- length(employed)

data <- list(gnp=gnp,employed=employed,n=n)

modfile <- tempfile()

#Write model to file

writeLines("

model{

  #Likelihood

  for (i in 1:n){

    employed[i] ~ dnorm(mu[i], tau)

    mu[i] <- alpha + beta*gnp[i]

  }

  #Priors

  alpha ~ dnorm(0, 0.00001)

  beta ~ dnorm(0, 0.00001)

  sigma ~ dunif(0,1000)

  tau <- pow(sigma,-2)

}

", con=modfile)

params <- c('alpha','beta','sigma')

out <- jags(data = data,

            inits = NULL,

            parameters.to.save = params,

            model.file = modfile,

            n.chains = 3,

            n.adapt = 100,

            n.iter = 1000,

            n.burnin = 500,

            n.thin = 2)

#Convert rjags mcmc.list to array format (iterations x parameters x chains)

arr <- as.array(out$samples)

#convert array to iterations x chains x parameters

arr <- aperm(arr, c(1,3,2))

#Calculate split-chain Rhat and other parameters

#no warmup/burn-in iterations saved by jagsUI, so set warmup to 0

mon <- rstan::monitor(arr, warmup=0)

I'd like to implement this in jagsUI at some point.

Ken

[Denis Valle]

The idea behind split chain Rhat seems very similar to Geweke's diagnostic test...

https://www.math.arizona.edu/~piegorsch/675/GewekeDiagnostics.pdf

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2020)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/d998f09e-d862-4434-b827-2ae4ab3794f3o%40googlegroups.com.

[Kery Marc]

Dear Ken,

thanks for the extremely informative and clear explanation (as usual)

Best regards  --- Marc

---

From: hmec...@googlegroups.com [hmec...@googlegroups.com] on behalf of kkellner [ken.k...@gmail.com]
Sent: 28 October 2020 17:11
To: hmecology: Hierarchical Modeling in Ecology
Subject: Re: Never blindly trust Rhat for convergence in a Bayesian analysis !

--

[Tomas T.]

Hello all!

thanks Marc for the example! And what about psrf (potential scale
reduction factor), as reported by gelman.diag() function from the R
package coda? I am using this instead of Rhat. Do you (or anyone else)
know if this statistics can suffer of similar problem or is it more
reliable?

Tomas

> --
> *** Three hierarchical modeling email lists ***
> (1) unmarked: for questions specific to the R package unmarked
> (2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
> (3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2020)
> ---
> You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/F801018864E940468323F405DF95D2DD76110FEE%40MailNT.vogelwarte.ch.

[kkellner]

Hi Tomas,

"Rhat" and "potential scale reduction factor" are two names for the same statistic, so they will have similar problems though they may be calculated slightly differently in different packages. 

The Rhats reported in Marc's example were calculated using coda::gelman.diag, assuming they came from jagsUI.

Ken

On Wednesday, November 11, 2020 at 8:25:57 AM UTC-5, Tomas T. wrote:

Hello all!

thanks Marc for the example! And what about psrf (potential scale
reduction factor), as reported by gelman.diag() function from the R
package coda? I am using this instead of Rhat. Do you (or anyone else)
know if this statistics can suffer of similar problem or is it more
reliable?

Tomas


On Mon, 26 Oct 2020 at 18:26, Kery Marc <marc...@vogelwarte.ch> wrote:
>
> Dear all,
>
> introducing a new "Cautionary Picture" feature: here is a particularly striking example of how Rhat can sometimes tell you an MCMC algorithm has converged ..... but it hasn't ! OK, this may be a pathological example, but nevertheless it's *really* good practice to always give a quick visual check to the chains for at least some key parameters in a model fit with MCMC.
>
> Best regards  --- Marc
>
> --
> *** Three hierarchical modeling email lists ***
> (1) unmarked: for questions specific to the R package unmarked
> (2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
> (3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2020)
> ---
> You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to hmec...@googlegroups.com.