Addressing R-hat=NA
Trinh Dinh Hoang
Dear all,
I am currently running N-mixture models (Binomial/ Poisson; see code below) to estimate population size. However, the results seem unusual — the R-hat values are 0 for all sites with no detections (see below).
Could you please advise me on how to address this issue?
Many thanks,
Hoang
1. Code:
modeltext <- "
model{
for(i in 1:nSites){
# Biological model Poisson distribution
N[i] ~ dpois(lambda[i])
log(lambda[i]) <- alpha.lam + 1 * logA[i] +
eps[i] # we used the log of the size of the sampling unit (logA) as an
additive offset
eps[i] ~ dnorm(0, tau) # random effect
# Observation model: Binomial distribution
for(j in 1:nOcc) {
C[i,j] ~ dbin(p, N[i])
# GOF
FT.obs1[i,j] <- (sqrt(C[i,j]) - sqrt(lambda[i]
* p))^2
C.sim[i,j] ~ dbin(p, N[i])
FT.sim1[i,j] <- (sqrt(C.sim[i,j]) -
sqrt(lambda[i] * p))^2
}
}
# Priors
alpha.lam
~ dunif(-10, 10)
tau <- pow(sigma, -2)
sigma ~ dunif(0, 10)
#psi ~ dbeta(4.375, 5.625)
# centered on 0.4375, naive occupation
p
~ dbeta(5.63, 4.37) # centered on 0.563,
naive average detection probablity
#p ~ dbeta(1,1) # use a
uniform Beta(1,1) prior for the probability of detection P
# GOF
FT.obs
<- sum(FT.obs1)
FT.sim
<- sum(FT.sim1)
# derived variable
totalN <- sum(N[]) #Estimate of total abundance across all sites
mean.abundance <-
mean(lambda[]) #mean expected abundance per occupied
site
}"
2. Results
MCMC values from jagsUI
object ‘jagsOut’
The object has 38 nodes with 3e+05 draws for each of 6 chains.
l95 and u95 are the limits of a 95% Highest Density Credible Interval.
Rhat is the estimated potential scale reduction factor:
largest is 2.00; 1 (3%) are greater than 1.10; 20
(53%) are NA.
MCEpc is the Monte Carlo standard error as a percentage of the posterior SD:
largest is 0.4%; NONE are greater than 5%.
mean
sd median l95 u95 Rhat MCEpc
totalN 40.107 8.291 38.000 28.000
55.000 1.000 0.349
mean.abundance 1.254 0.326 1.207 0.700 1.893 1.000
0.290
p 0.460 0.082 0.463
0.297 0.616 1.000 0.304
FT.obs 36.081 8.895 34.944 20.530
53.968 1.000 0.353
FT.sim 37.766 8.074 36.889 23.309
54.083 1.000 0.323
N[1] 0.036 0.202 0.000
0.000 0.000 NA 0.144
N[2] 0.027 0.172 0.000
0.000 0.000 NA 0.139
N[3] 0.041 0.215 0.000
0.000 0.000 NA 0.143
N[4] 0.029 0.182 0.000
0.000 0.000 NA 0.141
N[5] 0.034 0.195 0.000
0.000 0.000 NA 0.139
N[6] 0.034 0.196 0.000
0.000 0.000 NA 0.142
N[7] 0.037 0.205 0.000
0.000 0.000 NA 0.146
N[8] 1.059 0.262 1.000
1.000 2.000 NA 0.152
N[9] 0.035 0.198 0.000
0.000 0.000 NA 0.141
N[10] 0.035 0.198 0.000
0.000 0.000 NA 0.140
N[11] 0.035 0.197 0.000
0.000 0.000 NA 0.146
N[12] 0.033 0.194 0.000
0.000 0.000 NA 0.145
N[13] 1.104 0.355 1.000
1.000 2.000 NA 0.160
N[14] 2.297 0.634 2.000
2.000 4.000 1.000 0.186
N[15] 0.033 0.193 0.000
0.000 0.000 NA 0.145
N[16] 0.036 0.201 0.000
0.000 0.000 NA 0.140
N[17] 1.116 0.377 1.000
1.000 2.000 2.000 0.158
N[18] 2.496 0.835 2.000
2.000 4.000 1.091 0.225
N[19] 0.021 0.150 0.000
0.000 0.000 NA 0.138
N[20] 1.210 0.516 1.000
1.000 2.000 1.000 0.175
N[21] 1.390 0.709 1.000
1.000 3.000 1.000 0.194
N[22] 0.028 0.176 0.000
0.000 0.000 NA 0.142
N[23] 0.029 0.179 0.000
0.000 0.000 NA 0.138
N[24] 1.197 0.500 1.000
1.000 2.000 1.000 0.172
N[25] 3.663 1.017 3.000
3.000 6.000 1.000 0.227
N[26] 4.981 1.697 5.000
3.000 8.000 1.000 0.265
N[27] 6.709 2.114 6.000
4.000 11.000 1.000 0.275
N[28] 0.028 0.179 0.000
0.000 0.000 NA 0.148
N[29] 3.751 1.467 3.000
2.000 6.000 1.000 0.255
N[30] 2.599 0.938 2.000
2.000 4.000 1.000 0.219
N[31] 5.955 1.810 6.000
4.000 9.000 1.000 0.264
N[32] 0.030 0.184 0.000
0.000 0.000 NA 0.140
deviance 119.980 6.483 119.037 109.240 132.597 1.000 0.223
#Posterior predictive check
plot
pp.check(jagsOut, observed = 'FT.obs', simulated = 'FT.sim')
# P-value = 0.64
DIC(jagsOut)
# DIC = 141.00
Matthijs Hollanders
--
*** 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, 2021) and Schaub & Kéry (2022)
---
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 visit https://groups.google.com/d/msgid/hmecology/3b0f2985-a486-4693-a79d-e00e9898bc53n%40googlegroups.com.
Trinh Dinh Hoang
Hi,
The issue with NA values has been resolved by using the MCMCvis package instead of mcmcOutput.
However, I’m facing another problem: although all R-hat values are equal to 1, the trace plots (attached) do not indicate good convergence—they show noticeable spikes.
Could you please advise what might be causing this issue and how I could address it?
Thank you very much,
Hoang
Marc Kery
Sent: Tuesday, November 11, 2025 05:29
To: hmecology: Hierarchical Modeling in Ecology <hmec...@googlegroups.com>
Subject: Re: Addressing R-hat=NA
Trinh Dinh Hoang
Hi,
The issue with NA values has been resolved by using the MCMCvis package instead of mcmcOutput.
However, I’m facing another problem: although all R-hat values are equal to 1, the trace plots (attached) do not indicate good convergence—they show noticeable spikes.
Could you please advise what might be causing this issue and how I could address it?
Thank you very much,
Hoang