Chain convergence on model

[Sebastian Evans]

Hi all,

I am currently trying to model species richness using several environmental predictors, the most important among which is distance from a retreating glacier (a proxy of time since deglaciation). For this model, I am having serious difficulty getting convergence among my environmental covariates, even when they are centered and scaled. So much so that Rhat is "INF" in my MCMC summary output. I wondered if anyone had any quick fixes on this front, here is the code I am working with and can supply the data set if needed. This analysis is part of my Graduate Thesis. Thanks!

#### Model1B, Species Richness: Random effect for transect, distance and other covariates ####

Model1B_code <- nimbleCode({
 
  # Priors
  beta0 ~ dnorm(0, .00001)   #vague prior distribution for beta0
  beta1 ~ dnorm(0, .00001)   #vague prior distribution for beta1
 
  sd ~ dunif(0, 100)    #vague prior distribution for sd of epsilon
  ep.sig <- 1 / (sd*sd)      #convert sd to precision (nimble default is precision)
 
  sd2 ~ dunif(0, 100)    #vague prior distribution for sd of alpha
  t.sig <- 1 / (sd2*sd2)    #convert sd2 to precision (nimble default is precision)
 
 
  # Likelihood:
  for (i in 1:n){
    rich[i] ~ dnorm(mu[i], ep.sig)
    mu[i] <- beta0 + beta1*dist[i] + beta2*twi[i] + beta3*tri[i] + beta4*gravfines[i] + alpha[t[i]]
   
    #For GOF/Bayesian p-value
    rich.sim[i] ~ dnorm(mu[i], ep.sig)
    dev[i] <- -2*(-log(sd)-.5*log(6.283)- ((rich[i]-mu[i])^2)/(2*sd*sd))
    dev.sim[i] <- -2*(-log(sd)-.5*log(6.283)- ((rich.sim[i]-mu[i])^2)/(2*sd*sd))
  }
  for (j in 1:nSites){
    alpha[j] ~ dnorm(0, t.sig)
  }
 
  dev2 <- sum(dev[1:n])
  dev2.sim <- sum(dev.sim[1:n])
 
})


nSites = length(unique(model_data$transect))

# Bundle data
data.feed <- list(rich = as.numeric(model_data$richness),
                  n = as.numeric(length(model_data)),
                  dist = as.numeric(model_data$distance),
                  t = as.numeric(model_data$transect),
                  twi = as.numeric(model_data$twi),
                  tri = as.numeric(model_data$tri),
                  gravfines = as.numeric(model_data$grav_fines),
                  nSites=as.numeric(nSites))


# Initial values
inits <- function(){
  list(beta0 = runif(1, -2, 2),
       beta1 = runif(1, -3, 3),
       beta2 = runif(1, -3, 3),
       beta3 = runif(1, -3, 3),
       beta4 = runif(1, -3, 3),
       sd = runif(1, 0,1),
       sd2 = runif(1, 0,1),
       alpha=rnorm(nSites, 0, 2))
}

# Parameters monitored
params <- c("beta0", "beta1", "beta2", "beta3", "beta4", "sd", "sd2", "dev2", "dev2.sim")

# MCMC settings
ni <- 575000  
nt <- 5  
nb <- 560000  
nc <- 3


# Call nimble from R
start <- Sys.time()
Model1B <- nimbleMCMC(
  code = Model1B_code,
  constants = data.feed,
  inits = inits,
  monitors = params,
  summary = TRUE,
  WAIC = FALSE,
  nchains = nc,
  thin = nt,
  niter = ni,
  nburnin = nb,
  progressBar = TRUE)
Sys.time() - start

# Assess MCMC chain convergence
# Assess MCMC chain convergence
chainsPlot(Model1B$samples, var=c("beta0", "beta1", "sd", "sd2"))
chainsPlot(Model1B$samples, var=c("dev2", "dev2.sim"))
MCMCsummary(Model1B$samples, round = 2)
Model1B_stacked <- rbind(Model1A$samples$chain1,Model1A$samples$chain2,Model1A$samples$chain3)
mean(Model1B_stacked[,"beta1"] > 0)

[Chris Paciorek]

Hi Sebastian,

That's a pretty straightforward model, so I wouldn't think you should be having such poor performance. And running for 560000/575000 iterations seems like much more than should be needed.

You should have seen a message suggesting that you should put `t` in the constants list. Otherwise, nimble's mcmc calculations don't "know" that `t` never changes and a bunch more calculation needs to be done (though I wouldn't think that even that would fully explain what is happening).

So please try putting `t` in constants and let me know how much that changes things. 

Also, if you didn't see a message about putting `t` in constants it might be that the message is hidden somehow by use of `nimbleMCMC`. We'd want to know about that to fix it, so please let us know if you don't see that message when running your current code (before the change I recommended above).

In addition, which of the parameters mix poorly? As a vague guess, it might be that `beta0` trades off with the mean of the `alpha` random effects. So you could try reparameterizing to have `alpha[j] ~ dnorm(beta0, t.sig)` and remove `beta0` from the linear predictor. 

-chris

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/nimble-users/f11a92aa-cf83-462f-8821-cd78533882c8n%40googlegroups.com.