[R-sig-ME] Side effects from mcmcsamp()

[John Maindonald]

Dear Douglas -
The following is at the very least a serious trap!  The output
from VarCorr() and print() changes remarkably after running
mcmcsamp() with the lmer object as argument.  The estimate
of sigma is always (for these data) high, but roughly in the
middle of the range of values that come out of ant111b.samp at sigma

I have also been able to reproduce this behaviour
(a) under lme4_0.999375-28 [below, I use 30]
(b) under Mac OSX.  In fact, I thought initially that this
was an OSX problem!

While mcmcsamp() is in view, it would be very helpful
to have a progress report on where it is at.

Many thanks
John Maindonald.

I run the following code:

library(DAAG)
library(lme4)
ant111b.lmer <- lmer(harvwt ~ 1 + (1 | site), data=ant111b)
ant111b.lmer
VarCorr(ant111b.lmer)
ant111b.samp <- mcmcsamp(ant111b.lmer, n=1000)

## Now, extract the VarCorr and summary correlations
VarCorr(ant111b.lmer)
summary(ant111b.lmer)

 > ant111b.lmer <- lmer(harvwt ~ 1 + (1 | site), data=ant111b)
 > ant111b.lmer
Linear mixed model fit by REML
Formula: harvwt ~ 1 + (1 | site)
   Data: ant111b
   AIC   BIC logLik deviance REMLdev
100.4 104.8 -47.21    95.08   94.42
Random effects:
Groups   Name        Variance Std.Dev.
site     (Intercept) 2.36773  1.53874
Residual             0.57754  0.75996
Number of obs: 32, groups: site, 8

Fixed effects:
            Estimate Std. Error t value
(Intercept)   4.2917     0.5603   7.659
 > library(DAAG)
 > library(lme4)
 > ant111b.lmer <- lmer(harvwt ~ 1 + (1 | site), data=ant111b)
 > ant111b.lmer
Linear mixed model fit by REML
Formula: harvwt ~ 1 + (1 | site)
   Data: ant111b
   AIC   BIC logLik deviance REMLdev
100.4 104.8 -47.21    95.08   94.42
Random effects:
Groups   Name        Variance Std.Dev.
site     (Intercept) 2.36773  1.53874
Residual             0.57754  0.75996
Number of obs: 32, groups: site, 8

Fixed effects:
            Estimate Std. Error t value
(Intercept)   4.2917     0.5603   7.659
 > ant111b.samp <- mcmcsamp(ant111b.lmer, n=1000)
 > VarCorr(ant111b.lmer)
$site
            (Intercept)
(Intercept)    8.363436
attr(,"stddev")
(Intercept)
   2.891961
attr(,"correlation")
            (Intercept)
(Intercept)           1

attr(,"sc")
sigmaREML
1.428289
 > summary(ant111b.lmer)
Linear mixed model fit by REML
Formula: harvwt ~ 1 + (1 | site)
   Data: ant111b
   AIC   BIC logLik deviance REMLdev
112.1 116.5 -53.04    105.7   106.1
Random effects:
Groups   Name        Variance Std.Dev.
site     (Intercept) 8.3634   2.8920
Residual             2.0400   1.4283
Number of obs: 32, groups: site, 8

Fixed effects:
            Estimate Std. Error t value
(Intercept)   4.2917     0.4171   10.29

 > sessionInfo()
R version 2.9.0 (2009-04-17)
i386-pc-mingw32

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States. 
1252;LC_MONETARY=English_United States. 
1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] lme4_0.999375-30   Matrix_0.999375-26 lattice_0.17-22
[4] DAAG_0.99-1        MASS_7.2-47

loaded via a namespace (and not attached):
[1] grid_2.9.0



John Maindonald             email: john.maindonald at anu.edu.au
phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
Centre for Mathematics & Its Applications, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.