[Rd] degrees of freedom in nlme() (PR#2384)

[jerome at hivnet.ubc.ca]

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 at hivnet.ubc.ca
Phone: 604 806-9112   Fax: 604 806-9044