[R-sig-ME] pMCMC and HPD in MCMCglmm
Jarrod Hadfield
j.hadfield at ed.ac.uk
Sat Aug 27 14:20:18 CEST 2011
Hi,
Without prior knowledge of the total variability you have to guess, or
if you can live with the guilt use the posterior estimate from a
preliminary model. In practice, I think the posterior distribution of
the fixed effects are likely to change little with reasonable choices
for the prior unless there is complete or near-complete separation.
Even then, with large data sets differences may be small - but check!
Cheers,
Jarrod
Quoting "m.fenati at libero.it" <m.fenati at libero.it> on Fri, 26 Aug 2011
19:58:22 +0200 (CEST):
> Hi,
> when I set the fixed effect priors for the first example with repeated
> measures I follow the indication of both the coursenote and past
> post (https:
> //stat.ethz.ch/pipermail/r-sig-mixed-models/2010q3/004415.html):
>
> V=diag(n)*((varUNITS+varCLASS)+pi^2/3)
>
> where
> n= number of fixed effects
> varUNITS= residual variance (var of e)
> varCLASS = variance of random effect (var of Zu)
>
> You suggest to assume 1 for the residual variance (varUNITS) + ~2 for the
> random variable variance (varCLASS).
>
> How do you decide to set ~2 for varCLASS? is there an indication about a
> possible range of this value?
>
> Thank a lot
>
> Massimo
>
> -----------------------
> Massimo Fenati
> DVM
> Padova - Italy
>
>
>
>
>
>
>
>> ----Messaggio originale----
>> Da: j.hadfield at ed.ac.uk
>> Data: 25/08/2011 18.17
>> A: "m.fenati at libero.it"<m.fenati at libero.it>
>> Cc: <r-sig-mixed-models at r-project.org>
>> Ogg: Re: R: Re: [R-sig-ME] pMCMC and HPD in MCMCglmm
>>
>> Hi,
>>
>> You have specified a strong prior correlation (~0.5) between the
>> intercept and sex effect:
>>
>>> prior$B
>> $mu
>> [1] 0 0
>>
>> $V
>> [,1] [,2]
>> [1,] 6.289868 3.289868
>> [2,] 3.289868 6.289868
>>
>> Since the data have 50 0's and one 1 the data give a lot of weight to
>> the intercept being negative. As you believe, the data are not very
>> informative about sex-differences in this example but with your prior
>> specification we expect a priori that the intercept and sexM effect
>> are positively correlated. Hence, the "signifcant" negative sex effect.
>>
>> Perhaps you had mistyped the prior specification, and had intended
>>
>> V=diag(2)*(1+pi^2/3) rather than V=diag(2)+pi^2/3 ?
>>
>> as in Section 2.6 of the CourseNotes. Treating the effects as
>> independent in the prior gives results closer to what you would hope
>> for.
>>
>> Jarrod
>>
>>
>>
>> Quoting "m.fenati at libero.it" <m.fenati at libero.it> on Thu, 25 Aug 2011
>> 17:54:21 +0200 (CEST):
>>
>>> Hi Jarrod,
>>> In the past example, where HPD and pMCMC were slightly different, I tested
> an
>>> extreme dataset: 1 positive event on 51 sample of 34 animals. In this
>>> circumstance, even assuming not repeated data, the posterior distribution
> of
>>> the sex beta coefficient MCMC estimates seems to suggest a possible effect
> of
>>> “sex” on the response “dis”. But if I perform the analysis on the
>>> same dataset
>>> under frequentist approach this fails (using glm for perfect separation) or
>>> returns with high p-value (using aalysis of frequency table via fisher
> exact
>>> test). See the following example:
>>>
>>> sex<-c(rep("F",21),rep("M",30))
>>> dis<-c(1,rep(0,50))
>>> dat<-data.frame(sex,dis)
>>> prior<-list(R=list(V=1,fix=1),G=list(G1=list(V=1,nu=0.002)),B=list(mu=c(rep
>>> (0,2)),V=diag(2)*3+pi^2/3))
>>> m.1<-MCMCglmm(dis~sex,slice=T,prior=priorS,data=dat,nitt=800000,thin=100,
>>> burnin=250000,family="categorical",verbose=FALSE)
>>> summary(m.1)
>>>
>>> fisher.test(dat$dis,dat$sex)
>>>
>>> summary(glm(dis~sex,data=dat,family=binomial))
>>>
>>>
>>> How can I interpret the differences between Bayesian (MCMCglmm) and
>>> Frequentist approaches in these circumstances?
>>>
>>> Sorry for the basic question, but I am new in Bayesian world!
>>>
>>> Thanks
>>>
>>> Massimo
>>>
>>>
>>>
>>>
>>> -----------------------
>>> Massimo Fenati
>>> DVM
>>> Padova - Italy
>>>
>>>
>>>
>>>
>>>> ----Messaggio originale----
>>>> Da: j.hadfield at ed.ac.uk
>>>> Data: 24/08/2011 18.11
>>>> A: "m.fenati at libero.it"<m.fenati at libero.it>
>>>> Cc: <r-sig-mixed-models at r-project.org>
>>>> Ogg: Re: [R-sig-ME] pMCMC and HPD in MCMCglmm
>>>>
>>>> Hi,
>>>>
>>>> pMCMC is the two times the smaller of the two quantities: MCMC
>>>> estimates of i) the probability that a<0 or ii) the probability that
>>>> a>0, where a is the parameter value. Its not a p-value as such, and
>>>> better ways of obtaining Bayesian p-values exist.
>>>>
>>>> HPDinterval finds the closest points (c and d) for which Fa(d)-Fa(c) =
>>>> 0.95 (If prob=0.95 in HPDinterval) and Fa is the empirical cumulative
>>>> distribution of a.
>>>>
>>>> Cheers,
>>>>
>>>> Jarrod
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> On 24 Aug 2011, at 16:14, m.fenati at libero.it wrote:
>>>>
>>>>> Hi Jarrod,
>>>>> thanks for your answer, but I have again a lot of confusion. If
>>>>> possible,
>>>>> could you explain to me the definition of pMCMC?
>>>>> Maybe, knowing the right definition of pMCMC I will be able to
>>>>> understand
>>>>> completely your answer.
>>>>>
>>>>> Thank a lot!
>>>>>
>>>>> Massimo
>>>>>
>>>>> -----------------------
>>>>> Massimo Fenati
>>>>> DVM
>>>>> Padova - Italy
>>>>>
>>>>>
>>>>>
>>>>>> ----Messaggio originale----
>>>>>> Da: j.hadfield at ed.ac.uk
>>>>>> Data: 24/08/2011 13.24
>>>>>> A: "m.fenati at libero.it"<m.fenati at libero.it>
>>>>>> Cc: <ndjido at gmail.com>, <r-sig-mixed-models at r-project.org>
>>>>>> Ogg: Re: [R-sig-ME] pMCMC and HPD in MCMCglmm
>>>>>>
>>>>>> Hi Massimo,
>>>>>>
>>>>>> They only need to be slightly skewed (even up to Monte Carlo error
>>>>>> probably) - conclusions drawn from HPDinterval and pMCMC are in
>>>>>> reality almost identical in your example, it is the consequences of
>>>>>> the (arbitrary) distinction between <0.05 and >0.05 that makes them
>>>>>> "feel" different. Lets say we used the cutoff <0.06 and >0.06. Does
>>>>>> HPDinterval(m1$Sol[,3], prob=0.94) overlap zero? If not then
>>>>>> HPDinterval and pMCMC "agree" with respect to which side of the
>>>>>> cutoff
>>>>>> the probability lies ? It may make us happier, but it shouldn't.
>>>>>>
>>>>>> Jarrod
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> On 24 Aug 2011, at 11:45, m.fenati at libero.it wrote:
>>>>>>
>>>>>>> The posterior distribution seem to be only slightly skewed.
>>>>>>> However the question remains: what is the sense of the discrepancy
>>>>>>> between HPD
>>>>>>> and pMCMC?
>>>>>>>
>>>>>>> Thanks
>>>>>>>
>>>>>>> Massimo
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> ----Messaggio originale----
>>>>>>> Da: ndjido at gmail.com
>>>>>>> Data: 24/08/2011 11.43
>>>>>>> A: "m.fenati at libero.it"<m.fenati at libero.it>
>>>>>>> Cc: <r-sig-mixed-models at r-project.org>
>>>>>>> Ogg: Re: [R-sig-ME] pMCMC and HPD in MCMCglmm
>>>>>>>
>>>>>>> Check your posterior distributions, the one corresponding to GENDER
>>>>>>> seems to
>>>>>>> be skewed.
>>>>>>> Ardo.
>>>>>>> On Wed, Aug 24, 2011 at 11:33 AM, m.fenati at libero.it <m.fenati at libero.
> it
>>>>>>>>
>>>>>>> wrote:
>>>>>>> As suggested by Ben Bolker, I re-post the following question in this
>>>>>>> list.
>>>>>>> Thanks
>>>>>>>
>>>>>>>> Dear R users,
>>>>>>>> I’d like to pose aquestion about pMCMC and HDP.
>>>>>>>> I have performed a mixed logistic regression by MCMCglmm (a very
>>>>>>>> good
>>>>>>> package)
>>>>>>>> obtaining the following results:
>>>>>>>>
>>>>>>>> Iterations = 250001:799901
>>>>>>>> Thinning interval = 100
>>>>>>>> Sample size = 5500
>>>>>>>>
>>>>>>>> DIC: 10.17416
>>>>>>>>
>>>>>>>> G-structure: ~ID_an
>>>>>>>>
>>>>>>>> post.mean l-95% CI u-95% CIeff.samp
>>>>>>>> ID_an 0.7023 0.0001367 3.678 2126
>>>>>>>>
>>>>>>>> R-structure: ~units
>>>>>>>>
>>>>>>>> post.mean l-95% CIu-95% CI eff.samp
>>>>>>>> units 1 1 1 0
>>>>>>>>
>>>>>>>> Location effects: febbreq~ as.factor(sex)
>>>>>>>>
>>>>>>>> post.mean l-95% CIu-95% CI eff.samp pMCMC
>>>>>>>> (Intercept) -3.6332 -5.6136 -1.7719 3045 <2e-04 ***
>>>>>>>> as.factor(sex)M -2.9959 -6.0690 0.1969 2628 0.0455 *
>>>>>>>> ---
>>>>>>>> Signif. codes: 0 ‘***’0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>>>>>>>
>>>>>>>>
>>>>>>>> As you can see, pMCMC for gender is just less than 5%, but the
>>>>>>>> credible
>>>>>>>> interval (HPD) is wide and includes the 0 value.
>>>>>>>> How can I interpret these different results?
>>>>>>>>
>>>>>>>> Thank you in advance
>>>>>>>>
>>>>>>>> Massimo
>>>>>>>>
>>>>>>>> -----------------------
>>>>>>>> Massimo Fenati
>>>>>>>> DVM
>>>>>>>> Padova - Italy
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> R-sig-mixed-models at r-project.org mailing list
>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> 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.
>>
>>
>>
>
>
>
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list