[R-sig-ME] naive question about output from mcmcsamp v. lmer estimates

[Douglas Bates]

On 2/7/07, Gibbons,James <afs417 at bangor.ac.uk> wrote:
> 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?

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.