Greetings,
I noticed a bug in the getVarCov function from nlme. I am posting here
because it is not currently possible to register and file a report through
https://bugs.r-project.org/. (If this is not the appropriate venue for
this, I'd be grateful if someone could point me to the right place.)
The issue can be seen by observing that getVarCov is sensitive to the order
in which the data are sorted, as demonstrated in the following:
library(nlme)
data(Ovary)
gls_raw <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), data = Ovary,
correlation = corAR1(form = ~ 1 | Mare),
weights = varPower())
Mares <- levels(gls_raw$groups)
V_raw <- lapply(Mares, function(g) getVarCov(gls_raw, individual = g))
Ovary_sorted <- Ovary[with(Ovary, order(Mare, Time)),]
gls_sorted <- update(gls_raw, data = Ovary_sorted)
V_sorted <- lapply(Mares, function(g) getVarCov(gls_sorted, individual = g))
all.equal(gls_raw$modelStruct, gls_sorted$modelStruct)
all.equal(V_raw, V_sorted)
See here for more details and a simple patch:
http://jepusto.github.io//Bug-in-nlme-getVarCov
Or here for just the R code:
https://gist.github.com/jepusto/5477dbe3efa992a3d42c2073ccb12ce4
James
[[alternative HTML version deleted]]