[R-sig-ME] MCMCglmm for binomial models?

[Whit Armstrong]

Bob,

Do you mind posting the winbugs code for your model?  I am curious to
see how pymc performs.  Additionally, I have a pure c++ MCMC package.
Which I would like to try too.

A link to some sample data would be useful as well.

-Whit


On Mon, Sep 20, 2010 at 8:17 AM, Bob Farmer <farmerb at gmail.com> wrote:
> Thanks for the reply.  Unfortunately, I'm not sure that that's the
> whole story.  Upping the iterations to achieve an effective sample
> size of 1000 (output pasted below) still leads to some differences
> between the two models (total runtime on a quad-core Windows box is
> under 2 minutes, in case you want to try).  In fact, running it a
> second time shooting for an n.eff of 2000 (thin = 325) leads to even
> more divergent estimates (e.g. mean.int == -72.96).  Again, assuming
> that the glmer() output in this case is a `gold standard', maybe
> there's another puzzle piece missing?
>
> Thanks again.
> --Bob
>
>
> #################now binomial from a reference dataset
> gm3 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
>              family = binomial, data = cbpp)
> mc3<-MCMCglmm(cbind(incidence, size - incidence) ~ period,
>  random = ~ herd,
>  family = "multinomial2",
>  data = cbpp, verbose = FALSE,
>  nitt = 750E3, thin = 650, burnin = 100E3
> )
> summary(gm3)@coefs
> summary(mc3)
>
> Leads to:
>
>> summary(gm3)@coefs
>              Estimate Std. Error   z value     Pr(>|z|)
> (Intercept) -1.3985351  0.2278906 -6.136871 8.416284e-10
> period2     -0.9923347  0.3053852 -3.249452 1.156274e-03
> period3     -1.1286754  0.3260491 -3.461673 5.368286e-04
> period4     -1.5803739  0.4288037 -3.685542 2.282169e-04
>
> Versus:
>
>> summary(mc3)
>
>  Iterations = 749351
>  Thinning interval  = 100001
>  Sample size  = 1000
> ....
>            post.mean l-95% CI u-95% CI eff.samp pMCMC
> (Intercept)   -1.6102  -2.2092  -0.8357     1000 0.002 **
> period2       -1.2476  -2.2146  -0.1284     1000 0.020 *
> period3       -1.3743  -2.3387  -0.2746     1000 0.008 **
> period4       -1.9677  -3.2814  -0.8065     1000 0.008 **
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>