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

[Jake Bowers]

Hi All, There are 8 groups in the second example versus 66 in the  
first. Are we seeing a consistency/central limit problem. Thus,  
mcmcsamp() works when the posterior is plausibly mvnormal but not  
otherwise and in this particular example the information content of  
the data is too low? How's that for a guess? Jake

http://jakebowers.org


On Dec 5, 2009, at 2:21 AM, John Maindonald  
<john.maindonald at anu.edu.au> wrote:

> 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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models