class: center, middle, inverse, title-slide .title[ # MODELOS LINEARES GENERALIZADOS STA13829 ] .subtitle[ ## Modelos Lineares Generalizados ] .author[ ### Nátaly A. Jiménez Monroy ] .institute[ ### LECON/DEST - UFES ] --- [//]: <> (https://pkg.garrickadenbuie.com/extra-awesome-xaringan/intro/index.html#1) [//]: <> (https://pkg.garrickadenbuie.com/xaringanthemer/articles/xaringanthemer.html) [//]: <> (https://www.biostatistics.dk/talks/CopenhagenRuseRs-2019/index.html#1) [//]: <> (https://rstudio-education.github.io/sharing-short-notice/#1) [//]: <> (https://www.kirenz.com/slides/xaringan-demo-slides.html#1) [//]: <> (https://github.com/yihui/xaringan/issues/26) [//]: <> (https://github.com/emitanaka/anicon) [//]: <> (https://github.com/mitchelloharawild/icons) [//]: <> (https://slides.yihui.org/2020-genentech-rmarkdown.html#1) [//]: <> (https://github.com/gadenbuie/xaringanExtra) [//]: <> (class: center, middle, animated, slideInRight) class: animated, slideInRight <style> body {text-align: justify} </style> <!-- Justify text. --> # Introdução - I - Nelder e Wedderburn (1972) mostraram que uma série de técnicas estatísticas, comumente estudadas separadamente, podem ser formuladas, de uma maneira unificada, como uma classe de modelos de regressão. A essa teoria unificadora de modelagem estatística, uma extensão dos modelos clássicos de regressão, deram o nome de **<span style="color:orange">modelos lineares generalizados</span>** (MLG). Esses modelos envolvem uma variável resposta univariada, variáveis explanatórias e uma amostra aleatória de `\(n\)` observações independentes, sendo que: -- - a variável resposta, **<span style="color:orange">componente aleatório</span>** do modelo, tem uma distribuição pertencente a família de distribuições apresentada na Equação (3) que engloba as distribuições normal, gama e normal inversa para dados contínuos; binomial para proporções; Poisson para contagens; -- - as variáveis explanatórias entram na forma de uma estrutura linear, constituindo o **<span style="color:orange">componente sistemático</span>** do modelo; -- - a ligação entre os componentes aleatório e sistemático é feita através de uma função adequada como, por exemplo, logarítmica para os modelos log-lineares, chamada **<span style="color:orange">função de ligação</span>**. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Introdução - II O componente sistemático é estabelecido durante o planejamento (fundamental para a obtenção de conclusões confiáveis) do experimento, resultando em modelos de regressão (linear simples, múltipla, não linear, etc). O componente aleatório é estabelecido assim que são definidas as medidas a serem feitas, que podem ser contínuas ou discretas, exigindo o ajuste de distribuições diferentes. A partir de um mesmo experimento podem ser obtidas medidas de diferentes tipos, como por exemplo, dados de altura de plantas, número de lesões por planta e proporção de plantas doentes. -- No modelo de regressão clássico, tem-se $$ \boldsymbol{Y = \mu + \epsilon}, $$ sendo `\(\boldsymbol{Y}\)` o vetor, de dimensão `\(n \times 1\)`, da variável resposta, `\(\boldsymbol{\mu} = E(\boldsymbol{Y}) = \boldsymbol{X \beta}\)`, o componente sistemático, `\(\boldsymbol{X}\)` a matriz do modelo, de dimensões `\(n \times p\)`, `\(\boldsymbol{\beta} = (\beta_1, \; \beta_2, \; \cdots, \; \beta_p)^T\)`, o vetor dos parâmetros, `\(\boldsymbol{\epsilon} = (\epsilon_1, \; \cdots, \; \epsilon_n)^T\)` , o componente aleatório com `\(\epsilon_i \sim N(0, \sigma^2)\)`, `\(i = 1, ..., n\)`. Nesse caso, tem-se que `\(\boldsymbol{Y} \sim N_p(\boldsymbol{\mu}, \sigma^2 \boldsymbol{I})\)` e o vetor de médias `\(\boldsymbol{\mu}\)` da distribuição normal, que define o componente aleatório, é igual ao preditor linear que representa o componente sistemático. Essa é a forma mais simples de ligação entre esses dois componentes, sendo denominada **<span style="color:orange">função de ligação identidade</span>**. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Estruturação matemática de um MLG - I Os modelos lineares generalizados podem ser usados quando se tem uma única variável aleatória `\(Y\)` associada a um conjunto de variáveis explanatórias `\(x_1, \cdots, x_p\)`. Para uma amostra de `\(n\)` observações, `\((y_i; \boldsymbol{x}_i)\)` em que `\(\boldsymbol{x}_i = (x{i_1}, \cdots, x_{ip})^T\)` é o vetor coluna de variáveis explicativas, o MLG envolve os três componentes. -- - **Componente aleatório**: representado por um conjunto de variáveis aleatórias independentes `\(Y_1, \cdots, Y_n\)` provenientes de uma mesma distribuição que faz parte da família de distribuições na Equação (3) com médias `\(\mu_1, \cdots, \mu_n\)`, ou seja, $$ E(Y_i) = \mu_i, \quad i = 1, \cdots, n, $$ sendo `\(\phi> 0\)` um parâmetro de dispersão e o parâmetro `\(\theta_i\)` denominado parâmetro canônico. Então, a função de probabilidade de `\(Y_i\)` é dada por: -- `\begin{equation} f(y_i; \theta_i, \phi) = \exp[\phi \{ \theta_i y_i - b(\theta_i)\} + c(y_i, \phi)], \label{eq:fam_exp_uni_33} \end{equation}` em que `\(b(\cdot)\)` e `\(c(\cdot)\)` são funções conhecidas. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Estruturação matemática de um MLG - II Conforme visto, valem as seguintes relações: .center[ `\(E(Y_i) = \mu_i = b'(\theta_i), \;\;\; Var(Y_i) = \phi^{-1}V_i\)`, ] sendo `\(\phi^{-1} > 0\)` `\((\phi > 0)\)` o parâmetro de dispersão (precisão) e `\(V_i = V(\mu_i) = \partial\mu_i/\partial\theta_i\)` é a função de variância que depende unicamente da média `\(\mu_i\)`. -- O parâmetro natural `\(\theta_i\)` pode ser expresso como $$ \theta_i = \int V_i^{-1} \partial \mu_i = q(\mu_i), \;\; \mbox{pois} \;\;\; \frac{\partial \theta_i}{ \partial \mu_i} = \frac{1}{V_i} = q(\mu_i), $$ sendo `\(q(\mu_i)\)` uma função conhecida da média `\(\mu_i\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Estruturação matemática de um MLG - III - **Componente sistemático**: As variáveis explicativas entram na forma de uma soma linear de seus efeitos `$$\eta_i = \sum_{r = 1}^p x_{ir} \beta_r = \boldsymbol{x}_i^T \boldsymbol{\beta}\;\;\; \mbox{ou} \; \boldsymbol{\eta = X \beta},$$` sendo `\(\boldsymbol{X} = (\boldsymbol{x}_1, \cdots, \boldsymbol{x}_n)^T\)` a matriz do modelo, `\(\boldsymbol{\beta} = (\beta_1, \cdots, \beta_p)^T\)` o vetor de parâmetros e `\(\boldsymbol{\eta} = (\eta_1, \eta_2, \cdots, \eta_n)^T\)` o preditor linear. Se uma parâmetro tem valor conhecido, o termo correspondente na estrutura linear é chamado *offset*. -- - **Função de ligação**: uma função que **relaciona** o componente aleatório ao componente sistemático, ou seja, vincula a média ao preditor linear, isto é, `$$\eta_i = g(\mu_i),$$` sendo `\(g(\cdot)\)` uma função monótona e diferenciável. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Observações - Os parâmetros `\(\theta_i\)` da família de distribuições (\ref{eq:fam_exp_uni_33}) não são de interesse direto (pois há um para cada observação). -- - Temos interesse sim em um conjunto menor de parâmetros, `\(\beta_1, \cdots, \beta_p\)` tais que uma combinação linear dos `\(\beta's\)` seja igual a alguma função do valor esperado de `\(Y_i\)`. -- - Pelo menos em teoria, cada observação pode ter uma função de ligação diferente. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Funções de ligação - I Supondo `\(\phi\)` conhecido, o logaritmo da função de verossimilhança de um MLG com respostas independentes pode ser expresso na forma: `$$L(\boldsymbol{\beta}) = \sum_{i = 1}^n \phi \{ y_i \theta_i - b(\theta_i) \} + \sum_{i = 1}^n c(y_i, \phi).$$` -- Um caso particular importante ocorre quando o parâmetro canônico `\((\theta)\)` coincide com o preditor linear, isto é, quando `\(\theta_i = \eta_i = \sum_{j = 1}^p x_{ij} \beta_j\)`. Nesse caso, `\(L(\boldsymbol{\beta})\)` fica dado por: `$$L(\boldsymbol{\beta}) = \sum_{i = 1}^n \phi \left\{ y_i \sum_{j = 1}^p x_{ij} \beta_j - b\left(\sum_{j = 1}^p x_{ij} \beta_j \right) \right\} + \sum_{i = 1}^n c(y_i, \phi).$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Funções de ligação - II Definindo a estatística `\(S_j = \sum_{i = 1}^n Y_i x_{ij}\)`, temos que `\(L(\boldsymbol{\beta})\)` fica re-expresso por: `$$L(\boldsymbol{\beta}) = \phi \sum_{j = 1}^p s_j \beta_j - \phi \sum_{i = 1}^n b\left(\sum_{j = 1}^p x_{ij} \beta_j \right) + \sum_{i = 1}^n c(y_i, \phi).$$` Logo, pelo teorema de fatoração, a estatística `\(\boldsymbol{S} = (S_1, \cdots, S_p)^T\)` é suficiente para o vetor `\(\boldsymbol{\beta} = (\beta_1, \cdots, \beta_p)^T\)` e tem dimensão mínima `\(p\)` e, portanto, ocorre uma redução na dimensão das estatísticas suficientes de `\(n\)` (o número de observações) para `\(p\)` (o número de parâmetros a serem estimados). -- As estatísticas `\(S_1, \cdots, S_p\)` correspondem à maior redução que os dados podem alcançar, sem qualquer perda de informação relevante para se fazer inferência sobre o vetor de parâmetros desconhecidos `\(\boldsymbol{\beta}\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Observações - As ligações que correspondem a tais estatísticas suficientes `\(S_1, \cdots, S_p\)` são chamadas de ligações canônicas e desempenham um papel importante na teoria dos MLG's. -- - As ligações canônicas mais comuns são |Distribuição | Normal | Binomial | Poisson | Gama | N.Inversa| |:-----------:|:------------------:|:--------:|:-------:|:----:|:--------:| | Ligação | `\(\mu = \eta\)` | `\(\log\left( \frac{\mu}{1-\mu}\right) = \eta\)` | `\(\log \mu = \eta\)` | `\(\mu^{-1} = \eta\)` | `\(\mu^{-2} = \eta\)`| >**Importante**: As funções de ligação canônicas produzem propriedades estatísticas de interesse para o modelo, tais como, suficiência, facilidade de cálculo, unicidade das estimativas de máxima verossimilhança e, em alguns casos, interpretação simples. -- - Deve ser lembrado, porém, que embora as funções de ligação canônicas conduzam a propriedades estatísticas desejáveis para o modelo, principalmente, no caso de amostras pequenas, não há nenhuma razão a priori para que os efeitos sistemáticos do modelo devam ser aditivos na escala dada por tais funções. -- - Para ligações não-canônicas, Wedderburn (1976) discute condições à existência da concavidade de `\(L(\boldsymbol{\beta}; y)\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Outras ligações - Potência: `\(\eta = \mu^{\lambda}\)`, onde `\(\lambda\)` é um número real, principalmente para dados com média positiva; -- - Probit: `\(\eta = \Phi^{-1}(\mu)\)`, em que `\(\mu\)` é a proporção de sucessos de uma distribuição binomial e `\(\Phi(\cdot)\)` é a função de distribuição acumulada da distribuição normal padrão; -- - Logito: `\(\eta = \log[\mu/(1-\mu)]\)`, é a proporção de sucessos de uma distribuição binomial; -- - Complemento log-log: `\(\eta = \log[-\log(1-\mu)]\)`, é a proporção de sucessos de uma distribuição binomial; -- - Logaritmo: `\(\eta = \log[\mu]\)`, para dados com média positiva; -- - Box-Cox: `\(\eta = (\mu^{\lambda}-1)/\lambda\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Ligação de Aranda-Ordaz - I Uma outra transformação importante foi proposta por Aranda-Ordaz (1981) para dados binários. A transformação é dada por `$$\eta = \log\left[ \frac{(1-\mu)^{-\alpha} - 1}{\alpha} \right],$$` em que `\(0 < \mu < 1\)` e `\(\alpha\)` é uma constante desconhecida positiva. -- **Casos particulares**: - Quando `\(\alpha = 1\)` temos a função de ligação logito. -- - Quando `\(\alpha \rightarrow 0\)` temos a função de ligação complemento log-log. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Ligação de Aranda-Ordaz - II A figura mostra o comportamento de `\(\mu\)` (probabilidade) para alguns valores de `\(\alpha\)`. Em muitas situações práticas o interesse pode ser testarmos se o modelo logístico é apropriado, `\(H_0: \alpha = 1\)`, contra a necessidade de uma transformação na ligação, `\(H_1: \alpha \neq 1\)`. <div class="figure" style="text-align: center"> <img src="images/Aran.png" alt=" " width="40%" /> <p class="caption"> </p> </div> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Função Desvio - I Admitindo-se uma combinação satisfatória da distribuição da variável resposta e da função de ligação, o objetivo é determinar quantos termos são necessários na estrutura linear para uma descrição razoável dos dados. Um número grande de variáveis explanatórias pode conduzir a um modelo que explique bem os dados, mas com um aumento de complexidade na interpretação. Por outro lado, um número pequeno de variáveis explanatórias pode conduzir a um modelo de interpretação fácil, porém, que se ajuste pobremente aos dados. O que se deseja na realidade é um modelo intermediário, entre um modelo muito complicado e um modelo pobre em ajuste. A esse modelo denominamos **<span style="color:orange">modelo parcimonioso</span>**. -- > **Modelo Saturado** >Segundo Agresti (2007), “Let Ls denote the maximized log-likelihood value for the most complex model possible. This model has a separate parameter for each observation, and it provides a perfect fit to the data. The model is said to be saturated.” -- O modelo saturado tem, então, `\(n\)` parâmetros, um para cada observação, e as estimativas de MV das médias são `\(\widetilde{\mu_i} = y_i\)`, para `\(i = 1, ..., n\)`. O *til* é colocado para diferir das estimativas de MV do modelo sob investigação, cuja matriz modelo `\(\boldsymbol{X}\)` tem dimensões `\(n \times p\)`, com `\(p < n\)`. -- > Na prática, o modelo saturado é **pouco informativo** pois não sumariza os dados, mas, simplesmente, os repete. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Função Desvio - II A fim de discriminar entre modelos, medidas de discrepância devem ser introduzidas para medir o ajuste de um modelo e identificar a utilidade de um parâmetro extra no modelo sob pesquisa ou, então, verificar a falta de ajuste induzida pela omissão dele. -- Nelder e Wedderburn (1972) propuseram, como medida de discrepância, a **<span style="color:orange">deviance</span>** (traduzida como *desvio* por Cordeiro (1986)), que é uma distância entre o logaritmo da função de verossimilhança do modelo saturado (com `\(n\)` parâmetros) e do modelo sob investigação (com `\(p\)` parâmetros) avaliado na estimativa de máxima verossimilhança `\(\boldsymbol{\hat{\beta}}\)`. -- Sem perda de generalidade, vamos supor que o logaritmo da função de verossimilhança seja agora definido por `$$L(\boldsymbol{\mu}, \boldsymbol{y}) = \sum_{i = 1}^n L(\mu_i; y_i),$$` em que `\(\mu_i = g^{-1}(\eta_i)\)` e `\(\eta_i = \boldsymbol{x}_i^T \boldsymbol{\beta}\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Função Desvio - III Para o modelo saturado `\((p = n)\)`, a função `\(L(\boldsymbol{\mu}, \boldsymbol{y})\)` é estimada por `$$L(\boldsymbol{y}, \boldsymbol{y}) = \sum_{i = 1}^n L(y_i; y_i).$$` Ou seja, a estimativa de máxima verossimilhança de `\(\mu_i\)` fica nesse caso dada por `\(\mu_i = y_i\)`. Quando `\(p < n\)`, denotamos a estimativa de `\(L(\boldsymbol{\mu}, \boldsymbol{y})\)` por `\(L(\boldsymbol{\hat{\mu}}, \boldsymbol{y})\)`. Aqui, a estimativa de máxima verossimilhança de `\(\mu_i\)` será dada por `\(\hat{\mu_i} = g^{-1}(\hat{\eta}_i)\)`, em que `\(\hat{\eta}_i = \boldsymbol{x}_i^T \hat{\boldsymbol{\beta}}\)`. -- A qualidade do ajuste de um MLG é avaliada através da função desvio escalonado `$$D^*(\boldsymbol{y}; \hat{\boldsymbol{\mu}}) = \phi D(\boldsymbol{y}; \hat{\boldsymbol{\mu}}) = 2 \{ L(\boldsymbol{y}, \boldsymbol{y}) - L(\hat{\boldsymbol{\mu}}, \boldsymbol{y}) \},$$` que é uma distância entre o logaritmo da função de verossimilhança do modelo saturado (com `\(n\)` parâmetros) e do modelo sob investigação (com `\(p\)` parâmetros) avaliado na estimativa de máxima verossimilhança `\(\boldsymbol{\hat{\beta}}\)`. `\(D(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\)` é denominado *desvio não escalonado* ou, simplesmente, **<span style="color:orange">desvio</span>**. -- > Um valor pequeno para a função desvio indica que, para um número menor de parâmetros, obtemos um ajuste tão bom quanto o ajuste com o modelo saturado. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Função Desvio - IV Denotando por `\(\hat{\theta}_i = \theta_i(\hat{\mu}_i)\)` e `\(\tilde{\theta}_i = \theta_i(\tilde{\mu}_i)\)` as estimativas de máxima verossimilhança de `\(\theta\)` para os modelos com `\(p\)` parâmetros `\((p < n)\)` e saturado `\((p = n)\)`, respectivamente, temos que a função `\(D(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\)` fica, alternativamente, dada por: `$$D(\boldsymbol{y}; \hat{\boldsymbol{\mu}}) = 2 \sum_{i = 1}^n \{y_i(\tilde{\theta}_i - \hat{\theta}_i) + (b(\hat{\theta}_i) - b(\tilde{\theta}_i)) \}.$$` -- - Denotamos `\(D(\boldsymbol{y}; \hat{\boldsymbol{\mu}}) = \sum_{i = 1}^n d^2(y_i, \hat{\mu}_i)\)` em que `\(d^2(y_i, \hat{\mu}_i)\)` será denominado componente do desvio não-escalonado. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Casos particulares - I - **Normal** Aqui `\(\theta_i = \mu_i\)`, logo `\(\tilde{\theta}_i = y_i\)` e `\(\hat{\theta}_i = \hat{\mu}_i\)`. O desvio fica portanto dado por: <div class="math"> \[\begin{align*} D(\boldsymbol{y}; \hat{\boldsymbol{\mu}}) &= 2 \sum_{i = 1}^n \left\{y_i (y_i - \hat{\mu}_i) + \frac{\hat{\mu}_i^2}{2} - \frac{y_i^2}{2}\right\} \\ &= \sum_{i = 1}^n (y_i - \hat{\mu}_i)^2. \end{align*}\] </div> que coincide com a soma de quadrado de resíduos. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Casos particulares - II - **Poisson** Aqui `\(\theta_i = \log(\mu_i)\)`, logo `\(\tilde{\theta}_i = \log(y_i)\)`, para `\(y_i > 0\)` e `\(\hat{\theta}_i = \log(\hat{\mu}_i)\)`. O desvio fica portanto dado por: <div class="math"> \[\begin{align*} D(\boldsymbol{y}; \hat{\boldsymbol{\mu}}) = 2 \sum_{i = 1}^n \left\{ y_i \log \left( \frac{y_i}{\hat{\mu}_i} \right) - (y_i - \hat{\mu}_i) \right\}. \end{align*}\] </div> -- Se `\(y_i = 0\)`, o i-ésimo termo de `\(D(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\)` vale `\(2 {\hat{\mu}_i}\)`. Resumindo, em relação ao componente do desvio do modelo Poisson, temos que: <div class="math"> \[\begin{align*} d^2(y_i, \hat{\mu}_i) = \left\{\begin{array}{ll} 2 \left[ y_i \log \left(\frac{y_i}{\hat{\mu}_i}\right) - (y_i - \hat{\mu}_i) \right], & \mbox{se} \quad y_i > 0;\\ 2 \; {\hat{\mu}_i}, & \mbox{se} \quad y_i = 0. \end{array} \right. \end{align*}\] </div> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Casos particulares - III - **Binomial** Aqui assumimos `\(Y_i \sim B(n_i, \mu_i)\)`, `\(i = 1, ..., k\)`. Lembrando que na família exponencial consideramos que `\(y_i^* = \frac{y_i}{n_i}\)` obtemos que - `\(\theta_i = \log \left(\frac{\mu_i}{1-\mu_i}\right)\)` e `\(b(\theta_i) = \log (1 + e^{\theta_i})\)` -- - `\(\hat{\theta}_i = \log \left(\frac{\hat{\mu}_i}{1-\hat{\mu}_i}\right)\)` e `\(b( \hat{\theta}_i) = - \log(1-\hat{\mu}_i)\)` -- - `\(\tilde{\theta}_i = \log \left( \frac{y_i/n_i}{1-y_i/n_i}\right)\)` e `\(b(\tilde{\theta}_i) = -\log(1-y_i/n_i)\)`. -- Logo, o desvio assume a forma: <div class="math"> \[\begin{align*} D(\boldsymbol{y}; \hat{\boldsymbol{\mu}}) = 2 \sum_{i = 1}^n y_i \log \left( \frac{y_i}{n_i \hat{\mu}_i} \right) + (n_i - y_i) \log \left( \frac{1 - y_i/n_i}{1 - \hat{\mu}_i} \right), \quad \mbox{se} \quad 0 < y_i < n_i. \end{align*}\] </div> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Casos particulares - IV No entanto, nos casos em que `\(y_i = 0\)` ou `\(y_i = n\)`, o i-ésimo termo do `\(D(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\)` pode ser reescrito de uma maneira mais simples, sendo assim, <div class="math"> \[\begin{align*} d^2(y_i, \hat{\mu}_i) = \left\{\begin{array}{ll} 2 y_i \log \left( \frac{y_i}{n_i \hat{\mu}_i} \right) + (n_i - y_i) \log \left( \frac{1 - y_i/n_i}{1 - \hat{\mu}_i} \right), & \mbox{se} \quad 0 < y_i < n_i; \\ -2 \; n_i \log(1-{\hat{\mu}_i}), & \mbox{se} \quad y_i = 0;\\ -2 \; n_i \log({\hat{\mu}_i}), & \mbox{se} \quad y_i = n_i. \end{array} \right. \end{align*}\] </div> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Observações: - Um valor pequeno para a função desvio indica que, para um número menor de parâmetros, obtém-se um ajuste tão bom quanto o ajuste com o modelo saturado. -- - Embora seja usual compararmos os valores observados da função desvio com os percentis da distribuição qui-quadrado com `\(n\)` − `\(p\)` graus de liberdade, em geral `\(D(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\)` não segue assintoticamente uma `\(\chi_{n-p}^2\)`. -- - No caso binomial, quando `\(k\)` é fixo e `\(n_i \rightarrow \infty\)`, para cada `\(i\)`, `\(D(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\)` segue, sob a hipótese de que o modelo é verdadeiro, uma `\(\chi_{k-p}^2\)`. Isso não vale quando `\(n \rightarrow \infty\)` e `\(n_i \mu_i (1 - \mu_i)\)` permanece limitado. -- - Em geral, para os casos em que `\(D^*(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\)` depende do parâmetro de dispersão, `\(\phi^{-1}\)`, o seguinte resultado (Jørgensen, 1987) para a distribuição nula da função desvio pode ser utilizado: <div class="math"> \[\begin{align} D^*(\boldsymbol{y}; \hat{\boldsymbol{\mu}}) \sim \chi_{n-p}^2, \;\;\; \mbox{quando} \;\;\; \phi \rightarrow \infty \end{align}\] </div> Isto é, quando a dispersão é pequena, fica razoável compararmos os valores observados de `\(D^*(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\)` com os percentis da `\(\chi_{n-p}^2\)`. -- - Em particular, para o modelo normal linear, o resultado acima nos diz `\(\sum_{i = 1}^n \frac{(y_i -\hat{\mu}_i)^2}{\sigma^2}\)`, quando `\(\sigma^2 \rightarrow 0\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Análise do desvio - I Vamos supor para o vetor de parâmetros `\(\boldsymbol{\beta}\)` a partição `\(\boldsymbol{\beta} = (\boldsymbol{\beta}_1^T, \boldsymbol{\beta}_2^T)^T\)`, em que `\(\boldsymbol{\beta}_1\)` é um vetor `\(q\)`-dimensional enquanto `\(\boldsymbol{\beta}_2\)` tem dimensão `\(p-q\)` e `\(\phi\)` é conhecido (ou fixo). Pdemos ter interesse em testar as hipóteses `\(\mbox{H}_0: \;\; \boldsymbol{\beta}_1 = 0\)` contra `\(\mbox{H}_1: \;\; \boldsymbol{\beta}_1 \neq 0\)`. As funções desvio correspondentes aos modelos sob `\(\mbox{H}_0\)` e `\(\mbox{H}_1\)` serão denotadas por `\(D(\boldsymbol{y}; \hat{\boldsymbol{\mu}}^0)\)` e `\(D(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\)`, respectivamente, em que `\(\hat{\boldsymbol{\mu}}^0\)` é a estimativa de máxima verossimilhança sob `\(\mbox{H}_0\)`. -- A estatística da razão de verossimilhanças fica nesse caso dada por: <div class="math"> \[\begin{align} \xi_{RV} = \phi \{D(\boldsymbol{y}; \hat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\}, \label{eq:TRV} \end{align}\] </div> isto é, a diferença entre dois desvios. Como é conhecido, sob a hipótese nula, `\(\xi_{RV} \sim \chi_{n-p}^2\)`, quando `\(n \rightarrow \infty\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Análise do desvio - II De forma similar, podemos definir a estatística <div class="math"> \[\begin{align} F = \frac{\{D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0)-D(\boldsymbol{y};\hat{\boldsymbol{\mu}})\}/q}{D(\boldsymbol{y}; \hat{\boldsymbol{\mu}})/(n-p)}, \label{eq:Ftest} \end{align}\] </div> cuja distribuição nula assintótica é uma `\(F_{q,(n-p)}\)` quando o denominador de (\ref{eq:Ftest}) é uma estimativa consistente de `\(\phi^{-1}\)` (ver, por exemplo, Jørgensen, 1987). --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Análise do desvio - III **Esquema da análise do desvio com dois fatores na parte sistemática** | Modelo | Desvio | Diferença | G.L. | Testando | |:----------:|:-----------:|:-----------------:|:----------------------------:|:------------------------:| | Constante | `\(D_0\)` | | | | | | | `\(D_0-D_A\)` | `\(n(A)-1\)` | `\(A\)` ignorando `\(B\)` | | | | `\(D_0-D_B\)` | `\(n(B)-1\)` | `\(B\)` ignorando `\(A\)` | | `\(+A\)` | `\(D_A\)` | | | | | | | `\(D_A-D_{A+B}\)` | `\(n(B)-1\)` | `\(B\mid A\)` ignorando `\(AB\)` | | `\(+B\)` | `\(D_B\)` | | | | | | | `\(D_B-D_{A+B}\)` | `\(n(A)-1\)` | `\(A\mid B\)` ignorando `\(AB\)` | | `\(+A+B\)` | `\(D_{A+B}\)` | | | | | | | `\(D_{A+B}-D_{AB}\)` | `\(\{n(A)-1\}\times\{n(B)-1\}\)`| `\(AB\mid A + B\)` | | `\(+A+B+AB\)` | `\(D_{AB}\)` | | | . | --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Quando utilizar o teste da razão de verossimilhança, quando utilizar o teste F? - A vantagem de utilizarmos (\ref{eq:Ftest}) em relação a (\ref{eq:TRV}) é que a estatística F não depende do parâmetro de dispersão. -- - O resultado (\ref{eq:Ftest}) também é verificado quando `\(\phi \rightarrow \infty\)` e `\(n\)` é arbitrário. -- - Quando `\(\phi\)` é desconhecido a estatística da razão de verossimilhanças assume uma expressão diferente de (\ref{eq:TRV}). --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo de Análise do desvio - I Os dados se referem a um estudo de caso-controle realizado no Setor de Anatomia e Patologia do Hospital Heliópolis em São Paulo, no período de 1970 a 1982 (Paula e Tuder, 1986). Um total de 175 pacientes com processo infecioso pulmonar atendido no hospital no período acima foi classificado segundo as seguintes variáveis: - Tipo, tipo de tumor (1: maligno, 0: benigno); - IDADE, idade em anos; - SEXO (0: masculino, 1: feminino); - HL, intensidade da célula histiócitos-linfócitos (1: ausente, 2: discreta, 3: moderada, 4: intensa); - FF, intensidade da célula fibrose-frouxa (1: ausente, 2: discreta, 3: moderada, 4: intensa). --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo de Análise do desvio - II
--- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo de Análise do desvio - III ``` r fit0 <- glm(Tipo ~ 1, family = binomial, data = dados_cancer) summary(fit0) ``` ``` ## ## Call: ## glm(formula = Tipo ~ 1, family = binomial, data = dados_cancer) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -0.3817 0.1539 -2.479 0.0132 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 236.34 on 174 degrees of freedom ## Residual deviance: 236.34 on 174 degrees of freedom ## AIC: 238.34 ## ## Number of Fisher Scoring iterations: 4 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo de Análise do desvio - IV ``` r fit1 <- glm(Tipo ~ Sexo, family = binomial, data = dados_cancer) summary(fit1) ``` ``` ## ## Call: ## glm(formula = Tipo ~ Sexo, family = binomial, data = dados_cancer) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -0.8484 0.4654 -1.823 0.0683 . ## Sexo 0.3629 0.3400 1.067 0.2858 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 236.34 on 174 degrees of freedom ## Residual deviance: 235.21 on 173 degrees of freedom ## AIC: 239.21 ## ## Number of Fisher Scoring iterations: 4 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn, scrollable-slide ### Exemplo de Análise do desvio - V ``` r fit2 <- glm(Tipo ~ Sexo + Idade, family = binomial, data = dados_cancer) summary(fit2) ``` ``` ## ## Call: ## glm(formula = Tipo ~ Sexo + Idade, family = binomial, data = dados_cancer) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -4.68726 0.88711 -5.284 1.27e-07 *** ## Sexo 0.65520 0.40205 1.630 0.103 ## Idade 0.06954 0.01201 5.791 7.01e-09 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 236.34 on 174 degrees of freedom ## Residual deviance: 188.22 on 172 degrees of freedom ## AIC: 194.22 ## ## Number of Fisher Scoring iterations: 4 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo de Análise do desvio - VI ``` r fit3 <- glm(Tipo ~ Sexo + Idade + factor(HL), family = binomial, data = dados_cancer) summary(fit3) ``` ``` ## ## Call: ## glm(formula = Tipo ~ Sexo + Idade + factor(HL), family = binomial, ## data = dados_cancer) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -2.94989 1.21054 -2.437 0.01482 * ## Sexo 0.93736 0.45098 2.078 0.03767 * ## Idade 0.06534 0.01307 5.000 5.75e-07 *** ## factor(HL)2 -1.06661 0.95171 -1.121 0.26240 ## factor(HL)3 -2.64465 0.95431 -2.771 0.00558 ** ## factor(HL)4 -3.98225 1.48109 -2.689 0.00717 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 236.34 on 174 degrees of freedom ## Residual deviance: 162.55 on 169 degrees of freedom ## AIC: 174.55 ## ## Number of Fisher Scoring iterations: 5 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo de Análise do desvio - VII ``` r fit4 <- glm(Tipo ~ Sexo + Idade + factor(HL) + factor(FF), family = binomial, data = dados_cancer) summary(fit4) ``` ``` ## ## Call: ## glm(formula = Tipo ~ Sexo + Idade + factor(HL) + factor(FF), ## family = binomial, data = dados_cancer) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -2.63437 1.22220 -2.155 0.0311 * ## Sexo 0.78392 0.47066 1.666 0.0958 . ## Idade 0.06513 0.01328 4.906 9.31e-07 *** ## factor(HL)2 -0.86926 0.94750 -0.917 0.3589 ## factor(HL)3 -2.24903 0.97178 -2.314 0.0206 * ## factor(HL)4 -3.29561 1.47997 -2.227 0.0260 * ## factor(FF)2 -0.68718 0.50397 -1.364 0.1727 ## factor(FF)3 -1.02472 0.52721 -1.944 0.0519 . ## factor(FF)4 0.43149 1.12452 0.384 0.7012 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 236.34 on 174 degrees of freedom ## Residual deviance: 157.40 on 166 degrees of freedom ## AIC: 175.4 ## ## Number of Fisher Scoring iterations: 5 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo de Análise do desvio - VIII ANODEV ``` r anova(fit0,fit1) ``` ``` ## Analysis of Deviance Table ## ## Model 1: Tipo ~ 1 ## Model 2: Tipo ~ Sexo ## Resid. Df Resid. Dev Df Deviance ## 1 174 236.34 ## 2 173 235.21 1 1.1354 ``` ``` r anova(fit1,fit2) ``` ``` ## Analysis of Deviance Table ## ## Model 1: Tipo ~ Sexo ## Model 2: Tipo ~ Sexo + Idade ## Resid. Df Resid. Dev Df Deviance ## 1 173 235.21 ## 2 172 188.22 1 46.982 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo de Análise do desvio - IX ``` r anova(fit2,fit3) ``` ``` ## Analysis of Deviance Table ## ## Model 1: Tipo ~ Sexo + Idade ## Model 2: Tipo ~ Sexo + Idade + factor(HL) ## Resid. Df Resid. Dev Df Deviance ## 1 172 188.22 ## 2 169 162.55 3 25.674 ``` ``` r anova(fit3,fit4) ``` ``` ## Analysis of Deviance Table ## ## Model 1: Tipo ~ Sexo + Idade + factor(HL) ## Model 2: Tipo ~ Sexo + Idade + factor(HL) + factor(FF) ## Resid. Df Resid. Dev Df Deviance ## 1 169 162.55 ## 2 166 157.40 3 5.1476 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo de Análise do desvio - X Em resumo: | Modelo | Desvio | Diferença | G.L. | Testando | |:----------:|:-----------:|:-----------------:|:--------:|:----------------------------:| | Constante | 236,34 | | | | | +SEXO | 235,20 | 1,14 | 1 | SEXO | | +IDADE | 188,22 | 46,98 | 1 | IDADE `\(\mid\)` SEXO | | +HL | 162,55 | 25,67 | 3 | HL `\(\mid\)` SEXO + IDADE | | +FF | 157,40 | 5,15 | 3 | FF `\(\mid\)` SEXO + IDADE + HL| --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Estimação dos parâmetros - I Considere a partição `\(\boldsymbol{\theta} = (\boldsymbol{\beta}^T, \phi)^T\)` e denote o logaritmo da função de verossimilhança por `\(L(\theta;y)\)`. Para obter a função escore para o parâmetro `\(\beta\)`, inicialmente calculamos as derivadas <div class="math"> \[\begin{align} \frac{\partial L(\boldsymbol{\theta})}{\partial \beta} = \sum_{i=1}^n\phi \left\{y_i\frac{\partial \theta_i}{\partial \beta_j} - \frac{\partial b(\theta_i)}{\partial \beta_j}\right\}. \end{align}\] </div> -- No entanto, note que `$$L(\boldsymbol{\theta}) = f(\theta_1, \theta_2,...,\theta_i,...,\theta_n)$$` `$$\downarrow$$` `$$\theta_i = \int V^{-1} _i d \mu_i = q(\mu_i)$$` `$$\downarrow$$` $$\mu_i = g^{-1}(\eta_i) $$ `$$\downarrow$$` `$$\eta_i = \sum^p _{j=1} x_{ij}\beta_j$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Estimação dos parâmetros - II Para obtermos a função escore para o parâmetro `\(\boldsymbol{\beta}\)` calculamos inicialmente as derivadas <div class="math"> \[\begin{align*} U_j = \frac{\partial L(\boldsymbol{\theta})}{\partial\beta_j} &= \sum^n_{i=1} \phi\left\{y_i\frac{d\theta_i}{d\mu_i}\frac{d\mu_i}{d\eta_i}\frac{\partial\eta_i}{\beta_j} -\frac{db(\theta_i)}{d\theta_i}\frac{d\theta_i}{d\mu_i}\frac{d\mu_i}{d\eta_i}\frac{\partial\eta_i}{\partial\beta_j}\right\}\\ &= \sum^n_{i=1}\phi \left\{y_i V^{-1}_i \left(\frac{d\mu_i}{d \eta_i}\right)x_{ij} - \mu_i V^{-1}_i\left(\frac{d \mu_i}{d \eta_i}\right)x_{ij}\right\}\\ &=\sum^n_{i=1} \phi \left\{\sqrt{\frac{\omega_i}{V_i}}(y_i - \mu_i)x_{ij}\right\}, \end{align*}\] </div> em que `\(\omega_i = (d\mu_i/d\eta_i)^2/V_i\)` e `\(j = 1, ..., p\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Estimação dos parâmetros - III Então, a função escore na forma matricial fica na forma <div class="math"> \[\begin{align*} U_\beta(\theta) = \frac{\partial L(\theta)}{\partial \beta} = \phi X^T W^{1/2}V^{-1/2}(y-\mu), \end{align*}\] </div> onde - `\(X\)` é uma matriz `\(n \times p\)` de posto completo, cujas linhas serão denotadas por `\(x^T_i\)`, $i = 1,\ldots,n $, - `\(W\)` = `\(diag \left\{\omega_i,\ldots,\omega_n\right\}\)` é a matriz de pesos, - `\(V\)` = `\(diag \left\{V_i,\ldots,V_n\right\}\)`, - `\(\boldsymbol{y}\)` = `\((y_1,\ldots,y_n)^T\)` e - `\(\boldsymbol{\mu}\)` = `\((\mu_1,\ldots,\mu_n)^T\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Estimação dos parâmetros - IV Na busca da matriz de informação de Fisher precisamos das derivadas <div class="math"> \[\begin{align*} \frac{\partial^2 L(\boldsymbol{\theta})}{\partial \beta_j \partial\beta_l} &= \phi \sum^n_{i=1}(y_i -\mu_i)\frac{d^2\theta_i}{d\mu^2_i}\left(\frac{d\mu_i}{d \eta_i}\right)^2 x_{ij}x_{il}\\ &+ \phi \sum^n_{i=1}(y_i - \mu_i)\frac{d\theta_i}{d\mu_i}\frac{d^2 \mu_i}{d \eta^2 _i} x_{ij}x_{il}\\ &- \phi \sum^n _{i=1}\frac{d\theta_i}{d\mu_i}\left(\frac{d\mu_i}{d \eta_i}\right)^2 x_{ij}x_{il}, \end{align*}\] </div> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Estimação dos parâmetros - V cujos valores esperados são dados por <div class="math"> \[\begin{align*} E \left\{\frac{\partial^2 L(\boldsymbol{\theta})}{\partial\beta_j\partial\beta_l}\right\} &= {-\phi}\sum^n_{i=1}\frac{d\theta_i}{d\mu_i}\left(\frac{d\mu_i}{d\eta_i}\right)^2 x_{ij}x_{il}\\ &= {-\phi}\sum^n _{i=1}\frac{\left(d\mu_i/d\eta_i\right)^2}{V_i} x_{ij}x_{il}\\ &= {-\phi}\sum^n _{i=1}\omega_i x_{ij}x_{il}. \end{align*}\] </div> -- Na forma matricial, a **<span style="color:orange">informação de Fisher</span>** é denotada por <div class="math"> \[\begin{align*} K_{\beta\beta}(\boldsymbol{\theta}) = E \left\{-\frac{\partial^2 L(\theta)}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^T}\right\} = \phi X^T WX. \end{align*}\] </div> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Estimação dos parâmetros - VI Para ligação canônica, `\(\theta_i=\eta_i\)` e `$$V_i = V(\mu_i) = \frac{\partial\mu_i}{\partial\theta_i} = \frac{\partial\mu_i}{\partial\eta_i}.$$` Assim, `$$\omega_i = \frac{(d\mu_i/d\eta_i)^2}{V_i} = \frac{\partial\mu_i}{\partial\eta_i} = V_i.$$` -- E, dessa forma, temos que `$$U_\beta = \phi X^T (y - \mu)$$` e `$$K_{\beta\beta} = \phi X^T VX.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Estimação de `\(\phi\)` - I A função escore para o parâmetro `\(\phi\)` fica dada por <div class="math"> \[\begin{align*} U_\phi &= \frac{\partial L(\boldsymbol{\theta})}{\partial \phi}\\ &= \sum^n _{i=1}\left\{y_i \theta_i - b(\theta_i)\right\} + \sum^n _{i=1} c'(y_i,\phi) \end{align*}\] </div> em que `\(c'(y_i,\phi) = \partial c(y_i, \phi)/\partial \phi.\)` -- Para obtermos a informação de Fisher para `\(\phi\)` temos que calcular `$$\frac{\partial^2 L(\boldsymbol{\theta})}{\partial \phi^2} = \sum^n _{i=1} c''(y_i,\phi),$$` em que `\(c''(y_i,\phi) = \partial^2 c(y_i, \phi)/\partial \phi^2.\)` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Estimação de `\(\phi\)` - II Assim, a informação de Fisher para `\(\phi\)` fica dada por `$$K_{\phi\phi} = - \sum^n _{i=1} E\{c''(Y_i,\phi)\}.$$` -- **Ortogonalidade** Temos que <div class="math"> \[\begin{align*} \frac{\partial^2 L(\boldsymbol{\theta})}{\partial\boldsymbol{\beta}\partial \phi} = \sum^n _{i=1} \phi \left\{\sqrt{\frac{\omega_i}{V_i}}(y_i - \mu_i){\mbox x}_{i}\right\}. \end{align*}\] </div> Portanto, verificamos que `\(\boldsymbol{\beta}\)` e `\(\phi\)` são ortogonais, isto é, `$$K_{\boldsymbol{\theta}\phi} = E\left[− \frac{\partial^2 L(\boldsymbol{\theta})}{\partial\boldsymbol{\beta}\partial\phi}\right] = 0.$$` Logo, segue que a matriz de informação de Fisher para `\(\boldsymbol{\theta}\)` é bloco diagonal, sendo dada por `\(K_{\boldsymbol{\theta}\boldsymbol{\theta}} = diag\{K_{\beta\beta},K_{\phi\phi}\}\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Estimação de `\(\beta\)` - I O processo iterativo de Newton-Raphson para a obtenção da estimativa de máxima verossimilhança de `\(\boldsymbol{\beta}\)` é definido expandindo a função escore `\(U_{\boldsymbol{\beta}}\)` em torno de um valor inicial `\(\boldsymbol{\beta}^{(0)}\)`, tal que <div class="math"> \[\begin{align*} U_{\boldsymbol{\beta}} \cong {U}_{\boldsymbol{\beta}}^{(0)} + {U'}_{\boldsymbol{\beta}}^{(0)} (\boldsymbol{\beta} - \boldsymbol{\beta}^{(0)}), \end{align*}\] </div> em que `\({U'}_{\boldsymbol{\beta}}^{(0)}\)` denota a primeira derivada de `\(U_{\boldsymbol{\beta}}\)` com respeito a `\(\boldsymbol{\beta}^T\)`, sendo `\({U'}_{\boldsymbol{\beta}}^{(0)}\)` e `\({U}_{\boldsymbol{\beta}}^{(0)}\)`, respectivamente, essas quantidades avaliadas em `\(\boldsymbol{\beta}^{(0)}\)`. Assim, repetindo o procedimento acima, chegamos ao processo iterativo <div class="math"> \[\begin{align*} \boldsymbol{\beta}^{(m+1)} = \boldsymbol{\beta}^{(m)} + \{ {({-U'}_{\boldsymbol{\beta}}}^{-1}) \}^{(m)} {U}_{\boldsymbol{\beta}}^{(m)}, \end{align*}\] </div> `\(m = 0, 1, \ldots\)`. Como a matriz `\({-U'}_{\boldsymbol{\beta}}\)` pode não ser positiva definida, a aplicação do método escore de Fisher substituindo a matriz `\({-U'}_{\boldsymbol{\beta}}\)` pelo correspondente valor esperado `\(K_{\beta\beta}\)` pode ser mais conveniente. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Estimação de `\(\beta\)` - II Isso resulta em <div class="math"> \[\begin{align*} \boldsymbol{\beta}^{(m+1)} = \boldsymbol{\beta}^{(m)} + \{ {(K_{\beta\beta}}^{-1}) \}^{(m)} {U}_{\boldsymbol{\beta}}^{(m)}, \end{align*}\] </div> `\(m = 0, 1, \ldots\)`. Se trabalharmos um pouco o lado direito da expressão acima, chegaremos a <div class="math"> \[\begin{align*} \boldsymbol{\beta}^{(m+1)} = \boldsymbol{\beta}^{(m)} + \left({X^{'}} W^{(m)}X\right)^{-1}{X'} \left(W^{(m)}\right)^{1/2}\left(V^{(m)}\right)^{-1/2} (\boldsymbol{y} - \boldsymbol{\mu}^{(m)}), \end{align*}\] </div> `\(m = 0, 1,\ldots\)`. É possível ainda provar que o processo acima pode ser escrito como <div class="math"> \[\begin{align} \boldsymbol{\beta}^{(m+1)} = \left({X'}W^{(m)}X\right)^{-1} {X'} W^{(m)} \boldsymbol{z}^{(m)}, \label{eq:repond} \end{align}\] </div> `\(m = 0, 1,\ldots\)`, em que `\(\boldsymbol{z} = \boldsymbol{\eta} + W^{-1/2} V^{-1/2} (\boldsymbol{y} - \boldsymbol{\mu})\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Estimação de `\(\beta\)` - III A equação matricial (\ref{eq:repond}) é válida para qualquer MLG e mostra que a solução das equações de MV equivale a calcular repetidamente uma regressão linear ponderada de uma variável dependente ajustada `\(\boldsymbol{z}\)` sobre a matriz `\(X\)` usando uma função de peso `\(W\)` que se modifica no processo iterativo. As funções de variância e de ligação entram no processo iterativo através de `\(W\)` e `\(\boldsymbol{z}\)`. -- A convergência de (\ref{eq:repond}) ocorre em geral em um número finito de passos, independente dos valores iniciais utilizados (Wedderburn, 1976). -- É usual iniciar (\ref{eq:repond}) com `\(\eta_i^{(0)} = g(y_i)\)`, para `\(i = 1,\ldots, n\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Estimação de `\(\phi\)` - III Igualando a função escore `\(U_{\phi}\)` a zero chegamos à seguinte solução: <div class="math"> \[\begin{align*} \sum^n _{i=1} c'(y_i,\hat{\phi}) = \frac{1}{2} D(\boldsymbol{y}; \hat{\boldsymbol{\mu}}) - \sum^n_{i=1}\left\{y_i \tilde{\theta_i} - b(\tilde{\theta_i})\right\}, \end{align*}\] </div> -- em que `\(D(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\)` denota o desvio do modelo sob investigação. Verificamos que a estimativa de máxima verossimilhança para `\({\phi}\)` nos casos normal e normal inversa, igualando `\(U_{\phi}\)` a zero, é dada por <div class="math"> \[\begin{align*} \hat{\phi} = \frac{n}{D(\boldsymbol{y}; \hat{\boldsymbol{\mu}})}. \end{align*}\] </div> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Propriedades dos Estimadores de um MLG - Por se tratarem de estimadores de máxima verossimilhança, os estimadores dos parâmetros de um MLG compartilham as propriedades gerais dos EMVs. -- - Algo que será recorrente deste ponto em diante é que boa parte dos resultados são assintóticos, o que quer dizer que são rigorosamente válidos quando `\(n \rightarrow \infty\)`, mas aproximadamente válidos para grandes amostras. -- - Condições adicionais para a validade dos resultados serão oportunamente relatadas, caso se apliquem. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Distribuição assintótica dos estimadores - A distribuição dos estimadores é fundamental para se dimensionar o erro das estimativas, construir intervalos de confiança e testar hipóteses de interesse. -- - Seja `\(\hat{\boldsymbol{\beta}}\)` o estimador de máxima verossimilhança do vetor de parâmetros `\(\boldsymbol{\beta}\)` de um MLG. A distribuição assintótica de `\(\hat{\boldsymbol{\beta}}\)` é dada por: `$$\hat{\boldsymbol{\beta}} \sim N(\boldsymbol{\beta}, K^{-1}_{\beta\beta}),$$` sendo `\(K_{\beta\beta} (\boldsymbol{\theta}) = \phi X^T W X\)` e `\(W\)` a matriz diagonal com elementos `\(diag \left\{\omega_i,\ldots,\omega_n\right\}\)`. -- - Na prática, a matriz `\(K_{\beta\beta}\)` é estimada, usando-se alguma estimativa para `\(\phi\)` e a matriz `\(W\)` é avaliada com base nas estimativas obtidas no último passo do algoritmo de estimação. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Propriedades dos Estimadores de Máxima Verossimilhança - I Por se tratarem de estimadores de máxima verossimilhança, as seguintes propriedades são garantidas para cada `\(\hat{\beta}_j\)`, `\(j = 0, 1, \ldots, p\)`. - São assintoticamente não viciados: <div class="math"> \[\begin{align*} \lim_{n \rightarrow \infty} E(\hat{\beta}_j) = \beta_j. \end{align*}\] </div> -- - São estimadores consistentes: <div class="math"> \[\begin{align*} P( \mid \hat{\beta}_j - \beta_j \mid < \epsilon) \xrightarrow{n \rightarrow \infty} 0, \quad \text{para qualquer } \epsilon > 0. \end{align*}\] </div> -- - Apresentam distribuição assintótica normal. -- - São eficientes (apresentam variância mínima) na classe dos estimadores (assintoticamente) não viciados. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Propriedades dos Estimadores de Máxima Verossimilhança - II Como `\(K_{\boldsymbol{\theta}\boldsymbol{\theta}} = diag\{K_{\beta\beta},K_{\phi\phi}\}\)` (bloco diagonal), então assintoticamente segue que -- - `\(\hat{\boldsymbol{\beta}} \sim N_{p+1}(\beta, K^{-1} _{\beta\beta})\)`, -- - `\(\hat{\phi} \sim N(0,K^{-1} _{\phi\phi})\)` e -- - `\(\hat{\boldsymbol{\beta}}\)` e `\(\hat{\phi}\)` são independentes. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Teste de hipóteses Em MLGs, sob a hipótese nula `\(H_0\)`, três estatísticas são utilizadas para testar as hipóteses relativas aos parâmetros `\(\beta's\)`: - A estatística da razão de verossimilhança; -- - A estatística de Wald; -- - A estatística escore; -- que assintoticamente são equivalentes. >Sob a hipótese nula `\(H_0\)` e supondo que o parâmetro de dispersão `\(\phi\)` é conhecido, as três estatísticas convergem para uma variável aleatória com distribuição `\(\chi^2_{p}\)`, sendo a razão de verossimilhanças o critério que define um teste uniformemente mais poderoso. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Teste da razão de verossimilhanças Supondo `\(\phi\)` conhecido em um MLG, considere o teste da hipótese simples `$$H_0: \boldsymbol{\beta} = \boldsymbol{\beta}^0\quad \text{vs} \quad H_1: \boldsymbol{\beta} \neq \boldsymbol{\beta}^0,$$` onde `\(\boldsymbol{\beta}^0\)` é um vetor específico do vetor `\(\boldsymbol{\beta}\)`, `\(p-\)`dimensional e assumido conhecido. -- O TRV é usualmente definido por <div class="math"> \[\begin{align*} \xi_{RV} = 2 \left\{ L(\hat{\boldsymbol{\beta}}) - L(\boldsymbol{\beta}^0) \right\}, \end{align*}\] </div> -- e pode ser obtido como uma diferença entre duas funções desvio como segue: <div class="math"> \[\begin{align*} \xi_{RV} = \phi \left\{D(\boldsymbol{y};\widehat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\widehat{\boldsymbol{\mu}})\right\}, \end{align*}\] </div> em que, `\(\widehat{\boldsymbol{\mu}}^0 = g^{-1} (\widehat{\boldsymbol{\eta}}^0)\)`, `\(\widehat{\boldsymbol{\eta}}^0 = X\boldsymbol{\beta}^0\)`. -- Assintoticamente e sob `\(H_0\)` temos `\(\xi_{RV} \sim \chi^2 _p\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Modelos encaixados Queremos testar apenas um subconjunto do vetor de parâmetros `\(\boldsymbol{\beta}\)`. Considere a partição `\(\boldsymbol{\beta} = (\boldsymbol{\beta}^T _1\;\boldsymbol{\beta}^T _2)^T\)`, em que `\(\boldsymbol{\beta}_1\)` é um vetor `\(q-\)`dimensional e `\(\boldsymbol{\beta}_2\)` tem dimensão `\(p-q\)` dimensional. Supondo `\(\phi\)` conhecido em um MLG, considere as seguintes hipóteses: `$$H_0: \beta_1 = \beta^0_{1}\quad \text{ vs }\quad H_1: \beta_1 \neq \beta^0_{1}$$` -- e, teremos que, <div class="math"> \[\begin{align*} \xi_{RV} = \phi [D(\boldsymbol{y};\widehat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\widehat{\boldsymbol{\mu}})], \end{align*}\] </div> em que `\(\widehat{\boldsymbol{\mu}}^0\)` é a estimativa de máxima verossimilhança de MLG com parte sistemática `\(\boldsymbol{\eta} = \widehat{\boldsymbol{\eta}}^0 _1 + \boldsymbol{\eta}_2\)`, em que `\(\widehat{\boldsymbol{\eta}}^0 _1 = \sum^q _{j=1} x_j \beta^0 _j\)` e `\(\boldsymbol{\eta}_2 = \sum^p _{j=q+1} x_j \beta_j\)`. -- A quantidade `\(\widehat{\boldsymbol{\eta}}^0 _1\)` desempenha o papel de um *offset* (parte conhecida no preditor linear), conforme nomenclatura dos modelos lineares generalizados. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Modelos encaixados - Exemplo 1 Para ilustrarmos a utilização do *offset*, vamos supor um modelo de Poisson com ligação log-linear, resposta resp, covariáveis cov1 e cov2 e *offset* dado por logt0. -- Para ajustarmos o modelo e armazenarmos os resultados em fit1.poisson devemos fazer ```` fit1.poisson = glm(resp ~ cov1 + cov2 + offset(logt0), family= poisson) ```` Esse tipo de recurso é muito utilizado em estudos de seguimento em que cada indivíduo é observado durante um tempo diferente. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Modelos encaixados - Exemplo 2 Vamos supor um MLG com distribuição normal inversa, ligação canônica e preditor linear dado por $$\eta = \beta_1 + \beta_2 \text{cov2} + \beta_3 \text{cov3} $$ e que o interesse é testarmos `\(H_0 : \beta_2 = b\)`, em que `\(b\)` é uma constante diferente de zero, contra `\(H_1 : \beta_2 \neq b\)`. -- Sob `\(H_0\)` ```` fit0.ni = glm(resp ~ cov3 + offset(b cov2), family= inverse.gaussian). ```` -- Sob `\(H_1\)` ```` fit1.ni = glm(resp ~ cov 2 + cov3, family= inverse.gaussian). ```` -- Como `\(\phi\)` deve ser estimado, a estatística `\(F\)` deve ser usada para testarmos `\(H_0 : \beta_2 = b\)` versus `\(H_1 : \beta_2 \neq b\)`. ```` anova(fit0.ni, fit1.ni, test = "F"). ```` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Teste de Wald - I Para testarmos `\(H_0\)`, a estatística de Wald fica expressa na forma <div class="math"> \[\begin{align*} \xi_W = [\widehat{\boldsymbol{\beta}}_1 - \boldsymbol{\beta}^0_1]^T\; Var^{-1}(\widehat{\boldsymbol{\beta}}_1) [\widehat{\boldsymbol{\beta}}_1 - \boldsymbol{\beta}^0_1], \end{align*}\] </div> -- em que `\(\widehat{\boldsymbol{\beta}}_1\)` sai do vetor `\(\widehat{\boldsymbol{\beta}} = (\widehat{\boldsymbol{\beta}}_1^T, \widehat{\boldsymbol{\beta}}_2^T)^T\)`, em que, assintoticamente e sob `\(H_0\)`, temos `\(\xi_{RV} \sim \chi^2 _p\)`. -- **Observações** - Embora possam ser aplicados ao teste de hipóteses de dois ou mais parâmetros, o uso mais frequente do teste de Wald contempla apenas um parâmetro por vez. Em situações envolvendo mais parâmetros, é mais usual aplicar o teste da razão de verossimilhanças. -- - O teste de Wald baseia-se na distribuição assintótica normal dos estimadores de máxima verossimilhança dos parâmetros do modelo. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Teste de Wald - II - Seja `\(\hat{\beta}_j\)` o estimador de máxima verossimilhança de `\(\beta_j\)`, um particular parâmetro de um MLG. Conforme discutido anteriormente, para `\(n \rightarrow \infty\)`, $$ \hat{\beta}_j \sim N\left(\beta_j, Var(\hat{\beta_j})\right),$$ em que `\(Var(\hat{\beta_j})\)` é estimada através do correspondente termo da diagonal da matriz de covariâncias `$$\widehat{Var}(\widehat{\boldsymbol{\beta}}) = \hat{\phi} \; (X' \hat{W} X)^{-1}.$$` Vamos denotar por `\(\widehat{EP}(\hat{\beta}_j) = \sqrt{\widehat{Var}(\hat{\beta_j})}\)`, o erro padrão de `\(\hat{\beta_j}\)`. -- - O teste de Wald baseia-se na distribuição assintótica normal dos estimadores de máxima verossimilhança dos parâmetros do modelo. -- >A estatística e o teste de Wald são apresentados no próprio summary de um MLG. -- - A função *waldtest*, do pacote **lmtest** permite aplicar o teste de Wald para hipóteses envolvendo `\(p > 2\)` parâmetros, baseado numa distribuição assintótica `\(\chi^2_p\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Seleção de variáveis - I Em modelos de regressão normal, estudamos diversas técnicas para seleção de variáveis, entre elas: Cp de Mallows, Coeficiente de determinação múltiplo, `\(R^2\)` ajustado, quadrado médio residual que envolvem a análise de todas as possíveis regressões. -- - Técnicas iterativas para seleção de variáveis também foram vistas: forward, backward e stepwise. Esses métodos iterativos podem ser facilmente estendidos para os MLGs. -- **Método do Akaike** O método proposto por Akaike (1974) basicamente se diferencia dos procedimentos anteriores por ser um processo de minimização que não envolve testes estatísticos. -- - A ideia básica é selecionarmos um modelo que seja parcimonioso, ou em outras palavras, que esteja bem ajustado e tenha um número reduzido de parâmetros. -- - O critério consiste em encontrarmos o modelo tal que a quantidade abaixo seja minimizada: $$AIC = D^*(\boldsymbol{y}; \hat{\boldsymbol{\mu}}) + 2p, $$ em que `\(D^*(\boldsymbol{y}; \hat{\boldsymbol{\mu}})\)` denota o desvio do modelo e p o número de parâmetros. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Seleção de variáveis - II - Um modelo parcimonioso concilia bom ajuste (elevado valor para `\(l({\hat{\boldsymbol{\mu}}} ;y)\)` e baixa complexidade (reduzido número de parâmetros, `\(p\)`). Sob esta ótica, quanto menor o valor do AIC, melhor. -- - Uma medida de informação alternativa é o critério de informação de Schwarz (ou critério de informação Bayesiano), definido por: `$$BIC = D^*(\boldsymbol{y}; \hat{\boldsymbol{\mu}}) + p \log{n},$$` que penaliza mais fortemente o modelo por sua complexidade, em relação ao AIC. -- >No R: o AIC e o BIC podem ser extraídos usando AIC `\((\cdot)\)` (ou diretamente do summary) e BIC `\((\cdot)\)`. --- class: animated, hide-logo, bounceInDown ## Política de proteção aos direitos autorais > <span style="color:grey">O conteúdo disponível consiste em material protegido pela legislação brasileira, sendo certo que, por ser o detentor dos direitos sobre o conteúdo disponível na plataforma, o **LECON** e o **NEAEST** detém direito exclusivo de usar, fruir e dispor de sua obra, conforme Artigo 5<sup>o</sup>, inciso XXVII, da Constituição Federal e os Artigos 7<sup>o</sup> e 28<sup>o</sup>, da Lei 9.610/98. A divulgação e/ou veiculação do conteúdo em sites diferentes à plataforma e sem a devida autorização do **LECON** e o **NEAEST**, pode configurar violação de direito autoral, nos termos da Lei 9.610/98, inclusive podendo caracterizar conduta criminosa, conforme Artigo 184<sup>o</sup>, §1<sup>o</sup> a 3<sup>o</sup>, do Código Penal. É considerada como contrafação a reprodução não autorizada, integral ou parcial, de todo e qualquer conteúdo disponível na plataforma.</span> .pull-left[ <img src="images/logo_lecon.png" width="50%" style="display: block; margin: auto;" /> ] .pull-right[ <img src="images/logo_neaest.png" width="50%" style="display: block; margin: auto;" /> ] .center[ [https://lecon.ufes.br](https://lecon.ufes.br/) ] <font size="2"><span style="color:grey">Material elaborado pela equipe LECON/NEAEST: Alessandro J. Q. Sarnaglia, Bartolomeu Zamprogno, Fabio A. Fajardo, Luciana G. de Godoi e Nátaly A. Jiménez.</span></font>