- Fitting Normal mixture models
It is straightforward to use different variances in a normal mixture model:
for example, in the eyes example
{
for (i in 1:N){
y[i] ~ dnorm(lambda[T[i]],tau[T[i]]);
T[i] ~ dcat(P[])
}
lambda[1] ~ dnorm(0,1.0E-6);
lambda[2] <- lambda[1] + c;
c ~ dnorm(0,1.0E-6)I(0,);
}
There are various alternatives prior specifications for the precisions:
(1) tau[k] ~ dgamma(0.01,0.01);
(2) Make the variances similar, as Richardson and Green do in their
recent RSS paper. Then you could have
tau[k] ~ dgamma(a,b);
a and b are given appropriate priors.
(3)
log(tau[1]) <- (a+b)/2;
log(tau[2]) <- (a-b)/2;
appropriate priors are specified for a and b. b repesents the log ratio
of precisions and so a reasonable prior could be one such that
1/3 < exp(b) < 3. Similarly `a' represents the log product of the
precisions.
We have found that (3) gives the most stable results but BE VERY CAREFUL
WITH CONVERGENCE! These mixture problems are notoriously unstable.
Return to the top
- Mixture responses
Adapting the Hearts example, try the following
rather simpler construction that has a bivariate
structure for ep.
{
for (i in 1:I) {
group[i] ~ dcat(pn[]);
for (j in off[i]+1:off[i+1]) {
y[j] ~ dbin(ep[group[i],j],1);
ep[1,j] <- 0;
logit(ep[2,j]) <- a[i]+b.age*(age[j]-age.bar);
}
a[i] ~ dnorm(mu, tau);
}
}
Because of the line
logit(ep[2,j]) <- a[i]+b.age*(age[j]-age.bar);
a[i], and hence mu and tau, will only get likelihood
terms from y[j]'s that are currently in group 2. Otherwise
there will be no link between y[j] and a[i].
y[j] ~ dbin(ep[group[i],j],1);
The term ep[group[i],j] essentially acts as a `pick' function that
dynamically creates or removes the link between
y[j] and a[i].
Return to the top
- Fitting Bivariate Normals - an alternative to the Wishart distribution
In the `birats' example beta[i,] is assumed to have a bivariate normal
distribution:
for(i in 1:N){
beta[i,] ~ dmnorm(mu.beta[], Omega.beta[,]);}
where Wishart prior is then specified for the precision matrix
Omega.beta[,] ~ dwish(R[,],2);
This parameterisation is notoriously tricky - the following provides an
attractive alternative:
for(i in 1:N){
beta[i,1] ~ dnorm(mu.beta1,tau1);
beta[i,2] ~ dnorm(mu.beta2[i], tau2);
mu.beta2[i] <- lambda0 + lambda1*(beta[i,1] - mean(beta[i,1]));}
`Noninformative' priors can be specified for the lambda's and tau's.
This is equivalent to fitting a bivariate normal to beta[i,] with
mean (mu.beta1, lambda0 + lambda1*(mu.beta1 - mean(beta[i,1])), variance
(sigma1, sigma2 + sigma1*lambda1^2) and correlation
lambda1*sqrt(sigma1)/sqrt(sigma2 + sigma1*lambda1^2).
Return to the top
- An example of a hierarchical tobit model
This example of a hierarchical tobit model arises from a study in
which repeated measurements of aflatoxin exposure (a potential risk factor for
liver cancer) were taken, some of which fall below the detection limit
of the assay. I want to calculate the posterior expected mean level of
exposure for each subject.
This approach requires some manipulation of the data before using BUGS.
First a binary variable (is.below.cutoff) is created which indicates
whether the measurement (measured.afla) was below the detection limit.
Then the measured values which fall below the limit are set to missing.
In the BUGS model a bounded distribution is used for the measurements
with an upper bound equal to either the detection limit or a very large
(arbitrary) number, depending on whether the measurement fell below the
detection limit. The Gibbs sampler samples the observations we have set to
missing from the tail of the distribution.
Two minor points which might help you read this model
(1) Instead of log-transforming the data I'm using the lognormal distribution
(2) I'm using nested indexing because there is an unequal number of
observations on each subject. Subjects are identified by the variable "index".
model aflatoxin;
const
NOBS = 198, # number of observations
NSUB = 115, # number of subjects
CUTOFF = 5, # detection limit for assay
UPPERLIM = 10000; # very large upper limit for observations
var
mu, # Population mean on log scale
tau, # Inter-subject variation
sigma, # Intra-subject variation (measurement error)
mean.afla[NSUB], # Subject-specific mean
log.mean.afla[NSUB],
upper.lim[NOBS],
# -------- Data ----------
index[NOBS], # Subject identifier
measured.afla[NOBS], # Assay result
is.below.cutoff[NOBS];# Is assay result is below CUTOFF (1=yes,0=no)
data in "aflatoxin.dat";
inits in "aflatoxin.in";
{
mu ~ dnorm(0.0, 1.0E-4);
tau ~ dgamma(1.0E-3,1.0E-3);
sigma ~ dgamma(1.0E-3,1.0E-3);
for (i in 1:NSUB) {
log.mean.afla[i] ~ dnorm(mu, tau);
log(mean.afla[i]) <- log.mean.afla[i]
}
for (j in 1:NOBS) {
upper.lim[j] <- CUTOFF * is.below.cutoff[j]
+ UPPERLIM * (1 - is.below.cutoff[j]);
measured.afla[j] ~ dlnorm(log.mean.afla[index[j]], sigma)
I(,upper.lim[j]);
}
}
(Thanks to Martyn Plummer for this example.)
Return to the top
- Regression on the mean of a negative binomial
Suppose you have a set of counts, Y that you want to regress on X, according
to:
log(mu.Y[i]) <- a + bX;
where Y iself is ~ negative binomially with mean mu.Y and shape-parameter r.
Y[i] ~ dnegbin(p[i], r[i])
p[i] will need to be expressed as a function of r, a, b and
X[i]. This will presumably be a horrible non-linear expression,
so you will need to use Metropolis in Version 0.6 - making sure a and
b are bounded.
Return to the top
- How does BUGS deal with missing covariates?
As mentioned in Section 5.6 of the manual, BUGS has no automatic
way of dealing with missing data. If missing covariates are introduced
they must be given a prior distribution to fully specify
the probability model. The easiest thing is to give
each covariate a normal prior with unknown mean and
variance, and then the observed covariates can
provide information about these prior parameters.
The missing covariate is then estimated, automatically
taking into account the prior distribution and the observed
response. The cervix example of Examples Volume 2 shows this,
although with a discrete covariate.
Return to the top
- How can I use BUGS for simulation studies?
BUGS will easily simulate a single dataset doing a single update
and using `save', which then can be used
as a data file for a BUGS run. Alternatively, if multiple datasets
are wanted (ie many analyses) then provided the analysis is a closed
form expression then it can be put in as a logical function in BUGS
and will produce an answer at each iteration. However, if the analysis
itself requires iterative or simulation solution, there is no option but running BUGS many times. For helper code to do this using other software, see the Remote running page for WinBUGS v1.4 or consider using OpenBUGS resources.
Return to the top
- Is there a discrete uniform distribution available?
You need to construct your own distribution
as a categorical variable with a uniform prior
and which can take on the necessary integer
values. See the blockerht example in the first part
of the manual.
Return to the top
- How can I use a truncated sampling distribution?
BUGS will give the wrong answer if the I(,) construct is used,
since the parameters of a truncated distribution are also
present in the normalising constant. In theory, if this
constant can be written as an explicit function of the
parameters, BUGS could use the Metropolis algorithm
to learn about these parameters using the method of
arbitrary sampling distributions described below.
Another (and sometimes better) option is to implement the truncated distribution yourself using the WBDev tools. In particular, a truncated normal distribution is available as a simple add-in for WinBUGS.
Return to the top
- How can I use an arbitrary sampling distribution that isn't
in the current list that BUGS handles?
See the `Tricks' section in the manual.
Return to the top
- How can I generate a random position of a cell in a matrix?
Suppose the matrix is 3 rows by 4 columns, and you sample CELL using dcat and
get a number between 1 and 12, which you want to convert to a row and a
column. There are two solutions:
(1) do not use a matrix but use an array of size 12.
(2) set up two arrays
j[] = (1,1,1,2,2,2,3,3,3,4,4,4)
j[] = (1,1,1,1,2,2,2,2,3,3,3,3)
Then you can use
y[ i[cell],j[cell] ]
Return to the top
- How can I deal with interval censored data?
Interval censoring view can be seen as an extreme version of rounding, in which
one only knows the observation lies between two limits.
The methods used for classical rounding error
described in the manual (Section 9.6) should work OK.
Return to the top
- Zero denominators in a structured binomial
Consider the following part of a BUGS model:
for(i in 1:n) {
for(j in 1:(i-1)) {
logit(p[i,j]) <- beta[i]-beta[j];
N[i,j] <- X[i,j] + X[j,i];
X[i,j] ~ dbin(p[i,j],N[i,j]); (*)
}
}
For some (i,j), N[i,j] = X[i,j] + X[j,i] = 0 and this result
gives the "Invalid value for node expected positive value" error,
when statement (*) is "executed".
Is there some way of including statement (*) only under the
condition that N[i,j] is positive?
A clumsy way of dealing with this would be to place the X[i,j]'s with
non-zero N[i,j]'s into a single array indexed by k , thus leaving off
the tricky cases. Include as factors the row[k] and column[k] from
which they were taken, and then fit beta[row[k]] and beta[column[k]] in
the logistic model.
An alternative method of dealing with zero denominators (suggested by Michael Lavine, Duke University)
Return to the top
- How do I deal with data which is a random sum?
Generally speaking, you cannot give distributions to functions of random variables, so for example
model{
z <- x+y
x~dnorm(mu1,1)
y~dnorm(mu2,1)
}
#data
list(z <- 10)
}
is not possible in the BUGS language. However, if you can work out what the implied marginal distribution is for your data, you can use this. If the implied distribution is non-standard, or not available in closed form, you may need to use the WBDev tools, or the `Tricks' section of the manual.
Return to the top
- How do I program if...then statements?
These are not really part of the graphical model setup behind the BUGS language, and so are not generally available. However, simple if...then structures can be incorporated using step and equals commands, for example
y~dnorm(mu,tau)
mu <- mu0 + beta*(step(x)*x + (1-step(x))*pow(x,2))
fits a linear model for x>0 and a quadratic response for x<0
Return to the top
![[logo]](/bugs/images/bugslogo.gif)
