Royle-Nichols models for occupancy

In the simplest model, the number available for detection at site \( i \), \( N_i \), is modelled as being drawn from a Poisson distribution with parameter \( \lambda \). This is the biological model. The observation model assumes that detection of each individual is a Bernoulli trial, individuals being detected with probability \( r \). The data do not show which individuals are detected, just whether the species was detected, and that is recorded if at least one individual is observed.

To get from the probability of detecting one individual, \( r \), to the probability of detecting the species, \( p \), think of the probability of non-detection: The probability of not detecting an individual is \( 1 - r \) and if \( N_i \) individuals are present, the probability of detecting none of them is \( (1 - r)^{N_i} \). Detecting "at least one" is equivalent to detecting "not none", and the probability of that is \( p = 1 - (1 - r)^{N_i} \).

The model parameters to be estimated are \( \lambda \) and \( r \). The probability of occupancy, \( \psi \), is derived from \( \lambda \):

\[ \psi = {\rm I\!P}(N_i > 0) = 1 - e^{-\lambda} \]

As usual, we can model \( \lambda_i \) as a function of site covariates, and \( r_{it} \) as a function of site and survey covariates, but we assume that detection probability is the same for all individuals at site \( i \) and time \( t \) and detections are independent.

Maximum likelihood estimation for the Royle-Nichols occupancy model is available with the occSSrn function in the wiqid package and occuRN in unmarked . It was difficult to make the Bayesian version work in WinBUGS - it could be done by working with \( r' = 1 - r \)  and \( y' = 1 - y \) (ie, \( y' = 1 \) if the species is not detected, and \( y' = 0 \) if it is detected).

Implementation in JAGS

Simulating data

Let's begin by simulating a data set to work with. That process itself will clarify how the model works.

# Parameters for simulated data:
nSite <- 200   # number of sites
nRep <- 10     # number of replicate surveys
lambda <- 1.5  # mean number of individuals at each site
r <- 0.1       # probability of detection of one individual

set.seed(2017)
# Biological process:
N <- rpois(nSite, lambda)  # Latent number available for detection
mean(N) # Mean N for this sample (not same as lambda)
[1]  1.36
mean(N > 0) # True occupancy for this sample
[1]  0.725

# Observation process:
p <- 1 - (1-r)^N
y <- rbinom(nSite, nRep, p) # number of occasions the species was detected at each site
head(y)
[1]  3 1 2 2 2 2
mean(y > 0)  # naive estimate of occupancy 
[1] 0.555

Maximum likelihood estimation

We can now fit the Royle-Nichols model with the occSSrn0 function in wiqid ; occSSrn0 requires only the vector y and the number of replicates at each site, not the full detection/nondetection matrix, and runs quickly:

library(wiqid)
occSSrn0(y, n=rep(nSite, nRep))
Call: occSSrn0(y = y, n = rep(nSite, nRep))

Real values:
             est  lowCI  uppCI
psiHat    0.6830 0.5755 0.7856
lambdaHat 1.1488 0.8569 1.5401
rHat      0.1143 0.0859 0.1506

AIC: 1424.452 

The parameter estimates \( \lambda \) and \( \psi \) are a bit low for this realisation. We can compare with the standard psi(.)p(.) model:

occSS0(y, n=rep(nSite, nRep))
Call: occSS0(y = y, n = rep(nSite, nRep))

Real values:
          est  lowCI  uppCI
psiHat 0.6247 0.5402 0.7023
pHat   0.1969 0.1717 0.2247

AIC: 1431

The AIC value for this model is a lot higher, indicating that the Royle-Nichols model would give better predictions than the standard model.

The JAGS model

Here is the JAGS model I used. Note that in this model z[i] is only used to calculate the proportion of occupied sites in the sample, it does not appear in the calculation of probability of detection.

model {
  # Likelihood
  for (i in 1:nSite) {
    # Process model
    N[i] ~ dpois(lambda)
    z[i] <- step(N[i] - 1) # z=1 if N>0, ie. site is occupied

    # Observation model
    p[i] <- 1 - pow(1-nRep, N[i])
    y[i] ~ dbin(p[i], n)
  }
  # Priors
  lambda ~ dunif(0, 10)
  r ~ dunif(0, 1)
	
  # Derived quantities
  psi.sample <- mean(z[])   # Proportion of occupied sites in the sample
  psi <- 1 - exp( -lambda)  # Prob(N > 0 | lambda), probability that a random site is occupied.
}

