[R-sig-ME] MCMCglmm_2.13

[Jarrod Hadfield]

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.