# GLM standard errors

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?

• Can you provide a reproducible example? Fit a model with glm() and show how you get the FI matrix? Commented Mar 6, 2018 at 20:09
• 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. Commented Mar 6, 2018 at 20:26

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.

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"