7
$\begingroup$

I have a question about how to get the standard errors of the coefficients in my GLM model. I have the fisher information matrix which I calculated by hand, but it is unscaled. How can I scale the fisher information matrix so that I get the same standard errors from the GLM function?

$\endgroup$
  • $\begingroup$ Can you provide a reproducible example? Fit a model with glm() and show how you get the FI matrix? $\endgroup$ – atiretoo Mar 6 '18 at 20:09
  • $\begingroup$ I calculated the fisher information matrix using $(X^{T}WX)^{-1}$, where $W$ is the weight matrix. When I look at my results, they match exactly with what I get when I look at 'summary(m1)$cov.unscaled$. However this doesn't give me the standard errors of each coefficient. Do I have to calculate some scale parameter $\phi$? The dispersion parameter for my model is 1. $\endgroup$ – T. Nestor Mar 6 '18 at 20:26
6
$\begingroup$

You're very close! The standard errors of the coefficients are the square roots of the diagonal of your matrix, which is the inverse of the Fisher information matrix. Here is an example.


data <- caret::twoClassSim()
model <- glm(Class~TwoFactor1*TwoFactor2, data = data, family="binomial")
# here are the standard errors we want
SE <- broom::tidy(model)$std.error

X <- model.matrix(model)
p <- fitted(model)
W <- diag(p*(1-p))
# this is the covariance matrix (inverse of Fisher information)
V <- solve(t(X)%*%W%*%X)
all.equal(vcov(model), V)
#> [1] "Mean relative difference: 1.066523e-05"
# close enough

# these are the standard errors: take square root of diagonal 
all.equal(SE, sqrt(diag(V)))
#> [1] "names for current but not for target"  
#> [2] "Mean relative difference: 4.359204e-06"
$\endgroup$
5
$\begingroup$

How can I scale the fisher information matrix so that I get the same standard errors from the GLM function?

Time your unscaled co-variance matrix the dispersion paramter as done in summary.glm. The relevant code from summary.glm is

if (is.null(dispersion)) 
    dispersion <- if (object$family$family %in% c("poisson", 
        "binomial")) 
        1
    else if (df.r > 0) {
        est.disp <- TRUE
        if (any(object$weights == 0)) 
            warning("observations with zero weight not used for calculating dispersion")
        sum((object$weights * object$residuals^2)[object$weights > 
            0])/df.r
    }
    else {
        est.disp <- TRUE
        NaN
    }
# [other code...]
if (p > 0) {
    p1 <- 1L:p
    Qr <- qr.lm(object)
    coef.p <- object$coefficients[Qr$pivot[p1]]
    covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
    dimnames(covmat.unscaled) <- list(names(coef.p), names(coef.p))
    covmat <- dispersion * covmat.unscaled
    # [more code ...]

The chol2inv(Qr$qr[p1, p1, drop = FALSE]) computes $(R^\top R)^{-1}=(X^\top WX)^{-1}$ which you make a comment about. Here, $R$ is the upper triangular matrix from the QR decomposition $QR=\sqrt{W}X$.


atiretoo answers only holds when the dispersion paramter is one as with the Poisson and Binomial distribution.

$\endgroup$

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.