Figura 1: Funções de ligação logit, probit, log-log e complementar log-log para dados binomiais.
Modelos Lineares Generalizados
Devido originalmente a Nelder and Wedderburn (1972), os modelos lineares generalizados são uma síntese e extensão notáveis de modelos de regressão familiares, como os modelos lineares. O presente texto começa com uma consideração da estrutura geral e do alcance de aplicação dos modelos lineares generalizados; passa a examinar em mais detalhes modelos lineares generalizados para dados de contagem, incluindo tabelas de contingência; esboça brevemente a teoria estatística subjacente aos modelos lineares generalizados e conclui com a extensão dos diagnósticos de regressão para modelos lineares generalizados.
Modelos lineares generalizados tornaram-se tão centrais para a análise de dados estatísticos eficazes, entretanto, que vale a pena o esforço adicional necessário para adquirir um conhecimento básico do assunto.
Como suporte computacional utilizamos a linguagem de programação e ambiente de desenvolvimento integrado para cálculos estatísticos e gráficos R, versão 4.0.2 (2020-06-22) "Taking Off Again", especialmente a função glm e o pacote gamlss. Estas notas estão baseadas no livro de John Fox (2016).
Um modelo linear generalizado (ou GLM) consiste em três componentes:
| Ligação | \(\eta_i=g(\mu_i)\) | \(\mu_i=g^{-1}(\eta_i)\) | ||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Identidade | \(\mu_i\) | \(\eta_i\) | ||||||||||||||
| Log | \(\log_e\big(\mu_i\big)\) | \(e^{\eta_i}\) | ||||||||||||||
| Inversa | \(\mu_i^{-1}\) | \(\eta_i^{-1}\) | ||||||||||||||
| Inversa quadrada | \(\mu_i^{-2}\) | \(\eta_i^{-1/2}\) | ||||||||||||||
| Raiz quadrada | \(\sqrt{\mu_i}\) | \(\eta_i^2\) | ||||||||||||||
| Logit | \(\log_e\Big(\dfrac{\mu_i}{1-\mu_i}\Big)\) | \(\dfrac{1}{1+\exp\big(-\eta_i\big)}\) | ||||||||||||||
| Probit | \(\Phi^{-1}\big(\mu_i\big)\) | \(\Phi\big(\eta_i\big)\) | ||||||||||||||
| Log-log | \(-\log_e\big(-\log_e\big(\mu_i\big)\big)\) | \(\exp\big(-\exp\big(-\eta_i\big)\big)\) | ||||||||||||||
| Complementar log-log | \(\log_e\big(-\log_e\big(1-\mu_i\big)\big)\) | \(1-\exp\big(-\exp\big(\eta_i\big)\big)\) |
NOTA: \(\mu_i\) é o valor esperado da resposta; \(\eta_i\) é o preditor linear e \(\Phi(\cdot)\) é a função de distribuição normal padrão.
As funções de ligação comumente empregadas e seus inversos são mostrados na Tabela 1. Observe que a ligação identidade simplesmente retorna seu argumento inalterado, \begin{equation} \eta_i \, = \, g(\mu_i) \, = \, \mu_i \qquad \mbox{ e, portanto, } \qquad \mu_i \, = \, g^{−1}(\eta_i) \, = \, \eta_i\cdot \end{equation}
As últimas quatro funções de ligação na Tabela 1 são para dados binomiais, onde \(Y_i\) representa a proporção observada de sucessos em \(n_i\) tentativas binárias independentes; assim, \(Y_i\) pode assumir qualquer um dos valores \(0,1/n_i,2/n_i,\cdots,(n_i-1)/n_i,1\). Recordemos que os dados binomiais também abrangem dados bin´rios, onde todas as observações representam \(n_i = 1\) tentativas e, consequentemente, \(Y_i\) é 0 ou 1. A esperança da resposta \(\mu_i = \mbox{E}(Y_i)\) é então a probabilidade de sucesso. As ligações logit, probit, log-log e complementar log-log estão representados graficamente na Figura 1. Em contraste com as ligações logit e probit que, como observamos anteriormente, são quase indistinguíveis uma vez que as variâncias das distribuições normal e logística subjacentes são equacionadas, as ligações log-log e complementar log-log se aproximam das assíntotas de 0 e 1 assimetricamente.
Além do desejo geral de selecionar uma função de ligação que torne a regressão de \(Y\) nos \(X\)s linear, uma ligação promissora removerá as restrições no intervalo da resposta esperada. Esta é uma ideia familiar dos modelos logit e probit, onde o objetivo é modelar a probabilidade de sucesso, representada por \(\mu_i\) em nosso atual geral notação. Como probabilidade, \(\mu_i\) está confinada ao intervalo unitário [0,1]. As ligações logit e probit mapeiam este intervalo para toda a linha real, de \(-\infty,\infty\). Da mesma forma, se a resposta \(Y\) for uma contagem, assumindo apenas valores inteiros não negativos, \(0, 1, 2,\cdots\) e, conseqüentemente, \(\mu_i\) é uma contagem esperada que, embora não necessariamente um número inteiro também não é negativa, a ligação logaritmo mapeia \(μ_i\) para toda a linha real. Isso não quer dizer que a escolha da função de ligação seja inteiramente determinada pelo intervalo da variável resposta.
Figura 1: Funções de ligação logit, probit, log-log e complementar log-log para dados binomiais.
Uma propriedade conveniente das distribuições na famílias exponencial é que a variância condicional de \(Y_i\) é uma função de sua média \(\mu_i\), digamos, \(V(\mu_i)\) e, possivelmente, um parâmetro de dispersão \(\phi\). As funções de variância para as famílias exponenciais comumente usadas aparecem na Tabela 2. A variância condicional da resposta na família Gaussiana é uma constante \(\phi\), que é simplesmente uma notação alternativa para o que anteriormente denominamos variância do erro \(\sigma_\epsilon^2\). Nas famílias binomial e Poisson, o parâmetro de dispersão é definido com o valor fixo \(\phi=1\).
A Tabela 2 também mostra a faixa de variação da variável resposta em cada família e a função de elo chamada canônica ou natural associada a cada componente na família. A ligação canônica simplifica o GLM, mas outras funções de ligação também podem ser usadas. Na verdade, um dos pontos fortes do GLM - em contraste com as transformações da variável resposta na regressão linear - é que a escolha da transformação linearizante é parcialmente separada da distribuição da resposta, e a mesma transformação não precisa normalizar a distribuição de \(Y\) e fazer sua regressão linear nos \(X\).
| Família | Ligação canônica | Intervalo da resposta | \(\mbox{Var}(Y_i \, | \, \eta_i)\) | |||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Gaussiana | Identidade | \(-\infty.+\infty)\) | \(\phi\) | |||||||||||||||||||
| Binomial | Logit | \(0,1/n_i,\cdots,n_i/n_i\) | \(\mu_i(1-\mu_i)/n_i\) | |||||||||||||||||||
| Poisson | Log | \(0,1,2,\cdots\) | \(\mu_i\) | |||||||||||||||||||
| Gama | Inversa | \((0,+\infty)\) | \(\phi\mu_i^2\) | |||||||||||||||||||
| Normal inversa | Inversa quadrada | \((0,+\infty)\) | \(\phi\mu_i^3\) |
NOTA: \(\phi\) é o parâmetro de dispersão, \(\eta_i\) é o preditor linear e \(\mu_i\) é a esperança de \(Y_i\) (a resposta). Na família binomial, \(n_i\) é o número de tentativas independentes.
Há também esta diferença mais sutil: quando transformamos \(Y\) e regredimos a resposta transformada nos \(X\), nós estamos modelando a esperança da resposta transformada, \begin{equation} \mbox{E}\big( g(Y)\big) \, = \, \alpha+\beta_1 X_{i1}+\beta_2 X_{i2}+\cdots +\beta_k X_{ik}\cdot \end{equation} Em um modelo linear generalizado, em contraste, modelamos a esperança transformada da resposta, \begin{equation} g\big( \mbox{E}(Y)\big) \, = \, \alpha+\beta_1 X_{i1}+\beta_2 X_{i2}+\cdots +\beta_k X_{ik}\cdot \end{equation} Embora semelhante em espírito, isso não é exatamente a mesma coisa quando a função de ligação \(g(\cdot)\) é não linear.
As funções de ligação específicas que podem ser usadas variam de uma família para outra. Por exemplo, não seria promissor usar as ligações identidade, log, inversa, inversa quadrado ou raiz quadrada com dados binomiais, nem seria sensato usar o logit, probit, log-log ou complementar log-log com dados não binomiais.
Presumimos que o leitor esteja geralmente familiarizado com as famílias gaussiana e binomial e simplesmente apresentamos suas distribuições aqui para referência. As distribuições Poisson, gama e Gaussiana inversa são talvez menos familiares, então fornecemos mais alguns detalhes:
Figura 2: Distribuições Poisson para vários valores do parâmetro de taxa \(\mu\).
Figura 3: Várias distribuições gama para a escala \(\omega = 1\) e vários valores do parâmetro de forma \(\psi\).
Figura 4: Distribuições gaussianas inversas para várias combinações de valores da média \(\mu\) e dispersão \(1/\lambda\).
Os modelos lineares generalizados são ajustados aos dados pelo método de máxima verossimilhança, fornecendo não apenas estimativas dos coeficientes de regressão, mas também erros padrão assintóticos estimados, ou seja, em amostras grandes dos coeficientes. Para testar a hipótese nula \(H_0 \, : \, \beta_j\,= \, \beta_j^{(0)}\) podemos calcular a estatística de Wald \begin{equation} Z_0 \, = \, \dfrac{\widehat{\beta}_j - \beta_j^{(0)}}{SE(\widehat{\beta}_j)}, \end{equation} onde \(SE(\widehat{\beta}_j)\) é o erro padrão assintótico do coeficiente estimado \(\widehat{\beta}_j\). Sob a hipótese nula, \(Z_0\) segue uma distribuição normal padrão.
Conforme explicado, algumas das famílias exponenciais nas quais os modelos lineares generalizados são baseados incluem um parâmetro de dispersão desconhecido \(\phi\). Embora este parâmetro possa, em princípio, ser estimado por máxima verossimilhança é mais comum usar um estimador pelo método dos momentos, que denotaremos \(\widetilde{\phi}\).
A ANOVA para modelos lineares tem um análogo próximo na análise de desvio para modelos lineares generalizados. No contexto mais geral atual, o desvio residual para um modelo linear generalizado é \begin{equation} D_m \, = \, 2\big(\log_e(L_s)-\log_2(L_m) \big), \end{equation} onde \(L_m\) é a verossimilhança maximizada sob o modelo em questão e \(L_s\) é a verossimilhança maximizada sob um modelo saturado, que dedica um parâmetro a cada observação e conseqüentemente ajusta os dados o mais próximo possível. O desvio residual é análogo e, de fato, é uma generalização da soma residual dos quadrados para um modelo linear.
Em modelos lineares generallizados para os quais o parâetro de dispersão é fixado em 1, ou seja, binomial e Poisson, a estatística de teste da razão de verossimilhanças é simplesmente a diferença nos desvios residuais para modelos aninhados. Suponha que o Modelo 0, com \(k_0+1\) coeficientes, esteja aninhado no Modelo 1, com \(k_1+1\) coeficientes onde, então, \(k_0 < k_1\); mais comumente, o Modelo 0 simplesmente omitiria alguns dos regressores no modelo 1.
Testamos a hipótese nula de que as restrições no Modelo 1 representado pelo Modelo 0 estão corretas calculando a estatística de teste de razão de verossimilhança \begin{equation} G_0^2 \, = \, D_0 \, - \, D_1\cdot \end{equation} Sob a hipótese, \(G_0^2\) é assintoticamente distribuído como qui-quadrado com \(k_1 - k_0\) graus de liberdade.
Os testes de razão de verossimilhanças podem ser invertidos para fornecer intervalos de confiança para coeficientes; os testes e intervalos baseados na estatística da razão de verossimilhanças tendem a ser mais confiáveis do que aqueles baseados na estatística Wald. Por exemplo, o intervalo de confiança de 95% para \(\beta_j\) inclui todos os valores \(β_j'\) para os quais a hipótese \(H_0 \, : \, β_j = β_j'\) é aceitável no nível 0.05, ou seja, todos os valores de \(β_j\) para os quais \(2\big(\log_e(L_1)-\log_2(L_0) \big) \leq \chi^2_{0.05,1} = 3.84\), onde \(\log_e(L_1)\) é o logaritmo da verossimilhança maximizado para o modelo completo e \(log_e(L_0)\) é o logaritmo da verossimilhança maximizado para um modelo no qual \(\beta_j\) é restrito ao valor \(\beta_j'\). Este procedimento é computacionalmente intensivo porque exige a verossimilhança perfilada - reajustar o modelo para vários valores fixos \(\beta_j'\) de \(\beta_j\).
Para um modelo linear generalizado em que há um parâmetro de dispersão para estimar, gaussiano, gamma e gaussiana inversa, podemos, em vez disso, comparar modelos aninhados por um teste \(F\), \begin{equation} F_0 \, = \, \dfrac{\frac{D_0-D_1}{k_1-k_0}}{\widetilde{\phi}}, \end{equation} onde a dispersão estimada \(\widetilde{\phi}\), análoga à variância do erro estimado para um modelo linear é retirada do maior modelo ajustado aos dados, que não é necessariamente o Modelo 1. Se o maior modelo tem \(k+1\) coeficientes, então, sob a hipótese de que as restrições no Modelo 1 representadas pelo Modelo 0 está correto, \(F_0\) segue uma distribuição \(F\) com \(k_1−k_0\) e \(n-k-1\) graus de liberdade. Aplicado a um modelo linear generalizado gaussiano, este é simplesmente o familiar teste \(F\) incremental. O desvio residual dividido pela dispersão estimada, \(D^* ≡ D/\widetilde{\phi}\), é chamado de desvio escalonado. O uso não é totalmente uniforme aqui, e tanto o desvio residual quanto o desvio escalonado costumam ser simplesmente denominados de desvio.
Podemos basear um análogo nos modelos lineares generalizados da correlação múltipla ao quadrado no desvio residual: Seja \(D_0\) o desvio residual para o modelo incluindo apenas a constante de regressão \(\alpha\) - denominado desvio nulo - e \(D_1\) o desvio residual para o modelo em questão. Então, \begin{equation} R^2 \, = \, 1-\dfrac{D_1}{D_0}, \end{equation} representa a proporção do desvio nulo contabilizado pelo modelo.
O modelo linear generalizado básico para dados de contagem é o modelo de Poisson com ligação de logaritmo (log). Considere, a título de exemplo, os dados de Michael Ornstein sobre diretorias interligadas entre 248 empresas canadenses dominantes, exemplo a seguir. O número de interligamento para cada empresa é o número de laços que uma empresa manteve em virtude de seus membros do conselho e os principais executivos também atuarem como membros do conselho ou executivos de outras empresas no conjunto de dados. Ornstein estava interessado na regressão do número de bloqueios em outras características das empresas - especificamente, em seus ativos (medidos em bilhões de dólares), nação de controle (Canadá, Estados Unidos, Reino Unido ou outro país) e o principal setor de operação da empresa (10 categorias, incluindo bancos, outras instituições financeiras, manufatura pesada, etc.).
A descrição deste exemplo está disponível em:
Contêm as variáveis:
Figura 5: Distribuição do número de intertravamentos entre 248 corporações canadenses dominantes.
O exame da distribuição do número de intertravamentos, na Figura 5, revela que a variável resposta é altamente enviesada positivamente e que há muitas contagens zero. Embora a distribuição condicional de intertravamentos dadas as variáveis explicativas possa diferir de sua distribuição marginal, a extensão em que a distribuição marginal de intertravamentos se afasta da simetria é um mau presságio para a regressão de mínimos quadrados. Além disso, nenhuma transformação espalhará os zeros.
Os resultados da regressão Poisson do número de intertravamentos nos ativos, na nação de controle e no setor estão resumidos na saída da ANOVA acima. Todos os termos do modelo são, portanto, altamente significativos do ponto de vista estatístico.
Como o modelo usa a ligação logaritmo, podemos interpretar os coeficientes exponenciados, ou seja, o \(e^{\beta_j}\) como efeitos multiplicativos no número esperado de intertravamentos. Assim, por exemplo, mantendo a nação de controle e o setor constantes, aumentando os ativos em 1 bilhão de dólares, a unidade da variável de ativos, multiplica o número estimado estimado de intertravamentos por \(e^{0.020851} = 1.02107\) - ou seja, um aumento de pouco mais de 2%. Da mesma forma, o número estimado esperado de intertravamentos é \(e^{0.8259} = 2.283\) vezes mais alto em uma empresa controlada pelo Canadá do que em uma empresa comparável controlada pelos EUA.
Conforme mencionado, o desvio residual para o ajuste do modelo completo aos dados de Ornstein é \(D_1 = 1887.402\); o desvio para um modelo que se ajusta apenas à constante, ou seja, o desvio nulo é \(D_0 = 3737.010\). Consequentemente, \(R^2 = 1 - 1887.402/3737.010 = 0.495\), revelando que o modelo é responsável por quase metade do desvio no número de intertravamentos.
 |
 |
