Looking for help with the performance of an occupancy model

[Ryan Batt]

Hi all,

I decided to try Stan because I was starting to find computation time in JAGS to be limiting for the models I was running, and I had plans to increase model complexity even more. I've put together my first Stan model, and I'm a little nervous because the run time is quite long even for a small fraction of the full data set. I've read quite a bit of the manual, and found some wonderfully helpful examples on this mailing list, too.

My model is a minimal example that I think retains some of the essential elements of its complexity. I've tried making some changes to improve performance, but I'm still running rather slowly (taking about 15 minutes for ~1/1E4 of the full data set, and before introducing additional complexities/ parameters). So I was hoping I could get some pointers on where to focus my efforts before I begin hacking away.

The model

The model is similar to what can be found here: https://groups.google.com/forum/#!searchin/stan-users/occupancy/stan-users/Jc0Wej7XCKI/qalSXkVdUacJ

It's also an extension of this paper by Zipkin et al., which extends the previous example in a few ways, most notable of which is by adding covariates.

The presence/ absence data come in a 4-dimensional array, X, with the following dimensions:

1. year (time, t)
2. site (a 1º grid, j)
3. sample (a repeated visit to a site, k)
4. species (s, or n ... sorry in advance for any confusion)
   

There are two main components to the model:

1. Ψ (Psi) is the probability of a species being present
2. Θ (Theta) is the probability of detecting a species, given that it's present

Each of these probabilities is a linear function of covariates and associated coefficients (parameters). The coefficients vary among species, in a manner such that all species-specific parameters are drawn from the same parent distributions.

