fitList error after parLapply estimation of occuMulti models

[Andrew D]

Hello,

I am running into a strange error when trying to coerce a list of named occuMulti( ) model objects to a fitList( ) object (e.g., for use in modSel( ) and other post-processing functions). The code and error is listed below without the input data, but note that 8 total species are involved (i.e. surveyed across 14 sites) and no species-level interactions are estimated. Each model converges on identifiable MLEs, thanks to Nelder-Mead and many iterations.

Any thoughts on the fitList( ) error? Did I find a bug?

Thanks,

Andrew

#########################

library(unmarked)

library(parallel)

# define function to run unmarked on each core

run_um <- function(x){

  unmarked::occuMulti(detformulas=as.character(x$detformulas),

                     stateformulas=as.character(x$stateformulas),

                     data=x$data,

                     maxOrder=1,                       #Species interactions ignored

                     method="Nelder-Mead",   #To avoid "non-finite finite-difference value" errors

                     se=TRUE,

                     engine="C",

                     control=list(maxit=1000000))

}

# Define 17 candidate models

args <- list(

  fm0 = list(detformulas=rep('~1', 8), 

             stateformulas=rep('~1', 8), 

             data=umf2005),

  fmd = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

             stateformulas=rep('~1', 8), 

             data=umf2005),

  fmcp = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

             stateformulas=rep('~1+CP', 8), 

             data=umf2005),

  fafo = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

              stateformulas=rep('~1+CP+FA+FO', 8), 

              data = umf2005),

  fapg = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

              stateformulas=rep('~1+CP+FA+PG', 8), 

              data = umf2005),  

  fafbgs = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

                stateformulas=rep('~1+CP+FA+PF+PG+(FA*PF)+(FA*PG)', 8), 

                data = umf2005),

  f50m_1 = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

              stateformulas=rep('~1+CP+OC50_1', 8), 

              data = umf2005),

  f50m_2 = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

                stateformulas=rep('~1+CP+OC50_2', 8), 

                data = umf2005),

  f100m_1 = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

                stateformulas=rep('~1+CP+OC100_1', 8), 

                data = umf2005),

  f100m_2 = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

                stateformulas=rep('~1+CP+OC100_2', 8), 

                data = umf2005),

  f200m_1 = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

                 stateformulas=rep('~1+CP+OC200_1', 8), 

                 data = umf2005),

  f200m_2 = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

                 stateformulas=rep('~1+CP+OC200_2', 8), 

                 data = umf2005),

  f400m_2 = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

               stateformulas=rep('~1+CP+OC400_2', 8), 

               data = umf2005),

  f800m_1 = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

                 stateformulas=rep('~1+CP+OC800_1', 8), 

                 data = umf2005),

  f800m_2 = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

                 stateformulas=rep('~1+CP+OC800_2', 8), 

                 data = umf2005),

  f1600m_1 = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

                stateformulas=rep('~1+CP+OC1600_1', 8), 

                data = umf2005),

  f3200m_2 = list(detformulas=rep('~1+FS+JD+I(JD^2)+(FS*JD)+(FS*I(JD^2))', 8), 

                stateformulas=rep('~1+CP+OC3200_1+OC3200_2', 8), 

                data = umf2005)

)

# Run models in parallel

start_time <- Sys.time()

      nt <- 17

      cl <- makeCluster(nt, timeout=5184000)

      out2005 <- parLapply(cl=cl, X=args, fun=run_um)

      clusterExport(cl, c("args", "run_um", "umf2005"))

end_time <- Sys.time()

total_time<-end_time-start_time

total_time

stopCluster(cl)

#Coerce output list to fitList 

fitList(fits=out2005)

Error in .local(umf, ...) : 

          255 formulas are required in stateformulas list
##Note: 255 formulas / 17 models = 15, and each candidate model includes 16 total formulas total, 8 detformulas and 8 stateformulas

[Ken Kellner]

Hi Andrew,

Looks like a bug. Can you try this patch?

https://gist.github.com/kenkellner/3262e96fa28084b212523ca24913e47a

If you run this line of code and it should patch the  bug in one step (make sure unmarked is loaded first):

