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

Fri Mar 25 10:28:53 CET 2011

Hi Ned,

That sounds fine  - I thought it was a terminating error or lack of  
convergence. For assessing convergence I usually plot the multiple  
chains using mcmc.list, but if you want something more formal  
gelman.diag in coda implements the test in Gelman and Rubin 1992.

Cheers,

Jarrod

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

> Hi Jarrod,
>
> Thanks for getting back to me so promptly.
> The warning message was:
>
> In MCMCglmm(cbind(P.TimeFront, P.TimeBack, M.TimeFront, M.TimeBack) ~  :
>   good starting values not obtained: using Norm(0,1)
>
> This only shows up with start=list(QUASI=FALSE) which I was trying in order
> to run the Gelman and Rubin 1992 diagnostics as my understanding is that
> those diagnostics assume overdispersed starting values. I'm guessing that
> since I'm dealing with proportional data (the actual values are 0:300) the
> distribution starting values are being pulled from in this case isn't
> appropriate. Looking at the autocorrelation and diagnostic plots associated
> with individual chains using the default starting values suggests no problem
> with convergence but I'm trying to be thorough.
>
> With multiple chains, whether from overdispersed starting values or the
> default heuristic approach am I correct that section 3 of Gelman and Rubin
> (1992, Statistical Science) would be the place to go for estimating the
> correlation and its distribution from multiple chains? Nothing in coda looks
> appropriate. While posterior.cor is great for the individual chains, since
> I'm running many it seems reasonable to use those different chains given
> Gelman and Rubin's cautions.
>
> Regarding the "How to" paper it really came down more to my
> misreading/misinterpreting than anything else.
>
> Thanks again,
> Ned
>
>
>
>
> -----Original Message-----
> From: Jarrod Hadfield [mailto:j.hadfield at ed.ac.uk]
> Sent: Thursday, March 24, 2011 4:04 PM
> 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,
>
> 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.
>
>
>
>
>

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