[R] random effects with lme() -- comparison with lm()

[Douglas Bates]

Jerome Asselin <jerome at hivnet.ubc.ca> writes:

> On Thu, 2004-01-15 at 16:30, Douglas Bates wrote:
> <...snip...>
> > (BTW, I wouldn't say that this is equivalent to a fixed effects
> > model.  It is still a random effects model with two variance
> > components.  It just doesn't have well-defined estimates for those two
> > variance components.)
> 
> Agreed.
> 
> <...snip...>
> > You should find that intervals() applied to your fitted model produces
> > huge intervals on the variance components, which is one way of
> > diagnosing an ill-defined or nearly ill-defined model.
> 
> Following your suggestion, I got:
> > intervals(lme(Y~1,data=simdat,random=~1|A))
> Error in intervals.lme(lme(Y ~ 1, data = simdat, random = ~1 | A)) :
>         Cannot get confidence intervals on var-cov components:
> Non-positive definite approximate variance-covariance
> 
> This led me to:
> > lme(Y~1,data=simdat,random=~1|A)$apVar
> [1] "Non-positive definite approximate variance-covariance"
> 
> As a new feature suggestion for lme(), would it be appropriate to use
> "apVar" as a warning flag in this case?

Certainly.

You may know that we are doing a major revision of the lme
computational methods based on the ability to calculate both the
gradient and the Hessian of the profiled log-likelihood, as described
in
        http://www.stat.wisc.edu/~bates/reports/MixedComp.pdf

I think that when we have both the gradient and the Hessian we will be
in a much better situation to diagnose ill-defined estimates.  The
apVar component in the current lme objects is an approximate
variance-covariance matrix from numerical derivatives.  Working with
an exact Hessian should be much more reliable.