source("https://gist.githubusercontent.com/kenkellner/3262e96fa28084b212523ca24913e47a/raw/fc0585caf4d3b68cb9deaaaa7577375fa3c4dc0c/occuMulti_fitList_patch.R")

You shouldn't need to re-run your models, the fitList command should work again after applying the patch.

Ken

[Andrew D]

Hi Ken,

Thank you for your prompt reply, including code for the patch. I tested it on my machine, and it worked as planned. An "unmarkedFitList" returns, and I can also use it with the modSel( ) function.

I will let you know if I encounter any further issues on this front, but I think that I have what I need. I sincerely appreciate your time and assistance.

All the best,

Andrew

[Andrew D]

Hi Ken,

Unfortunately, I did encounter the same error from other functions following activation of the fitList( ) patch; although, modSel( ) worked as expected and without error. Please see below.

Do each of these (and other post-processing functions) require similar patches for occuMulti( ) models excluding species interactions (i.e. specifying maxOrder=1)?

Thanks,
Andrew

#########################

fmList <- fitList(fits=out2005)

modSel(fmList, nullmod="fm0")

         nPars    AIC delta                 AICwt cumltvWt  Rsq

fm0         16 525.56  0.00 0.9999997012598793589     1.00 0.00

fmcp        64 556.36 30.80 0.0000002053319522057     1.00 0.99

fmd         56 560.31 34.74 0.0000000285767077444     1.00 0.96

f200m_2     72 560.45 34.88 0.0000000266224140482     1.00 1.00

f800m_2     72 560.80 35.23 0.0000000223325079522     1.00 1.00

fafo        80 562.38 36.81 0.0000000101514179519     1.00 1.00

f100m_1     72 564.83 39.27 0.0000000029703107834     1.00 0.99

fapg        80 566.82 41.25 0.0000000011010065162     1.00 1.00

f50m_2      72 568.06 42.50 0.0000000005910835245     1.00 0.99

f400m_2     72 569.06 43.50 0.0000000003589975714     1.00 0.99

f800m_1     72 569.25 43.69 0.0000000003263202653     1.00 0.99

f3200m_2    80 570.32 44.76 0.0000000001907710119     1.00 1.00

f100m_2     72 572.86 47.30 0.0000000000535831695     1.00 0.99

f200m_1     72 572.96 47.40 0.0000000000509908495     1.00 0.99

f1600m_1    72 573.03 47.46 0.0000000000494254893     1.00 0.99

f50m_1      72 573.86 48.29 0.0000000000326314704     1.00 0.99

fafbgs     104 607.77 82.21 0.0000000000000000014     1.00 1.00

 
#Testing a subset of post-processing functions

lapply(residuals(fmList@fits$fm0), function(x) { head(x, 8) })

Error in .local(umf, ...) : 

          255 formulas are required in stateformulas list

ranef(fmList@fits$fm0, species='sp1')

Error in .local(umf, ...) : 

          255 formulas are required in stateformulas list

parboot(fmList@fits$fm0, nsim=5)

Error in .local(umf, ...) : 

          255 formulas are required in stateformulas list

predict(fmList@fits$fm0, type="state")                              #Cannot use backTransform( ) with occuMulti( ) models

Error in .local(umf, ...) : 

          255 formulas are required in stateformulas list

predict(fmList@fits$fm0, type="det")                                  #Cannot use backTransform( ) with occuMulti( ) models

Error in .local(umf, ...) : 

          255 formulas are required in stateformulas list

[Ken Kellner]

Hi Andrew,

These other issues with maxOrder have already been fixed in the dev version of unmarked but not in the CRAN version yet. It can be kind of frustrating to install the dev version so I made another patch that should fix your issues for now:

https://gist.github.com/kenkellner/dab98ebf204234a8b8558a15c1ddce87

source("https://gist.githubusercontent.com/kenkellner/dab98ebf204234a8b8558a15c1ddce87/raw/fd6b172e44c985d2fd883d62e00e3636a05b1152/occuMulti_maxOrder_patch.R")

Let me know if you find other things that don't work. You may also want to just avoid using the maxOrder argument and fill in your interaction formulas with '0' instead if you fit future models. That approach shouldn't have any issues.

