I would like to document my findings (with a potential FIX) regarding the
issue of calculation of the degrees of freedom with nlme().
The program given at the bottom of this email generates and fit 20 data
sets with a mixed-effects LINEAR model, but using the function nlme()
instead of lme(). In each case, the correct number of degrees of freedom
for the intercept parameter is 12. However, in 7 of the 20 sets, round-off
errors cause the calculated number of degrees of freedom to be 2 instead
of 12. Moreover, in 3 cases, the algorithm failed to converge (represented
by NA).
> DF
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
2 NA 2 NA 12 12 12 NA 12 12 2 2 12 2 12 12 2 12 2 12
A potential fix for this issue would be to replace line 7282 of
R/library/nlme/R/nlme from this expression:
if (any(notIntX <- apply(X, 2, function(el) any(el != el[1])))) {
to this expression:
if (any(notIntX <- apply(X, 2, function(el) any(abs(el -
el[1])>10*sqrt(.Machine$double.eps))))) {
The constant 10*sqrt(.Machine$double.eps) is rather arbitrary and I am not
sure which value would be best to use. However, this correction allows to
calculate correct degrees of freedom in the simulations, although it does
not correct convergence failures, as shown below. I believe those failures
are happening in the C code, which I haven't deeply examined.
> DF
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
12 NA 12 NA 12 12 12 NA 12 12 12 12 12 12 12 12 12 12 12 12
I agree that in principle linear mixed-effects models should be fit with
lme(), not nlme(). However, nlme() may return the wrong number of degrees
of freedom without warning and because of that I believe this issue should
be considered seriously.
Sincerely,
Jerome Asselin
nlme package version: 3.1-43
R version: 1.7.1 running on Red Hat Linux 7.2
###########
library(nlme)
nseeds <- 20
DF <- numeric(nseeds)
for(i in 1:nseeds)
{
set.seed(i)
a <- 2
x <- rep(rnorm(3),rep(5,3))
id <- rep(c("a","b","c"),rep(5,3))
y <- a+x+rnorm(15)
data <- data.frame(y=y,id=id)
initx <-
matrix(x[c(1,6,11)],dimnames=list(c("a","b","c"),"x"))
tryerror <- try(fit.nlme <-
nlme(y ~ a + x, fixed= a~1, random=x~1|id,
data=data, start=list(fixed=c(a=2),
random=list(id=initx)),method="ML"))
DF[i] <-
ifelse(attr(tryerror,"class")[1]=="try-error",NA,fit.nlme$fixDF$X)
}
names(DF) <- 1:nseeds
DF
###########
--
Jerome Asselin (Jérôme), Statistical Analyst
British Columbia Centre for Excellence in HIV/AIDS
St. Paul's Hospital, 608 - 1081 Burrard Street
Vancouver, British Columbia, CANADA V6Z 1Y6
Email: jerome@hivnet.ubc.ca
Phone: 604 806-9112 Fax: 604 806-9044
