Issues setting priors for elemets in an array

[Jordan Dinardo]

Hi all,

I am VERY new to nimble and currently trying to convert my multivariate state-space model in jags to nimble code and I'm running into what is probably rudimentary issues. 

I am trying to set the priors on the first timestep of each population's (location) state vector within an array. 

for ( l in 1:Nloc){

X[1,1,l] ~dnorm(1,5)

}

where Nloc equals the number of states, which in this case represents populations from different locations. X will ultimately be an array of matrices, where each matrix has the dimensions of Nyears by 1, representing the population state vector through time.   

The code above outputs the error: 

Error in array(as.numeric(NA), dim = varInfo[[iV]]$maxs) : 

  negative length vectors are not allowed

In addition: Warning message:

In array(as.numeric(NA), dim = varInfo[[iV]]$maxs) :

  NAs introduced by coercion

I have tried to set the dimensions of X in the nimbleModel with no luck. Might anyone know what could be the issue and any solutions? I've come across nimbleFunctions as a solution, but would appreciate some direction.

Thanks,

Jordan 

[Perry de Valpine]

Hi Jordan,

Glad to hear of your application.

Can you please send code that reproduces the error?

If I run it as below, it works fine.

-Perry

nc <- nimbleCode({

for ( l in 1:Nloc){
 X[1,1,l] ~dnorm(1,5)
 }

})

m <- nimbleModel(nc, constants = list(Nloc = 3))

--
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/c8de1530-5794-4d9b-ab3e-ed49ee220973n%40googlegroups.com.

[Jordan Dinardo]

Hi Perry,

Thanks for getting back to me so quickly. Here is the nimble code. After some tinkering I realize that the error code most likely stems from the observation model later in the model (highlighted in yellow). Here is the full model:

