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

covariancerregression

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

Best Answer

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