NA warning with compile

[Beth Ross]

Hi all,
I'm new to Nimble (just attended the recent workshop!). Post-workshop, I wanted to convert my recent JAGS code to Nimble to see if it would help with convergence issues. I've been working on a problem with two levels of variable selection - one for the covariate itself and one for the "best of 3" spatial or temporal scales ("scale" below) given the covariate is selected. I've tried to implement the code using RJMCMC instead of variable selection which generally seems to work well - up to a certain number of covariates. If I use more than 10 covariates (out of 14), I get the following warning message when compiling the model:

In if (is.na(code)) stop("Character NA is not supported for compilation.") :
the condition has length > 1 and only the first element will be used

I understand what this message means in R, but I've been stymied trying to chase down the problem in this code. I've provided initial values for everything I can, and the calculated values from btpi.model (code below) all seem reasonable for multiple different initial values. I've tried different values for covariates which also works, so long as there are fewer than 10 (and triple checked all data for NAs!). I also don't get the warning message if I use the full model with the other form of RJMCMC without indicator variables, but given the quadratic terms in the model, I wanted to keep these in. I've also tried constraining the indicator variables based on a previous post discussing RJMCMC and interaction terms, but the constraint did not change the warning message. I get considerably different results if I use 10 or all 14 covariates, so want to be sure the warning message isn't indicative of a major problem.

Thank you for any help or advice! And for the workshop last week, it was incredibly helpful!

Best,
Beth


> btpi_code <- nimbleCode({
+ #priors
+
+ beta0 ~ dnorm(0,.1)
+ size ~ dunif(0,10)
+ psi <- 0.5
+
+ for(s in 1:ncov){
+ indA[s] ~ dbern(psi)
+ beta[s] ~ dnorm(0,.1)
+ }
+
+ prob.c[1:3] <- c(1/3, 1/3, 1/3)
+ for(i in 1:7){
+ scale[i] ~ dcat(prob.c[1:3])}
+
+ #likelihood
+ for(i in 1:nsite){
+ log.mu[i] <- beta0 +
+ indA[1]*beta[1]*covs[i,scale[1],1] +
+ indA[1]*indA[2]*beta[2]*(covs[i,scale[1],1]^2) +
+ indA[3]*beta[3]*covs[i,scale[2],2] +
+ indA[3]*indA[4]*beta[4]*(covs[i,scale[2],2]^2) +
+ indA[5]*beta[5]*covs[i,scale[3],3] +
+ indA[5]*indA[6]*beta[6]*(covs[i,scale[3],3]^2) +
+ indA[7]*beta[7]*covs[i,scale[4],4] +
+ indA[7]*indA[8]*beta[8]*(covs[i,scale[4],4]^2) +
+ indA[9]*beta[9]*covs[i,scale[5],5] +
+ indA[9]*indA[10]*beta[10]*pow(covs[i,scale[5],5],2) +
+ indA[11]*beta[11]*covs[i,scale[6],6] +
+ indA[11]*indA[12]*beta[12]*(covs[i,scale[6],6]^2) +
+ indA[13]*beta[13]*covs[i,scale[7],7] +
+ indA[13]*indA[14]*beta[14]*(covs[i,scale[7],7]^2)
+
+ log.mu.lim[i] <- min(250, max(-250, log.mu[i]))
+ mu[i] <- exp(log.mu.lim[i])
+
+ y[i] ~ dnegbin(prob = p[i],size = size)
+ p[i] <- size/(mu[i]+size)
+
+ #fit stats
+
+ E[i] <- pow(sqrt(y[i])-sqrt((size*(1-p[i]))/p[i]),2)
+ y.pred[i] ~ dnbinom(p[i],size)
+ E.pred[i] <- pow(sqrt(y.pred[i])-sqrt((size*(1-p[i]))/p[i]),2)
+ }
+
+ fit <- sum(E[1:nsite])
+ fit.new <- sum(E.pred[1:nsite])
+
+ })
> dim(cov.scale)
[1] 1932 31
> covs = array(unlist(cov.scale[,-c(1,17:19,26:31)]),
+ dim=c(nsite,3,7))
> ncov = dim(covs)[3]*2
> win.data <- list(y = counts,covs = covs)
> inits <- list(beta0 = runif(1,-1,1),
+ size = 1,
+ scale = rep(1,7),
+ indA = rep(0,14),
+ beta = rep(0,14),
+ y.pred = counts)
> btpi.model <- nimbleModel(code = btpi_code,
+ constants = list(nsite=nsite,
+ ncov=14),
+ data = win.data,
+ inits = inits)
defining 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...
model building finished.
> btpi.config.mcmc <- configureMCMC(btpi.model)
===== Monitors =====
thin = 1: beta0, size, beta, indA, scale
===== Samplers =====
RW sampler (16)
- beta0
- size
- beta[] (14 elements)
posterior_predictive sampler (1932)
- y.pred[] (1932 elements)
categorical sampler (7)
- scale[] (7 elements)
binary sampler (14)
- indA[] (14 elements)
> btpi.config.mcmc$addMonitors(c('indA', 'fit','fit.new'))
thin = 1: beta0, size, beta, indA, scale, fit, fit.new
> configureRJ(btpi.config.mcmc,
+ targetNodes = 'beta',
+ indicatorNodes = 'indA',
+ control = list(mean=0, scale = .2))
> btpi.mcmc <- buildMCMC(btpi.config.mcmc)
> btpi.cmodel <- compileNimble(btpi.model)
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
compilation finished.
Warning message:
In if (is.na(code)) stop("Character NA is not supported for compilation.") :
the condition has length > 1 and only the first element will be used
> btpi.model$mu[1:10]
[1] 1.546792 1.546792 1.546792 1.546792 1.546792 1.546792
[7] 1.546792 1.546792 1.546792 1.546792
> btpi.model$p[1:10]
[1] 0.3926509 0.3926509 0.3926509 0.3926509 0.3926509
[6] 0.3926509 0.3926509 0.3926509 0.3926509 0.3926509

