6
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ Welcome Kati. You use the dollar sign before and after your equation for it to be displayed as math. $\endgroup$ – Antoine Vernet Jul 18 '16 at 11:57
  • $\begingroup$ before apply $X'X$ you need to subtract the mean for each column, otherwise it will not be covariance matrix $\endgroup$ – Haitao Du Jul 18 '16 at 18:51
3
+50
$\begingroup$

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

$\endgroup$
1
3
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ 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! $\endgroup$ – Tom Wenseleers May 13 '19 at 15:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.