Error in runjags when fitting a multinomial model: "is partly observed and partly missing"

Question

I am new to Jags and I'm trying to fit a multinomial model to my data. When I run the code I get the following error: "positive.counts[1,1:9] is partly observed and partly missing". I googled it and I found that this is due to the fact that a node cannot have observed and missing values at the same time. This is because in my data (see code below) there are values and NA in the same row. If I substitute the NA with 0 values the model works properly. Does anyone have a solution to this? Below you can find the data and the code!

Many thanks in advance,

Elisa

##########################################################################
# load jags
library(runjags)

# define the data: 
data <- list("N_cases"=c(978,737, 737, 1189,  270,  268), "positive.counts" = matrix(c(649 ,567 ,464 ,821,  98, 117,203 , 133,  81, 290,  41,  26,3, 7, 4, 6, 5, 0,NA, NA, NA, NA, 20, 19,24, 15,  3, NA, 21, 15,NA,  NA, 184,  NA,  17,  23, NA, NA, NA, NA, 26, 17,99, 15,  1, 72, 14, 25,NA, NA, NA, NA, 28, 26),6,9), "n"= 6,"n_responses" = 9)

# define the model
model {   mu.w <- 0
          prec <- 0.5
          for (s in 1:n_responses) {
            w[s] ~ dnorm(mu.w, prec);
            a[s] <- exp(w[s]); # positive parameter
          }
          for (i in 1:n){
            positive.counts[i,1:n_responses] ~ dmulti(p[i,1:n_responses], N_cases[i])
          }  
            for(i in 1:n){  
            for (s in 1:n_responses) {
              delta[i, s] ~ dgamma(a[s], 1)
            }
          }
          for(i in 1:n){  
            for (s in 1:n_responses) {
              p[i, s] <- delta[i,s] / sum(delta[i,1:n_responses])
            }
          }
        }  

# run the model
n.adapt=1000
n.burn=5000
n.iters=10000
n.chains <- 5;
n.total.samples <- 10000
n.samples.per.chain <- (n.total.samples %/% n.chains) 
n.thin <- n.iters %/% n.samples.per.chain; if(n.thin==0) n.thin <- 1;

tomonitor <- c("a","p")

mcmc.post <- run.jags(model="multi_model.jags",
                                  data=data,
                                   method="parallel", 
                                  sample=n.samples.per.chain, 
                                  burnin=n.burn,
                                  adapt=n.adapt,
                                  n.chains=n.chains,
                                  thin=n.thin,
                                  monitor= tomonitor);

Answer 1

As you've discovered, JAGS doesn't let you model a partially observed multinomial distribution. That's a limitation. If you replace the NAs with zeroes, you're no longer saying "I have no data for the combination of row and column variables" but rather "I have observed no events for this combination of row and column variables". So the model will (probably) run and produced some output. But you're no longer modelling your observed dataset. So the model's results don't apply to your data any more. Conclusions you draw from it will be invalid.

So, you need a way of modelling multinomial data with incomplete observations. One way of doing that is to transform your data into a series of related binomial variables. For each category c[i] in your multinomial data with k categories, you create a binomial variable Y[i] ~ Bin(p[i], n) where n is the number of observations (N_cases in your data?) with that combination of predictors (rows of positive.counts in your data?) and p[i] = Prob(x >= c[i]). In other words, p[i] is the probability that an observation with the relevant set of predictor variables is in category c[i] or above. (So, by definition, p[1] = 1.)

I'm sorry if that explanation is a little abstract, but you've given us your code but no context. I can't explain it in terms of your actual model, because I don't know what it is.

Two ways of fitting such a model in JAGS are given in this post. You can trust the author of the answer. He is Martyn Plummer. He wrote JAGS.