For this trial run I used a Uniform(0, 10) prior and will check to see if the upper limit constrains the posterior. It would be better to use a gamma prior with sensible parameters for a run with real data that you wished to report.

This is saved in the file named "modelRN.txt" in the current working directory.

Running the model from R

library(jagsUI)
# The data
str(JAGSdata <- list(y = y, nRep = nRep, nSite = nSite))
List of 3
 $ y    : int [1:200] 3 1 2 2 2 2 0 2 1 0 ...
 $ nRep : num 10
 $ nSite: num 200
# Starting values
Nst <- as.numeric(y > 0)  # N must be >0 if y > 0, so start at 1
inits <- function() list(N = Nst, lambda=runif(1,0,3)) 
# Parameters to monitor
params <- c("lambda", "r", "psi", "psi.sample", "N")
# Takes < 2 mins on my laptop:
out <- jags(JAGSdata, inits, params, "modelRN.txt", DIC=FALSE,
  n.chains = 3, n.adapt=1000, n.burnin=0, n.iter=100000, parallel=TRUE)

Here are the first few rows of the output:

            mean    sd  2.5%   50% 97.5% overlap0 f Rhat  n.eff
lambda     1.183 0.184 0.884 1.163 1.600    FALSE 1    1  60856
r          0.114 0.016 0.083 0.114 0.146    FALSE 1    1  31678
psi        0.689 0.054 0.587 0.687 0.798    FALSE 1    1  36368
psi.sample 0.689 0.041 0.620 0.685 0.780    FALSE 1    1  58580
N[1]       2.110 0.932 1.000 2.000 4.000    FALSE 1    1  80386
N[2]       1.390 0.638 1.000 1.000 3.000    FALSE 1    1 300000
N[3]       1.703 0.807 1.000 2.000 4.000    FALSE 1    1 255018

For this data set, the estimate of \( \lambda \) and \( \psi \) are a bit low, as we saw for the ML estimates. I'd like to compare the estimates of \( N \) with the true values for each of the 200 sites, but plotting  \( \hat{N} \) versus  \( N \) (graph not shown) is not very informative because integers are involved. Instead, we'll tabulate the mean estimate  \( \hat{N} \) for each of the true values of \( N \):

Nhat <- unlist(out$mean)[-(1:4)]
plot(jitter(N), Nhat, ylim=range(0, Nhat))
abline(0, 1, col='red')
abline(v=mean(N), h=mean(Nhat), lty=3)

# Get mean(Nhat) for each value of N:
table(N)
N
 0  1  2  3  4  6  7 
55 70 40 23 10  1  1 
tapply(Nhat, N, mean)
        0         1         2         3         4         6         7 
0.3663258 1.0717541 1.6843473 1.9255339 2.2463930 2.1042233 4.2685067 

There are only 2 sites with >4 individuals, so we shouldn't take too much notice of those. The estimate for \( N = 0 \) is inevitably biased high (as any value >0 in the MCMC chain will immediately produce an estimate >0), those for \( N > 1 \) appear to be biased low, at least for this simulated sample.

Investigating accuracy

Probability of occupancy

I ran the simulation 1000 times (with different seeds) and estimated occupancy using both the RN model and the null model (the psi(.)p(.) model). Comparing the estimates with the data-generating value, \( \psi = 1 - e^{-\lambda} = 0.777 \), the null model showed a bias of -0.077 and the RN model was almost unbiased (0.006).

Root mean squared errors (RMSEs) calculated using the data-generating value were 0.087 for the null model and 0.053 for the RN model.

Some of the error here is due to sampling error, the fact that the proportion of occupied sites in the sample is not the same as the data-generating \( \psi \). Comparing estimates with the proportion in sample, psiSample , seems a fairer way to test the two methods. Defining "error" as the difference between the estimate and psiSample gives RMSEs of 0.083 for the null model and 0.044 for the RN model.

While doing this, I used the loo package ( Vehtari et al 2016 ) to estimate looaic and waic model selection criteria. There were a few issues with this which I'll look at in a separate post, but in all 1000 simulations the RN model was preferred over the null model by both criteria.

I also ran the ML analyses, using the occSS0 and occSSrn0 functions in wiqid. These gave very similar results, as one would expect as we are running the Bayesian analysis with flat priors. AIC for the RN model was lower than the null model for 86% of the samples.

