Addressing R-hat=NA

11 views

Skip to first unread message

Trinh Dinh Hoang

unread,

Nov 5, 2025, 3:01:15 PM (2 days ago) Nov 5

to hmecology: Hierarchical Modeling in Ecology

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

unread,

Nov 5, 2025, 3:33:47 PM (2 days ago) Nov 5

to Trinh Dinh Hoang, hmecology: Hierarchical Modeling in Ecology

Can you attach some traceplots of the problematic parameters?

--
*** 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.

Reply all

Reply to author

Forward

0 new messages