[R-sig-ME] MCMCglmm output 1) on right scale? 2) produces huge deviance spread?

Ryan King c.ryan.king at gmail.com
Sun Jan 15 20:35:40 CET 2012


Hi, I have an MCMCglmm run predicting a binary outcome with a single
fixed effect and two groups of random effects like so

testprior4<-list( R=list(V=1, fix=1), G=list(G1=list(V=1, nu=1,
alpha.mu=0, alpha.v=.5^2 ), G2=list(V=1, nu=1, alpha.mu=0, alpha.v=1 )
), B=list(mu=c(0,0), V=diag(2)*4 ) )

testmcmc5b<-MCMCglmm(fixed=myy2~myx,random=~idv(myz)+idv(newz3)   ,
family = "ordinal", prior= testprior4, data=mydata, nitt=6000,
thin=10, burnin=5000 , pr=TRUE)

I can provide the full simulation details, but hope that this is a
general rather than specific problem.

The first set of random effects have large effects, but the second set
is just noise.

I'm setting up a bridge sampler to get a bayes factor, which means I
need to get a log-likelihood + log-prior for my proposal.

I tested my log-likelihood calculation with  MCMC output and got
atrociously spread out log-likelihoods (sd of 15-20). That seems like
rather a lot on the log scale, and means that numerical problems are
likely for the bridge sampler. GLM produces the same deviance up to a
constant (regressing the outcome versus the linear predictor), so that
function is fine. This problem is much smaller when using only the
true predictors (sd = 6) .  The estimated variance component for the
noise predictors never goes above .015, so that they cause this
problem is surprising. Is there a handy explanation / work around?

More troubling, a calibration exercise in the deviance calculation
showed that the groups of linear predictors were not on the right
scale.
That is, multiplying through the two sets of random effects and design
matrices to get batches of linear predictors, and running a glm of the
binary outcome versus the linear predictors using a probit (or logit)
link gets coefficients consistently not 1. The true random effects get
a coefficient of about .65, and the noise random effects between .7
and .4 depending on how I set them up. I don't think it's a bad mixing
problem, because 1) I have gotten the same thing on independent
chains, 2) effectiveSize() reports large effective samples 3) there is
quite a bit of spread in the deviances.  This problem persists if I
run with family="categorical".

When I plug in the re-calibrated linear predictor, the sd of the
log-likelihood goes down a lot, from 15 to 6 in the
two-predictor-group case and from 6 to 3 when using only the true
predictors. That's much better, but still quite a premium for having a
group of noise predictors.

Thanks,
Ryan King
Dept Health Studies
University of Chicago




More information about the R-sig-mixed-models mailing list