[prev in list] [next in list] [prev in thread] [next in thread] 

List:       r-sig-mixed-models
Subject:    Re: [R-sig-ME] simulate for glmmADMB
From:       Ben Bolker <bbolker () gmail ! com>
Date:       2016-08-10 21:01:31
Message-ID: eb5167d3-0fdb-e461-076d-bc2b083bedf0 () mcmaster ! ca
[Download RAW message or body]


 [cc'ing to r-sig-mixed-models]

On 16-08-10 09:44 AM, Chris wrote:
> I'm hoping to do some posterior predictive simulation with a ZIP
> model fit in GLMMADMB  (to better assess model fit and to compare its
> predictive ability with models fit with glmer).  I understand there
> is no simulate method implemented yet for glmmADMB, and the "predict"
> method (and code that I've been able to locate, such as on your
> wikidot GLMM page) does not appear to include random effects.  My
> model has a single random effect (intercept only) that I would like
> to include (so that I have output analogous to that of the "simulate"
> method in lme4).  Given that it is not an incredibly complex model,
> I'm wondering if this is something that I can tackle myself;  I'm
> assuming I would have to adapt the "predict" code to include the
> estimates of the random intercepts and their variability (which
> appear to be "U" and "sd_U" in the fitted glmmADMB model object) and
> perhaps as well  the zero-inflation coefficient (pz" and "sd_pz"?).
> Can you help to point me in the right direction as to how I might
> achieve this?

  I think you're on the right track.  I would look at
glmmADMB:::predict.glmmadmb (or just use it) for constructing the
fixed-effects part of the simulation.  Then for a single random
intercept, you can either simulate from the estimated distribution of
conditional modes *unconditionally* (which is what lme4 does),

  uvals <- rnorm(length(fm$U),mean=0,sd=sqrt(c(VarCorr(fm)[[1]]))

where fm is the fitted model, or simulate *conditionally* on the
observed values

  uvals <- rnorm(length(fm$U),mean=fm$U,sd=fm$sd_U)

Then add uvals to the predictions (pred_vals <- pred_vals +
uvals[grpvar]); transform back to the response scale (e.g. exp-transform
for NB models); add zero-inflation if necessary

val <- ifelse(runif(nobs)<fm$pz,0,rnbinom(nobs,mu=pred_vals,size=fm$alpha))

if you use NB1 or some more exotic distribution it may be harder to
simulate ...

  cheers
    Ben Bolker

> 
> Many thanks, Chris
>

_______________________________________________
R-sig-mixed-models@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
[prev in list] [next in list] [prev in thread] [next in thread]