occupancy model using RStan giving funky results

[Quresh Latif]

I am trying to run occupancy models using RStan. I simulated some data for 240 sites and using a randomly generated covariate of occupancy. This is a hidden data model where the occupancy status (Z) is latent but partially observed. If the species is detected, Z = 1. If the species is not detected, Z could equal either with detection probability (p) governing whether the species is detected when Z = 1. The complete code for simulating a dataset, defining what I think should be the right model for RStan, and executing the model are pasted below. The model compiles and runs, but I get really funky results. Parameters have either really high or really low n effectives and R_hats are all extremely high. The lp__ parameter (I think this is the same as deviance right?) is not calculated successfully at all. I think I must be doing something wrong presumably in the order in which I declare/define variables. Appreciate any insight.

 

library(rstan)
set_cppo("fast")  # for best running speed
#set_cppo("debug")  # for debugging
p <- 0.7   #Values for p at intercept (mean of covariate values)
psi <- 0.5 #Values for psi at intercept (mean of covariate value(s))
n.site <- 240 #no. sites
n.visit <- 2 #no. visits (will try n.visit<-c(5,10,20,40) for n.site = 30 and p = 0.1 later)
B.x1 <- 1   #Slope of linear relationship btwn psi and covariate 1
B.x2 <- 0  #Slope of linear relationship btwn psi and covariate 2
A.x1 <- 0   #Slope of linear relationship btwn p and covariate 1
cov.x <- 0    #Correlation between x1 and x2

 

# Function for simulating data
expit <-
  function(x){
    exp(x)/(1+exp(x))
  }

logit <-
  function(x){
    log(x/(1-x))
  }

Dat.sim <- function(p0,psi0,n.sts,n.vsts,b.x1,b.x2,a.x1,covx) {
  Y<-matrix(NA,nrow=n.sts,ncol=n.vsts)
  mns<-c(0,0)
  Sgma<-as.matrix(rbind(c(1,covx),c(covx,1)))
  x<-mvrnorm(n.sts,mns,Sgma)
  psi<-expit(logit(psi0)+b.x1*x[,1]+b.x2*x[,2])
  p<-expit(logit(p0)+a.x1*x[,1])
  for(ii in 1:nrow(Y)) {
    z<-rbinom(1,1,psi[ii])
    Y[ii,]<-rbinom(n.vsts,1,z*p[ii])
  }
  Dat<-list(Y=Y,psi=psi,p=p,x1=x[,1],x2=x[,2])
  return(Dat)
}

 

# Simulate data
dat <- Dat.sim(p,psi,n.site,n.visit,B.x1,B.x2,A.x1,cov.x)
R <- nrow(dat$Y)
J <- ncol(dat$Y)
Data <- list(Y=dat$Y,R=R,J=J,x1=dat$x1,x2=dat$x2)

 

# Set initial values for HMC sampling
z.init <- apply(dat$Y,1,max)

inits <- function()
  list(Z=z.init,psi_mean=runif(1),p_mean=runif(1),B_x1=rnorm(1),B_x2=rnorm(1),A_x1=rnorm(1))

 

#Define which parameters to save.
parameters <- c("psi","p","B0","B_x1","B_x2","A_x1")

 

# Model code
mod_occupancy_simple <- '
data {
int<lower=0> R; // number of sites
int<lower=0> J; // number of visits
int<lower=0,upper=J> Y[R,J]; // detection-nondetection data
real x1[R]; // covariate x1
real x2[R];
}

parameters {
real<lower=0,upper=1> psi_mean;
real<lower=0,upper=1> p_mean;
real B_x1;   
real B_x2;   
real A_x1;   
}

transformed parameters {
real B0;
real A0;
real psi[R];
real p[R,J];
B0 <- logit(psi_mean);
A0 <- logit(p_mean);
for (r in 1:R) {
  psi[r] <- inv_logit(B0 + B_x1*x1[r] + B_x2*x2[r]);
  for (j in 1:J) {
    p[r,j] <- inv_logit(A0 + A_x1*x1[r]);
    }
  }
}

