dimension error in spatiotemporal model

[Insang Song]

Hello, thank you for maintaining the nimble package!

I am trying to fit a spatiotemporal hierarchical model with nimble code as below:

intIV <- nimbleCode({

        ### Prior ####

        # betas

        beta0 ~ dflat()

        for (k in 1:p)

            betas[k] ~ dflat()

        # sigma

        sigma.U ~ dflat()

        sigma.v ~ dflat()

        sigma.delta ~ dflat()

        sigma.sigma ~ dflat()

        sigma.delta1it ~ dflat()

        sigma.betaU ~ dflat()

        sigma.betav ~ dflat()

        sigma.gamma ~ dflat()

        sigma.betag ~ dflat()

        # precision

        prec.U <- pow(sigma.U, -2)

        prec.v <- pow(sigma.v, -2)

        prec.delta <- pow(sigma.delta, -2)

        tau.sigma <- pow(sigma.sigma, -2)

        tau.delta1it <- pow(sigma.delta1it, -2)

        tau.U <- pow(sigma.U, -2)

        tau.betaU <- pow(sigma.betaU, -2)

        tau.v <- pow(sigma.v, -2)

        tau.betav <- pow(sigma.betav, -2)

        tau.gamma <- pow(sigma.gamma, -2)

        tau.betag <- pow(sigma.betag, -2)

        for (i in 1:N){

            U[i] ~ dnorm(0, tau.U)

            betaU[i] ~ dnorm(0, tau.betaU)

        }

        v[1:N] ~ dcar_normal(adj = tm.adj[1:L], num = tm.num[1:N], tau = tau.v)

        betav[1:N] ~ dcar_normal(adj = tm.adj[1:L], num = tm.num[1:N], tau = tau.betav)

        gamma[1] ~ dnorm(0, 1)

        betag[1] ~ dnorm(0, 1)

        for (t in 2:T){

            gamma[t] ~ dnorm(gamma[t-1], tau.gamma)

            betag[t] ~ dnorm(betag[t-1], tau.betag)

        }

        

        for (t in 2:T) {

            delta_1it[1:N, t-1] ~ dcar_normal(adj = tm.adj[1:L], num = tm.num[1:N], tau = tau.delta1it)

            }

        for (i in 1:N){

            for (t in 2:T){

                epsilonit[i,t] ~ dflat()

                delta_1it[i,t] <- delta_1it[i,t-1] + epsilonit[i,t]

                

            }

            for (t in 1:T){

                epsilon[i,t] ~ dflat()

                beta1[i,t] <- betaU[i] + betav[i] + betag[t] + delta_1it[i, t]

                y[i,t] ~ dnorm(beta0 + (beta1[i, t] * x[(N*(t-1))+i, 1]) + (betas %*% x[(N*(t-1))+i, 2:(p+1)]) + U[i] + v[i] + gamma[t] + epsilon[i,t], tau.sigma)

                

            }

            

        }

})

When I build the model, then it crashes with error Error in genVarInfo3() : 

  Dimension of betas is 1, which does not match its usage in 'y[i, t] ~ dnorm(mean = beta0 + (beta1[i, t] * x[(3108 * (t - 1)) + i, 1]) + (betas %*% x[(3108 * (t - 1)) + i, 2:9]) + U[i] + v[i] + gamma[t] + epsilon[i, t], sd = lifted_d1_over_sqrt_oPtau_dot_sigma_cP, .tau = tau.sigma, lower_ = -Inf, upper_ = Inf, .var = lifted_d1_over_sqrt_oPtau_dot_sigma_cP * lifted_d1_over_sqrt_oPtau_dot_sigma_cP)'.

The dimensions of x is (3108*14, 9), y is (3108, 14), betas are (8), beta1 is (3108, 14). As far as I understood, the dimensions are compatible in the expression (betas %*% x[(3108 * (t - 1)) + i, 2:9]). Could anyone give any idea for the error?

Thank you very much.

[Perry de Valpine]

Hi Insang,

Thanks for your question.  You've landed in some of the issues of types in nimble.  Because nimble generates C++ and C++ is statically typed, nimble requires stricter coding than R to ensure correct handling of types.  Types here refers to the number of indices.  As a result of static typing, matrix multiplication always returns a matrix, but you need a scalar.  In R, a 1x1 matrix may be dynamically treated as a scalar, but in nimble this must be done explicitly.  For example, if you need the result of A[,] %*% B[,] to be a vector, you can do (A[,] %*% B[,] ) [1:nrow, j], and for a scalar you can do (A[,] %*% B[,] )[1,1].  nimble is sometimes able to convert a vector (like betas) to a one-column or one-row matrix, although in your case I am wondering if you wanted x[(3108 * (t - 1)) + i, 2:9] %*% betas?  However, I think there is an easier solution for your need, which is the inprod (inner product) function.  Can you please try inprod(betas[1:8], x[(3108 * (t - 1)) + i, 2:9])?  This will return a scalar in nimble's type management system.  Also, in model code (but not nimbleFunctions), it is required to include square brackets on non-scalars.  That means you need "betas[]" or "betas[1:8]".  In the former usage, it must either be the case that nimble can determine the size of betas from other model code (which appears to be the case in your code) or that the size is provided in the dimensions argument to nimbleModel.

Here are some links that I hope are helpful:

