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:
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