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