model {
int Z[R];   // declare latent variable indicating occupancy status of site.
real mu_p[R,J];

psi_mean ~ uniform(0,1);
p_mean ~ uniform(0,1);
B_x1 ~ normal(0,10);   
B_x2 ~ normal(0,10);
A_x1 ~ normal(0,10);

//loop over all sites R
for (r in 1:R) {
  Z[r] ~ bernoulli(psi[r]);     // define latent variable indicating occupancy status of site.
  for (j in 1:J) {
    mu_p[r,j] <- Z[r]*p[r,j];
    Y[r,j] ~ bernoulli(mu_p[r,j]);
    }
  }
}
'

# Settings for HMC chains
nc <- 6
nb <- 1000
ni <- 2000
nt <- 10

#Run model
out <- stan(model_code = mod_occupancy_simple, data = Data, pars = parameters, init = inits,
            iter = ni, chains = nc, warmup = nb, thin = nt)

 

[Quresh Latif]

I tried simplifying further by removing the covariates. Now there are only two parameters to estimate, occupancy and detection probability. Estimates coming out of RStan seem pretty off, n_effective = the number of samples drawn, and R_hat is huge. At the very least, there seems to be something wrong with how the convergence diagnostics are being calculated, so I can't tell if I'm sampling the posterior distributions adequately. Longer chains don't seem to be bringing down the R hat or making the estimates look any better. Code for simplified model without covariates is below. Any ideas on what I am doing wrong would be appreciated.

 

library(rstan)
set_cppo("fast")  # for best running speed
#set_cppo("debug")  # for debugging
p <- 0.7   #Values for p at intercept (mean of covariate values)
psi <- 0.5 #Values for psi at intercept (mean of covariate value(s))

n.site <- 1000 #no. sites

n.visit <- 2 #no. visits (will try n.visit<-c(5,10,20,40) for n.site = 30 and p = 0.1 later)

 

# Functions for simulating data
expit <-
  function(x){
    exp(x)/(1+exp(x))
  }

logit <-
  function(x){
    log(x/(1-x))
  }

Dat.sim <- function(p,psi,n.sts,n.vsts) {
  Y <- matrix(NA,nrow=n.sts,ncol=n.vsts)
  for(ii in 1:nrow(Y)) {
    z<-rbinom(1,1,psi)
    Y[ii,]<-rbinom(n.vsts,1,z*p)
  }
  Dat<-list(Y=Y,psi=psi,p=p)
  return(Dat)
}

 

# Simulate data
dat <- Dat.sim(p,psi,n.site,n.visit)

R <- nrow(dat$Y)
J <- ncol(dat$Y)

Data <- list(Y=dat$Y,R=R,J=J)

# Set initial values for HMC sampling
z.init <- apply(dat$Y,1,max)

inits <- function()
  list(Z=z.init,psi=runif(1),p=runif(1))

#Define which parameters to save.

parameters <- c("psi","p")

# Model code
mod_occupancy_simple <- '
data {
int<lower=0> R; // number of sites
int<lower=0> J; // number of visits
int<lower=0,upper=J> Y[R,J]; // detection-nondetection data
}

parameters {
real<lower=0,upper=1> psi;
real<lower=0,upper=1> p;
}

model {
int Z[R];   // declare latent variable indicating occupancy status of site.
real mu_p[R,J];

psi ~ uniform(0,1);
p ~ uniform(0,1);

//loop over all sites R
for (r in 1:R) {

  Z[r] ~ bernoulli(psi);     // define latent variable indicating occupancy status of site.

  for (j in 1:J) {

    mu_p[r,j] <- Z[r]*p;

    Y[r,j] ~ bernoulli(mu_p[r,j]);
    }
  }
}
'

# Settings for HMC chains
nc <- 6
nb <- 1000

ni <- 5000
nt <- 1

#Run model
out <- stan(model_code = mod_occupancy_simple, data = Data, pars = parameters, init = inits,
            iter = ni, chains = nc, warmup = nb, thin = nt)

 

out
Inference for Stan model: mod_occupancy_simple.
6 chains, each with iter=5000; warmup=1000; thin=1; 
post-warmup draws per chain=4000, total post-warmup draws=24000.

     mean se_mean  sd 2.5%  25%  50%  75% 97.5% n_eff         Rhat
