[R-sig-ME] R: Re: pMCMC and HPD in MCMCglmm

Jarrod Hadfield j.hadfield at ed.ac.uk
Thu Aug 25 18:17:03 CEST 2011 

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.

More information about the R-sig-mixed-models mailing list