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.