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

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}$$.

• I'm not entirely sure mathematically what matrix you're asking for. The covariance matrix of the observation-level fitted values incorporating the random effects? Something else? Can you please clarify? Commented Jun 26, 2019 at 16:04
• Hello, I need an estimation of the covariance matrix $n^{-1}\Sigma^{-1}$ where for a regression parameter estimate $\hat{\alpha}$, $\sqrt{n}(\hat{\alpha} - \alpha) \rightarrow N(0, \Sigma^{-1})$ and $\Sigma$ is the limiting value of the partial likelihood information matrix normalized through division by $n$. So, I need the observed inverse information matrix evaluated at $\hat{\alpha}$. Commented Jun 27, 2019 at 7:54

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.)

• Thanks a lot for your answer. It is exactly what I needed. I found the merDerivpackage with the function vcov(..., option = FULL) which gives the same result than your procedure. I am afraid I did not understand your warning : my model has two independant random effects. I need their variance estimation. What do you mean by " if you need this in variance-covariance parameterization" ? Commented Jun 27, 2019 at 12:18
• If the random effects are both scalar (i.e. something like (1|A) + (1|B) then the parameters for which you're computing the variance-covariance matrix are their standard deviations. If you have random-slopes models it gets more complicated. Commented Jun 27, 2019 at 22:34
• My random parameters are written : (0 + A | B ) + (1 | B ). So, I suppose I have a random-slopes model. It is written as $\lambda_i = \lambda_0 \exp(b_j + a_j X_i + \beta Z)$. I need the observed inverse information matrice for the random and fixed effects. The random parameters are independant. What should I do in this case ? Commented Jun 28, 2019 at 9:07

In the GLMMadaptive package the vcov() method returns the covariance matrix of the maximum likelihood estimates for both the fixed effects coefficients and the parameters of the variance-covariance matrix of the random effects (the later in the log-Cholesky factor scale).

For an example, check here.

• Thanks for your answer. I was using the GLMMadaptive package till recently. But I've had estimation problems with it ("Model convergence problem; non-positive-definite Hessian matrix") which was not the case while applying glmer from lme4 on the same set of data! I also get very high values for the covariance of the random effects with the vcov() function on glmmTMB which let me think there is some problem with it. Commented Jun 27, 2019 at 7:48
• Because I do not have access to your data, it is difficult to see why you get the warning. But in any case, you could try changing some of the defaults. For example, increasing the number of EM iterations via iter_EM or changing the starting values. Commented Jun 27, 2019 at 12:41
• I'm simulating survival data and I do not understand why GLMMadaptative does not work (sometimes) where other packages get there. About 15% of my simulations fail because of a definite positive Hessian matrix. When I use another way to simulate data, the problem disappears. But my main problem is that I do not trust the vcov (.., full = TRUE) function on glMM because I get very high values for the random covariance estimate. Commented Jun 27, 2019 at 13:18