Bayesian Version of Consensus Means Using the BUGS Program (Blaza Toman)

Introduction

The problem of determining a consensus mean based on data from multiple labs or methods has been fully addressed at NIST from the classical point of view. The result is encoded in the DATAPLOT command CONSENSUS MEAN. A Bayesian solution to this problem can be obtained through the application of hierarchical Bayes models via the Markov Chain Monte Carlo (MCMC) simulation software called BUGS. This report gives the BUGS code for two common hierarchical models and applies it to an example from SRM 1946, Lake Superior Fish Tissue.

Model 1

Multiple labs: Normal distributions with different means and variances

This model assumes that the multiple labs data comes from Normal distributions which have different means and different variances, that is that Y_ij is distributed as

and that the means delta(i)

and the variances sigma(i)

are related. That is,

and

where

tau ~ IGamma(0.0001,0.0001)

a ~ exp(1)

b ~ IGamma(0.0001,0.0001)

Borrowing and Vague Priors

This model allows for borrowing in the estimation of both the means delta(i)

, and the variances sigma(i)

. This means that even when some of the labs have extremely small sample sizes (even n=1), the lab variances can be estimated through the pooling of the hierarchical model. The prior distributions above are "vague", which is appropriate when real prior information in the form of expert opinion or prior data is not available. Analysis performed using vague prior distributions can be considered objective and is generally preferred by classical statisticians. When real prior information is available in the form of a mean and variance of

, it can be included by simply changing the mean and variance of the Normal distribution. If more robustness is required, a t-distribution with small degrees of freedom can be substituted for the Normal.

BUGS code for a two-stage Normal hierarchical model with vague priors and borrowing for both means and variances

The following is the BUGS code which will carry out the MCMC simulation to estimate the consensus mean.

model model1;
const
N=13, k=7;
var
theta[k], sigma[k],mu,tau, Y[N], lab[N], a, b, vw[k], vb;
data lab, Y in "fat.dat";
inits in "fat.in";

{
mu    ~ dnorm(0, 1.0E-4);
 tau   ~ dgamma(1.0E-4, 1.0E-4);
 a     ~ dexp(1.0);
 b     ~ dgamma(1.0E-4, 1.0E-4);
 vb    <- 1/tau;

 for(i in 1:k) {
   sigma[i]  ~dgamma(a,b);
   theta[i]  ~ dnorm(mu, tau);  
   vw[i] <- 1.0/(sigma[i]);  } 

for (i in 1:N) {
   
       Y[i] ~ dnorm( theta[lab[i]],sigma[lab[i]]);
      }
        }

Model 2

In some cases, the assumption that the variances are related may not be appropriate. In that case the following model should be used.

delta(i) ~ N(mu,tau)

and

where

tau ~ IGamma(0.0001,0.0001)

This model does not have the property of model 1 that allows for pooling of data in the estimation of the lab variances. For that reason, when sample sizes are small for some of the labs, the precision of the posterior estimates of the means and variances will be smaller than for model 1.

BUGS code for a two-stage Normal hierarchical model with vague priors and borrowing for means only

The following is the BUGS code for this model:

model model2;
const
N= 13, k=7;
var
theta[k], sigma[k],mu,tau, Y[N], lab[N], vw[k], vb;
data lab, Y in "fat.dat";
inits in "fatonest.in";

{
 mu    ~ dnorm(0, 1.0E-4);
 tau   ~ dgamma(1.0E-4, 1.0E-4);
 
 for(i in 1:k) {
   sigma[i]  ~dgamma(1.0E-4, 1.0E-4);
   theta[i]  ~ dnorm(mu, tau);  
   vw[i] <- 1.0/(sigma[i]);  } 

for (i in 1:N) {
   
       Y[i] ~ dnorm( theta[lab[i]],sigma[lab[i]]);
      }
        }

Example

Data

Data from an experiment which weighed the amount of solids in a sample of dried fish tissue is given in the following table.

**Table 1. Amount of solids measured per sample per lab.**
Lab	solids


1	28.58
1	28.98
2	28.41
2	28.58
3	28.86
3	28.72
4	29.3
4	28.7
5	28.6
6	28.64
6	28.75
7	28.78
7	29.31
8	28.71
8	28.89

9	28.5
9	28.6

Results

Both models were applied with the results given in the following table.

**Table 2. Consensus means, posterior standard deviations, and 95% HPDs for the two models.**
Model	Consensus Mean	Standard Deviation	95% HPD

1	2.871E+1	6.116E-2	(2.861E+1,2.884E+1)
2	2.868E+1	7.077E-2	(2.858E+1,2.884E+1)

Interpretation of the Results

It is clear from the table that the estimates of the consensus means are very close and that there is a slight increase in the size of the posterior HPD due to the reduced pooling of Model 2.

It is interesting to note that lab 5 had only one observation. This causes most classical consensus mean methods to fail because they require at least two data points to estimate each labs variance. The Bayesian models can handle this situation, even in the case of model 2 where there is no pooling for variance estimation. This is due to the fact that a proper prior (probability distribution) is used and so the variance estimate is based on the prior together with any data that is available. For a comparison, the DATAPLOT Consensus Mean procedure was run on this data and produced the following two estimates:

Method	Consensus Mean	95% Confidence Interval

The Mean of Means	2.875E+1	(2.860E+1, 2.889E+1)
Grand Mean	2.875E+1	(2.863E+1, 2.888E+1)

Date created: 9/20/2001
Last updated: 9/20/2001
Please email comments on this WWW page to sedwww@nist.gov.