Hi Hugh, Yes, I've been able to reproduce this behavior (on OS X R 2.8.1, same lme4 version as you). I don't see any reason why this should happen and it violates R's functional-programming modus operandi. I suggest you report it to the R-sig-ME list as Doug Bates will probably be grateful.
Best Roger Hugh Rabagliati wrote: > Hi all, > > This is a bit of a general question, but it arose out of using the > pvals.fnc function so I figured this forum might have a few ideas about it. > > The issue is whether there's a very odd bug in the mcmcsamp package of > lme4? (or whether I still don't understand mcmc methods well enough). > > If I create a mer object using lmer, use it as an argument for mcmcsamp > (sampling > 1 times), assign the output to a new mermcmc object 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 code is only supposed 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. > > Thanks very much, > > Hugh Rabagliati > > > ######## > # > #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-lang mailing list > [email protected] > http://pidgin.ucsd.edu/mailman/listinfo/r-lang -- Roger Levy Email: [email protected] Assistant Professor Phone: 858-534-7219 Department of Linguistics Fax: 858-534-4789 UC San Diego Web: http://ling.ucsd.edu/~rlevy _______________________________________________ R-lang mailing list [email protected] http://pidgin.ucsd.edu/mailman/listinfo/r-lang
