# Royle and Nichols

Dear people,

I am trying to implement the Royle and Nicholson model for abundance-induced heterogeneity in detection probability.

There are i = 1,2,…,M sites sampled J times and the apparent presence/absence of a species is recorded during each visit.

Assume that the population size exposed to sampling varies among sites. Let N[i] be the population size at site i.

If individuals are detected independently of one another with probability r, then site-specific detection probability is

p[i] = 1 - (1 - r)^N[i]

To implement this in Stan, I follow

and write:

``````data {
int<lower=1> M;                 // total number of sites
int<lower=1> J;                 // Number of replicates at each site
int<lower=0,upper=1> Y[M, J];   // response variable
int<lower=0,upper=1> yz[M];     // at least one detection
int<lower=0> n_max;             // Upper bound of population size
}

parameters {
real<lower=0> lambda;       // mean population size
real<lower=0,upper=1> r;    // per capita detection probability
}

model {
r ~ beta(2,2);
lambda ~ cauchy(0, 10);

for (n in 1:M){
vector[n_max - yz[n] + 1] lp;
for (j in 1:(n_max - yz[n] + 1)){
lp[j] = poisson_lpmf(yz[n] + j - 1 | lambda)
+ binomial_lpmf(Y[n] | 1, 1 - (1 - r)^(yz[n] + j - 1) );
}
target += log_sum_exp(lp);
}
}
``````

And try to run a simple example:

``````
M = 10
J = 5

stan_dat <- list(
M = M,
Y = matrix(0, M, J),
yz = numeric(M),
J = J,
n_max = 10
)

fit <- stan(file = 'RoyleNicholson.stan',
data = stan_dat,
iter = 100,
chains = 1)
``````

But I get an error message saying that the iGradient evaluated at the initial value is not finite. I played around a bit and found that the problem seem to be with

binomial_lpmf(Y[n] | 1, 1 - (1 - r)^(yz[n] + j - 1) );

I have tried for example

binomial_lpmf(Y[n] | 1, 1 - (1 - 0.8)^0 );

and it works. But I get an error if I do

binomial_lpmf(Y[n] | 1, 1 - (1 - r)^0 );

Any help would be very much appreciated!

2 Likes

You have to handle `yz=0;j=1` case separately.

``````    vector[n_max - yz[n] + 1] lp;
if (yz[n] == 0) {
lp = poisson_lpmf(0 | lambda)
+ binomial_lpmf(Y[n] | 1, 0. );
} else {
lp = poisson_lpmf(yz[n] | lambda)
+ binomial_lpmf(Y[n] | 1, 1 - (1 - r)^(yz[n]) );
}
for (j in 2:(n_max - yz[n] + 1)) {
lp[j] = poisson_lpmf(yz[n] + j - 1 | lambda)
+ binomial_lpmf(Y[n] | 1, 1 - (1 - r)^(yz[n] + j - 1) );
}
target += log_sum_exp(lp);
``````
1 Like

Thanks Nico!

First of all, it’s NICHOLS not Nicholson. Sorry…

I tried your suggestion and it now works. There’s probably a more efficient way to do it, but here’s what I did:

``````data {
int<lower=1> M;                 // total number of sites
int<lower=1> J;                 // Number of replicates at each site
int<lower=0,upper=1> Y[M, J];   // response variable
int<lower=0,upper=1> yz[M];     // at least one detection
int<lower=0> n_max;             // Upper bound of population size
}

parameters {
real<lower=0> lambda;       // mean population size
real<lower=0,upper=1> r;    // per capita detection probability
}

model {
r ~ beta(2,2);
lambda ~ cauchy(0, 10);

for (n in 1:M){
vector[n_max - yz[n] + 1] lp;
if(yz[n] == 0){
lp = poisson_lpmf(0 | lambda) + 1;
}
else lp = poisson_lpmf(1 | lambda) + binomial_lpmf(Y[n] | 1, r);

for (j in 2:(n_max - yz[n] + 1)){
lp[j] = poisson_lpmf(yz[n] + j - 1 | lambda)
+ binomial_lpmf(Y[n] | 1, 1 - (1 - r)^(yz[n] + j - 1) );
}
target += log_sum_exp(lp);
}
}
``````

Not the crux of the problem, but if yz[n] = 0 (no detections), then shouldn’t the binomial contribution to lp be 0 instead 1, since it’s a probability of 1 but in the log scale?

1 Like