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}}