## Monday, June 23, 2008

### 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 $y_{ij}$ 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. $y_{ij}$ is N($\mu + b_i, \sigma_y^2$)

2. $b_1, ..., b_N$ are iid N(0, $\sigma^2_b)$

3. $(\sigma_y^2, \sigma_b^2)$ assigned a uniform prior

In this situation, the focus is on the marginal posterior distribution of $(\mu, \sigma_y^2, \sigma_b^2)$ . It is possible to analytically integrate out the random effects $b_1, ..., b_N$, resulting in the marginal posterior
density

$\prod_{i=1}^N (\exp\{-\frac{1}{2\sigma_y^2} S_i\} \frac{1}{\sigma_y})$ $\prod_{i=1}^N (\exp\{-\frac{1}{2(\sigma^2_y/n+\sigma^2_b)} (\bar y_i - \mu)^2\} \frac{1}{\sqrt{\sigma^2_y/n+\sigma^2_b}} )$

where $S_i$ 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
$(\mu, \log \sigma_y, \log \sigma_b)$ 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 ($\log \sigma_y, \log \sigma_b$) is (3.80, 3.79).

> library(LearnBayes)
> source("logpostnorm1.R")