Royle and Nichols

Modeling
1

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
2

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

    vector[n_max - yz[n] + 1] lp;
    if (yz[n] == 0) {
      lp[1] = poisson_lpmf(0 | lambda) 
      + binomial_lpmf(Y[n] | 1, 0. );
    } else {
      lp[1] = 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
3

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[1] = poisson_lpmf(0 | lambda) + 1;
    }
    else lp[1] = 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);
  }
}
4

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

1 Like