psi   0.3       0 0.2    0  0.2  0.2  0.4   0.8 24000          Inf
p     0.7       0 0.3    0  0.5  0.8  0.9   1.0 24000 8.020037e+15
lp__ -Inf     NaN NaN -Inf -Inf -Inf -Inf  -Inf 24000          NaN

Samples were drawn using NUTS(diag_e) at Thu Jan 02 11:58:49 2014.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

 

out2 <- stan(fit=out,data=Data,init=inits,iter=5000,chains=nc,warmup=nb,thin=nt)

 

out2
Inference for Stan model: mod_occupancy_simple.
6 chains, each with iter=50000; warmup=1000; thin=1; 
post-warmup draws per chain=49000, total post-warmup draws=294000.

     mean se_mean  sd 2.5%  25%  50%  75% 97.5%  n_eff         Rhat
psi   0.5       0 0.3  0.1  0.3  0.5  0.8   0.8 294000          Inf
p     0.5       0 0.3  0.0  0.2  0.6  0.8   0.9 294000 9.833529e+15
lp__ -Inf     NaN NaN -Inf -Inf -Inf -Inf  -Inf 294000          NaN

Samples were drawn using NUTS(diag_e) at Thu Jan 02 11:50:26 2014.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

[Jiqiang Guo]

stan does not support your model as the latent parameters are discrete (which should be declared parameters though). 

After all, n_eff and rhat should be combined as evidence for convergence diagnostics.

Jiqiang 

On Thu, Jan 2, 2014 at 2:00 PM, Quresh Latif <latif...@gmail.com> wrote:

I tried simplifying further by removing the covariates. Now there are only two parameters to estimate, occupancy and detection probability. Estimates coming out of RStan seem pretty off, n_effective = the number of samples drawn, and R_hat is huge. At the very least, there seems to be something wrong with how the convergence diagnostics are being calculated, so I can't tell if I'm sampling the posterior distributions adequately. Longer chains don't seem to be bringing down the R hat or making the estimates look any better. Code for simplified model without covariates is below. Any ideas on what I am doing wrong would be appreciated.

library(rstan)
set_cppo("fast")  # for best running speed

#set_cppo("debug")  # for best running speed

# Set values for simulating data

p <- 0.7   #Values for p at intercept (mean of covariate values)
psi <- 0.5 #Values for psi at intercept (mean of covariate value(s))

n.site <- 1000 #no. sites

n.visit <- 2 #no. visits (will try n.visit<-c(5,10,20,40) for n.site = 30 and p = 0.1 later)

# Function for simulating data
expit <-
  function(x){
    exp(x)/(1+exp(x))
  }

logit <-
  function(x){
    log(x/(1-x))
  }

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.

[Bob Carpenter]

The key is that sampling statements, like x ~ foo(theta), don't
actually set the value of x. Instead, they increment
the total model log probability with log foo(x|theta).
This means x must already be set.

We say it right in the introduction of the manual, but I fear that people
don't read any of the manual, but rather read the
R getting started guide and then they work by analogy to BUGS,
which then leads them astray as in this example.


We initialize local integer values to 0, so that when Quresh did this:

