[R-sig-ME] opposite MCMC results with two lmer versions

[robert espesser]

Dear R Users,

Looking for a teaching  example about the interest of mixed models,
I tried the "shoes data" from the "SMpracticals" package.
It's about the amount of wear in a paired comparison of two materials
used for soling the shoes of 10 boys.
The materials were allocated randomly to the left and right feet.
The dataframe is:	

 >shoe

    material boy foot    y
1         A   1    L 13.2
2         B   1    R 14.0
3         A   2    L  8.2
4         B   2    R  8.8
5         A   3    R 10.9
6         B   3    L 11.2
7         A   4    L 14.3
8         B   4    R 14.2
9         A   5    R 10.7
10        B   5    L 11.8
11        A   6    L  6.6
12        B   6    R  6.4
13        A   7    L  9.5
14        B   7    R  9.8
15        A   8    L 10.8
16        B   8    R 11.3
17        A   9    R  8.8
18        B   9    L  9.3
19        A  10    L 13.3
20        B  10    R 13.6

I obtained very different HPD intervals and pMCMC values according to 
these 2 sessionInfo :

############
A)   Last available versions of lme4 and languageR

 > sessionInfo()
R version 2.8.0 (2008-10-20)
i386-pc-mingw32

locale:
LC_COLLATE=French_France.1252;LC_CTYPE=French_France.1252;LC_MONETARY=French_France.1252;LC_NUMERIC=C;LC_TIME=French_France.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] languageR_0.953    zipfR_0.6-0        lme4_0.999375-27 
Matrix_0.999375-16
[5] lattice_0.17-15

 > lmer(y ~ material +(1|boy),data=shoe) ->  mat.lmer
 > summary(mat.lmer)
Linear mixed model fit by REML
Formula: y ~ material + (1 | boy)
    Data: shoe
    AIC   BIC logLik deviance REMLdev
  62.94 66.92 -27.47    53.82   54.94
Random effects:
  Groups   Name        Variance Std.Dev.
  boy      (Intercept) 6.100889 2.47000
  Residual             0.074944 0.27376
Number of obs: 20, groups: boy, 10

Fixed effects:
             Estimate Std. Error t value
(Intercept)  10.6300     0.7858  13.527
materialB     0.4100     0.1224   3.349

 > pvals.fnc(mat.lmer)
$fixed
             Estimate MCMCmean HPD95lower HPD95upper  pMCMC Pr(>|t|)
(Intercept)    10.63  10.6379      9.126     12.175 0.0001   0.0000
materialB       0.41   0.4098     -1.593      2.345 0.6542   0.0036

$random
     Groups        Name Std.Dev. MCMCmedian MCMCmean HPD95lower HPD95upper
1      boy (Intercept)   2.4700     0.8747   0.8570     0.0000     1.8015
2 Residual               0.2738     2.1350   2.1818     1.2191     3.1653

The MCMC results show materialB is here clearly unsignificant.


############
B)  Previous versions  of lme4 and languageR
same model & data, with :
R version 2.6.2 (2008-02-08)
i386-pc-mingw32

locale:
LC_COLLATE=French_France.1252;LC_CTYPE=French_France.1252;LC_MONETARY=French_France.1252;LC_NUMERIC=C;LC_TIME=French_France.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] languageR_0.4     coda_0.13-1       lme4_0.99875-9    Matrix_0.999375-6
[5] zipfR_0.6-0       lattice_0.17-6

 > lmer(y ~ material +(1|boy),data=shoe) ->  mat_old.lmer
 > summary(mat_old.lmer)
Linear mixed-effects model fit by REML
Formula: y ~ material + (1 | boy)
    Data: shoe
    AIC   BIC logLik MLdeviance REMLdeviance
  60.94 63.92 -27.47      53.82        54.94
Random effects:
  Groups   Name        Variance Std.Dev.
  boy      (Intercept) 6.073101 2.46437
  Residual             0.075285 0.27438
number of obs: 20, groups: boy, 10

Fixed effects:
             Estimate Std. Error t value
(Intercept)  10.6300     0.7841  13.557
materialB     0.4100     0.1227   3.341

 > pvals.fnc(mat_old.lmer)
$fixed
             Estimate MCMCmean HPD95lower HPD95upper  pMCMC Pr(>|t|)
(Intercept)    10.63  10.6182      8.818     12.357 0.0001   0.0000
materialB       0.41   0.4118      0.119      0.684 0.0102   0.0036

$random
          MCMCmean HPD95lower HPD95upper
sigma      0.2928     0.1829     0.4869
boy.(In)   2.5999     1.6610     4.3302

The MCMC results show  materialB is clearly significant


The both versions of  lme4 showed a similar LRT result between the
"null" model and the regular test :
 > anova(empty.lmer,mat.lmer)
Data: shoe
Models:
empty.lmer: y ~ 1 + (1 | boy)
mat.lmer: y ~ material + (1 | boy)
           Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
empty.lmer  3  67.937  70.924 -30.968
mat.lmer   4  61.817  65.800 -26.909 8.1197      1   0.004379 **

I know it's anti-conservative, but it's quite opposite to the
results of MCMC values of A).

The old classical paired t.test  method gave results similar to B)



So, which are the  results I can trust in ?
Thank you very much

Robert Espesser
Laboratoire Parole and Langage, CNRS UMR 6057
Université de Provence
Aix en Provence, France