question about dinterval() and dimension of cutpoints

[Rob Erhardt]

Hello,

I greatly enjoyed the NIMBLE workshop last month, and am struggling with what I think is an imprecise definition of the dimension of an object in NIMBLE.

I am trying to use  ~ dinterval() along with a vector "cut" to handle an ordinal response variable Y, defined from binning a latent variable Z using the vector "cut".  NIMBLE code is shown below, with the relevant bits shown in blue and the resulting error message in red.

DEcode <- nimbleCode({

  #Drought Variable
  for(i in 1:I){
    for(t in 1:J){
      Y[i,t] ~ dinterval(Z[i,t], cut)
    }
  }
  #Latent Gaussian Variable
  for(i in 1:I){
    Z[i,1] ~ dnorm(beta[1] + beta[2]*A[i,1] + beta[3]*P[i,1] + beta[4]*S[i,1] + U[i], tau.z)
    for(t in 2:J){
      Z[i,t] ~ dnorm(beta[1] + beta[2]*A[i,t] + beta[3]*P[i,t] + beta[4]*S[i,t] + rho.z[i]*(Z[i,(t-1)] - (beta[1] + beta[2]*A[i,(t-1)] + beta[3]*P[i,(t-1)] + beta[4]*S[i,(t-1)])) + U[i], tau.z)
    }
  }
 
  # ICAR prior for the spatial random effects
  U[1:I] ~ dcar_normal(adj[], weights[], num[], tau.U, zero_mean=1)
 
  #Covariate Models
  for(i in 1:I){
    A[i,1] ~ dnorm(G[1,i] + delta[1,1]*w[1,1] + delta[1,2]*w[2,1], tau.x[1])
    P[i,1] ~ dnorm(G[2,i] + delta[2,1]*w[1,1] + delta[2,2]*w[2,1], tau.x[2])
    S[i,1] ~ dnorm(G[3,i] + delta[3,1]*w[1,1] + delta[3,2]*w[2,1], tau.x[3])
    for(t in 2:J){
      A[i,t] ~ dnorm(G[1,i] + delta[1,1]*w[1,t] + delta[1,2]*w[2,t] + rho.x[1]*(A[i,(t-1)] - (G[1,i] + delta[1,1]*w[1,(t-1)] + delta[1,2]*w[2,(t-1)])), tau.x[1])
      P[i,t] ~ dnorm(G[2,i] + delta[2,1]*w[1,t] + delta[2,2]*w[2,t] + rho.x[2]*(P[i,(t-1)] - (G[2,i] + delta[2,1]*w[1,(t-1)] + delta[2,2]*w[2,(t-1)])), tau.x[2])
      S[i,t] ~ dnorm(G[3,i] + delta[3,1]*w[1,t] + delta[3,2]*w[2,t] + rho.x[3]*(S[i,(t-1)] - (G[3,i] + delta[3,1]*w[1,(t-1)] + delta[3,2]*w[2,(t-1)])), tau.x[3])
    }
  }
  #Random Effects
  for(i in 1:I){
    for(r in 1:3){G[r,i] ~ dnorm(0, tau.G[r])}
  }
  for(i in 1:I){rho.z[i] ~ T(dnorm(rho.mean, tau.rho),-1,1)}
 
  #Priors
  for(r in 1:3){
    tau.G[r] ~ dgamma(0.01, 0.01)
    tau.x[r] ~ dgamma(0.01, 0.01)
  }
  tau.rho ~ dgamma(0.01, 0.01)
  tau.z ~ dgamma(0.01, 0.01)
 
  sigma ~ dunif(0, 100)   # prior for variance components based on Gelman (2006)
  tau.U <- 1 / sigma^2
 
  for(i in 1:3){rho.x[i] ~ dunif(-1,1)}
  rho.mean ~ dunif(-1,1)
  for(i in 1:3){
    delta[i,1] ~ dnorm(0, 0.01)
    delta[i,2] ~ dnorm(0, 0.01)
  }
  for(i in 1:4){beta[i] ~ dnorm(0, 0.01)}
  for(i in 1:4){d[i] ~ dnorm(0, 0.01)}
 
  cut <- c(0, rep(exp(d[1]), 4)) + c(0, 0, rep(exp(d[2]), 3)) + c(0, 0, 0, rep(exp(d[3]), 2)) + c(0,0,0, 0,exp(d[4]))
 
  #Predict Covariates, Generate Z's, Predict Drought
  for(i in 1:I){
    A.p[i,1] ~ dnorm(G[1,i] + delta[1,1]*w[1,(J+1)] + delta[1,2]*w[1,(J+1)] + rho.x[1]*(A[i,J] - (G[1,i] + delta[1,1]*w[1,J] + delta[1,2]*w[2,J])), tau.x[1])
    P.p[i,1] ~ dnorm(G[2,i] + delta[2,1]*w[1,(J+1)] + delta[2,2]*w[1,(J+1)] + rho.x[2]*(P[i,J] - (G[2,i] + delta[2,1]*w[1,J] + delta[1,2]*w[2,J])), tau.x[2])
    S.p[i,1] ~ dnorm(G[3,i] + delta[3,1]*w[1,(J+1)] + delta[3,2]*w[1,(J+1)] + rho.x[3]*(S[i,J] - (G[3,i] + delta[3,1]*w[1,J] + delta[1,2]*w[2,J])), tau.x[3])
    Z.p[i,1] ~ dnorm(beta[1] + beta[2]*A.p[i,1] + beta[3]*P.p[i,1] + beta[4]*S.p [i,1] + rho.z[i]*(Z[i,J] - (beta[1] + beta[2]*A[i,J] + beta[3]*P[i,J] + beta[4]*S[i,J])), tau.z)
    Y.p[i,1] ~ dinterval(Z.p[i,1], cut)
    for(t in 2:p){
      A.p[i,t] ~ dnorm(G[1,i] + delta[1,1]*w[1,(J+t)] + delta[1,2]*w[2,(J+t)]  + rho.x[1]*(A.p[i,(t-1)] - (G[1,i] + delta[1,1]*w[1,(J+t-1)] + delta[1,2]*w[2,(J+t-1)])), tau.x[1])
      P.p[i,t] ~ dnorm(G[2,i] + delta[2,1]*w[1,(J+t)] + delta[2,2]*w[2,(J+t)]  + rho.x[2]*(P.p[i,(t-1)] - (G[2,i] + delta[2,1]*w[1,(J+t-1)] + delta[2,2]*w[2,(J+t-1)])), tau.x[2])
      S.p[i,t] ~ dnorm(G[3,i] + delta[3,1]*w[1,(J+t)] + delta[3,2]*w[2,(J+t)]  + rho.x[3]*(S.p[i,(t-1)] - (G[3,i] + delta[3,1]*w[1,(J+t-1)] + delta[3,2]*w[2,(J+t-1)])), tau.x[3])
      Z.p[i,t] ~ dnorm(beta[1] + beta[2]*A.p[i,t] + beta[3]*P.p[i,t] + beta[4]*S.p[i,t] + rho.z[i]*(Z.p[i,(t-1)] - (beta[1] + beta[2]*A.p[i,(t-1)] + beta[3]*P.p[i,(t-1)] + beta[4]*S.p[i,(t-1)])), tau.z)
      Y.p[i,t] ~ dinterval(Z.p[i,t], cut)
    }
  }
})

