[R-sig-ME] glmADMB mcmc output

[Dean Castillo]

Ben, Thank you for your response and the code you included. I think
this is what I was looking for.

I have tried the zero-inflated model in another package and the fit
for the zeroinflated negative binomial is better than negative
binomial alone.
I am having trouble fitting the zeroinflated model in the glmmADMB
package and I am getting an error I have seen come up in the list
archives. I have not had time to trouble shoot this just yet.

Dean

On Tue, Mar 13, 2012 at 4:02 PM, Ben Bolker <bbolker at gmail.com> wrote:
> 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))
> }
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models