Figura 6: Número de intertravamentos preditos ou estimados segundo o número de intertravamentos observados. Esquerda: de acordo com a nação de controle e direita: de acordo com o setor industrial.
O modelo de regressão de Poisson é um modelo não linear para a resposta esperada e, portanto, geralmente é mais simples interpretar o modelo graficamente usando exibições de efeito do que examinar os coeficientes estimados diretamente. Os princípios de construção de telas de efeito para modelos lineares generalizados são essencialmente o mesmo que para modelos lineares e para modelos logit e probit. Normalmente construímos uma exibição para cada termo de ordem superior no modelo, permitindo que as variáveis explicativas nesse termo variem sobre seus valores, enquanto mantemos outras variáveis explicativas no modelo com valores típicos. Em um GLM, é vantajoso plotar os efeitos na escala do preditor linear estimado \(\eta\), um procedimento que preserva a estrutura linear do modelo. Em um modelo de Poisson com a ligação logaritmo, o preditor linear est&aaute; na escala logarítmica. Podemos, no entanto, tornar a exibição mais fácil de interpretar, rotulando novamente o eixo vertical na escala da resposta esperada \(\mu\), mais informativamente, fornecendo um segundo eixo vertical no lado direito do gráfico. Para um modelo de Poisson, a resposta esperada é uma contagem.
O desvio residual para o ajuste do modelo de regressão de Poisson aos dados da diretoria de intertravamento, \(D = 1887.4\), é muito maior do que os 234 graus de liberdade residuais do modelo. Se o modelo de Poisson se ajusta aos dados razoavelmente, esperaríamos que o desvio residual fosse aproximadamente igual aos graus de liberdade residuais. Ou seja, a razão entre o desvio residual e os graus de liberdade pode ser tomada como uma estimativa do parâmetro de dispersão \(\phi\) que, em um modelo de Poisson, é fixado em 1. Deve-se notar, entretanto, que este estimador baseado em desvio da dispersão pode ter um desempenho fraco. Um estimador pelo método de momentos geralmente preferível é fornecido na Seção III.
O fato de o desvio residual ser tão grande sugere que a variação condicional do número esperado de intertravamentos excede a variação de uma variável com distribuição Poisson, para a qual a variância é igual à média. Essa ocorrência comum na análise de dados de contagem é chamada de superdispersão.
Embora seja muito menos comum, também é possível que os dados de contagem sejam subdispersos, ou seja, que a variação condicional da resposta seja menor do que a média. A solução para dados de contagem subdispsidos é a mesma que para dados superdispersos; por exemplo, podemos ajustar um modelo quase-Poisson com um parâmetro de dispersão, conforme descrito imediatamente abaixo.
Na verdade, a superdispersão é tão comum em modelos de regressão para dados de contagem e suas consequências são potencialmente tão graves, que modelos como os GLMs de quase Poisson e binomial negativo discutidos em esta seção devem ser empregados como uma coisa natural.
Um remédio simples para dados de contagem superdispersos é introduzir um parâmetro de dispersão no modelo Poisson, de forma que a variância condicional da resposta seja agora \(\mbox{Var}(Y_i \, | \, \eta_i) = \phi \mu_i\). Se \(\phi> 1\), portanto, a variância condicional de \(Y\) aumenta mais rapidamente do que sua média. Não há família exponencial correspondente a esta especificação e o GLM resultante não implica uma distribuição de probabilidade específica para a variável de resposta. Em vez disso, o modelo especifica a méia condicional e variância de \(Y_i\) diretamente. Como o modelo não fornece uma distribuição de probabilidade para \(Y_i\), ele não pode ser estimado por máxima verossimilhança. No entanto, o procedimento usual para estimação por máxima verossimilhança de um GLM produz os chamados estimadores de quase verossimilhança dos coeficientes de regressão, que compartilham muitas das propriedades dos estimadores de máxima verossimilhança.
Acontece que as estimativas de quase verossimilhança dos coeficientes de regressão são idênticas às estimativas de máxima verossimilhança para o modelo Poisson. Os erros padrão dos coeficientes estimados diferem, entretanto: Se \(\widetilde{\phi}\) for a dispersão estimada para o modelo, então os erros padrão dos coeficientes para o modelo de quase-Poisson é \(\widetilde{\phi}^{1/2}\) vezes aquela para o modelo Poisson. No caso de sobredispersão, portanto, onde \(\phi> 1\), o efeito de introduzir um parâmetro de dispersão e obter estimativas de quase-verossimilhança é realisticamente inflar os erros padrão do coeficiente. Da mesma forma, os testes \(F\) para termos no modelo refletirão o parâmetro de dispersão estimado, produzindo estatísticas de teste menores e \(p\)-valores maiores.
Conforme explicado na seção seguinte, usamos o estimador obtido pelo método dos momentos para o par&acir;metro de dispersão. No modelo quase-Poisson, o estimador de dispersão assume a forma \begin{equation} \widetilde{\phi} \, = \, \dfrac{1}{n-k-1}\sum_{i=1}^n \dfrac{\big(Y_i-\widehat{\mu}_i \big)^2}{\widehat{\mu}_i}, \end{equation} onde onde \(\widehat{\mu}_i = g^{−1}(\widehat{\eta}_i)\) é a esperança ajustada de \(Y_i\). Aplicado à regressão de direção interligada de Ornstein, por exemplo, obtemos \(\widetilde{\phi} = 7.943873\) e, portanto, os erros padrão dos coeficientes de regressão para o modelo Poisson na Tabela 3 são cada um multiplicados por \(\sqrt{7.943873} = 2.818488\).
Observo de passagem que existe um modelo quase-binomial semelhante para proporções superdispersas, substituindo o parâmetro de dispersão fixo de 1 na distribuição binomial por um parâmetro de dispersão \(\phi\) a ser estimado a partir dos dados. Dados binomiais superdispersos podem surgir, por exemplo, quando diferentes indivíduos que compartilham os mesmos valores das variáveis explicativas diferem em sua probabilidade \(\mu\) de sucesso, uma situação que é denominada heterogeneidade não modelada. Da mesma forma, a superdispersão pode ocorrer quando as observações binomiais não são independentes, conforme exigido pela distribuição binomial - por exemplo, quando cada observação binomial é para indivíduos relacionados, como membros de uma família.
Existem várias rotas para modelos de contagens com base na distribuição binomial negativa ver, por exemplo, Long (1997) e McCullagh and Nelder (1989). Uma abordagem, seguindo McCullagh and Nelder (1989) é adotar um modelo Poisson para a contagem \(Y_i\), mas supor que a contagem esperada \(\mu_i^*\) é em si uma variável aleatória não observável que é distribuída segundo a distribuição gama com média \(\mu_i\) e parâmetro de escala constante \(\omega\), implicando que o parâmetro de forma gama é \(\psi_i = \mu_i/\omega\). Então, a contagem observada \(Y_i\) segue uma distribuição binomial negativa \begin{equation} P(Y_i=y_i) \, = \, \dfrac{\Gamma(y_i+\omega)}{y_i!\Gamma(\omega)}\dfrac{\mu_i^{y_i}\omega^\omega}{(\mu_i+\omega)^{\mu_i+\omega}}, \end{equation} com valor esperado \(\mbox{E}(Y_i)=\mu_i\) e variância \(\mbox{Var}(Y_i)=\mu_i+\mu_i^2/\omega\). A menos que o parâmetro \(\omega\) seja grande, portanto, a variância de \(Y\) aumenta mais rapidamente com a média do que a variância de uma variável Poisson. Tornar o valor esperado de \(Y_i\) uma variável aleatória incorpora variação adicional entre as contagens observadas para observações que compartilham os mesmos valores das variáveis explicativas e, conseqüentemente, têm o mesmo preditor linear \(\eta_i\).
Com o parâmetro da escala gama \(\omega\) fixado em um valor conhecido, a distribuição binomial negativa é um elemento da família exponencial e um GLM baseado nesta distribuição pode ser ajustado por mínimos quadrados ponderados iterados, conforme desenvolvido na próxima seção. Se em vez disso - e normalmente é o caso - o valor de \(\omega\) é desconhecido e deve, portanto, ser estimado a partir dos dados, os métodos padrão para GLMs baseados em famílias exponenciais não se aplicam. Podemos, entretanto, obter estimativas tanto dos coeficientes de regressão quanto de \(\omega\) pelo método da máxima verossimilhança. Aplicado à regressão de direção interligada de Ornstein e usando a ligação logaritmo, o GLM binomial negativo produz resultados muito semelhantes aos do modelo quase Poisson. O parâmetro de escala estimado para o modelo binomial negativo é \(\widehat{\omega} = 1.2953\), com erro padrão \(SE(\widehat{\omega}) = 0.4627\); temos, portanto, fortes evidências de que a variâcia condicional do número de intertravamentos aumenta mais rapidamente do que seu valor esperado.
Modificando o número de iterações do algoritmo obtemos a seguinte resultado melhorado:Um tipo particular de superdispersão é obtido quando há mais zeros nos dados do que é consistente com uma distribuição Poisson ou binomial negativa, uma situação que pode surgir quando apenas certos membros da população estão em risco de uma contagem diferente de zero. Imagine, por exemplo, que estamos interessados em modelar o número de filhos nascidos de uma mulher. Podemos esperar que esse número seja uma função parcial de variáveis explicativas como estado civil, idade, etnia, religião e uso de anticoncepcionais. Também é provável, no entanto, que algumas mulheres ou seus parceiros sejam inférteis e sejam diferentes das mulheres férteis que, embora em risco de ter filhos, acontecem não ter nenhum. Se soubéssemos quais mulheres são inférteis, poderíamos simplesmente excluí-las da análise, mas suponhamos que não seja esse o caso. Para reiterar, existem duas fontes de zeros nos dados que não podem ser perfeitamente distinguidas: as mulheres que não podem ter filhos e as que não querem ter filhos.
Vários modelos estatísticos têm sido propostos para dados de contagem com excesso de zeros, incluindo o modelo de regressão de Poisson inflado por zeros ou ZIP, devido a Lambert (1992). O modelo ZIP consiste em dois componentes:
Seja \(\pi_i\) representam a probabilidade de que a resposta \(Y_i\) para o \(i\)-ésimo indivíduo seja necessariamente 0. Então \begin{equation} \log_e\left( \dfrac{\pi_i}{1-\pi_i}\right) \, = \, \gamma_0+\gamma_1 Z_{i1}+\gamma_2 Z_{i2}+\cdots+\gamma_p Z_{ip} \end{equation} onde os \(Z_{ij}\) são regressores para prever a adesão &\grave; primeira classe latente; e \begin{equation} \log_e \big( \mu_i\big) \, = \, \alpha+\beta_1 X_{i1}+\beta_2 X_{i2}+\cdots+\beta_p X_{ik} \end{equation} e \begin{equation} p(y_i \, | \, x_1,\cdots,x_k) \, = \, \dfrac{\mu_i^{y_i}e^{-\mu_i}}{y_i!}, \qquad \mbox{para} \qquad y_1=0,1,2,\cdots, \end{equation} onde \(\mu_i = \mbox{E}(Y_i)\) é a contagem esperada para um indivíduo na segunda classe latente e o \(X_{ij}\) são regressores para o submodelo de Poisson. Nas aplicações, os dois conjuntos de regressores - os \(X\)s e os \(Z\)s - costumam ser iguais, mas não é necessariamente o caso. De fato, um caso especial particularmente simples surge quando o submodelo logístico é \begin{equation} \log_e\left( \dfrac{\pi_i}{1 - \pi_i}\right) = \gamma_0, \end{equation} uma constante, implicando que a probabilidade de pertencer à primeira classe latente é idêntica para todas as observações.
A probabilidade de observar uma contagem 0 (zero) é \begin{equation} p(0) \, = \, P(Y_i=0) \, = ,\ \pi_1+(1-\pi_i)e^{-\mu_i} \end{equation} e a probabilidade de observar qualquer contagem diferente de zero \(y_i\) é \begin{equation} p(y_i) \, = \, (1-\pi_i)\times \dfrac{\mu_i^{y_i}e^{-\mu_i}}{y_i!}\cdot \end{equation}
A esperança condicional e a variância de \(Y_i\) são \begin{array}{rcl} \mbox{E}(Y_i) & = & (1-\pi_i)\mu_i \\[0.4em] \mbox{Var}(Y_i) & = & (1-\pi_i)\mu_i(1+\pi_i\mu_i), \end{array} com \(\mbox{Var}(Y_i) > \mbox{E}(Y_i)\) para \(\pi_> 0\), ao contrário de uma distribuição de Poisson pura, para a qual \(\mbox{Var}(Y_i) = \mbox{E}(Y_i) = \mu_i\). Embora esta forma do modelo de contagem inflacionada de zero seja a mais comum, Lambert (1992) também sugeriu o uso de outros GLMs binários para associação na classe latente zero, ou seja, probit, log-log e modelos complementar log-log e o uso alternativo da distribuição binomial negativa para o submodelo de contagem.
A estimação do modelo ZIP seria simples se soubéssemos a qual classe latente cada observação pertence, mas, como indicado, isso não é verdade. Em vez disso, devemos maximizar a log-verossimilhança combinada um pouco mais complexa para os dois componentes do modelo ZIP \begin{array}{rcl} \log_e\Big(L(\beta,\gamma;y) \Big) & = & \displaystyle \sum_{y_i=0} \log_e\left( \exp(Z_i^\top \gamma)+\exp\big( -\exp(X_i^\top \beta)\big)\right) + \displaystyle \sum_{y_i>0} \left( y_iX_i^\top \beta-\exp(X_i^\top \beta)\right) \\[0.4em] & & -\displaystyle \sum_{i=1}^n \log_e\big(1+\exp(Z_i^\top \gamma) \big)-\sum_{y_i>0} \log_e(y_i!), \end{array} onde \(Z_i^\top =(1,Z_{i1},\cdots,Z_{ip})\) , \(X_i^\top=(1,X_{i1},\cdots,X_{ik})\), \(\gamma=(\gamma_0,\gamma_1,\cdots,\gamma_p)^\top\) e \(\beta=(\alpha,\beta_1,\cdots,\beta_p)^\top\).
Fonte: Institute for Digital research & Education. UCLA: Statistical Consulting Group.
https://stats.idre.ucla.edu/ (acessado 12 de maio de 2021)
Os biólogos estaduais da vida selvagem querem modelar quantos peixes estão sendo capturados pelos pescadores em um parque estadual. Os visitantes são questionados sobre quanto tempo permaneceram, quantas pessoas estavam no grupo, se havia crianças no grupo e quantos peixes foram pescados. Alguns visitantes não pescam, mas não há dados sobre se uma pessoa pescou ou não. Alguns visitantes que pescaram não apanharam nenhum peixe, por isso existem zeros em excesso nos dados por causa das pessoas que não pescaram.
Como dizemos a regressão Poisson inflacionada com zeros é usada para modelar dados de contagem que têm um excesso de contagens zero. Além disso, a teoria sugere que os zeros em excesso são gerados por um processo separado dos valores de contagem e que os zeros em excesso podem ser modelados independentemente. Assim, o modelo ZIP tem duas partes, um modelo de contagem de poisson e o modelo logit para prever zeros em excesso.
Temos dados de 250 grupos que foram a um parque. Cada grupo foi questionado sobre quantos peixes pescaram count, quantas crianças estavam no grupo child, quantas pessoas estavam no grupo pearsons e se trouxeram ou não um trailer para o parque camper.
Além de prever o número de peixes capturados, há interesse em prever a existência de zeros em excesso, ou seja, a probabilidade de que um grupo tenha capturado zero peixe. Usaremos as variáveis child, pearsons e camper em nosso modelo.
Embora possamos executar uma regress&aatilde;o de Poisson em R usando a função glm em um dos pacotes principais, precisamos de outro pacote para executar o modelo de poisson inflado de zeros. Usamos o pacote pscl.
A saída se parece muito com a saída de duas regressões OLS em R. Abaixo da chamada do modelo, você encontrará um bloco de saída contendo os coeficientes da regressão de Poisson para cada uma das variáveis, juntamente com erros padrão, escores \(z\) e os \(p\)-valores para os coeficientes. Segue-se um segundo bloco que corresponde ao modelo de inflação. Isso inclui coeficientes logit para prever zeros em excesso junto com seus erros padrão, escores \(z\) e \(p\)-valores.
Todos os preditores nas partes de contagem e inflação do modelo são estatisticamente significativos. Este modelo se ajusta aos dados significativamente melhor do que o modelo nulo, ou seja, o modelo somente de intercepto. Para mostrar que esse é o caso, podemos comparar com o modelo atual a um modelo nulo sem preditores usando o teste qui-quadrado da diferença de log-verossimilhanças.
Como temos três variáveis preditoras no modelo completo, os graus de liberdade para o teste qui-quadrado são 3. Isso resulta em um \(p\)-valor significativo alto; portanto, nosso modelo geral é estatisticamente significativo.
Observe que a saída do modelo acima não indica de forma alguma se nosso modelo inflado de zero é uma melhoria em relação a uma regressão de Poisson padrão. Podemos determinar isso executando o modelo de Poisson padrão correspondente e, em seguida, executando um teste de Vuong dos dois modelos.
O teste de Vuong compara o modelo inflado de zero com um modelo de regressão de Poisson comum. Neste exemplo, podemos ver que nossa estatística de teste é significativa, indicando que o modelo inflado de zero é superior ao modelo de Poisson padrão.
Existem algumas críticas e debates na literatura sobre o uso e mau uso do teste de Vuong para modelos não aninhados no teste de inflação zero em modelos de contagem. O objetivo aqui é mais focado em como implementar várias análises em R. Verifique as referências abaixo para maiores informações:
Podemos obter intervalos de confiança para os parâmetros e os parâmetros exponenciados de diversas maneiras, aqui usando bootstrapping. Para o modelo de Poisson, esses seriam os índices de risco incidente; para o modelo de inflação zero, os odds ratios.
Usamos o pacote boot. Primeiro, obtemos os coeficientes de nosso modelo original para usar como valores iniciais no modelo e acelerar o tempo que leva para estimar. Em seguida, escrevemos uma função curta que recebe dados e índices como entrada e retorna os parâmetros nos quais estamos interessados. Por fim, passamos isso para a função boot e fazemos 1200 replicações. Para resultados finais, pode-se desejar aumentar o número de replicações para ajudar a garantir resultados estáveis.
Os resultados são estimativas alternadas dos parâmetros e erros padrão. Ou seja, a primeira linha contém a primeira estimativa de parâmetro de nosso modelo. O segundo tem o erro padrão do primeiro parâmetro. A terceira coluna contém os erros padrão inicializados, que são consideravelmente maiores do que aqueles estimados por zeroinfl.
Agora podemos obter os intervalos de confiança para todos os parâmetros. Começamos na escala original com intervalos de confianças (ICs) ajustados por percentil e viés. Também comparamos esses resultados com os intervalos de confiança regulares com base nos erros padrão.
Os intervalos de confiança bootstrapped são consideravelmente maiores do que a aproximação baseada na normal. Agora podemos estimar a razão de risco de incidente (TIR) para o modelo de Poisson e a razão de chances (OR) para o modelo logístico (inflação zero). Isso é feito usando um código quase idêntico ao anterior, mas passando uma função de transformação para o argumento \(h\) de boot.ci, neste caso, \(\exp\) para exponenciar.
Para entender melhor nosso modelo, podemos calcular o número esperado de peixes capturados para diferentes combinações de nossos preditores. Na verdade, como estamos trabalhando com preditores essencialmente categóricos, podemos calcular os valores esperados para todas as combinações usando a função expand.grid para criar todas as combinações e, em seguida, a função de previsão para fazer isso. Também removemos todas as linhas em que o número de filhos excede o número de pessoas, o que não faz sentido logicamente, usando a função subset. Finalmente criamos um gráfico.
Outras referências:
A distribuição conjunta de várias variáveis categóricas define uma tabela de contingência. Se uma das variáveis em uma tabela de contingência for tratada como a variável de resposta, podemos ajustar um modelo logit ou probit, ou seja, para uma resposta dicotómica, um GLM binomial à tabela. Os modelos loglineares, em contraste, que são modelos para as associações entre as variáveis em uma tabela de contingência, tratam as variáveis simetricamente - eles não distinguem uma variável como a resposta. Há, entretanto, uma relação entre os modelos loglinear e os modelos logit que desenvolveremos posteriormente nesta seção. Como veremos também, os modelos loglineares têm a estrutura formal de modelos ANOVA de dois fatores e fatores superiores e podem ser ajustados aos dados por regressão de Poisson.
Os modelos loglineares para tabelas de contingência têm muitas aplicações especializadas nas ciências sociais - por exemplo, para quadrar tabelas, como tabelas de mobilidade, em que as variáveis na tabela têm as mesmas categorias. O tratamento de modelos loglineares nesta seção apenas arranha a superfície. Relatos mais extensos estão disponíveis em muitas fontes, incluindo Agresti (1990), Fienberg (1980) e Powers and Xie (2008).
Examinaremos as tabelas de contingência para duas variáveis com algum detalhe, pois este é o caso mais simples e os principais resultados que estabelecemos aqui se estendem diretamente a tabelas de dimensão superior. Considere a tabela bidirecional ilustrativa mostrada abaixo, construída a partir de dados relatados no American Voter (Campbell, Converse, Miller & Stokes, 1960). A tabela relaciona a intensidade da preferência partidária à participação na vota¸ão nas eleições presidenciais dos EUA em 1956. Para antecipar a análise, os dados indicam que o comparecimento às urnas está positivamente associado à intensidade da preferência partidária.
| Participação | eleitoral | |||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Intensidade de preferência | Votou | Não votou | Total | |||||||||||||||||||
| Fraca | 305 | 126 | 431 | |||||||||||||||||||
| Média | 405 | 125 | 530 | |||||||||||||||||||
| Forte | 265 | 49 | 314 | |||||||||||||||||||
| Total | 975 | 300 | 1275 |
De forma mais geral, duas variáveis categóricas com \(r\) e \(c\) categorias, respectivamente, definem uma tabela de contingência \(r\times c\), onde \(Y_{ij}\) é a contagem de frequência observada na célula \((i,j)\)-ésima da tabela. Assim, \(Y_{i+} = \sum_{j=1}^c Y_{ij}\) é a frequência marginal na \(i\)-ésima linha; \(Y_{+j} = \sum_{i=1}^r Y_{ij}\) é a frequência marginal na \(j\)-ésima coluna e \(n = Y_{++} = \sum_{i=1}^c \sum_{j=1}^r Y_{ij}\) é o número de observações na amostra.
Supomos que as \(n\) observações na tabela são independentemente amostradas de uma população com proporção \(\pi_{ij}\) na célula \((i,j)\) e, portanto, que a probabilidade de amostragem de uma observação individual nesta célula é \(\pi_{ij}\). As distribuições de probabilidade marginais \(\pi_{i+}\) e \(\pi_{+j}\) podem ser definidas como acima; observe que \(\pi_{++} = 1\). Se as variáveis de linha e coluna são estatisticamente independentes na população, então a probabilidade conjunta \(\pi_{ij}\) é o produto das probabilidades marginais para todo \(i\) e \(j\): \(\pi_{ij} = \pi_{i+}\times \pi_{+j}\).
Como as frequências observadas \(Y_{ij}\) resultam do desenho de uma amostra aleatória, elas são variáveis aleatórias que geralmente assumem valores diferentes em amostras diferentes. Então, a frequência esperada na célula \((i,j)\) é \(\mu_{ij}=\mbox{E}(Y_{ij}) = n\pi_{ij}\).
Se as variáveis são independentes, temos que \(\mu_{ij} = n\pi_{i+}\times \pi_{+j}\). Além disso, como \(\mu_{i+} = \sum_{j = 1}^c n\pi_{ij} = n\pi_{i+}\) e \(\mu_{+j} = \sum_{i=1}^r n\pi_{ij} = n\pi_{+j}\), podemos escrever \(\mu_{ij} = \mu_{i+}\times \mu_{+j}/n\).
Tomando o logaritmo de ambos os lados desta última equação produz \begin{equation*} \eta_{ij}=\log_e(\mu_{ij})=\log_e(\mu_{i+})+\log_e(\mu_{+j})-\log_e(n)\cdot \end{equation*}
Ou seja, sob independência, o logaritmo das frequências esperadas \(\eta_{ij}\) depende aditivamente dos logaritmos das frequências marginais esperadas da linha, das frequências marginais esperadas das colunas e do tamanho da amostra. A expressão acima é uma reminiscência de um modelo ANOVA bidirecional de efeitos principais, onde − \(\log_e(n)\) desempenha o papel da constante, \(\log_e(\mu_{i+})\) e \(\log_e(\mu_{+j})\) são análogos a parâmetros de “efeito principal” e \(\eta_{ij}\) aparece no lugar da média da variável de resposta.
Se impormos restrições ao modelo ANOVA, podemos reparametrizar a equação acima da seguinte forma: \begin{equation*} \eta_{ij}=\mu+\alpha_i+\beta_j, \end{equation*} onde \(\alpha_+=\sum \alpha_i = 0\) e \(\beta_+ =\sum\beta_j = 0\). A equação acima é o modelo loglinear para independência na tabela bidirecional. Resolvendo para os parâmetros do modelo, obtemos \begin{array}{rcl} \mu & = & \dfrac{\eta_{++}}{r c}, \\ \alpha_i & = & \dfrac{\eta_{i+}}{c}-\mu, \\ \beta_j & = & \dfrac{\eta_{+j}}{r}-\mu\cdot \end{array}
É importante ressaltar que embora o modelo loglinear seja formalmente semelhante a um modelo ANOVA, o significado dos dois modelos difere muito: Na análise de variância, \(\alpha_i\) e \(\beta_j\) são parâmetros de efeito principal, especificando a relação parcial da variável de resposta (quantitativa) com cada variável explicativa. O modelo loglinear, em contraste, não distingue uma variável de resposta e, por ser um modelo para independência, especifica que as variáveis de linha e coluna na tabela de contingência não estão relacionadas; para este modelo, o \(\alpha_i\) e o \(\beta_j\) apenas expressam a relação das frequências logarítmicas esperadas das células com os marginais da linha e da coluna.
O modelo de independência descreve \(rc\) frequências esperadas em termos de \begin{equation*} 1+(r-1)+(c-1)=r+c-1 \end{equation*} parâmetros independentes.
Por analogia ao modelo ANOVA bidirecional, podemos adicionar parâmetros para estender o modelo loglinear aos dados para os quais as classificações de linha e coluna não são independentes na população, mas sim relacionadas de maneira arbitrária: \begin{equation*} \eta_{ij} = \mu+\alpha_i+\beta_j+\gamma_{ij}, \end{equation*} onde \(\alpha_+ = \beta_+ = \gamma_{i+} = \gamma_{+j} = 0\) para todos \(i\) e \(j\). Como antes, podemos escrever os parâmetros do modelo em termos do logaritmo de contagens esperadas \(\eta_{ij}\). De fato, a solução para \(\mu\), \(\alpha_i\) e \(\beta_j\) é a mesma que acima encontrada e \begin{equation*} \gamma_{ij}=\eta_{ij}-\mu-\alpha_i-\beta)j \end{equation*}
Nesta seção, revisitamos com maior rigor e mais detalhes muitos dos pontos levantados nas seções anteriores. A exposição aqui deve ao Capítulo 2 de McCullagh e Nelder (1989), que se tornou a fonte padrão sobre modelos lineares generalizados e ao tratamento mais breve e notavelmente lúcido e perspicaz do tópico por Firth (1991).
Como muito mais nas estatísticas modernas, a perspicácia de que muitas das distribuições mais importantes na estatística poderiam ser expressas na seguinte forma linear-exponencial comum foi devido a R.A. Fisher: \begin{equation} p(y;\theta,\phi) \, = \, \exp\left( \dfrac{y\theta-b(\theta)}{a(\phi)}+c(y,\phi)\right) \end{equation} onde
Considere, por exemplo, a distribuição normal ou gaussiana com média \(\mu\) e variância \(\sigma^2\). Colocar a distribuição normal na forma da equação acima requer alguma manipulação algébrica, eventualmente produzindo \begin{equation} p(y;\theta,\phi) \, = \, \exp\left( \dfrac{y\theta-\theta^2/2}{\phi}-\dfrac{1}{2}\left(\dfrac{y^2}{\phi}+\log_e\big(2\pi\phi\big)\right)\right), \end{equation} onde \(\theta=g_c(\mu)=\mu\), \(\phi=\sigma^2\), \(a(\phi)=\phi\), \(b(\theta)=\theta^2/2\) e \(c(y,\phi)=-\frac{1}{2}\Big( y^2/\phi+\log_e\big(2\pi\phi\big)\Big)\).
Agora considere a distribuição binomial, onde \(Y\) é a proporção de sucessos em \(n\) tentativas binárias independentes e \(\mu\) é a probabilidade de sucesso em uma tentativa individual. Escrito depois da ginástica algébrica como uma família exponencial \begin{equation} p(y;\theta,\phi) \, = \, \exp\left( \dfrac{y\theta-\log_e\big(1+e^\theta\big)}{1/n}+\log_e{n \choose ny} \right), \end{equation} onde \(\theta=g_c(\mu)=\log_e\big(\mu/(1-\mu)\big)\), \(\phi=1\), \(a(\phi)=1/n\), \(b(\theta)=\log_e\big(1+e^\theta\big)\) e \(c(y,\phi)=\log_e{n \choose ny}\).
Da mesma forma, as famílias Poisson, gama e Gaussiana inversa podem ser colocadas na forma linear-exponencial acima, usando os resultados dados na Tabela 9.
| Família | \(a(\phi)\) | \(b(\theta)\) | \(c(y,\phi)\) | |||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Gaussiana | \(\phi\) | \(\theta^2/2\) | \(-\frac{1}{2}\Big( y^2/\phi+\log_e\big(2\pi\phi\big)\Big)\) | |||||||||||||||||||
| Binomial | \(1/n\) | \(\log_e\big(1+e^\theta\big)\) | \(\log_e{n \choose ny}\) | |||||||||||||||||||
| Poisson | \(1\) | \(e^\theta\) | \(-\log_e(y!)\) | |||||||||||||||||||
| Gama | \(\phi\) | \(-\log_e(-\theta)\) | \(\phi^{-2}\log_e(y/\phi)-\log_e(y)-\log_e\big(\Gamma(\phi^{-1})\big)\) | |||||||||||||||||||
| Gaussiana invera | \(\phi\) | \(-\sqrt{-2\theta}\) | \(-\frac{1}{2}\Big( \log_e\big(\pi\phi y^3\big)+(\phi y)^{-1}\Big)\) |
NOTA: Nesta tabela, \(n\) é o número de observações binomiais e \(\Gamma(\cdot)\) é a função gama.
A vantagem de expressar diversas famílias de distribuições na forma exponencial comum é que as propriedades gerais das famílias exponenciais podem então ser aplicadas aos casos individuais. Por exemplo, é verdade em geral que \begin{equation} b'(\theta) \, = \, \dfrac{\mbox{d}b(\theta)}{\mbox{d}\theta} \, = \, \mu \end{equation} e que \begin{equation} V(Y) \, = \, a(\phi)b''(\theta) \, = \, a(\phi)\dfrac{\mbox{d}^2b(\theta)}{\mbox{d}\theta^2} \, = \, a(\phi)\nu(\mu) \end{equation} levando aos resultados na Tabela 2. Observe que \(b(\cdot)\) é o inverso da função de ligação caônica. Por exemplo, para a distribuição normal, \begin{array}{rcl} b'(\theta) & = & \dfrac{\mbox{d}(\theta^2/2)}{\mbox{d}\theta} \, = \, \theta \, = \, \mu, \\ a(\phi)b''(\theta) & = & \phi\times 1 \, = \, \sigma^2, \\ \nu(\mu) & = & 1 \end{array} e para a distribuição binomial, \begin{array}{rcl} b'(\theta) & = & \dfrac{\mbox{d}\log_e\big(1+e^\theta\big)}{\mbox{d}\theta} \, = \, \dfrac{e^\theta}{1+e^\theta} \, = \, \dfrac{1}{1+e^{-\theta}} \, = \, \mu, \\ a(\phi)b''(\theta) & = & \dfrac{1}{n}\left(\dfrac{e^\theta}{1+e^\theta}-\Big(\dfrac{e^\theta}{1+e^\theta}\Big)^2\right) \, = \, \dfrac{\mu(1-\mu)}{n}, \\ \nu(\mu) & = & \mu(1-\mu)\cdot \end{array}
O logaritmo da função de verossimilhança para uma observação individual \(Y_i\) assume a forma \begin{equation} \log_e\big( L(\theta_i,\phi; y_i)\big) \, = \, \dfrac{y_i\theta_i-b(\theta_i)}{a_i(\phi)} + c(y_i,\phi)\cdot \end{equation} Para \(n\) observações independentes, temos \begin{equation} \log_e\big( L(\theta,\phi; y)\big) \, = \, \sum_{i=1}^n \left(\dfrac{y_i\theta_i-b(\theta_i)}{a_i(\phi)} + c(y_i,\phi)\right), \end{equation} onde \(\theta=(\theta_1,\theta_2\cdots,\theta_n)=\{\theta_i\}_{i=1}^n\) e \(y=(y_1,y_2,\cdots,y_n)=\{y_i\}_{i=1}^n\).
Suponha que um modelo linear generalizado use a função de ligação \(g(\cdot)\), é notacionalmente conveniente escrever \(\beta_0\) para a constante de regressão \(\alpha\); de modo que, \begin{equation} g(\mu_i) \, = \, \eta_i \, = \, \beta_0+\beta_1 X_{i1}+\beta_2 X_{i2}+\cdots+\beta_k X_{ik}\cdot \end{equation}
O modelo, portanto, expressa os valores esperados das \(n\) observações em termos de um número muito menor de parâmetros de regressão. Para obter equações de estimação para os parâmetros de regressão, temos que diferenciar o logaritmo da verossimilhança em relação a cada coeficiente por sua vez. Consideremos que \(\ell_i\) represente a \(i\)-ésima componente da log verossimilhança. Então, pela regra da cadeia, \begin{equation} \dfrac{\partial\ell_i}{\partial\beta_j} \, = \, \dfrac{\partial \ell_i}{\partial\theta_i}\times \dfrac{\mbox{d}\theta_i}{\mbox{d}\mu_i} \times \dfrac{\mbox{d}\mu_i}{\mbox{d}\eta_i}\times \dfrac{\partial\eta_i}{\partial\beta_j}, \qquad \mbox{para} \qquad j=0,1,\cdots,k\cdot \end{equation}
Depois de algum trabalho, podemos reescrever a equação acima como \begin{equation} \dfrac{\partial\ell_i}{\partial\beta_j} \, = \, \dfrac{y_1-\mu_i}{a_i(\phi)\nu(\mu_i)}\times \dfrac{\mbox{d}\mu_i}{\mbox{d}\eta_i}\times x_{ij}\cdot \end{equation} Somando as observações e definindo a soma sendo zero, produz as equações de estimação de máxima verossimilhança para os modelos lineares generalizados, \begin{equation} \sum_{i=1}^n \dfrac{y_i-\mu_i}{a_i\nu(\mu_i)}\times \dfrac{\mbox{d}\mu_i}{\mbox{d}\eta_i} \times x_{ij} \, = \, 0, \qquad \mbox{para} \qquad j=0,1,\cdots,k, \end{equation} onde \(a_i=a_i(\phi)/\phi\) não depende do parâmetro de dispersão, que é constante nas observações. Por exemplo, em um modelo linear generalizado gaussiano, \(a_i = 1\), enquanto em um modelo linear generalizado binomial, \(a_i = 1/n_i\).
Simplificação adicional pode ser alcançada quando \(g(\cdot)\) é a ligação canônica. Neste caso, as equações de estimação de máxima verossimilhança tornam-se \begin{equation} \sum_{i=1}^n \dfrac{y_i x_{ij}}{a_i} \, = \, \sum_{i=1}^n \dfrac{\mu_i x_{ij}}{a_i}, \end{equation} definindo a soma observada à esquerda da equação para a soma esperada à direita. Observamos esse padrão nas equações de estimação para modelos de regressão logística anteriormente. No entanto, mesmo aqui as equações de estimação são, exceto no caso da família Gaussiana, emparelhados com o ligação de identidade, funções não lineares dos parâmetros de regressão e geralmente requerem métodos iterativos para sua solução.
Seja \begin{equation} Z_i \, = \, \eta_i +(y_i-\mu_i)\dfrac{\mbox{d}\eta_i}{\mbox{d}\mu_i} \, = \, \eta_i+(y_i-\mu_i)g'(\mu_i)\cdot \end{equation} Então \begin{equation} \mbox{E}(Z_i) \, = \, \eta_i \, = \, \beta_0+\beta_1 X_{i1}+\beta_2 X_{i2}+\cdots+\beta_k X_{ik}, \end{equation} e \begin{equation} \mbox{Var}(Z_i) \, = \, \big( g'(\mu_i)\big)^2 a_i \nu(\mu_i)\cdot \end{equation}
Se, portanto, pudéssemos calcular o \(Z_i\), seríamos capazes de ajustar o modelo por regressão de mínimos quadrados ponderados de \(Z\) nos \(X\)s, usando os inversos da \(\mbox{Var}(Z_i)\) como pesos. Claro, este não é o caso porque não sabemos os valores de \(\mu_i\) e \(\eta_i\), que, de fato, dependem dos coeficientes de regressão que desejamos estimar - ou seja, o argumento é essencialmente circular. Esta observação sugeriu a Nelder and Wedderburn (1972) a possibilidade de estimar os modelos lineares generalizados por mínimos quadrados ponderados iterativos (IWLS), habilmente transformando a circularidade em um procedimento iterativo:
Aplicado à ligação canônica, o IWLS é equivalente ao método Newton-Raphson, como descobrimos para um modelo logit anteriormebe; mais geralmente, IWLS implementa o método de Escore de Fisher.
Observe que não exigimos uma estimativa do parâmetro de dispersão para estimar os coeficientes de regressão em um modelo linear generalizado. Embora seja, em princípio, possível estimar \(\phi\) também por máxima verossimilhança, isso raramente é feito. Em vez disso, lembre-se de que \(\mbox{Var}(Y_i) = \phi a_i \nu(\mu_i)\). Resolvendo para o parâmetro de dispersão, obtemos \(\phi = \mbox{Var}(Y_i)/a_i\nu(\mu_i)\), sugerindo o método do estimador de momentos \begin{equation} \widetilde{\phi} \, = \, \dfrac{1}{n-k-1}\sum_{i=1}^n \dfrac{\big(y_i-\widehat{\mu}_i \big)^2}{a_i \nu(\widehat{\mu}_i)}\cdot \end{equation}
A matriz de covariância assintótica estimada dos coeficientes é então obtida a partir da última iteração IWLS como \begin{equation} \widehat{\mbox{Var}}(\widehat{\beta}) \, = \, \widetilde{\phi}\big({\bf X}^\top {\bf W}{\bf X} \big)^{-1}\cdot \end{equation} Como o estimador de máxima verossimilhança \(\widehat{\beta}\) é normalmente distribuído assintoticamente, \(\widehat{\mbox{Var}}(\widehat{\beta})\) pode ser usada como base para os testes Wald dos parâmetros de regressão.
O argumento que leva à estimativa de IWLS repousa apenas na linearidade da relação entre \(\eta = g(\mu)\) e os \(X\)s e na suposição de que \(\mbox{Var}(Y)\) depende de uma maneira particular do parâmetro de dispersão e \(\mu\). Contanto que possamos expressar a média transformada de \(Y\) como uma função linear dos \(X\)s, e podemos escrever uma função de variância para \(Y\), expressando a variância condicional de \(Y\) como uma função de sua média e um parâmetro de dispersão, podemos aplicar as equações de estimação de máxima verossimilhança e obter estimativas por IWLS - mesmo sem nos comprometermos com uma distribuição condicional particular para \(Y\).
Este é o método de estimação de quase-verossimilhança, introduzido por Wedderburn (1974), e mostrou reter muitas das propriedades da estimação de máxima-verossimilhança: Embora o estimador de quase-verossimilhança possa não ser maximamente assintoticamente eficiente, é consistente e tem a mesma distribuição assintótica que o estimador de máxima verossimilhança de um modelo linear generalizado em uma família exponencial. Podemos pensar na estimativa de quase verossimilhança de modelos lineares generalizados como análoga à estimativa de mínimos quadrados de modelos de regressão lineares com erros potencialmente não normais: Lembre-se de que, enquanto a rela¸ão entre \(Y\) e os \(X\)s for linear, a variância do erro é constante e as observações são amostradas de forma independente, a teoria subjacente à estimativa de OLS se aplica - embora o estimador de OLS possa não ser mais maximamente eficiente.
Originalmente, escrevemos o logaritmo da verossimilhança para um modelo linear generalizado como uma função \(\log_e\big(L(\theta,\phi; y)\big)\) do parâmetro canônico \(\theta\) para as observações. Como \(\mu_i = g_c^{−1}(\theta_i)\), para o elo canônico \(g_c(\cdot)\), podemos igualmente pensar na probabilidade logarítmica como uma função da resposta esperada e, portanto, podemos escrever a verossimilhança maximizada como \(\log_e\big(L(\widehat{\mu},\phi; y)\big)\). Se, então, dedicarmos um parâmetro a cada observação, de modo que \(\widehat{\mu}_i = y_i\), por exemplo, removendo a constante do modelo de regressão e definindo um regressor fictício para cada observação, o log-verossimilhança torna-se \(\log_e\big(L(y,\phi;y)\big)\).
O desvio residual sob o modelo inicial é duas vezes a diferença nessas log-verossimilhanças: \begin{array}{rcl} D(y;\widehat{\mu}) & = & \displaystyle 2\left( \log_e\big(L(y,\phi;y)\big)-\log_e\big(L(\widehat{\mu},\phi; y)\big)\right) \\[0.6em] & = & \displaystyle 2\sum_{i=1}^n \left( \log_e\big(L(y_i,\phi;y_i)\big)-\log_e\big(L(\widehat{\mu},\phi; y_i)\big)\right) \\[0.6em] & = & \displaystyle 2\sum_{i=1}^n \dfrac{1}{a_i}\Big(y_i\big(g(y_i)-g(\widehat{\mu}_i)\big)-b\big(g(y_i)\big)+b\big(g(\widehat{\mu}_i)\big) \Big)\cdot \end{array} A divisão do desvio residual pelo parâmetro de dispersão estimado produz o desvio escalonado, \begin{equation} D^*(y;\widehat{\mu}) = \dfrac{1}{\widetilde{\phi}}D(y;\widehat{\mu})\cdot \end{equation} Conforme explicado na Seço I.1, os desvios são os blocos de construção do teste da razão de verossimilhanças e dos testes \(F\) para GLMs.
Aplicando a equação acima à distribuição Gaussiana, onde \(g_c(\cdot)\) é o elo de identidade, \(a_i = 1\) e \(b(\theta) = \theta^2/2\), produz, após alguma simplificação \begin{equation} D(y;\widehat{\mu}) \, = \, \sum_{i=1}^n \big(y_i-\widehat{\mu}_i \big)^2 , \end{equation} ou seja, a soma residual dos quadrados do modelo. Da mesma forma, aplicando a mesma equação à distribuição binomial, onde \(g_c(\cdot)\) é a ligação logit, \(a_i = n_i\) e \(b(\theta) = \log_e\big( 1 + e^\theta\big)\), obtemos, após um pouco de simplificação \begin{equation} D(y;\widehat{\mu}) \, = \, \sum_{i=1}^n n_i\left(y_i\log_e\Big(\dfrac{y_i}{\widehat{\mu}_i}\Big)+(1-y_i)\log_e\Big(\dfrac{1-y_i}{1-\widehat{\mu}_i}\Big) \right)\cdot \end{equation} Deixamos como exercício para o leitor o desenvolvimento de fórmulas para o desvio nos modelos de Poisson, gama e Gaussiano inverso.
Como no caso dos modelos lineares, podemos formular um teste para a hipótese linear geral \begin{equation} H_0 \, : \, {\bf L}_{q\times (k+1)}\beta_{(k+1)\times 1} \, = \, {\bf c}_{q\times 1}, \end{equation} onde a matriz de hipótese \({\bf L}\) e o vetor \({\bf c}\) do lado direito contêm constantes pré-especificadas; geralmente, \({\bf c = 0}\). Para um modelo linear generalizado, a estatística Wald \begin{equation} Z_0^2 \, = \, ({\bf L}\widehat{\beta}-{\bf c})^\top \big({\bf L}\widehat{\mbox{Var}}(\widehat{\beta}){\bf L}^\top\big)^{-1} ({\bf L}\widehat{\beta}-{\bf c}) \end{equation} segue uma distribuição assintótica qui-quadrado com \(q\) graus de liberdade sob a hipótese. A aplicação mais simples desse resultado é a estatística Wald \(Z_0 = \widehat{\beta}_j/SE(\widehat{\beta}_j)\), testando se um coeficiente de regressão individual é zero. Aqui, \(Z_0\) segue uma distribuição normal padrão sob \(H_0 \, : \, \beta_j = 0\) ou, equivalentemente, \(Z_0^2\) segue uma distribuição qui-quadrado com um grau de liberdade.
Alternativamente, quando o parâmetro de dispersão é estimado a partir dos dados, podemos calcular a estatística de teste \begin{equation} F_0 \, = \, \dfrac{1}{q}Z_0^2, \end{equation} que é distribuída como \(F_{q, n−k−1}\) sob \(H_0\). Aplicado a um coeficiente individual, \(t_0 = \pm \sqrt{F_0} = \widehat{\beta}_j/SE(\widehat{\beta}_j)\) produz um teste \(t\) com \(n-k-1\) graus de liberdade.
Ocasionalmente, é interessante testar uma hipótese ou construir um intervalo de confiança para uma função não linear dos parâmetros de um modelo linear ou linear generalizado. Se a função não linear em questão é uma função diferenciável dos coeficientes de regressão, um erro padrão assintótico aproximado pode ser obtido pelo método delta.
O método delta (Rao, 1973) emprega uma aproximação da série de Taylor de primeira ordem, isto é, linear para a função não linear. O método delta é apropriado aqui porque as estimativas de máxima verossimilhança ou quase verossimilhança dos coeficientes de um GLM são normalmente distribuídas assintoticamente.
Na verdade, o procedimento descrito nesta seção é aplicável sempre que os parâmetros de um modelo de regressão sejam normalmente distribuídos e podem, portanto, ser aplicados em uma ampla variedade de contextos. Em pequenas amostras, no entanto, a aproximação do método delta para o erro padrão pode não ser adequada e os procedimentos de bootstrapping geralmente fornecerão resultados mais confiáveis.
Suponha que estejamos interessados na função \begin{equation} \gamma \, = \, f(\beta) \, = \, f(\beta_0,\beta_1,\cdots,\beta_k), \end{equation} onde, por conveniência de notação, usamos \(\beta_0\) para denotar a constante de regressão. A função \(f(\beta)\) não precisa usar todos os coeficientes de regressão, veja o exemplo abaixo. O estimador de máxima verossimilhança de \(\gamma\) é simplesmente \(\widehat{\gamma} = f(\widehat{\beta})\) que, como um estimador de máxima verossimilhança, também é assintoticamente normal e a variância amostral aproximada de \(\gamma\) é então \begin{equation} \widehat{\mbox{Var}}(\widehat{\gamma}) \approx \sum_{j=0}^k \sum_{l=0}^k \widehat{\sigma}_{jl} \dfrac{\partial \widehat{\gamma}}{\partial\widehat{\beta}_j} \dfrac{\partial \widehat{\gamma}}{\partial\widehat{\beta}_l}, \end{equation} onde \(\widehat{\sigma}_{jl}\) é o \((j,l)\)-ésimo elemento da matriz de covariância assintótica estimada dos coeficientes, \(\widehat{\mbox{Var}}(\widehat{\beta})\).
Para ilustrar a aplicação desse resultado, imagine que estamos interessados em determinar o valor máximo ou mínimo de uma regressão parcial quadrática. Enfocando a relação parcial entre a variável resposta e um \(X\) em particular, temos uma equação da forma \begin{equation} \mbox{E}(Y) \, = \, \cdots + \beta_1 X+\beta_2 X^2+\cdots\cdot \end{equation} A aplicação do método delta para encontrar o mínimo ou máximo de uma curva quadrática foi sugerida por Weisberg (2005).
Diferenciando esta equação em relação a \(X\), obtemos \begin{equation} \dfrac{\mbox{d}\mbox{E}(Y)}{\mbox{d}X} \, = \, \beta_1+2\beta_2 X\cdot \end{equation} Definindo a derivada como 0 e resolvendo para \(X\) produz o valor no qual a função atinge um mínimo, se \(\beta_2\) for positivo ou um máximo se \(\beta_2\) for negativo, \begin{equation} X \, = \, -\dfrac{\beta_1}{2\beta_2}, \end{equation} que é uma função não linear dos coeficientes de regressão \(\beta_1\) e \(\beta_2\).
Por exemplo, usando dados da Pesquisa Canadense de Dinâmica de Trabalho e Renda ou SLID ajustamos uma regressão de mínimos quadrados do logaritmo base 2 da taxa de salário segundo a idade quadrática, um regressor fictício para sexo e o quadrado dos anos de estudo, obtendo:
Imagine que estamos interessados na idade \(\gamma = −\beta_1/(2\beta_2)\) em que os salários estão no máximo, mantendo o sexo e a educação constantes. As derivadas necessárias são \begin{array}{rcl} \left.\dfrac{\mbox{d}\widehat{\gamma}}{\mbox{d}\beta_1}\right|_{\beta=\widehat{\beta}} & = & -\dfrac{1}{2\widehat{\beta}_2} \, = \, -\dfrac{1}{2(-0.001230)} \, = \, 406.5 \\[0.6em] \left.\dfrac{\mbox{d}\widehat{\gamma}}{\mbox{d}\beta_2}\right|_{\beta=\widehat{\beta}} & = & -\dfrac{\widehat{\beta}_1}{2\widehat{\beta}_2^2} \, = \, -\dfrac{0.1198}{2(-0.001230)^2} \, = \, 39.593 \end{array} Nossa estimativa pontual de \(\gamma\) é \begin{equation} \widehat{\gamma} \, = \, -\dfrac{\widehat{\beta}_1}{2\widehat{\beta}_2} \, = \, -\dfrac{0.1198}{2\times 0.001230} \, = \, 48.70 \; \mbox{anos}\cdot \end{equation}
A variância amostral estimada do coeficiente de idade é \(\widehat{\mbox{Var}}(\widehat{\beta}_1) = 2.115\times 10^{−5}\), e do coeficiente de idade ao quadrado, \(\widehat{\mbox{Var}}(\widehat{\beta}_2) = 3.502\times 10^{−9}\); a covariância amostral estimada para os dois coeficientes é \begin{equation} \widehat{\mbox{Cov}}(\widehat{\beta}_1, \widehat{\beta}_2) \, = \, −2.685\times 10^{−7}\cdot \end{equation} A variância estimada aproximada de \(\gamma\) é então \begin{array}{cl} \widehat{\mbox{Var}}(\widehat{\gamma}) & \approx & \big(2.115\times 10^{-5}\big)\times 406.5^2-\big(2.685\times 10^{-7}\big)\times 406.5 \times 39.593 \\[0.4em] & & \qquad -\big(2.685\times 10^{-7}\big)\times 406.5\times 39.593+\big(3.502\times 10^{-9}\big)\times 39.593^2 \, = \, 0.3419\cdot \end{array} Consequentemente, o erro padrão aproximado de \(\gamma\) é \(SE(\gamma) \approx \sqrt{0.3419} = 0.5847\), e um intervalo de confiança de aproximadamente 95% para a idade em que a renda é mais alta em média é \begin{equation} \gamma \in 48.70 \pm 1.96 (0.5847) = (47.55, 49.85)\cdot \end{equation}
Vamos escrever o modelo linear generalizado em forma matricial, com preditor linear \begin{equation} \eta_{n\times 1} \, = \, {\bf X}_{n\times (k+1)}\beta_{(k+1)\times 1} \end{equation} e função de ligação \(g(\mu) = \eta\), onde \(\mu\) é a esperança do vetor de resposta \(Y\).
Conforme descrito anteriormente, calculamos o estimador de máxima verossimilhança \(\widehat{\beta}\) de \(\beta\), juntamente com o estimador da matriz de covariância assintótica \(\widehat{\mbox{Var}}(\widehat{\beta})\) de \(\widehat{\beta}\).
Considere que as linhas de \({\bf X}^*\) incluam regressores correspondentes a todas as combinações de valores de variáveis explicativas que aparecem em um termo de ordem superior do modelo ou, para uma variável explicativa contínua, valores abrangendo o intervalo da variável, junto com valores típicos das restantes regressores. A estrutura de \({\bf X}^*\) com respeito às interações, por exemplo, é a mesma que a da matriz do modelo \({\bf X}\). Então os valores ajustados \begin{equation} \widehat{\eta}^* \, = \, {\bf X}^*\widehat{\beta}, \end{equation} representam o termo de ordem superior em questão e uma tabela ou gráfico destes valores - ou, alternativamente, dos valores ajustados transformados para a escala da variável resposta, \(g^{−1}(\widehat{\eta}^*)\) - é uma exibição de efeito. Os erros padrão de \(\widehat{\eta}^*\), disponíveis como entradas diagonais de raiz quadrada de \({\bf X}^*\widehat{\mbox{Var}}(\widehat{\beta}){{\bf X}^*}^\top\), podem ser usados para calcular intervalos de confiança pontuais para os efeitos, cujos pontos finais também podem ser transformados na escala da resposta.
Os diagnósticos de regressão são métodos para determinar se um modelo de regressão ajustado representa adequadamente os dados. A maioria dos diagnósticos para modelos lineares se estendem de forma relativamente diretamente aos GLMs. Essas extensões normalmente tiram vantagem do cálculo dos estimadores de máxima verossimilhança e máxima quase-verossimilhança para GLMs obtidos por mínimos quadrados ponderados iterados, conforme descrito na Seção III.2. O ajuste final por mínimos quadrados ponderados lineariza o modelo e fornece uma aproximação quadrática para o logaritmo da verossimilhança. Os diagnósticos aproximados são baseados diretamente na solução WLS ou são derivados de estatísticas facilmente calculadas a partir desta solução. O trabalho sobre a extensão de diagnósticos de mínimos quadrados lineares para GLMs foi feito por Pregibon (1981), Landwehr, Pregibon e Shoemaker (1984), Wang (1985, 1987) e Williams (1987).
Modelos lineares ajustados por mínimos quadrados fazem suposições fortes e às vezes irrealistas sobre a estrutura dos dados. Quando essas premissas são violadas, as estimativas de mínimos quadrados podem se comportar mal e podem até representar os dados de maneira completamente incorreta. Os diagnósticos de regressão podem revelar esses problemas e, muitas vezes, apontar o caminho para as soluções.
Todos os métodos discutidos estão disponíveis nas funções R padrão ou são implementados no pacote do car. Um dos objetivos do pacote car é fazer diagnósticos para modelos lineares e GLMs prontamente disponíveis em R. Nossa experiência mostra que os métodos de diagnóstico são muito mais prováveis de serem usados quando são convenientes. Por exemplo, gráficos de variáveis adicionadas são construídos regredindo um regressor específico e a resposta em todos os outros regressores, computando os resíduos dessas regressões auxiliares e plotando um conjunto de resíduos contra o outro. Isso não é difícil de fazer em R, embora as etapas sejam um pouco mais complicadas quando há fatores, interações ou termos polinomiais ou de spline de regressão no modelo. A função avPlots no pacote car constrói todos os gráficos de variáveis adicionadas para um modelo linear ou GLM e adiciona melhorias, como uma linha de mínimos quadrados e identificação de ponto.
Os valores \(h_i\), da matriz chapéu, para um modelo linear generalizado podem ser obtidos diretamente da iteração final do procedimento de mínimos quadrados ponderados interados para ajustar o modelo, e têm a interpretação usual - exceto que, ao contrário de um modelo linear, os valores \(h_i\) em um modelos linear generalizado dependem da variável de resposta \(Y\), bem como na configuração dos \(X\).
A matriz chapéu \(H\) é \begin{equation*} H \, = \, W^{1/2}X(X^\top WX)^{-1}X^\top W^{1/2}, \end{equation*} onde \(W\) é a matriz de peso da iteração final do procedimento de estimação.
As observações que estão relativamente longe do centro do espaço do regressor, levando em consideração o padrão correlacional entre os regressores, têm uma influência potencialmente maior nos coeficientes de regressão de mínimos quadrados; tais pontos são considerados como tendo alta alavancagem. A medida mais comum de alavancagem é o \(h_i\) ou hat-values.
O nome hat-values vem da relação entre o vetor de respostas observado e os valores ajustados. O vetor de valores ajustados é dado por \(\widehat{y} = X\widehat{\beta} = Hy\), onde \(H\), definida acima e chamada de matriz hat, projeta \(y\), os valores observados da variável resposta \(Y\), no subespaço estendido pelas colunas da matriz do modelo \(X\). Como \(H = H^\top H\), os valores hat \(h_i\) são simplesmente as entradas diagonais da matriz chapéu.
Os \(h_i\) são limitados entre 0 e 1; em modelos com um intercepto, eles são limitados entre 1/n e 1 e sua soma \(\sum_i h_i\) é sempre igual ao número de coeficientes no modelo, incluindo o intercepto.
Situações nas quais há alguns \(h_i\) muito grandes podem ser problemáticas: em particular, a normalidade de grandes amostras de algumas combinações lineares dos regressores tende a falhar e as observações de alta alavancagem podem exercer influência indevida sobre os resultados (veja abaixo).
A função hatvalues funciona para modelos lineares e GLMs. Uma maneira de examinar os valores de \(h_i\) e outras estatísticas de diagnóstico de observação individual é construir gráficos de índice, representando graficamente as estatísticas em comparação com os índices de observação correspondentes.
Por exemplo, o comando a seguir usa a função influenciaIndexPlot, no pacote car, para produzir a figura acima, que inclui gráficos de índice de resíduos estudentizados, os p-valores Bonferroni correspondentes para o teste de outlier, os hat-values e as distâncias de Cook (discutida em Medidas de influência) para a regressão de intertravamentos:
Observe que as duas primeiras observações se destacam na distância de Cook e no hat-value, o assets, ou seja, os ativos em milhões de dólares na primeira observação é o máximo do conjunto de dados e o interlocks, ou seja, o número de cargos de diretor e executivo interligados compartilhados com outras empresas importantes assim como o assets da segunda observação correspondem a uns dos maiores valores em ambas variáveis.
Os resíduos no podem ser encontrados utilizando a função genérica R residuals e podem calcular-se vários tipos de resíduos. O padrão para um modelo linear é retornar os resíduos ordinários, mesmo se houver pesos. Definir o argumento type = "pearson", retorna os resíduos de Pearson, que produzem resíduos corretamente ponderados se houverem pesos e resíduos ordinários se não houverem pesos. Resíduos de Pearson são o padrão quando os resíduos são usados com um GLM. As funções rstandard e rstudent retornam os resíduos padronizados e estudentizados, respectivamente. A função hatvalues retorna os hat-values.
Gráficos de resíduos em relação aos valores ajustados e em relação a cada um dos preditores, por sua vez, são os gráficos de diagnóstico mais básicos. Se um modelo linear for especificado corretamente, os resíduos de Pearson são independentes dos valores ajustados e dos preditores, e esses gráficos devem ser gráficos nulos, sem características sistemáticas - no sentido de que a distribuição condicional dos resíduos, no eixo vertical do gráfico, não deve ser alterada com os valores ajustados ou com um preditor, no eixo horizontal. A presença de características sistemáticas geralmente implica uma falha de uma ou mais suposições do modelo. De interesse nesses gráficos são as tendências não lineares, as tendências de variação no gráfico e os pontos isolados.
A plotagem de resíduos em relação aos valores ajustados e preditores é útil para revelar problemas, mas menos útil para determinar a natureza exata do problema. Consequentemente, vamos empregar outros gráficos de diagnóstico para sugerir melhorias para um modelo.
Uma variação no gráfico de resíduos básico é o gráfico do modelo marginal, proposto por Cook and Weisberg (1997):
Os gráficos da resposta versus preditores individuais exibem a distribuição condicional da resposta dado cada preditor, ignorando os outros preditores; estes são gráficos marginais no sentido de que mostram a relação marginal entre a resposta e cada preditor contínuo. O gráfico em relação aos valores ajustados é um pouco diferente, pois exibe a distribuição condicional da resposta de acordo com o ajuste do modelo.
Podemos estimar uma função de regressão para cada um dos gráficos marginais ajustando uma suavização aos pontos do gráfico. A função marginalModelPlots usa uma suavização inferior, conforme mostrado pela linha sólida no gráfico.Agora imagine um segundo gráfico que substitui o eixo vertical com os valores ajustados do modelo. Se o modelo for apropriado para os dados, então, sob condições bastante suaves, o ajuste suave para este segundo gráfico também deve estimar a esperança condicional da resposta dado o preditor no eixo horizontal. A segunda suavização também é desenhada no gráfico do modelo marginal, como uma linha tracejada. Se o modelo se ajusta bem aos dados, então as duas suavizações devem corresponder em cada um dos gráficos do modelo marginal; se algum par de alisamentos não corresponder, então temos evidências de que o modelo não se ajusta bem aos dados.
Uma característica interessante dos gráficos do modelo marginal é que, embora o modelo que ajustamos aos dados especifique relações parciais lineares entre interlocks e assets, ele é capaz de reproduzir relações marginais não lineares para esse preditor. Na verdade, o modelo, conforme representado pelas linhas tracejadas, faz um trabalho bastante bom em combinar as relações marginais representadas pelas linhas sólidas, embora as falhas sistemáticas descobertas nos gráficos de resíduos sejam discerníveis aqui.
Uma aproximação da medida de influência da distância de Cook é \begin{equation*} D_i \, = \, \dfrac{R_{P_i}^2}{\widetilde{\phi}(k+1)}\times \dfrac{h_i}{1-h_i}\cdot \end{equation*} Esta é essencialmente a definição de Williams, exceto que dividimos pela dispersão estimada \(\widetilde{\phi}\) para escalar \(D_i\) como uma estatística \(F\) em vez de uma estatística qui-quadrado.
A função InfluencePlot no pacote car fornece uma alternativa para o gráficos de índice de estatísticas de diagnóstico:
Este comando produz um gráfico de bolhas, mostrado na figura acima, que combina a exibição de resíduos estudentizados, valores de \(h_i\) e distâncias de Cook, com as áreas dos círculos proporcionais à distância de Cook.
Gráficos de variáveis adicionadas são um diagnóstico útil para encontrar pontos potencialmente influentes no conjunto, que corresponderão a conjuntos de pontos que estão fora de linha com o resto dos dados e estão no extremo esquerdo ou direito do eixo horizontal.
 |
 |
