[R] data order affects glmmPQL

[Jack Tanner]

Is it to be expected that the way a data frame is sorted should affect the 
model fit by glmmPQL?

Example:

library(MASS)
library(nlme)

fit.model <- function(il, model.family) {
  cs <- Initialize(corSymm(form=~1|id), data=il)
  glmmPQL(score~test+coder, random=~1|id, # 
control=lmeControl(msMaxIter=100),
          family=model.family, data=il, correlation=cs)
}

score <- 
c(1,8,1,3,4,4,2,5,3,6,0,3,1,5,0,5,1,11,1,2,4,5,2,4,1,6,1,2,8,16,5,16,3,15,3,12,4,9,2,4,1,8,2,6,4,11,2,9,3,17,2,6)
id <- rep(1:13, rep(4, 13))
test <- gl(2, 1, length(score), labels=c("pre", "post"))
coder <- gl(2, 2, length(score), labels=c("two", "three"))

foo <- data.frame(id, score, test, coder) # Define data frame
print(summary(fit.model(foo, poisson)))

bar <- foo[order(id, score, test, coder),] # Reorder data frame
print(summary(fit.model(bar, poisson)))

The two summaries are clearly different. Is this to be expected? Is there a 
canonical way one should order a data frame before passing it to glmmPQL?