model {
int Z[R]; // declare latent variable indicating occupancy status of site.
real mu_p[R,J];
psi ~ uniform(0,1);
p ~ uniform(0,1);
//loop over all sites R
for (r in 1:R) {

Z[r] ~ bernoulli(psi); *** PROBLEM HERE, Z[r] not set ***

for (j in 1:J) {

mu_p[r,j] <- Z[r]*p; *** PROBLEM: mu_p[r,j] can be 0 ***

Y[r,j] ~ bernoulli(mu_p[r,j]);
}
}

rather than getting an error message, it just behaves unexpectedly, because
each of the Z[r] is set to 0.


Usually, you can solve this problem the way Jiqiang suggests,
by moving the local variable into a parameter.

But for integers, you need to marginalize them out of the log
probability function. We discuss how to do this in the manual.


Quresh's model also brings up another issue with inits for local
variables. For real-valued local variables, we initialize to NaN.
There's no equivalent for integer values, and I think we're
initializing to 0 (I'd have to check the generated code).

One thing we could do is try to catch uses of variables
before they're ever set. I'm adding an issue for that:

https://github.com/stan-dev/stan/issues/477

I assigned myself, but don't know when I'll be able to
get around to it.


Also, it's not necessary to set variables BUGS-style by defining
locals for whole loops. The for loop over j could be replaced
with a single vectorized statement with no locals:

Y[r] ~ bernoulli(Z * p);

at least if Z were defined as a vector rather than as an array.

- Bob

[Quresh Latif]

Thanks for your responses. To be honest, I do not have a strong programming background or deep knowledge of Bayesian statistics, so a lot of what's in the manual went over my head. For example, I'm not yet clear on what is meant by "...increment the total model log probability with log foo(x|theta)." and why this means that "x must already be set." But, I will look back at the manual and continue working on my understanding.

 

One thing I want to be clear on: if I do marginilize 'Z' out of the log probability function, will this model work or do I need an alternative that does not involve explicitly defining 'Z'? You mentioned that local integer parameter values are initialized at 0, which would definitely not work with this model, because when a species is observed at a site, Z must equal 1 at that site. So, the question is, will marginilization of Z as described in the manual take care of this problem, or is this model simply a no-go in the current version of STAN?

 

Thanks again for all your help.

[Bob Carpenter]

On 1/3/14, 5:44 PM, Quresh Latif wrote:
> Thanks for your responses. To be honest, I do not have a strong programming background or deep knowledge of Bayesian
> statistics, so a lot of what's in the manual went over my head. For example, I'm not yet clear on what is meant by
> "...increment the total model log probability with log foo(x|theta)."

Stan just defines the log probability function for p(data,parameters),
which is proporitional to the posterior p(params|data) given Bayes's rule.

> and why this means that "x must already be set."

Stan is more like a programming language than like BUGS. So it's
just like in R --- if you use a variable in R, it must already be set.

You can't write an R program that looks like this:

Z <- rep(0,R);

for (r in 1:R)

Z[r] <- rbinom(1,1,psi[r]);

because psi hasn't been set to any value at the point you
execute the rbinom(1,1,psi[r]) expression.

> But, I will look back at the manual and continue working on my understanding.
> One thing I want to be clear on: if I do marginilize 'Z' out of the log probability function, will this model work

Yes, it should.

> I need an alternative that does not involve explicitly defining 'Z'?

Hmm. I fyou marginalize out Z, you won't have an explicit Z.
For instance, if b is discrete and has two possible values, 0 and 1,
then to marginalize b out of p(a,b), you have

p(a) = SUM_(b in {0,1}) p(a,b)

= p(a,0) + p(a,1)

So b doesn't appear in the final expression for p(a).

> You mentioned that local integer parameter values
> are initialized at 0, which would definitely not work with this model, because when a species is observed at a site, Z
> must equal 1 at that site. So, the question is, will marginilization of Z as described in the manual take care of this
> problem, or is this model simply a no-go in the current version of STAN?

Yes, it should work. But it'll require some algebra and
won't look at all like the model you have now.

> Thanks again for all your help.

Here's the critical bit where you need to marginalize Z[r] out:

Z[r] ~ bernoulli(psi[r]); // define latent variable indicating occupancy status of site.

for (j in 1:J) {

mu_p[r,j] <- Z[r]*p[r,j];

Y[r,j] ~ bernoulli(mu_p[r,j]);
}

If Z[r] = 1, which happens with probability psi[r], you get:

for (j in 1:J) {

mu_p[r,j] <- 1 * p[r,j];

Y[r,j] ~ bernoulli(mu_p[r,j]);
}

or equivalently

for (j in 1:J)

Y[r,j] ~ bernoulli(p[r,j]);

which can be vectorized in Stan along with the probability psi[r]
to produce log probability

log(psi[r]) + bernoulli_log(Y[r], p[r]);

Now if Z[r] = 0, which happens with probability 1 - psi[r], you get:

for (j in 1:J)

Y[r,j] ~ bernoulli(0);

This can be vectorized to log probability

log1m(psi[r]) + bernoulli_log(Y[r],0);

where log1m(u) =def= log(1 - u), but it's faster and more aritmetically stable.

So writing out the marginalized probability in full, you have:

increment_log_prob(log_sum_exp(log(psi[r]) + bernoulli_log(Y[r], p[r]),
log1m(psi[r]) + bernoulli_log(Y[r], 0)));

where log_sum_exp(u,v) = log(exp(u) + exp(v)), but is also faster and more stable.

Next, we recognize that

bernoulli_log(Y[r],0) = 0 if all Y[r,j] = 0
= -infinity otherwise,

Because exp(-infinity) = 0, you can just reduce this whole thing to
a much more efficient form as:

if (sum(Y[r]) == J)
increment_log_prob(log_sum_exp(log(psi[r]) + bernoulli_log(Y[r], p[r]),
log1m(psi[r]));
else
increment_log_prob(log(psi[r]) + bernoulli_log(Y[r], p[r]);

It's ugly, but it should be much much faster than actually sampling
a discrete Z[r].

You could make it even faster by pulling out an indicator for sum(Y[r]) == J
in the transformed data block so you don't have to keep redundantly calculating
the sum.

You should check my work on the algebra, too!

- Bob

[Bob Carpenter]

Speaking of occupancy models, I just did a quick search and found
this simple (no predictors) example in BUGS:

http://www.fisheriesstockassessment.com/TikiWiki/tiki-index.php?page=BUGS+Occupancy

I'll try to add this to the manual. Maybe I'll add a whole chapter on
ecology models. They require some care in coding in Stan, as you saw.

I'm about to add the Cormack-Jolly-Seber model for population size estimation;
it requires some algebra, but works out very neatly and runs very efficiently.

- Bob

On 1/3/14, 5:44 PM, Quresh Latif wrote:

[Quresh Latif]

Great thanks for all your help and suggestions! I am understanding better now, and I'm sure will understand even better once I've gone through and revised my model according to your suggestions. Once I do have a revised version that works correctly, I will post it here, in case it would be of use to anyone.

 

Thanks,

Quresh

[Quresh Latif]

I followed Bob's example and used his formulation (see post entitled "Simple occupancy model...") to analyze simulated data with true parameter values (psi, p, no. sites, no. visits) set to levels in this post. It worked great! The complete code for simulating the data, setting up the model in R, and the results are below. I will work on incorporating covariates and post the result here as well in case anyone is interested.

 

Incidentally, there is some material in the Link and Barker book (Bayesian Inference with Ecological Applications) on complete data likelihoods versus observed data likelihoods that seems relevant to the difference between how occupancy models are typically compiled in BUGS (i.e., the Royle and Kery formulation) versus how they need to be compiled in STAN. It appears to me that STAN requires us to use the observed data likelihood, which is actually the original formulation provided by MacKenzie (2002) for use with maximum likelihood estimation. By going to the observed data likelihood (i.e., not explicitly sampling the posterior for 'Z'), one loses some of the major advantages for shifting from the maximum likelihood to Bayesian framework when running occupancy models. For example, I have been working with multi-species occupancy models (see Dorazio and Royle 2005) from which I am interested in estimating species richness for individual sites. Doing so is easy having estimated 'Z'. I can just add up Z across species to get an estimate of species richness. I imagine there is a way to get at species richness from estimates of occupancy probabilities, but doing so would take some additional thought and would probably be more complicated.

 

Regardless, appreciate all the assistance and the availability of this new tool.

 

Working code and results in R (no covariates (yet)):

library(rstan)
set_cppo("fast")  # for best running speed

# Set values for simulating data
p <- 0.7   #Values for p at intercept (mean of covariate values)
psi <- 0.5 #Values for psi at intercept (mean of covariate value(s))

n.site <- 240 #no. sites

n.visit <- 2 #no. visits (will try n.visit<-c(5,10,20,40) for n.site = 30 and p = 0.1 later)

 

# Functions for simulating data
expit <-
  function(x){
    exp(x)/(1+exp(x))
  }

logit <-
  function(x){
    log(x/(1-x))
  }

Dat.sim <- function(p,psi,n.sts,n.vsts) {
  Y <- matrix(NA,nrow=n.sts,ncol=n.vsts)
  for(ii in 1:nrow(Y)) {
    z<-rbinom(1,1,psi)
    Y[ii,]<-rbinom(n.vsts,1,z*p)
  }
  Dat<-list(Y=Y,psi=psi,p=p)
  return(Dat)
}

 

# Simulate data
dat <- Dat.sim(p,psi,n.site,n.visit)
R <- nrow(dat$Y)
J <- ncol(dat$Y)
Data <- list(Y=dat$Y,R=R,J=J)

 

# Set initial values for HMC sampling

inits <- function()
  list(psi=runif(1),p=runif(1))

 

#Define which parameters to save.
parameters <- c("psi","p")

 

# Model code
mod_occupancy_simple <- '
data {
int<lower=0> R; // number of sites
int<lower=0> J; // number of visits
int<lower=0,upper=J> Y[R,J]; // detection-nondetection data
}

parameters {
real<lower=0,upper=1> psi;
real<lower=0,upper=1> p;
}

model {
// local variables to avoid recomputing log(psi) and log(1-psi)
real log_psi;
real log1m_psi;
log_psi <- log(psi);
log1m_psi <- log1m(psi);

// priors

psi ~ uniform(0,1);
p ~ uniform(0,1);

// likelihood
   for (r in 1:R) {
     if (sum(Y[r]) > 0)
       increment_log_prob(log_psi + bernoulli_log(Y[r],p));
     else
       increment_log_prob(log_sum_exp(log_psi + bernoulli_log(Y[r],p),log1m_psi));
   }
}
'

# Settings for HMC chains
nc <- 6
nb <- 1000
ni <- 5000
nt <- 1

 

#Run model
out <- stan(model_code = mod_occupancy_simple, data = Data, pars = parameters, init = inits,
            iter = ni, chains = nc, warmup = nb, thin = nt)

Results:

 

TRANSLATING MODEL 'mod_occupancy_simple' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'mod_occupancy_simple' NOW.
cygwin warning:
  MS-DOS style path detected: C:/PROGRA~1/R/R-30~1.1/etc/x64/Makeconf
  Preferred POSIX equivalent is: /cygdrive/c/PROGRA~1/R/R-30~1.1/etc/x64/Makeconf
  CYGWIN environment variable option "nodosfilewarning" turns off this warning.
  Consult the user's guide for more details about POSIX paths:
    http://cygwin.com/cygwin-ug-net/using.html#using-pathnames
C:/Program Files/R/R-3.0.1/library/rstan/include//stansrc/stan/agrad/rev/var_stack.hpp:49:17: warning: 'void stan::agrad::free_memory()' defined but not used [-Wunused-function]
SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 1).
Iteration: 5000 / 5000 [100%]  (Sampling)
Elapsed Time: 0.533 seconds (Warm-up)
              1.963 seconds (Sampling)
              2.496 seconds (Total)

SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 2).
Iteration: 5000 / 5000 [100%]  (Sampling)
Elapsed Time: 0.503 seconds (Warm-up)
              1.755 seconds (Sampling)
              2.258 seconds (Total)

SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 3).
Iteration: 5000 / 5000 [100%]  (Sampling)
Elapsed Time: 0.518 seconds (Warm-up)
              1.747 seconds (Sampling)
              2.265 seconds (Total)

SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 4).
Iteration: 5000 / 5000 [100%]  (Sampling)
Elapsed Time: 0.503 seconds (Warm-up)
              1.734 seconds (Sampling)
              2.237 seconds (Total)

SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 5).
Iteration: 5000 / 5000 [100%]  (Sampling)
Elapsed Time: 0.52 seconds (Warm-up)
              1.873 seconds (Sampling)
              2.393 seconds (Total)

SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 6).
Iteration: 5000 / 5000 [100%]  (Sampling)
Elapsed Time: 0.505 seconds (Warm-up)
              1.763 seconds (Sampling)
              2.268 seconds (Total)

> out
Inference for Stan model: mod_occupancy_simple.
6 chains, each with iter=5000; warmup=1000; thin=1; 
post-warmup draws per chain=4000, total post-warmup draws=24000.

       mean se_mean sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
psi     0.5       0  0    0.5    0.5    0.5    0.6    0.6 16893    1
p       0.7       0  0    0.6    0.7    0.7    0.7    0.8 16380    1
lp__ -287.5       0  1 -290.1 -287.9 -287.2 -286.8 -286.6 10545    1

Samples were drawn using NUTS(diag_e) at Sun Jan 05 10:29:05 2014.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

[Quresh Latif]

I was able to get a version of the above model with covariates of occupancy and detection to work. I simulated data with a single covariate of occupancy (x1) and no covariates of detection. I fitted a model with two occupancy covariates (x1, x2) and one covariate for detection, so I expected estimates for B_x1, B_x2, and A_x1 of 1,0, and 0, respectively. I got roughly that.

 

