Question on Error in chol.default(model$taur[1:2, 1:2])

[MPrzys]

Dear Nimble users,

I try to fit a model for the system x environment x variety data. The model can be written as 

 y[i]~dnorm(mu[i],tau[sys[i]]) , where mu[i] is equal to: 

           mu[i]<-b[env[i],rep[i],sys[i]]+va[v[i],sys[i]]+ge[env[i],v[i],sys[i]].

I fitted the model in WinBUGS (see the code below)

library(R2WinBUGS)
library(coda)


trial<-read.table("D:/starykomp/programy/old/convvsorgZIK/proba/proba.txt",h=TRUE)
#trial
#length(trial$yi)

N<-length(trial$yield)
n<-length(trial$yield)/2
ne<-6
nv<-11
nr<-3
#### organic
y1<-trial$yi[1:n]/10
rep1<-trial$rep[1:n]
si1<-trial$si[1:n]
v1<-trial$var[1:n]
va1<-var(y1);va1


#### conventional
y2<-trial$yi[(n+1):N]/10
rep2<-trial$rep[(n+1):N]
si2<-trial$si[(n+1):N]
v2<-trial$var[(n+1):N]
va2<-var(y2);va2

R<-matrix(nrow=2,ncol=2)
R[1,1]<-va1/2
R[1,2]<-0
R[2,1]<-0
R[2,2]<-va2/2
Rinv<-solve(R)

R2<-diag(1,2)

y<-trial$yield/10
sys<-trial$sys
env<-trial$env
rep<-trial$rep
v<-trial$var

model<-function(){
          for(i in 1:N){
           y[i]~dnorm(mu[i],tau[sys[i]])
           mu[i]<-b[env[i],rep[i],sys[i]]+va[v[i],sys[i]]+ge[env[i],v[i],sys[i]]
           }
for(l in 1:2){a[l]~dnorm(0,0.00001)}
for(l in 1:2){tau[l]~dgamma(2,10)
sig2[l]<-1/tau[l]}
taur[1:2,1:2]~dwish(R2[1:2,1:2],2)
tauv[1:2,1:2]~dwish(R2[1:2,1:2],4)
tauge[1:2,1:2]~dwish(R2[1:2,1:2],2)
mu1[1]<-0
mu1[2]<-0
for(j in 1:ne){
  for(r in 1:nr){
   b[j,r,1:2]~dmnorm(a[1:2],taur[1:2,1:2])#taur[1:2,1:2]
   }
}
for(k in 1:nv){
  va[k,1:2]~dmnorm(mu1[1:2],tauv[1:2,1:2])
 }
covV[1:2,1:2]<-inverse(tauv[1:2,1:2])

for(j in 1:ne){
  for(k in 1:nv){
    ge[j,k,1:2]~dmnorm(mu1[1:2],tauge[1:2,1:2])
  }
 }
covGE[1:2,1:2]<-inverse(tauge[1:2,1:2])
covRep[1:2,1:2]<-inverse(taur[1:2,1:2])
}


model.file<-file.path("D:/starykomp/programy/old/convvsorgZIK","model.txt")
write.model(model,model.file)

data<-list("ne","nv","nr","N","rep","sys","env","v","y","R2")#"Rinv",","T"
parameters<-c("a","covV","covRep","covGE","sig2")#"b","v",


model<-bugs(data,inits=NULL,parameters,model.file,
                n.chains=3,n.iter=1000000,n.burnin=500000,n.thin=25,
                codaPkg=FALSE,debug=TRUE,DIC=TRUE,
                bugs.directory = "D:/Programy/winbugs143/winbugs14_full_patched/WinBUGS14/",
                working.directory="D:/starykomp/programy/old/convvsorgZIK")#,bugs.seed=1231

print(model,digits=3)

I tried to fit the same model in nimble (the code below) 

library(nimble)
library(coda)
library(R2WinBUGS)



trial<-read.table("D:/starykomp/programy/old/convvsorgZIK/proba/proba.txt",h=TRUE)
trial
#length(trial$yi)

### definition of factor in the data set
trial$sys<-as.factor(trial$sys)
#trial$system<-as.factor(trial$system)
#trial$year<-as.factor(trial$year)
#trial$site<-as.factor(trial$site)
trial$var<-as.factor(trial$var)
trial$rep<-as.factor(trial$rep)
trial$env<-as.factor(trial$env)
#trial$s<-as.factor(trial$s)


