Verificando a normalidade dos resíduos em modelos lineares generalizados

Este artigo utiliza modelos lineares generalizados (distribuições binomiais e negativas de erros binomiais negativos) para analisar dados. Mas então, na seção de análise estatística dos métodos, há esta declaração:

... e segundo modelando os dados de presença usando Modelos de Regressão Logística e os dados de tempo de forrageamento usando um Modelo Linear Generalizado (GLM). Uma distribuição binomial negativa com uma função de link de log foi usada para modelar os dados do tempo de forrageamento (Welsh et al. 1996) e a adequação do modelo foi verificada pelo exame de residentes (McCullagh & Nelder 1989). Os testes de Shapiro – Wilk ou Kolmogorov – Smirnov foram utilizados para testar a normalidade, dependendo do tamanho da amostra; os dados foram transformados em log antes das análises para aderir à normalidade.

Se eles assumem distribuições binomiais e negativas de erros binomiais negativos, certamente não devem verificar a normalidade dos resíduos?

Nota: os resíduos de desvio (ou Pearson) não devem ter uma distribuição normal, exceto para um modelo gaussiano. Para o caso de regressão logística, como diz @Stat, resíduos deviance para o th observação y i são dadas por$i$$y_i$

se e$y_i=0$

se , onde ^ π i é a probabilidade ajustada de Bernoulli. Como cada um pode ter apenas um dos dois valores, fica claro que sua distribuição não pode ser normal, mesmo para um modelo especificado corretamente:$y_i=1$$\hat{\pi_i}$

#generate Bernoulli probabilities from true model
x <-rnorm(100)
p<-exp(x)/(1+exp(x))

#one replication per predictor value
n <- rep(1,100)
#simulate response
y <- rbinom(100,n,p)
#fit model
glm(cbind(y,n-y)~x,family="binomial") -> mod
#make quantile-quantile plot of residuals
qqnorm(residuals(mod, type="deviance"))
abline(a=0,b=1)

$n_i$$i$

$y_i$$n_i$$n_i$

#many replications per predictor value
n <- rep(30,100)
#simulate response
y<-rbinom(100,n,p)
#fit model
glm(cbind(y,n-y)~x,family="binomial")->mod
#make quantile-quantile plot of residuals
qqnorm(residuals(mod, type="deviance"))
abline(a=0,b=1)

As coisas são semelhantes para Poisson ou GLMs binomiais negativos: para contagens previstas baixas, a distribuição de resíduos é discreta e distorcida, mas tende à normalidade para contagens maiores em um modelo especificado corretamente.

Não é comum, pelo menos não no meu pescoço, realizar um teste formal de normalidade residual; se o teste de normalidade é essencialmente inútil quando seu modelo assume a normalidade exata, a fortiori é inútil quando não. No entanto, para modelos não saturados, o diagnóstico gráfico de resíduos é útil para avaliar a presença e a natureza da falta de ajuste, normalizando com uma pitada ou um punhado de sal, dependendo do número de repetições por padrão preditivo.

O que eles fizeram está correto! Vou lhe dar uma referência para verificar novamente. Consulte a Seção 13.4.4 em Introdução à análise de regressão linear, 5ª ediçãopor Douglas C. Montgomery, Elizabeth A. Peck, G. Geoffrey Vining. Em particular, veja os exemplos na página 460, onde eles se encaixam em um binômio glm e verifique duas vezes a suposição de normalidade dos "Deviance Residuals". Como mencionado na página 458, isso ocorre porque "os resíduos de desvio se comportam como os resíduos comuns em um modelo de regressão linear da teoria normal padrão". Portanto, faz sentido se você os plotar na escala normal de plotagem de probabilidade, bem como vs valores ajustados. Novamente, consulte a página 456 da referência acima. Nos exemplos que eles forneceram nas páginas 460 e 461, não apenas para o caso binomial, mas também para o Poisson glm e o Gamma com (link = log), eles verificaram a normalidade dos resíduos de desvio.
Para o caso binomial, o desvio residual é definido como:

> attach(npk)

> #Fitting binomila glm
> fit.1=glm(P~yield,family=binomial(logit))
> 
> #Getting deviance residuals directly
> rd=residuals(fit.1,type = c("deviance"))
> rd
         1          2          3          4          5          6          7 
 1.1038306  1.2892945 -1.2912991 -1.1479881 -1.1097832  1.2282009 -1.1686771 
         8          9         10         11         12         13         14 
 1.1931365  1.2892945  1.1903473 -0.9821829 -1.1756061 -1.0801690  1.0943912 
        15         16         17         18         19         20         21 
-1.3099491  1.0333213  1.1378369 -1.2245380 -1.2485566  1.0943912 -1.1452410 
        22         23         24 
 1.2352561  1.1543163 -1.1617642 
> 
> 
> #Estimated success probabilities
> pi.hat=fitted(fit.1)
> 
> #Obtaining deviance residuals directly
> rd.check=-sqrt(2*abs(log(1-pi.hat)))
> rd.check[P==1]=sqrt(2*abs(log(pi.hat[P==1])))
> rd.check
         1          2          3          4          5          6          7 
 1.1038306  1.2892945 -1.2912991 -1.1479881 -1.1097832  1.2282009 -1.1686771 
         8          9         10         11         12         13         14 
 1.1931365  1.2892945  1.1903473 -0.9821829 -1.1756061 -1.0801690  1.0943912 
        15         16         17         18         19         20         21 
-1.3099491  1.0333213  1.1378369 -1.2245380 -1.2485566  1.0943912 -1.1452410 
        22         23         24 
 1.2352561  1.1543163 -1.1617642 
> 

Verifique aqui também o caso Poisson.