library(rstan)
library(MASS)

set_cppo("fast")  # for best running speed

# Set values for simulating data
p <- 0.7   #Values for p at intercept (mean of covariate values)
psi <- 0.5 #Values for psi at intercept (mean of covariate value(s))
n.site <- 240 #no. sites
n.visit <- 2 #no. visits (will try n.visit<-c(5,10,20,40) for n.site = 30 and p = 0.1 later)

B.x1 <- 1   #Slope of linear relationship btwn psi and covariate 1
B.x2 <- 0  #Slope of linear relationship btwn psi and covariate 2
A.x1 <- 0   #Slope of linear relationship btwn p and covariate 1
cov.x <- 0    #Correlation between x1 and x2

# Functions for simulating data
expit <-
  function(x){
    exp(x)/(1+exp(x))
  }

logit <-
  function(x){
    log(x/(1-x))
  }

Dat.sim <- function(p0,psi0,n.sts,n.vsts,b.x1,b.x2,a.x1,covx) {
  Y<-matrix(NA,nrow=n.sts,ncol=n.vsts)
  mns<-c(0,0)
  Sgma<-as.matrix(rbind(c(1,covx),c(covx,1)))
  x<-mvrnorm(n.sts,mns,Sgma)
  psi<-expit(logit(psi0)+b.x1*x[,1]+b.x2*x[,2])
  p<-expit(logit(p0)+a.x1*x[,1])

  for(ii in 1:nrow(Y)) {

    z<-rbinom(1,1,psi[ii])
    Y[ii,]<-rbinom(n.vsts,1,z*p[ii])
  }
  Dat<-list(Y=Y,psi=psi,p=p,x1=x[,1],x2=x[,2])

  return(Dat)
}

 

