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

[Douglas Bates]

I agree that mcmcsamp in current versions of the lme4 package provide
suspect results.  I still have problems formulating an MCMC update for
the variance components when there is a possibility the value getting
close to zero.

In the development branch of the lme4 package I have added profiling
of the deviance with respect to the parameters as a method of
determining the precision of the parameter estimates.  This package is
called lme4a in the R packages tab of
http://lme4.r-forge.r-project.org/  Unfortunately it looks like the
Windows binary is not available at R-forge (and the error message
looks peculiar because it is an error from an R declarations file - it
may be that I don't have the correct sequence of include files).

I usually look at profiles of the change in the deviance by taking the
signed square root transformation.  Results for your models are
enclosed.  In these cases there aren't problems with the variance of
the random effects approaching zero.  In the case of the Dyestuff data
we do get such behavior and the splom plot shows some rough edges as a
result.



On Sat, 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
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Rplots.pdf
Type: application/pdf
Size: 175874 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20091207/e62db89c/attachment.pdf>