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

Thu Aug 25 17:54:21 CEST 2011

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.
>
>