using NIMBLE with jagam

[dbowler]

I keep running into the same error message when trying to fit a model using code developed by jagam in NIMBLE. The model runs in JAGS but not NIMBLE. Its probably simple but I am new to NIMBLE...

Help appreciated!

#simulate data
require(mgcv)

set.seed(2) 
n <- 400
dat <- gamSim(1,n=n,dist="normal",scale=2)

#get spline code
jd <- jagam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,
            file="jagam.txt",
              sp.prior="gamma")
 
#first fit using jags as standard
library(rjags)
library(jagsUI)
out <- jags(data = jd$jags.data,
            parameters.to.save = c("mu"),
            model.file = "jagam.txt",
            n.thin=2, n.chains=3, n.burnin=1000,n.iter=2000,
            parallel=T)
#works ok

#also using NIMBLE
library(nimble)
library(igraph)
library(coda)

#sort data
nimbleConstants <- jd$jags.data[c("n","zero")]
nimbleData <- jd$jags.data[c("y","X","S1","S2","S3","S4")]

#model code - using jagam.txt above - except adding in indices
mynimbleCode <- nimbleCode({
 
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i],tau)
    mu[i] <- X[i,1:37] * b[1:37]
  }
  tau ~ dgamma(.05,.005)
 
  ## Parametric effect
  for (i in 1:1) {
    b[i] ~ dnorm(0,0.01)
  }
 
  ## prior for s(x0)...
  K1[1:9,1:9] <- S1[1:9,1:9] * lambda[1]  + S1[1:9,10:18] * lambda[2]
  b[2:10] ~ dmnorm(zero[2:10],K1[1:9,1:9])
  ## prior for s(x1)...
  K2[1:9,1:9] <- S2[1:9,1:9] * lambda[3]  + S2[1:9,10:18] * lambda[4]
  b[11:19] ~ dmnorm(zero[11:19],K2[1:9,1:9])
  ## prior for s(x2)...
  K3[1:9,1:9] <- S3[1:9,1:9] * lambda[5]  + S3[1:9,10:18] * lambda[6]
  b[20:28] ~ dmnorm(zero[20:28],K3[1:9,1:9])
  ## prior for s(x3)...
  K4[1:9,1:9] <- S4[1:9,1:9] * lambda[7]  + S4[1:9,10:18] * lambda[8]
  b[29:37] ~ dmnorm(zero[29:37],K4[1:9,1:9])
 
  ## smoothing parameter priors
  for (i in 1:8) {
    lambda[i] ~ dgamma(.05,.005)
    rho[i] <- log(lambda[i])
  }
 
})

#fit model
Rmodel <- nimbleModel(mynimbleCode, constants = nimbleConstants)
#error here:
#Error in chol.default(model$K1[1:9, 1:9]) :
#  the leading minor of order 1 is not positive definite

Thanks for any pointers

Diana

[Daniel Turek]

Diana, no problem at all.  The main concern here is about initialization of variables in your model, but I also noted one thing in your model code which might a mistake.

Since you don't provide initial values for lambda[1:8], they are all taken to be NA.  That causes the calculation of

K1[1:9,1:9] <- S1[1:9,1:9] * lambda[1]  + S1[1:9,10:18] * lambda[2]

and similar terms to be 9x9 matrices of NA's.

Then, nimble re-parameterizes the multivariate normal distributions:

b[2:10] ~ dmnorm(zero[2:10],K1[1:9,1:9])

to (internally) use the cholesky factorization of the precision matrix K1[1:9,1:9], which attempts to calculate this matrix factorization of the 9x9 matrix of NA's.  That's what leads to the error you're seeing.

You can easily fix this by providing any positive initial values for lambda[1:8]:

nimbleInits <- list(lambda = rep(1, 8))
Rmodel <- nimbleModel(mynimbleCode, constants = nimbleConstants, data = nimbleData, inits = nimbleInits)

This gets your model mostly working fine, but then warnings about the dimensionality of mu appear.  I took a look at the model code, and this line was concerning:

mu[i] <- X[i,1:37] * b[1:37]

because the left-hand-side is a *vector* of length 37, which is then assigned into the *scalar* model variable mu[i].  So I think you might want to change that line to:

mu[i] <- sum(X[i,1:37] * b[1:37])

With that change, the model builds fine.

Also noting the model is not fully initialized, and you might consider adding reasonable initial values for tau, and also b[1:37].

Hope this helps, thanks for sending the question,

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/a84141f2-e3ac-4eef-9188-73a2f6c30a25n%40googlegroups.com.

[Perry de Valpine]

I'll just add, in case it helps, that a recent workshop module on spatial models in nimble included a jagam example: https://github.com/nimble-training/nimble-lisbon-2022/tree/main/content/9_spatial_models.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CAKbe0hreQutznQKPugEAZgFcm2M5Br3V6s47DoWY5DxYh2h9rA%40mail.gmail.com.