Now the acid test: how does the RN model perform when probability of detection does not depend on number available for capture, merely on whether that is > 0? I did 1000 simulations with data generated with the same biological model and a detection model where \( y_i ~ Binomial(p * z_i, n_i) \) where \( z_i = 1 \) when \( N > 0 \) and 0 otherwise, and \( p = 0.179 \), the expected value of \( 1 - (1 - r)^{N_i} \) for \( N > 0 \). The exact same code was used for the analysis. Now the null model has lower bias (0.003 vs 0.124) and smaller RMSE (0.031 vs 0.130) when compared with psiSample . But both looaic and waic indicated that the RN model was better for all except 1 of the data sets. MLE analysis again gave very similar results, except that AIC was lower for the null model for 88% of the samples.

Based on these simulations, choosing the correct model gives more accurate results and lower bias. AIC picks the right model >80% of the time, but not always. LOO and WAIC were almost always lower for the RN model, even when the data came from a null model.

Numbers at each site

For the same 1000 simulations of the RN model, I checked the estimates of the number of individuals available for detection at each site as I did above. All 1000 had sites with values of \( N \) between 0 and 4, and 931 simulations had sites with \( N =5 \). Higher values were uncommon, so are not included here. The mean estimates across all the simulations are in the table below:

For this particular scenario, it seems that the estimates are more or less unbiased for \( N = 2 \), but are positively biased for \( N < 2 \) and negatively biased for \( N > 2 \). The results would certainly be different with other values for \( \lambda \) and \( r \) and for different numbers of replicates. But the take home message would likely be similar: (a) don't regard the output as unbiased estimates of abundance, or even "number of animals available for detection", at each location; but (b) the estimates may still be useful to investigate trends in abundance with habitat covariates.

The Royle-Nichols model represents an important extension of standard occupancy modelling, linking abundance and detection probability through a simple yet powerful framework. Rather than treating occupancy as a binary state, this approach recognises that detection probability is influenced by the number of individuals present at a site. When more animals are available to be detected, the overall detection probability increases naturally. This insight allows researchers to estimate both occupancy and relative abundance from detection-nondetection data alone, without requiring individual identification. The model has proven especially valuable for studies using camera traps or sign surveys where direct counts are difficult but detection patterns carry information about underlying population size.

Bayesian methods offer a natural framework for fitting occupancy models, particularly when dealing with complex hierarchical structures or sparse data. Using Markov chain Monte Carlo sampling, analysts can obtain full posterior distributions for all model parameters, including derived quantities such as the probability of occupancy at unsampled locations. Software packages like JAGS and BUGS make these analyses accessible to wildlife researchers, allowing them to specify models in a straightforward manner. The Bayesian approach also facilitates the incorporation of prior information from previous studies or expert knowledge, which can be especially helpful when detection probabilities are low or sample sizes are limited.

Covariates play a central role in modern occupancy analyses, enabling researchers to investigate the ecological factors that influence both occupancy and detection. Habitat characteristics such as vegetation cover, distance to water, or elevation might affect where a species chooses to live. Meanwhile, survey-specific covariates like temperature, time of day, or observer experience can influence detection probability. By including these variables in the model, analysts can account for heterogeneity that would otherwise bias estimates. This flexibility makes occupancy models far more useful than simple presence-absence comparisons, as they can reveal the underlying ecological processes driving species distributions across landscapes.

Spatial capture-recapture methods extend the principles of occupancy modelling to estimate population density directly from detection data. When animals can be identified individually through natural markings or artificial tags, researchers can fit spatial models that relate detection probability to distance from an animal's home range centre. These models produce estimates of density per unit area, which is often more meaningful than simple occupancy for conservation and management decisions. The spatial component also allows investigators to map the distribution of activity centres across the landscape, providing insights into habitat use and movement patterns that would be impossible to obtain from nonspatial approaches alone.

The practical implementation of occupancy models requires careful attention to study design and data collection protocols. Repeated surveys at each site are essential to separate nondetection from true absence, and the timing of surveys should reflect the target species' activity patterns. Modern software tools have made these analyses increasingly accessible, but understanding the underlying assumptions remains crucial. Violations of assumptions such as population closure or independent detections can lead to biased estimates if not properly addressed. Researchers should conduct goodness-of-fit tests and consider alternative model structures when appropriate, ensuring that their conclusions about species occupancy and distribution are supported by robust statistical evidence.