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