[R-sig-ME] Question about autocorrelation in nlme

[Paul Buerkner]

Thank you both for your helpful comments and suggestions!

This indeed answers may question thoroughly! :-)

best
 paul

2015-09-16 16:04 GMT+02:00 Dan Wright <dbrookswr at gmail.com>:

> 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]]
>> >
>> > _______________________________________________
>> > R-sig-mixed-models at r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>>
>> --
>> Kevin Wright
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
>

	[[alternative HTML version deleted]]