# Solved – lme4::glmer : Get the covariance matrix of the fixed and random effect estimates

covariance-matrixlme4-nlmemixed modelrandom-effects-model

My problem may seem easy but I have found no satisfactory solution. I am stuck on this problem for a few days already.
How to obtain the covariance matrix of the fixed AND random effects estimates while using the glmer function in thelme4library.

I tried vcov(.., full = TRUE) without success.

Is there a function or a way to calculate this matrix of variance-covariance ?

Edit

The covariance matrix I need is $$n^{−1}\Sigma^{−1}$$. For a regression parameter estimate $$\hat{\alpha}$$, $$\sqrt{n}(\hat{\alpha}−\alpha)\rightarrow N(0,\Sigma^{−1})$$. $$\Sigma$$ is the limiting value of the partial likelihood information matrix normalized through division by $$n$$. In brief, I need the observed inverse information matrix evaluated at $$\hat{\alpha}$$.

Unless you've gone out of your way to not compute the Hessian, it's hiding in the output model structure. You can look in lme4:::vcov.merMod to see where these computations come from (what's there is more complicated because it handles a bunch of edge cases; it also extracts just the part of the covariance matrix relevant to the fixed effects ...)

Example:

library(lme4)
object <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp,
family = binomial)


This extracts the Hessian, inverts it, and doubles it (since the Hessian is computed on the (-2 log likelihood) scale. The h+t(h) is a clever way to improve symmetry while doubling (if I recall correctly ...)

h <- object@optinfo$$derivs$$Hessian
h <- solve(h)
v <- forceSymmetric(h + t(h))


Check that the fixed-effect part agrees (random-effect parameters come first):

all.equal(unname(as.matrix(vcov(object))),
unname(as.matrix(v)[-1,-1])) ## TRUE


Warning: the random effects are parameterized on the Cholesky scale (i.e., the parameters are the lower triangle, in column-major order, of the Cholesky factor of the random effect covariance matrix) ... if you need this in variance-covariance parameterization, or in standard deviation-correlation parameterization, it's going to take more work. (If you only have a single scalar random effect, then the parameter is the standard deviation.)