fitting a model with nimbleMCMC for matrix data

[MPrzys]

Dear Nimble users

I'm trying to fit a very simple model in nimble (see the code below)

library(agridat)

library(nimble)

library(coda)

data(cornelius.maize)

dat <- cornelius.maize

#env<-dat$env

#gen<-c(dat$gen)

#gen

yield<-dat$yield

Y<-matrix(c(yield),nrow=20,ncol=9)

yn<-t(Y)

mm<-mean(y)

#Y;Yn

nv<-nrow(yn);nv

ne<-ncol(yn);ne

sortR <- nimbleRcall(prototype   = function(x = double(1)) {},

                     Rfun       = "sort",

                     returnType = double(1))

mcode<-nimbleCode({

mu~dnorm(0,tau=0.000001)

for(i in 1:nv){alpha[i]~dnorm(0,tau=1/nv)}

for(j in 1:nv){beta[j]~dnorm(0,tau=1/ne)}

for(i in 1:nv){

  for(j in 1:ne){

    y[i,j]~dnorm(m[i,j],tau)

    m[i,j]<-mu+alpha[i]+beta[j]

  }

}

tau~dgamma(0.1,0.1)

})

mConts<-list(ne=ncol(yn),nv=nrow(yn))#,m.mu=mean(y)

mData<-list(y=matrix(c(yn),nrow=nrow(yn))

nimbleModel(code=mcode,data=mData)

mcmc <- nimbleMCMC(code=mcode,constants=mConts,data=mData,

niter=10000,nburnin=2500,

summary=TRUE, samples=FALSE,

monitors=c('mu','alpha','beta'))

mcmc.

The data used in the code is type matrix with 9 row and 20 columns. I've the model code in OpenBUGS and is correct. But when I've run nimbleMCMC in the summary I've obtained NA. I suspect that I didn't correctly define the data.  

For this reason, I'd like to ask if someone could help me?

Best regards

Marcin

[Chris Paciorek]

Hi Marcin,

It looks like you are not setting up priors for all the elements of beta.

You have 

for(j in 1:nv){beta[j]~dnorm(0,tau=1/ne)}

but I think you mean to have "for j in 1:ne".  The issue is that nimble can't initialize all the beta values because you don't have priors for all the elements. 

I'm surprised that code would work in BUGS.

Also, it's often a good idea to provide good initial values and not rely on simulating from the priors. 

-chris

--
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/360ae477-a053-4468-b31b-c2ed4fb40eb9n%40googlegroups.com.

[MPrzys]

Hi Chris,

Thanks for pointing the mistake. Now it works perfectly.

I extended the model to the AMMI model with 2 principal components described in Josse et (2014) JABES. For this purpose I used the nimbleRcall function according to Sally's suggestion but when I run the nimbleMCMC I obtain the following message

 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 model$lam[1:2] <<- sortR(model$lamb[1:2]) : 

  replacement has length zero

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

checking model calculations...

Error in model$lam[1:2] <<- sortR(model$lamb[1:2]) : 

  replacement has length zero

Error in model$lam[1:2] <<- sortR(model$lamb[1:2]) : 

  replacement has length zero

Note: cannot calculate logProb for node  lam[1:2] .

NAs were detected in model variables: mu, logProb_mu, alpha, logProb_alpha, beta, logProb_beta, lamb, U, logProb_U, V, logProb_V, tau, logProb_tau, m, lam, lifted_d1_over_sqrt_oPtau_cP, se2, D, logProb_y.

Infinite values were detected in model variables: logProb_lamb, logProb_U.

model building finished.

compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.

Warning, in eigenizing model_D[1] the [ is still there but nDim is not 0 (not a scalar).

Error: Failed to create the shared library. Run 'printErrors()' to see the compilation errors.

In addition: Warning message:

In model$D[1] <<- diag(model$lam[1:2]) :

  number of items to replace is not a multiple of replacement length

Could You help and explain me what mistake I've made? Below I enclose the modified code.

 library(agridat)

library(nimble)

library(coda)

data(cornelius.maize)

dat <- cornelius.maize

yield<-dat$yield

Y<-matrix(c(yield),nrow=20,ncol=9)

yn<-t(Y)/1000

mm<-mean(y)

nv<-nrow(yn);nv

ne<-ncol(yn);ne

sortR <- nimbleRcall(prototype   = function(x = double(1)) {},

                     Rfun       = "sort",

                     returnType = double(1))

#### nimble code for fitting the AMMI model

#### according to Josse et al. (2014) JABES 19, 240-257.

mcode<-nimbleCode({

mu~dnorm(4.85,sd=1)

for(i in 1:nv){alpha[i]~dnorm(0,tau=1/nv)}

for(j in 1:ne){beta[j]~dnorm(0,tau=1/ne)}

for(q in 1:2){lamb[q]~T(dnorm(0,sd=13),0,)}#### Singular nvalues

lam[1:2]<-sortR(lamb[1:2])

for(q in 1:2){

U[1,q]~T(dnorm(0,sd=1),0,) #### Left singular vectors

for(i in 2:nv){

U[i,q]~dnorm(0,sd=1)

   }

}

for(q in 1:2){    #### Right singular vectors

  for(j in 1:ne){

    V[j,q]~dnorm(0,sd=1)  

  }

}

#for(q in 1:2){

# for(i in 1:nv){

#  W[i,q]<-U[i,q]*lam[q]

# }

#}

D<-diag(lam[1:2])

W<-U%*%D

inte<-W%*%t(V)

for(i in 1:nv){

  for(j in 1:ne){

    y[i,j]~dnorm(m[i,j],tau=tau)

    m[i,j]<-mu+alpha[i]+beta[j]+inte[i,j]

  }

}

tau~dgamma(2,0.1)

se2<-1/tau

})

mConts<-list(ne=ncol(yn),nv=nrow(yn))#,m.mu=mean(y)

mData<-list(y=yn)#matrix(c(yn),nrow=nrow(yn)

nimbleModel(code=mcode,data=mData)

mcmc <- nimbleMCMC(code=mcode,constants=mConts,data=mData,

niter=10000,nburnin=2500,

summary=TRUE, samples=FALSE,

monitors=c('mu','alpha','beta','se2','lam'))

mcmc

Best regards

Marcin

[Chris Paciorek]

hi Marcin,

Nimble requires that you provide brackets when using vectors, matrices, arrays (and commas indicating dimensionality for matrices and arrays).

We also often need explicit indexing. So, for example, in the lines of code where you define D, W, and inte, you'll need to add this information, e.g.,

    D[1:2,1:2]<-diag(lam[1:2])
    W[1:nv,1:2]<-U[1:nv,1:2]%*%D[1:2,1:2]
    inte[1:nv,1:ne]<-W[1:nv,1:2]%*%t(V[1:ne,1:2])

There might be more going on, but that is a good place to start.

-Chris

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/d47fa8b6-3e81-43f0-aa4c-7bd8d02b4f02n%40googlegroups.com.