From deville.yves at alpestat.com  Tue Jan 22 13:09:40 2013
From: deville.yves at alpestat.com (Yves Deville)
Date: Tue, 22 Jan 2013 13:09:40 +0100
Subject: [Rd] ordering in 'gnls' with 'corCompSymm' corStruct
Message-ID: <50FE8184.5060309@alpestat.com>

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
##----------------------------------------------------------------------------------



[Rd] ordering in 'gnls' with 'corCompSymm' corStruct
Yves Deville