kenkellner / ubms Goto Github PK

Models with random effects take a while to run, it may need to be precomputed. Potentially useful:

https://ropensci.org/technotes/2019/12/08/precompute-vignettes/
https://github.com/ropensci/eia/blob/master/vignettes/

Hi Ken,
when using a very long model formula, check_formula throws this warning:

In if (grepl(":|/", char)) { :
  the condition has length > 1 and only the first element will be used

It is because

deparse1() instead of deparse() would avoid that, but that's only for R4.0 and higher.

char <- paste(deparse(formula_long, width.cutoff = 500L))

seems to work though.

Altenratively, if char is a vector of length > 1 (this might print unwanted output):

sapply(char, FUN = function(x){
  if(grepl(":|/", x)){
    stop("Nested random effects (using / and :) are not supported",
         call.=FALSE)
  }
  })

Here's a reproducible example:

formula_long <- formula("~ (longcovariate1 + longcovariate + longcovariate3 + longcovariate4 || group)")

char <- paste(deparse(formula_long))

if(grepl(":|/", char)){
  stop("Nested random effects (using / and :) are not supported",
       call.=FALSE)
}

Note: I didn't check how a long formula would affect get_reTrms a few lines below

Like unmarked::gdistremoval (probably without the temporary emigration part though)

Hi Ken,

I am trying to run ubms model to my data but because I did not obtain similar estimates to the models with unmarked, I tested both with example data from unmarked.

The results with ubms sometimes do not have a good fit and sometimes the R-hat is not appropriate.

data(crossbill)

site_covs <- crossbill[,c("id", "ele", "forest")]

y <- crossbill[,c("det991","det992","det993")]

date <- crossbill[,c("date991","date992","date993")]

umf <- unmarkedFrameOccu(y=y, siteCovs=site_covs, obsCovs=list(date=date))

stan_global <- stan_occu(~scale(date)~scale(forest)+scale(ele), data=umf, chains=4)

um_global <- occu(~scale(date)~scale(forest)+scale(ele), data=umf)

cbind(unmarked=coef(um_global), stan=coef(stan_global))

psi(Int) -0.7434025 1.4468469
psi(scale(forest)) 0.9782374 2.4063905
psi(scale(ele)) 0.5898850 1.4289875
p(Int) -0.6738667 -0.9996905
p(scale(date)) 0.5505878 0.5685348

Any suggestion on how to reach similar results?

Based on AICcmodavg

Hi!

Your package looks awesome! I'm excited to try it but couldn't find any code examples on how to do multi-season occupancy (non-stacked). Is this possible with ubms and if so do you have a code example I could try? If not, do you have plans to add this feature soon?

Thank you!

-Carrie

Started in this branch, but it doesn't recover parameters.

Hi Ken,

I'm not sure whether the use of a bayes R2 would be broadly suitable for occupancy models given the limitations with its application in logistic regression and often low values. However I had a crack at getting a bayes R2 for the state, detection submodels (and their combined value) using the residual method (https://avehtari.github.io/bayes_R2/bayes_R2.html). A bit of a messy function and example below.

# R2 example for logistic occupancy model
    library(ubms)
#> Loading required package: unmarked
#> Loading required package: lattice
#> 
#> Attaching package: 'ubms'
#> The following objects are masked from 'package:unmarked':
#> 
#>     fitList, projected
#> The following object is masked from 'package:base':
#> 
#>     gamma
    
    # R2 function
    bayes_R2_res_ubms <- function(fit, re.form = NULL, draws = draws) {
        
        # Get the observed occupancy at each site for each observation period
        y <- ubms::getY(fit)
        
        # Get the observed occupancy at each site
        y_mod <- matrix(apply(y, 1, FUN = function(x) max(x, na.rm = TRUE), simplify = TRUE), ncol = 1)
        
        # Get the linear predictor for the 'state'
        ypred <- ubms::posterior_linpred(fit, transform = TRUE, submodel = "state", re.form = re.form, draws = draws)
        yp <- ubms::posterior_linpred(fit, transform = TRUE, submodel = "det", re.form = re.form, draws = draws)
        
        # Assume perfect detection (no effect of detection submodel)
        yp1 <- yp
        yp1[!is.na(yp1)] <- 1
        
        # For each draw obtain the probability for each site that a detection will occur (0-1).
        # Probably a less messy way to do this
        ys <- list("state" = yp,
                   "detection" = yp1)
        
        ypred_mod <- list()
        r2 <- list()
        
        for(j in 1:length(ys)) {
            
            # detection * state
            ypred_prob <- ys[[j]]
            
            # draws x sites x obs
            ypred_cond <- array(data = NA, dim = c(draws, dim(ypred)[2], fit@response@max_obs))
            ypred_mod[[j]] <- ypred
            
            # loop over draws
            for(i in 1:draws) {
                # repeat the latent state by observation periods
                ypred_vec <- rep(ypred[i,], each = fit@response@max_obs)
                # det * state
                ypred_prob[i,] <- ypred_prob[i,]*ypred_vec
                # force into matrix with nsites X nobs
                ypred_cond[i,,] <- matrix(ypred_prob[i,], ncol = fit@response@max_obs, byrow = TRUE)
                # 1 minus the probability that all observations are zero = at least 1 detection
                ypred_mod[[j]][i,] <- 1-apply(1-ypred_cond[i,,], 1, function(x) {
                    prod(x, na.rm = TRUE)
                })
                
            }
            
            if (fit@response@z_dist == "binomial" && NCOL(y_mod) == 2) {
                trials <- rowSums(y_mod)
                y_mod <- y_mod[, 1]
                ypred_mod[[j]] <- ypred_mod[[j]] %*% diag(trials)
            }
            e <- -1 * sweep(ypred_mod[[j]], 2, y_mod)
            var_ypred <- apply(ypred_mod[[j]], 1, var)
            var_e <- apply(e, 1, var)
            r2[[j]] <- var_ypred / (var_ypred + var_e)
        }
        r2[[3]] <- r2[[1]] - r2[[2]]
        names(r2) <- c("both", "state", "det")
        return(r2)
    }
    
    # Data simulation
    
    set.seed(123)
    dat_occ <- data.frame(x1=rnorm(500))
    dat_p <- data.frame(x2=rnorm(500*5))
    y <- matrix(NA, 500, 5)
    z <- rep(NA, 500)
    b <- c(0.4, -0.5, 0.3, 0.5)
    re_fac <- factor(sample(letters[1:26], 500, replace=T))
    dat_occ$group <- re_fac
    re <- rnorm(26, 0, 1.2)
    re_idx <- as.numeric(re_fac)
    idx <- 1
    for (i in 1:500){
        z[i] <- rbinom(1,1, plogis(b[1] + b[2]*dat_occ$x1[i] + re[re_idx[i]]))
        for (j in 1:5){
            y[i,j] <- z[i]*rbinom(1,1,
                                  plogis(b[3] + b[4]*dat_p$x2[idx]))
            idx <- idx + 1
        }
    }
    
    # unmarked frame
    umf <- unmarkedFrameOccu(y=y, siteCovs=dat_occ, obsCovs=dat_p)
    
    # model
    options(mc.cores=3) #number of parallel cores to use
    fm <- stan_occu(~x2 ~x1 + (1|group), umf, chains=3)
    
    bayes_R2_res_ubms(fm, draws = 100)
#> $both
#>   [1] 0.2436708 0.1813600 0.1861829 0.1968885 0.2093023 0.2065249 0.2262449
#>   [8] 0.1605566 0.2356042 0.2089506 0.2313770 0.2569806 0.2476316 0.2080628
#>  [15] 0.1884270 0.2458283 0.2353730 0.2410739 0.2087004 0.2352926 0.2182499
#>  [22] 0.2195304 0.2197247 0.2299345 0.1642092 0.2512100 0.2537433 0.1866575
#>  [29] 0.2420921 0.1616609 0.1850477 0.2040931 0.2660747 0.1736265 0.2113426
#>  [36] 0.2345318 0.2077254 0.1940378 0.2166815 0.2828791 0.2590453 0.1380645
#>  [43] 0.2494774 0.2094812 0.1596296 0.2125061 0.2406926 0.2306142 0.2238410
#>  [50] 0.1424783 0.2422054 0.2399446 0.2433104 0.2083924 0.2163885 0.2143759
#>  [57] 0.2868863 0.2545129 0.2662583 0.2061908 0.2272771 0.1880253 0.2036795
#>  [64] 0.2259326 0.2396712 0.1560981 0.2029403 0.2575136 0.2447248 0.2600542
#>  [71] 0.1503401 0.1612022 0.2184497 0.2333725 0.2235703 0.2745752 0.1902804
#>  [78] 0.2021373 0.1704586 0.2658361 0.2470988 0.2298406 0.2281244 0.2116103
#>  [85] 0.2407505 0.2564699 0.2642025 0.2530184 0.2964650 0.1884810 0.1804507
#>  [92] 0.2371928 0.2304418 0.2121007 0.2515888 0.2200955 0.1956262 0.2275631
#>  [99] 0.2210517 0.2139585
#> 
#> $state
#>   [1] 0.17605019 0.10263418 0.10580724 0.13545364 0.14482887 0.13800919
#>   [7] 0.17355460 0.08628623 0.17066872 0.14583832 0.16249959 0.20375155
#>  [13] 0.18365678 0.13735974 0.11535767 0.17160635 0.16365889 0.18999839
#>  [19] 0.16342425 0.17257228 0.16156030 0.15115874 0.15180402 0.17367579
#>  [25] 0.09177282 0.18924085 0.20799994 0.13745090 0.18211136 0.09588723
#>  [31] 0.11718993 0.14892238 0.22422619 0.13182738 0.15381415 0.17530606
#>  [37] 0.14920495 0.12388629 0.15416240 0.24336614 0.20041741 0.08053604
#>  [43] 0.20827089 0.14264055 0.09731327 0.14357953 0.18287715 0.16933439
#>  [49] 0.15418602 0.07828137 0.17981013 0.17819720 0.19682582 0.15552595
#>  [55] 0.15226781 0.15547727 0.24004102 0.20505126 0.22142930 0.14618668
#>  [61] 0.15626259 0.12546194 0.14349652 0.16704033 0.18849389 0.08944449
#>  [67] 0.12046083 0.22172548 0.19268544 0.20941069 0.07355600 0.11754425
#>  [73] 0.15150645 0.17893636 0.16290501 0.22485436 0.12711429 0.15046677
#>  [79] 0.11693493 0.22099586 0.18770485 0.17481347 0.16442692 0.14398650
#>  [85] 0.17886823 0.20548550 0.20694618 0.20829883 0.24947301 0.11531789
#>  [91] 0.10101353 0.18807278 0.18736563 0.14776186 0.19111258 0.16002263
#>  [97] 0.11816607 0.17083967 0.16663188 0.15840369
#> 
#> $det
#>   [1] 0.06762062 0.07872586 0.08037562 0.06143486 0.06447340 0.06851572
#>   [7] 0.05269027 0.07427036 0.06493549 0.06311229 0.06887738 0.05322905
#>  [13] 0.06397480 0.07070302 0.07306933 0.07422198 0.07171414 0.05107554
#>  [19] 0.04527614 0.06272035 0.05668959 0.06837161 0.06792069 0.05625875
#>  [25] 0.07243638 0.06196911 0.04574332 0.04920663 0.05998069 0.06577362
#>  [31] 0.06785781 0.05517070 0.04184849 0.04179915 0.05752846 0.05922575
#>  [37] 0.05852047 0.07015146 0.06251909 0.03951300 0.05862791 0.05752848
#>  [43] 0.04120654 0.06684066 0.06231634 0.06892656 0.05781542 0.06127982
#>  [49] 0.06965496 0.06419698 0.06239527 0.06174742 0.04648460 0.05286642
#>  [55] 0.06412072 0.05889858 0.04684531 0.04946166 0.04482896 0.06000413
#>  [61] 0.07101450 0.06256336 0.06018296 0.05889223 0.05117732 0.06665366
#>  [67] 0.08247942 0.03578808 0.05203937 0.05064353 0.07678410 0.04365799
#>  [73] 0.06694323 0.05443615 0.06066532 0.04972085 0.06316610 0.05167058
#>  [79] 0.05352365 0.04484027 0.05939398 0.05502711 0.06369745 0.06762376
#>  [85] 0.06188222 0.05098439 0.05725632 0.04471958 0.04699198 0.07316309
#>  [91] 0.07943715 0.04912006 0.04307613 0.06433879 0.06047627 0.06007289
#>  [97] 0.07746012 0.05672347 0.05441978 0.05555479

^{Created on 2021-11-11 by the reprex package (v2.0.1)}

• predict
• posterior_predict, simy, simz, etc.
• gof

Hi Ken,
would it please be possible to allow a named list of models as input for fitList, as in unmarked::fitList() via its "fits" argument?
I believe it would allow for much easier (and automated) model comparison if one wouldn't have to manually add and name each model in the fitList.

# ... data preparation from github readme file

(fm <- stan_occu(~x2 ~x1 + (1|group), umf, chains=3, cores=3))

The following works as expected:

 fl1 <- fitList(mod1 = fm, 
                mod2 = fm)
  modSel(fl1)

Providing a list as input fails (and seems to use unmarked::fitList, despite asking for ubms::fitList):

modList <- list(modelname1 = fm, 
                modelname2 = fm)
fl2 <- ubms::fitList(modList)

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘getData’ for signature ‘"ubmsFitOccu"’
In addition: Warning message:
In unmarked::fitList(...) :
  If supplying a list of fits, use fits = 'mylist'

I am trying to install ubms but am struggling with installing rstan. Any suggestions?

• installing source package 'ubms' ...
  ** using staged installation
  Error: .onLoad failed in loadNamespace() for 'rstan', details:
  call: get_rng(0)
  error: function 'Rcpp_precious_remove' not provided by package 'Rcpp'
  Execution halted
  ERROR: configuration failed for package 'ubms'
• removing 'C:/Users/serge/Documents/R/win-library/4.0/ubms'

Hi Ken,
I am interested in exploring the spatial patterns of residuals of some occupancy models. However, I am somewhat lost when it comes to understanding the matrix created by the residuals() function of an ubm object.
I fit a model with covariates to a matrix with 30 sites and 7 survey occasions, all well so far.
Once I try to extract the residuals with
res <- residuals(fit_ubms, submodel= "state")
I get a res matrix of [1:3000, 1:9]. This is where I'm confused, I don't know if the rows represent my 30 sites100 times, and I figure the columns will be the draws?
When I extract the residuals for the detection submodel I now get a matrix of [1:3000, 1:63]. Am I correct in assuming that is my 30sites X 100 and my 7 survey occasionsX 9?

If so, to proceed to construct a Moran I correlogram plot I should use the mean of all repetitions per site, right?

An apology for how silly the question may be and thank you very much for your time.

Hi Ken,
first of all, this is a fabulous package and overall works very well already. Really exciting. Thank you very much for your efforts.

Here's a little problem I encountered. I ran a big model and wished to save disk space by doing:

stan_occu(..., save_warmup = FALSE)

When I tried to use predict(), SD was all NA. I get values for SD though when I don't set save_warmup = FALSE. It turned out that nsamples() (called by predict) returns 0 because:

sum(object@stanfit@sim$n_save - object@stanfit@sim$warmup)
with n_save = 2500 2500 2500
and
warmup = 2500

I didn't follow the issue further (how exactly it impacts sim_lp etc) but it seems to be the cause of SD being NA in the prediction.
Best regards,
Jürgen

An issue was exposed when I was trying to upload a new major version of StanHeaders to CRAN, but actually the issue is present (albeit not being caught) with the ubms / StanHeaders on CRAN currently,

If you add lines to the effect of

  if (nb != cols(pars)) {
    print("pars = ");
    for (n in 1:rows(pars)) print(pars[n, ]);
    print ("dist = ", dist);
    print("beta = ");
    print(beta);
    print("idx = ", idx);
    print("nb = ", nb);
}

in the lp_priors function in ubms/inst/stan/include/functions_priors.stan , your examples and vignettes will run but output things like

pars = 
[0,1]
[1,1]
[0,0]
dist = [4,0,5]
beta = 
[-4.00449]
idx = 2
nb = 1

In other words, idx:nb evaluates to 2:1, which is not a valid index in the Stan language and is not correct logic for subsetting beta and pars anyways. The new StanHeaders catches this, but unfortunately this issue is one of the few things preventing the StanHeaders from making its way to CRAN.

If you could fix this as soon as possible, that would be great.

Hi Ken,

Thought I might raise a suggested enhancement for an expansion to the list of prior distributions available in ubms. I'm interested in using some of the priors listed here (e.g. hs() and laplace())

Thanks
Justin

Sometimes returns a list instead of a vector

At least make it possible to change SD of normal priors.

Hey Ken,

When I fit a model with the following nested random effect:
(1|vegetation_type/site)

I get

Error in .subset2(x, i, exact = exact) : no such index at level 1

Is this just because nested random effects or not supported? Or am I doing something wrong?

Cheers,
Matt.

Hi Ken,
saving a large stan_occu() model (>100k sites, 6 occasions, 2500 iterations without warmup) results in an .RData file of almost 10Gb size, a bit unwieldy on an average computer and also very slow to process (even printing the model summary can take several minutes). It is caused by the posterior samples for log_lik in each chain within model@stanfit@sim$samples.
There are around 100 actual model parameters in my example, but almost 90,000 log_lik parameters.

Returning log_lik is currently hard-coded and one cannot pass the pars argument to stan() via ... since stan_occu() defines it already. Removing log_lik from the posterior samples manually breaks print() and possibly other methods, so that's not an option either.

So, could computing / returning log_lik be made optional please to avoid such huge model files? From what I understand Stan does not automatically compute and store the log-likelihood, so it might be nice if ubms would also give that option (returning log_lik should be on by default though in my opinion).

Thank you!

See CRAN checks.

The more top-level .stan models exist, the longer compilation takes (currently about 5 min 42 seconds on my laptop). I should be able to combine some models and use an if/else structure to minimize the number of models. Possible combinations of existing models:

1. Basic single season models (occu, occuRN, pcount)
2. Dynamic occupancy (colext)

Hi Ken,
I have very large and complex model for which gof() takes over 5 minutes per posterior sample. Below I suggest a way to speed up gof considerably for this large models.

The model contains >50000 sites and >50 occasions. There are about 16000 unique cohorts for the gof test. More than half the cohorts had only one site, and mean cohort size was 3.5.

I found that in my example, this line in mb_chisq():

took about 10ms for one iteration, there were 16000 iterations, and mb_chisq is called twice during each posterior sample. So that really adds up.

If the line is replaced with:

 p_sub <- p[Reduce(c,  lapply(inds[[i]], FUN = function(k){
      seq((((k - 1) * object@response@max_obs) + 1) , (k * object@response@max_obs))
    }    ))]  

it only took about 40 microseconds instead of 10 milliseconds in my example (and it makes the object site_idx obsolete). Of course the speed gains don't fully translate to gof() because of the other functions (get_exp_counts / get_obs_counts: 110/130 microseconds), but overall gof() was still much faster (10 posterior samples: 4 minutes vs. >1 hour).

Please also note that it is slower than the current code in the example in your readme (20s instead of 12s for 1000 draws). When I added some NAs to the detection matrix (resulting in 8 cohorts) it was 24s instead of 17s.

So it's not a universal solution, but may help in large models (in my case, site_idx had over 3 million entries, so naturally "%in%" got a bit slow).

Also, before that I also wrote another little hack of sim_gof() which allows parallel computation of gof across the posterior samples. But the speed gains were less noticable than the modification above and the parallelization overhead was immense for the huge model (exporting a >1Gb model to multiple cores). If it is of interest I can also share it though.

• stan_occu
• stan_occuRN
• stan_pcount
• stan_multinomPois
• stan_distsamp
• stan_colext

Hello,

I am trying to use the ubms package to fit an occupancy model with random effects. for model diagnosis, the LOO function does not seem to work anymore. I keep getting this error:

' Error in UseMethod("loo") :
no applicable method for 'loo' applied to an object of class "c('ubmsFitOccu', 'ubmsFit')" '

This even happens with the example code on the ubms github page. Has this issue been encountered before?

See Gelman et al. 2021, Regression and Other Stories, pg. 179

Best regards Ken,
First I wanted to congratulate and thank you for your work on unmarked as well as this new package. I am working with basicoccupancy models for camera trap data. I have 67 sites with 27 sampling events, 110 detections, 1038 no-detections and 661 NAs. Everything works perfect when fitting the models with ubms and I have very similar results to unmarked.
My best model was this

(mg5_stan <- stan_occu(~ Slope~ 1, data=gf, chains=3, iter=600)) # Only 600 iterations for the purposes of the exercise.

I have no problem generating the residuals plots for the detection parameter

plot_residuals(mg5_stan, submodel="det")

plot_residuals(mg5_stan, submodel="det", covariate="Slope")

But the problem comes when I want to make the residuals diagnostic plots for the occupancy parameter

plot_residuals(mg5_stan, submodel="state")
Error: Couldn't find working breakpoints

I am new to Bayesian approach and the Wright et al 2019 definition of residuals. which doesn't give me much clue as to what might be generating this error. is it a bug or am I missing something?

I really like that this package is bundling so many models in such an accessible way! Thank for your for efforts in that direction.

Eventually, users will want to build upon the included models (I do!). The packaged stan code seems to be written to be able to accommodate all the models at once, making it difficult to extract the code specific to a given model, understand said code, and build upon it.

I understand that modifying the code so that only code relevant to the chosen model is included at runtime would require a lot of efforts. As a first step, would it be possible to comment the stan code so that it is easier to understand? It would add great value to those who like me, have found using ubms a much needed shortcut from coding models directly from the literature.

Requires adding some auxiliary data slots into ubmsResponse objects

When I read ubms paper in Methods in Ecology and Evolution, it was clear for that it is possible to estimate random slopes, not only random intercepts, even though I just found examples with random intercepts.

When I try the model with only random intercet, it runs perfectly:
stan_occu(~method+hort+wash+(1|parasite)+(1|id)~hort+type+(1|parasite)+(1|id),input,chains=4,iter=4000,cores=4,control=list(adapt_delta=.90))

But when I include any random slope (I want to see effect to very among parasites), it fails:
stan_occu(~method+(1+method|parasite)+(1|id)~hort+(1+hort|parasite)+(1|id),input,chains=4,iter=4000,cores=4,control=list(adapt_delta=.90))

Error: Failed to create random effects model matrix.
Possible reasons:
(1) Correlated slopes and intercepts are not supported. Replace | with || in your formula(s).
(2) You have specified random slopes for an R factor variable. Try converting the variable to a series of indicator variables instead.

Could you gently guide me. I appreciate any advice. G.

See fix here: https://github.com/rbchan/unmarked/blob/e128480ed012131a46ffa50587064742659c3946/R/mixedModelTools.R#L111

Tests are taking too long on Travis. Run models with fewer sites/iterations for majority of tests, run longer tests locally only.

Hi Ken,

Sorry for posting this here, but I am just wondering how you fixed the Clang-UBSAN error in your 1.0.2 package release. I encountered a similar error in my R package (which depends on RStan) and I just couldn't find a way to fix it.

Thank you so much!
James

Single-season only for now

See draft here: https://kenkellner.com/blog/ubms-vignette.html

Hello Ken,

Many thanks for the hard work on this package, and just wanted to report on a potential issue. After running the goodness-of-fit test function, I've been unable to create the posterior predictive check plot. I get the following error when running the vignette code:

fm_fit <- gof(fm, draws=500)
plot(fm_fit)

Error in as.double(y) : 
  cannot coerce type 'S4' to vector of type 'double'

I also see this error when running my own occupancy models. Here's my session info if that's helpful:

> sessionInfo()
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ubms_1.2.6     unmarked_1.3.2

Thanks!

Jay

Hi Ken,

I'm predicting estimates from a few models to make some custom plots, but having trouble removing the effect of another explanatory variable - it's categorical so I can't just take the mean etc.

Just wanting to confirm there isn't currently a simple way to exclude particular effects (e.g. an argument in the predict function like there is for random effects or by specifying a NA - if not can I put in a suggestion for an enhancement 😉), and ask if you have any recommendations of a work-around?

Cheers,
Matt.

Need to change how observed data are represented here:

So that output says e.g. "running model occu" rather than "running model single_season".

• stan_occu
• stan_occuRN
• stan_pcount
• stan_multinomPois
• stan_distsamp
• stan_colext

ranef should return just random part, and coef should return the full random slope/intercept.
Also I should add a fixef function.

When the true value of a coefficient is relatively large (>3-4 for a standardized covariate), the estimated value might be pulled too much towards zero. Need to set even weaker default priors e.g. normal(0,10).