Dear R-devel members,
While writing a new correlation structure similar to 'corCompSymm' and
intended to be used with 'gnls', I got puzzled with the
'Initialize' method.
Using 'Initialize' before 'gnls' may be regarded as a mean to
set an
initial value for the corStruct parameter. However 'gnls' does not work
properly with a 'corCompSymm' correlation structure when the data frame
is not suitably ordered by group and when 'Initialize' is used for the
corStruct before calling 'gnls'. Maybe the documentation could warn
about this problem?
As an example, the results of 'fit1' and 'fit2' below are
different, and
only 'fit1' gives the good results. The problem remains if we drop the
0.2 (for the 'value' formal) when calling Initialize.
My explanation is that when 'Initialize' is used with a
'corCompSymm'
object and a data frame which is not sorted by the group variable, the
'start' indices computed by the 'Dim' method are misleading.
These
'start' values will remain attached to the object in 'gnls' even
though
the data frame is sorted therein, since the call to 'Initialize' within
'gnls' will not then change them. If this explanation is true,
'Initialize.corCompSymm' should not accept an unordered 'data'
formal,
or should sort the data before computing 'start' with 'match'.
Best
Yves Deville, statistical consultant. France
##-----------------------------------------------------------------------------------
## generate grouped data
M <- 20
set.seed(123)
Ni <- 1 + rpois(M, lambda = 4); N <- sum(Ni)
grp <- factor(rep(1:M, times = Ni))
x <- lapply(Ni, function(n) 1:n)
y <- lapply(x, function(x) 1 + 2*x + rnorm(1) + rnorm(length(x)))
dataOrd <- data.frame(grp = grp, x = unlist(x), y = unlist(y))
## change order
ind <- sample(1:N); dataUnOrd <- dataOrd[ind, ]
library(nlme)
cS1 <- corCompSymm(0.2, form = ~ 1 | grp)
fit1 <- gnls(model = y ~ a + b*x, data = dataUnOrd, start = c(a = 0.1, b
= 0.3),
correlation = cS1)
## the same with 'Initialize'
cS2 <- corCompSymm(0.2, form = ~ 1 | grp)
cS2 <- Initialize(cS2, data = dataUnOrd)
fit2 <- gnls(model = y ~ a + b*x, data = dataUnOrd, start = c(a = 0.1, b
= 0.3),
correlation = cS2)
## Dim does not give the good 'start' on 'cS2'
cS3 <- corCompSymm(0.2, form = ~ 1 | grp)
cS3 <- Initialize(cS3, data = dataOrd)
Dim(cS2)$start
Dim(cS3)$start
##----------------------------------------------------------------------------------
