Skip to content
Prev 11271 / 20628 Next
Original Raw

Fit negative binomial glmm credible intervals from the predict function?

Ben Bolker 
Sheryn Olson <sheryn.olson at ...> writes:
In case it helps, what you're assuming is Normal is the *sampling
distribution of the parameters*, not the data themselves ...

  The bottom line is that what you're doing seems reasonable, with
the caveats already expressed in ?predict.glmmadmb ...
Well, technically it's a discrete (count) variable, not continuous ...

  Section 2 of the glmmADMB vignette [vignette("glmmADMB")] has
a bit of information about running post-fit MCMC.   There should
be an example of how to generate confidence intervals (basically,
you would generate predictions for each sampled set of coefficients
and then compute the quantiles or credible intervals of the
predictions), but there isn't (yet ...)
Good question.  Not easy, you may need to take some shortcuts
and hope for the best.
It depends what values you want to predict for.  The data frame
typically includes "typical" values of the parameters; you're welcome
to compute predictions for the medians instead if you prefer, or
for a range of values.
??  Do you mean the 'mcmc' argument of glmmadmb?
[Update: I just wrote a bunch of stuff about how se.fit doesn't
work for glmmadmb, when it indeed does! However, I had already
written all the stuff below, which describes more generally how
you would do it by hand, so I'm just going to send it anyway]

  The basic recipe for constructing confidence intervals, as laid out
in http://glmm.wikidot.com/faq in the "Predictions and/or confidence
(or prediction) intervals on predictions" section, is

 (1) construct the model matrix (X)
  * If you want predictions on the original model (not this case),
the model.matrix() accessor may work to extract it from the fit for you
  * If not, or if you want to construct predictions for new data,
you need to call model.matrix() with the fixed-effect formula (you
may or may not be able to extract the fixed-effect part of the
formula only from the model, but it's usually fairly easy just
to respecify it

  (2) extract the variance-covariance matrix of the fixed effect
predictors, V (usually just vcov(model))

  (3) extract the fixed-effect coefficients (usually just fixef(model)
or coef(model))

 then the predictions (on the linear predictor scale, at the population
level [i.e. ignoring all random effects]) are

 X %*% beta

(plus an offset if there is one)

and the standard errors of prediction, *CONDITIONAL ON THE
RANDOM EFFECTS* [i.e. ignoring all uncertainty due to uncertainty
of the random effects] are

  sqrt(diag(X %*% V %*% t(X)))

 Then to get confidence intervals you should compute

  inverse_link(fit +/- 1.96*se.fit) exactly as you have below

This assumes further that the sampling distribution of the
coefficients is Normal (in which case the linear computation
we did above will also lead to a Normal distribution, on 
the linear predictor scale).  This might not be true for
small or badly behaved data sets (of course it is never
exactly true except asymptotically ...), but there's not
much you can do about it without working much harder.
Do you mean you're suspicious of the code, or that the numbers
it produces seem wrong?

Thread (2 messages)

Sheryn Olson Fit negative binomial glmm credible intervals from the predict function? Jan 1 Ben Bolker Fit negative binomial glmm credible intervals from the predict function? Jan 6