help with an n mixture model

[Brian Hatfield]

Hello all,

I have hit a wall here, might be a simple error, but i can't seem to find it, maybe someone can help. 

I am attempting to run an n-mixture on model on camera trap data.  the data was collected at a different study area each year, so I want to loop over years, then sites/reps for each year. 

Another small wrinkle is that I am using a long form (vector) for the detection data, and then supplying indexing variables as constants. I have removed NAs from the data set so the capture histories have uneven numbers of replicates by site (dealt with by the indexing variables).

the model includes some covariates  on detection (days since lure applied, Julian date, julian date sq), I get the same error below if I remove the covs from the model.

I have fairly meticulously checked the indexing vars, this is the messiest  part of the model, but seems to be correct as far as i can tell. 

I can create and configure the model, but when running the samples I get the following error:

running chain 1... Error in if (nrow > rows) { : missing value where TRUE/FALSE needed In addition: Warning message: In container$resize(as.integer(k)) : NAs introduced by coercion to integer range

Thank you in advance, I would be grateful for any suggestions

I can post the data if it would be helpful

Model:

fm.site.Nmix <- nimbleCode({
 
  # Priors
  alpha0 ~ dnorm(0, 0.1)

  # alphas - det covs
  for(l in 1: n.alpha){               
    alpha[l] ~ dnorm(0, 0.1)  
  }
 
  for(k in 1:n_years){               ## loop over years
  
    lambda[k] ~ dgamma(0.001, 0.001)
   
    for (i in 1:n.sites.by.year.index[k]) {              # Loop over sites (for year K )
      N[site.by.year.index[k,i]] ~ dpois(lambda[k])
     
      for (j in 1:n.reps.by.year.site.index[k,i]) {        # Loop over replicates (by site and year)
       
        y[rep.by.year.site.index[k,i,j]] ~ dbin(p[rep.by.year.site.index[k,i,j]] , N[site.by.year.index[k,i]])
       
        logit(p[rep.by.year.site.index[k,i,j]]) <- alpha0 +
          (alpha[1:n.alpha] %*% det.covs.on.reps[rep.by.year.site.index[k,i,j], 1:n.alpha])[1,1]
 
      }
    }
  }
})

indexing vars:

n.sites.by.year.index[k]:  vector of number of sites by year

site.by.year.index[k,i]:  ragged matrix of individual sites by year

n.reps.by.year.site.index[k,i]: ragged matrix of number of reps for each site by year

rep.by.year.site.index[k,i,j]:  3d array of individual reps by site and year

Data lists:

dat<-list(y = y,
          det.covs.on.reps = rep.cov)

const<-list(
            n_years = 5,
            n.sites.by.year.index = index$yr.nsite,
            rep.by.year.site.index = index$yr.site.rep,
            n.reps.by.year.site.index = index$yr.nreps,
            site.by.year.index = index$yr.site,
            n.alpha = ncol(rep.cov)
            )

init<- list( N = n.start,   #total detections by site plus 1 (from ahm -> nimble on git)
             lambda = rep(.5, 5),
             alpha = rep(0,ncol(rep.cov)),
             p = runif(length(y), 0, 1)
            
)

params <- c('N', 'alpha0', 'alpha', 'lambda', 'p')

model call:

run.model<-nimbleModel(
  code = fm.site.Nmix ,
  constants = const,
  data = dat,
  inits = init 
)

model.config<-configureMCMC(run.model, monitors = params, onlySlice = F)

mcmc.model<-buildMCMC(model.config, enableWAIC = T)

samples<-runMCMC(
                 mcmc.model,
                 niter = 100,
                 nburnin = 0,
                 thin = 0,
                 summary = T,
                 WAIC = T)

[Ben Goldstein]

Hi Brian,

If you're able to provide the data to me (if you want to reply off-list) I'd be happy to take a look at debugging.

Incidentally you might be interested in the dNmixture() custom distribution provided in the nimbleEcology package. I can whip up an example of that with the data you provide, too.

Ben

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/71998b5a-faf6-49b4-b8ec-5293d017bcd8n%40googlegroups.com.

[Brian Hatfield]

Hi Ben,

Here are rds files of the model input lists. Will that work?

also included is the full long capture history with covariates (weekly.cap.hist.wtj.rds).  fyi i have the data set up with weekly replicates.

I am interested in the nimbleEcology nmixture model (and occupancy.. etc),  was actually planning on trying to run and compare the nimble ecology model to this model if i could get it running.  Can I index the dNmixture distribution with a long form capture history similar to the model above? 

Thank you for your help,

Brian

[Ben Goldstein]

Bringing this back to the list after some off-list conversation.

The immediate issue turned out to be the use of thin=0 instead of thin=1 in the call to runMCMC. This took me debugging to notice when I didn't immediately spot it--maybe a more informative error message could be useful here.

Brian, I also prepared a version of your model code that uses dNmixture from nimbleEcology (below). The main advantage of this is that the distribution marginalizes over values of N, which is no longer a sampled node in the model, which may improve MCMC runtime. The main downside is that if you were planning to directly monitor N, this is no longer possible.

To marginalize over N, dNmixture takes a vector of observations from a single site and (if they're visit-variant) a vector of probabilities. This means you have to be a little thoughtful with your indexing. I got a version working with vectors you already provided in constants, but it's not pretty. My go-to is to provide a vector (or in your case maybe a matrix) of start/end indices, i.e. y[start[k,i]:end[k,i]] ~ dNmixture(...).

Hope this is helpful.

library(nimbleEcology)

fm.site.Nmix <- nimbleCode({
 
  # Priors
  alpha0 ~ dnorm(0, 0.1)
 
  # alphas - det covs
  for(l in 1:n.alpha){              
    alpha[l] ~ dnorm(0, 0.1)  
  }
 
  for(k in 1:n_years){               ## loop over years

    lambda[k] ~ dgamma(0.001, 0.001)

    for (i in 1:n.sites.by.year.index[k]) {              # Loop over sites (for year K )

      y[(rep.by.year.site.index[k,i,1]):(rep.by.year.site.index[k,i,n.reps.by.year.site.index[k,i]])] ~
        dNmixture_v(
          lambda = lambda[k],
          prob = p[(rep.by.year.site.index[k,i,1]):(rep.by.year.site.index[k,i,n.reps.by.year.site.index[k,i]])],
          Nmin = -1, Nmax = -1,
          len = n.reps.by.year.site.index[k,i]

        )
     
      for (j in 1:n.reps.by.year.site.index[k,i]) {        # Loop over replicates (by site and year)

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/f8e681a2-059f-4c9c-950b-dda95f667a6cn%40googlegroups.com.

[Brian Hatfield]

Thank you Ben!   I'm very happy to have my model working.

Ill try the dNmixture version, and the start:end indexing (might be better/ less messy in my original model as well).

for a covariate on lambda would you include something like the following at the bottom of the i loop?

lambda[index[k,i]] <- beta0 + exp(betas[ :  ] %*% beta.covs[ index[k,i]  ,  :  ])     

Brian