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

[espesser]

I think it's due to the same bug discussed in:

http://finzi.psych.upenn.edu/R/Rhelp08/2009-February/002001.html

R. Espesser
Laboratoire Parole et Langage
CNRS/Université Aix-Marseille


John Maindonald a écrit :
> 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.
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models