Wang (1985) sugere uma extensão de gráficos de variáveis adicionadas para GLMs que funciona da seguinte maneira: suponha que o regressor focal seja \(X_j\). Reajuste o modelo com \(X_j\) removido, extraindo os resíduos de trabalho deste ajuste. Em seguida, regredir \(X_j\) nos outros \(X\)s por mínimos quadrados ponderados, usando os pesos da última etapa de mínimos quadrados ponderados iterados, obtendo os resíduos. Finalmente, plote os resíduos de trabalho da primeira regressão contra os resíduos de \(X_j\) da segunda regressão.
Em vez de resumir a influência observando todos os coeficientes simultaneamente, poderíamos criar \(k + 1\) medidas de influência observando as diferenças individuais: \begin{equation} DFBETA_{ij} \, = \, \widehat{\beta}_{(-i)j}-\widehat{\beta}_j, \qquad \mbox{para} \quad j=0,\cdots,k, \end{equation} onde \(\widehat{\beta}_j\) é o coeficiente calculado usando todos os dados e \(\widehat{\beta}_{(-i)j}\) é o mesmo coeficiente calculado com a observação \(i\) omitida. Tal como acontece com \(D_i\), o cálculo de \(DFBETA_{ij}\) pode ser realizado de forma eficiente, sem a necessidade de reajustar o modelo. Os \(DFBETA_{ij}\) são expressos na métrica, unidades de medida, do coeficiente \(\widehat{\beta}_j\).
Uma versão padronizada, \(DFBETAS_{ij}\), divide \(DFBETA_{ij}\) por uma estimativa do erro padrão de \(\widehat{\beta}_j\) calculado com a observação \(i\) removida.
A função dfbeta em R leva um modelo linear ou objeto GLM como seu argumento e retorna todos os valores de \(DFBETA_{ij}\); da mesma forma, dfbetas calcula o \(DFBETAS_{ij}\), como no exemplo a seguir para a regressão de intertravamentos:
Poderíamos examinar cada coluna da matriz dfbetas separadamente, por exemplo, por meio de um gráfico de índices:
A falta de relação entre os valores \(DFBETAS_{ij}\) para os dois regressores reflete a falta de correlção dos próprios regressores.
Os gráficos de componente mais residual e CERES também se estendem diretamente aos GLMs. A suavização não paramétrica dos diagramas de dispersão resultantes pode ser importante para a interpretação, especialmente em modelos para variáveis de resposta binárias, onde a discrição da resposta torna os diagramas difíceis de examinar. Efeitos semelhantes, se normalmente menos extremos, podem ocorrer para dados binomiais e de contagem.
Componente mais residual e CERES traçam um modelo linearizado da última etapa do ajuste IWLS. Por exemplo, o resíduo parcial para \(X_j\) adiciona o resíduo de trabalho \(\widehat{\beta}_j X_{ij}\) o componente mais residual o gráfico então representa o resíduo parcial em relação a \(X_j\). Ao suavizar um gráfico de componente mais residual para um GLM não gaussiano, geralmente é preferível usar um suavizador não robusto.
Um gráfico de componente mais residual para ativos na regressão quase Poisson para os dados da diretoria interligados é mostrado na figura acima. Os ativos são tão inclinados positivamente que o gráfico é diferente de examinar, mas é, no entanto, aparente que a relação parcial entre o número de intertravamentos e ativos é não linear, com uma inclinação muito mais acentuada à esquerda do que à direita. Como a protuberância aponta para a esquerda, podemos tentar endireitar essa relação transformando os ativos na escada de potências e de raízes. Tentativas e erro sugerem a transformação logaritmo de ativos, após a qual um gráfico de componente mais residual para o modelo modificado (Figura abaixo) não é digno de nota.
Finalmente, é importante notar a relação entre os problemas de influência e não linearidade neste exemplo: A observação 1 foi influente na regressão original porque é muito grande e os ativos deram-lhe alta alavancagem e, como a não linearidade não foi modelada, colocou a observação abaixo do ajuste erroneamente linear para os ativos, puxando a superfície de regressão em sua direção. Os ativos com a transformação logaritmo corrigem esses dois problemas.
Vamos utilizar duas formas diferentes de implementação do modelo de regressão Beta.
O pacote R betareg é específico para este modelo de regressão enquanto o segundo, pacote gamlss, é muito mais geral no sentido de conter este modelo como uma situação a mais. Maiores informações acerca dos muitos modelos de regressão implementados neste pacote podem ser consultadas aqui.
Utilizando a pacote gamlss.
Estimação numérica
Vamos entender estes resultados. Obtemos como estimativa do intercepto dos coeficientes do modelo de regressão -0.8064; lembrando que a função de ligação é logística significa que, para obtermos \(\widehat{\mu}\), a estimativa de \(\mu\), a estimativa da esperança da distribuição é
Significa que a estimativa obtiva com 200 observações independentes é bem próxima à teórica 0.3. Uma forma automática de obtermos estas estimativas é mostrada a seguir:
Nos dedicamos agora a entender diversas forma da densidade Beta.
Ainda temos a chamada densidade Beta inflacionada de zeros.
Vejamos o seguinte exemplo. Considere a proporção de óleo cru convertido em gasolina depois do proceso de destilação fracionada. Os resultados consistem em 32 observações e 5 variáveis:
A resposta neste estudo é a variável yield e fazemos com o comando abaixo o correlograma entre todas as variáveis disponíveis. Percebemos que conforme a temperatura à qual toda a gasolina é vaporizada (temp) aumenta, aumenta a proporção de óleo cru convertido em gasolina depois da destilação fracionada (yield). A influência das outras variáveis não está clara.
Modelos de regressão Beta. Primeiro vamos escolher dentre as diversas funções de ligação. A opção control = gamlss.control(trace = FALSE) permite ocultar as diferentes iterações do algoritmo de maximização.
Observamos que a funções de ligação escolhida para o modelo de regressão é probit.
Permanecem as variáveis temp10 e temp.
Mostremos o resultado do ajuste através das estimativas do yield segunda cada uma das variáveis explicativas resultantes no modelo.
Informações de diversos países aparecem a continuação. O objetivo é tentar explicar o Produto Interno Bruto (GNP) per-capita segundo informações populacionais e de escolaridade.
Dispomos das seguintes variáveis:
Mostramos a continuação como realizar a leitura dos dados. Os dados podem ser observados completamente utilizando comando showData, o qual habilita uma janela especifica com a informação contida no argumento, como mostrado abaixo. De outra forma, o comando head mostra somente as primeiras linhas do argumento.
No estudo descritivo abaixo observamos que a relação do número de estudantes no ensino superior por cada 100,000 habitantes com a resposta é relativamente linear positiva, enquanto que as variáveis porcentagem da população letrada, habitantes por médico e mortalidade infantil mostram relação não linear inversa. As outras variáveis não ficam claras a forma de explicarem a resposta, caso isso aconteça.
Escolhemos primeiro a função de ligação.
Selecionamos então a ligação inversa. Verifiquemos se é possível identificar alguma transformação adequada para as variáveis esplicativas.
Prefere-se manter todas as variáveis na forma original e utilzar o Critério AIC para escolher às influintes.
Observamos na figura a seguir os seis diferentes gráficos de resíduos programados. Destaca-se, fundamentalmente, a observação 13.
A observção 13 corresponde aos Estados Unidos. Retiramos a observação identificada e ajustamos o modelo à nova base de dados.
O modelo encontrado é \begin{equation} \widehat{\mu} \, = \, \dfrac{1}{\mbox{4.714e-03}+\mbox{1.859e-05}\times\mbox{Infd}-\mbox{4.257e-05}\times\mbox{Lit}}, \end{equation} e a resposta estimada ou a esperança estimada do Produto Interno Bruto (GNP) segundo as variáveis explicativas mostra-se a continuação.
Este conjunto de dados foi coletado pela Duke University Cardiovascular Disease Databank em 2002 e consiste em 3504 pacientes e seis variáveis. Os pacientes foram encaminhados para a Duke University Medical Center com dor no peito. O arquivo de dados acath com as informações deste estudio está disponível no pacote Hmisc.
Algumas análises interessantes incluem prever a probabilidade de doença coronária (>= 75% de estreitamento do diametro em pelo menos uma arteria coronária importante) e prever a probabilidade de doença coronária grave, uma vez que é uma doença significativa. A primeira análise usaria sigdz como variável de resposta e a segunda usaria tvdlm no subconjunto de pacientes com sigdz = 1. Doença coronária grave é definida como três vasos doentes ou doença principal na esquerda e é denotada por tvdlm = 1.
Observação: o sexo do paciente é codificado como:
Prevendo a probabilidade de doença coronária. A variável resposta é sigdz e é construída como:
A questão agora é saber se as outras variáveis disponíveis explicam a resposta:
 |
 |
