[R-lang] Fwd: [R-sig-ME] lmer & mcmcsamp bug?
Hugh Rabagliati
hugh.rabagliati at gmail.com
Thu Feb 26 17:44:58 PST 2009
Hi all,
So it turns out to be a bug. I've attached Douglas Bates' response
below. It doesn't seem like this should be a hard problem to fix, so
hopefully there'll be an update soonish.
Thanks for the replies,
Hugh
Begin forwarded message:
> From: Douglas Bates <bates at stat.wisc.edu>
> Date: February 26, 2009 5:51:58 PM EST
> To: Hugh Rabagliati <hugh.rabagliati at gmail.com>
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] lmer & mcmcsamp bug?
>
> On Thu, Feb 26, 2009 at 4:21 PM, Hugh Rabagliati
> <hugh.rabagliati at gmail.com> wrote:
>> Hi all,
>>
>> Yesterday I noticed what I take to be a bug in the current version
>> of lme4,
>> and the folks over at the R-lang forum suggested checking in about
>> it here.
>>
>> If I create a mer object using lmer, use it as an argument for
>> mcmcsamp
>> (sampling > 1 times) and then examine my mer object again, I
>> notice a rather
>> peculiar thing. In particular, all of the variance/standard error
>> terms
>> change, as do the associated t values for the fixed effects. The
>> estimated
>> coefficients are unaffected. I figure this is a bug, because I
>> can't see any
>> reason why mcmcsamp would want to do this. I took a look through
>> the code
>> for mcmcsamp, but I don't speak C and nothing jumped out at me.
>> Certainly
>> the function looks like its only meant to return a mermcmc object,
>> so I have
>> no idea why its messing with my mer.
>>
>> I've got this same result on a couple of different computers (all
>> macs). It
>> does, however, seem to be specific to either the version of lmer (
>> 0.999375-28) or of R (2.8.1) that I'm using. I tried it on an old
>> PC version
>> of R (2.5.1) using lme4 version 0.99875-9, and the same problems
>> don't
>> happen then. I've included the output from both the PC and mac
>> versions
>> below.
>>
>> Sage advice, comments, or confirmation that this is a known bug (I
>> couldn't
>> find mention of it elsewhere) would be welcome.
>
> Well, it's a bug.
>
> The C code called by mcmcsamp does something naughty - it changes the
> value of slots of the fitted model object in place. I plead the usual
> excuse for such inexcusable behavior: efficiency. If one copies the
> whole fitted model object at each step in the MCMC iterations the
> function would only be applicable to small models and would take
> forever, even on those. (Actually I guess such a statement further
> makes me guilty of Knuth's "root of all evil" - premature optimization
> - because I haven't actually checked that.)
>
> The code at the end of the C function mer_MCMCsamp is supposed to
> restore the original values. I guess some "infelicities", as Bill
> Venables calls them, must have crept in. I'll take a look.
>
>> ########
>> #
>> #PC Code & output - R v. 2.5.1 & lme4 v. 0.99875-9.
>> #
>> (fm1 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject),
>> sleepstudy))
>> #Linear mixed-effects model fit by REML
>> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
>> # Data: sleepstudy
>> # AIC BIC logLik MLdeviance REMLdeviance
>> # 1752 1764 -871.8 1752 1744
>> #Random effects:
>> # Groups Name Variance Std.Dev.
>> # Subject (Intercept) 627.508 25.0501
>> # Subject Days 35.858 5.9881
>> # Residual 653.590 25.5654
>> #number of obs: 180, groups: Subject, 18; Subject, 18
>>
>> #Fixed effects:
>> # Estimate Std. Error t value
>> #(Intercept) 251.405 6.885 36.51
>> #Days 10.467 1.560 6.71
>>
>> #Correlation of Fixed Effects:
>> # (Intr)
>> #Days -0.184
>>
>> fm1 -> fm1.old
>> samp0 <- mcmcsamp(fm1, n = 1000)
>> fm1
>> #Linear mixed-effects model fit by REML
>> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
>> # Data: sleepstudy
>> # AIC BIC logLik MLdeviance REMLdeviance
>> # 1752 1764 -871.8 1752 1744
>> #Random effects:
>> # Groups Name Variance Std.Dev.
>> # Subject (Intercept) 627.508 25.0501
>> # Subject Days 35.858 5.9881
>> # Residual 653.590 25.5654
>> #number of obs: 180, groups: Subject, 18; Subject, 18
>>
>> #Fixed effects:
>> # Estimate Std. Error t value
>> #(Intercept) 251.405 6.885 36.51
>> #Days 10.467 1.560 6.71
>>
>> #Correlation of Fixed Effects:
>> # (Intr)
>> #Days -0.184
>>
>> #
>> # As you can see, the estimates don't change here
>> #
>> ###########
>>
>>
>> #########
>> #
>> # Mac R v. 2.8.1, ‘lme4’ version 0.999375-28
>> #
>> # I make two mers, which should be identical (I make two because
>> there are
>> some assignment weirdnesses going on here too, #which I haven't yet
>> understood)
>>
>> (fm1.to_mcmc <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|
>> Subject),
>> sleepstudy))
>>
>> #Linear mixed model fit by REML
>> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
>> # Data: sleepstudy
>> # AIC BIC logLik deviance REMLdev
>> # 1754 1770 -871.8 1752 1744
>> #Random effects:
>> # Groups Name Variance Std.Dev.
>> # Subject (Intercept) 627.568 25.0513
>> # Subject Days 35.858 5.9882
>> # Residual 653.584 25.5653
>> #Number of obs: 180, groups: Subject, 18
>>
>> #Fixed effects:
>> # Estimate Std. Error t value
>> #(Intercept) 251.405 6.885 36.51
>> #Days 10.467 1.559 6.71
>>
>> #Correlation of Fixed Effects:
>> # (Intr)
>> #Days -0.184
>>
>>
>> (fm1.not_to_mcmc <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|
>> Subject),
>> sleepstudy))
>>
>> #Linear mixed model fit by REML
>> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
>> # Data: sleepstudy
>> # AIC BIC logLik deviance REMLdev
>> # 1754 1770 -871.8 1752 1744
>> #Random effects:
>> # Groups Name Variance Std.Dev.
>> # Subject (Intercept) 627.568 25.0513
>> # Subject Days 35.858 5.9882
>> # Residual 653.584 25.5653
>> #Number of obs: 180, groups: Subject, 18
>>
>> #Fixed effects:
>> # Estimate Std. Error t value
>> #(Intercept) 251.405 6.885 36.51
>> #Days 10.467 1.559 6.71
>>
>> #Correlation of Fixed Effects:
>> # (Intr)
>> #Days -0.184
>>
>>
>> samp0 <- mcmcsamp(fm1.to_mcmc, n = 1000)
>> fm1.to_mcmc
>>
>> #Linear mixed model fit by REML
>> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
>> # Data: sleepstudy
>> # AIC BIC logLik deviance REMLdev
>> # 1763 1779 -876.7 1761 1753
>> #Random effects:
>> # Groups Name Variance Std.Dev.
>> # Subject (Intercept) 868.398 29.469
>> # Subject Days 49.619 7.044
>> #Residual 904.398 30.073
>> #Number of obs: 180, groups: Subject, 18
>>
>> #Fixed effects:
>> # Estimate Std. Error t value
>> #(Intercept) 251.405 5.260 47.79
>> #Days 10.467 1.518 6.90
>>
>> # All the variances etc are different from before
>>
>> #Correlation of Fixed Effects:
>> # (Intr)
>> #Days -0.343
>>
>> fm1.not_to_mcmc
>>
>> #Linear mixed model fit by REML
>> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
>> # Data: sleepstudy
>> # AIC BIC logLik deviance REMLdev
>> # 1754 1770 -871.8 1752 1744
>> #Random effects:
>> # Groups Name Variance Std.Dev.
>> # Subject (Intercept) 627.568 25.0513
>> # Subject Days 35.858 5.9882
>> # Residual 653.584 25.5653
>> #Number of obs: 180, groups: Subject, 18
>>
>> #Fixed effects:
>> # Estimate Std. Error t value
>> #(Intercept) 251.405 6.885 36.51
>> #Days 10.467 1.559 6.71
>>
>> # Variances are unaffected here
>>
>> #Correlation of Fixed Effects:
>> # (Intr)
>> #Days -0.184
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://pidgin.ucsd.edu/pipermail/r-lang/attachments/20090226/d5d758d4/attachment-0001.htm>
More information about the R-lang
mailing list