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

[Jarrod Hadfield]

Hi,

MCMCglmm cannot fit AR1 models but it can fit 1st-order antedependence 
models (an ar1 process sampled from time t=0 rather than from the 
stationary process). That being said, it only really works when you have 
short-time series where all observations have been sampled at the same 
time. I'm not sure if this is your case? Its also not clear how you 
would set up the AR1 for a 5-category multinomial  - something like a 
vector autoregressive might make sense but this is not implemented in 
MCMCglmm except for Gaussian traits and at the observation-level rather 
than the residual level. For non-Gaussian traits, particularly 
multinomial, you're probably going to have to code the model yourself in 
stan/JAGS if you want autocorrelation at the observation-level.

Cheers,

Jarrod



(except for Gaussian traits and at the observation-level).


On 15/08/2017 10:39, Kim, S. wrote:
> 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]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.