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