Side effects from mcmcsamp()
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