[R-sig-ME] Discrepancy between mcmcsamp value and variance component esitamate

[Alan Bergland]

Hi all,

	I am getting a funny result using the mcmcsamp function, perhaps it  
is a bug or perhaps there are convergence errors.  I get this result  
using the latest stable version of lme4 on CRAN and the alpha version  
on R-forge, on both a Mac and linux machines.  I apologize for not  
being able to construct a self contained example, so perhaps my  
output will suffice.

My data is structured in the following way:
 > str(ov)
'data.frame':   7958 obs. of  7 variables:
  $ ovn  : num  27 16 29 22 13 18 17 12 16 18 ...
  $ geno : Factor w/ 214 levels "2","3","10","12",..: 52 99 25 25 41  
32 24 32 99 207 ...
  $ food : Factor w/ 4 levels "0.2","0.4","0.6",..: 3 3 3 3 3 3 3 3 3  
3 ...
  $ block: Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
  $ gf   : Factor w/ 856 levels "2:0.8","2:0.6",..: 206 394 98 98 162  
126 94 126 394 826 ...
  $ gb   : Factor w/ 642 levels "2:1","2:2","2:3",..: 154 295 73 73  
121 94 70 94 295 619 ...
  $ bf   : Factor w/ 12 levels "1:0.8","1:0.6",..: 2 2 2 2 2 2 2 2 2  
2 ...

where gf is the interaction of 'geno' and 'food'; gb is the  
interaction of 'geno' and 'block'; and 'bf' is the interaction of  
'block' and 'food'.  There is no geno:food:block interaction.


I fit the following model.  I should also note that lmer2 produces  
the same results.

 > full.ovn<-lmer(ovn~1+(1|block)+(1|gb)+(1|gf)+(1|bf)+(1|geno)+(1| 
food), ov)
 > summary(full.ovn)
Linear mixed-effects model fit by REML
Formula: ovn ~ 1 + (1 | block) + (1 | gb) + (1 | gf) + (1 | bf) + (1  
|      geno) + (1 | food)
    Data: ov
    AIC   BIC logLik MLdeviance REMLdeviance
  49097 49146 -24542      49087        49083
Random effects:
  Groups   Name        Variance Std.Dev.
  gf       (Intercept)  6.43894 2.53751
  gb       (Intercept)  6.00979 2.45149
  geno     (Intercept)  3.37417 1.83689
  bf       (Intercept)  0.84951 0.92169
  food     (Intercept) 19.30673 4.39394
  block    (Intercept)  0.35833 0.59861
  Residual             22.38789 4.73158
number of obs: 7958, groups: gf, 787; gb, 527; geno, 214; bf, 12;  
food, 4; block, 3

Fixed effects:
             Estimate Std. Error t value
(Intercept)   21.009      2.249   9.342


And, then run mcmcsamp:
 > full.mcmc<-mcmcsamp(full.ovn, n=10000)
 > summary(full.mcmc)

Iterations = 1:10000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000

1. Empirical mean and standard deviation for each variable,
    plus standard error of the mean:

                    Mean       SD  Naive SE Time-series SE
(Intercept)     21.0016  3.82131 0.0382131      0.0350391
log(sigma^2)     3.1106  0.01734 0.0001734      0.0003096
log(gf.(In))     1.7228  0.08343 0.0008343      0.0019563
log(gb.(In))     1.8292  0.09193 0.0009193      0.0019939
log(geno.(In))   1.6577  0.13624 0.0013624      0.0015901
log(bf.(In))     0.5115  0.62757 0.0062757      0.0119208
log(food.(In))   3.2965  1.01289 0.0101289      0.0154601
log(blck.(In)) -50.0119 43.12205 0.4312205      5.3344899

2. Quantiles for each variable:

                     2.5%       25%      50%      75%   97.5%
(Intercept)      13.7807  19.34548  21.0349  22.6964 27.9320
log(sigma^2)      3.0763   3.09890   3.1105   3.1222  3.1446
log(gf.(In))      1.5619   1.66684   1.7226   1.7782  1.8891
log(gb.(In))      1.6491   1.76648   1.8300   1.8926  2.0060
log(geno.(In))    1.3956   1.56445   1.6588   1.7480  1.9267
log(bf.(In))     -0.6863   0.08243   0.5012   0.9256  1.7489
log(food.(In))    1.6572   2.58508   3.1775   3.8774  5.6251
log(blck.(In)) -152.0792 -61.58789 -38.9542 -19.1747 -0.4419

Warning message:
In glm.fit(x = X, y = Y, weights = weights, start = start, etastart =  
etastart,  :
   algorithm did not converge

and then, the HPD interval:

 > ovn.hpd
                      lower     upper
(Intercept)      13.744164 27.856893
log(sigma^2)      3.075533  3.143706
log(gf.(In))      1.564024  1.890378
log(gb.(In))      1.650405  2.006905
log(geno.(In))    1.399007  1.929362
log(bf.(In))     -0.721077  1.706577
log(food.(In))    1.518141  5.352884
log(blck.(In)) -142.946749  1.565227
attr(,"Probability")
[1] 0.95
 >

As you can see, the HPD interval for the 'geno' variance component do  
no enclose the estimate for the variance component itself.  Further,  
the mean value of 'geno' from the MCMC run does not equal the  
variance component for 'geno'.  For other variance components there  
is no problem.

Is the warning message at the end of the summary(full.mcmc) a clue  
that there is just poor convergence?  Or, is this a bug?

Thanks for the help,
Alan