[Rd] ordering in 'gnls' with 'corCompSymm' corStruct
Yves Deville
deville.yves at alpestat.com
Tue Jan 22 13:09:40 CET 2013
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
##----------------------------------------------------------------------------------
More information about the R-devel
mailing list