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

Question

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

score 4 · Answer 1 · edited Jul 08 '19 at 08:57

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

score 3 · Answer 2 · edited May 13 '19 at 13:22

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

Many thanks for this - get it now! Oh yes and would you also know if what I did here to calculate 95% confidence intervals on nnls estimates by using bbmle's mle2 is correct? https://stats.stackexchange.com/questions/373253/calculating-the-p-values-in-a-constrained-least-squares/405606#405606 — Tom Wenseleers, May 13 '19 at 15:44

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

2 Answers2