Using different design matrix for different time periods

[Andrew Greene]

Hi all,

I am trying to define a model that uses all available covariates at each time period. For example in time period t=1, there are two available and complete variables, x1 and x2. In time period t=2, x1, x2, and x3 are available.

I am trying to use the following model to accomplish this (only the highlighted bits are relevant to the issue):

model_code=nimbleCode({
  for(t in 1:T){ # Loop through time periods
    for(s in 1:S){ # Loop through different areas
      y[((t-1)*S+s)] ~ dZIP(lambda[((t-1)*S+s)],zeroProb=p[s])
     
      if (t==1) {
        log(lambda[((t-1)*S+s)]) <- off[((t-1)*S+s)] + StInt[states[s]] + mu[t] + inprod(X1[((t-1)*S+s),],beta1[]) + u[s]
      } else {
        log(lambda[((t-1)*S+s)]) <- off[((t-1)*S+s)] + StInt[states[s]] + mu[t] + inprod(X2[((t-1)*S+s),],beta2[]) + u[s]
      }
    }
   
    mu[t] ~ dflat()
   
  }
  #### define spatial random effect with exponential cov function
  cov[1:S, 1:S] <- expcov(dists[1:S, 1:S], rho, sigma)
  u[1:S] ~ dmnorm(zeros[1:S], cov = cov[1:S, 1:S])
 
  #### define prior distributions
  for(i in 1:bp1){
    beta1[i] ~ dflat()
  }
  for(i in 1:bp2){
    beta2[i] ~ dflat()
  }
  for(i in 1:numStates) {
    StInt[i] ~ dflat()
  }
  rho ~ dunif(0,5)
  sigma ~ dunif(0, 100)
  for (s in 1:S) {
    p[s] ~ dunif(0,1)
  }
})


constants <- list(S = n, T=T, dists = dist.mat, zeros = rep(0, n),
                  X1 = X1, X2 = X2, off=log(Pop),bp1=dim(X1)[2], bp2=dim(X2)[2],
                  states = States, numStates = max(States))
data <- list(y = VFdata$cases)
inits <- list(beta1 = rep(0,ncol(X1)), beta2 = rep(0,ncol(X2)), StInt = rep(0,max(States)), sigma = 1, rho = 0.2,mu=rep(-10,T),p=rep(.5, n))

set.seed(1)

## setup initial spatially-correlated latent process values
inits$cov <- cExpcov(dist.mat, inits$rho, inits$sigma)
inits$u <-  as.vector(t(chol(inits$cov)) %*% rnorm(n))
#inits$u <- inits$s[ , 1]  # so can give nimble a vector rather than one-column matrix
inits$y <- y
inits$y[which(is.na(y)==TRUE)]=0

model <- nimbleModel(model_code, constants = constants, data = data, inits = inits)

I think this approach makes sense, but when I try to define the nimble model, an error is thrown:

Defining model

Error in t == 1 :

   comparison (1) is possible only for atomic and list types

Error in codeProcessIfThenElse(code[[i]], constants, envir) :

   Cannot evaluate condition of 'if' statement: t == 1.

Condition must be able to be evaluated based on values in 'constants'.

I understand that if-else statements are only used to create model variants at the time of definition, but is there any way to do what I am trying to accomplish in nimble?

Thank you

[Daniel Turek]

Andrew, I can't run your model without the relevant functions and data, but replacing the highlighted part of the model (the if(t == 1) ... else ...) with the following:

log(lambda[((t-1)*S+s)]) <-
    equals(t,1) * (off[((t-1)*S+s)] + StInt[states[s]] + mu[t] + inprod(X1[((t-1)*S+s),],beta1[]) + u[s]) +
    (1 - equals(t,1)) * (log(lambda[((t-1)*S+s)]) <- off[((t-1)*S+s)] + StInt[states[s]] + mu[t] + inprod(X2[((t-1)*S+s),],beta2[]) + u[s])

should accomplish the same.  Please double-check the code, but it's just making basic use of the equals(t,1) function.

Another option would be pushing this code into a custom-written (deterministic) nimbleFunction, but if the approach above suffices, then there's no need to do that.

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/c6cbb6d3-ba87-4b59-9940-a09ba146db3dn%40googlegroups.com.

[Chris Paciorek]

Just to be explicit, one cannot use 'if' statements as part of the code of a single model. "if" can only be used to define multiple models from the same model code, and the "if" is stripped out when defining a given model from the original code.

Hence Daniel's suggested syntax to avoid the use of "if".

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