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

[John Maindonald]

I wonder whether you have considered simulation from the full
variance-covariance matrix.  For some samples, one or more
variance parameter estimates will be allowed to go negative, 
but the full variance-covariance matrix should still be positive 
definite.  This avoids any problem with variance close to zero,
but it does mean (and I think this a useful side-effect) that 
credible intervals will sometimes have a negative lower bound,
or even lie entirely on the negative real line!

Why a useful side effect?  I've encountered cases where scientists, 
in a field experiment, have chosen blocks that are e.g., long strip 
at right angles to a river bank, thus spanning about as wide a range 
of variation as is possible.  (I've been told about the river bank 
example, the examples in my own experience were a bit different.)
A model in which there is a negative component of variance may 
then give a formally correct variance-covariance structure.  It is a bit 
like using a complex number representation to solve some problems 
on the real line.

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 9:30 AM, John Maindonald wrote:

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