[R-sig-ME] lmer & mcmcsamp bug?
Hugh Rabagliati
hugh.rabagliati at gmail.com
Thu Feb 26 23:21:21 CET 2009
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.
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
More information about the R-sig-mixed-models
mailing list