[R-lang] lmer objects and mcmcsamp

[Roger Levy]

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
> R-lang at ling.ucsd.edu
> http://pidgin.ucsd.edu/mailman/listinfo/r-lang


-- 

Roger Levy                      Email: rlevy at ling.ucsd.edu
Assistant Professor             Phone: 858-534-7219
Department of Linguistics       Fax:   858-534-4789
UC San Diego                    Web:   http://ling.ucsd.edu/~rlevy