If you are coming from BUGS or JAGS, you may find this conversion guide helpful (indexing brackets are one of the points covered): https://r-nimble.org/quick-guide-for-converting-from-jags-or-bugs-to-nimble

Information about required indexing formats in model code is in Chapter 5 of the nimble User Manual: https://r-nimble.org/html_manual/cha-writing-models.html

Information on types in nimbleFunctions (much of which is relevant for model code) is in Section 11.3 of the User Manual: https://r-nimble.org/html_manual/cha-RCfunctions.html#sec:how-nimble-handles

Good luck and let us know how it goes.

-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/cfbd4012-5b2a-452d-9953-7be67e988e5dn%40googlegroups.com.

[Insang Song]

Thank you for the prompt response!

Using inprod and adding brackets after betas solved my problem. But the other problem occurs when I try to compile and configure a model. When I try to run the code as below, it shows me an error message:

Rmodel <- nimbleModel(intIV, data = dat, constants, inits, returnDef = TRUE)

Error: Some items provided for compilation do not have types that can be compiled: Rmodel.  The types provided were: unknown. Be sure only specialized nimbleFunctions are provided, not nimbleFunction generators.

When I remove the returnDef argument, it compiles without errors. But in this case, I got a different error:

R3 <- compileNimble(Rmodel)

mcmc.init <- configureMCMC(R3)

Error in distInfo$discrete : $ operator is invalid for atomic vectors

I searched for the message in this forum and other sources, but the only thing I got was the link to the source code in GitHub. Maybe I am missing lots of fundamentals in nimble. I hope I could get any helps to deal with the problem.

Thank you again!

2020년 12월 1일 화요일 오후 1시 5분 4초 UTC-8에 pdevalpine님이 작성:

[Perry de Valpine]

Typically the default of returnDef = FALSE is what one wants.  returnDef = TRUE requests that a model definition but not a model itself be returned.

The distInfo$discrete error is not immediately obvious.  Would you be willing to share reproducible code, which you might prefer to do off-list since it would include your data?  That would be the fastest way for one of us to look at what's happening.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/21f75545-b2dc-47b6-b526-e81daa86f663n%40googlegroups.com.

[Insang Song]

Dear Perry:

Thank you for the response again.

I attached the files and the code. I hope they work well!

[Perry de Valpine]

Hi Insang,

The problem appears to be that you declare elements of delta_1it twice, first as vector nodes following a CAR distribution and next as scalar nodes in an autoregressive structure.  This won't work.  Nodes should be declared only once.  This case got past error-trapping for duplicate declarations, so that's something we'll look at improving. 

I'll offer a few additional comments that might be helpful:

You can use calculate = FALSE in nimbleModel to save time, of course without the benefit of initialized calculations, but MCMC will initialize calculations later anyway (in compiled code).

We typically pass the uncompiled model (rather than the compiled model) to configureMCMC, but either should be ok.

nimbleMCMC will call nimbleModel and compileNimble for you.  You might realize that and were breaking down the steps in a debugging effort.  One advantage of going through the handful of steps yourself is to gain control such as using calculate=FALSE in nimbleModel.

I hope you are able to get up and running now.

-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/83bf9c4d-cefe-4dc9-b0ec-e5938c2b9563n%40googlegroups.com.

[Insang Song]

Sorry for the late response.

I tested several cases following your instruction (and it took 10+ hours per trial). MCMC ran properly, but at last all models crashed. I changed the number of chains (4->3->2), but none of them gave me any result. It is astonishing to see the ~1MB data resulted in 100+GB of memory (including paged files, I have a machine with 64GB RAM) after running MCMC. I think the MCMC crashed at last because of the shortage of memory. Could you let me know there is a rule of thumb to estimate how much memory would be consumed? Then I can prepare other resources (e.g., AWS, Azure, or adding more RAMs in my machine) to run my model.

Meanwhile, I saw a message before the MCMC began to run: "The logProb of y[3108, 14] is NA or NaN" sort of thing. I am not sure the message means, as the input data has no NA or NaN values (even zeros).

Thank you.

2020년 12월 1일 화요일 오후 6시 15분 43초 UTC-8에 pdevalpine님이 작성:

[Perry de Valpine]

Hi Insang,

There are various reasons memory use can explode for models like this.  Let me refer you to two previous threads for some of the ideas on how you can manage this: https://groups.google.com/u/2/g/nimble-users/c/8vaTHQoJz2M, https://groups.google.com/u/2/g/nimble-users/c/920URfnFiTY. 

It is also important to realize that you are monitoring (in the example code you sent previously) about 87000 nodes.  That will take about 2/3 Mb per stored iteration, or about 3.25 Gb per 5000 iterations, per chain (the number that was in your example code).  That doesn't account for all the memory use, but it's worth noting.

The NA or NaN messages are saying that the log probability calculation for a node is not valid.  That usually means not that the input data for the node is invalid but that some initial values for its parent nodes are giving an invalid probability.  This may or may not be a problem. The MCMC does try to initialize everything.  A good way to diagnose issues is to run the MCMC for 0 iterations (which will trigger the initialization) and then inspect the compiled model object.  For example, if it is called Cmodel, you can do:

Cmodel$calculate("y[5, 10]") # perhaps this is NA

# Then inspect parent nodes to see how this is happening:

Cmodel$beta1[5, 10] # etc.

I hope that helps.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/4a12254d-6f5a-4953-8c9b-1ccdba467f57n%40googlegroups.com.