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

[jerome@hivnet.ubc.ca]

I believe that a round-off error is causing nlme() to provide a wrong number 
of degrees of freedom. In the example below, getFixDF() returns 2 degrees
of freedom, but it should be 12 degrees of freedom according to lme().

R 1.6.1 on RedHat Linux 7.2
Package: nlme
Version: 3.1-36
Date: 2002-12-28

library(nlme)
set.seed(14)

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"))
rm(a,x,id,y)

summary(fit.lme <-
              lme(y ~ 1,data=data,random=~1|id,method="ML"))
debug(getFixDF)
summary(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"))

# From within getFixDF(), do the following check...
Browse[1]> apply(X, 2, function(el) any(el != el[1])) #Round-off error here!!
   a
TRUE
Browse[1]> X
   a
1  1
2  1
3  1
4  1
5  1
6  1
7  1
8  1
9  1
10 1
11 1
12 1
13 1
14 1
15 1

# If I force X[ ] <- X[1] in this case, then the result of nlme() is 
concordant with that of lme().


# Below is what I get with the current version of the nlme package.

> summary(fit.lme <-
+               lme(y ~ 1,data=data,random=~1|id,method="ML"))
Linear mixed-effects model fit by maximum likelihood
 Data: data
       AIC      BIC    logLik
  46.09706 48.22121 -20.04853

Random effects:
 Formula: ~1 | id
        (Intercept)  Residual
StdDev:   0.5019741 0.8299741

Fixed effects: y ~ 1
               Value Std.Error DF  t-value p-value
(Intercept) 1.671345 0.3730901 12 4.479736   8e-04

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max
-2.1912914 -0.5408147 -0.1515690  0.7275933  1.3373572

Number of Observations: 15
Number of Groups: 3

> summary(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"))
Nonlinear mixed-effects model fit by maximum likelihood
  Model: y ~ a + x
 Data: data
       AIC      BIC    logLik
  46.09706 48.22121 -20.04853

Random effects:
 Formula: x ~ 1 | id
                x  Residual
StdDev: 0.5019711 0.8299748

Fixed effects: a ~ 1
     Value Std.Error DF  t-value p-value
a 1.671345 0.3730888  2 4.479752  0.0464

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max
-2.1912919 -0.5408166 -0.1515696  0.7275931  1.3373593

Number of Observations: 15
Number of Groups: 3