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

[Bob Farmer]

seems I replied to Dr. Hadfield without adding the mailing list!
Sorry about that.  See below.
-------------------

Thank you very much, Drs. Hadfield and Onkelinx, for your comments.
I've refit the glmer() model with a quasibinomial link, but it was my
understanding that overdispersion changes the estimated error ranges
of the parameters, and not the values of the parameters themselves
(Gelman and Hill, "Data Analysis Using Regression and
Multilevel/Hierarchical Models 2007, p 115); this is borne out by a
comparison of the equivalent glmer() fits with the data here.  So I
don't see how accounting for overdispersion could change the parameter
estimates, just the inferences arising from them.

I've refit the MCMCglmm model using the priors supplied, and while I
don't see a fantastic alignment with glmer()'s parameter values, there
is good congruence with the bugs() (package "arm") equivalent (at
least, as I've coded it (below), albeit with low n.eff values, and
almost-adequate convergence statistics (Rhat <= 1.1 in most
cases)...).  On the whole, this third hat in the ring suggests to me
that the fits of all 3 model approaches must all be reasonably close
to what the data can say.  So many thanks again for the help --
although I welcome further thoughts.

My only lingering request for Dr. Hadfield is a more systematic
explanation of the priors arguments in the MCMCglmm helpfile -- even
in Chapter 8 of the course notes, the list structure and parameter
values seem to be assumed knowledge, rather than provided as an
exhaustive list of options, along with descriptions of each option
(e.g. alpha.V -- what is it explicitly?).  Presently, I can't derived
these priors myself, because I don't know their context well enough to
make informed choices -- but I would prefer a bit more up-front
context to use these existing prior sets with confidence.
Nonetheless, MCMCglmm seems like a wonderful addition to the set of
free modeling tools available!

--Bob

Code:

#################now binomial from a reference dataset
rm(list=ls())
library(MCMCglmm); library(lme4)
####lmer version
(gm3 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
             family = quasibinomial, data = cbpp))

####MCMCglmm version
#priors in email from Jarrod Hadfield
 prior2=list(R=list(V=1, nu=0.002), G=list(G1=list(V=1, nu=1,
alpha.mu=0, alpha.V=100)))
st<-Sys.time()
mc3<-MCMCglmm(cbind(incidence, size - incidence) ~ period,
 random = ~ herd,
 family = "multinomial2",
 data = cbpp, verbose = FALSE,
 nitt = 750E3, thin = 650, burnin = 100E3,
 prior = prior2
)
runt<-Sys.time() - st
runt
summary(gm3)@coefs
summary(mc3)

####now a WinBUGS equivalent
library(arm)
n<-nrow(cbpp)
n.herd<-length(unique(cbpp$herd))
model.data<-list(
 "incidence" = cbpp$incidence,
 "size" = cbpp$size,
 "period2" = as.numeric(model.matrix(~cbpp$period - 1)[,2]),
 "period3" = as.numeric(model.matrix(~cbpp$period - 1)[,3]),
 "period4" = as.numeric(model.matrix(~cbpp$period - 1)[,4]),
 "herd" = as.numeric(cbpp$herd),
 "n" = n,
 "n.herd" = n.herd
)
model.inits<-function(){
 list(
   B.0 = rnorm(1),
   B.period2 = rnorm(1),
   B.period3 = rnorm(1),
   B.period4 = rnorm(1),
   b.herd = rnorm(15),
   overdisp = rnorm(56),
   sigma.overdisp = runif(1),
   sigma.b.herd = runif(1)
 )
}
params<-c("B.0", "B.period2", "B.period3", "B.period4",
"sigma.overdisp", "sigma.b.herd")
model<-function(){
 for(i in 1:n){
   incidence[i] ~ dbin(phi2[i], size[i])
   logit(phi[i]) <- B.0 + B.period2*period2[i] + B.period3*period3[i]
+ B.period4*period4[i] +
     b.herd[herd[i]] + overdisp[i]
   phi2[i] <- max(0.00001, min(phi[i], 0.9999999))
   overdisp[i] ~ dnorm(0,tau.overdisp)
 }
 B.0 ~ dnorm(0,0.001)
 B.period2 ~ dnorm(0, 0.001)
 B.period3 ~ dnorm(0, 0.001)
 B.period4 ~ dnorm(0, 0.001)
 tau.overdisp <- pow(sigma.overdisp, -2)
 sigma.overdisp ~ dunif(0, 1000)
 for(j in 1:n.herd){
   b.herd[j] ~ dnorm(0, tau.b.herd)
 }
 tau.b.herd <- pow(sigma.b.herd, -2)
 sigma.b.herd ~ dunif(0, 100)
}
write.model(model, con="mcmcTest.bug")
starT<-Sys.time()
bug3<-bugs(model.data, model.inits, params, "mcmcTest.bug",
 n.iter = 500E3, debug = FALSE)
runt<-Sys.time() - starT
runt

summary(gm3)@coefs
summary(mc3)
bug3$summary[which(rownames(bug3$summary) %in% params), c("mean",
"sd", "2.5%", "97.5%", "Rhat", "n.eff")]