NA in node

[cwliu007]

Hi jagsUI developer,

Say I have a monitored parameter b which is a 2*3 matrix like c(1, 2, 3 ,4, NA, NA) I found the JAGS can finish the estimation, but jagsUI failed to compute the summary statistics for b. A workaround to avoid this issue is to give zeros to b[2,2] and b[2,3] in the BUGS syntax. Anyway, I think jagsUI should be able to deal with this situation by default.

Thank you.

[slozan0]

I think I ran into the same problem. For my particular case there is no workaround because the data is censored, and JAGS uses NA to process censored data.

model { for (i in 1:N) { y[i] ~ dinterval(x[i], lim[i,]) x[i] ~ dweib(a, b) } a ~ dunif(0,10) b ~ dunif(0,10) }

Thanks in advance,

[kenkellner]

Hi @genwei007 and @slozan0, sorry for my slow response on this. I'm getting started in a new job and plan to tackle these outstanding issues as soon as I get my feet under me. This one's at the top of my list.

[slozan0]

Thank you @kenkellner and congratulations. Would posting a complete example help you? JAGS sets this model in a peculiar way.

[kenkellner]

@slozan0 I think I have a working fix but if you could post a (hopefully small) complete example for me to test that would indeed be helpful.

[slozan0]

@kenkellner, here is a working example. Thank you for your help.

Set the model.

modelString <- "
model
{
  for (i in 1:n) {
    y[i] ~ dinterval(t[i], c[i, ])
    t[i] ~ dweib(v, lambda)
  }
  v      ~ dunif(0, 10)
  lambda ~ dunif(0, 10)
}
writeLines(modelString, con = "model.txt")

The data.

time,time2,event
5,  15, 3
7.5,15, 3
8,  15, 3
9,  15, 3
8,  10, 3
8,  10, 3
8,  10, 3

Fitting the data.

# read time to event
mySurvDf <- read.csv('temp1.csv', header = T, stringsAsFactors = FALSE)

# following JAGS notation
# y ~ dinterval(t, c[1:M])
# where dinterval will return a likelihood P( y| t, c ) for the parameters
# c and t that takes the value 1 when t lies within an interval defined by
# the cut points c1 .. cM and 0 otherwise

c <- cbind(mySurvDf$time, mySurvDf$time2)  # time intervals for each data point
n <- nrow(c)
t <- rep(NA, n) # everything is censored in the bioassay therefore NA for st
y <- rep(1, n)  # everything is censored in the bioassay therefore 1 for indicator variable

dataList <- list(
  n = n,
  t = t,
  c = c,
  y = y)
#Run the chains

parameters         <- c("v", "lambda")
adaptSteps         <- 1000         # Number of steps to "tune" the samplers.
burnInSteps        <- 1000        # Number of steps to "burn-in" the samplers.
numSavedSteps  <- 1000       # Total number of steps in chains to save.
nChains               <- 3               # Number of chains to run.
thinSteps             <- 5               # Number of steps to "thin" (1=keep every step).

# Create, initialize, and adapt the model:
jagsModel <- jags.model( "model.txt", data = dataList,
                         n.chains = nChains, n.adapt=adaptSteps )
# Burn-in:
cat( "Burning in the MCMC chain...\n" )
update( jagsModel, n.iter = burnInSteps )
# The saved MCMC chain:
cat( "Sampling final MCMC chain...\n" )
mcmcCoda <- coda.samples( jagsModel, variable.names = parameters,
                         n.iter = nPerChain, thin = thinSteps )

# EXAMINE THE RESULTS
# Convert coda-object codaSamples to matrix object for easier handling.
mcmcChain <- as.matrix(mcmcCoda)
chainLength <- nrow(mcmcChain)

median(mcmcChain[ , "lambda"])
median(mcmcChain[ , "v"])

The resulting values should be around

> median(mcmcChain[ , "lambda"])
[1] 0.003485706
> median(mcmcChain[ , "v"])
[1] 2.510426

[kenkellner]

Hi @slozan0:

Thanks. Using your example, I made a few modifications/additions:

#Can leave t out because it's all NAs
dataList <- list(
  n = n,
  c = c,
  y = y)

library(jagsUI)

out <- jags(data = dataList,
                 parameters.to.save = parameters,
                 model.file ="model.txt" ,
                 n.chains = nChains,
                 n.adapt = adaptSteps,
                 n.iter = numSavedSteps+burnInSteps,
                 n.burnin = burnInSteps,
                 n.thin = thinSteps)

Got this as output from jagsUI:

print(out,digits=5)

JAGS output for model 'model.txt', generated by jagsUI.
Estimates based on 3 chains of 2000 iterations,
adaptation = 1000 iterations (sufficient),
burn-in = 1000 iterations and thin rate = 5,
yielding 600 total samples from the joint posterior.
MCMC ran for 0.003 minutes at time 2018-01-17 14:33:23.

            mean      sd    2.5%     50%   97.5% overlap0 f    Rhat n.eff
v        2.85536 0.88507 1.47659 2.69840 5.03008    FALSE 1 1.09430    29
lambda   0.00563 0.01041 0.00001 0.00205 0.03470    FALSE 1 1.30988    23
deviance 0.00000 0.00000 0.00000 0.00000 0.00000    FALSE 1     NaN     1

Which is close to your results. It seems like c8c643f fixes the problem. I'm working on an adjustment to jags() which will dump raw rjags output when processing output fails, instead of just crashing. Once that's done I'll package this up and submit it.