[R-sig-ME] Intelligent use of mcmcsamp()

[John Maindonald]

Below are two examples of the use of mcmcsamp() with output
from lmer().  The first set of results are believable.  
For the second set of results, both variance estimates
(for 'site' and for 'Residual' or scale) lie outside of the
credible intervals that are obtained.  Assuming I have
used the functions correctly, it seems surprising that
mcmcsamp() would 'fail' on the second example, which is
balanced and where both variances seem well away from
zero.  These results are consistent over repeated 
simulations.

I'd be interested to hear from list members who make regular use of
mcmcsamp(), as well as maybe from Douglas. Is there any advance on
the current routines in the development branch? 

Questions:

(1) Are instances of less obvious failoure common? How does one check?  
(2) Is there are more direct way to get the credible intervals?
(3) What insight is avaiable on why the second example fails?

John Maindonald.


## EXAMPLE 1:

library(DAAG)
science1.lmer <- lmer(like ~ sex + PrivPub + (1 | school:class),
                     data = science)
> science1.lmer
> ...
Random effects:
Groups       Name        Variance Std.Dev.
school:class (Intercept) 0.321    0.566   
Residual                 3.052    1.747   
Number of obs: 1383, groups: school:class, 66
> ...

science1.mcmc <- mcmcsamp(science1.lmer , n=1000)
z <- VarCorr(science1.mcmc, type="varcov")
colnames(z) <- c("school:class", "Residual")

> ## The following are believable!
> t(apply(z,2, function(x)quantile(x, prob=c(.025, .975))))
              2.5%  97.5%
school:class 0.1442 0.4334
Residual     2.8601 3.3427


## EXAMPLE 2:
ant111b.lmer <- lmer(harvwt ~ 1 + (1 | site), data=ant111b)
> ant111b.lmer
...
Random effects:
Groups   Name        Variance Std.Dev.
site     (Intercept) 2.368    1.54    
Residual             0.578    0.76
Number of obs: 32, groups: site, 8    
> ...

ant111b.mcmc <- mcmcsamp(ant111b.lmer, n=1000)
z <- VarCorr(ant111b.mcmc, type="varcov")
colnames(z) <- c("site", "Residual")
> t(apply(z,2, function(x)quantile(x, prob=c(.025, .975))))
          2.5% 97.5%
site     0.2376 1.878
Residual 0.6370 2.488


John Maindonald             email: john.maindonald at anu.edu.au
phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
Centre for Mathematics & Its Applications, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.
http://www.maths.anu.edu.au/~johnm