N<-length(trial$yield)
n<-N/2
ne<-6
nv<-11
nr<-3
me1<-aggregate(yield/10~var,data=trial[1:n,],na.rm=TRUE,FUN="mean")
va1<-var(me1[,2])
me2<-aggregate(yield/10~var,data=trial[(n+1):N,],na.rm=TRUE,FUN="mean")
va2<-var(me2[,2])

R<-matrix(nrow=2,ncol=2)
R[1,1]<-va1/2
R[1,2]<-0
R[2,1]<-0
R[2,2]<-va2/2
Rinv<-solve(R)

R2<-diag(1,2)

y<-trial$yield/10;y
sys<-trial$sys;sys
env<-trial$env
rep<-trial$rep
v<-trial$var;v

model<-nimbleCode({
          for(i in 1:N){
           y[i]~dnorm(mu[i],tau[sys[i]])
           mu[i]<-b[env[i],rep[i],sys[i]]+va[v[i],sys[i]]+ge[env[i],v[i],sys[i]]
           }
for(l in 1:2){a[l]~dnorm(0,0.00001)}
for(l in 1:2){tau[l]~dgamma(2,10)
sig2[l]<-1/tau[l]}
cc[1]<-1
cc[2]<-2
R2[1:2,1:2]<-diag(cc[1:2])
taur[1:2,1:2]~dwish(R=R2[1:2,1:2],df=3)
tauv[1:2,1:2]~dwish(R=R2[1:2,1:2],df=4)
tauge[1:2,1:2]~dwish(R=R2[1:2,1:2],df=4)
mu1[1]<-0
mu1[2]<-0
for(j in 1:ne){
  for(r in 1:nr){
   b[j,r,1:2]~dmnorm(m=a[1:2],prec=taur[1:2,1:2])#taur[1:2,1:2]
   }
}
for(k in 1:nv){
  va[k,1:2]~dmnorm(mu1[1:2],prec=tauv[1:2,1:2])
 }
covV[1:2,1:2]<-inverse(tauv[1:2,1:2])

for(j in 1:ne){
  for(k in 1:nv){
    ge[j,k,1:2]~dmnorm(mu1[1:2],tauge[1:2,1:2])
  }
 }
covGE[1:2,1:2]<-inverse(tauge[1:2,1:2])
covRep[1:2,1:2]<-inverse(taur[1:2,1:2])
}
)

modelConsts<-list(N=N,ne=ne,nv=nv,nr=nr,rep=c(rep),env=c(env),sys=c(sys),v=c(v))
modelData<-list(y=c(y))

mod <- nimbleModel(code = model, name = "model", constants = modelConsts,
data = modelData)
myMCMC <- buildMCMC(mod)
myModel <- compileNimble(mod)

mcmc.out <- nimbleMCMC(code = model, constants = modelConsts,
data = modelData,
nchains = 2, niter = 10000,nburnin=2000,
summary = TRUE, WAIC = TRUE,
monitors = c('a','sig2'),
samples=TRUE,samplesAsCodaMCMC=TRUE).

When I've run the model I've obtained the following error:

 "> mod <- nimbleModel(code = model, name = "model", constants = modelConsts,

+ data = modelData)
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) ... Error in chol.default(model$taur[1:2, 1:2]) :
  the leading minor of order 1 is not positive definite

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."

I assumed Wishart distribution for all precision matrices in my model. Could You tell me where I've made a mistake or what is wrong in my model?

I enclosed my toy data set. 

I would be grateful for any help in solving my problem.

Best regards 

Marcin

[Daniel Turek]

Marcin, I think the issue comes down to a lack of initialization of the multivariate normal precision matrices, specifically:

taur[1:2,1:2]

tauv[1:2,1:2]

tauge[1:2,1:2]

Try adding initial values for these, which would be any positive definite 2x2 matrices.  The identity matrix would be fine.  This will give them valid initial states, but they will still undergo MCMC sampling.  Notice that in NIMBLE, fixed "data" is different from "initial values".  You can add these initial values as shown below (new code is in blue):

modelConsts<-list(N=N,ne=ne,nv=nv,nr=nr,rep=c(rep),env=c(env),sys=c(sys),v=c(v))
modelData<-list(y=c(y))

modelInits <- list(taur = diag(2), tauv = diag(2), taure = diag(2))
mod <- nimbleModel(code = model, name = "model", constants = modelConsts, data = modelData, inits = modelInits)

You could also consider providing neutral initial values for the other model parameters, as well, such as "a", "tau", and others.

Give this a go, and see if it works now.

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/697259c9-de1e-4a8d-b8e0-a356b1c1bb72n%40googlegroups.com.

