[R] Perhaps Off-topic lme question

[Douglas Bates]

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"))
Read 82 items
 > 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.