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 the`lme4`

library.

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

**Contents**hide

#### Best Answer

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 <- [email protected]$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.)

### Similar Posts:

- Solved – How are the p-values in the glmer (lme4) output calculated
- Solved – How to calculate the robust standard error of predicted y from a linear regression model in R?
- Solved – Using Singular Value Decomposition to Compute Variance Covariance Matrix from linear regression model
- Solved – Mixed effects negative binomial with robust standard errors (Huber-white) in R
- Solved – Variance-covariance matrix of survival model