[R-lang] lmer objects and mcmcsamp

[Hugh Rabagliati]

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