[R-sig-ME] Side effects from mcmcsamp()
John Maindonald
john.maindonald at anu.edu.au
Sat May 16 00:26:05 CEST 2009
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
More information about the R-sig-mixed-models
mailing list