Percebemos que a doença coronária se manifesta mais no sexo masculino e está clara a influência da idade mas não do colesterol e nem da duração dos sintomas.
Percebemos que a doença coronária se manifesta principalmente em pessoas maiores de 35 anos e os níveis de colesterol apresentam-se legeiramente superiores em àqueles com resposta positiva à doença.
O gráfico acima mostra navamente que os homens tiveram, fundamentalmente, resposta positiva à doença coronária e que a não está clara a influência da duração dos sintomas da doença arterial coronariana, variável Cad.dur.
Lembremos que a resposta é dicotômica, logo faz sentido pensarmos na distribuição Binomial e acima consideramos o modelo completo. Escolhemos, a seguir, a função de ligação.
Escolhendo a função de ligação mais adequada segundo o menos AIC.
De maneira automática descartamos somente a variável Cad.dur, a duração dos sintomas da doença arterial coronariana.
A proporção do desvio nulo contabilizado pelo modelo é calculado abaixo e, como pode ser observado, é muito baixo.
Utilizamos o pacote
> library(rms)
> R2=lrm(formula(modelo04), x=T, y=T)$stats[10]
> R2
R2
0.2981746
o qual é também um valor baixissimo!
Fica claro no gráfico dos resíduos acima que os resíduos não seguem a distribuição normal padrão. Uma primeira tentativa de melhorar o ajuste &eacue; investirmos mais nas variáveis explicativas e, nesse sentido, fizemos diversas tentativas de transformações de cada varível e interações de duplas delas. Àquela que melhor apresentou-se como explicativa da resposta foi a transformação \(\log\big( Colesterol\times Idade\big)\), vejamos isto aseguir:
Percebemos que conforme os valores da vari´vel Transformada aumentam a proporção de Resposta afirmativa à doença coronária aumenta. Propomos agora o seguinte modelo:
Um pouco melhor, pouco mesmo, mas conseguimos perceber mais simetria e comportamento aleatório no gráfico a esquerda superior, assim como também o gráfico Q-Q Plot mais acertivo quanto à normalidade.
O teste Hosmer-Lemeshow de bondade de ajuste aceita o ajuste do modelo, porém no gráfico de valores estimados vs observados percebemos que há muito erro nas estimativas da resposta. Seria esperado nesse gráfico que para pequenos valores preditos obtenhamos muitos zeros como observados da resposta e para valores altos dos preditos teriamos muitos valores observados 1 da resposta ou vice-versa. Vamos tentar melhorar as estimativas no nosso modelo a seguir.
Para quais valores de \(\widehat{P(Y=1)}\) assumir que o valor estimado da resposta é \(Y=1\) e para quais \(Y=0\)?
As tabelas acima mostram como mudam as estimativas segundo a forma como decidimos escoler quando o valor estimado é 0 ou 1. Assim, por exemplo, a primeira tabela mostra que utilizamos no modelo um total de 768+1490=2258 observções; delas 768 correspondem a resposta codificada como 0 e 1490 correspondem a resposta codificada como 1.
Caso decidirmos que, se \(\widehat{P(Y=1)}>0.5\) então o valor \(\widehat{Y}=1\), padrão no R, obtemos 166 observações estimadas erradamente como sendo \(\widehat{Y}=0\) quando foram observadas sendo \(Y=0\) e 1324 foram estimadas corretamente. Observe que se escolhemos \(\widehat{P(Y=1)}>0.4\) como critério para decidir que \(\widehat{Y}=1\), obtemos 110 observações estimadas erradamente como sendo \(\widehat{Y}=0\) quando foram observadas sendo \(Y=0\) e 1380 foram estimadas corretamente. Isto é conhecido como sensitividade do modelo ou o percentual de verdadeiros positivos; que nesta última situação foi de 1380/1490=0.9261745 (92%).
Um outro conceito importante é a especificidade ou percentual de verdadeiros negativos. Nos dados utilizados para modelar foram observados 768 verdadeiros negativos, enquanto se \(\widehat{P(Y=1)}>0.5\) consideramos \(\widehat{Y}=1\) estimamos 369 verdadeiros negativos; caso \(\widehat{P(Y=1)}>0.4\) como critério para decidir que \(\widehat{Y}=1\), estimamos 293 verdadeiros negativos. Pode-se observar que caso a sensitividade seja muito acurada a especificidade será comprometida.
Nesse sentido a Curva ROC permite-nos calibrar entre sensitividade a especificidade, dessa maneira podemos escolher qual será o critério para escolher se \(\widehat{P(Y=1)}>\delta\) então \(\widehat{Y}=1\).
Observemos que a sensitividade (verdadeiros positivos) obtida com esta escolha é de 70% e a especificidade (verdadeiros negativos) é de 74%. Maiores informações acerca da Curva ROC podem ser encontradas aqui.
Prevendo a probabilidade de doença coronária grave. Resposta: tvdlm caso sigdz=1.
Percebemos que a doença coronária grave se manifesta na mesma proporção nos sexos e não está clara a influencia da idade, colesterol e da duração dos sintomas.
Escolhendo a função de ligação:
Selecionamos a função de ligação cauchit (menor valor de AIC):
Este modelo acerta 80% aproximadamente a probabilidade do indivíduo apresentar doença coronária grave, porém estima muitos falsos positivos.
Dados da alteração do peso em pacientes jovens do sexo feminino com anorexia. Contên 72 linhas de observações e três variás;veis:
Percebemos que existem relação entre o peso ante do estudo e depois em duas situações: tratamento CBT e FT.
Escolhendo a função de ligação masi adequada:
Assim, selecionamos a função de ligação inversa (menor valor de AIC).
De maneira automática não descartamos nenhuma variável explicativa.
A figura acima resume o ajuste de mosso modelo, apresentamos as curvas dos pesos estimados depois do estudo para cada um dos tratamentos.
No arquivo olhos.txt são apresentados dados referentes a 78 famílias com pelo menos seis filhos cada uma. Este archivo contêm as seguintes informações:
Seja \(Y_i\) o número de filhos com olhos claros pertencentes à \(i\)-ésima família. Sugere-se então considerar um modelo logístico linear. Existem pontos aberrantes? Há indícios de superdispersão? Verifique se o fator ''cor dos olhos dos avôs'' deve permanecer no modelo.
Todos as variáveis neste exemplo são categóricas, significa que a realização do estudo descritivo é mais complicado. Utilizamos a seguinte código.
 |
 |
No gráfico a esquerda observamos que os códigos 1, 4 e 5 da classificação dos olhos dos pais segundo a cor destacam-se por serem mais comuns. Aparentemente, a cor dos olhos dos avôs influencia pouco na ocorrência de olhos claros nos netos. No gráfico à direita percebemos que as categorias 4 e 9 da classificação dos olhos dos avôs segundo a cor influenciam na ocorrência de olhos claros nos netos.
Selecionamos o modelo com a função de ligação canônica (menor valor de AIC) e escolhemos quais variáveis devem permanecer no modelo.
De maneira automáica descartamos a cor dos olhos dos avôs.
Percebemos na análise de resíduos que existem três pontos aberrantes: 18, 37 e 47, fundamentalmente. Observamos abaixo suas características e percebemos que são situações nas quais os pais tiveram poucos filhos com olhos claros.
Decidimos então eliminar os pontos aberrantes/influentes e continuar com a análise.
Decide-se por manter o fator da classificação dos olhos dos pais modelo. Posteriormente tentaremos ajustar modelos considerando a possibilidade de superdispersão nesses dados.
O cálculo descrito aqui é a base do pacote qvcalc de Firth (descrito em Firth, 2003) para o ambiente de programação estatística R.