• Introdução
• Definição e construção
• Resultados teóricos
• Estimação
• Inferência
• Exemplos
• Considerações finais

Neste seminário vamos apresentar os Modelos Lineares Generalizados (MLG), ou no inglês Generalized Linear Models, introduzidos em 1972 por John Nelder and Robert Wedderburn.

Este modelos recebem este nome pois são uma generalização dos modelos lineares tradicionais de regressão. Através de uma função de ligação, pode-se permitir que a variável resposta tenha uma distribuição (relativamente) arbitrária. Tais modelos tem sido utilizados intensivamente e tem se demonstrado muito úteis.

O presente trabalho está dividido em 6 partes. Primeiramente veremos os conceitos básicos com a definição e construção do modelo, logo depois resultados teóricos, seguiremos com a estimação, depois com detalhes de inferência, uma seção com alguns exemplos e, por fim, as considerações finais.

O modelo de regressão simples tem como objetivo predizer a esperança de uma variável aleatória através de uma combinação linear de suas covariáveis. Na notação matricial comum:

\[ {\bf y} = {\bf X}\beta + {\bf \epsilon}\]

Para isso, algumas suposições fortes são assumidas:

• o vetor de erros tem distribuição normal de média 0 e variância \(\sigma^2\)

• Efeitos fixos e lineares.

• Homoscedasticidade: variância constante.

Porém, em muitas situações, tais suposições não são razoáveis. Por exemplo, a primeira suposição, dos erros, \(\epsilon \sim N(0, \sigma^2)\), implica que \(Y \sim N({\bf X}\beta, \sigma^2)\). Entretanto, muitas vezes a variável resposta não tem uma distribuição normal.

Considere os seguintes estudos:

• Estudos sobre o \(defaults\) de empréstimos;
• número de bactérias sobreviventes em amostras de um produto alimentício exposto.
• nível de doêncas;

Em geral, nas análises de regressão, procurava-se algum tipo de transformação que levasse à normalidade, tais como a transformação de Box-Cox (1964) dada abaixo:

\[z = \begin{cases} \frac{y^\lambda -1}{\lambda}, & \text{se } \lambda \neq 0, \\ log(y), & \text{se } \lambda = 0. \end{cases} \] onde \(y>0\) e \(\lambda\) é uma constante desconhecida.

Sendo assim, um dos métodos propostos são os MGLs. Os MLGs foram criados com o objetivo de reunir numa mesma família vários modelos estatísticos que eram tratados separadamente.

A construção do modelo é dada pela suposição:

\[g(\mathbb{E}(Y|{\bf X})) = {\bf X}\beta \] Ou seja, o modelo assume que a média de uma função da esperança é uma combinação linear dos parâmetros. Outra definição equivalente é:

\[\mathbb{E}(Y|{\bf X}) = g^{-1}({\bf X}\beta) \]

O modelo consiste dos seguintes elementos:

• Variável resposta: \({\bf Y}\).
• Preditor linear, composto das preditoras e um vetor de parâmetros: \({\bf X}\beta\). Muitas vezes referido como \({\bf \eta}\).
• Função de ligação: \(g(\cdot)\).

Em geral, escolhe-se uma distribuição para \({\bf Y}\). Preferencialmente da família exponencial.

A família exponencial é uma família de funções de densidade de probabilidade \(f\) que satisfazem:

\[f(y, \theta, \phi ) = exp\Big\{ \frac{y \theta - b(\theta)}{a(\phi)} + c(y, \phi) \Big\}\] Costuma-se chamar o parâmetro \(\theta\) de \({\bf natural}\) e \(\phi\) de parâmetro de dispersão. Tal nomenclatura se deve ao fato de que densidades da família exponencial escritas nessa formulação apresentam vantagens na estimação.

Algumas das distribuições que fazem parte da família exponencial:

