[R-sig-ME] Covariance between two traits in MCMCglmm

[Jarrod Hadfield]

Hi Ned,

I've made comments in-line below.

Cheers,

Jarrod



Quoting Ned Dochtermann <ned.dochtermann at gmail.com>:

> Jarrod,
> Thank you very much, this clarifies some things both in regards to this
> analysis and some other projects in which I'm calculating covariance
> matrices from repeated measures. Using the flat prior I'm getting reasonable
> agreement among multiple (5) chains and reasonable estimates as well as full
> (2000) effective sample sizes. Unfortunately I've had some trouble with
> getting multiple chain results from overdispersed starting values
> (QUASI=TRUE seems to produce values that don't work for these data, or
> that's my guess for basis of the error). Once I get that taken care of
> hopefully it will be reasonably simple to get estimates and credibility
> ranges from across chains.

When you say "QUASI=TRUE seems to produce values that don't work" what  
do you actually mean - that is fails with an error, or that the chain  
doesn't converge?  If the latter I would worry - perhaps you could use  
plot(mcmc.list()) to show us the traces?




>
> Regarding the degree of informativeness of the prior, I guess I
> misunderstood a comment in the "How to" paper's supplementary
> material--there's a comment about setting the certainty to the dimensions of
> G which I took to mean something other than it does. As mentioned my
> understanding of prior specification is still a work in progress.

I have asked that the MCMCglmm section of Wilson's "How to" paper be  
rectified or removed and I will hassle them again to do so.

>
> Thanks for your help, it is greatly appreciated.
> Ned
>
> --
> Ned Dochtermann
> Department of Biology
> University of Nevada, Reno
>
> ned.dochtermann at gmail.com
> http://wolfweb.unr.edu/homepage/mpeacock/Ned.Dochtermann/
> http://www.researcherid.com/rid/A-7146-2010
> --
>
>
>
> -----Original Message-----
> From: Jarrod Hadfield [mailto:j.hadfield at ed.ac.uk]
> Sent: Thursday, March 24, 2011 4:41 AM
> To: Ned Dochtermann
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] Covariance between two traits in MCMCglmm
>
> Hi Ned,
>
> You should probably interact your fixed effects with trait, because at
> the moment you are assuming things like Sex have the same effect on
> P.TimeFront and M.TimeFront  - I'm not sure if this in intentional?
>
> Regarding priors - the ones you have used are quite informative with
> respect to some parameters. For example, specifying  nu=2 as a prior
> for a 2x2 covariance matrix means that the marginal prior for each of
> the variances is inverse-Wishart with nu=1. This can have a strong
> effect if there are not much data.
>
>  From experience I find
>
> list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*a))
>
> where a is something large (e.g. 1000, depending on the scale of the
> data) works well for the two standard deviations and the correlation,
> in terms of informativeness.
>
> You can't use parameter expanded priors for the residual term yet, so
> you will have to stick with the standard inverse-Wishart (or use
> another program). Generally, data are highly informative for the
> residual part so often the posterior is not very sensitive to the
> prior specification. Nevertheless, you should check alternatives:
>
> V=diag(2), nu=1.002 gives the inverse-gamma prior for the variances
> with shape=scale=0.001
> V=diag(2)*1e-6, nu=3 is flat for the correlation from -1 to 1
>
> Cheers,
>
> Jarrod
>
>
>
>
> On 23 Mar 2011, at 19:59, Ned Dochtermann wrote:
>
>> List members,
>>
>> I'm currently trying to estimate the covariance and correlation
>> between two
>> variables for which I have repeated measurements at the level of
>> individuals. I'm getting output but would appreciate feedback
>> regarding my
>> coding, if I've learned anything it is that getting output isn't the
>> same as
>> getting output because you've done things correctly!
>>
>> BACKGROUND:
>> I had conducted univariate analyses using glmer and the binomial
>> family, the
>> variables are the amount of time spent essentially in one area versus
>> another (based on some other validation work these variables equate to
>> tolerance of predation risk by individuals and aggressiveness towards
>> conspecifics). I then extracted the BLUPS from glmer, recognizing
>> that they
>> were neither unbiased nor linear in this case, and calculated
>> Spearman's
>> correlation (rs ~ 0.4). For well known reasons this is not a
>> preferred way
>> to go so I wanted to calculate the correlation directly from a
>> multi-response model. I know you can do this in glmer by including
>> one of
>> the two variables as a covariate but I wanted to do it from a single
>> model.
>>
>> MCMCglmm MODEL:
>> I chose a generally uninformative prior because I don't have any a
>> priori
>> expectations for the variances but used a lengthy burnin. However,
>> to be
>> honest, I don't have a good grasp of priors despite having read all
>> the
>> MCMCglmm documentation and some other readings besides those.
>>
>> "mass.bar" and "mass.dev" represent mass, obviously, but split to
>> distinguish between within versus among individual effects. I'm not
>> actually
>> interested in the fixed effects from a multi-response model though,
>> just the
>> correlation. Autocorrelation and other diagnostics look ok. Anyway,
>> the
>> code:
>>
>> ####
>> multi.prior<-
>> list(R=list(V=diag(2),nu=2),G=list(G1=(list(V=diag(2),nu=2))))
>>
>> corr.mcmc<-
>> MCMCglmm(cbind(P.TimeFront,P.TimeBack,M.TimeFront,M.TimeBack)
>> ~(trait-1)+Site+Sex+mass.bar+mass.dev, random=~us(trait):Sub_ID,
>> rcov=~us(trait):units,
>> family=c("multinomial2","multinomial2"),
>> data=Compiled.3,
>> nitt=1080000,thin=480,burnin=120000,prior=multi.prior, verbose=FALSE)
>>
>> xx<-matrix(c(posterior.mode(corr.mcmc$VCV)
>> [1],posterior.mode(corr.mcmc$VCV)[
>> 2],
>> posterior.mode(corr.mcmc$VCV)[3],posterior.mode(corr.mcmc$VCV)[4]),2)
>> cov2cor(xx)
>> ####
>>
>> This produces an estimate (0.9) of the correlation between
>> P.TimeFront and
>> M.TimeFront (TimeBack's are linearly dependent on TimeFronts) which
>> is the
>> value I want. Unfortunately I can't find an example of a multi-
>> response
>> binomial model in the documentation or in the archives of this
>> listserv so
>> I'm a bit concerned that my specification of the response variables in
>> particular, or other model terms, is incorrect. Any feedback on this
>> would
>> be appreciated.
>>
>> Thanks a lot,
>> Ned
>>
>> --
>> Ned Dochtermann
>> Department of Biology
>> University of Nevada, Reno
>>
>> ned.dochtermann at gmail.com
>> http://wolfweb.unr.edu/homepage/mpeacock/Ned.Dochtermann/
>> http://www.researcherid.com/rid/A-7146-2010
>> --
>>
>> _______________________________________________
>> 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.
>
>
>
>



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