Simulation from a model using nimble

[Bill Chen]

Hi Nimble expert,

I am trying to use nimble to simulate 100 datasets based on different simulated variables. However, when I replace different variables, I always get the same simulated dataset, I guess it might be because I did not change the initials...I have attached the code and could anyone please let me know if I messed something up? Thanks in advance!

simCode  <- nimbleCode(
  {
    for (i in 1:N){
     
      O[i] ~ dpois(mu[i])  
      log(mu[i]) <- log(expected[i])+alpha+ice[i]*betai[i]+u[i]+v[i]+b0*ice[i]
       
     
     
        betai[i]<-inprod(beta_sm[], X_sm[i,])+inprod(beta_pm[],X_pm[i,])+
        inprod(beta_un[],X_un[i,]) #this is from spline function
     
        v[i] ~ dnorm(0, tau = tau.v)  
     
    }
   
    u[1:N]~dcar_normal(adj[1:L], weights[1:L], num[1:N], tau1, zero_mean = 1)
   
   
    # Priors:
    alpha ~dnorm(0,sd=100)                                  
   
    b0~dnorm(0,sd=100)

    tau1~dgamma(1,0.01)
    tau.v~dgamma(1,0.01)
   
    for(j in 1:7){
      beta_sm[j] ~ dnorm(0, sd=100)
      beta_pm[j] ~ dnorm(0, sd=100)
      beta_un[j] ~ dnorm(0, sd=100)

    }
   
 
  }
)



**each simulated dataset should be based on different values of X_sm, X_pm and X_un.

    bsMat.s, bsMat.p and bs.Mat.u are the predefined matrix**

simModel <- nimbleModel(code = simCode,

                       data=list(X_sm=bsMat.s, X_pm=bsMat.p, X_un=bsMat.u),

                        constants = list(expected =n,ice=mergeus$ICEwbinc,
                                         N = Nareas, L=length(adj), adj=adj, weights=weights,
                                         num=num),
                        inits = list(alpha=0.01, betai=rep(0.1,3108),beta_sm=rep(0,7),beta_pm=rep(0,7),
 beta_un=rep(0,7),u=rep(1,3108),v= rep(1,3108),tau1=3, tau.v=1, b0=1))
       
simModel$calculate('mu')


nodesToSim <- simModel$getDependencies(c("alpha", "betai","beta_sm","beta_pm","beta_un","u","v",
                                         "tau1","tau.v","b0"),self = F, downstream = T)

simModel$simulate(nodesToSim,includeData = FALSE)

simModel$mu

I tried to replace different values for X_sm, X_pm, X_un, but always get the same mu.

Many thanks!

Bill

[Perry de Valpine]

Hi Bill,

The example isn't fully reproducible (I don't have your constants or data) so I'm just going on code reading.

I'm wondering if you want "self=TRUE". It looks like the nodes argument to simModel$getDependences includes all of the stochastic nodes. But if you put self=FALSE, then they will be omitted from nodesToSim. (Take a look at nodesToSim to confirm). Without any new stochastic draws, I don't think you'll get different values of mu, since mu is deterministic given some of the random variables, and none of those are included in nodesToSim. If that's not the right answer, please provide a fully reproducible example. HTH.

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/CAJsVCKQiYeHc-o3eEhkHYw_SMmL_a5_Ama-sJpTsR_BubAWO7Q%40mail.gmail.com.

[Bill Chen]

Hi Perry,

Thanks for the tips! However, I still could not get the simulated datasets...Please see the attached document. I appreciate any suggestions.

Best,

Bill

[Perry de Valpine]

Hi Bill,

Thanks for sharing the full setup. I was able to modify the file names according to what you provided. I commented out the u part of the model because I didn't seem to have weights, adj or num set up from your code.

It looks like it works for me. With the numbers you have, the values of mu can be enormous, so I look at log(mu).

old_log_mu <- log(simModel$mu)

simModel$simulate(nodesToSim,includeData = FALSE)

new_log_mu <- log(simModel$mu)

plot(old_log_mu, new_log_mu) # different numbers.

I didn't move on to the second setup in your code because the first one seemed to work. Does this solve your problem, or can you clarify what you're seeing as not working?

Cheers,

Perry