Sorry for the trouble,

Ken

[Andrew D]

Understood; thanks again, Ken. I did not know that the dev version had addressed these issues previously, and I will also try using '0' to ignore interactions in future models.

I will let you know if I have any further trouble, but I am confident that your suggestions will help solve any lingering problems.

Best,

Andrew

[Andrew D]

Looks like there's a new maxOrder-related error inhibiting occuMulti( ) estimation with parLapply( ). I should also note that, because of the number of computing resources necessary, I am processing my R scripts on a high-powered computing cluster with a SLURM job submission system (i.e. in case that information matters to the problem at hand). Any thoughts?

Here's the error: 

Error in terms.formula(object, data = data) : invalid power in formula

Calls: unmarkedFrameOccuMulti ... model.matrix -> model.matrix.default -> terms -> terms.formula

Here's a snippet of relevant code:

umf2005 <- unmarkedFrameOccuMulti(y=y2005, siteCovs=site.covs.2005,

                                  obsCovs=obs.covs.2005, maxOrder=0)

# define function to run unmarked in each core

run_um <- function(x){

  set.seed(999)

  unmarked::occuMulti(detformulas=as.character(x$detformulas),

                     stateformulas=as.character(x$stateformulas),

                     data=umf2005,

                     maxOrder=0,

                     method="Nelder-Mead",

                     se=TRUE,

                     engine="C",

                     control=list(maxit=1000))

}

# define one of 29 total candidate models

args <- list(

  fm0 = list(detformulas=rep('~1', 8), 

             stateformulas=rep('~1', 8))
)

# run one of 29 total candidate models in parallel

start_time <- Sys.time()

nt <- 29

cl <- makeCluster(nt, timeout=5184000)

clusterExport(cl, "umf2005")

out_2005 <- parLapply(cl=cl, X=args, fun=run_um)

                      Error in terms.formula(object, data = data) : invalid power in formula

                      Calls: unmarkedFrameOccuMulti ... model.matrix -> model.matrix.default -> terms -> terms.formula

end_time <- Sys.time()

total_time<-end_time-start_time

total_time

stopCluster(cl)

[Ken Kellner]

Hi Andrew,

The minimum value of maxOrder is 1, indicating no species interactions, as you had it originally, but you have 0.

I think maybe the confusion comes from my previous message where I suggested putting '0' in your formulas. If so, to clarify - 

In the simplest case with two species, you will have three occupancy formulas to provide - species 1, species 2, and the interaction of species 1 & 2. So that would look something like:

occ_formulas <- c("~1", "~1", "~1")

If you don't want the interaction, set the appropriate formula(s) to "0" :

occ_formulas <- c("~1", "~1", "0")

This allows you to ignore interactions without setting maxOrder. For three species, ignoring just the three-way interaction would look like this:

occ_formulas <- c("~1", "~1", "~1", "~1", "~1", "~1", "0")

Ignoring both pairwise and three-way interactions:

occ_formulas <- c("~1", "~1", "~1", "0", "0", "0", "0")

Are you eventually planning to include some species interaction terms in your models? Because if not, it would probably be faster to fit 8 single-species occupancy models with occu() rather than a single multispecies model with no interactions.

Ken

[Andrew D]

Hi Ken,

Thank you for clarifying. You're absolutely right that I misunderstood your previous suggestion; my sincerest apologies there. I also appreciate your detailed explanation.

One thing to point out:
My model setup currently looks like the following (i.e. at least for a null model), and I had decided to replicate each function 8 times (i.e. for 8 species) because occuMulti( ) would formerly throw errors concerning "not enough state/detection formulas" when I only included one formula per argument. Perhaps this has something to do with parLapply( ) and its involvement.

# define one of 29 total candidate models

args <- list(

  fm0 = list(detformulas=rep('~1', 8), 

                     stateformulas=rep('~1', 8))
)

How would the above code change, based on your latest suggestions? Would maxOrder also remain at 1?

To answer your question, I was hoping to test some models with interactions in the future, and these initial tests were for evaluating parameter identifiability; however, I worry that low detection estimates may jeopardize models with interactions entirely (e.g., the inputs comprise data on rare grassland birds). Indeed, there may come a time where occu( ) is the only way to go.

