[R-sig-ME] glmADMB mcmc output

[Ben Bolker]

Dean Castillo <dmcastil at ...> writes:

> I have run the ML and mcmc for my model. I am wondering if there is a
> simple transformation between the the ML coefficients and mcmc
> coefficients. Or in general how do I interpret the mcmc coefficients
> (I know how to interpret the ml coefficients). For example e^4.03 is
> not in the range of the HPD for beta1 of the mcmc
> 
> nb1<-glmmadmb(D2_eggs~D2_female*(Bacteria+Enviro)+(1|Pop), data=data,
> family="nbinom1")
> > summary(nb1)
> 
> Call:
> glmmadmb(formula = D2_eggs ~ D2_female * (Bacteria + Enviro) +
>   (1 | Pop), data = data, family = "nbinom1")
> 

  [snip]

> Negative binomial dispersion parameter: 243.68 (std. err.: 48.422)

  (This very large NB dispersion parameter is a little bit worrying.
Have you examined your data for weirdness/zero-inflation etc.?  Do
you results look reasonable?)


> >nb_mcmc<-glmmadmb(D2_eggs~D2_female*(Bacteria+Enviro)+(1|Pop), data=data,
> family="nbinom1",mcmc=T, mcmc.opts=mcmcControl(mcmc=10000))
> > m<-as.mcmc(nb_mcmc$mcmc)
> > head(HPDinterval(m))
>           lower     upper
> beta.1 50.4190302 55.231164
> beta.2  6.3920598  9.033473
> beta.3  0.5165269  4.145737
> beta.4  1.7327276  5.167999
> beta.5  1.6289002  3.791544
> beta.6 -3.9886068 -2.057463

  I haven't fully implemented this yet, but here is an embryonic
function for transforming an mcmc object (m) based on a glmmADMB fit
(fit).  At the moment it only does some of the easy stuff
(transforming fixed-effect parameters to their original scale, and
variance parameters to the standard deviation scale)

  There will be a bit more detail about this in the vignette in
the next release.

mcmc_transform <- function(m,fit) {
  if (missing(fit)) {
    fit0 <- fit
    m <- fit$mcmc
    fit <- fit0
  }
  if (!is(m,"mcmc")) stop("m must be an 'mcmc' object")
  if (!is(fit,"glmmadmb")) stop("fit must a 'glmmadmb' object")
  ## zero-inflation
  pz <- m[,"pz",drop=FALSE]
  t_pz <- pz  ## (not transformed)
  ## fixed effects
  fixed <- m[,grep("^beta",colnames(m)),drop=FALSE]
  t_fixed <- as.mcmc(fixed %*% fit$phi)
  colnames(t_fixed) <- names(fixef(fit))
  ## variance parameters: log std dev
  theta <- m[,grep("^tmpL",colnames(m)),drop=FALSE]
  t_theta <- exp(theta)
  ## corr parameters ("offdiagonal elements of cholesky-factor of correlation
matrix")
  corr <- m[,grep("^tmpL1",colnames(m)),drop=FALSE]
  t_corr <- corr
  ## scale/overdispersion parameter
  logalpha <- m[,grep("^log_alpha",colnames(m)),drop=FALSE]
  t_alpha <- matrix(exp(logalpha),dimnames=list(NULL,"alpha"))
  ## random effects
  re <- m[,grep("^u\\.[0-9]+",colnames(m)),drop=FALSE]
  t_re <- re
  mcmc(cbind(t_pz,t_fixed,t_theta,t_corr,t_alpha,t_re),
       start=start(m),end=end(m),thin=frequency(m))
}