• Gaussiana: \(\mathbb{R}\)
• Multinomial: Categórico
• Bernoulli: Binário \(\{0, 1\}\)
• Binomial: Contagens de sucesso/fracasso
• Gama: \(\mathbb{R}\)
• Poisson: \(\mathbb{N}\)
• Laplace: \(\mathbb{R}\)
• Exponencial: \(\mathbb{R}^+\)
• Beta: \((0, 1)\)
• Dirichlet: ∆ (Simplex)
• Weibull: \(\mathbb{R}\)
• Weishart: matrizes simétricas positivas definidas

A função geradora de momentos de \(Y_i\) é dada por:

\[M_{Y_i}(t) = \mathbb{E}(e^{tY_i}) = \int_A e^{tY_i}f(y_i, \theta,\phi)dy_i \] Integrando e utilizando o fato que a primeira derivada avaliada em \(t = 0\) é a esperança, temos:

\[ \mathbb{E}(Y_i) = \frac{d}{dt}M_{Y_i}(t) \Bigr|_{\substack{t=0}}= exp\Big( \frac{b(ta(\phi) +\theta_i) - b(\theta_i)}{a(\phi)} \Big)\frac{b'(ta(\phi) +\theta_i)a(\phi)}{a(\phi)}\Bigr|_{\substack{t=0}}\] Ou seja:

\[\mathbb{E}(Y_i) = \mu_i = b'(\theta_i)\]

Sob as condições de regularidade usuais é possível provar que:

\[Var\Big(\frac{\partial l_i}{\partial\theta_i} \Big) = - \mathbb{E}\Big(\frac{\partial^2 l_i}{\partial\theta_i^2} \Big) \]

Agora, por outro lado, percebendo que:

\[ \frac{\partial l_i}{\partial\theta_i} = \frac{y_i - b'(\theta_i)}{a(\phi)} \ \ \ \ \ ; \ \ \ \ \ \frac{\partial^2 l_i}{\partial\theta_i^2} = \frac{ - b''(\theta_i)}{a(\phi)}\]

Pode-se provar um resultado fundamental:

\[ Var(Y_i) = a(\phi)b''(\theta_i)\] Observação: a quantidade \(b''(\theta)\) é conhecida como função de variância. Ela auxilia tanto no cálculo da variância quanto da identificação da distribuição.

As funções de ligação têm como objetivo relacionar o preditor linear com a esperança \(\mu\). Ou seja:

\[ \sum_{i = 1}^p x_{ij} \beta_j = \eta_i = g(\mu_i)\] A função de ligação que transforma \(\mu_i\) para \(\theta_i\) é denominada função de ligação canônica.

Considere \(Y \sim P(\mu)\). Então:

\[\begin{align*} f(y;\mu) = \frac{e^{-\mu}\mu^y}{y!} &= e^{log\big(\frac{e^{-\mu}\mu^y}{y!}\big) }\\ &= e^{-\mu + ylog(\mu) - log(y!) } \end{align*}\]

Porém, para adequar a definição utilizada, precisamos transformar para o parâmetro natural. A substituição pode ser feita por:

\[ \mu = exp(\theta) \iff log(\mu) = \theta\] Garantindo que:

\[\begin{align*} f(y;\mu) &= exp \Big\{-exp\{\theta\} + y\theta - log(y!) \Big\} \end{align*}\]

Ou seja, a Poisson(\(\mu\)) pertence à família exponencial com parâmetros:

• \(b(\theta) = exp\{\theta\}\)
• \(a(\phi) = 1\)
• \(c(y, \phi) = - log(y!)\)

Logo, a função de ligação canônica é \(log(\cdot)\). Pois:

\[log(\mu) = \theta \]

A escolha do uso das ligações canônicas se deve ao fato que estas funções mantém a concavidade da log-verossimilhança. Ou seja, garantimos a unicidade da estimativa por máxima verossimilhança e temos resultados assintóticos fortes.

Outro resultado importante se deve a Nelder e Wedderburn (1972). Eles garantem que a estimação por \(Newton-Raphson\) e pelo \(Escore\) \(de\) \(Fisher\) são idênticas para as funções canônicas.

Existem vários métodos de estimação existentes. Os mais usuais são: máxima verossimilhança, mínima deviância e mínimos quadrados. Em geral, o primeiro é o mais utilizado. O segundo método é facilmente empregado no modelo normal, porém sabemos (pelo Teorema de Gauss-Markov) da equivalência entre os dois métodos.

É importante comentar que na prática a minimização da deviância funciona quase como um critério de ajuste.

Estes métodos necessitam algoritmos numéricos, por exemplo: Newton-Raphson, Escore de Fisher, BFGS, Nelder-Mead, etc.

Na função \(glm\) implementada no \(r\) usa-se o Escore de Fisher.

Considerando somente a família exponencial, a log-verossimilhança do modelo é dada por:

\[\mathcal{l} ({\bf \beta}, \phi) = \sum_{i=1}^n \Bigg(\frac{y_i\theta_i - b(\theta_i)}{a(\phi} + c(y_i, \phi)\Bigg) \] Em que \(\theta_i = q(\mu_i)\).

Considerando que a relação da função de ligação: \(\eta_i = g(\mu_i)\), temos:

\[ \frac{\partial l_i}{\partial \beta_j} = \frac{\partial l_i}{\partial \theta_i} \frac{\partial \theta_i}{\partial \mu_i} \frac{\partial \mu_i}{\partial \eta_i} \frac{\partial \eta_i}{\partial \beta_j} \] O que implica que a verossimilhança é:

\[ \frac{\partial l({\bf \beta}, \phi)}{\partial \beta_j} = \sum_{i=1}^n \frac{((y_i-\mu_i))x_{ij}}{Var(Y_i)}\frac{\partial \mu_i}{\eta_i}\]

Maximizando as equações, temos o desejado.

É possível provar que a distribuição assintótica de \({\bf \beta}\) sob as condições de regularidade é \(\sim N_p({\bf \beta}, ({\bf X'WX})^{-1})\). Onde \(p\) é o número de covariáveis e \({\bf W}\) é a diagonal da matriz da informação esperada de Fisher.

Testar a qualidade geral do ajuste pode ser feito utilizando a função desvio. Comparando a diferença entre a verossimilhança de um modelo saturado e um modelo não saturado, pode-se ter uma ideia do ajuste. Isto é, desvio pequenos indicam que o modelo proposto se ajusta bem aos dados.

Para isso, considere a log-verossimilhança do modelo:

\[ l({\bf \mu, y}) = \sum_{i = 1}^n \mathcal{l} ({\bf \mu_i}; y_i) \]

E a log-verossimilhança saturada. Ou seja, onde a média é estimada pela própria observação (também podemos pensar como \(p = n\)):

\[ l({\bf y; y}) = \sum_{i = 1}^n \mathcal{l} (y_i; y_i) \]

A estatística então é dada da forma:

\[ D({\bf y; \hat{\mu}}) = 2[l({\bf y; y}) - l({\bf \hat{\mu}; y})]\] Onde \(D({\bf y; \hat{\mu}})\) tem \(n-p\) graus de liberdade.

Embora o desvio muitas vezes não tenha uma distribuição conhecida (mesmo assintoticamente), na prática costuma-se a aceitar o teste:

\[ \frac{D({\bf y; \hat{\mu}})}{\phi} \leq \chi^2_{n-p, \alpha}\]

Para considerar que existem evidências, a um nível de significância \(\alpha\), do bom ajusto do modelo proposto.

Em especial este teste deriva do resultado de Jφrgensen (1987).

Outra estatística utilizada é a estatística de Pearson generalizada, da forma:

\[ X^2 = \sum_i \frac{w_i(y_i- \hat{\mu_i})^2}{V(\hat{\mu_i})}\]

Onde \(V(\hat{\mu}_i)\) representa a função de variância.

Em geral, tal estatística tem distribuição qui-quadrado.

Considere o teste:

\[ H_{0} : \mathop{\beta}\limits^{\sim} = 0 \ \ vs\ H_{1} : \mathop{\beta}\limits^{\sim} \neq 0\] onde \(\mathop{\beta}\limits^{\sim}\) representa algum subconjunto dos parâmetros.

Realizar este teste pode ser feito pensando em modelos encaixados. Onde queremos verificar se alguns dos parâmetros não tem efeitos significativos.

A estatística utilizada são os próprios desvios:

\[D({\bf y; \hat{\mu_0}})-D({\bf y; \hat{\mu_1}})\] Assintóticamente e conhecida a dispersão, temos que a estatística tem distribuição \(\chi^2_{p_1 - p_0}\)

Onde \(p_1 - p_0\) é a diferença da quantidade de parâmetros do modelo 1 (completo) para o modelo proposto.

Outros testes estão disponíveis e são muito utilizados em geral:

• Razão de verossimilhanças: considerando \(L_0\) o modelo relativo a hipótese testada:

\[-2l\Big(\frac{L_0}{L_1}\Big)\sim \chi^2_p \]

• Wald:

\[W = (\hat{\beta} -\beta_0)'[\hat{Cov}(\hat{\beta})]^{-1}(\hat{\beta} -\beta_0) \sim \chi^2_p\] No caso do teste para um único efeito, \(\sqrt{W} \sim N(0,1)\)

• Escore: usa a função Escore e a curvatura esperada da função de log-verossimilhança, avaliada no valor da hipótese nula:

\[U'_\beta(\beta_0)(I(\beta_0))^{-1}U_\beta(\beta_0)\sim \chi^2_p \] onde \(U(\cdot)\) representa o vetor \(score\) e \(I(\cdot)\) a matriz de informação esperada. .

Na função \(glm\) o teste utilizado é o de \(Wald\)*.

Na regressão linear consi dera-se normalmente o resíduo:

\[y_i - \hat{y_i}\]

Porém, em GLM tal resíduo não faz muito sentido e não goza de propriedades interessantes. Existem diversos tipos de resíduos propostos para os GLM, entre eles: Pearson, desvio, padronizado, quantílico, aleatorizado etc.

Como na regressão linear, o estudo dos resíduos é uma ferramenta para identificar onde o ajuste de um MLG é razoável. Este estudo também permite a identificação de outliers.

Alguns resíduos para algumas distribuições tem aproximações aproximadamente normais. A bibliografia sobre o assunto é extensa: Williams 1984, McCullagh, 1987, Dunn e Smyth 1996.

No artigo Certainty vs. Severity of Punishment, Jeffrey Grogger examina os efeitos de punições criminais, emprego e renda na redução da criminalidade. Ele utiliza um grande conjunto de dados sobre históricos criminais e de mercado de trabalho de jovens homens presos. O estudo avalia o impacto da certeza e severidade das punições, além de efeitos como incapacitação e capital humano criminoso, e como emprego e renda influenciam a atividade criminosa. Os resultados conciliam achados conflitantes de pesquisas anteriores.

Para modelar o número de prisões, Grogger usa uma função não linear das variáveis explicativas, como raça, etnia e ano de nascimento. A teoria econômica prevê que a probabilidade de condenação e a severidade das penas tenham efeitos negativos, assim como emprego e renda, já que o desemprego aumenta a atração pelo crime.

Como as prisões são uma variável de contagem, ele aplica uma regressão de Poisson, ideal para esse tipo de dado, garantindo que as previsões sejam inteiros não negativos, semelhante ao que os modelos \(probito\) e \(logito\) fazem para dados binários.

Foram coletados as seguintes covariáveis:

• \(narr86_i\): Número de vezes preso em 1986 (variável dependente)
• \(pcnv_i\): Proporção de condenações anteriores
• \(avgsen_i\): Duração média da sentença (meses) -\(tottime_i\): Tempo total na prisão desde os 18 anos (meses) -\(ptime86_i\): Meses na prisão em 1986
• \(qemp86_i\): Número de trimestres empregado em 1986
• \(inc86_i\): Renda legal em 1986 (em centenas de dólares)
• \(black_i\): Dummy para raça negra
• \(hispan_i\): Dummy para hispânico
• \(born60_i\): Dummy para nascidos em 1960 -\(\varepsilon_i\): Erro aleatório no modelo OLS

Logo, os modelos propostos são:

\[ \begin{multline*} narr86_i = \beta_0 + \beta_1 pcnv_i + \beta_2 avgsen_i + \beta_3 tottime_i + \beta_4 ptime86_i + \beta_5 qemp86_i\\ + \beta_6 inc86_i + \beta_7 black_i + \beta_8 hispan_i + \beta_9 born60_i + \varepsilon_i \end{multline*} \] E, com MLG: \[ \begin{equation*} \mathbb{E}(narr86_i \mid X_i) = \exp(\textbf{X'}\beta) \end{equation*} \]

ols_model <- lm(narr86 ~ pcnv + avgsen + tottime + ptime86 +
                  qemp86 + inc86 + black + hispan + born60, data = crime1)
poisson_model <- glm(narr86 ~ pcnv + avgsen + tottime + ptime86 +  qemp86 + inc86 + black + hispan + born60,  family = poisson(link = "log"), data = crime1)
res <- stargazer(ols_model,poisson_model) 

Agora vamos trabalhar com um modelo binomial. Vamos considerar o conjunto de dados \(frogs\).

Nestes dados temos o estudo da distribuição dos sapos da espécie \(Southern\) \(Corroboree\), na área de \(Snowy\) \(Mountains\), na Austrália.

Sapinho

• sapo: indica de há sapo ou não (0-1);
• norte: distância ao norte em relação ao ponto de referência;
• leste: distância ao leste em relação ao ponto de referência;
• alt: altura em metros;
• dis: distância até a população mais próxima;
• lagos: #n de potenciais criadores;
• reg: n de potenciais criadores em 2km;
• chuva: chuva média na primavera;
• c_min: chuva mínima na primavera;
• c_max: chuva máxima na primavera;

m1 <- glm(formula = sapo ~ alt + dis + lagos + reg + chuva,family = binomial(), data = frogs)
summary(m1)

## 
## Call:
## glm(formula = sapo ~ alt + dis + lagos + reg + chuva, family = binomial(), 
##     data = frogs)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.9560229  2.3568358   1.679 0.093243 .  
## alt         -0.0122487  0.0025902  -4.729 2.26e-06 ***
## dis         -0.0005789  0.0001932  -2.996 0.002731 ** 
## lagos        0.0248106  0.0083524   2.970 0.002973 ** 
## reg          0.0190144  0.1019500   0.187 0.852047    
## chuva        0.0977021  0.0275511   3.546 0.000391 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 279.99  on 211  degrees of freedom
## Residual deviance: 210.56  on 206  degrees of freedom
## AIC: 222.56
## 
## Number of Fisher Scoring iterations: 6

Então perceba o teste de hipótese (feito individualmente), não rejeita a hipótese de que o efeito da variável “\(reg\)” seja nulo. Logo podemos descartá-la do modelo:

## 
## Call:
## glm(formula = sapo ~ alt + dis + lagos + chuva, family = binomial(), 
##     data = frogs)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.1064847  2.2131372   1.856 0.063524 .  
## alt         -0.0121824  0.0025646  -4.750 2.03e-06 ***
## dis         -0.0005954  0.0001734  -3.434 0.000596 ***
## lagos        0.0249273  0.0083352   2.991 0.002784 ** 
## chuva        0.0964965  0.0267310   3.610 0.000306 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 279.99  on 211  degrees of freedom
## Residual deviance: 210.60  on 207  degrees of freedom
## AIC: 220.6
## 
## Number of Fisher Scoring iterations: 6

Perceba que o \(deviance\) não altera significativamente entre os dois modelos. Ou seja, aparentemente a qualidade dos ajustes dos modelos são iguais. Por parcimônia, o segundo é preferível, porém, para garantir a validade devemos realizar um teste de hipótese.

No \(r\), neste caso, dado a construção do modelo, podemos realizar o teste dos modelos encaixados:

des_dif <- deviance(m2) - deviance(m1)
p_valor <- 1-pchisq(des_dif, df = 1) 
p_valor

## [1] 0.8521861

Sendo assim, podemos concluir que o modelo 2 é de fato preferível.

A interpretação pode ser realizada considerando o modelo:

\[log(\mu) = 4.1-0.012alt -0.0005954dis +0.0249lagos+0.09chuva \]

Uma observação útil: a covariável \(dis\) aparenta ter efeito 0, mas é importante lembrar que seus valores são numericamente altos. Um teste de hipótese confirma a avaliação:

m3<- glm(formula = sapo ~ alt  + lagos + chuva,family = binomial(), data = frogs)
des_dif <- deviance(m3) - deviance(m2)
p_valor <- 1-pchisq(des_dif, df = 1) 
p_valor

## [1] 5.513202e-07

É interessante sempre realizar a seleção das covariáveis utilizando métodos como o \(step-wise\).

Por fim, uma das ferramentas para avaliar a qualidade geral do ajuste pode ser a análise dos resíduos. Neste caso:

Embora tenhamos falado da modelagem. Não falamos de dois tópicos muito importantes: previsão e interpretação. Considerando o último exemplo, estimamos a média, que ainda é um valor real entre 0 e 1. Neste caso, para realizar a previsão, estratégias como a curva \(ROC\) são eficientes e de fáceis implementações no \(r\).

Outro detalhe importante é notar que as análises feitas dependem da função de ligação utilizada. Por exemplo, escolhendo a função de ligação complemento log-log, no exemplo dos sapos, o efeito fixo deixa de ser significativo e o \(AIC\) abaixa, nos fornecendo indícios de um melhor modelo. Porém, é útil perceber das limitações do uso de funções de ligação não canônicas.

A interpretação do modelo, em geral, depende da construção.

No exemplo dos sapos, a interpretação do modelo pode ser feita considerando a razão de chances. Ilustrando: o que acontece quando o número de potenciais criadores aumenta de 3 para 4, mantendo as outras covariáveis iguais?

\[ \begin{align*} \frac{\mu_i}{\mu_j} &= \frac{exp\{4.1-0.012alt -0.0005954dis +0.0249\cdot 4 +0.09chuva\}}{ exp\{4.1-0.012alt -0.0005954dis +0.0249\cdot 3+0.09chuva\}}\\ &= exp\{0.0249\}\\ &= 1.0252 \end{align*}\]

Ou seja, há um aumento de 2,52% nas chances de encontrar sapos.

Um ponto de grande importância que não foi abordado é a estimação do parâmetro de dispersão \(\phi\). Nos exemplos apresentados o parâmetro sempre foi estipulado como 1, porém para algumas distribuições também é necessário realizar sua estimação.

Resultados assintóticos e testes de hipóteses sofrem uma pequena alteração. Até mesmo a estimação de \(\beta\) se altera.

Os modelos lineares generalizados não tem uma construção relativamente complicada e sua implementação, em si, não é muito diferente de outros métodos. Logo, a flexibilização proporcionada por estes modelos e seus resultados geralmente satisfatórios, são fatores que justificam a boa aceitação do modelo.

Por outro lado, mesmo que a família exponencial contenha diversas distribuições, seu uso é limitado. Na prática, os modelos lineares tem problemas com outliers e dados de alta dimensão.

AGRESTI, A. Foundations of Linear and Generalized Linear Models. 2012. Wiley series in probability and statistics.

AVILLA, L. Notas de aula. 2023.

DOBSON, A.J, BARNETT, A.G. An Introduction to Generalized Linear Models. 2008. CRC Press.

GENTLE, J.E., HARDLE, W.K, MORI, Y. Handbook of Computational Statistics. Berlin, Heidelberg: Springer Berlin Heidelberg. Disponível: https://doi.org/10.1007/978-3-642-21551-3_1.

GROGGER, J. CERTAINTY VS. SEVERITY OF PUNISHMENT. Economic Inquiry, 29: 297-309, 1991. Disponível em: https://doi.org/10.1111/j.1465-7295.1991.tb01272.x

NELDER, J.A., WEDDERBURN, R.W.M. (1972). Generalized Linear Models. Journal of the Royal Statistical Society: Series A (General), 135(3), 370–384. doi:10.2307/2344614

PAULA, G.A. Modelos de Regressão com Apoio Computacional. 2013.

Agradecimento especial para a Prof. Dra. Larissa Ávilla pelo excelente material e curso em GLM.

A D.L.B.