Best,

Andrew

[Ken Kellner]

Andrew,

I'm not completely sure what you mean by "replicate each function 8 times", do you mean replicate each *formula* 8 times?

Regardless, this code should work, after applying the patch I sent previously:

run_um <- function(x){

  set.seed(999)

  unmarked::occuMulti(detformulas=as.character(x$detformulas),

                     stateformulas=as.character(x$stateformulas),

                     data=umf2005,

                     maxOrder=1,

                     method="Nelder-Mead",

                     se=TRUE,

                     engine="C",

                     control=list(maxit=1000))

}

# define one of 29 total candidate models

args <- list(

  fm0 = list(detformulas=rep('~1', 8), 

             stateformulas=rep('~1', 8))
)

If you don't want to apply the patch, this is the alternative that doesn't use maxOrder:

run_um <- function(x){

  set.seed(999)

  unmarked::occuMulti(detformulas=as.character(x$detformulas),

                     stateformulas=as.character(x$stateformulas),

                     data=umf2005,

                    #don't set maxOrder

                     method="Nelder-Mead",

                     se=TRUE,

                     engine="C",

                     control=list(maxit=1000))

}

# define one of 29 total candidate models

# You should have 8 detection formulas: 1 per species

# You should have 255 state formulas: (2^S - 1), where S=number of species = 8

# This is 8 first order terms + all the possible 2-way interactions + plus all the possible 3 way interactions and so on up to the 8-way interaction term = 255 total

# In your example with no interactions, the first 8 of these 255 will be "~1", and the  remaining 247 will be "0"

# To check this look at

umf2005@design

# Each column of the matrix above needs 1 formula; length(stateformulas) == ncol(umf2005@design)

