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

[Jarrod Hadfield]

Hi Bob,

Quoting Bob Farmer <farmerb at gmail.com>:

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

There are two ways of handling over-dispersion. If you use a  
multiplicative model the fixed effect estimates should stay the same  
and the standard error should increase (if the data are  
over-dispersed). I am not sure what the state of play is with  glmer's  
quasi models but with my version at least, it does not look right. The  
alternative is to fit an additive model by fitting a level for each  
observation, which is what I did in the last post. You should not  
expect the fixed effect coefficients to remain the same with additive  
models unless they are zero.

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

Could you post the output from

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

I'm on a Mac so I can't run this (unless things have changed?)


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

Further info can be found in section's 1.3, 1.5, 2.7, 3.6 of the  
CourseNotes and Gelman's  2006 Bayesian Analysis 1(3) 515-533 is a  
good reference for parameter expanded priors, at least for single  
variance components. If I get more time I will make the notes more  
complete and indexed!

Cheers,

Jarrod


>
> --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")]
>
> _______________________________________________
> 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.