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

[jerome@hivnet.ubc.ca]

Full_Name: Jerome Asselin
Version: 1.6.1
OS: RedHat Linux 7.2
Submission from: (NULL) (142.103.173.179)



There is something very queer happening with the degrees
of freedom in nlme(). In the example below, I am fitting
the same model with lme() and nlme(). I am using the
nlme package version 3.1-34.


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"))

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


The number of degrees of freedom for the intercept parameter
is 12 degrees using lme(), but only 2 degrees using nlme().
The same problem is noticed with seeds 15, 16 and 17.

But if I change the seed to 18 and try again, then I get 12
degrees of freedom for both lme() and nlme().


set.seed(18)

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"))

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


If I set the seed to 5 or 8, then nlme() fails despite the
simplicity of the model. The error messages are:
- for set.seed(5): 
    Step halving factor reduced below minimum in PNLS step
- for set.seed(8):
    Singularity in backsolve at level 0, block 1
This issue of degrees of freedom may be related to the
convergence problems stated in another bug report (PR#2369).

Sincerely,
Jerome Asselin

set.seed(5)

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"))

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


set.seed(8)

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"))

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