[R-sig-ME] naive question about output from mcmcsamp v. lmer estimates
Douglas Bates
bates at stat.wisc.edu
Wed Feb 7 17:26:27 CET 2007
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.
More information about the R-sig-mixed-models
mailing list