Many rows of NA values in MCMC output

[Athul Sudheesh]

In what circumstances does one get MCMC output like this (this is just a subset of the mcmc output where  Rhat where not < 1.1):

I have seen Rhat being NA values before, but this is the first time I see entire rows of MCMC output being NAs. I know probably I have a wrong model, but I am interested in knowing what might be causing this behavior. 

[Perry de Valpine]

Dear Athul,

NAs in the output typically mean something was not initialized in a valid way or the model has a problem.

Below is a version of your code (from your previous message) that works. Please check if I understood what you wanted to express in model code.

I could not follow entirely your intention from your previous message. Lists can not be used in nimble models.

I also could not follow entirely your intended model, because you had an "m_A[k]" inside of a for-loop without k as an index. No worries that you were not sure how to express your model, but I'm just saying please check below to see if I made good guesses. I think you want eff_ability[i,j] to be the minimum over k of m_A[k] * skills[i, k], considering only k such that Q[j, k] == 1. That is what I have implemented below to try to give you something that works to build from. I hope that helps get you launched.

This looks related to Item Response Theory models. In case it is helpful, we have a paper on IRT models in nimble here: https://doi.org/10.3102/10769986221136105

Also in case it is helpful, we have a user manual at r-nimble.org.

Best wishes,

Perry

library(CDM)
library(nimble)
data( data.fraction1 )
X = data.fraction1$data
Q = data.fraction1$q.matrix



fullModel <- nimbleCode({
  #coeff. for the skills
  for(k in 1:K){
    m_A[k] ~ dbeta(1,1)
  }


  # Prior for item difficulty
  for(j in 1:J){
    D[j] ~ dbeta(1,1)
  }

  for (i in 1:N){


    for(k in 1:K){
      lambda[i,k] ~ dbeta(1,1)
      skills[i,k] ~ dbern(lambda[i,k])
    }

  }


  for(i in 1:N){
    for(j in 1:J){

      # Modeling item response probability ========
      # Logic to compute effective ability ========
#      q_j[1:K] <- which(Q[j,1:K] == 1)
      eff_ability[i,j] <- indexed_min(m_A[1:K], skills[i,1:K], Q[j, 1:K])
      # =========================================
      theta[i,j] <- eff_ability[i,j] - D[j]
      pi[i,j] <- ilogit(1.7*theta[i,j])

      # Item-Response Likelihood
      X[i,j] ~ dbern(pi[i,j])
      # ===============================================
    }
  }
})

indexed_min <- nimbleFunction(
  run = function(m = double(1), skills = double(1), Q = double(1)) {
    returnType(double())
    K <- length(m)
    result <- Inf
    for(k in 1:K) {
      if(Q[k] == 1) {
        result <- min(result, m[k] * skills[k])
      }
    }
    return(result)
  }
)

modelConstants <- list(
  N = dim(X)[1],
  J = dim(X)[2],
  K = dim(Q)[2]
)
data <- list(
  X = X,
  Q = Q
)
initial_values <- list() # your value of NULL generates a warning

M0 <- nimbleModel(
  code = fullModel,
  data = data,
  inits = initial_values,
  constants = modelConstants,
  name = "baseline model"
)


dataNodes <- M0$getNodeNames(dataOnly = TRUE)
parentNodes <- M0$getParents(dataNodes, stochOnly = TRUE)
simNodes <- M0$getDependencies(parentNodes, self = FALSE)

M0$calculate() #NA
# This may simply mean the model is not fully initialized.
# It is better to provide a full set of initial values.
# If you don't want to, the MCMC will initialize from the priors.
# I'll build the mcmc manually to see if it works
mcmc <- buildMCMC(M0, monitors = params.to.save)
cM0 <- compileNimble(M0)
cmcmc <- compileNimble(mcmc, project = M0)

samples <- runMCMC(cmcmc, niter = 100, nburnin = 0)
dim(samples)
sum(is.na(samples)) # 0
# It looks superficially like it worked.


# I will reduce these run sizes by 10x for a quicker and smaller check.
n.iter <- 1500
n.burnun <- 200
n.chain <- 1
params.to.save <- parentNodes

mcmc_fit<- nimbleMCMC( model = M0,
                      monitors = params.to.save,
                      niter = n.iter,
                      nburnin = n.burnun,
                      nchains = n.chain,
                      samplesAsCodaMCMC=TRUE, WAIC=TRUE)

summary(mcmc_fit[[1]])

# It looks superficially like it worked.

--
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/76189bbf-5eb9-4c18-bc7e-c349058a1c3dn%40googlegroups.com.

[Athul Sudheesh]

Thank you so much for your solution Dr. Valpine. What you described for eff_ability[i,j] was what I was trying to achieve.

Sorry, you are correct it should have been m_A[1:K] instead of m_A[k]. I am really sorry for the confusion. 

Once again, THANK YOU SO MUCH for taking the time and effort from your busy schedule to help me fix my NIMBLE model! :)