[R-sig-ME] MCMCglmm with a AR(1)
Kim, S.
S.Kim at tue.nl
Tue Aug 15 14:09:56 CEST 2017
Thanks Jarrod,
Unfortunately the data I have spans a few months with repeated measurements and each individual participates the survey at random time points.
Yes. I would consider a vector autoregressive (VAR(p)).
And sorry for the messy data. The example data and code looks like follows. Also it would be greatly appreciated if you comment whether there is any mistakes in the code (without AR(1)).
> head(sig)
ind.id P1 P2 T1 T2 T3 T4 date.time response
1 1 3 2 410 2.8 0 3 1/26/2017 9:35 A
2 1 3 2 603 11.9 0 1 1/26/2017 9:48 C
3 2 1 3 410 13.2 0 3 1/25/2017 12:28 B
4 2 1 3 639 8.9 0 1 1/25/2017 12:50 C
5 2 1 3 538 11.8 50.5 3 1/25/2017 17:20 E
6 2 1 3 628 14 0 3 1/25/2017 18:08 D
> 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= diag(k-1), 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)
Best regards,
Seheon
-----Original Message-----
From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Jarrod Hadfield
Sent: Tuesday, August 15, 2017 12:33 PM
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] MCMCglmm with a AR(1)
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.
_______________________________________________
R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list