Dear BUGS users,
I am interested in whether anyone has coded the product binomial
distribution as described in McCullagh & Nelder (1989, pg 338). It is
sometimes called the joint binomial distribution, although M&N didn't
give it a name.
I have written code for the distribution since it is not included as one
of the standard BUGS distribution. Below is a description of the
distribution and the BUGS code that Ive written for a particular model,
which uses the distribution. Can someone point out for me where I have
gone wrong?
The distribution :
Suppose we have G 2x2 tables, where for the g'th table, the (1,1)'th
cell value is distributed so that N_{11g} ~ Binomial (N_{1 \bullet g},
\pi_{1g}) and similarly, N_{21g} ~ Binomial(N_{2 \bullet g}, \pi_{2g}).
Here N_{1 \bullet g} and N_{2 \bullet g} are the marginal frequencies of
the two row categories, and \pi_{1g} and \pi_{2g} are the parameters of
interest. It ends up that the likelihood function is written as a
summation so, where the product of the two binomials is summed between
the Frechet bounds of N_{11g}.
What we want to do is model \pi_{1g} and \pi_{2g} in terms of the single
parameter \theta_g and the marginal frequencies, where \theta_g is given
a normal prior with known mean and variance.
The BUGS code (for five groups) :
==================
Model {
########################################################################
# * Estimates the values of \pi_{1g} and \pi_{2g}, which are parameters
of the product #
# binomial distribution of McCullagh & Nelder (1989, pg
338). #
# * The paremeters \pi_{1g} and \pi_{2g} are parameterised in terms of
theta_g and the #
# known marginal proportions of G groups of 2x2 tables.
#
# * theta_g is assumed to be normally distributed, centred at 1 and
with a variance #
# depending on the marginal proportions #
########################################################################
for (g in 1:G){
X1[g] <- NX[g]/N[g] # X1=Marginal proportion for Row 1
Y1[g] <- NY[g]/N[g] # Y1 = Marginal proportion for Column 1
a[g] <- max(0, NY[g]-(N[g]-NX[g])) # Lower Frechet Bound of
N_{11g}
b[g] <- min(NX[g], NY[g]) # Upper Frechet Bound for N_{11g}
d[g] <- b[g] - a[g] + 1
}
for (g in 1:G){
ones[g] <-1
ones[g] ~ dbern(f [g])
tildeX[g] <- X1[g]/(1-X1[g])
tildeY[g] <- Y1[g]/(1-Y1[g])
########################################################
# The parameter theta_g is generated so that it lies within these
limits #
########################################################
lower[g] <- max(0, 1/Y1[g] + 1/X1[g] - 1/(X1[g]*Y1[g]) )
upper[g] <- min( 1/Y1[g], 1/X1[g] )
vartheta[g] <- 1/(tildeX[g]*tildeY[g]) # Variance of theta[g]
dependent on marginal proportions
######################################################################
# theta_g is a parameter whose prior distribution is normal centred
about 1, with variance #
# depending on the marginal proportions. This parameter is generated
to be censored #
# between two known quantities (above), so that estimates \pi_{1g}
and $\pi_{2g} #
# are valid #
######################################################################
theta[g] ~ dnorm(1, vartheta[g]) I(lower[g], upper[g])
####################################################################
# The parameterisation of pi_{1g} and \pi_{2g} so that
#
# \pi_{1g} = f _1(theta_g, X1_g, Y1_g) and \pi_{2g} = f_2 (theta_g,
X1_g, Y1_g) #
####################################################################
pi1[g] <- Y1[g] * theta[g]
pi2[g] <- (Y1[g] - pi1[g] * X1[g])/(1- X1[g])
#######################################
# The Product, or Joint, Binomial Distribution #
#######################################
for (j in 1 : d[g]){
X1C[g, j] <- exp(logfact(NX[g]))/(exp(logfact(a[g] + j - 1)) *
exp(logfact(NX[g]-(a[g] + j - 1))))
X0C[g, j] <- exp(logfact(N[g]-NX[g]))/(exp(logfact(NY[g]-(a[g] + j -
1))) *
exp(logfact(N[g] - NX[g] - NY[g] +(a[g] + j - 1))))
ff[g, j] <- X1C[g, j] * X0C[g, j] * pow(pi1[g], a[g] + j - 1)*
pow(1-pi1[g], NX[g] - (a[g] + j - 1)) * pow(pi2[g], NY[g] -
(a[g] + j - 1)) * pow(1 - pi2[g],
N[g] - NX[g] -NY[g] + (a[g] + j - 1)) }
f[g] <- ff[g, ]
}
}
DATA
list(N=c(707, 408, 668, 293, 583), NX=c(584, 273, 451, 220, 445),
NY=c(173, 99, 177, 79, 171), G=5)
INITIAL VALUES
list( theta = c(0, 0, 0, 0, 0), pi1 = c(0.5, 0.5, 0.5, 0.5, 0.5), pi2 =
c(0.5, 0.5, 0.5, 0.5, 0.5))
==========
BUGS tells me there is a problem with d[g] in the for loop, as the
distribution is being defined. I can understand this as X1C[g, d[g]]
doesn't make sense. Is this how BUGS interprets the problem, and how
does BUGS deal with this situation ??
Hope someone can help.
Cheers
Eric
Dr Eric J. Beh,
Research Fellow,
School of Mathematics and Applied Statistics,
University of Wollongong,
Northfields Avenue,
Wollongong, NSW, 2522,
Australia
Email : [log in to unmask]
Web: http://www.uow.edu.au/informatics/maths/staff/beh.html
Phone : (+61) 2 42215662 Fax : (+61) 2 42214845
-------------------------------------------------------------------
To mail the BUGS list, mail to [log in to unmask]
You can search old messages at www.jiscmail.ac.uk/lists/bugs.html
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
This list is for discussion of modelling issues and the BUGS
software. For help with crashes and error messages, first mail
[log in to unmask]