# Simulate data
dat <- Dat.sim(p,psi,n.site,n.visit,B.x1,B.x2,A.x1,cov.x)

R <- nrow(dat$Y)
J <- ncol(dat$Y)

Data <- list(Y=dat$Y,R=R,J=J,x1=dat$x1,x2=dat$x2)

 

# Set initial values for HMC sampling
inits <- function()

  list(psi_mean=runif(1),p_mean=runif(1),B_x1=rnorm(1),B_x2=rnorm(1),A_x1=rnorm(1))

 

#Define which parameters to save.

parameters <- c("psi_mean","p_mean","B0","B_x1","B_x2","A0","A_x1")

 

# Model code
mod_occupancy_simple <- '
data {
int<lower=0> R; // number of sites
int<lower=0> J; // number of visits

int<lower=0,upper=1> Y[R,J]; // detection-nondetection data
vector[R] x1; // 1st covariate
vector[R] x2; // 2nd covariate
}

parameters {
real<lower=0,upper=1> psi_mean;
real<lower=0,upper=1> p_mean;
real B_x1;   
real B_x2;   
real A_x1;   
}

transformed parameters {
real B0; // intercept for linear logistic model for occupancy
real A0; // intercept for linear logistic model for detection
vector[R] psi;
matrix[R,J] p;

B0 <- logit(psi_mean);
A0 <- logit(p_mean);

for (r in 1:R) {

psi[r] <- inv_logit(B0 + B_x1*x1[r] + B_x2*x2[r]);

for (j in 1:J) p[r,j] <- inv_logit(A0 + A_x1*x1[r]);
}
}

