[R-sig-ME] Intelligent use of mcmcsamp()
Jake Bowers
jwb1970 at gmail.com
Sat Dec 5 16:00:22 CET 2009
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
More information about the R-sig-mixed-models
mailing list