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

