Side effects from mcmcsamp()
I have uploaded a new version of the lme4 package to CRAN's incoming directory. This version (0.999375-31) forces a copy of the fitted model object in the mcmcsamp function before passing it to the C code. I enclose a copy of John's example run with the new release of lme4. On Fri, May 15, 2009 at 5:26 PM, John Maindonald
<john.maindonald at anu.edu.au> wrote:
I believe that assignment to a new name creates, in the first place, a promise. ?A new object is created only when obmix0 is changed. ?John Chambers, in "Software for Data Analysis", has a lot to say about promises, though not in this particular connection as far as I could see at a quick check. 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. On 16/05/2009, at 7:52 AM, Ben Bolker wrote:
Juan Pedro Steibel wrote:
Hello, I noticed this a wile ago. One thing that occurred to me was to copy to another object before running mcmcsamp, example: obmix0<-obmix where obmix is the lmer object. then do: mcmcsamp(obmix).... For some reason obmix0 gets modified too.
Yeah, this aspect of the whole thing mystified me.
I understand that mcmcsamp does nasty things internally,
but I don't understand why it would also do them to a *copy* of the
object. ?I poked around in R manuals (language definitions etc.) but
didn't manage to find anything useful for understanding this. ?I can
imagine there's some slick optimization where R says to itself "gee,
this object hasn't been modified from its original yet, so I won't
actually make a new copy ..." but is fooled by the internal
manipulation. ?("Don't anthromorphize computers -- they get really mad"
:-) )
gc()
? ? ?used (Mb) gc trigger (Mb) max used (Mb) Ncells 209832 ?5.7 ? ? 467875 12.5 ? 350000 ?9.4 Vcells 476846 ?3.7 ? ?1165368 ?8.9 ?1028740 ?7.9 ## make a big object -- 70 Mb
z <- runif(1e7) gc()
? ? ? ?used (Mb) gc trigger (Mb) max used (Mb) Ncells ? 209821 ?5.7 ? ? 467875 12.5 ? 350000 ?9.4 Vcells 10476828 80.0 ? 11888118 90.7 10477151 80.0 ## make a copy
z2 <- z
## Vcells only increases a little
gc()
? ? ? ?used (Mb) gc trigger (Mb) max used (Mb) Ncells ? 209826 ?5.7 ? ? 467875 12.5 ? 350000 ?9.4 Vcells 10476829 80.0 ? 12562523 95.9 10477151 80.0 ## ditto
z3 <- z2 gc()
? ? ? ?used (Mb) gc trigger ?(Mb) max used (Mb) Ncells ? 209831 ?5.7 ? ? 467875 ?12.5 ? 350000 ?9.4 Vcells 10476830 80.0 ? 13270649 101.3 10477151 80.0 ## tweak z3
z3[1] <- z3[1]+1 gc()
? ? ? ?used ?(Mb) gc trigger ?(Mb) max used ?(Mb) Ncells ? 209838 ? 5.7 ? ? 467875 ?12.5 ? 350000 ? 9.4 Vcells 20476831 156.3 ? 22913123 174.9 20476845 156.3 ## now memory usage increases I'm sure this is obvious to hardened R programmers, but I can't find any explicit discussion in the R manuals. So I would imagine that making some trivial modification of the object (like ?"class(object) <- class(object)", for example) ?would mark it as a "new" object ... cheers ?Ben
It's been a while since I checked this so it may have been fixed. Anyways, a way around this is to extract all the info needed from your lmer object to other variables before running mcmcsamp, that's what I do now and it works fine. Thanks JP John Maindonald wrote:
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
-- Ben Bolker Associate professor, Biology Dep't, Univ. of Florida bolker at ufl.edu / www.zoology.ufl.edu/bolker GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-------------- next part -------------- R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R.
library(DAAG)
Loading required package: MASS Attaching package: 'DAAG' The following object(s) are masked from package:MASS : hills
library(lme4)
Loading required package: Matrix Loading required package: lattice Attaching package: 'Matrix' The following object(s) are masked from package:stats : xtabs The following object(s) are masked from package:base : rcond
(ant111b.lmer <- lmer(harvwt ~ 1 + (1 | site), data=ant111b))
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) 2.367733
attr(,"stddev")
(Intercept)
1.538744
attr(,"correlation")
(Intercept)
(Intercept) 1
attr(,"sc")
sigmaREML
0.759959
summary(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
sessionInfo()
R version 2.9.0 (2009-04-17) x86_64-pc-linux-gnu locale: LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] lme4_0.999375-31 Matrix_0.999375-27 lattice_0.17-25 DAAG_0.99-1 [5] MASS_7.2-47 loaded via a namespace (and not attached): [1] grid_2.9.0
proc.time()
user system elapsed 14.096 0.224 14.291