[R-sig-ME] Failure of mcmcsamp()

[Douglas Bates]

I also saw your later message about getting a warning in the earlier
version of lme4.

You have encountered the reason that the development version of the
lme4 package is not yet the released version.  The development version
handles singular variance-covariance matrices more cleanly than does
the release version because of the way that the log-likelihood is
computed.  The release version uses the older-style calculation based
on the precision matrix (i.e. the inverse of the variance-covariance
matrix) of the random effects.  As you might imagine one encounters
difficulties when that variance-covariance matrix approaches
singularity as it can.  The ML or REML estimates of individual
variance components can be zero or, more generally, the
variance-covariance matrix of the random effects can be singular as in
the example below.  In the development version of the lme4 package the
lmer function uses a parameterization of the variance-covariance
matrix and a method of evaluating the log-likelihood at the parameter
values that behaves smoothly as the matrix approaches singularity.
The current implementation of the mcmcsamp function is not as
well-suited to near-singular or singular variances in the random
effects.  There are ways around this.  Gelman et al. (J. Comp. Graph.
Statist., 2008) describe one such way.  I think I have another way
that fits into this framework more cleanly but I still need to
implement and test it.


On Fri, May 2, 2008 at 8:59 PM, John Maindonald
<john.maindonald at anu.edu.au> wrote:
> Dear Douglas -
> The following happens under the versions of lme4 and Matrix from R-forge:
>
>> Orthodont$logdist <- log(Orthodont$distance)
>> keep <- !(Orthodont$Subject%in%c("M04","M13"))
>> orthodont <- subset(Orthodont, keep)
>> orthodont$Subject <- factor(orthodont$Subject)
>> orthdiff.lmer <- lmer(logdist ~ Sex * I(age-11) + (I(age-11) |
> +                                                    Subject),
> +                       data=orthodont,  method="ML")
>> orthdiffr.lmer <- update(orthdiff.lmer, method="REML")
>> orth.mcmc <- mcmcsamp(orthdiffr.lmer, n=1000)
> Error in .local(object, n, verbose, ...) :
>  crossproduct matrix 1 is not positive definite
>
>> VarCorr(orthdiffr.lmer)
> $Subject
>            (Intercept) I(age - 11)
> (Intercept)   1.746e-02   1.555e-07
> I(age - 11)   1.555e-07   1.384e-12
> attr(,"stddev")
> (Intercept) I(age - 11)
>  1.322e-01   1.177e-06
> attr(,"correlation")
>            (Intercept) I(age - 11)
> (Intercept)           1           1
> I(age - 11)           1           1
>
> attr(,"sc")
> sigmaREML
>  0.05238
>
> R version 2.7.0 Patched (2008-05-01 r45578)
> i386-apple-darwin8.10.1
>
> locale:
> C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] DAAG_0.97          MASS_7.2-41        lme4_0.999375-14
> Matrix_0.999375-10
> [5] lattice_0.17-6
>
> loaded via a namespace (and not attached):
> [1] grid_2.7.0  tools_2.7.0
>
> This proceeds without error when I use the CRAN versions (albeit under
> Windows)
>
> other attached packages:
> [1] coda_0.13-1       MEMSS_0.2-4       lme4_0.99875-9    Matrix_0.999375-9
> lattice_0.17-6
>
> If I do not remove the "outliers" the calculations proceed without
> complaint.  The issue is,
> maybe, that the variance component associated with the I(age-11) slope is
> rather small?
> Regards
>
> 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.
>
>