model {
// local variables to avoid recomputing log(psi) and log(1-psi)

vector[R] log_psi;
vector[R] log1m_psi;

for (r in 1:R) {

  log_psi[r] <- log(psi[r]);
  log1m_psi[r] <- log1m(psi[r]);
}

// priors

psi_mean ~ uniform(0,1);
p_mean ~ uniform(0,1);
B_x1 ~ normal(0,10);
B_x2 ~ normal(0,10);
A_x1 ~ normal(0,10);

// likelihood
   for (r in 1:R) {
     if (sum(Y[r]) > 0)

       increment_log_prob(log_psi[r] + bernoulli_log(Y[r],p[r]));
     else
       increment_log_prob(log_sum_exp(log_psi[r] + bernoulli_log(Y[r],p[r]),log1m_psi[r]));
   }
}
'

# Settings for HMC chains
nc <- 6
nb <- 1000

ni <- 2000
nt <- 10

#Run model
out <- stan(model_code = mod_occupancy_simple, data = Data, pars = parameters, init = inits,
            iter = ni, chains = nc, warmup = nb, thin = nt)

Results:

 

TRANSLATING MODEL 'mod_occupancy_simple' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'mod_occupancy_simple' NOW.
cygwin warning:
  MS-DOS style path detected: C:/PROGRA~1/R/R-30~1.1/etc/x64/Makeconf
  Preferred POSIX equivalent is: /cygdrive/c/PROGRA~1/R/R-30~1.1/etc/x64/Makeconf
  CYGWIN environment variable option "nodosfilewarning" turns off this warning.
  Consult the user's guide for more details about POSIX paths:
    http://cygwin.com/cygwin-ug-net/using.html#using-pathnames

C:/Program Files/R/R-3.0.1/library/rstan/include//stansrc/stan/agrad/rev/var_stack.hpp:49:17: warning: 'void stan::agrad::free_memory()' defined but not used [-Wunused-function]
SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 1).
Iteration: 2000 / 2000 [100%]  (Sampling)
Elapsed Time: 2.604 seconds (Warm-up)
              2.148 seconds (Sampling)
              4.752 seconds (Total)

SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 2).
Iteration: 2000 / 2000 [100%]  (Sampling)
Elapsed Time: 2.554 seconds (Warm-up)
              2.452 seconds (Sampling)
              5.006 seconds (Total)

SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 3).
Iteration: 2000 / 2000 [100%]  (Sampling)
Elapsed Time: 2.764 seconds (Warm-up)
              2.238 seconds (Sampling)
              5.002 seconds (Total)

SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 4).
Iteration: 2000 / 2000 [100%]  (Sampling)
Elapsed Time: 2.553 seconds (Warm-up)
              2.429 seconds (Sampling)
              4.982 seconds (Total)

SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 5).
Iteration: 2000 / 2000 [100%]  (Sampling)
Elapsed Time: 2.585 seconds (Warm-up)
              2.431 seconds (Sampling)
              5.016 seconds (Total)

SAMPLING FOR MODEL 'mod_occupancy_simple' NOW (CHAIN 6).
Iteration: 2000 / 2000 [100%]  (Sampling)
Elapsed Time: 2.606 seconds (Warm-up)
              2.291 seconds (Sampling)
              4.897 seconds (Total)

> out
Inference for Stan model: mod_occupancy_simple.
6 chains, each with iter=2000; warmup=1000; thin=10; 
post-warmup draws per chain=100, total post-warmup draws=600.

           mean se_mean  sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
psi_mean    0.5     0.0 0.0    0.4    0.5    0.5    0.5    0.6   600    1
p_mean      0.7     0.0 0.0    0.6    0.7    0.7    0.7    0.8   600    1
B0          0.0     0.0 0.2   -0.3   -0.1    0.0    0.1    0.3   600    1
B_x1        1.0     0.0 0.2    0.6    0.8    1.0    1.1    1.4   600    1
B_x2       -0.1     0.0 0.2   -0.4   -0.2   -0.1    0.0    0.2   600    1
A0          1.0     0.0 0.2    0.6    0.8    0.9    1.1    1.4   600    1
A_x1       -0.1     0.0 0.2   -0.4   -0.2   -0.1    0.0    0.3   592    1
lp__     -262.5     0.1 1.5 -266.3 -263.3 -262.2 -261.4 -260.6   600    1

Samples were drawn using NUTS(diag_e) at Sun Jan 05 12:08:04 2014.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).