[R-sig-ME] lmer & mcmcsamp bug?
S Ellison
S.Ellison at lgc.co.uk
Fri Feb 27 00:47:17 CET 2009
This behaviour has also been implicated - correctly or otherwise I can't say - in the 'negative PWRSS' question reported previously (http://www.nabble.com/Problem-with-aovlmer.fnc-in-languageR-td20706128.html, http://www.nabble.com/Re%3A-Problem-with-aovlmer.fnc-in-languageR-p21365322.html), so a fix may kill two birds with one stone. That would be most welcome!
Steve Ellison
>>> Douglas Bates <bates at stat.wisc.edu> 02/26/09 10:51 PM >>>
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
>
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
*******************************************************************
This email and any attachments are confidential. Any use...{{dropped:8}}
More information about the R-sig-mixed-models
mailing list