[R-sig-ME] MCMCglmm with a AR(1)

[Kim, S.]

Dear list member,

I was wondering if anyone would able to give me a guidance with the following:

I would like to consider both random effects and temporal autocorrelation in my model. So I have used MCMCglmm (developed by Jarrod Hadfield).
But I have a difficulty with considering temporal autocorrelation (in my case, first-order autoregressive process, AR(1)).

At the most of previous studies, the autocorrelation is evaluated on the covariance structure of residuals.

Corr(eit, eit')~rho^|t-t'|
where, e is the residual, i is the individual, and t is the time.

However, in my multinomial response, this is not the case because the residuals do not follow the normal distribution anymore in logistic regression model (So I set R-structure as fix=1).
I have a feeling that the autocorrelation should be evaluated on the realization of response rather than residuals, therefore,

Corr(yit, yit')~rho^|t-t'|
where, e is the residual, i is the individual, and t is the time.

Could anyone help me to take temporal autocorrelation into account in my model?

Belows are the short version of the data.
I have a longitudinal data with repeated measures. Observation spans several months with the choice of individuals out of 5 choice alternatives (A,B,C,D,E).

ind.id

P1

P2

T1

T2

T3

T4

date.time

response

1

3

2

410

2.8

0

3

1/26/2017 9:35

A

1

3

2

603

11.9

0

1

1/26/2017 9:48

C

2

1

3

410

13.2

0

3

1/25/2017 12:28

B

2

1

3

639

8.9

0

1

1/25/2017 12:50

C

2

1

3

538

11.8

50.5

3

1/25/2017 17:20

B

2

1

3

628

14

0

3

1/25/2017 18:08

B

2

1

3

239

7.9

0

3

1/25/2017 18:25

B

2

1

3

561

12.7

50.5

2

1/25/2017 18:50

D

2

1

3

150

9.4

50.5

2

1/26/2017 13:14

D

3

1

3

476

13

50.5

3

1/26/2017 18:45

B

3

1

3

1146

12

0

3

1/26/2017 21:14

E

3

1

3

1038

11.5

50.5

1

1/26/2017 21:37

D

3

1

3

2230

1.8

18

3

1/27/2017 8:06

B

4

1

3

1916

6.2

18

2

1/27/2017 17:08

C

4

1

3

603

10.6

50.5

1

1/27/2017 19:29

D

4

1

3

627

3.6

0.7

3

1/27/2017 20:52

B

4

1

3

918

2.2

0

3

1/28/2017 9:16

B

5

1

3

382

9.6

50.5

3

1/28/2017 9:56

B

5

1

3

4940

0.9

0.7

3

1/29/2017 9:53

E

5

1

3

383

10.6

50.5

3

1/29/2017 12:14

B

5

1

3

232

8.1

0

3

1/29/2017 12:22

A

5

1

3

302

5.4

50.5

3

1/29/2017 12:30

B



The code below is the model without AR(1) but with random effects grouped by individual.

# Set variable type
sig$ind.id<- as.factor(sig$ind.id)
sig$P1<- as.factor(sig$P1)
sig$P2<- as.factor(sig$P2)
sig$T1<- as.factor(sig$T1)
sig$T2<- as.factor(sig$T2)
sig$T3<- as.factor(sig$T3)
sig$T4<- as.factor(sig$T4)
sig$response<- as.factor(sig$response)

library(MCMCglmm)

# Set priors
k <- length(levels(sig$response))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))

prior <- list(
                  G = list(G1 = list(V = diag(k-1), n = k-1)),
                  R = list(fix=1, V= (1/2) * (I + J), n = k-1))

# Fit MCMCglmm
m <- MCMCglmm(fixed = response ~ P1+P2+T1+T2+T3+T4+trait-1,
                                     random = ~ us(trait):ind.id,
                                     rcov = ~us(trait):units,
                                     prior = prior,
                                     burnin =10000,
                                     nitt = 110000,
                                    thin = 10,
                                    pr = TRUE,
                                     family = "categorical",
                                     data = sig,
                                     verbose = T)

Many thanks in advance.
Seheon


	[[alternative HTML version deleted]]