[R-sig-ME] Question about autocorrelation in nlme
Kevin Wright
kw.stat at gmail.com
Wed Sep 16 15:46:21 CEST 2015
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]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
--
Kevin Wright
More information about the R-sig-mixed-models
mailing list