args <- list(

  fm0 = list(detformulas=rep('~1', 8), 

             stateformulas=c(rep('~1', 8), rep("0", 255-8))
)

[Andrew D]

Yeah, I meant to say "formula." My apologies again; it's been a frustrating few days.

Thank you for providing those examples as well as for the detailed explanations concerning them. They each give me a much better idea about what is going on behind the scenes for occuMulti( ). Thank you very much for your generous time and assistance, Ken.

[Andrew D]

Hi Ken et al.,

I'm not sure if the following predict(occuMulti( ), "state") errors concern the maxOrder issue discussed above (or perhaps it relates more to this particular thread for occuMS( )), but I thought that I would share them here. Perhaps one or more of the errors relate to lacking co-occurrences (see summary of detection histories below).

Following completion of one of my (candidate) single-season multi-species occupancy models, I get the following results and errors with respect to predicting site-level occupancy rates for one or more species:

summary(fit)

Call: occuMulti(detformulas = as.character(x$detformulas), stateformulas = as.character(x$stateformulas),
    data = umf, maxOrder = 1, method = "Nelder-Mead", se = TRUE, engine = "C", control = list(maxit = 1e+05))

Occupancy (logit-scale):
                   Estimate     SE      z P(>|z|)
[sp1] (Intercept)    1.248  0.655  1.904  0.0569
[sp1] AREA_MN_400   -0.410  0.593 -0.692  0.4890
[sp2] (Intercept)   -2.525  1.614 -1.565  0.1176
[sp2] AREA_MN_400   -7.222  4.190 -1.724  0.0848
[sp3] (Intercept)    2.739  5.745  0.477  0.6336
[sp3] AREA_MN_400    2.115  6.365  0.332  0.7397
[sp4] (Intercept)    1.731 14.143  0.122  0.9026
[sp4] AREA_MN_400   -5.661 21.805 -0.260  0.7952
[sp5] (Intercept)   -0.177  0.811 -0.218  0.8272
[sp5] AREA_MN_400   -1.948  1.668 -1.168  0.2429
[sp6] (Intercept)    6.683  7.957  0.840  0.4009
[sp6] AREA_MN_400   -0.589  4.824 -0.122  0.9029
[sp7] (Intercept)    2.317  2.326  0.996  0.3192
[sp7] AREA_MN_400    2.056  3.783  0.544  0.5868

Detection (logit-scale):
                   Estimate    SE       z  P(>|z|)
[sp1] (Intercept)  1.01256 0.290  3.4926 0.000478
[sp1] FS          -0.35481 0.274 -1.2939 0.195687
[sp1] JD          -0.08730 0.309 -0.2827 0.777387
[sp2] (Intercept)  0.40252 0.618  0.6514 0.514758
[sp2] FS           1.49672 0.933  1.6049 0.108515
[sp2] JD          -0.09781 0.345 -0.2834 0.776872
[sp3] (Intercept) -1.62137 0.502 -3.2302 0.001237
[sp3] FS           0.33226 0.330  1.0081 0.313411
[sp3] JD          -0.00367 0.347 -0.0106 0.991549
[sp4] (Intercept) -1.83220 0.822 -2.2290 0.025811
[sp4] FS          -1.38306 1.084 -1.2758 0.202032
[sp4] JD           0.18443 0.348  0.5300 0.596094
[sp5] (Intercept) -0.22650 0.350 -0.6464 0.517990
[sp5] FS          -0.14937 0.536 -0.2786 0.780572
[sp5] JD           0.35482 0.331  1.0708 0.284245
[sp6] (Intercept)  3.89184 1.020  3.8169 0.000135
[sp6] FS           2.72284 1.148  2.3714 0.017722
[sp6] JD           0.11730 0.424  0.2764 0.782249
[sp7] (Intercept) -1.45085 0.418 -3.4735 0.000514
[sp7] FS           0.07592 0.329  0.2308 0.817445
[sp7] JD           0.33828 0.335  1.0113 0.311879

AIC: 530.5768
Number of sites: 14
optim convergence code: 0
optim iterations: 34622
Bootstrap iterations: 0

#marginal occupancy

predict(fit, 'state', species='sp1') #nsite long
Bootstrapping confidence intervals with 100 samples
 Error in mvrnorm(1, coef(object), Sigma) :
  'Sigma' is not positive definite 

#conditional occupancy & detection

predict(fit, 'state', species=c('sp1','sp2','sp3','sp4','sp5','sp6','sp7'))
Bootstrapping confidence intervals with 100 samples
 Error in mvrnorm(1, coef(object), Sigma) :
  'Sigma' is not positive definite

predict(fit, 'det', species=c('sp1','sp2','sp3','sp4','sp5','sp6','sp7'))
Error in out[[species]] : recursive indexing failed at level 3

#detection histories summarized
AICcmodavg::detHist(fit)
Summary of detection histories:

Note:  Detection histories exceed 80 characters and are not displayed

Species-specific detection histories:

sp1
          000000 000101 001011 100110 101101 110110 111000 111110 111111
Frequency      3      1      1      1      1      1      1      2      3
--------

sp2
          000000 000001 000010 001011 001110 010110 011100
Frequency      8      1      1      1      1      1      1
--------

sp3
          000000 000010 000100 001011 010000 010011 011000
Frequency      6      2      1      1      2      1      1
--------

sp4
          000000 000001 000100 000101 001010 010000 010001 110111
Frequency      7      1      1      1      1      1      1      1
--------

sp5
          000000 000001 000010 000011 000110 001110 110101 111111
Frequency      7      1      1      1      1      1      1      1
--------

sp6
          011010 110111 111101 111110 111111
Frequency      1      1      2      1      9
--------

sp7
          000000 000010 000100 001000 100000 100001 100101 111000
Frequency      5      1      1      2      2      1      1      1
--------

Frequency of co-occurrence among sites:
(sp1(a), sp2(b), sp3(c), sp4(d), sp5(e), sp6(f), sp7(g))
          0000000 000000g 00000f0 00000fg 0000ef0 000d0f0 00c00f0 00c0ef0 0b00000 0b000f0
Frequency       2       2      14       2       3       2       5       1       2       2
          0b0d0f0 0b0d0fg a000000 a0000f0 a0000fg a000ef0 a00d0f0 a00d0fg a00def0 a00defg
Frequency       1       1       1      13       6       8       3       1       2       1
          a0c00f0 a0c00fg a0cd0f0 a0cdef0 ab000f0 abc0ef0 abcdef0
Frequency       1       1       1       1       5       2       1

Proportion of sites with at least one detection:
sp1 sp2 sp3 sp4 sp5 sp6 sp7
0.79 0.43 0.57 0.50 0.50 1.00 0.64

Frequencies of sites with detections:
     sampled detected
sp1      14       11
sp2     14        6
sp3      14        8
sp4      14        7
sp5      14        7
sp6      14       14
sp7      14        9

The current patch I am working with (i.e. source("https://gist.githubusercontent.com/kenkellner/dab98ebf204234a8b8558a15c1ddce87/raw/fd6b172e44c985d2fd883d62e00e3636a05b1152/occuMulti_maxOrder_patch.R")) does not resolve these issues, but perhaps a different patch exists (e.g., as does one for occuMS( ) here). What are some recommendations for working around the aforementioned problems, especially to generate species-level occupancy and co-occupancy/-detection rates and their 95% CIs?

Thanks,

Andrew

[Ken Kellner]

Andrew,

I think the first error related to bootstrapping is a result of issues with optimization. Since you aren't fitting any interaction terms, degree of co-occurrence should not be a major issue. The model apparently converged (after quite a large number of iterations), but occupancy parameter estimates are very imprecise, especially for species 4. This is causing issues with the generation of confidence intervals around predictions. You probably could get predict to work if you set argument se.fit=FALSE, but I would not be very confident in using results from a model with some parameter SEs that are that large. Might be worth trying a different optim algorithm like BFGS or L-BFGS-B, although I suspect you are using Nelder-Mead for a reason.

In your specific case here I'm not sure it makes a lot of sense to fit this as a multi-species model (given my limited info on your situation). You have a pretty small sample size (14 sites) and are trying to handle a relatively large number of species (likely resulting in these poor estimates), but also don't have any species interactions included. I'd at least start by trying to fit single-species occupancy models with covariates for interest for each species, to see if it's even possible in this simpler case. If species 4, for example, is problematic even with a single-species occupancy model then I would consider not attempting to include it in the multispecies model. There may also be better multi-species frameworks to use with your data (see e.g. "Multi‐species occupancy models: review, roadmap, and recommendations", Deverajan et al. 2020).

For the second issue with predicting detection probabilities, I'm not able to replicate that problem myself. If you send me your fitted model saved as an Rdata or Rds I can take a look and see if there's something unique about your model causing issues.

Ken

[Andrew D]

Hi Ken,

My apologies for previously listing "conditional occupancy & detection" in my code annotations above. I meant to list "co-occupancy & co-detection," instead.

Thank you for taking the time to answer my questions today. I also appreciate your insights concerning the aforementioned model. Indeed, I am working with species that (as a group) are rarely detected well, or in sufficient number. Such conditions thus lead to imprecise parameter estimates and related issues with bootstrapping. I did start fitting candidate multispecies models using BFGS or L-BFGS-B, but later encountered numerous "non-finite finite-difference value" errors; each one likely due to imprecise parameter estimates (e.g., boundary violations such as p or psi > 1). Nelder-Mead helped alleviate those errors yet, as you can see above, still led to imprecise parameter estimates for a couple of the rarer species in the group. Additionally, given the sparsity of the data (as well as limited co-occupancy between sites), I was precluded from considering any species-level interactions, unfortunately. Despite there being only 14 sites, each was surveyed 5-6 times per season, for what that's worth.

Thank you also for your suggestions regarding single-species models and for pointing out Deverajan et al. (2020). I will revisit my preliminary single-species models, particularly to check for similar imprecisions, as well as revisit Deverajan et al. (2020).

For the second issue with predicting detection probabilities, I suspect that that problem stems from those issues discussed above, including (1) the lack of species-level interactions in the model, (2) limited species co-occurrences between sites, and (3) imprecise model parameter estimates due to reduced sample size and other related constraints. I was able to identify less-aberrant models in my candidate set (i.e. those still outcompeting the null [intercept-only] model), for which I can generate more precise parameter estimates as well as predictions of occurrence and detection (e.g., marginal, conditional, and co-occupancy and -detection).

Best,

Andrew