[R-sig-ME] MCMCglmm output 1) on right scale? 2) produces huge deviance spread?
Jarrod Hadfield
j.hadfield at ed.ac.uk
Thu Jan 19 12:29:02 CET 2012
Hi Ryan,
Sorry for not replying earlier - I've been away for two months. I'm
having trouble understanding the issues but my gut feeling (especially
for the second issue) is that they arise because of the error variance
being set to one in MCMCglmm rather than zero as in glm. For example,
id<-gl(100,2) # 100 individuals observed twice
u<-rnorm(100) # individual random effects ~ N(0,1)
lp<-u[id]+rnorm(200) # linear predictor includes error ~ N(0,1)
y<-rbinom(200, 1, plogis(lp)) # response (with logit link)
glm(y~lp, family="binomial")$coef[2]
# regression on linear predictor should have coefficient of 1 on
average as you suggest
glm(y~u[id], family="binomial")$coef[2]
# regression on linear predictor omitting noise has a coefficient < 1
Section 2.5 has ways of obtaining a rescaled linear predictor for the
logit link, but the probit link may be a bit easier to work with in
this respect as you can set the variance of the normal in the cdf
calculation to 2 rather than 1.
Alternatively you could extract the latent variables (pl=TRUE) stored
under Liab, which are the linear predictors including the noise term
Note sure this helps for the first problem though - perhaps you could
provide the simulation code together with comments on what you expect
to happen?
Cheers,
Jarrod
On 15 Jan 2012, at 19:35, Ryan King wrote:
> 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
>
> _______________________________________________
> 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.
More information about the R-sig-mixed-models
mailing list