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

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