[R-sig-ME] Side effects from mcmcsamp()
Douglas Bates
bates at stat.wisc.edu
Wed May 20 19:56:00 CEST 2009
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
More information about the R-sig-mixed-models
mailing list