[R-sig-ME] MCMCglmm_2.13

[Malcolm Fairbrother]

Jarrod,

This is indeed very interesting, and plenty fancy enough. What I'm particularly intrigued by is the possibility of fitting models with *both* time-series and spatial correlation structures. The code below represents my attempt to do that. Does this make sense?


library(nlme); library(MCMCglmm)
cs1AR1 <- corAR1(0.2, form = ~ Time | Rat) # the value chosen here doesn't seem to affect the estimates much
cs1AR1 <- Initialize(cs1AR1, BodyWeight)
corMatrix(cs1AR1)[[1]] # shows the correlation structure for each Rat (all 16 are the same, since each Rat is measured on the same occasions)

# a revised version of what you proposed (which, if I've got this right, proves about 25% faster for this dataset, and is simpler to set up):

smat <- Matrix(band(solve(corMatrix(cs1AR1)[[1]]), -1, 1), sparse=T)
rownames(smat) <- unique(BodyWeight$Time)
BodyWeight$Time2 <- BodyWeight$Time # distinguishing between two different "Time" variables seems to prevent R from crashing when calling MCMCglmm...

prior=list(R=list(V=1, nu=0.02), G=list(G1=list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*1000), G2=list(V=1, nu=1, alpha.mu=0, alpha.V=1000)))

MC0 <- MCMCglmm(weight ~ as.numeric(Time), random= ~ us(1+Time):Rat + Time, data=BodyWeight, prior=prior)
MC1 <- MCMCglmm(weight ~ as.numeric(Time), random= ~ us(1+Time):Rat + Time2, ginverse=list(Time2=smat), data=BodyWeight, prior=prior)
# this approach takes advantage of the fact that the correlation structure is the same for each Rat
# it generates the same results as your m3.mcmc, though it reallocates some of the variance from your "timeRat" level to the units level
# but I think that's not a problem, because the units/residual variance can still be identified using va*mat+vb*I
# make sense?

# anyway, here's the (maybe) interesting bit, combining time-series and spatial autocorrelation:

BodyWeight$x <- sample(1:length(unique(BodyWeight$Rat)))[BodyWeight$Rat]
BodyWeight$y <- sample(1:length(unique(BodyWeight$Rat)))[BodyWeight$Rat]
plot(unique(BodyWeight$x), unique(BodyWeight$y)) # a map

cE <- corExp(form = ~ x + y | Time)
cE <- Initialize(cE, BodyWeight)
corMatrix(cE)[[1]] # a spatial correlation structure (albeit for data without any real spatial pattern)

smat2 <- Matrix(solve(corMatrix(cE)[[1]]), sparse=T) # not actually sparse, unfortunately
rownames(smat2) <- unique(BodyWeight$Rat)
BodyWeight$Rat2 <- BodyWeight$Rat # again, this seems necessary to prevent R from crashing when calling MCMCglmm

MC2 <- MCMCglmm(weight ~ as.numeric(Time), random= ~ us(1+Time):Rat + Time2, ginverse=list(Time2=smat, Rat2=smat2), data=BodyWeight, prior=prior)

From the code MC2 *looks like* it's taking spatial correlation into account... but is it? I'm not sure whether I've got this right, and I'm a bit perplexed about how one might go about simulating data to test this (obviously the data here have no spatial autocorrelation).

Basically, is this an AR1 + spatial corExp model? Your clarifications are always appreciated.

Thanks again,
Malcolm




On 28 Jul 2011, at 10:20, Jarrod Hadfield wrote:

> Hi Malcolm,
> 
> I don't think it's as fancy as you think/hope. The covariance structure has to be seperable into a matrix of prameters (V) and a fixed matrix (A), whereas in many spatial/time-series models A is also a function of some parameter(s). It is flexible in the sense that A can be any positive-definite matrix.
> 
> I've edited your code below to get MCMCglmm to fit the AR1 model (with fixed phi).
> 
> Cheers,
> 
> Jarrod
> 
> 
> On 27 Jul 2011, at 19:18, Malcolm Fairbrother wrote:
> 
>> Dear Jarrod,
>> 
>> Can I please just confirm what you meant here by "covariance structures for the random effects"? If I understand correctly, this technique might be very useful for quite a lot of things--though I'm also trying to think it through.
>> 
>> Specifically, the nlme package for example allows for error structures like AR1, etc., in the *residuals* of a mixed model. Are you saying that MCMCglmm's "ginverse" function allows for similar sorts of structure to be modelled in the *random effects*? Consider:
>>