# Occupancy model error with K value (number of visits) – SOLVED

**URL:** <https://discourse.mc-stan.org/t/occupancy-model-error-with-k-value-number-of-visits-solved/5411>\
**Category:** Modeling\
**Tags:** ecology\
**Created:** [August 28, 2018, 6:33pm UTC](https://discourse.mc-stan.org/t/occupancy-model-error-with-k-value-number-of-visits-solved/5411 "2018-08-28T18:33:41Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [August 28, 2018, 6:33pm UTC](https://discourse.mc-stan.org/t/occupancy-model-error-with-k-value-number-of-visits-solved/5411/1 "2018-08-28T18:33:41Z")

</div>

Afternoon all,

I am slowing working my way through different ecology models using Stan. I stumbled across this occupancy model in the case studies.  
[http://mc-stan.org/users/documentation/case-studies/dorazio-royle-occupancy.html](http://mc-stan.org/users/documentation/case-studies/dorazio-royle-occupancy.html)

I ran the example as is on:  
R version 3.4.4 (2018-03-15)  
Platform: x86\_64-apple-darwin15.6.0 (64-bit)  
Running under: OS X El Capitan 10.11.6  
rstan\_2.17.3

That runs just fine. When I roll my data into it (attached here as a csv) I get the following error:  
Error in new\_CppObject\_xp(fields$.module, fields$.pointer, …) :  
Exception: model346a164a8a02\_f46f730e64151416119f700ecfb76ff2\_namespace::model346a164a8a02\_f46f730e64151416119f700ecfb76ff2: x[k0\_\_][k1\_\_] is 5, but must be less than or equal to 3 (in ‘model346a164a8a02\_f46f730e64151416119f700ecfb76ff2’ at line 34)  
failed to create the sampler; sampling not done

We visit our sites three times a year. It’s possible for us to have hundreds of a single species in a trap.  
Our n = 28, J = 20, K = 3, S = 1000. I set the super community high since at one site we 889 ants of a single species caught. I’ve also set S as high as 10000 but got the same error.

Based on the error I am guessing the model code has a part where it’s expecting a certain number of visits based on the counts. However I haven’t figured out where that is in the model.

Model from the case study

```stan
model functions {
matrix cov_matrix_2d(vector sigma, real rho) {
matrix[2,2] Sigma;
Sigma[1,1] = square(sigma[1]);
Sigma[2,2] = square(sigma[2]);
Sigma[1,2] = sigma[1] * sigma[2] * rho;
Sigma[2,1] = Sigma[1,2];
return Sigma;
}

real lp_observed(int x, int K, real logit_psi, real logit_theta) {
return log_inv_logit(logit_psi)
+ binomial_logit_lpmf(x | K, logit_theta);
}

real lp_unobserved(int K, real logit_psi, real logit_theta) {
return log_sum_exp(lp_observed(0, K, logit_psi, logit_theta),
log1m_inv_logit(logit_psi));
}

real lp_never_observed(int J, int K, real logit_psi, real logit_theta,
real Omega) {
real lp_unavailable = bernoulli_lpmf(0 | Omega);
real lp_available = bernoulli_lpmf(1 | Omega)
+ J * lp_unobserved(K, logit_psi, logit_theta);
return log_sum_exp(lp_unavailable, lp_available);
}
}
data {
int<lower=1> J; // sites within region
int<lower=1> K; // visits to sites
int<lower=1> n; // observed species
int<lower=0, upper=K> x[n,J]; // visits when species i was detected at site j
int<lower=n> S; // superpopulation size
}
parameters {
real alpha; // site-level occupancy
real beta; // site-level detection
real<lower=0, upper=1> Omega; // availability of species

real<lower=-1,upper=1> rho_uv; // correlation of (occupancy, detection)
vector<lower=0>[2] sigma_uv; // sd of (occupancy, detection)
vector[2] uv[S]; // species-level (occupancy, detection)
}
transformed parameters {
vector[S] logit_psi; // log odds of occurrence
vector[S] logit_theta; // log odds of detection
for (i in 1:S)
logit_psi[i] = uv[i,1] + alpha;
for (i in 1:S)
logit_theta[i] = uv[i,2] + beta;
}
model {
// priors
alpha ~ cauchy(0, 2.5);
beta ~ cauchy(0, 2.5);
sigma_uv ~ cauchy(0, 2.5);
(rho_uv + 1) / 2 ~ beta(2, 2);
uv ~ multi_normal(rep_vector(0, 2), cov_matrix_2d(sigma_uv, rho_uv));
Omega ~ beta(2,2);

// likelihood
for (i in 1:n) {
1 ~ bernoulli(Omega); // observed, so available
for (j in 1:J) {
if (x[i,j] > 0)
target += lp_observed(x[i,j], K, logit_psi[i], logit_theta[i]);
else
target += lp_unobserved(K, logit_psi[i], logit_theta[i]);
}
}
for (i in (n + 1):S)
target += lp_never_observed(J, K, logit_psi[i], logit_theta[i], Omega);
}
generated quantities {
real<lower=0,upper=S> E_N = S * Omega; // model-based expectation species
int<lower=0,upper=S> E_N_2; // posterior simulated species
vector[2] sim_uv;
real logit_psi_sim;
real logit_theta_sim;

E_N_2 = n;
for (i in (n+1):S) {
real lp_unavailable = bernoulli_lpmf(0 | Omega);
real lp_available = bernoulli_lpmf(1 | Omega)
+ J * lp_unobserved(K, logit_psi[i], logit_theta[i]);
real Pr_available = exp(lp_available
- log_sum_exp(lp_unavailable, lp_available));
E_N_2 = E_N_2 + bernoulli_rng(Pr_available);
}

sim_uv = multi_normal_rng(rep_vector(0,2),
cov_matrix_2d(sigma_uv, rho_uv));
logit_psi_sim = alpha + sim_uv[1];
logit_theta_sim = beta + sim_uv[2];
}

```

[ant\_species\_by\_site\_2016.csv](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/3/3f7d57a8602b77140d1c9c25a702bb2f3336adca.csv) (1.7 KB)

---

<div class="post-metadata">

**Author:** ![mbjoseph](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mbjoseph/32/870_2.png) [@mbjoseph](https://discourse.mc-stan.org/u/mbjoseph)\
**Post date:** [August 28, 2018, 10:34pm UTC](https://discourse.mc-stan.org/t/occupancy-model-error-with-k-value-number-of-visits-solved/5411/2 "2018-08-28T22:34:26Z")

</div>

Hi @Ara_Winter - great to see more occupancy models on this forum! The values in the observation matrix x should be between 0 and K, where K is the number of visits to each site. That’s why the error is raised - it’s telling you that a particular element in the matrix x is 5, but must be less than or equal to K=3.

Looking at the attached data, the values in x are pretty large, in the hundreds, which along with this:

> [@Ara\_Winter](#):
>
> at one site we 889 ants of a single species caught

makes me suspect that maybe the values in x are counts of individual ants, rather than integers that represent the number of surveys that each species was observed at each site, which would be between 0 and K. In the case study when x is described, it’s introduced as:

- x: n \times J matrix where x\_{i, j} \in 1:K is the number of visits in which species i was observed at site j

Now that I read that, I see there’s a typo (it should be x\_{i, j} \in 0:K), but the error you’re getting is because you have values x\_{i, j} that are greater than K. In other words, the way the data are encoded now imply that species 14 at site 2 was detected on 889 out of 3 surveys, which is impossible. This species could have been detected at this site at most 3 times out of 3 surveys.

To fix this, you’ll probably want to go back to the original data and construct the matrix x based on the number of surveys that each species was seen at each site, so that the values are all between 0 and K.

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [August 28, 2018, 10:38pm UTC](https://discourse.mc-stan.org/t/occupancy-model-error-with-k-value-number-of-visits-solved/5411/3 "2018-08-28T22:38:32Z")

</div>

@mbjoseph Funny I was reading the original paper and it dawned on me that I had the matrix setup wrong :(

Thanks for clarifying what I was doing wrong. That solves the problem.

Hopefully once I get my head wrapped around this model a bit better I’ll have some fancified version to contribute.

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [August 30, 2018, 1:41pm UTC](https://discourse.mc-stan.org/t/occupancy-model-error-with-k-value-number-of-visits-solved/5411/4 "2018-08-30T13:41:37Z")

</div>

Just for completeness the model is up and running.

 ![ant%20detection%20of%20species%20at%20sites%20over%20visits](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/b/b6238a04c17e61fdd5b3fbdc29b156269e13180c.png)

 ![ants_traceplot](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/b/b2793a3f0b0e1dee3317a7b0e303cfbb7f35ce74.png)

 ![Site-Level%20Parameter%20Posterior](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/e/e7d0ebc686e61fb8bd8ecb2dbc384e079bf35bba.png)

With the attached data file here – in case anyone at home wants to play with ant data.  
[ant\_species\_by\_site\_2016.csv](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/d/de198740046667d372e1b2d4266bfc9fe0256ca2.csv) (1.5 KB)