[R-sig-ME] prior specification in MCMCglmm

[Celine Teplitsky]

Ned,

thanks a lot, I'll try that and compare results with the different priors.

What I just begin to understand is that nu is not really a degree of  
belief you can play with but a more specific parameter of a  
distribution, right?

All the best

Celine

> There's a mistake in my email below.
> The generalization of an inverse-gamma prior for multi-response models
> should, I think, just be:
> list(V=diag(x), nu=1.002)
>
> Sorry about that,
> Ned
>
>
>
>
> -----Original Message-----
> From: Ned Dochtermann [mailto:ned.dochtermann at gmail.com]
> Sent: Tuesday, April 26, 2011 1:12 PM
> To: 'r-sig-mixed-models at r-project.org'; 'bonamy at horus.ens.fr'
> Subject: Re: [R-sig-ME] prior specification in MCMCglmm
>
> Celine and Pierre,
>
> I too am still very unclear on priors however a month ago Jarrod replied to
> some questions of mine regarding multi-response models with this:
>
> 	"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"
>
> I think that this would generalize to G1=(list(V=diag(x), nu=1.00x)) for x
> response variables for an inverse-gamma prior for G but I'm not entirely
> positive that such is the case. I'm also not positive how to generalize the
> flat prior (mainly would nu stay at three?), the first one I assume
> generalizes to (V=diag(x), nu=x, alpha.mu=c(0,0...0x),alpha.V=diag(x)*a).
>
> The whole thread for this discussion starts here:
> https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q1/005694.html
>
> The whole prior issue is still a mystery to me so I continue to be
> uncomfortable with those approaches where they're necessary; nonetheless for
> some data I'm analyzing they're basically necessary (e.g. multi-response
> generalized mixed models). By and large the univariate analyses I've
> conducted produce highly concordant results regardless of whether fitting
> via REML, ML, Laplace or MCMC and regardless of priors for the latter. This
> gives me a bit of confidence in the quantitative results and quite a bit of
> confidence in the inferences. Unfortunately such comparisons aren't as
> feasible (for me) with multi-response models.
>
> 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
> --
>
> Message: 5
> Date: Tue, 26 Apr 2011 09:41:58 +0200
> From: "Pierre B. de Villemereuil" <bonamy at horus.ens.fr>
> To: Celine Teplitsky <teplitsky at mnhn.fr>
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] prior specification in MCMCglmm
> Message-ID: <4DB67746.4000607 at horus.ens.fr>
> Content-Type: text/plain
>
> Dear Ciline,
>
> I'm not very comfortable with covariance priors, but my guess is that,
> is this case, you've got to really specify a inverse-Wishart as a prior.
>
> You should check into Hadfield's article introducing MCMCglmm, they use
> something like V=diag(dimV)/4 and n=dimV. Why ? I have no idea. If your
> data are not standardized, I don't think you should divide by 4 (but
> then, you specify your prior as if your guess for variance components is
> that they all equal 1), but for the rest...
>
> Sorry I can not help further. Maybe somebody else would be able to help
> on this subject ?
>
> Pierre.
>
>
> Le 26/04/2011 09:17, Celine Teplitsky a icrit :
>> Dear Pierre,
>>
>> thanks a lot! It does help, but I will need time to fully understand
>> the paper.
>>
>> Just one further question if I may, what prior would you use for a
>> covariance then?
>>
>> Many thanks again,
>>
>> All the best
>>
>> Celine
>>
>>> Dear Ciline !
>>>
>>> One usual "non informative" prior on variance component is V=1, and
>>> nu=0.002, which correspond to a inverse-Gamma(0.001,0.001). This is
>>> usual, but that is not to say that it is really non informative.
>>> Indeed, inverse-Gamma(e,e) is weakly informative, since the posterior
>>> can depend on the choice of e.
>>>
>>> Concerning the WAMwiki suggestion to use the phenotypic variance to
>>> set the prior, this quite not orthodox since no information coming
>>> from your current dataset should be used to define the prior (but you
>>> could use previous data to parametrize your prior).
>>>
>>> I would suggest you to refit your model considering V=1 and nu=0.002
>>> as a (so-called) non informative prior. Other solutions exist like
>>> using the parameter expansion and a chi2 distribution by setting
>>> V=1,nu=1000,alpha.mu=0 and alpha.V=1, which is also weakly
>>> informative (it has more weight in variance values less than 10).
>>>
>>> For more information about priors on variance component, and
>>> parameter expansion, it would suggest you to read :
>>>
>>>    1. A. Gelman, + Prior distributions for variance parameters in
>>>       hierarchical models ;, /Bayesian analysis/ 1, n^o . 3 (2006):
>>>       515--533.
>>>
>>> In the hope I'm helping.
>>>
>>> Pierre de Villemereuil.
>>>
>>>
>>> Le 25/04/2011 13:26, Celine Teplitsky a icrit :
>>>> Dear all,
>>>>
>>>> I realise that Jarrod is doing field work, but I'm really hoping
>>>> someone can answer my question while he's not around.
>>>>
>>>> I am running animal models estimating covariances between life
>>>> history traits, and I'm having trouble knowing which prior to use.
>>>>
>>>> Thing is, if I use a prior as described on the Wam wiki site with
>>>> V=PhenotypicVar/4 (as I have 3 random effects + residual), I have
>>>> very nice results, with some significant genetic correlations
>>>> between some life history traits.
>>>>
>>>> However, one reviewer asked about prior sensitivity because CI were
>>>> pretty large, so I went back to MCMCglmm course notes and saw that
>>>> non informative prior were supposed to be V=diag(nbDimV)*0 and
>>>> n=nbDimV-3. This led to an error message about G being ill
>>>> conditioned, so I tried with diag(nbDimV)*0.001 and
>>>> diag(nbDimV)*0.01 instead of diag(nbDimV)*0, and diag(nbDimV)*0.01
>>>> worked... But then I have the posterior of additive genetic variance
>>>> collapsing on 0 for some trait. So my guess would be that I should
>>>> use those latest priors, and believe my nice results did not exist.
>>>> But as Hadfield et al paper and the Wam wiki website do not
>>>> recommend those priors, I am a bit confused. Could someone help me
>>>> figure out what would be the right thing to do?
>>>>
>>>> All my apologies if this is a silly question, but I'm feeling a bit
>>>> lost here
>>>>
>>>> Thanks a lot in advance
>>>>
>>>> Celine
>>>>
>>>
>>
>>
>> --
>>
>> Celine Teplitsky
>> Dipartement Ecologie et Gestion de la Biodiversiti UMR 7204
>> Uniti Conservation des Esphces, Restauration et Suivi des Populations
>> Case Postale 51
>> 55 rue Buffon 75005 Paris
>>
>> Webpage :http://www2.mnhn.fr/cersp/spip.php?rubrique96
>> Fax : (33-1)-4079-3835
>> Phone: (33-1)-4079-3443
>
>
> 	[[alternative HTML version deleted]]
>
>
>
> ------------------------------
>
> Message: 6
> Date: Tue, 26 Apr 2011 11:53:02 +0200
> From: peter dalgaard <PDalgd at gmail.com>
> To: Junqian Gordon Xu <xjqian at gmail.com>
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] interaction term in null hypothesis
> Message-ID: <E1153437-757B-42E9-9A72-89EF9A368101 at gmail.com>
> Content-Type: text/plain; charset="us-ascii"
>
>
> On Apr 26, 2011, at 00:05 , Junqian Gordon Xu wrote:
>
>> I have a quick question for a simple model as below:
>>
>>> Fix + (1 | Rand) + (1 | Rand : Fix)
>>
>> Which one is the null hypothesis:
>>
>>> 1 + (1 | Rand)
>>
>> or
>>
>>> 1 + (1 | Rand) + (1 | Rand : Fix)
>>
>> To me the interaction term (1 | Rand : Fix) does not make much sense if
>> no fixed effect term is present in the model, but I'm not sure.
>>
>
> It does make sense, at least sometimes. For one thing, such interactions are
> often aliased to aspects of the experimental design; e.g., if you have
> randomized different treatment to left side and the right side of test
> subjects, then the random interaction is equivalent to the (random)
> difference between the two sides within the same subject. Also, you could
> conceivably have a kill-or-cure drug with positive effects for some and
> negative effects for others, with the question being whether the total
> effect is positive or negative.
>
>> Regards
>> Gordon
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> --
> Peter Dalgaard
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com
>
>
>
> ------------------------------
>
> _______________________________________________
> R-sig-mixed-models mailing list
> R-sig-mixed-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
> End of R-sig-mixed-models Digest, Vol 52, Issue 47
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>



-- 

Celine Teplitsky
UMR 5173 MNHN-CNRS-P6 'Conservation des espèces, restauration et suivi des
populations'
Muséum National d'Histoire Naturelle
CRBPO, 55, Rue Buffon, CP51, 75005 Paris, France

Fax : (33-1)-4079-3835
Phone: (33-1)-4079-3443