I'm facing a problem when simulating with the function gls() in the R package nlme. I'd like to impose an AR(1) correlation structure on linear models fitted with gls(). When I do this in a simulation framework (e.g., repeating the whole thing 10,000 times) R usually (but not always) crashes out of the blue at some point. I don't get an error message in R but R abruptly terminates and Windows reports "R for Windows GUI front-end has stopped working". I am using OS Windows 7 and R version 3.0.0 as well as the latest version of nlme (3.1-109).
This is my (simplified) code:
library(nlme)
group <- rep(c("A", "B", "C", "D"), each=5)
person <- as.factor(rep(c(1:5), times=4))
value <- c(8.50, 11.69, 9.07, 8.82, 10.63, 10.87, 10.02, 9.84, 9.82, 9.68,
7.55, 12.19, 9.16, 8.59, 9.98, 11.16, 9.63, 9.73, 10.15, 9.06)
dat <- data.frame(group, person, value)
xyz <- numeric(1)
runs <- 1E4
for(r in 1:runs){
model <- gls(value ~ group - 1, data=dat, corAR1(form=~1|person))
if(summary(model)$corBeta[1,2] < 0.9){xyz <- xyz + 1}
}
xyz
The problem occurs with different datasets as well, not only with the particular one in this example. Moreover, I ran the above code repeatedly on different (but identical) workstations and in different R sessions. R crashes in about 50% of cases. Unfortunately, I have no idea what makes this happen (or not happen), and I would be very grateful if you could help me.