## ends nimble code

DEconstants <- list(I=I,J=J,p=p, adj=adj, num=num, weights=weights)

DEmodel <- nimbleModel(DEcode,
                       constants = DEconstants)

defining model... Adding adj, num, weights as data for building model. building model... setting data and initial values... running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... checking model sizes and dimensions... This model is not fully initialized. This is not an error. To see which variables are not initialized, use model$initializeInfo(). For more information on model initialization, see help(modelInitialization). model building finished.

Warning message: In model$cut[1] <<- nimC(0, nimRep(exp(model$d[1]), 4)) + nimC(0, : number of items to replace is not a multiple of replacement length

"cut" is a vector whose first element is fixed at zero, but elements 2 through 5 will be updated in the MCMC.  If it helps, in JAGS the problem doesn't arise with a simpler version of the model.  However, I think this is because I am able to give a modified argument in the dinterval() command for JAGS, which is unavailable in NIMBLE:

(From JAGS)

  lambda <- rep(exp(d[1]), 4) + c(0, rep(exp(d[2]), 3)) + c(0, 0, rep(exp(d[3]), 2)) + c(0,0,0, exp(d[4]))

  Y[i,t] ~ dinterval(Z[i,t], c(0,lambda)) ## <-- multivariate argument c(0, lambda) allowed in JAGS.

Am I being careless about the dimension of "cut" someplace in the NIMBLE code?  Or perhaps I've broken a NIMBLE rule by trying to update only some elements of the vector "cut" so it isn't a proper node?  Many thanks for any thoughts you can share.  Cheers,

Rob

[Perry de Valpine]

Hi Rob,

Thanks for posting the question. 

nimble models always require square brackets for non-scalars, so please try

cut[1:5] <- c(0, rep(exp(d[1]), 4)) + c(0, 0, rep(exp(d[2]), 3)) + c(0, 0, 0, rep(exp(d[3]), 2)) + c(0,0,0, 0,exp(d[4]))

or

cut[] <- c(0, rep(exp(d[1]), 4)) + c(0, 0, rep(exp(d[2]), 3)) + c(0, 0, 0, rep(exp(d[3]), 2)) + c(0,0,0, 0,exp(d[4]))

along with including the dimensions argument to nimbleModel, as

nimbleModel(<other stuff>, dimensions = list(cut = 5))

If that doesn't work,  I suggest making a tiny toy model with just the cut and dinterval to have a simpler case for debugging.

-Perry

--
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/8a91830c-c971-43e0-8878-f139562566cdn%40googlegroups.com.