How to obtain multinomial probabilities in WinBUGS with multiple regression

Question

In WinBUGS, I am specifying a model with a multinomial likelihood function, and I need to make sure that the multinomial probabilities are all between 0 and 1 and sum to 1.

Here is the part of the code specifying the likelihood:

e[k,i,1:9] ~ dmulti(P[k,i,1:9],n[i,k]) 

Here, the array P[] specifies the probabilities for the multinomial distribution.

These probabilities are to be estimated from my data (the matrix e[]) using multiple linear regressions on a series of fixed and random effects. For instance, here is the multiple linear regression used to predict one of the elements of P[]:

P[k,1,2] <- intercept[1,2]  +  Slope1[1,2]*Covariate1[k]  +
Slope2[1,2]*Covariate2[k]  +  Slope3[1,2]*Covariate3[k]  
+ Slope4[1,2]*Covariate4[k] +  RandomEffect1[group[k]]  +  
RandomEffect2[k]

At compiling, the model produces an error:

elements of proportion vector of multinomial e[1,1,1] must be between zero and one

If I understand this correctly, this means that the elements of the vector P[k,i,1:9] (the probability vector in the multinomial likelihood function above) may be very large (or small) numbers. In reality, they all need to be between 0 and 1, and sum to 1.

I am new to WinBUGS, but from reading around it seems that somehow using a beta regression rather than multiple linear regressions might be the way forward. However, although this would allow each element to be between 0 and 1, it doesn't seem to get to the heart of the problem, which is that all the elements of P[k,i,1:9] must be positive and sum to 1.

It may be that the response variable can very simply be transformed to be a proportion. I have tried this by trying to divide each element by the sum of P[k,i,1:9], but so far no success.

Any tips would be very gratefully appreciated!

(I have supplied the problematic sections of the model; the whole thing is fairly long.)

Answer 1

The usual way to do this would be to use the multinomial equivalent of a logit link to constrain the transformed probabilities to the interval (0,1). For example (for a single predictor but it is the same principle for as many predictors as you need):

Response[i, 1:Categories] ~ dmulti(prob[i, 1:Categories], Trials[i])

phi[i,1] <- 1
prob[i,1] <- 1 / sum(phi[i, 1:Categories])
for(c in 2:Categories){
    log(phi[i,c]) <- intercept[c] + slope1[c] * Covariate1[i]
    prob[i,c] <- phi[i,c] / sum(phi[i, 1:Categories])
}

For identifibility the value of phi[1] is set to 1, but the other values of intercept and slope1 are estimated independently. When the number of Categories is equal to 2, this collapses to the usual logistic regression but coded for a multinomial response:

log(phi[i,2]) <- intercept[2] + slope1[2] * Covariate1[i]
prob[i,2] <- phi[i, 2] / (1 + phi[i, 2])
prob[i,1] <- 1 / (1 + phi[i, 2])

ie:

logit(prob[i,2]) <- intercept[2] + slope1[2] * Covariate1[i]
prob[i,1] <- 1 - prob[i,2]

In this model I have indexed slope1 by the category, meaning that each level of the outcome has an independent relationship with the predictor. If you have an ordinal response and want to assume that the odds ratio associated with the covariate is consistent between successive levels of the response, then you can drop the index on slope1 (and reformulate the code slightly so that phi is cumulative) to get a proportional odds logistic regression (POLR).

---

Edit

Here is a link to some example code covering logistic regression, multinomial regression and POLR from a course I teach:

http://runjags.sourceforge.net/examples/squirrels.R

Note that it uses JAGS (rather than WinBUGS) but as far as I know there are no differences in model syntax for these types of models. If you want to quickly get started with runjags & JAGS from a WinBUGS background then you could follow this vignette:

http://runjags.sourceforge.net/quickjags.html