The observations come from X, but the true presence/ absences are in Z, a latent parameter matrix of 1's and 0's (I used Bob Carpenter's code as an example of how to marginalize this out). There's also a weird parameter, w, which is a binary (1, 0) indicating whether a given species is "real" (available to be sampled).

The model looks something like this:

​

U and V have dimensions K x nV and K x nU, where K is the maximum number of samples taken from any site in any year, and nV and nU are the number of covariates. The parameters

matrices α and β have dimensions of nV x S and nU x S, respectively; the S parameters in each row are drawn from Normal distribution, giving rise to a different parameter for each species. The parameter array  is of length T (number of years), and for each t, it has a length (S) with definition .

In words, the model uses covariates to predict both presence/ absence and detectability. The presence/ absence parameters vary hierarchically among species. The covariates are sample-specific (they are things like "day of year" for detection, and "temperature" for presence/absence). In addition, there's a parameter that determines whether or not a species is even real, and this parameter is drawn from a parent distribution whose sole parameter is a year-specific average probability of any given species being part of the supercommunity.

The Code

I've attached an R script which contains the data, Stan model, and R code to run the model. This data set is a tiny subset of what I'll be working with, but I think it's substantial enough to illustrate some of the challenges with the data.

In particular, I had to use a weird workaround to handle sites that weren't sampled in some years, and the fact that each site-year combination has a different number of samples (K). I still presented the data and covariates in an array, but I turned the NA's to 0's. Then, I have a matrix called "nK" that indicates the number of samples in a given site-year. When an element of nK is 0, it means that that site wasn't sampled (I don't think this happens in the example data set). When looping through the data, I don't loop through the full extent of the covariate/ response arrays, instead stopping the loop at the sampled extent.

It is worth noting that the sites are not arbitrary, as the same site might be sampled repeatedly in different years. Thus, I may introduce site-specific parameters in the future. On the other hand, the distinction between the k=1 and k=K samples is arbitrary. Thus, it might be possible to recode the samples using the ragged array format (although I can't figure out how to do this with the degree of nesting I have), and it could be possible to recode the sites using a data base format, but not the ragged array format.

The Problem

My code runs slowly. I'm pretty sure it's because it isn't coded efficiently. I tried using matrix algebra and vectorized operations when I could, but I ran into trouble here because of the raggedness of the samples/ sites (see above). I also wasn't sure how to code the ragged array examples given that my situation is a fair bit more complex than any examples I could find. I was thinking of reparameterizing the distributions to, e.g., use normal(0,1) for alpha and beta, and re-scale them elsewhere. I didn't do this yet because it is my impression that this will help things converge more quickly, but not run more quickly. Right now it seems that regardless of how fast this code converges, it'll still be kinda slow (I'm probably showing my ignorance here).

I tried to be careful in how I arranged my dimensions, and my choice of containers. I read through the manual, and tried to optimize the order in which I looped through arrays and matrices.

Lastly, I have a hunch that there might be something I did wrong when repurposing some of Bob Carpenter's code. He had written that code for a case where the samples didn't have covariates, so the response variable could be Binomial instead of Bernoulli, and where the parameters (e.g., even Psi and Theta) didn't change among sites, so some of those operations could have been vectorized. Some of the statistical/ coding strategies don't really apply to my data, and while I tried to adapt that code properly, I may have done so in an very inefficient way.

Any help with how to improve this model (in the context of these unwieldy data) would be appreciated. The mode may seem complicated, but it's actually going to get a fair bit more complex. It's also working right now, so maybe it was a fair starting point. Thank you for reading through this description, and I apologize for the errors I undoubtedly made in specifying the model here.

Attached is the code.

Best,

Ryan

[Ryan Batt]

To correct a small Typo, the piece of the code that runs the model should read

# =============
# = Fit Model =
# =============
ebs_msom <- stan(
 model_code=model_file,
 data=staticData,
 control=list(stepsize=0.05, adapt_delta=0.95),
 chains=4, iter=500, refresh=1, seed=1337, cores=4
)

Additionally, a corrected file is attached.

-Ryan

[Bob Carpenter]

The model looks OK, and if you've validated it with fake data, you should
be good to go in terms of what the model does. So what follows is just about
the comptutation.

Do you need to save all those transformed parameters? If
not, you can define them as local variables in the model block.

The bigger speedup gains would come from vectorizing more operations
and/or aggregating observation data.

I think it'd make sense to put Omega on the logit scale and give
it a hierarchical prior (you could do that on the (0,1) scale
with a hierarchical beta).

By transposing alpha and making it an array of vectors, you will avoid
local copies (see the chapter in the manual on arrays vs. vectors)
and you should be able to get away with replacing

for(u in 1:nU){
alpha[u,s] ~ normal(alpha_mu[u], alpha_sd[u]);
}

with

alpha[s] ~ normal(alpha_mu, alpha_sd);

and same for the beta[v,s] loop.

This is redundan

t:


if(tK>0){ // if >0 sampling events, proceed
for(k in 1:tK){// loop through sampling events


if tK <= 0, then the loop will have zero iterations and you save
a level of indentation, some clarity, and a conditional evaluation.

This

1 ~ bernoulli(Omega[t])

is nested a couple levels deep. You could pull it out and then multiply
by the number of times its evaluated.

There's no point in defining a local here:

int K;
K <- nK[t,j];
if(K>0){ // if K is not 0
for(k in 1:K){

This should all be replaced with

for (k in 1:nK[t, j]) { ...

If it makes sense, tighter priors than cauchy(0, 2.5) can help,
either in tails (Say student-t with 4 dof or even normal), or
in reducing the scale from 2.5. I don't know where things liveh
ehre, but these look like logistic regression parameters.

Along the same lines as earlier examples, this loop can be
vectorized:

for(n in 1:N){// loop through observed species
1 ~ bernoulli(Omega[t]); // observed, so available
if(X[t,j,k,n] > 0){ // if present
// increment likelihood
increment_log_prob(
lp_obs(X[t,j,k,n], logit_psi[t,j,k,n], logit_theta[t,j,k,n])

The goal is to reduce the number of calls to the density functions
with vectorization, which in turn reduces the size of the expression graph for
autodiff.

- Bob

P.S. Some minor formatting.

I would get rid of these comments given that the
very next line of code says the same thing.

// ==========================
// = Transformed Parameters =
// ==========================

or things that are totally obvious from looking at the code like

// ============================================
// = Begin Looping down to Point Observations =
// ============================================

For a start, I wouldn't call a string "model_file". But I would
recommend actually pulling the model out and referring to it by path.
That lets you use line-number based diagnostics.

I really don't like 4 indents (see how out of hand the cut-and-pastes
got above) or the quirky R style of rendering loops as
"for(t in 1:nT){" because I find it hard to read, but consistency's
even more important. My main goal is to leave the code uncluttered,
and I find comments like this hurt more than they help:

}else{ // if absent
// increment likelihood
increment_log_prob(
lp_unobs(logit_psi[t,j,k,n], logit_theta[t,j,k,n])
);

You should assume people reading the program know the programming
language, so they know what "else" does and what "increment_log_prob"
does. The comments add clutter and make it hard to read the code
itself.

> On Nov 24, 2015, at 12:47 PM, Ryan Batt <bat...@gmail.com> wrote:
>
> Hi all,
>
> I decided to try Stan because I was starting to find computation time in JAGS to be limiting for the models I was running, and I had plans to increase model complexity even more. I've put together my first Stan model, and I'm a little nervous because the run time is quite long even for a small fraction of the full data set. I've read quite a bit of the manual, and found some wonderfully helpful examples on this mailing list, too.
>
> My model is a minimal example that I think retains some of the essential elements of its complexity. I've tried making some changes to improve performance, but I'm still running rather slowly (taking about 15 minutes for ~1/1E4 of the full data set, and before introducing additional complexities/ parameters). So I was hoping I could get some pointers on where to focus my efforts before I begin hacking away.
>
> The model
>
> The model is similar to what can be found here: https://groups.google.com/forum/#!searchin/stan-users/occupancy/stan-users/Jc0Wej7XCKI/qalSXkVdUacJ
>
> It's also an extension of this paper by Zipkin et al., which extends the previous example in a few ways, most notable of which is by adding covariates.
>
> The presence/ absence data come in a 4-dimensional array, X, with the following dimensions:
> • year (time, t)
> • site (a 1º grid, j)
> • sample (a repeated visit to a site, k)
> • species (s, or n ... sorry in advance for any confusion)
> There are two main components to the model:
> • Ψ (Psi) is the probability of a species being present
> • Θ (Theta) is the probability of detecting a species, given that it's present
> Each of these probabilities is a linear function of covariates and associated coefficients (parameters). The coefficients vary among species, in a manner such that all species-specific parameters are drawn from the same parent distributions.
>
> The observations come from X, but the true presence/ absences are in Z, a latent parameter matrix of 1's and 0's (I used Bob Carpenter's code as an example of how to marginalize this out). There's also a weird parameter, w, which is a binary (1, 0) indicating whether a given species is "real" (available to be sampled).
>
> The model looks something like this:
>
>
>
>
>
>
> ​

> U and V have dimensions K x nV and K x nU, where K is the maximum number of samples taken from any site in any year, and nV and nU are the number of covariates. The parameters

> matrices α and β have dimensions of nV x S and nU x S, respectively; the S parameters in each row are drawn from Normal distribution, giving rise to a different parameter for each species. The parameter array is of length T (number of years), and for each t, it has a length (S) with definition .

>
> In words, the model uses covariates to predict both presence/ absence and detectability. The presence/ absence parameters vary hierarchically among species. The covariates are sample-specific (they are things like "day of year" for detection, and "temperature" for presence/absence). In addition, there's a parameter that determines whether or not a species is even real, and this parameter is drawn from a parent distribution whose sole parameter is a year-specific average probability of any given species being part of the supercommunity.
>
> The Code
> I've attached an R script which contains the data, Stan model, and R code to run the model. This data set is a tiny subset of what I'll be working with, but I think it's substantial enough to illustrate some of the challenges with the data.
>
> In particular, I had to use a weird workaround to handle sites that weren't sampled in some years, and the fact that each site-year combination has a different number of samples (K). I still presented the data and covariates in an array, but I turned the NA's to 0's. Then, I have a matrix called "nK" that indicates the number of samples in a given site-year. When an element of nK is 0, it means that that site wasn't sampled (I don't think this happens in the example data set). When looping through the data, I don't loop through the full extent of the covariate/ response arrays, instead stopping the loop at the sampled extent.
>
> It is worth noting that the sites are not arbitrary, as the same site might be sampled repeatedly in different years. Thus, I may introduce site-specific parameters in the future. On the other hand, the distinction between the k=1 and k=K samples is arbitrary. Thus, it might be possible to recode the samples using the ragged array format (although I can't figure out how to do this with the degree of nesting I have), and it could be possible to recode the sites using a data base format, but not the ragged array format.
>
> The Problem
> My code runs slowly. I'm pretty sure it's because it isn't coded efficiently. I tried using matrix algebra and vectorized operations when I could, but I ran into trouble here because of the raggedness of the samples/ sites (see above). I also wasn't sure how to code the ragged array examples given that my situation is a fair bit more complex than any examples I could find. I was thinking of reparameterizing the distributions to, e.g., use normal(0,1) for alpha and beta, and re-scale them elsewhere. I didn't do this yet because it is my impression that this will help things converge more quickly, but not run more quickly. Right now it seems that regardless of how fast this code converges, it'll still be kinda slow (I'm probably showing my ignorance here).
>
> I tried to be careful in how I arranged my dimensions, and my choice of containers. I read through the manual, and tried to optimize the order in which I looped through arrays and matrices.
>
> Lastly, I have a hunch that there might be something I did wrong when repurposing some of Bob Carpenter's code. He had written that code for a case where the samples didn't have covariates, so the response variable could be Binomial instead of Bernoulli, and where the parameters (e.g., even Psi and Theta) didn't change among sites, so some of those operations could have been vectorized. Some of the statistical/ coding strategies don't really apply to my data, and while I tried to adapt that code properly, I may have done so in an very inefficient way.
>
> Any help with how to improve this model (in the context of these unwieldy data) would be appreciated. The mode may seem complicated, but it's actually going to get a fair bit more complex. It's also working right now, so maybe it was a fair starting point. Thank you for reading through this description, and I apologize for the errors I undoubtedly made in specifying the model here.
>
> Attached is the code.
>
> Best,
>
> Ryan
>
>

> --
> 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <minimal_msomStatic.R>

[Ryan Batt]

Thank you so much for those extensive comments! I've tested single-year versions of this model pretty extensively on simulated data, and it worked pretty well. Like I said, I'm building up complexity, so this exact model hasn't been completely validate (let that be a fair heads-up for others who find this code).

Along those lines, anyone should feel free to use this code/ data. I'm working it into an R package that will go along with a manuscript, and in the package's DESCRIPTION file I've listed it as GPL-3. The package is currently a subdirectory on a private repo, but later will be transferred into its own public repo.

Again, thank you for the comments. Below I've included some responses to your suggestions. I'll make some changes, and post back with results.

-Ryan

On Tuesday, November 24, 2015 at 3:00:32 PM UTC-5, Bob Carpenter wrote:

The model looks OK, and if you've validated it with fake data, you should
be good to go in terms of what the model does.  So what follows is just about
the comptutation.

Do you need to save all those transformed parameters?   If
not, you can define them as local variables in the model block. 

You're referring to `logit_psi` and `logit_theta`, right? I could live without having those saved, but I would be trying recalculate them in R later. Also, when you say I should define them as local variables in the model block, are you suggested that I define temporary variables of lower dimensionality? For example, I currently refer to those transformed parameters 1 element at a time in the model block as `logit_psi[t,j,k,n]`; so instead I could define a `logit_psi` that is overwritten for every t,j,k,n. Is that what you had in mind? I can see how this would save on memory, but it seems like this would reduce efficiency (having to redefine a variable and give it content, as opposed to defining a variable of pre-allocated size then filling), or perhaps be a wash (less indexing is needed, but that's fast in Stan). <shrug> I only have a vague idea of what I'm talking about, whereas you know much more :)

The bigger speedup gains would come from vectorizing more operations
and/or aggregating observation data.   

To your second suggestion there, you're simply saying that maybe I don't have to define the model process to such a fine index. So I could have an average temperature for each sample in a site-year, so instead of having that U matrix indexed as [t,j,k] it could be just [t,j]. That what you mean?

I think it'd make sense to put Omega on the logit scale and give
it a hierarchical prior (you could do that on the (0,1) scale
with a hierarchical beta).

Due to methodological concerns, I wanted to be a bit conservative with how I defined the relationship of parameters among years. That aside, is it generally true that having 3 parameters drawn from a hierarchical prior will be more efficient than having 3 parameters defined 'independently' of each other?

So I could do

\\ in parameter block
Omega_mu ~ beta(2,2);
Omega_sd ~ cauchy(0, 1);
for (t in 1:nT) {
  Omega_logit[t] ~ normal(logit(Omega_mu), Omega_sd);
}

​

I get the feeling that what I just wrote is wrong, but I'm a bit confused by the convention you're suggesting.

 

By transposing alpha and making it an array of vectors, you will avoid
local copies (see the chapter in the manual on arrays vs. vectors)
and you should be able to get away with replacing

                for(u in 1:nU){
                        alpha[u,s] ~ normal(alpha_mu[u], alpha_sd[u]);
                }

with  

  alpha[s] ~ normal(alpha_mu, alpha_sd);

and same for the beta[v,s] loop.

I'd had a similar thought myself. I originally wrote it that way so I could use matrix algebra here:

    for (t in 1:nT) {
        for (j in 1:Jmax) {
            int K;
            K <- nK[t,j];
            if (K>0) { // if K is not 0
                // if the K matrix gets too large, 
                // it is in here that I should add
                // a loop, wherein I subset to U[t,j,k], and
                // vectorize over spp only
                for (k in 1:K) {
                    logit_psi[t,j,k] <- U[t,j,k]*alpha;
                    logit_theta[t,j,k] <- V[t,j,k]*beta;
                }

            }
        }
    }

​

As you can see from my comments, I changed my mind and decided to do the looping. Originally I wanted `U[t,j] * alpha;`. The downside is that with the potential for many missing K, I wasn't sure the matrix algebra would be a savings. 

If I wanted to go back to the matrix algebra, I suppose I could do what you suggest, then also transpose the matrix portion

of U and V (currently declared as matrix[Kmax,nU] U[nT,Jmax], switch to matrix[nU,Kmax] U[nT,Jmax]), then switch toalpha * U[t,j. I'd have to check to make sure that'd work out.

Regardless, I remember thinking that I couldn't keep the current dimensions and vectorize across the column (in R, e.g., do a alpha[,1] <-  statement), so I understand why the transposition is needed for this vectorization. I'll try it.

 

This is redundan

t:


                        if(tK>0){ // if >0 sampling events, proceed
                                for(k in 1:tK){// loop through sampling events
                                        

if tK <= 0, then the loop will have zero iterations and you save
a level of indentation, some clarity, and a conditional evaluation.

Oh, OK — so if I had for (k in 1:0) {, nothing would happen. That sounds familiar from something I read in the manual, because Stan doesn’t loop backwards. Gotcha. I hadn’t thought of the implications of the rule in this context, thank you.

​

 

This

   1 ~ bernoulli(Omega[t])

is nested a couple levels deep.  You could pull it out and then multiply
by the number of times its evaluated.  

For multiplying out the probability increments, can I still use the ~ notation? Or would I just do something like 1 ~ bernoulli(Omega[t]) * N * sum(nK[t])? Note that nK[t] is a row_vector, and if that notation works, I could move this out of 3 for loops. Good call on moving this; however it’s done, it’ll help.

​

There's no point in defining a local here:

                        int K;
                        K <- nK[t,j];
                        if(K>0){ // if K is not 0
                                for(k in 1:K){

This should all be replaced with

                        for (k in 1:nK[t, j]) { ...

Yep, same as with the previous example that used tK. Thanks again.

 

If it makes sense, tighter priors than cauchy(0, 2.5) can help,
either in tails (Say student-t with 4 dof or even normal), or
in reducing the scale from 2.5.  I don't know where things liveh
ehre, but these look like logistic regression parameters.

I'm unclear about the implications of your last sentence. In your example occupancy model (linked in the original message) you used cauchy on similar parameters. So that's why I used them, but there might be some context I'm missing as to why it's a good idea in the example, but not in my model. But I'll trying reducing tail thickness. 

Along the same lines as earlier examples, this loop can be
vectorized:

                                        for(n in 1:N){// loop through observed species
                                                1 ~ bernoulli(Omega[t]); // observed, so available
                                                if(X[t,j,k,n] > 0){ // if present
                                                        // increment likelihood
                                                        increment_log_prob(
                                                                lp_obs(X[t,j,k,n], logit_psi[t,j,k,n], logit_theta[t,j,k,n])

The goal is to reduce the number of calls to the density functions
with vectorization, which in turn reduces the size of the expression graph for
autodiff.

I think I’m running into at least 2 potential problems vectorizing these pieces. First, the lp_obs function can only be applied to the first N elements of that fourth dimenion; i.e., to X[t,j,k, 1:N] … which is R syntax, and which I am unsure of how to do in Stan. On page 353 of the manual I see the vector and array slicing operations, but I didn’t understand them right away. Taking a second look, I see that maybe segment(X[t,j,k], 1, N) would give the desired result for vectorization. I’ve read elsewhere that in many cases the copy command induced by segment is a small loss in performance compared to the performance gained by vectorizing. So I see how this could be handy.

Second, this entire process could be greatly simplified if I could condense the lp_obs and lp_unobs functions into 1. These are your functions, which I simplified to my needs. After simplification, I didn’t see how to easily condense them into 1 function. It’s because of that log1m_inv_logit(logit_psi) term in the lp_unobs function —- I have to cancel it out if x is 1, and I can’t think of a way to do that, unless something like value * (x==0) can take advantage of TRUE being 1, and FALSE being 0, and then working for algebra (like in R). Any ideas?

​

 

- Bob

P.S.  Some minor formatting.  

I would get rid of these comments given that the
very next line of code says the same thing.

// ==========================
// = Transformed Parameters =
// ==========================

My text editor doesn't have syntax highlighting for Stan (TextMate2), so I'm current using the C++ highlighting (I saw this: https://github.com/dougalsutherland/sublime-stan, but unfortunately TM2 only recently made it easy to export entire bundles, and this repo only has certain grammars, and I haven't gotten around to manually installing them). All of this is to say that I still need obvious visual pointers in the code. But point taken.

or things that are totally obvious from looking at the code like

        // ============================================
        // = Begin Looping down to Point Observations =
        // ============================================

For a start, I wouldn't call a string "model_file".  But I would
recommend actually pulling the model out and referring to it by path.
That lets you use line-number based diagnostics.

Ah, yes --- that what I did originally, so "model_file" is left over from when it was actually just a path. My copy is saved as a .stan file. I just put it all in 1 place for convenience.  

I really don't like 4 indents (see how out of hand the cut-and-pastes
got above) or the quirky R style of rendering loops as
"for(t in 1:nT){" because I find it hard to read, but consistency's
even more important.  My main goal is to leave the code uncluttered,
and I find comments like this hurt more than they help:

                                                }else{ // if absent
                                                        // increment likelihood
                                                        increment_log_prob(
                                                                lp_unobs(logit_psi[t,j,k,n], logit_theta[t,j,k,n])
                                                        );

My apologies for this -- I actually read the manual section on style, and appreciated it a lot. Good to know the convention. I had intended to sweep out tab for spaces, but I forgot. As for the loops, that's a habit that'll be hard to reverse, but maybe in time :) And I agree about my comments. Thanks for the pointers here.

 

You should assume people reading the program know the programming
language, so they know what "else" does and what "increment_log_prob"
does.  The comments add clutter and make it hard to read the code
itself. 

Interesting philosophy on reading the code. Being new to the language, I go into a more detailed mode of having the comments almost serve as a syntax tutorial. I usually pull that back quite a bit after my first few lines, once I get the hang of it. But sometimes I leave things in, because I try to share my code with other people who I know would like to run and understand it, but definitely don't know the language. So I sometimes adopt the philosophy of "if this code is the first time someone sees this language, is it going to be too confusing?". But I'm not consistent in that regard, so maybe I should adopt your philosophy.

Thanks again for all of the comments. I'll start making some of these changes, and testing the effects. I'll post more soon. It'll be fun to see the changes in performance! 

[Maxwell Joseph]

TFPU - always glad to see more occupancy models here. Just a quick comment on using segment() in-line below:

It looks like log_inv_logit(logit_psi) + bernoulli_logit_log(x, logit_theta) is incremented for each detection, and log_sum_exp(lp_obs(0, logit_psi, logit_theta), log1m_inv_logit(logit_psi)) for each nondetection. But, the marginalization over z I think you're looking for is: 

If species n seen at time t, site, j:

so that log(\psi_{ntj}) gets incremented once instead of k_{jt} times, and

If species n not seen at time t, site j

where again the occupancy parameters are incremented only once per combination of n, t, and j.

With the current formulation, there are potentially instances where species n was seen in time t at site j, such that X[t, j, , n] > 0, but there were some surveys (values of k) with X[t, j, k, n] = 0. In these cases, if we assume that there are no false-detections, then the species is known to be present and any zeros must come from the detection process, so the component of the likelihood that corresponds to species absence log1m_inv_logit(logit_psi) should probably not be incremented. In Bob's original example, the detection histories are modeled as binomial rather than Bernoulli observations, resulting in the correct marginalization over z. 

I think you can get around this either by aggregating your Bernoulli observations into binomial form, or by using segment() to feed a vector of observations X[t, j, , n] in R speak, into your likelihood functions (the log version of the above likelihood expressions). In conventional occupancy model formulations, the occupancy states and probabilities do not vary among surveys (e.g., see top of page 481 in the Zipkin paper: \psi is not indexed by k, but p is). So in your case \psi should be indexed by t, j, and n, but not k. Hope that makes sense!

[Ryan Batt]

Maxwell,

Thanks for those comments. You've homed in on the big hangup here. The covariates change among the K. So there aren't K replicates, there are K sampling locations within a site (j).

To put it in ecological context, we're looking at marine critters in the Eastern Berring Sea. A "site" is a 1º grid (100+ km^2). A "sample" (k) is a tow of the net by a boat. This tow is very small compared to the size of the site, and the tows within the site can occur months apart in the same year. So tows can be dozens of kilometers and scores of days apart. However, given the size of the region, and given that there's almost 100 sites in the area, the different K are still really close together compared to the full context.

All this is to say that you could argue that the K are different (not really replicates, have different covariate values [both for detection and presence]), but the samples within a site are still more similar than among sites.

Does this change your perspective? It's why I didn't use a binomial. But I'm considering taking the average of the covariates within a site-year, then using binomial. But as it stands, having a species present for 1 K in a site-year doesn't mean that all the other K truly have the species present, too. Necessarily. See where I'm coming from?

-Ryan

[Bob Carpenter]

> On Nov 24, 2015, at 4:47 PM, Ryan Batt <bat...@gmail.com> wrote:
>
> Thank you so much for those extensive comments!

I'm very interested in seeing more ecology applications. More replies below.

> I've tested single-year versions of this model pretty extensively on simulated data, and it worked pretty well. Like I said, I'm building up complexity, so this exact model hasn't been completely validate (let that be a fair heads-up for others who find this code).
>
> Along those lines, anyone should feel free to use this code/ data. I'm working it into an R package that will go along with a manuscript, and in the package's DESCRIPTION file I've listed it as GPL-3. The package is currently a subdirectory on a private repo, but later will be transferred into its own public repo.

Thanks.

> Again, thank you for the comments. Below I've included some responses to your suggestions. I'll make some changes, and post back with results.
>
> -Ryan
>
> On Tuesday, November 24, 2015 at 3:00:32 PM UTC-5, Bob Carpenter wrote:
> The model looks OK, and if you've validated it with fake data, you should
> be good to go in terms of what the model does. So what follows is just about
> the comptutation.
>
> Do you need to save all those transformed parameters? If
> not, you can define them as local variables in the model block.

> You're referring to `logit_psi` and `logit_theta`, right? I could live without having those saved, but I would be trying recalculate them in R later.

OK, then you probably want to save them in terms of total work
required.

> Also, when you say I should define them as local variables in the model block, are you suggested that I define temporary variables of lower dimensionality?

I just meant transfer all the code, but you could do that, too, and
it'd probably be faster.

> For example, I currently refer to those transformed parameters 1 element at a time in the model block as `logit_psi[t,j,k,n]`; so instead I could define a `logit_psi` that is overwritten for every t,j,k,n. Is that what you had in mind? I can see how this would save on memory, but it seems like this would reduce efficiency (having to redefine a variable and give it content, as opposed to defining a variable of pre-allocated size then filling), or perhaps be a wash (less indexing is needed, but that's fast in Stan). <shrug> I only have a vague idea of what I'm talking about, whereas you know much more :)

There's a chapter in the manual on arrays and vectors
and efficiency where I try to explain what's going on.
The only time it makes sense to define something on the
outside is if you're going to reuse it. Otherwise, it
should be local for both efficiency and clarity.

> The bigger speedup gains would come from vectorizing more operations
> and/or aggregating observation data.

> To your second suggestion there, you're simply saying that maybe I don't have to define the model process to such a fine index. So I could have an average temperature for each sample in a site-year, so instead of having that U matrix indexed as [t,j,k] it could be just [t,j]. That what you mean?

Yes, that's one way to do it. Sometimes that's exact
depending on your model. For instance, if it's normal,
you can use the sufficient stats. But we don't have
the sufficient stats distribution (though someone took
a partial stab at it in the past). An easier case is
converting bernoullis to binomials, replacing

for (i in 1:I)
y[i] ~ bernoulli(theta);

with

sum(y) ~ binomial(size(y), theta);

It's exactly the same up to a proportion but much more efficient.

Or, given that the variable's saved, you could do something like
this, which is more general:

increment_log_prob(sum(y) * bernoulli_log(1, theta));
increment_log_prob((size(y) - sum(y)) * bernoulli_log(0, theta));

> I think it'd make sense to put Omega on the logit scale and give
> it a hierarchical prior (you could do that on the (0,1) scale
> with a hierarchical beta).

> Due to methodological concerns, I wanted to be a bit conservative with how I defined the relationship of parameters among years. That aside, is it generally true that having 3 parameters drawn from a hierarchical prior will be more efficient than having 3 parameters defined 'independently' of each other?

It's a result of what Andrew calls the "Folk Theorem". If
you define a hierarchical model, then it figures out how
much partial pooling to do, which improves estimates, and
as a pleasant side effect, more highly constrains the sampler
so that efficiency is better, too. At 3 members per group,
it's dicy unless you have another layer of hierarchy, because
it's hard to estimate anything with 3 data points.

> So I could do
> \\ in parameter block
> Omega_mu ~ beta(2,2);
> Omega_sd ~ cauchy(0, 1);
> for (t in 1:nT) {
> Omega_logit[t] ~ normal(logit(Omega_mu), Omega_sd);
> }
>
> ​
> I get the feeling that what I just wrote is wrong, but I'm a bit confused by the convention you're suggesting.

What you wrote is one way to do it. I'd have probably stuck
to either the logit scale (note that I'm doing this in a
vectorized way):

real mu_Omega;
real<lower=0> sigma_Omega;
vector[?} Omega_logit;


mu_Omega ~ normal(0, 5);
sigma_Omega ~ cauchy(0, 2.5);
Omega_logit ~ normal(mu_Omega, sigma_Omega);

or on the (0, 1) scale:

real<lower=0, upper=1> theta_Omega;
real<lower=0.1> count_Omega;

Omega ~ beta(theta_Omega * count_Omega, (1 - theta_Omega) * count_Omega);

and then you need priors on theta_Omega and on count_Omega (the
manual suggests uniform or a light beta on theta_Omega and a pareto(0.1, 1.5)
on count_Omega, though I'd probably recommend pareto(2, 1.5) and constraining
count_Omega > 2.

As to the rest, we really need sparse matrix operations. It's on
my to-do list, but it's a huge huge project at this point given everything
else going on with Stan (in particular, all the linear algebra on the
one side and all the I/O in the interfaces and indxing in the language).
> ...

> For multiplying out the probability increments, can I still use the ~ notation? Or would I just do something like 1 ~ bernoulli(Omega[t]) * N * sum(nK[t])? Note that nK[t] is a row_vector, and if that notation works, I could move this out of 3 for loops. Good call on moving this; however it’s done, it’ll help.

See the example above using increment_log_prob.

>
> ​
>
> There's no point in defining a local here:
>
> int K;
> K <- nK[t,j];
> if(K>0){ // if K is not 0
> for(k in 1:K){
>
> This should all be replaced with
>
> for (k in 1:nK[t, j]) { ...
> Yep, same as with the previous example that used tK. Thanks again.
>
>
>
> If it makes sense, tighter priors than cauchy(0, 2.5) can help,
> either in tails (Say student-t with 4 dof or even normal), or
> in reducing the scale from 2.5. I don't know where things liveh
> ehre, but these look like logistic regression parameters.

> I'm unclear about the implications of your last sentence. In your example occupancy model (linked in the original message) you used cauchy on similar parameters. So that's why I used them, but there might be some context I'm missing as to why it's a good idea in the example, but not in my model. But I'll trying reducing tail thickness.

We're moving over time to recommending stronger and stronger
priors on both statistical and computational grounds (luckily,
they tend to move together, hence the "Folk Theorem").

> Along the same lines as earlier examples, this loop can be
> vectorized:
>
> for(n in 1:N){// loop through observed species
> 1 ~ bernoulli(Omega[t]); // observed, so available
> if(X[t,j,k,n] > 0){ // if present
> // increment likelihood
> increment_log_prob(
> lp_obs(X[t,j,k,n], logit_psi[t,j,k,n], logit_theta[t,j,k,n])
>
> The goal is to reduce the number of calls to the density functions
> with vectorization, which in turn reduces the size of the expression graph for
> autodiff.
> I think I’m running into at least 2 potential problems vectorizing these pieces. First, the lp_obs function can only be applied to the first N elements of that fourth dimenion; i.e., to X[t,j,k, 1:N] … which is R syntax, and which I am unsure of how to do in Stan. On page 353 of the manual I see the vector and array slicing operations, but I didn’t understand them right away. Taking a second look, I see that maybe segment(X[t,j,k], 1, N) would give the desired result for vectorization. I’ve read elsewhere that in many cases the copy command induced by segment is a small loss in performance compared to the performance gained by vectorizing. So I see how this could be handy.

Better slicing is coming in our next release, but we'll keep
segment around as Maxwell Joseph suggested using.

>
> Second, this entire process could be greatly simplified if I could condense the lp_obs and lp_unobs functions into 1. These are your functions, which I simplified to my needs. After simplification, I didn’t see how to easily condense them into 1 function. It’s because of that log1m_inv_logit(logit_psi) term in the lp_unobs function —- I have to cancel it out if x is 1, and I can’t think of a way to do that, unless something like value * (x==0) can take advantage of TRUE being 1, and FALSE being 0, and then working for algebra (like in R). Any ideas?

Algebra (including collapsing loops into sums and products)
can often greatly help reduce computation.

> Interesting philosophy on reading the code. Being new to the language, I go into a more detailed mode of having the comments almost serve as a syntax tutorial. I usually pull that back quite a bit after my first few lines, once I get the hang of it. But sometimes I leave things in, because I try to share my code with other people who I know would like to run and understand it, but definitely don't know the language. So I sometimes adopt the philosophy of "if this code is the first time someone sees this language, is it going to be too confusing?". But I'm not consistent in that regard, so maybe I should adopt your philosophy.

Not just mine. Almost everyone comes to this conclusion eventually.
It's even in the Google Coding style guidelines :-)

https://google.github.io/styleguide/cppguide.html#Implementation_Comments

Don'ts

Do not state the obvious. In particular, don't literally describe what
code does, unless the behavior is nonobvious to a reader who understands
C++ well. Instead, provide higher level comments that describe why the
code does what it does, or make the code self describing.

The problem with not adopting this philosophy is that you never
know where to stop.

> Thanks again for all of the comments. I'll start making some of these changes, and testing the effects. I'll post more soon. It'll be fun to see the changes in performance!

It's probably not going to be huge if most of the work
is down in the lower-level densities. But vectorizations
of distributions with shared parameters can speed things
up an order of magnitude, and computing sufficient statistics
can speed things up immensely, even to the point of making
some computations analytical.

I wish our compiler was smart enough to do a lot of this
automatically, but that's a huge challenge to even begin to
undertake.

- Bob

[Maxwell Joseph]

Hi Ryan, 

Thanks for clarifying! I assumed that the surveys were replicated within the sites a la traditional MacKenzie et al. occupancy setups. Given the large spatial scale & lack of boundaries between the "sites" and the temporal/spatial distance between tows, I see why occupancy is a tow rather than site level state variable. Seems alright so long as the occupancy and detection parameters are still identifiable (sounds like they are based on your simulations, which is great). 

Cheers,

Max

P.S. I'd upvote Bob's recommendation to tighten up some of the cauchy priors for the occupancy and detection coefficients. These can translate to U-shaped priors on the probability scale, which may or may not be desirable!

[Ryan Batt]

Hi all,

I've been working to make as many improvements as I could. Thanks to Max and Bob for some insights.

I boiled Bob's original suite of comments down to a to-do list of optimizations to try. I re-wrote the call to `stan()` to only do 50 iterations for each chain, so that the code would run quickly for testing. The attached R file includes all changes I've made since the original posting. At the end of the code I've included the calculation of metrics that I used to measure speed/ efficiency.

Before I made any changes, the code ran in 33.9 sec. These timings are generally from only 1 trial run. Typically stable within a second or so. Does not include compilation time (see attached code).

Please note: A 1 second savings here could end up saving much more time (hours? days?) when the model is run for more chain iterations (e.g., 5k instead of 50) and for the full data set (but making this extrapolation to full data set is tricky). Currently the array containing species abundance is 3x3x7x15 (for these tests). The full data set is 32x99x13x466 + (that last dimension gets padded with, say, 100+ never-seen species) (and I have similar data sets for 11 regions). So even a little savings in these tests will make a big difference in the end.

Here are the list of changes suggested, in the order that I made the changes, and the associated timings:

1. Don't use if() statements to avoid going into a loop when K=0; `for (k in 1:0) {` will skip the loop anyway.
• 34.4 sec 
• (I hadn't started doing seconds/ effective sample size yet)
2. Move `1 ~ Bernoulli(Omega[t])` out a few levels of for() loops, using `increment_log_prob(bernoulli_log(1, Omega[t]) * N * sum(nK[t]))` instead
• 33.4 sec 
• (didn't calculate efficiency)
3. vectorize `alpha[u,s] ~ normal(` and similar for beta; Thought this would require transpose, but it didn't
   
• 31.6 sec
• 0.89 sec per sample (denominator is average nEff across all parameters; see attached)
4. use tighter priors than cauchy(0, 2.5); I did cauchy(0,2) on the sd, and cauchy(0,1) on the mu; originally tried cauchy(0,1) for both, but that was slower/ less efficient than current
• 32.8 sec
• 0.87 sec per sample
5. Make Psi and Theta local variables that are no longer tracked or reported in end (I can recalculate on my own if needed)
• 31.6 sec
• 0.52 sec per sample
6. Use matrix algebra for calculating each Psi and Theta (no transpose needed). Used `block()` to adapt the size of covariate matrix (U or V) to the number of K in that site-year, so not wasting computations for when `Kmax` is much larger than most `nK[,]` (in some cases there are ~100 K, but only for 1 site-year, e.g.). Avoids many empty elements of Psi and Theta, too.
• 29.0 sec
• 0.48 sec per sample
7. Try vectorizing the calculation of log probabilities. To do this, I decided to aggregate my process covariates (U) over the sample (K) dimension. Then I tried creating vectorized form of the `lp_obs()` function (`lp_obs2()`), and created a new `lp_exists()` function that can works in both of the situations covered by `lp_obs()` and `lp_unobs()`. I struggled with this step, mostly due to matching of types and those cascading effects. I had made some intermediate changes when working towards this that I didn't save, but were a fair bit faster, but didn't include my attempt to vectorize the `lp_obs()` functions. But I forget exactly what I'd had.
1. 22.9 sec
2. 0.72 sec per sample
8. I have not implemented this yet: Make Omega hierarchical, and/or put it on logit scale (see Bob comments)

As you can see, things have improved. I struggled a lot with #7, and as a result I am hesitant to attempt a similar vectorization for the never-seen species.

As a side note, I found it cumbersome that I could not store integers in a matrix, because Bernoulli requires integer values to be supplied. I found my work arounds, but it was trying. Along those lines, many of the container formats I chose were the result of iterative trial and error; it was hard to get everything to match up correctly for the functions, but I'm learning.

I would greatly appreciate any further comments, but I understand that I may be asking a lot, and I'm very happy with the feedback I've gotten thus far. Thank you all.

-Ryan

[Ryan Batt]

Hi all,

Small update on timings.

I decided to let this model chug on a larger-sized data set (32x99x13x466 + 10 0's , so 476). Same model and stan call as in the previously-attached R script. It's been running for ~12 hours now, and the chains are 44%, 52%, 30%, and 28% done (of only 50 iterations; I'm running this on 4 cores, 1 for each chain).

I think that having 31 year-specific dummy variables in V, and 32 non-hierarchical Omega could be hurting things a bit. But this still seems a lot slower than what I would have hoped for. If it takes 24 hours for this model to run, I guess for 5000 iterations it would take a cool .... 100 days? Ouch.

I should also point out that the last dimension there (the 466) is the number of species. Most of them are observed only a handful of times, so this is quite sparse.

I point out these timings just to emphasize my motivation for speeding the code up. Suggestions are welcome.

Oh, and happy Thanksgiving everyone!

Ryan

...

[Bob Carpenter]

What's the effective sampling rate and how fast is convergence?
You probably don't need 5000 draws for inference.

It would be a huge speedup if rather than encoding sparse matrices,
you used a long form of the data. So if you have a sparse matrix x,
the data structure is a table like this:

x[m, n] | m | n

where first column is value, and second two columns are indices. Then
you need the number of entries. If you can recast the data in this
form, it can be a huge speedup.

The other thing to try other than optimizing the code is to
play with the parameterization. You can use a non-centered
parameterization where instead of

alpha ~ normal(mu, sigma);

you use

alpha_raw ~ normal(0, 1);

alpha <- sigma * alpha_raw + mu;

and then use alpha as before. If you have a hierarchical
model with small counts in the groups, this can also have
a very large effect, in some cases leading to models that
will converge where they wouldn't with the centered
parameterization.

You may also find that lowering the initial step size will
help it converge faster during warmup and increasing the target
acceptance rate may help with effective sample size per iteration.
Both cause slower iterations.

- Bob

> On Nov 26, 2015, at 10:26 AM, Ryan Batt <bat...@gmail.com> wrote:
>
> Hi all,
>

> so that log(\psi_{ntj}) gets incremented once instead of k_{jt} times, and
>
> If species n not seen at time t, site j
>

[Ryan Batt]

The effective sample size varies a lot among parameters. With the 50 iterations, some parameters have hundreds, while others are in the single digits. I could provide a more detailed summary of results, though.

From my initial tests, some parameters seemed to converge rather quickly. But I didn't look at all of them, just a smattering.

I'm honestly pretty confused by how to take advantage of the database representation that you suggest. Am I supposed to drop the 0's from that representation? Including the 0's definitely means that the object itself will take up much more memory, as it will repeat the information  contained in the array's dimension names. On the other hand, I'm less certain about what to do with a database representation that doesn't have all of the year-site-sample-species combinations encoded in the identifying columns. I remember the section in the manual. I remember thinking it would be challenging to code up correctly, but do-able. I would appreciate some elaboration on the topic, if you don't mind. Particularly, you seem very confident in the speed increase, but I don't understand why this would happen. My lack of understanding is probably because I don't have a clear vision of the desired implementation. Sorry for my slowness here.

I'll try the other suggestions about non-centered distributions and tweaking function arguments. Thank you again, very much.

I'll post the results from the full data set run when it's finished. Just in terms of convergence, etc. Right now the chains are between 40% and 90% complete. That's another question: why are the chains showing such different run times? This question clearly illustrates my lack of understanding of what's going on under the hood. I'd enjoy reading more about what leads to the chains sampling at such different rates.

Thanks!!

-Ryan

[Bob Carpenter]

Since you said you had *sparse* data (I take "sparse" to mean
a matrix with lots of zeroes that then get multiplying by other
things to produce more zeros), it led me to *guess*
that you could do it more efficiently than coding it up as
matrices. The problem with sparse matrix operations using
dense matrices is all the multiply-by-zero and add-zero operations.
Each of these adds to the expression graph, which then takes
time during automatic differentiation to compute gradients.

We're getting far too much mail now for me to have time to
go through everyone's code closely. Sorry about that.

The difference in chain speeds is almost always due to failures
to find the typical set (volume of high posterior mass) efficiently
during warmup. The chain will wander off into an area
of high curvature and have to take the step size down to a very
small value to move, which then requires lots of steps elsewhere.
The real problem is that in some cases, Stan can't find its
way back to the typical set (volume of high posterior mass).
You can usually see this in the marginal trace plots of the parameters.

- Bob

[Ryan Batt]

No worries at all, Bob. You've been quite gracious with your time and advice! I am very appreciative.

Your explanations have given me some insight. I'll a couple other strategies for optimizing/ reformulating the code.

For the (potential) benefit of others who may stumble upon this post, I'll be sure to update my code/ strategies. Eventually when this analysis is ready for publication, I'll share more content as well.

Thanks again.

-Ryan

[Michael Betancourt]

The non-centered parameterization is likely going to be a huge

effect on performance so I’d focus on that first.  And remember

that the ultimate measure of performance is ESS / time.

[Ryan Batt]

Thank you, Michael. I think that's a step I understand. I also intend to spend some time trying to vectorize those log prob functions. That second step is tricky b/c there are integer parameters that I have to marginalize over, so I'm glad you mentioned to not do that first.

More soon.

Thanks again, sincerely.

-Ryan

...

[Ryan Batt]

Update:

Doing the non-centered parameterization was a massive speed-up (~20-30 fold).
I noticed that I made an error in my most recent version of the model, after aggregating the presence process to the site level (still incrementing log_inv_logit(logit_psi) for each sample, instead of just 1 time per site).

For a similar reason, I think I also needed to change increment_log_prob(bernoulli_log(1, Omega[t]) * N * sum(nK[t])); to increment_log_prob(bernoulli_log(1, Omega[t]) * N * sum(nKJt]));; i.e., do 1~Omega[t] J*T times, not also incrementing for each sample K in each site J.

I’m currently working on correcting the above mistakes. I’m also going to aggregate the presence covariates, but then model them as ~normal(mean(covariate), sd(covariate)) so that I don’t lose all of the variability when I aggregate.

​

[Ryan Batt]

Hi all,

I've been working on the model. I think aspects of my original re-purposing of Bob's code were flawed. I believe I've fixed most of those errors, but now I'm struggling with how I can vectorize my loops. Without such vectorization, things are still running fairly slow.

I'm attaching the code I have right now, which is working. I've gained a lot of speed from past suggestions. I'm still working on the vectorization. In the Stan model I've commented out large section in the model block -- those sections are my failed attempts at vectorizing the preceding for() loops. The functions are all for the vectorized form, and I don't currently use them in the for() loops.

The whole problem for me can probably be summarized as the need to reformulate the following functions (Bob's) for use when psi and theta not just among species, but among species and sites:

real lp_observed(int x, int K, real logit_psi, real logit_theta) {
    return log_inv_logit(logit_psi) + binomial_logit_log(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;
    real lp_available;
    lp_unavailable <- bernoulli_log(0, Omega);
    lp_available <- bernoulli_log(1, Omega) + J * lp_unobserved(K, logit_psi, logit_theta);
    return log_sum_exp(lp_unavailable, lp_available);
}

​

Terms involving Omega should only be incremented once per species, I believe. So the first term in lp_available <- bernoulli_log(1, Omega) + J * lp_unobserved(K, logit_psi, logit_theta); must be separated from the second, and the second has to be first be summed among all T,J for a species before being added to the first term. That's why I use the lp_available_pt1 vector in my current iterative form of the Never Observed species.

I should mention that I think my vectorized version is wrong because I'm getting a different answer with the same seed. When I initially said I'd tested my model, I didn't mean the exact Stan encoding per se, but rather the general structure (and I'd tested what I had written in JAGS). Things have been reformulated quite a bit now, so it's less clear what's right or wrong.

I'll post back if I figure something out, but I've been struggling with this lately. 

[Bob Carpenter]

Unless you can get things into vectorized probability functions
or simple matrix operations, there's not much else you can do to
vectorize.

The other alternative is to compute sufficient statistics and/or
turn things like multiple Bernoullis into binomials.

Another avenue of optimization is to reuse shared expressions.
It's even worthwhile to build data structures up to cache the
intermediate results if they can be reused. For instance, if
there are two random effects with M and N levels, building
the M x N combined table and indexing into that rather than adding
each time.

If you can post a currently working version of the model, I might
have more concrete suggestions.

- Bob

[Ryan Batt]

Thanks Bob. The current model produces “reasonable” results when fed larger data sets (5 species isn’t very many!). I’ll work on adapting old methods for verifying these models.

If I have a statement like bernoulli_log(1, Omega)*D, is it faster/ correct to switch it to binomial_log(D, D, Omega) ? Similar question for increment_log_prob(bernoulli_log(1, Omega) * N).

I’ll report back once I’ve verified my current Stan model. Thank you so much.

​

-Ryan

[Bob Carpenter]

> On Nov 30, 2015, at 9:24 PM, Ryan Batt <bat...@gmail.com> wrote:
>
> Thanks Bob. The current model produces “reasonable” results when fed larger data sets (5 species isn’t very many!). I’ll work on adapting old methods for verifying these models.
>
> If I have a statement like bernoulli_log(1, Omega)*D, is it faster/ correct to switch it to binomial_log(D, D, Omega) ?

That'll actually be slower because binomial_log() will compute
the normalizing constant for the binomial.

> Similar question for increment_log_prob(bernoulli_log(1, Omega) * N).

That's the speediest approach.

- Bob

[Ryan Batt]

I ran some simulation tests that show pretty good correlations between estimated and true Psi (using posterior means). It's a little off when Psi is low (less than 20%, say, the model overestimates). The model is quite data-hungry, but it seems to be giving good estimates.

[Michael Betancourt]

You really want to compare posterior intervals to the true values in order to

incorporate posterior uncertainty — for example see how often the true value 

is, say, in between the 5% an 95% intervals.

[Ryan Batt]

Thanks for the suggestion, Michael. I'm working on this right now. In short, only about half of the parameters fall within the middle 90% of the posteriors. I'm currently working to figure out if the simulated data isn't really appropriate, or if the code is wrong. My hunch is that the simulation isn't providing enough context for the model, but I'm unsure. Just wanted to give an update.

[Ryan Batt]

OK, I've verified the model. I wasn't getting great results before because of a bug in a data-input formatting function I wrote. Oops.

Sorry this took so long, I wanted to do my best to make sure I could honestly say the model is correct and working.

I'm attaching an updated, minimal, all-in-one R script (contains data, Stan model, R code for running and timing). The data are real, but very small. If desired, I can share a larger data set of simulated data where we'll know the answer. I'd be happy to share some of the code I used to verify the model (although I'd have to clean that code up a bit, just to make it readable for others).

I mentioned this in my first message, but I'd like to reiterate that although I believe the results of this model will be ecologically meaningful, I intend to make 2 'enhancements' to the model:

1. year- and/or site-varying parameters (hierarchical). This depends on how good the covariates are at pulling out space-time differences in occupancy
2. make the model of a dynamic flavor, whereby colonization and persistence terms are modelled. I know that Bob has graciously shared a simple version of such a model in the past.

One thing I've noticed during my testing is that I don't have to run a Stan model for very many iterations to get close to reasonable results. Honestly, with simulated data and this particular model, I'm finding that ~50 iterations gets be in the ballpark. In fact, I can usually figure out if the model is 'working' if I do ~20. That's orders of magnitude fewer iterations than I would do in JAGS for similar feel of accuracy. 

Is there any standard expectation for the number of iterations I should run the Stan model? If I put in a paper that I ran my model for 200 iterations, a reviewer used to BUGS or JAGS would probably question the choice (if I didn't explain effective sample size, that I'm using Stan, etc). What should I consider to be in the ballpark of 'normal' for number of iterations? this point is germane to the performance of the occupancy model, because "performance" is really about "how long it takes for the model to finish", which obviously depends on # iterations. As Bob briefly alluded to, and is implied by extensive discussions here and elsewhere about "effective sample size", perhaps the easiest way to improve my definition of "performance" is to run the model for fewer iterations.

I would certainly appreciate more feedback or tips on how to improve the speediness/ accuracy (!)/ performance of this model. But in the past ~10 days (since original post), I've really become less worried after making the efficiency improvements thus far (thanks all!) and after realizing that I probably don't need 1000's of iterations.

-Ryan

p.s. When I do 20-50 iterations, there's usually a couple parameters here and there that weren't even close to converging; these are typically some of the Psi parameters (or Theta), and there's one of those for each site-year-species (a lot). Just thought I'd point out that there are obvious problems with so few iterations. Nonetheless, most parameters are pretty close to true values by that time.

[Michael Betancourt]

One thing I've noticed during my testing is that I don't have to run a Stan model for very many iterations to get close to reasonable results. Honestly, with simulated data and this particular model, I'm finding that ~50 iterations gets be in the ballpark. In fact, I can usually figure out if the model is 'working' if I do ~20. That's orders of magnitude fewer iterations than I would do in JAGS for similar feel of accuracy. 

Yup, HMC is fast.  One thing to be careful about is that while HMC does converge

really fast, our adaptation may not be as fast.  So if your parameters have very

different scales then it’s better to run longer warmup to learn a good metric to

standardize all of those scales.  But if your parameters are all similarly scaled

then you can get away with less warmup.

Is there any standard expectation for the number of iterations I should run the Stan model? If I put in a paper that I ran my model for 200 iterations, a reviewer used to BUGS or JAGS would probably question the choice (if I didn't explain effective sample size, that I'm using Stan, etc). What should I consider to be in the ballpark of 'normal' for number of iterations? this point is germane to the performance of the occupancy model, because "performance" is really about "how long it takes for the model to finish", which obviously depends on # iterations. As Bob briefly alluded to, and is implied by extensive discussions here and elsewhere about "effective sample size", perhaps the easiest way to improve my definition of "performance" is to run the model for fewer iterations.

For convergence you want to quote something like R-hat.

For statistical power you want to quote ESS.

For performance you want to quote ESS / time (or time / ESS).

Raw numbers of iterations is completely meaningless and so

is not helpful in any discussion.  If editors give you crap about

this don’t hesitate to contact us and we can help push back

against poor MCMC understanding.

p.s. When I do 20-50 iterations, there's usually a couple parameters here and there that weren't even close to converging; these are typically some of the Psi parameters (or Theta), and there's one of those for each site-year-species (a lot). Just thought I'd point out that there are obvious problems with so few iterations. Nonetheless, most parameters are pretty close to true values by that time.

Did you full non-center the model?  This could be a big

improvement.

[Ryan Batt]

Thanks Michael, your comments are making me feel good and giving me some confidence in moving forward with my approach of fitting and diagnosing convergence.

As far as I know, I used non-centered methods everywhere. I can repost the model with a message or as a .stan or .txt (for easier browser viewing) for convenience.