# How does R function summary.glm calculate the covariance matrix for glm model?

I would like to know how the covariance matrix of estimated coefficients is actually calculated. The code uses QR-decomposition and inversion of some sort. I have an idea that it would go something like this:

$(X'X)^{-1}=[(QR)'QR]^{-1}=(R'R)^{-1}=\Sigma$

Could someone explain the code?

p <- object$rank p1 <- 1L:p Qr <- qr.lm(object) covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
covmat <- dispersion * covmat.unscaled

• Welcome Kati. You use the dollar sign before and after your equation for it to be displayed as math. – Antoine Vernet Jul 18 '16 at 11:57
• before apply $X'X$ you need to subtract the mean for each column, otherwise it will not be covariance matrix – Haitao Du Jul 18 '16 at 18:51

## 2 Answers

?chol2inv says:

Invert a symmetric, positive definite square matrix from its Choleski decomposition. Equivalently, compute (X'X)^(-1) from the (R part) of the QR decomposition of X.

and if you want to dig deeper into the numerics:

This is an interface to the LAPACK routine ‘DPOTRI’

and the documentation for DPOTRI says:

DPOTRI computes the inverse of a real symmetric positive definite matrix A using the Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF. [here **T denotes transposition]

The point is that once you have the $$R$$ component, you don't need to take the crossproduct and then invert - you can much more easily compute the inversion directly.

Hence, the chol2inv(Qr\$qr[p1, p1, drop = FALSE]) computes $$(R^\top R)^{-1}=(X^\top WX)^{-1}$$, where $$R$$ is the upper triangular matrix from the QR decomposition $$QR(\sqrt{W}X)$$.

I was confused about this several months ago. After reading some material, I have a clue about the problem. But I am not sure that the following is correct.

Suppose the distribution of response $$Y$$ is an exponential family, $$X^TX>0$$. $$p(Y|\eta)\sim \exp(\phi^{-1}\eta^T T(Y)-A(\eta)+B(Y))$$, $$\eta=g(E(Y|X))=X\beta$$, $$g$$ is the canonical link, $$\phi$$ is the dispersion parameter of the distribution.

Fitting GLM uses MLE, and we have the Wald-type CI for MLE. The variance-covariance matrix is $$I(\hat\beta)^{-1}$$. $$I(\hat\beta)=-E(\dfrac{\partial^2}{\partial\beta^2}\log p(Y|\eta)|\hat\beta)=\dfrac{\partial^2}{\partial\beta^2}A(X\beta)|\hat\beta=X^T(\dfrac{\partial^2}{\partial\eta^2}A(\eta)|\hat\eta )X=X^TWX$$. Here $$W$$ is a diagonal matrix with weights. The variance-covariance matrix is $$(X^TWX)^{-1}$$.

Edit: An easier way is to use $$I(\beta)=(\dfrac{d \eta}{d \beta})^TI(\eta)\dfrac{d \eta}{d \beta}$$.

reference:

exponential family

Fisher information

• Many thanks for this - that's really helpful - I think this is right! Checked the other answer as the correct one as it addressed my question more closely, but +1! – Tom Wenseleers May 13 '19 at 15:46