Erros padrão do GLM

7

Eu tenho uma pergunta sobre como obter os erros padrão dos coeficientes no meu modelo GLM. Eu tenho a matriz de informações de fisher que calculei manualmente, mas não é dimensionada. Como posso escalar a matriz de informações do fisher para obter os mesmos erros padrão da função GLM?

r generalized-linear-model covariance T. Nestor
fonte
Você pode fornecer um exemplo reproduzível? Ajuste um modelo com glm () e mostre como você obtém a matriz FI?
Atiretoo - restabelecer monica
Calculei a matriz de informações de Fisher usando , em que é a matriz de peso. Quando olho para os meus resultados, eles correspondem exatamente ao que recebo quando olho para o resumo (m1) . No entanto, isso não me dá os erros padrão de cada coeficiente. Preciso calcular algum parâmetro de escala ? O parâmetro de dispersão para o meu modelo é 1.(XTWX)−1Wcov.unscaledϕ
T. Nestor

Respostas:

6

Você está muito perto! Os erros padrão dos coeficientes são as raízes quadradas da diagonal da sua matriz, que é o inverso da matriz de informações de Fisher. Aqui está um exemplo.


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"
atiretoo - restabelecer monica
fonte
5

Como posso escalar a matriz de informações do fisher para obter os mesmos erros padrão da função GLM?

Cronometre sua matriz de co-variância não dimensionada no parâmetro de dispersão, conforme feito em summary.glm. O código relevante de summary.glmé

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 ...]

Os chol2inv(Qr$qr[p1, p1, drop = FALSE])cálculos quais você faz um comentário. Aqui, é a matriz triangular superior do Code decomposição .(R⊤R)−1=(X⊤WX)−1RQR=WX

as respostas do atiretoo só são válidas quando o parâmetro de dispersão é um, como ocorre com a distribuição de Poisson e Binomial.

Benjamin Christoffersen
fonte