MRC Biostatistics UnitHow to use this page
Below are some frequently asked WinBUGS questions. If you don’t find an appropriate solution to your problem here, we advise that you
- Read the `health warning’ on our homepage and in the manual. MCMC can be less stable than other methods – this is not WinBUGS’ fault! – and there is no automatic protection against fitting inappropriate models.
- Read the manual – particularly the `Tips and Troubleshooting’ sections
- Work through some examples – understanding what’s happening in something similar to your code will help you identify coding and modelling problems
- New users: if possible, ask a local expert for advice. Remember WinBUGS is free; we try to help where possible but do not have time or resources to debug all users’ code `online’.
- Get in touch: See the list of contacts for the appropriate place to send your question
We apologise for not keeping these FAQs up-to-date.
FAQs regarding DIC can be found here.
Choose a topic
- Modelling Issues – generic tips for complex Bayesian modelling using the BUGS language
- WinBUGS Issues – strange things and bugs in WinBUGS. This section also contains information about implementation of WinBUGS on the various Windows systems that are available.
- Or search the archive of messages to the discussion list for past queries and solutions on a wide range of issues.
Modelling Issues: Questions
- Fitting Normal mixture models
- Mixture responses
- Fitting Bivariate Normals – an alternative to the Wishart distribution
- An example of a hierarchical tobit model
- Regression on the mean of a negative binomial
- How does BUGS deal with missing covariates?
- How can I use BUGS for simulation studies?
- Is there a discrete uniform distribution available?
- How can I use a truncated sampling distribution?
- How can I use an arbitrary sampling distribution that isn’t in the current list that BUGS handles?
- How can I generate a random position of a cell in a matrix?
- How can I deal with interval censored data?
- Zero denominators in a structured binomial
- How do deal with data which is a random sum?
- How do I program an if...then statement?
Modelling Issues: Questions and Answers
- Fitting Normal mixture modelsIt 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.
- Mixture responsesAdapting 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].
- Fitting Bivariate Normals – an alternative to the Wishart distributionIn 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[,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[,1])), variance (sigma1, sigma2 + sigma1*lambda1^2) and correlation lambda1*sqrt(sigma1)/sqrt(sigma2 + sigma1*lambda1^2).
- An example of a hierarchical tobit modelThis 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.)
- Regression on the mean of a negative binomialSuppose 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.
- 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] ]
- 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 binomialConsider 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.
- 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.
- 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
WinBUGS Issues: Questions
- A.1: What if WinBUGS crashes with a Trap full of gibberish?
- A.2: Why does WinBUGS just hang when updating?
- A.3: I want to have a vector of probabilities for a categorical variable, which is a row in a 4 dimensional array. WinBUGS says ‘multiple definitions of node’?
- A.4: What does the error message ‘logical expression too complex’ mean?
- B.2: I’ve drawn a new doodle – how do I add the data lists and other texts?
- B.3: Can WinBUGS run under Linux?
- C.1: What happens to virtual memory under Windows Pro2K?
- C.2: Can WinBUGS be run on a Mac?
- C.3: Is there a size restriction in WinBUGS?
WinBUGS Issues: Questions and Answers
A.1: What if WinBUGS crashes with a Trap full of gibberish?
The usual reason is that the sampling algorithm has failed due to the prior being incomparable with the data, or because therer is too little information concerning a parameter given a “vague” prior. Try to check what was happening just before the crash. Ensure initial values are appropriate, or use stronger priors. It is sometimes possible to carry on sampling by closing the Trap window and clicking on update a couple of times. However, if you repeatedly crash, it may be best to quit and start WinBUGS again.
A.2: Why does WinBUGS just hang when updating?
Two reasons have been found:
- There is an occasional problem with Windows NT – if the system is rebooted it goes away.
- WinBUGS may be struggling to update a variable. Check initial values are appropriate
A.3: I want to have a vector of probabilities for a categorical variable, which is a row in a 4 dimensional array. WinBUGS says ‘multiple definitions of node’?
WinBUGS complains about something like
a ~ dcat(p.a[1:2])
b ~ dcat(p.b[1:2])
c ~ dcat(p.c[1:2])
d ~ dcat(p.d[a,b,c,1:2])
There is a restriction in WinBUGS that is not in the manual. p.d[a,b,c,1:2] is a “mixture node”, in the sense that it identifies an appropriate vector to use as a parameter in a distribution. Currently mixture nodes are only allowed with one or two indices ie
p[i, 1:2] or p[i, j, 1:2]
say. There is a simple work round: for example write
d ~ dcat(p[k, 1:2])
k <- 8 * (a - 1) + 4 * (b - 1) + c
This `calculated index’ trick is useful in many circumstances.
A.4: What does the error message ‘logical expression too complex’ mean?
WinBUGS can only handle logical expressions of a certain length. Break the expression into two or more stages.
B.2: I’ve drawn a new doodle – how do I add the data lists and other texts?
You need to highlight the doodle, then either drag it into a new document, or paste it into the clipboard and copy into a new document. Then you should be able to add text.
B.3: Can WinBUGS run under Linux?
Running WinBUGS v1.4.x under Linux using Wine is described by Martyn Plummer on his web site. Others have reported successful and stable running under VMware.
Locally at the MRC Biostatistics Unit (home of WinBUGS) we are using CrossOver Office from Codeweavers to run WinBUGS under Linux.
Far more extensive resources for other versions of BUGS which operate under Linux are available on the OpenBUGS website.
C.1: What happens to virtual memory under Windows Pro2K?
When using extremely large models, your computer may resort to using virtual memory. Under Windows Pro2K, it appears that this will cause a crash. The problem does not seem to occur under Windows NT
C.2: Can WinBUGS be used on a Mac?
You can run WinBUGS on a Mac under emulators such as Virtual PC. Please note we have no plans to produce a dedicated Macintosh version
C.3: Is there a size restriction in WinBUGS?
Once you have downloaded the key the only restriction is the size of your computer.