[R-sig-ME] MCMCglmm with a AR(1)
Kim, S.
S.Kim at tue.nl
Tue Aug 15 11:39:51 CEST 2017
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]]
More information about the R-sig-mixed-models
mailing list