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

[John Maindonald]

This looks promising, thanks for your efforts.  Fortunately the MacOS X 
version runs fine, under R-devel (2.11.0 Under development).

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

On 08/12/2009, at 12:41 AM, Douglas Bates wrote:

> 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
>> 
> <john.Rout><Rplots.pdf>