[Daniel Turek]

Beth, good question, and I'm glad you were concerned by this warning.  Unfortunately, this relates to some internal code processing that's going on inside the NIMBLE compiler, which isn't working as it should.  The problem is triggered by the lengthy model declaration line that defines log.mu[i].  Essentially, it's "too long" in terms of the sheer number of characters, and somewhere inside NIMBLE processing, this gets split into 2 chunks, which is then causing the warning you're seeing - and perhaps also other undesired behaviour.

I'm going to flag this for the development team to look into.  This is something we should address and fix.  But in the meantime, you can easily work around this by just splitting that line into two, for example as:

        temp1[i] <- beta0 +

            indA[1]*beta[1]*covs[i,scale[1],1] +

            indA[1]*indA[2]*beta[2]*(covs[i,scale[1],1]^2) +

            indA[3]*beta[3]*covs[i,scale[2],2] +

            indA[3]*indA[4]*beta[4]*(covs[i,scale[2],2]^2) +

            indA[5]*beta[5]*covs[i,scale[3],3] +

            indA[5]*indA[6]*beta[6]*(covs[i,scale[3],3]^2) +

            indA[7]*beta[7]*covs[i,scale[4],4] +

            indA[7]*indA[8]*beta[8]*(covs[i,scale[4],4]^2)

        temp2[i] <- indA[9]*beta[9]*covs[i,scale[5],5] +

            indA[9]*indA[10]*beta[10]*pow(covs[i,scale[5],5],2) +

            indA[11]*beta[11]*covs[i,scale[6],6] +

            indA[11]*indA[12]*beta[12]*(covs[i,scale[6],6]^2) +

            indA[13]*beta[13]*covs[i,scale[7],7] +

            indA[13]*indA[14]*beta[14]*(covs[i,scale[7],7]^2)

        log.mu[i] <- temp1[i] + temp2[i]

This change won't affect the effective model structure or the RJMCMC at all, and will have a truly minimal effect on the speed of model calculations - which won't be noticeable.  Give it a try, and see if it works.  And it was just pure coincidence that the threshold for when this became "too long" was after the first 10 covariates.

Hope this helps.

Cheers,

Daniel

--
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/4c839ef7-23f9-4730-bc58-71dbc7ec667fn%40googlegroups.com.