# Re: [R] Perhaps Off-topic lme question

From: Douglas Bates <bates_at_stat.wisc.edu>
Date: Thu 14 Apr 2005 - 00:18:10 EST

Berton Gunter wrote:
> A question on lme() :
>
> details: nlme() in R 2.1.0 beta or 2.0.1
>
> The data,y, consisted of 82 data value in 5 groups of sizes 3 9 8 28 34
> .
>
> I fit a simple one level random effects model by:
>
> myfit <- lme( y~1, rand = ~1|Group)
>
> The REML estimates of between and within Group effects are .0032 and .53,
> respectively; the between group component is essentially zero as is clearly
> evident from a plot of the data. So, thus far, no problem.
>
> However, the confidence interval for the between Groups sd that I get from
> intervals(myfit) goes from essentially 0 to infinity (6 x 10^13, actually).
> I assume that this is because the between component estimate is too close to
> the boundary of 0 so that the likelihood approximations with default
> control values fail, but I would appreciate a more definitive comment from
> someone who knows what they're talking about.
>
> If anyone cares to try this, the data in group order are below (i.e., first
> 3 from Group 1, next 9 from Group 2, etc.).

The REML and ML estimates for the variance component associated with groups in these data are zero but the way they are estimated in the lme function will always provide non-zero estimates. As you have seen the intervals constructed in such cases are essentially [0, infinity).

The lmer function in the lme4 package does somewhat better in that it shows that the estimates are on the boundary of the parameter space. For technical reasons at present that boundary is not at zero but at a very small value - a relative variance of 10^{-10}. (The optimization uses an analytic gradient and I haven't worked out what to give the optimizer for the gradient when the relative variance is zero so I use a small positive value for the boundary instead.) However, if you use a value greater than 2 for the msVerbose control option you will see that the optimizer has converged on the boundary.

> bert <- data.frame(grp = factor(rep(1:5, c(3, 9, 8, 28, 34))), resp = scan("/tmp/bert.txt"))
> library(lme4)
> (fm1 <- lmer(resp ~ 1 + (1|grp), bert, control = list(msV=3))) N = 1, M = 5 machine precision = 2.22045e-16 At X0, 0 variables are exactly at the bounds At iterate 0 f= 132.39 |proj g|= 0.0084942

iterations 1
function evaluations 2
segments explored during Cauchy searches 1 BFGS updates skipped 0
active bounds at final generalized Cauchy point 1 norm of the final projected gradient 0
final function value 132.1

F = 132.1
final value 132.099870
converged
Linear mixed-effects model fit by REML
Formula: resp ~ 1 + (1 | grp)

Data: bert

```       AIC      BIC    logLik MLdeviance REMLdeviance
138.0999 145.3200 -66.04993   128.2635     132.0999
Random effects:
Groups   Name        Variance   Std.Dev.
grp      (Intercept) 2.8325e-11 5.3221e-06
Residual             2.8325e-01 5.3221e-01
```
# of obs: 82, groups: grp, 5

Fixed effects:

Estimate Std. Error DF t value Pr(>|t|) (Intercept) 15.352683 0.058773 81 261.22 < 2.2e-16

The important information from the optimizer is that there is 1 active bound at the final point. Also the estimate for the variance component for grp is exactly 1e-10 times the variance component from the residual.

R-help@stat.math.ethz.ch mailing list