naive question about output from mcmcsamp v. lmer estimates

Dear List,

I am  attempting to estimate the quantitative trait statistic Qst for a
colleague. Essentially I am fitting the model

X[ijk] = mu + alpha[i] + beta[ij] + gamma[k] + epsilon[ijk]

to data from measured from trees planted in an unbalanced, incomplete
block experiment, where X[ijk] = individual observation in the ith
provenance, ijth family and kth block; mu, alpha[i], beta[ij], and
gamma[k] are the effects of population mean, provenance, family within
provenance and block respectively. Qst is estimated using variance
components for provenance and family within provenance. Because of the
problems inherent in using point estimates to estimate Qst I would like
to use a Bayesian approach as proposed by Waldmann et al. (Heredity
(2005) 94, 623?629). I have implemented this model in winBugs and get
similar summary output to fitting the following lmer model (from lme4,
using R 2.4.1 on windows):

 >summary(lme.HT02)

Linear mixed-effects model fit by REML
Formula: HT02 ~ (1 | Block) + (1 | Provenance/Family)
    Data: quant
    AIC   BIC logLik MLdeviance REMLdeviance
  12615 12636  -6303      12611        12607
Random effects:
  Groups            Name        Variance Std.Dev.
  Family:Provenance (Intercept)   3.0350  1.7421
  Provenance        (Intercept) 191.7605 13.8478
  Block             (Intercept)  16.1534  4.0191
  Residual                      593.3983 24.3598
number of obs: 1359, groups: Family:Provenance, 335; Provenance, 18;
Block, 14

Fixed effects:
             Estimate Std. Error t value
(Intercept)  123.497      3.549    34.8

As winBugs is very slow for this problem I was hoping to use mcmcsamp
instead. I am being naive in thinking that I can use mcmcsamp to
estimate posterior variance components densities? E.g.

 >colMeans(mcmcsamp(lme.HT02,10000,trans=FALSE))

(Intercept)     sigma^2   Fm:P.(In)   Prvn.(In)   Blck.(In)
   122.57740   571.59628    60.77002    61.59250    17.56870

I was expecting these numbers to be similar, although not identical, to
the lmer output. This is all quite new ground for me so, I readily
accept I have probably gone wrong somewhere. But where?

I was expecting these numbers to be similar, although not identical, to
the lmer output. This is all quite new ground for me so, I readily
accept I have probably gone wrong somewhere. But where?

As a first step I would suggest that you examine the mcmcsamp result
to see if the chains are stable.  Save the output from mcmcsamp as,
say, samp.HT02 and plot the chains as parallel time series

library(coda)
xyplot(samp.HT02)

You might also want to examine empirical density estimates derived
from the chains

densityplot(samp.HT02, plot.points = FALSE)

or normal probability plots

qqmath(samp.HT02, type = c("g", "p"))

If the distribution of the MCMC sample for a parameter is badly skewed
or otherwise unstable then  the sample mean will not be near to the
estimate.