Question about autocorrelation in nlme

There are also examples on
http://www.ats.ucla.edu/stat/r/examples/alda/ch7.htm (and examples from
their other chapters may also be useful ... Singer and Willet's book is
highly recommended). We go through a similar example using gls in
https://www.researchgate.net/publication/23134911_Multilevel_modeling_Beyond_the_basic_applications

On Wed, Sep 16, 2015 at 8:46 AM, Kevin Wright <kw.stat at gmail.com> wrote:

This might help:

http://stats.stackexchange.com/questions/6469/simple-linear-model-with-autocorrelated-errors-in-r

Adapting that to your example:

### autocor with random effects
set.seed(12345)
b <- 0.3  # effect of predictor x
x <- sample(1:20, 1000, TRUE)  # predictor x (time)
ID <- rep(1:100, each = 10)  # person IDs
r <- rnorm(100, sd = 1)  # random Intercepts with sd = 1
e <- arima.sim(model=list(ar=0.5), n=1000, sd=2) # .5 = autocor
plot(e)
y <- b*x + r[ID] + e

data <- data.frame(y = y, x = x, ID = ID, time = rep(1:10, 100))
head(data)

library(lattice)
xyplot(y~x|ID, data)

### fit with AR1 effect
library(nlme)
fit_lmeA <- lme(y ~ x, data = data, random = ~ 1 | ID, cor = corAR1(form =
~1|ID))
summary(fit_lmeA)

The model recovers phi and b nicely:

Correlation Structure: AR(1)
 Formula: ~1 | ID
 Parameter estimate(s):
      Phi
0.5522512
Fixed effects: y ~ x
                 Value  Std.Error  DF   t-value p-value
(Intercept) -0.1189931 0.18745596 899 -0.634779  0.5257
x            0.3076964 0.01017406 899 30.243235  0.0000


Kevin Wright


On Wed, Sep 16, 2015 at 3:19 AM, Paul Buerkner <paul.buerkner at gmail.com>
wrote:

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

--
Kevin Wright