[MPrzys]

Hi Daniel 

Thanks for the tip. Now the model works. I've got further question on Error in chol.default(model$taur[1:2, 1:2]). Can this error be avoided by defining the prior for precision matrices in terms of Cholesky decomposition? Does nimble do this automatically? 

I have further question which maybe silly. How to define initial values for more than one chain. I tried to use list of list of initial values as in WinBUGS

### initial values
inits1<-list(taur = diag(2), tauv = diag(2), tauge = diag(2),tau=c(0.01,0.01))
inits2<-list(taur= diag(c(0.01,1000),2),tauv=diag(c(100,0.1),2),tauge=diag(c(1000,1000),2),tau=c(10,100))
inits3<-list(taur=X,tauv=diag(c(0.001,1000),2),tauge=diag(c(1000,0.001),2),tau=c(0.001,0.001))

modelInits <- list(inits1,inits2,inits3)

mod <- nimbleModel(code = model, name = "model", constants = modelConsts,
data = modelData,inits = modelInits )

but I get the following error

"Defining model
Error in if (initName %in% names(dL)) { : argument is of length zero}

Best regards

Marcin

[Daniel Turek]

Yes, briefly, under the hood NIMBLE transforms model distributions to a "canonical" parmeterization.  For the dmnorm distribution, this canonical parameterization is always in terms of the Cholesky decomposition of either the covariance matrix, or the precision matrix (depending on which was provided in the original dmnorm declaration.

That is, if you specify precision, then NIMBLE changes this line:

... ~ dmnorm(mu[1:k], Q[1:k,1:k])

into this:

... ~ dmnorm(mu[1:k], cholesky = chol(Q[1:k,1:k]), prec_param = 1)   ## prec_param = 1 indicates that the "cholesky" argument is the factorization of the precision

Or, if you specify dmnorm in terms of covariance, then NIMBLE changes this line:

... ~ dmnorm(mu[1:k], cov = Sigma[1:k,1:k])

into this:

... ~ dmnorm(mu[1:k], cholesky = chol(Sigma[1:k,1:k]), prec_param = 0)   ## prec_param = 0 indicates that the "cholesky" argument is the factorization of the covariance

So, yes, if you want  you can specify the dmnorm model declaration directly using either canonical parameterization, just provide the mean vector, the (upper) triangular cholesky factor matrix as the "cholesky" parameter, and then provide prec_param = 0 or 1, depending on whether the "cholesky" argument is the factorization of the precision (1) or the covariance (0) matrix.

As for multiple initial values: the list you supply to nimbleModel must only be a single list of initial values - for example, your "inits1" object.  Nimble will create the (single) model using these initial values.  For supplying different initial values for multiple chains, you would subsequently pass your "modelInits" (a list of lists of initial values) as the "inits" argument to either nimbleMCMC() or runMCMC() function, also specifying the nchains argument.  A good workflow for you might be to use nimbleModel to create the model, then buildMCMC to build an MCMC, then using compileNimble() to compile both the model and the MCMC, and finally runMCMC (along with the list of different initial values) to run the compiled MCMC.

Hope this helps, keep at it,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/acbda5d0-377e-48bd-b92e-42955ab6695bn%40googlegroups.com.

[MPrzys]

Hi Daniel,

Thanks for the tips. For the same data set, I've rewrote the R code as You suggested. I have a quick question about runMCMC. I defined my model using nimbleModel function with one set of initial values. But in runMCMC, I've put inits=Inits which is list of three lists of initial values for each chain. Does runMCMC use the latter set of initial values?

Further, instead of using invWishart distribution, I tried to use the LKJ distribution  to define the covariance matrix for some random effect in the model. I used the from User manual, but  I've obtained the following error

===== Monitors =====
thin = 1: b, covEV, sdv, tau, Ustar
===== Samplers =====
RW_block_lkj_corr_cholesky sampler (1)
  - Ustar[1:2, 1:2]
RW_block sampler (68)
  - va[]  (11 multivariate elements)
  - ge[]  (57 multivariate elements)
RW sampler (2)
  - sdv[]  (2 elements)
posterior_predictive sampler (9)
  - ge[]  (9 multivariate elements)
conjugate sampler (39)
  - tau[]  (2 elements)
  - b[]  (36 elements)
  - covEV[1:2, 1:2]
Error in samplerFunction(model = model, mvSaved = mvSaved, target = target,  :
  RW_block_lkj_corr_cholesky sampler requires target node dimension to be at least 2x2.

In the attachments I enclosed the R-code. I'm thinking  of using scaled-inverse-Wishart. Should it be defined in the same manner?

Finally, I was also experimenting with runCrossValidate function for model with invWishart for covV. But I've obtained the following error

dropping data fold 1
Checking model sizes and dimensions
  [Note] 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).

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
Error in as.matrix(C.modelMCMC$mvSamples)[, leaveOutNames, drop = FALSE] :
  subscript out of bounds

Can You or anyone explain me what mistake I've made?

Best regards

Marcin

[Chris Paciorek]

I'll let Daniel respond to the other issues, but the LKJ error is one that someone pointed out recently and we've just fixed it on the devel branch on GitHub. We'll have that fix in our next release, which should be soon.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/2da57b66-5a75-47f8-96d6-74fb2649df26n%40googlegroups.com.

[Daniel Turek]

Marcin, thanks for the questions.

As for initial values, those (single set of initial values) that you supply at the time of model building (nimbleModel) are used as the initial values in the model object.  Subsequently, regardless of what was supplied at the time of model building, if you use runMCMC with nchains = N, and provide the argument (to runMCMC) inits = initsList, where initsList is a list (of length N) of lists of initial values, then at the beginning of chain i, runMCMC will set the values into the model to be initsList[[i]], overwriting whatever previous values were in the model (which were either provided at model build time, or otherwise would be where the previous MCMC chain ended at).

As for the runCrossValidate, it's difficult to say with just that small bit of output.  Are you able to provide a fully reproducible example, demonstrating the error from runCrossValidate?

Thanks,

Daniel

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

[MPrzys]

Dear Chris and Daniel, 

thanks for answering my questions. As for runCrossValidate, I've put my nimble code below and enclosed in the attachments the toy data set 

library(nimble)
library(coda)
library(R2WinBUGS)
library(MCMCpack)
#library(MCMCglmm)

trial<-read.table("D:/starykomp/programy/old/convvsorgZIK/proba/proba.txt",h=TRUE)
#trial
#length(trial$yi)

### definition of factor in the data set
trial$sys<-as.factor(trial$sys)
#trial$system<-as.factor(trial$system)
#trial$year<-as.factor(trial$year)
#trial$site<-as.factor(trial$site)
trial$var<-as.factor(trial$var)
trial$rep<-as.factor(trial$rep)
trial$env<-as.factor(trial$env)
#trial$s<-as.factor(trial$s)


N<-length(trial$yield)
n<-N/2
ne<-6
nv<-11
nr<-3
me1<-aggregate(yield/10~var,data=trial[1:n,],na.rm=TRUE,FUN="mean")

va1<-var(me1[,2]);va1

me2<-aggregate(yield/10~var,data=trial[(n+1):N,],na.rm=TRUE,FUN="mean")

va2<-var(me2[,2]);va2
s1<-2.5*sqrt(va1);
s2<-2.5*sqrt(va2)


S<-matrix(nrow=2,ncol=2)
S[1,1]<-va1/2
S[1,2]<-0
S[2,1]<-0
S[2,2]<-va2/2

#Rinv<-solve(R);Rinv


R2<-diag(1,2)

y<-trial$yield/10#;y
sys<-trial$sys#;sys
env<-trial$env
rep<-trial$rep
v<-trial$var#;v


uppertri_mult_diag <- nimbleFunction(
run = function(mat = double(2), vec = double(1)) {
returnType(double(2))
p <- length(vec)
out <- matrix(nrow = p, ncol = p, init = FALSE)
for(i in 1:p)
out[ , i] <- mat[ , i] * vec[i]
return(out)
})

#### Model definition

model<-nimbleCode({
          for(i in 1:N){

           y[i]~dnorm(mu[i],tau=tau[sys[i]])#dt(mu[i],tau=tau[sys[i]],df=7)
           mu[i]<-b[rep[i],env[i],sys[i]]+va[v[i],sys[i]]+ge[env[i],v[i],sys[i]]

           }
for(l in 1:2){tau[l]~dgamma(2,10)
sig2[l]<-1/tau[l]}
cc[1]<-1

cc[2]<-1
S2[1:2,1:2]<-diag(cc[1:2])
covV[1:2,1:2]~dinvwish(S=S[1:2,1:2],df=3)
covEV[1:2,1:2]~dinvwish(S=S2[1:2,1:2],df=4)

mu1[1]<-0
mu1[2]<-0

for(r in 1:nr){
  for(j in 1:ne){
    for(l in 1:2){
   b[r,j,l]~dnorm(0,tau=0.000001)
   }
}
}


for(k in 1:nv){
  va[k,1:2]~dmnorm(mu1[1:2],cov=covV[1:2,1:2])

 }

for(j in 1:ne){
  for(k in 1:nv){

    ge[j,k,1:2]~dmnorm(mu1[1:2],cov=covEV[1:2,1:2])
  }
 }
}
)

X<-riwish(4,diag(2))
Z1<-riwish(3,diag(2))
Z2<-riwish(3, diag(c(0.001,1000),2));Z2
Z3<-riwish(3,diag(c(500,0.1),2));Z3# diag(2)

#### Data and constants

modelConsts<-list(N=N,ne=ne,nv=nv,nr=nr,rep=c(rep),env=c(env),sys=c(sys),v=c(v))
modelData<-list(y=c(y),S=S)

### initial valuestauR = c(0.01,1000), tauR= c(100,0.1),tauR=c(1000,0.001),

inits1<-list(covV = diag(2),covEV = diag(2) ,tau=c(0.01,0.01))#
inits2<-list(covV=diag(c(100,0.1),2),covEV=Z2,tau=c(10,100))#
inits3<-list(covV=X,covEV=Z3,tau=c(0.001,0.001))#

Inits<- list(inits1,inits2,inits3)


#### MCMC setup
Rmodel <- nimbleModel(code = model, constants = modelConsts,data = modelData,inits=inits1)

myModel <- compileNimble(Rmodel)
myMCMC <- buildMCMC(Rmodel)
myMCMCconf<-configureMCMC(Rmodel,enableWAIC=TRUE,print=TRUE)
myMCMCconf$addMonitors(c("sig2","covV","covEV","tau","b","va","ge"))
CmyModel <- buildMCMC(myMCMCconf)
CmyMCMC <- compileNimble(CmyModel,project=myModel)

#### MCMC
mcmc.out <- runMCMC(CmyMCMC, niter = 20000, nburnin=10000,thin=40,
nchains = 3,inits = Inits,
summary = TRUE, progressBar=TRUE,
samples=TRUE,samplesAsCodaMCMC=TRUE)

#### MCMC summary statistics and WAIC

print(mcmc.out$summary$all.chains,digit=4)
calculateWAIC(CmyMCMC)
runCrossValidate(MCMCconfiguration=myMCMCconf,
                        k=4,
                        foldFunction="random",
                        lossFunction="predictive",
                        MCMCcontrol = list(niter = 5000,
                  nburnin = 500))

#calculateWAIC(mcmc.out,CmyMCMC)

###### Samples
mcmc1<-as.mcmc(mcmc.out$samples$chain1)#
mcmc2<-as.mcmc(mcmc.out$samples$chain2)
mcmc3<-as.mcmc(mcmc.out$samples$chain3)
mcmcb<-mcmc.list(mcmc1,mcmc2,mcmc3)

###### Diagnostic plots
gelman.diag(mcmcb,multivariate=FALSE)
gelman.plot(mcmcb,ask=TRUE)
plot(mcmcb,ask=TRUE)

When I've run the code I've obtained the following error

 

> runCrossValidate(MCMCconfiguration=myMCMCconf,

+ k=4,
+ foldFunction="random",
+ lossFunction="predictive",
+ MCMCcontrol = list(niter = 5000,
+                   nburnin = 500))

dropping data fold 1
Checking model sizes and dimensions
  [Note] 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).
Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
Error in as.matrix(C.modelMCMC$mvSamples)[, leaveOutNames, drop = FALSE] :
  subscript out of bounds

I would be grateful for explaining me where I've made the mistake

Best regards

Marcin  

 

[Daniel Turek]

Marcin, thanks for sending the reproducible example.  This revealed a bug in nimble's cross validation routine.  In the presence of multivariate model parameters, runCrossValidate was using the wrong form of the node names when trying to extract posterior samples from the cross validation folds.  Fortunately, this bug would not lead anyone to incorrect results, but rather an error (as you experienced).

I tried to create a fix for this, and put it into a new GitHub branch called "fixCrossVal".  Can you try re-installing nimble from this GitHub branch, and re-running your code?  You can install nimble from this fixed branch using:

remove.packages('nimble')
library(devtools)
install_github('nimble-dev/nimble', ref = 'fixCrossVal', subdir = 'packages/nimble')
library(nimble)

Please give this a try, and let me know if it works for you.  And also be aware, the cross validation does take a while to run.

Cheers,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/b4e8c672-76a2-40ef-91bc-aea8dc5c3913n%40googlegroups.com.