nc <- nimbleCode({ 

  #---------------

  # DEFINE PRIORS

  #---------------

  

  ## Prior for step parameter (U parameter in MARSS models) 

  #--------------------------------------------------------

  

  for (l in 1:Nloc){

    

    step[l] ~ dunif(-10,10) # time-invariant

    

  }

  

  ## Prior for Q parameter (Process error MARSS models)

  #--------------------------------------------------------

  

  # Q prior from wishart distribution: dwish ~ (omega, degrees of freedom)

  

  tauQ[1:Nloc, 1:Nloc] ~ dwish(Omega[1:Nloc,1:Nloc], Nloc)

  

  for (l in 1:Nloc) {

    

    Omega[l,l] <- l

    

  }

  

  for (i in 1:(Nloc-1)) {

    for (j in (i+1):Nloc) {

      

      Omega[i,j] <- 0

      Omega[j,i] <- 0

      

    }

  }

  

  # MVN Process error

  

  sigmaQ[1:Nloc, 1:Nloc] <- inverse(tauQ[1:Nloc, 1:Nloc]) # convert from precision to sd

  

  ## Prior for R parameter (observation error in MARRS model)

  #----------------------------------------------------------

  

  # Assuming observation error is different for each location

  

  for(l in 1:Nloc) {

    

    tauR[l] ~ dgamma(0.001,0.001)

    sigmaR[l] <- 1/sqrt(tauR[l]) 

    

  }

  

  ## MVN prior on initial state for each location

  #---------------------------------------------

  

  

  

  for (l in 1:Nloc){

    

    X[1,1,l] ~dnorm(1,5)

    

  }

  

  ## Priors for thresholds

  #------------------------  

  

  for(l in 1:Nloc){

    

    tau[l,1] ~ dnorm(0,0.01)

    

    for(j in 1:3){

      

      delta[l,j] ~ dexp(2)

      tau[l,j+1] <- tau[l,j] + delta[l,j]

      

    }

  }

  

  ## Prior for random efect of location 

  #---------------------------------

  

  for (l in 1:Nloc){

    

    loc.eta[l] ~ dnorm(loc.muhat,loc.tau)

    

  }

  

  loc.muhat ~ dunif(-10,10)

  loc.tau <- 1/pow(loc.sigma,2) ## convert stddev to precision (sigma defiend above)

  loc.sigma ~ dunif(0,5)

  

  ## Prior for random efect of site 

  #---------------------------------

  

  for (s in 1:Nsite){

    

    site.eta[s] ~ dnorm(site.muhat,site.tau)

    

  }

   

  site.muhat ~ dunif(-10,10)

  site.tau <- 1/pow(site.sigma,2) ## convert stddev to precision (sigma defiend above)

  site.sigma ~ dunif(0,5)

  

  ## Prior for random effect of diver 

  #-----------------------------------

  

  for (k in 1:Ndiver){

    

    diver.eta[k] ~ dnorm(diver.muhat,diver.tau)

    

  }

  

  diver.muhat ~ dunif(-10,10)

  diver.tau <- 1/pow(diver.sigma,2) ## convert stddev to precision (sigma defiend above)

  diver.sigma ~ dunif(0,5)

  

  #-------------------

  # OBSERVATION MODEL

  #------------------- 

  

  for (l in 1:Nloc){

    for(i in 1:Nobs[l]) {

      exp.mu[i,1,l] <- X[year[l,i],1,l]+loc.eta[loc[l,i]]+site.eta[site[l,i]]+diver.eta[diver[l,i]]

      

      #cumulative logistic probabilities

      logit(Q[i,1,l]) <- tau[l,1]-exp.mu[i,1,l]

      p[i,1,l] <- Q[i,1,l]

      

      for (j in 2:4){

        

        logit(Q[i,j,l]) <- tau[l,j]-exp.mu[i,1,l]

        p[i,j,l] <- Q[i,j,l] - Q[i,j-1,l]

      }

      

      p[i,5,l] <- 1 - Q[i,4,l]

      predY[l,i]  ~ dcat(p[i,1:5,l]) #p[i,] sums to 1 for each i 

      Y[l,i] ~ dnorm(predY[l,i],tauR[l])

      

    }

  }

  #-------------------

  # PROCESS MODEL

  #------------------- 

  for (t in 2:Nyear){

    for (l in 1:Nloc){

      predX[t,1,l] <- X[t-1,1,l]+step[l]

    }

    X[t,1,1:Nloc] ~ dmnorm(predX[t,1,1:Nloc],tauQ[1:Nloc, 1:Nloc])

  }

  

})

  

m <- nimbleModel(nc, constants = list(Nloc=6,Nsite=539,Ndiver=1572,Nyear=27,Nobs=c(2170,977,839,696,875,683),year=constants_list[[1]]$year,loc=constants_list[[1]]$loc,site=constants_list[[1]]$site,diver=constants_list[[1]]$diver,Y=constants_list[[1]]$Y)

 

I realize you cannot reproduce some of the constants, but thought post this in the mean time while I try to figure out how to reproduce the data. Maybe you can see something blaring that I have overlooked. 

Thanks,

Jordan 

[Perry de Valpine]

Please try:

options(error = recover)

# run code that leads to the error, and with luck you will be thrown into a browser at the error.  Then:

varInfo[[iV]]

That should display information about the variable causing the problem.

For a reproducible example, you can send it offline if sharing the data publicly is an issue.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/662c829e-e9ee-450b-a503-3e1430485e1fn%40googlegroups.com.

[Jordan Dinardo]

Hi Perry, 

I will try the code you suggested to pinpoint the variable that is causing the problem. It might be best to send a reproducible example offline so I will shoot you an email.

Best,

Jordan 

[Jordan Dinardo]

Hi Perry, 

I sent you an email offline with data and code. Even after running the suggested code above I am still having trouble deciphering the issue. I'd appreciate any insight on this.

Thanks, 

Jordan 

[Perry de Valpine]

Just reporting to the list that we determined the issue was how some of the inputs were set up.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/edefc800-8a5b-432b-ac12-c53e5a4b1c64n%40googlegroups.com.