Variance components model

Here is a simple illustration of an variance components model given by "Dyes" in the WinBUGS 1.4 Examples, volume 1:

******************************************************
Box and Tiao (1973) analyse data first presented by Davies (1967) concerning batch to batch variation in yields of dyestuff. The data (shown below) arise from a balanced experiment whereby the total product yield was determined for 5 samples from each of 6 randomly chosen batches of raw material.

Batch Yield (in grams)
_______________________________________
1 1545 1440 1440 1520 1580
2 1540 1555 1490 1560 1495
3 1595 1550 1605 1510 1560
4 1445 1440 1595 1465 1545
5 1595 1630 1515 1635 1625
6 1520 1455 1450 1480 1445
*******************************************************

Let denote the jth observation in batch i. To determine the relative importance of between batch variation versus sampling variation, we fit the multilevel model.

1. is N()

2. are iid N(0,

3. assigned a uniform prior

In this situation, the focus is on the marginal posterior distribution of . It is possible to analytically integrate out the random effects , resulting in the marginal posterior
density



where is the "within batch" sum of squares for the ith batch. To use the computational algorithms in LearnBayes, we consider the log posterior distribution of
that is programed in the function logpostnorm1:

logpostnorm1=function(theta,y)
{
mu = theta[1]; sigma.y = exp(theta[2]); sigma.b = exp(theta[3])
p.means=apply(y,1,mean); n=dim(y)[2]
like1=-(apply(sweep(y,1,p.means)^2,1,sum))/2/sigma.y^2-n*log(sigma.y)
like2=-(p.means-mu)^2/2/(sigma.y^2/n+sigma.b^2)-.5*log(sigma.y^2/n+sigma.b^2)
return(sum(like1+like2)+theta[2]+theta[3])
}

In the following R code, I load the LearnBayes package and read in the function logpostnorm1.R and the Dyes dataset stored in "dyes.txt".

Then I summarize the posterior by use of the laplace function -- the mode of () is (3.80, 3.79).

> library(LearnBayes)
> source("logpostnorm1.R")
> y=read.table("dyes.txt")
> fit=laplace(logpostnorm1,c(1500,3,3),y)
> fit$mode
[,1] [,2] [,3]
[1,] 1527.5 3.804004 3.787452