Hi all,
So it turns out to be a bug. I've attached Douglas Bates' response
below. It doesn't seem like this should be a hard problem to fix, so
hopefully there'll be an update soonish.
Thanks for the replies,
Hugh
Begin forwarded message:
From: Douglas Bates <[email protected]>
Date: February 26, 2009 5:51:58 PM EST
To: Hugh Rabagliati <[email protected]>
Cc: [email protected]
Subject: Re: [R-sig-ME] lmer & mcmcsamp bug?
On Thu, Feb 26, 2009 at 4:21 PM, Hugh Rabagliati
<[email protected]> 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
_______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________
R-lang mailing list
[email protected]
http://pidgin.ucsd.edu/mailman/listinfo/r-lang
