[R-sig-ME] About glmer.nb

[Ben Bolker]

For GLMMs, R computes the full Hessian of all of the top-level
(theta=var-cov parameters + beta=fixed-effect parameters); it is
available, with a little digging, in the fitted object.

Set up example:

library(lme4)
## from ?glmer.nb
set.seed(101)
dd <- expand.grid(f1 = factor(1:3),
                  f2 = LETTERS[1:2], g=1:9, rep=1:15,
                  KEEP.OUT.ATTRS=FALSE)
mu <- 5*(-4 + with(dd, as.integer(f1) + 4*as.numeric(f2)))
dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5)
m.nb <- glmer.nb(y ~ f1*f2 + (1|g), data=dd)


   Now pull out the full variance-covariance matrix:

## extract Hessian
hh <- m.nb at optinfo$derivs$Hessian
## invert
vv <- solve(hh)
## double/symmetrize (move from deviance to log-likelihood scale)
v2 <- unname(as.matrix(forceSymmetric(vv + t(vv))))
## compare ...
v3 <- unname(as.matrix(vcov(m.nb)))
all.equal(v2[-1,-1],v3)

Note that this is the covariance between the fixed-effect parameters
and the "theta" parameters, i.e. the lower triangle of the Cholesky
factor ... in the simple case of a scalar random effect, this reduces
to the standard deviation of the random effect. (In the LMM case you
also have to account for the fact that the thetas give the _scaled_
random effects var-cov, i.e. relative to the residual variance.)

  If you're doing lots of NB fitting you might also be interested in the
glmmTMB package (on github ...) -- there, vcov(.,full=TRUE) gives the
desired full variance-covariance matrix.

On 16-11-21 10:43 AM, Levine, Michael wrote:
> 
> While working on one of my projects, I encountered the function
> glmer.nb that is a part of the package lme4 that you maintain.  As
> far as I can see, one can only obtain a variance-covariance matrix of
> estimated fixed effects using this function...Is there any way for me
> to obtain a full variance-covariance matrix that includes covariances
> between fixed effects and estimated variance of the Gaussian random
> effect?
> 
> Thank you very much in advance for your help!
>