[R-sig-ME] lmer & mcmcsamp bug?

[Douglas Bates]

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
>