[R] pMCMC and HPD in MCMCglmm

[Ben Bolker]

m.fenati <at> libero.it <m.fenati <at> libero.it> writes:

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

[snip]
> 
>  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?

  My first suggestion would be to post this sort of question to
the r-sig-mixed-models list instead (I would forward it but I'm
replying via Gmane)

  I haven't been able to find much information on the definition
of the Bayesian p-value, but looking into the code of summary.MCMCglmm you 
can see:

solutions <- cbind(colMeans(object$Sol[, 1:nF, drop = FALSE]),   ## means
        coda::HPDinterval(object$Sol[, 1:nF, drop = FALSE]),     ## HPD
        effectiveSize(object$Sol[, 1:nF, drop = FALSE]), ## eff sample size
    2 * pmax(0.5/dim(object$Sol)[1], 
     pmin(colSums(object$Sol[,1:nF, drop = FALSE] > 0)/dim(object$Sol)[1], 
     1 - colSums(object$Sol[, 1:nF, drop = FALSE] >0)/dim(object$Sol)[1])))

  That is, it's basically the minimum of the proportion of samples that
are on one side or the other of zero.

  I would suggest looking at density plots to see what might be happening:

library(coda)
xyplot(mcmc(modelfit$Sol))