Dear list,
I have a question concerning autocorrelation models as estimated by nlme.
When I try to recover an autocorrelative (AR1) effect from simulated data,
I do not get the expected results, that is nlme gives different estimations
than I would expect based on the parameters used to simulate the data
(reproducible example below).
### autocor with random effects
set.seed(12345)
phi <- 0.5 # autocorrelation
b <- 0.3 # effect of predictor x
x <- sample(1:20, 1000, TRUE) # predictor x
ID <- rep(1:100, each = 10) # person IDs
r <- rnorm(100, sd = 1) # random Intercepts with sd = 1
e <- rnorm(1000, sd = 2) # residuals with sd = 2
y <- rep(0,1000) # initialize y
y[1] <- b*x[1] + r[ID[1]] + e[1]
for (i in 2:1000) {
if (ID[i] == ID[i-1])
y[i] <- phi * y[i-1] + b*x[i] + r[ID[i]] + e[i] # not the first
observation for this ID
else y[i] <- b*x[i] + r[ID[i]] + e[i] # first observation for this ID
}
data <- data.frame(y = y, x = x, ID = ID, time = rep(1:10, 100))
head(data)
p <- ggplot(data, aes(x=time, y=y, group=ID))
p + geom_line()
### fit with AR1 effect
library(nlme)
fit_lmeA <- lme(y ~ x, data = data, random = ~ 1 | ID, cor = corAR1(form =
~1|ID))
summary(fit_lmeA)
Accordingly, the model I think of as AR1 does not seem to be the one nlme
fits to the data. Does someone know what type of model nlme actually
assumes in this situation?
Thanks for your help!
Best
Paul
[[alternative HTML version deleted]]