[R-sig-ME] MCMCglmm_2.13

[Malcolm Fairbrother]

Hi Jarrod,

Thanks--that makes sense.

I'm not familiar with ASReml, but it would be difficult, if not impossible, to use nlme to fit the kind of model I have in mind. Bracketing the issue of estimating the variance structure parameters, the problem with nlme is that it can't handle distances of zero in the covariance structure. But where data consisting of repeated observations are located on a spatial grid that doesn't change over time, there will be zero distances for the spatial covariances (i.e., any observation location will have zero spatial distance from itself at another time).

The approach below, which I think takes into account your responses to my e-mails, allows for *both* temporal and spatial autocorrelation.

Cheers,
Malcolm


library(nlme); library(MCMCglmm)
cs1AR1 <- corAR1(0.8, form = ~ Time | Rat) # value here could be estimated roughly using a semivariogram
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)

smat <- band(solve(corMatrix(cs1AR1)[[1]]), -1, 1)
invC <- as(kronecker(diag(nlevels(BodyWeight$Rat)), smat),  "CsparseMatrix")
rownames(invC) <- 1:dim(BodyWeight)[1]
BodyWeight$timeRat <- as.factor(1:dim(BodyWeight)[1])

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

# add some simulated spatial structure to the data (imagining these were not rats but, say, spatial observation points):
BodyWeight$weight.xy <- BodyWeight$weight + matrix(rnorm(nlevels(BodyWeight$Rat), sd=25) %*% chol(corMatrix(cE)[[1]]), ncol=1)[BodyWeight$Rat]

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)))
MC3a <- MCMCglmm(weight.xy ~ as.numeric(Time), random= ~ us(1+Time):Rat + timeRat, ginverse=list(timeRat=invC), data=BodyWeight, prior=prior)
MC3b <- MCMCglmm(weight.xy ~ as.numeric(Time), random= ~ us(1+Time):Rat + timeRat, ginverse=list(timeRat=invC, Rat2=smat2), data=BodyWeight, prior=prior)

# MC3b treats the spatial covariance between observations at different times as the same as observations at the same time
# I can imagine objections to this, but I think that conditional on everything else this should in many instances be OK





On 29 Jul 2011, at 09:33, Jarrod Hadfield wrote:

> Hi Malcolm,
> 
> You can combine variance structures as you have done, but given some of the parameters of the variance structure (e.g. phi) have to be known, I think you're better off fitting these models in lme or ASReml where they are estimated.
> 
> The first Rat example using Time2 is not equivalent to the model I set up (or your lme model m3). In the original models there are covariances between observations within Rats, whereas in your specification you have covariances between observations within AND between Rats (i.e. You have COV(u_{11},u_{12}) =  COV(u_{11},u_{22}), where the first subscript is Rat and the second Time2, whereas before it COV(u_{11},u_{22})=0).
> 
> I'll take a look into the problem you were getting when time was used as a fixed covariate and a factor associated with a ginverse.
> 
> Cheers,
> 
> Jarrod
> 
> 
> 
> 
> 
> Quoting Malcolm Fairbrother <m.fairbrother at bristol.ac.uk> on Thu, 28 Jul 2011 22:17:46 +0100:
> 
>> 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:
>>>> 
>> 
>> 
>> 
> 
> 
> 
> -- 
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
> 
>