Modelos Lineares Generalizados, Notas de estudo de Estatística

Baixe Modelos Lineares Generalizados e outras Notas de estudo em PDF para Estatística, somente na Docsity! MODELOS DE REGRESSÃO com apoio computacional Gilberto A. Paula Instituto de Matemática e Estat́ıstica Universidade de São Paulo e-mail:giapaula@ime.usp.br home-page:www.ime.usp.br/∼giapaula ii Sumário Prefácio iii 1 Modelos Lineares Generalizados 1 1.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Definição . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 Casos particulares . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Ligações canônicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3.1 Outras ligações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Função desvio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4.1 Análise do desvio . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.5 Função escore e matriz de informação . . . . . . . . . . . . . . . . . . . . . 15 1.6 Estimação dos parâmetros . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.6.1 Estimação de β . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.6.2 Estimação do parâmetro de dispersão . . . . . . . . . . . . . . . . 18 1.7 Teste de hipóteses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 1.7.1 Hipóteses simples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 1.7.2 Modelos encaixados . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 1.7.3 Modelo de análise de variância . . . . . . . . . . . . . . . . . . . . . 26 1.7.4 Regressão linear simples . . . . . . . . . . . . . . . . . . . . . . . . 27 1.7.5 Hipóteses restritas . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.8 Técnicas de diagnóstico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.8.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.8.2 Pontos de alavanca . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 1.8.3 Reśıduos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 1.8.4 Influência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 1.8.5 Influência local . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 v vi 1.8.6 Gráfico da variável adicionada . . . . . . . . . . . . . . . . . . . . . 43 1.8.7 Seleção de modelos . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 1.8.8 Técnicas gráficas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 1.8.9 Bandas de confiança . . . . . . . . . . . . . . . . . . . . . . . . . . 47 1.9 Extensão para os MLGs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 1.9.1 Pontos de alavanca . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 1.9.2 Reśıduos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 1.9.3 Influência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 1.9.4 Influência local . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 1.9.5 Gráfico da variável adicionada . . . . . . . . . . . . . . . . . . . . . 53 1.9.6 Seleção de modelos . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 1.9.7 Técnicas gráficas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 1.9.8 Bandas de confiança . . . . . . . . . . . . . . . . . . . . . . . . . . 55 1.10 Aplicações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 1.10.1 Estudo entre escolaridade e renda . . . . . . . . . . . . . . . . . . . 56 1.10.2 Estudo comparativo de processo infeccioso pulmonar . . . . . . . . 60 1.10.3 Sobrevivência de bactérias . . . . . . . . . . . . . . . . . . . . . . . 61 1.10.4 Estudo seriado com ratos . . . . . . . . . . . . . . . . . . . . . . . . 64 1.10.5 Comparação de cinco tipos de turbina de avião . . . . . . . . . . . 66 1.10.6 Consumo de combust́ıvel . . . . . . . . . . . . . . . . . . . . . . . . 71 1.11 Exerćıcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 2 Modelos para Dados Binários 85 2.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 2.2 Métodos clássicos: uma única tabela 2 × 2 . . . . . . . . . . . . . . . . . . 85 2.2.1 Risco relativo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 2.2.2 Modelo probabiĺıstico não-condicional . . . . . . . . . . . . . . . . . 87 2.2.3 Modelo probabiĺıstico condicional . . . . . . . . . . . . . . . . . . . 88 2.2.4 Teste de hipóteses e estimação intervalar . . . . . . . . . . . . . . . 91 2.3 Métodos clássicos: k tabelas 2 × 2 . . . . . . . . . . . . . . . . . . . . . . . 94 2.3.1 Estimação da razão de chances comum . . . . . . . . . . . . . . . . 94 2.3.2 Testes de homogeneidade . . . . . . . . . . . . . . . . . . . . . . . . 95 2.4 Métodos clássicos: tabelas 2 × k . . . . . . . . . . . . . . . . . . . . . . . . 96 2.5 Aplicações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 vii 2.5.1 Influência do fungicida Avadex no desenvolvimento de tumor em ratos 98 2.5.2 Efeito de um tipo de extrato vegetal na morte de embriões . . . . . 100 2.6 Regressão loǵıstica linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 2.6.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 2.6.2 Regressão loǵıstica simples . . . . . . . . . . . . . . . . . . . . . . . 101 2.6.3 Regressão loǵıstica múltipla . . . . . . . . . . . . . . . . . . . . . . 104 2.6.4 Amostragem retrospectiva . . . . . . . . . . . . . . . . . . . . . . . 105 2.6.5 Seleção de modelos . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 2.6.6 Técnicas de diagnóstico e qualidade do ajuste . . . . . . . . . . . . 113 2.6.7 Modelos de dose-resposta . . . . . . . . . . . . . . . . . . . . . . . . 118 2.6.8 Modelos de dose-resposta de retas paralelas . . . . . . . . . . . . . 127 2.6.9 Superdispersão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 2.6.10 Modelo loǵıstico condicional . . . . . . . . . . . . . . . . . . . . . . 137 2.7 Exerćıcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 3 Modelos para Dados de Contagem 153 3.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 3.1.1 Métodos clássicos: uma única tabela 2 × 2 . . . . . . . . . . . . . . 154 3.1.2 Estratificação : k tabelas 2 × 2 . . . . . . . . . . . . . . . . . . . . 158 3.2 Modelos de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 3.2.1 Propriedades da Poisson . . . . . . . . . . . . . . . . . . . . . . . . 159 3.2.2 Modelos log-lineares . . . . . . . . . . . . . . . . . . . . . . . . . . 160 3.2.3 Relação com a exponencial . . . . . . . . . . . . . . . . . . . . . . . 161 3.2.4 Aplicação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162 3.2.5 Modelo log-linear geral . . . . . . . . . . . . . . . . . . . . . . . . . 164 3.2.6 Superdispersão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 3.3 Relação entre a multinomial e a Poisson . . . . . . . . . . . . . . . . . . . 180 3.3.1 Modelos log-lineares hierárquicos . . . . . . . . . . . . . . . . . . . 182 3.3.2 Exemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 3.4 Exerćıcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188 4 Modelos de Quase-Verossimilhança 195 4.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195 4.2 Respostas independentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197 2 Caṕıtulo 1 √ y no sentido de buscarmos a normalidade dos dados, podemos supor que a distribuição de Y é Poisson e que a relação funcional entre a média de Y e o preditor linear é dada por logµ = η. Essa relação funcional é conveniente, uma vez que garante para quaisquer valores dos parâmetros do preditor linear um valor positivo para µ. Similarmente, para proporções, pode-se pensar na distribuição binomial para a resposta e numa relação fun- cional do tipo log{µ/(1 − µ)}, em que µ é a proporção esperada de sucessos. Nelder e Wedderburn propuseram também um processo iterativo para a estimação dos parâmetros e introduziram o conceito de desvio que tem sido largamente utilizado na avaliação da qualidade do ajuste dos MLGs, bem como no desenvolvimento de reśıduos e medidas de diagnóstico. Inúmeros trabalhos relacionados com modelos lineares generalizados foram publica- dos desde 1972. Um aplicativo, GLIM (Generalized Linear Interactive Models) (vide Aitkin et al., 1989), foi desenvolvido para o ajuste dos MLGs e hoje outros aplicativos, tais como o S-Plus (http://www.insightful.com), R (http://www.r-project.org), SAS(http://www.sas.com), STATA (http://www.stata.com), SUDAAN (http://www.rti. org/sudaan) dentre outros apresentam procedimentos para o ajuste dos MLGs. Uma das extensões mais importantes dos MLGs foi apresentada por Wedderburn (1974), os mod- elos de quase-verossimilhan,̧ que estendem a idéia dos MLGs para situações mais gerais incluindo dados correlacionaods. Os modelos de dispersão (Jørgensen, 1983) ampliam o leque de opções para a distribuição da variável resposta. Liang e Zeger (1986) estendem os modelos de quase-verossimilhança propondo as equações de estimação generalizadas (EEGs) que permitem o estudo de variáveis aleatórias correlacionadas não-Gaussianas. Os modelos não-lineares de famı́lia exponencial (Cordeiro e Paula, 1989a e Wei, 1998) ad- mitem preditor não-linear nos parâmetros. Temos ainda os modelos aditivos generalizados (Hastie e Tibshirani, 1990) que supõem preditor linear formado também por funções semi- paramétricas e os modelos lineares generalizados mistos (Breslow e Clayton, 1993) que admitem a inclusão de efeitos aleatórios Gaussianos no preditor linear. Recentemente, Lee e Nelder (1996, 2001) estenderam o trabalho de Breslow e Clayton propondo modelos lineares generalizados hierárquicos em que o preditor linear pode ser formado por efeitos fixos e efeitos aleatórios não-Gaussianos. Muitos desses resultados são discutidos no livro de McCulloch e Searle (2001). Outras aplicações da estrutura dos MLGs podem ser en- contradas em diversos artigos e livros da literatura Estat́ıstica. Referências de texto no assunto são os livros de McCullagh e Nelder (1989) e Cordeiro (1986). Modelos Lineares Generalizados 3 1.2 Definição Suponha Y1, . . . , Yn variáveis aleatórias independentes, cada uma com densidade na forma dada abaixo f(y; θi, φ) = exp[φ{yθi − b(θi)} + c(y, φ)], (1.1) em que E(Yi) = µi = b ′(θi), Var(Yi) = φ −1Vi, V = dµ/dθ é a função de variância e φ −1 > 0 é o parâmetro de dispersão. A função de variância desempenha um papel importante na famı́lia exponencial, uma vez que a mesma caracteriza a distribuição. Isto é, dada a função de variância, tem-se uma classe de distribuições correspondentes, e vice-versa. Essa propriedade permite a comparação de distribuições através de testes simples para a função de variância. Para ilustrar, a função de variância definida por V (µ) = µ(1 − µ), caracteriza a classe de distribuições binomiais com probabilidades de sucesso µ ou 1 − µ. Uma propriedade interessante envolvendo a distribuição de Y e a função de variância é a seguinte: √ φ(Y − µ) →d N(0, V (µ)), quando φ→ ∞. Ou seja, para φ grande Y segue distribuição aproximadamente normal de média µ e variância φ−1V (µ). Esse tipo de abordagem assintótica, diferente da usual em que n é grande, foi introduzida por Jørgensen (1987). Os modelos lineares generalizados são definidos por (1.1) e pela componente sis- temática g(µi) = ηi, (1.2) em que ηi = x T i β é o preditor linear, β = (β1, . . . , βp) T , p < n, é um vetor de parâmetros desconhecidos a serem estimados, xi = (xi1, . . . , xip) T representa os valores de p variáveis explicativas e g(·) é uma função monótona e diferenciável, denominada função de ligação. Apresentamos a seguir as distribuições mais conhecidas pertencentes à famı́lia exponencial. 1.2.1 Casos particulares Normal Seja Y uma variável aleatória com distribuição normal de média µ e variância σ2, Y ∼ N(µ, σ2). A densidade de Y é expressa na forma 1 σ √ 2π exp{− 1 2σ2 (y − µ)2} = exp[{ 1 σ2 (µy − µ 2 2 ) − 1 2 {log2πσ2 + y 2 σ2 }], 4 Caṕıtulo 1 em que −∞ < µ, y < ∞ e σ2 > 0. Logo, para θ = µ, b(θ) = θ2/2, φ = σ−2 e c(y, φ) = 1 2 logφ/2π − φy2 2 tem-se (1.1). Verifica-se facilmente que a função de variância é dada por V (µ) = 1. Poisson No caso de Y ∼ P (µ), a densidade fica dada por e−µµy/y! = exp{ylogµ− µ− logy!}, em que µ > 0 e y = 0, 1, . . .. Fazendo logµ = θ, b(θ) = eθ, φ = 1 e c(y, φ) = −logy! tem-se (1.1). Segue portanto que V (µ) = µ. Binomial Seja Y ∗ a proporção de sucessos em n ensaios independentes, cada um com probabil- idade de ocorrência µ. Assumiremos que nY ∗ ∼ B(n, µ). A densidade de Y ∗ fica então expressa na forma ( n ny∗ ) µny ∗ (1 − µ)n−ny∗ = exp { log ( n ny∗ ) + ny∗log ( µ 1 − µ ) + nlog(1 − µ) } , em que 0 < µ, y∗ < 1. Obtém-se (1.1) fazendo φ = n, θ = log{µ/(1−µ)}, b(θ) = log(1+eθ) e c(y∗, φ) = log ( φ φy∗ ) . A função de variância aqui fica dada por V (µ) = µ(1 − µ). Gama Seja Y uma variável aleatória com distribuição gama de média µ e coeficiente de variação φ−1/2, denotaremos Y ∼ G(µ, φ). A densidade de Y é dada por 1 Γ(φ) ( φy µ )φ exp ( −φy µ ) d(logy) = exp [ φ { −y µ + log ( 1 µ )} − logΓ(φ) + φlog(φy)− logy ] , em que y ≥ 0, φ > 0, µ > 0 e Γ(φ) = ∫∞0 tφ−1e−tdt é a função gama. Logo, fazendo θ = −1/µ, b(θ) = −log(−θ) e c(y, φ) = (φ− 1)logy + φlogφ− logΓ(φ) tem-se (1.1). Para 0 < φ < 1 a densidade da gama tem uma pole na origem e decresce monotonicamente quando y → ∞. A exponencial é um caso especial quando φ = 1. Para φ > 1 a densidade assume zero na origem, tem um máximo em y = µ − µ/φ e depois decresce para y → ∞. A χ2k é um outro caso especial quando φ = k/2 e µ = k. A distribuição normal é obtida fazendo φ→ ∞. Isto é, quando φ é grande Y ∼ N(µ, φ−1V (µ)). Note que φ = E2(Y )/Var(Y ) é o inverso do coeficiente de variação de Y ao quadrado (φ = 1/(CV )2). A função de variância da gama é dada por V (µ) = µ2. Modelos Lineares Generalizados 7 em que −∞ < y < ∞. Dáı segue que a função de distribuição acumulada fica expressa na forma F (y) = ey/(1 + ey). O modelo loǵıstico binomial é obtido substituindo F (y) por µ e y por η na expressão acima. Como no caso binomial o parâmetro de interesse sempre é uma probabilidade, fica muito razoável que funções de distribuições acumuladas sejam utilizadas para gerarem novas ligações e consequentemente novos modelos. Na Figura 1.1 apresentamos a F (y) da distribuição loǵıstica e da distribuição do valor extremo para valores de y variando no intervalo [−3 , 3]. Note que, a curva loǵıstica é simétrica em torno de F (y) = 1/2, enquanto que a curva do valor extremo apresenta comportamentos distintos para F (y) ≤ 1/2 e F (y) > 1/2. y F( y) -3 -2 -1 0 1 2 3 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Logistica V.Extremo Figura 1.1: Função de distribuição acumulada das curvas loǵıstica e valor extremo. Ligação de Box-Cox Uma classe importante de ligações, pelo menos para observações positivas, são as ligações de Box-Cox, definidas por η = (µλ − 1)/λ, para λ 6= 0 e η = logµ para λ → 0. Note que a idéia agora é aplicar a transformação de 8 Caṕıtulo 1 0 2 4 6 8 10 0 10 20 30 λ = 0.5 λ = 0.6 λ = 0.8 µ η Figura 1.2: Transformação de Box-Cox para alguns valores de λ. Box-Cox, definida no Caṕıtulo 1, na média da variável resposta ao invés de transformar a própria variável resposta. Temos na Figura 1.2 o comportamento de µ para alguns valores de λ e para η variando no intervalo [0 , 10]. -3 -2 -1 0 1 2 3 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 α = 0.5 α = 1.0 α = 2.0 µ η Figura 1.3: Transformação de Aranda-Ordaz para alguns valores de α. Modelos Lineares Generalizados 9 Ligação de Aranda-Ordaz Uma outra transformação importante foi proposta por Aranda-Ordaz (1981) para dados binários. A transformação é dada por η = log { (1 − µ)−α − 1 α } , em que 0 < µ < 1 e α é uma constante desconhecida. Quando α = 1 tem-se a ligação logito η = log{µ/(1− µ)}. Quando α→ 0 tem-se {(1− µ)−α − 1}/α → log(1− µ)−1 de modo que η = log{−log(1− µ)}, obtendo-se portanto a ligação complemento log-log. Na Figura 1.3 temos o comportamento de µ para alguns valores de α. Em muitas situações práticas o interesse pode ser testar se o modelo loǵıstico é apropriado, H0 : α = 1, contra a necessidade de uma transformação na ligação, H1 : α 6= 1. Os MLGs são ajustados no aplicativo S-Plus através do comando glm. Para ilustrar uma aplicação, suponha que temos interesse em ajustar um modelo de Poisson com ligação canônica e que a variável resposta é denotada por resp com variáveis explicativas cov1 e cov2. Podemos mandar os resultados do ajuste para um arquivo (objeto no S-Plus), por exemplo com nome fit.poisson, através do comando fit.poisson < − glm( resp ∼ cov1 + cov2, family=poisson) Com o comando summary(fit.poisson) podemos obter um resumo dos resultados do ajuste. 1.4 Função desvio Sem perda de generalidade, suponha que o logaritmo da função de verossimilhança seja agora definido por L(µ;y) = n ∑ i=1 L(µi; yi), em que µi = g −1(ηi) e ηi = x T i β. Para o modelo saturado (p = n) a função L(µ;y) é estimada por L(y;y) = n ∑ i=1 L(yi; yi). Ou seja, a estimativa de máxima verossimilhança de µi fica nesse caso dada por µ̂ 0 i = yi. Quando p < n, denotaremos a estimativa de L(µ;y) por L(µ̂;y). Aqui, a estimativa de máxima verossimilhança de µi será dada por µ̂i = g −1(η̂i), em que η̂i = x T i β̂. 12 Caṕıtulo 1 acima diz que ∑n i=1(yi − µ̂i)2/σ2 ∼ χ2n−p quando σ2 → 0. No caso do modelo gama, o desvio estará bem aproximado por uma qui-quadrado com n − p graus de liberdade a medida que o coeficiente de variação ficar próximo de zero. 1.4.1 Análise do desvio Suponha para o vetor de parâmetros β a partição β = (βT1 ,β T 2 ) T , em que β1 é um vetor q-dimensional enquanto β2 tem dimensão p − q e φ é conhecido (ou fixo). Portanto, podemos estar interessados em testar as hipóteses H0 : β1 = 0 contra H1 : β1 6= 0. As funções desvio correspondentes aos modelos sob H0 e H1 serão denotadas por D(y; µ̂ 0) e D(y; µ̂), respectivamente, em que µ̂0 é a estimativa de máxima verossimilhança sob H0. A estat́ıstica da razão de verossimilhanças fica nesse caso dada por ξRV = φ{D(y; µ̂0) −D(y; µ̂)}, (1.3) isto é, a diferença entre dois desvios. Como é conhecido, sob a hipótese nula, ξRV ∼ χ2q quando n→ ∞. De forma similar, podemos definir a estat́ıstica F = {D(y; µ̂0) −D(y; µ̂)}/q D(y; µ̂)/(n− p) , (1.4) cuja distribuição nula assintótica é uma Fq,(n−p) quando o denominador de (1.4) é uma estimativa consistente de φ−1 (vide Jørgensen, 1987). A vantagem de utilizar (1.4) em relação a (1.3) é que a estat́ıstica F não depende do parâmetro de dispersão. O resultado (1.4) também é verificado quando φ → ∞ e n é arbitrário. Quando φ é desconhecido a estat́ıstica da razão de verossimilhanças assume uma expressão diferente de (1.3). A estat́ıstica F acima fica, no caso normal linear, reduzida à forma conhecida dada abaixo F = (qs2)−1{ n ∑ i=1 (yi − µ̂0i )2 − n ∑ i=1 (yi − µ̂i)2}, em que s2 = ∑n i=1(yi−µ̂i)2/(n−p) é o erro quadrático médio do modelo com p parâmetros. A forma da estat́ıstica F dada em (1.4) pode ser obtida, em particular, quando testamos uma hipótese de igualdades lineares num modelo de regressão normal linear. Para ilustrar, suponha o modelo y = Xβ + Wγ + , em que  ∼ N(0, σ2I), X é uma matriz n × p, W é aqui uma matriz n × q, ambas de posto completo, β = (β1, . . . , βp) T e γ = (γ1, . . . , γq) T . Vamos supor as hipóteses H0 : Cθ = 0 contra H1 : Cθ 6= 0, Modelos Lineares Generalizados 13 em que θ = (βT ,γT )T e C é uma matriz k × (p+ q) de posto completo. O acréscimo na soma de quadrados de reśıduos devido às restrições em H0 é dado por ASQ(Cθ = 0) = (Cθ̂)T{C(ZTZ)−1CT}−1(Cθ̂), em que θ̂ = (ZTZ)−1ZTy e Z = (X,W). A estat́ıstica F para testar H0 fica então dada por F = ASQ(Cθ = 0)/k D(y; µ̂)/(n− p− q) , em que D(y; µ̂) é o desvio do modelo completo com p+ q parâmetros e ASQ(Cθ = 0) = D(y; µ̂0) −D(y; µ̂), com D(y; µ̂0) sendo o desvio do modelo sob H0. Portanto, F toma a forma F = {D(y; µ̂0) −D(y; µ̂)}/k D(y; µ̂)/(n− p− q) , e segue, sob H0, uma distribuição Fk,(n−p−q). No caso de testarmos H0 : γ = 0 contra H1 : γ 6= 0 a matriz C tem dimensão q × (p + q) com a i-ésima linha tendo o valor 1 na posição p+ i e zeros nas demais posições. Essa formulação pode também ser aplicada quando testamos a inclusão de novas covariáveis num modelo de regressão normal linear. Tabela 1.2 Análise do desvio (ANODEV) supondo dois fatores na parte sistemática. Modelo Desvio Diferença G.L. Testando Constante D0 D0 −DA n(A) − 1 A ignorando B D0 −DB n(B) − 1 B ignorando A +A DA DA −DA+B n(B) − 1 B|A ignorando AB +B DB DB −DA+B n(A) − 1 A|B ignorando AB +A+B DA+B DA+B −DAB {n(A) − 1}× AB|A + B {n(B) − 1} +A+B+AB DAB Para ilustrar o uso das diferenças de desvios para testar hipóteses em modelos en- caixados, suponha um MLG com dois fatores, A e B. O fator A com n(A) ńıveis e o fator B com n(B) ńıveis. Descrevemos na Tabela 1.2 os posśıveis testes envolvendo os dois fatores. Note que, se o interesse é testar a inclusão do fator B dado que o fator A já está no modelo, devemos comparar a diferença φ{D(y; µ̂A) − D(y; µ̂A+B)} com os 14 Caṕıtulo 1 ńıveis cŕıticos da distribuição qui-quadrado com {n(B) − 1} graus de liberdade. Alter- nativamente, podemos comparar o valor observado da estat́ıstica F correspondente com os ńıveis da distribuição F com {n(B) − 1} e {n− n(A) − n(B) + 1} graus de liberdade. No caso normal linear a tabela ANOVA é constrúıda utilizando-se a estat́ıstica F no lugar da diferença entre desvios. A vantagem disso é o fato do parâmetro de dispersão φ−1 não precisar ser estimado. Através do comando anova() o S-Plus fornece uma tabela ANODEV para os ajustes colocados como objetos. Por exemplo, suponha que os objetos fit1.reg, fit2.reg e fit3.reg correspondam aos ajustes de um MLG com um, dois e três fatores, respectivamente. Então, o comando anova(fit1.reg,fit2.reg,fit3.reg) fornece uma tabela ANODEV comparando os três fatores. Tabela 1.3 Análise do desvio referente ao exemplo sobre processo infeccioso pulmonar. Modelo Desvio Diferença G.L. Testando Constante 236,34 - - - + SEXO 235,20 1,14 1 SEXO + IDADE 188,22 46,98 1 IDADE | SEXO + HL 162,55 25,67 3 HL | SEXO + IDADE + FF 157,40 5,15 3 FF | SEXO + IDADE + HL Como aplicação do ANODEV, vamos considerar o exemplo descrito na Seção 1.10.2 em que um modelo loǵıstico linear é ajustado para explicar a ocorrência ou não de câncer de pulmão em pacientes com processo infeccioso pulmonar. A parte sistemática do modelo é representada abaixo 1 + SEXO + IDADE + HL + FF, em que 1 denota a presença de intercepto no modelo, SEXO (1:feminino, 0:masculino), IDADE (em anos) e HL e FF são dois fatores com 4 ńıveis cada um representando a intensidade de dois tipos de célula. Na Tabela 1.3 resumimos alguns resultados. Para calcular os ńıveis descritivos das diferenças apresentadas na Tabela 1.3, usamos o comando pchisq(dv,q) do S-Plus. Por exemplo, para calcular o ńıvel descritivo referente ao efeito do fator SEXO, fazemos Modelos Lineares Generalizados 17 Gama Para o caso gama V (µ) = µ2. Logo, ω = µ2(dθ/dη)2. Em particular, para um modelo log-linear (logµ = η), temos dµ/dη = µ, o que implica em ω = 1. Assim, U(β) = φXTV−1/2(y−µ) e K(β) = φXTX, similarmente ao caso normal. Para ligação canônica, ω = µ2. Normal inversa Nesse caso a função de variância é dada por V (µ) = µ3. Assim, ω = µ3(dθ/dη)2. Pode ser muito razoável aplicar aqui um modelo log-linear, uma vez que as respostas são sempre positivas. Portanto, como ocorre nos modelos log-lineares com resposta de Poisson, os pesos seriam as próprias médias, isto é ω = µ. Em particular para ligação canônica, ω = µ3. 1.6 Estimação dos parâmetros 1.6.1 Estimação de β O processo iterativo de Newton-Raphson para a obtenção da estimativa de máxima verossimilhança de β é definido expandindo-se a função escore U(β) em torno de um valor inicial β(0), tal que U(β) ∼= U(β(0)) + U′(β(0))(β − β(0)), em que U′(β) denota a primeira derivada de U(β) com respeito a β. Assim, repetindo-se o procedimento acima, chega-se ao processo iterativo β(m+1) = β(m) + {−U′(β(m))}−1U(β(m)), m = 0, 1, . . .. Como a matriz −U′(β) pode não ser positiva definida, a aplicação do método de scoring de Fisher substituindo a matriz −U′(β) pelo correspondente valor esperado, pode ser mais conveniente. Isso resulta no seguinte processo iterativo: β(m+1) = β(m) + K−1(β(m))U(β(m)), m = 0, . . .. Se trabalharmos um pouco o lado direito da expressão acima, chegaremos a um processo iterativo de mı́nimos quadrados reponderados β(m+1) = (XTW(m)X)−1XTW(m)z(m), (1.5) 18 Caṕıtulo 1 m = 0, 1, . . ., em que z = η + W−1/2V−1/2(y − µ). Note que z desempenha o papel de uma variável dependente modificada, enquanto W é uma matriz de pesos que muda a cada passo do processo iterativo. A convergência de (1.5) ocorre em um número finito de passos, independente dos valores iniciais utilizados. É usual iniciar (1.5) com η(0) = g(y). Apenas para ilustrar, note que para o caso loǵıstico binomial, tem-se ω = nµ(1 − µ) e variável dependente modificada dada por z = η + (y − nµ)/nµ(1 − µ). Lembrando, para o modelo normal linear tradicional não é preciso recorrer ao processo iterativo (1.5) para a obtenção da estimativa de máxima verossimilhança. Nesse caso, β̂ assume a forma fechada β̂ = (XTX)−1XTy. Tem-se, sob condições gerais de regularidade (vide, por exemplo, Sen e Singer, 1993, Cap. 7), que β̂ é um estimador consistente e eficiente de β e que √ n(β̂ − β) →d Np(0, φ−1Σ−1(β)), conforme n→ ∞, em que Σ(β) = lim n→∞ K(β) n , sendo Σ(β) uma matriz positiva definida e K(β) não contém aqui o multiplicador φ. A demonstração da existência de Σ(β) nem sempre é simples, sendo necessário muitas vezes recorrer a condições suficientes que impliquem na existência de Σ(β). Para ilustrar um caso, vamos supor um MLG com respostas Yij, i = 1, . . . , g e j = 1, . . . , ni, tais que E(Yij) = µij e a parte sistemática é dada por g(µij) = x T i β. As condições suficientes para que Σ(β) exista e seja positiva definida são que ni n → ai > 0 quando n → ∞ e que ∑g i=1 xix T i seja de posto completo, em que n = n1 + · · ·+ng. Outra referência importante sobre as propriedades assintóticas dos estimadores de máxima verossimilhança dos MLGs é Fahrmeir e Kaufmann (1985). Mostra-se também sob certas condições de regularidade que √ n(φ̂− φ) →d N(0, σ2φ), conforme n→ ∞, em que σ2φ = limn→∞−n{ ∑n i=1 c”(yi, φ)}−1. Portanto, um estimador consistente para Var(φ̂) é dado por {∑ni=1 −c”(yi, φ)}−1. 1.6.2 Estimação do parâmetro de dispersão É interessante observar que os parâmetros β e φ são ortogonais, isto é, E[∂2L(β, φ;y)/∂β∂φ] = 0. Uma consequência desse fato é a independência assintótica entre φ̂ e β̂. Derivando o Modelos Lineares Generalizados 19 logaritmo da função de verossimilhança apenas com respeito ao parâmetro φ e igualando a zero, chega-se à seguinte solução: n ∑ i=1 c′(yi, φ̂) = 1 2 D(y; µ̂) − n ∑ i=1 {yiθ̂0i − b(θ̂0i )}, em que D(y; µ̂) denota o desvio do modelo sob investigação. Verifica-se facilmente que as estimativas de máxima verossimilhança para φ nos casos normal e normal inversa são dadas por φ̂ = n/D(y; µ̂). Para o caso gama, a estimativa de máxima verossimilhança de φ sai da equação 2n{logφ̂− ψ(φ̂)} = D(y; µ̂), em que ψ(φ) = Γ′(φ)/Γ(φ) é a função digama. A equação acima pode ser resolvida di- retamente pelo S-PLus através da library mass (Venables e Ripley, 1999). Para ilustrar suponha que os resultados do ajuste sejam guardados em fit.model. Então, para en- contrar a estimativa de máxima verossimilhança de φ com o respectivo desvio padrão aproximado deve-se usar os comandos library(mass) gamma.shape(fit.model) Cordeiro e McCullagh(1991) propõem uma solução em forma fechada para φ usando a expansão (φ grande) ψ(φ) ∼= logφ− 1/2φ− 1/12φ2, que leva ao seguinte resultado: φ̂ ∼= 1 + (1 + 2D̄/3) 1/2 2D̄ , (1.6) em que D̄ = D(y; µ̂)/n. Um problema com essa estimativa é que a mesma não é con- sistente quanda a suposição de distribuição gama é falsa. Um estimador preferido nesse caso, que é consistente, é baseado na estat́ıstica de Pearson φ̃−1 = n ∑ i=1 {(yi − µ̂i)/µ̂i}2/(n− p). A suposição aqui é que β̂ tem sido consistentemente estimado. O S-Plus solta a estimativa φ̂−1 = D(y; µ̂)/(n− p) que não é consistente para φ. 1.7 Teste de hipóteses 1.7.1 Hipóteses simples Buse (1982) apresenta de uma forma bastante didática a interpretação geométrica dos testes da razão de verossimilhanças, escore e Wald para o caso de hipóteses simples. 22 Caṕıtulo 1 em que L(β) = L(β;y). Se, em particular, estamos interessados num subconjunto β1 q-dimensional, a região assintótica de confiança utilizando as estat́ısticas de Wald e da razão de verossimilhanças ficam, respectivamente, dadas por [β; (β̂1 − β)T V̂ar(β̂1)(β̂1 − β) ≤ φ−1χ2q(1 − α)] e [β; 2{L(β̂) − L(β, β̂2(β))} ≤ χ2q(1 − α)], em que β é aqui q-dimensional e β̂2(β) é a estimativa de máxima verossimilhança de β2 dado β (vide Seber e Wild, 1989). 1.7.2 Modelos encaixados φ conhecido(ou fixo) Suponha novamente a partição β = (βT1 ,β T 2 ) T definida na Seção 1.4.1 e as seguintes hipóteses: H0 : β1 = β 0 1 contra H1 : β1 6= β01. Para esse caso temos ξRV = φ{D(y; µ̂0) −D(y; µ̂)}, em que µ̂0 é a estimativa de máxima verossimilhança do MLG com parte sistemática η = η̂01 + η2, em que η̂ 0 1 = ∑q j=1 xjβ 0 j e η2 = ∑p j=q+1 xjβj . A quantidade η̂ 0 1 desempenha o papel de um offset (parte conhecida no preditor linear), conforme a nomenclatura de modelos lineares generalizados. Para ilustrar a utilização do offset, suponha um modelo de Poisson com ligação log-linear, resposta resp, covariáveis cov1 e cov2 e offset dado por logt0. Para ajustar o modelo e armazenar 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 (vide Exemplo 1.10.4). Como ilustração, suponha um MLG com distribuição normal inversa, ligação canônica e preditor linear dado por η = β1 + β2cov2 + β3cov3 e que o interesse é testar H0 : β2 = b, em que b é uma constante diferente de zero, contra H1 : β2 6= b. Os ajustes correspondentes a H0 e H1 são, respectivamente, dados por fit1.ni < − glm( resp ∼ cov3 + offset(b*cov2), family=inverse.gaussian) fit2.ni < − glm( resp ∼ cov2+cov3, family=inverse.gaussian) Logo, de (1.4), a estat́ıstica F para testar H0 : β2 = b contra H1 : β2 6= b fica dada por Modelos Lineares Generalizados 23 F < − (deviance(fit1.ni) - deviance(fit2.ni))/(deviance(fit2.ni)/(n-3)) Note que o offset desaparece para b = 0. O ajuste, nesse caso, fica simplesmente dado por fit1.ni < − glm( resp ∼ cov3, family=inverse.gaussian) Teste de Wald Para testar H0, a estat́ıstica de Wald fica expressa na forma ξW = [β̂1 − β01]T V̂ar−1(β̂1)[β̂1 − β01], em que β̂1 sai do vetor β̂ = (β̂ T 1 , β̂ T 2 ) T . Usando resultados conhecidos de álgebra linear, mostra-se que a variância assintótica de β̂1 é dada por Var(β̂1) = φ −1[XT1 W 1/2M2W 1/2X1] −1, em que X1 sai da partição X = (X1,X2), sendo portanto n×q, X2 é n×(p−q), M2 = I− H2 e H2 = W 1/2X2(X T 2 WX2) −1XT2 W 1/2 é a matriz de projeção ortogonal de vetores do Rn no subespaço gerado pelas colunas da matriz W1/2X2. Em particular, no caso normal linear, temos as simplificações H2 = X2(X T 2 X2) −1X2 e Var(β̂1) = σ 2[XT1 (I− H2)X1]−1. Teste de escore A função escore pode ser expressa alternativamente na forma U(β) = φ1/2XTW1/2rP , em que rP = φ 1/2V−1/2(y − µ) é conhecido como reśıduo de Pearson. Note que rP tem a mesma distribuição de Y, no entanto, E(rP ) = 0 e Var(rP ) = I. O teste de escore é definido por ξSR = U1(β̂ 0 )T V̂ar0(β̂1)U1(β̂ 0 ), em que U1(β) = ∂L(β;y)/∂β1 = φX T 1 W 1/2V−1/2(y − µ), β̂0 = (β0T1 , β̂ 0T 2 ) T e β̂ 0 2 é a estimativa de máxima verossimilhança de β2 sob o modelo com componente sistemática η = η̂01 +η2, isto é, sob H0, em que η̂ 0 1 = X1β 0 1 e η2 = X2β2. Se trabalharmos um pouco mais a expressão para Var(β̂1), chegaremos ao seguinte: Var(β̂1) = φ −1(RTWR)−1, em que R = X1 − X2C e C = (XT2 WX2)−1XT2 WX1. Aqui C é uma matriz n × q cuja j-ésima coluna é o vetor de coeficientes da regressão linear (com pesos W) da j-ésima coluna de X1 sobre X2. Assim, R pode ser interpretado como sendo uma matriz n× q de 24 Caṕıtulo 1 reśıduos. A j-ésima coluna de R corresponde aos reśıduos ordinários da regressão linear (com pesos W) da j-ésima coluna de X1 sobre X2. Assim, o teste de escore fica reexpresso na forma (vide Cordeiro, Ferrari e Paula, 1993) ξSR = r̂ T P0 Ŵ 1/2 0 X1(R̂ T 0 Ŵ0R̂0) −1XT1 Ŵ 1/2 0 r̂P0, com as quantidades r̂P0, Ŵ0 e R̂0 sendo avaliadas em β̂ 0 . Para ilustrar o cálculo da estat́ıstica de escore, suponha um MLG com preditor linear dado por η = β1 +β2cov2 +β3cov3 +β4cov4 e que o interesse é testar H0 : β3 = β4 = 0. As matrizes X1 e X2 serão então dadas por X1 = [cov3 , cov4] e X2 = [1 , cov2]. Se temos um modelo de Poisson, por exemplo com ligação canônica, então como já vimos ω = µ. Logo, Ŵ0 = diag{µ̂01, . . . , µ̂0n}, em que µ̂01, . . . , µ̂0n são os pesos sob H0, ou seja, os pesos do modelo ajustado de Poisson com preditor linear η = β1 +β2cov2. Portanto, precisamos apenas fazer esse ajuste e dáı computarmos Ŵ0, R̂0, r̂P0 e finalmente ξSR. Chamando no S-Plus os pesos por w, Ŵ0 por W, r̂P0 por rp e R̂0 por R, os passos para o cálculo de ξSR são dados abaixo X1 < − cbind(cov3 , cov4) X2 < − cbind(1 , cov2) fit.poisson < − glm( resp ∼ cov2, family=poisson) rp < − resid(fit.poisson, type=‘‘pearson") w < − fit.poisson$weights W < − diag(w) A < − solve(t(X2)%*%W%*%X2) C1 < − A%*%t(X2)%*%W%*%cov3 C2 < − A%*%t(X2)%*%W%*%cov4 C < − cbind(C1 , C2) R < − X1 - X2%*%C SR < − solve(t(R)%*%W%*%R) SR < − t(rp)%*%sqrt(W)%*%X1%*%SR%*%t(X1)%*%sqrt(W)%*%rp Em particular, para o caso normal linear, C = (XT2 X2) −1XT2 X1 e rP = (y−µ)/σ. Logo, ξSR = σ −2(y−µ̂0)TX1(RTR)−1XT1 (y−µ̂0), em que R = X1−X2(XT2 X2)−1XT2 X1 = (I− H2)X1. Aqui, também as estat́ısticas da razão de verossimilhanças e de Wald coincidem com a estat́ıstica de escore. Isso em geral vale para o modelo normal linear. Modelos Lineares Generalizados 27 Tabela 1.4 Expressões para as estat́ısticas de escore e de Wald. Distribuição ξSR ξW Normal m 2σ2 (ȳ1 − ȳ2)2 m2σ2 β̂2 Poisson m 2ȳ (ȳ1 − ȳ2)2 mȳ1ȳ2(ȳ1+ȳ2) β̂ 2 Binomial 2m y(2m−y) (y1 − y2)2 β̂ 2 m y1(m−y1)y2(m−y2) y1(m−y1)+y2(m−y2) Gama φm 2ȳ2 (ȳ1 − ȳ2)2 φm(ȳ1ȳ2) 2 (ȳ21+ȳ 2 2) β̂2 Normal inversa φm 2ȳ3 (ȳ1 − ȳ2)2 φm(ȳ1ȳ2) 3 (ȳ31+ȳ 3 2) β̂2 1.7.4 Regressão linear simples Suponha agora um MLG com parte sistemática na forma linear simples g(µi) = α + βxi, i = 1, . . . , n, e as hipóteses H0 : β = 0 contra H1 : β 6= 0 com φ conhecido. Nesse caso obtemos Rj = (xj ∑n i=1 ωi− ∑n i=1 ωixi)/ ∑n i=1 ωi e R TWR = ∑n i=1 ωiR 2 i . Consequentemente, R̂0j = xj− x̄ e R̂T0 Ŵ0R̂0 = ω̂0 ∑n i=1(xi − x̄)2. Aqui, também temos µ̂0 = ȳ. A estat́ıstica de escore fica portanto dada por ξSR = φ V̂0 {∑ni=1 xi(yi − ȳ)}2 ∑n i=1(xi − x̄)2 , (1.9) em que V̂0 = V (ȳ). Similarmente, obtemos para a estat́ıstica de Wald ξW = φβ̂ 2 n ∑ i=1 ω̂iR̂ 2 i , (1.10) em que β̂ é a estimativa de β sob H1. 1.7.5 Hipóteses restritas Pode haver interesse, em algumas situações práticas, em testar hipóteses na forma de igualdades lineares, isto é, H0 : Cβ = 0 contra H1 : Cβ 6= 0, em que C é uma matriz k×p de posto completo. A estimativa de máxima verossimilhança sob a hipótese alternativa coincide com a estimativa de máxima verossimilhança irrestrita β̂, no entanto, obter a 28 Caṕıtulo 1 estimativa de máxima verossimilhança sob H0 pode ser mais complexo, requerendo o uso de algum processo iterativo. Nyquist (1991) propõe um processo iterativo para a obtenção da estimativa de máxima verossimilhança em MLGs com parâmetros restritos na forma Cβ = 0. O processo iterativo é dado abaixo β(m+1)c = β̃ (m+1) − (XTW(m)X)−1CT{C(XTW(m)X)−1CT}−1Cβ̃(m+1), m = 0, 1, . . ., em que β̃ (m+1) é (1.5) avaliado na estimativa restrita β(m)c . A matriz de variância-covariância assintótica de β̂c é dada por Var(β̂c) = φ −1(XTWX)−1[I −CT{C(XTWX)−1CT}−1C(XTWX)−1]. Os testes estat́ısticos tomam formas similares aos testes do caso irrestrito. Em particular, quando φ é conhecido, o teste da razão de verossimilhanças fica dado por ξRV = φ{D(y; µ̂0) −D(y; µ̂)}, em que µ̂0 denota aqui a estimativa de máxima verossimilhança de µ sob H0 : Cβ = 0. Já, o teste de escore, toma a forma ξSR = φ −1U(β̂c) T (XTŴ0X) −1U(β̂c), em que Ŵ0 é aqui avaliado em β̂c. Finalmente, o teste de Wald fica dado por ξW = [Cβ̂ − 0]T [V̂ar(Cβ̂)]−1[Cβ̂ − 0] = φβ̂ T CT [C(XTŴX)−1CT ]−1Cβ̂. Sob H0 e para grandes amostras, as estat́ısticas ξRV , ξW e ξSR seguem uma distribuição χ2k. A distribuição nula assintótica dos testes acima para o caso H0 : Cβ = 0 contra H1 − H0, em que H1 : Cβ ≥ 0, é uma mistura de distribuições do tipo qui-quadrado. Fahrmeir e Klinger (1994) discutem esse tipo de teste em MLGs ( vide também Paula, 1997). 1.8 Técnicas de diagnóstico 1.8.1 Introdução Uma etapa importante na análise de um ajuste de regressão é a verificação de posśıveis afastamentos das suposições feitas para o modelo, especialmente para a parte aleatória Modelos Lineares Generalizados 29 e para a parte sistemática do modelo, bem como a existência de observações extremas com alguma interferência desproporcional nos resultados do ajuste. Tal etapa, conhecida como análise de diagnóstico, tem longa data, e iniciou-se com a análise de reśıduos para detectar a presença de pontos extremos e avaliar a adequação da distribuição proposta para a variável resposta. Uma referência importante nesse tópico é o artigo de Cox e Snell (1968) em que é apresentada uma forma bastante geral de definir reśıduos, usada até os dias atuais. Belsley, Kuh e Welsch (1980) e Cook e Weisberg (1982) discutem a padronização de reśıduos para o caso normal linear. Pregibon (1981) propõe o componente do desvio como reśıduo na classe dos modelos lineares generalizados e sugere uma padronização que é mais tarde comprovada matematicamente por McCullagh (1987) que usa as aproximações propostas por Cox e Snell (1968). Nesse mesmo trabalho McCullagh apresenta uma outra forma de padronização para o componente do desvio em que procura corrigir os efeitos de assimetria e curtose. Atkinson (1981) propõe a construção por simulação de Monte Carlo de uma banda de confiança para os reśıduos da regressão normal linear, a qual denominou envelope, e que permite uma melhor comparação entre os reśıduos e os percentis da dis- tribuição normal padrão. Williams (1984,1987) discute, com base em estudos de simulação de Monte Carlo, a aproximação da forma padronizada proposta por Pregibon (1981) en- contrando fortes evidências de concordância entre a distribuição emṕırica do componente do desvio padronizado e a distribuição normal padrão para vários MLGs. Williams (1987) também discute a construção de envelopes em MLGs. Davison e Gigli (1989) estendem a proposta de Cox e Snell (1968) e definem uma forma geral de padronização para o com- ponente do desvio para distribuições cont́ınuas, mesmo quando a função de distribuição acumulada não é expressa em forma fechada. Fahrmeir e Tutz (1994) estendem o trabalho de McCullagh (1987) para modelos mais gerais, não pertencentes à famı́lia exponencial de distribuições. Paula (1995) apresenta uma forma padronizada para o componente do desvio em MLGs com parâmetros restritos na forma de desigualdades lineares Cβ ≥ 0 e verifica através de estudos de simulação forte concordância na maioria dos modelos estudados entre a distribuição emṕırica do reśıduo padronizado e a distribuição normal padrão, generalizando os resultados de Williams para parâmetros restritos. De Souza e Paula (2002) usam o método proposto por Davison e Gigli (1989) a fim de obterem uma forma padronizada para o componente do desvio em modelos de regressão von Mises, os quais têm sido aplicados na análise de dados circulares. A construção de envelopes com 32 Caṕıtulo 1 A matriz H é simétrica e idempotente e é conhecida como matriz hat, uma vez que faz µ̂ = Hy. Por ser idempotente, tem-se que posto(H) = tr(H) = ∑n i=1 hii = p. O elemento hii = x T i (X TX)−1xi desempenha um papel importante na construção de técnicas de diagnóstico. Mostra-se que 1 n ≤ hii ≤ 1c (vide Cook e Weisberg, 1982), em que c é o número de linhas de X idênticas a xTi . O i-ésimo valor ajustado fica então dado por ŷi = hiiyi + ∑ i6=j hjiyj, (1.11) e pelo fato da matriz H ser idempotente ∑ j 6=i h2ij = hii(1 − hii). Note que hii = 1 implica em ŷi = yi, todavia a rećıproca não é necessariamente verdadeira. Logo, para valores altos de hii predomina na expressão (1.11) a influência de yi sobre o correspondente valor ajustado. Assim, é muito razoável utilizar hii como uma medida da influência da i-ésima observação sobre o próprio valor ajustado. Note também que hii = ∂ŷi/∂yi, ou seja, hii corresponde à variação em ŷi quando yi é acrescido de um infinitésimo. Supondo que todos os pontos exerçam a mesma influência sobre os valores ajustados, pode-se esperar que hii esteja próximo de tr(H) n = p n . Convém então examinar aqueles pontos tais que hii ≥ 2pn , que são conhecidos como pontos de alavanca ou de alto leverage e geralmente estão localizados em regiões remotas no subespaço gerado pelas colunas da matriz X. Esses pontos podem ser também informativos com relação à estimativa β̂. Uma outra maneira de entender hii é construindo a matriz Jacobiana de leverages (vide, por exemplo, St. Laurent e Cook, 1993; Paula, 1999) quando a i-ésima observação é perturbada de modo que o novo valor observado seja dado por yi(b) = yi + b, em que b é uma constante real. O novo vetor de valores ajustados fica dado por ŷ(b) = X(XTX)−1XTy(b), em que y(b) = (y1, . . . , yi−1, yi + b, yi+1, . . . , yn) T . A matriz Jacobiana de leverages é definida por J(b) = lim b→0 1 b {ŷ(b) − ŷ}, e representa a variação no vetor de valores ajustados sob uma variação infinitesimal no i-ésimo valor observado. É fácil verificar que J(b) = X(XTX)−1XT f = Hf , Modelos Lineares Generalizados 33 em que f é um vetor n × 1 de zeros com o valor 1 na i-ésima posição. Portanto, prova- se que hii representa a variação no valor predito da i-ésima observação quando o valor observado é acrescido de um infinitésimo. Para ilustrar como obter os valores hii no S-Plus, suponha um modelo normal linear de variável resposta resp, fatores A e B e covariáveis cov1 e cov2. Supor ainda que os resultados do ajuste serão armazenadas em fit.model. Esse modelo pode ser ajustado de duas formas fit.model < − lm( resp ∼ A + B + cov1 + cov2) ou, alternativamente, como um MLG fit.model < − glm( resp ∼ A + B + cov1 + cov2, family=normal) É claro que a primeira maneira é mais simples. Para gerar a matriz modelo (incluindo a constante) fazemos X < − model.matrix( ∼ A + B + cov1 + cov2) Assim, temos em X a matriz modelo correspondente. O cálculo da matriz de projeção H pode ser feito seguindo os passos descritos abaixo H < − solve(t(X)% ∗ %X) H < − X% ∗ %H% ∗ %t(X) Logo, podemos obter hii extraindo os elementos da diagonal principal de H h < − diag(H) Outras maneiras mais fáceis de extrair os elementos h′iis de uma regressão linear são através dos comandos h < − lm.influence(fit.model)$hat h < − hat(X,T) Para construir um index plot de hii, a fim de detectar pontos de alavanca, fazemos plot(h, xlab=‘‘indice ’’, ylab= ‘‘leverage ’’) É importante que os comandos openlook() ou motif() tenham sido acionados na versão UNIX e win.graph() na versão Windows. 1.8.3 Reśıduos Dos resultados descritos na seção anterior segue que E(r) = (I−H)E(Y) = 0 e Var(r) = σ2(I−H). Isto é, ri tem distribuição normal de média zero e variância Var(ri) = σ2(1−hii). Além disso, a covariância entre ri e rj , i 6= j, fica dada por Cov(ri, rj) = −σ2hij. Como os r′is têm variâncias diferentes, é conveniente expressá-los em forma padronizada 34 Caṕıtulo 1 a fim de permitir uma comparabilidade entre os mesmos. Uma definição natural seria di- vidir ri pelo respectivo desvio padrão, obtendo-se o reśıduo studentizado ti = ri s(1 − hii)1/2 , i = 1, . . . , n, em que s2 = ∑n i=1 r 2 i /(n− p). No entanto, como ri não é independente de s 2, ti não segue uma distribuição t de Student como se poderia esperar. Mostra-se (vide Cook e Weisberg, 1982) que t2i /(n − p) segue uma distribuição beta com parâmetros 1 2 e (n − p − 1)/2. Logo, temos que E(ti) = 0, Var(ti) = 1 e Cov(ti, tj) = −hij/{(1 − hii)(1 − hjj)}1/2, i < j. O problema da dependência entre ri e s 2 pode ser contornado substituindo s2 por s2(i), o erro quadrático médio correspondente ao modelo sem a i-ésima observação. O ı́ndice (i) indica que a i-ésima observação foi exclúıda. Mostra-se usando (1.16) que (n− p)s2 σ2 = (n− p− 1)s2(i) σ2 + r2i σ2(1 − hii) , e dáı segue usando o teorema de Fisher-Cochran (vide, por exemplo, Rao, 1973, p.185) a independência entre s2(i) e r 2 i . Além disso, obtém-se (n− p− 1)s2(i) = n ∑ j=1 r2j − r2i (1 − hii) e dáı segue, após alguma álgebra, que s2(i) = s 2 ( n− p− t2i n− p− 1 ) . (1.12) Assim, fica fácil mostrar que o novo reśıduo studentizado t∗i = ri s(i){1 − hii}1/2 segue uma distribuição central tn−p−1. Se ainda substituirmos (1.12) na expressão acima mostramos que t∗i é uma transformação monótona de ti, t∗i = ti ( n− p− 1 n− p− t2i )1/2 . O reśıduo ti pode ser calculado pela sequência de comandos lms < − summary(fit.model) s < − lms$sigma Modelos Lineares Generalizados 37 que para o caso de p = 2 é um elipsóide no R2 centrado em β̂. Tal medida, conhecida como distância de Cook, é definida por Dδ = (β̂ − β̂δ)T (XTX)(β̂ − β̂δ) ps2 , (1.17) e mede quanto a perturbação δ = (δ1, . . . , δn) T afasta β̂δ de β̂, segundo a métrica M = XTX. Por exemplo, se Dδ > Fp,(n−p)(1−α), significa que a perturbação está deslocando o contorno do elipsóide para um contorno correspondente a um ńıvel de significância menor do que α. Em particular, quando o i-ésimo ponto é exclúıdo, a distância de Cook fica expressa na forma Di = (β̂ − β̂(i))T (XTX)(β̂ − β̂(i)) ps2 = { ri s(1 − hii)1/2 }2 hii (1 − hii) 1 p = t2i hii (1 − hii) 1 p . Portanto, Di será grande quando o i-ésimo ponto for aberrante (ti grande) e/ou quando hii for próximo de um. A distância de Cook pode ser calculada da seguinte maneira: di < − (ti^ 2)*h / (p*(1-h)) A distância Di poderá não ser adequada quando ri for grande e hii for pequeno. Nesse caso, s2 pode ficar inflacionado e não ocorrendo nenhuma compensação por parte de hii, Di pode ficar pequeno. Uma medida supostamente mais apropriada foi proposta por Belsley, Kuh e Welsch (1980), sendo definida por DFFITSi = |ri| s(i)(1 − hii)1/2 { hii (1 − hii) }1/2 = |t∗i | { hii (1 − hii) }1/2 . O DFFITSi é calculado conforme abaixo dfit < − abs(tsi)*(h/(1-h))^ .5 Como o valor esperado de hii é p n , é razoável dar mais atenção àqueles pontos tais que DFFITSi ≥ 2 { p (n− p) }1/2 . 38 Caṕıtulo 1 Aparentemente Di e DFFITSi seriam medidas de influência competitivas, uma vez que DFFITSi parece ser mais adequada para avaliar a influência nas estimativas dos coefi- cientes de um ponto aberrante com hii pequeno. No entanto, como mostram Cook, Peña e Weisberg (1988) Di e DFFITSi medem coisas diferentes. Ambas podem ser expressas a partir da medida mais geral de influência denominada afastamento da verossimilhança (likelihood displacement) proposta por Cook e Weisberg (1982). A medida Di mede essen- cialmente a influência das observações nos parâmetros de locação, enquanto DFFITSi tem o propósito de medir a influência das observações nos parâmetros de locação e escala. Como é pouco provável que um ponto com ri alto e hii pequeno seja influente nas estima- tivas dos coeficientes, o uso de Di não compromete a detecção de observações influentes. Cook, Peña e Weisberg observam também que DFFITSi não é um medida completa de in- fluência nos parâmetros de locação e escala simultaneamente, podendo falhar em algumas situações. Uma medida mais geral nesse caso é proposta pelos autores. Atkinson (1985) propôs uma outra medida de influência que é um aperfeiçoamento do DFFITSi, Ci = { (n− p) p hii (1 − hii) }1/2 |t∗i |. Aqui, quando o experimento for balanceado, isto é, todos os h′iis forem iguais, tem-se Ci = |t∗i |. A vantagem de Ci é que a mesma pode ser utilizada em gráficos normais de probabilidades. Ilustração As Figuras 1.4a-1.4d ilustram as diferenças entre pontos aberrantes, alavanca e influentes. Na Figura 1.4a temos os pontos alinhados sem nenhum tipo de perturbação. Na Figura 1.4b perturbamos o ponto #3 fazendo-o aberrante. Note que a exclusão do mesmo (reta pontilhada) altera apenas o intercepto, isto é, os valores ajustados. É um ponto que não está muito afastado dos demais, logo tem um valor para hii relativamente pequeno. Já na Figura 1.4c, perturbamos o ponto #5 de modo que o mesmo fique mais afastado no subespaço gerado pelas colunas da matriz X. É um ponto de alavanca, todavia a elim- inação do mesmo não muda praticamente nada nas estimativas dos parâmetros. Como é um ponto com hii relativamente alto, as variâncias dos valores ajustados dos pontos próximos ao mesmo serão maiores do que as variâncias dos valores ajustados correspon- dentes aos demais pontos. Finalmente, na Figura 1.4d, perturbamos novamente o ponto #5 fazendo-o agora influente e também alavanca. O mesmo, além de mudar a estimativa Modelos Lineares Generalizados 39 da inclinação da reta ajustada, continua mais afastado do que os demais. As posśıve is situações discutidas acima, quando detectadas num ajuste de regressão, devem ser examinadas cuidadosamente antes de qualquer decisão. Encontrar razões que expliquem o fato dos pontos terem um comportamento at́ıpico com relação aos demais pon- tos pode ajudar a entender melhor a relação entre as variáveis explicativas e o fenômeno sob investigação como também a traçar uma poĺıtica de utilização do modelo ajustado, que não necessariamente implica na eliminação de tais pontos que deve ser o último recurso a ser utilizado. Mudanças na distribuição postulada para a variável resposta, inclusão, eliminação ou mesmo transformação de variáveis explicativas podem ajudar a atenuar a influência de observações. O uso de métodos robustos (vide, por exemplo, Venables e Rip- ley, 1999, Cap.8) ou modelos robustos (vide, por exemplo, Galea, Paula e Uribe-Opazo, 2003) são outras opções a serem tentadas antes da eventual eliminação de pontos. x y 1 2 3 4 5 1 2 3 4 5 (a) x y 1 2 3 4 5 1 2 3 4 5 (b) 3 x y 1 2 3 4 5 6 7 1 2 3 4 5 6 7 (c) 5 x y 1 2 3 4 5 6 7 2 4 6 8 (d) 5 Figura 1.4: Ilustração de pontos aberrantes, influentes e alavanca. 42 Caṕıtulo 1 De forma similar, o autovetor correspondente padronizado e em valor absoluto é obtido com os comandos dmax < − eigen(A)$vec[,1] dmax < − dmax/sqrt(Lmax) dmax < − abs(dmax) Quando o interesse é verificar a influência local das observações num coeficiente partic- ular, Cook (1986) mostra que o autovetor dmax pode ser obtido de forma similar ao caso descrito acima. Esse autovetor contém a influência local das observações na estimativa do coeficiente sob estudo. Assim, particionando a matriz X tal que X = (X1,X2), em que X1 é um vetor n × 1 correspondente à variável explicativa sob estudo e X2 uma matriz n× (p− 1) correspondente às demais variáveis explicativas, o vetor dmax fica dado por dTmax = ( v1r1√ λmax , . . . , vnrn√ λmax ) , em que v1, . . . , vn são os reśıduos ordinários da regressão linear de X1 sobre as colunas de X2, ou seja, o vetor v = (v1, . . . , vn) T é dado por v = (I−H2)X1, H2 = X2(XT2 X2)−1XT2 . Aqui, a matriz A tem posto m = 1. Logo, há apenas um autovalor diferente de zero. Nesse caso, podemos tanto utilizar o procedimento descrito acima para calcular dmax como obtê- lo diretamente sem precisar calcular a matriz H2. Para ilustrar, suponha que os resultados do ajuste estão armazenados em fit.model. Para extrair o vetor r precisamos fazer r < − resid(fit.model) Se o modelo tem as covariáveis cov1 e cov2 além dos fatores A e B, o vetor dmax corre- spondente, por exemplo à covariável cov1, sai de fit < − lm( cov1 ∼ A + B + cov2 - 1) v < − resid(fit) dmax < − v*r tot < − t(dmax)%*%dmax dmax < − dmax/sqrt(tot) dmax < − abs(dmax) Uma outra maneira de interpretação do método de influência local que usa conceitos de curvatura pode ser encontrado em diversos artigos tais como Cook (1986, 1987), Thomas e Cook (1990) e Galea, Paula e Bolfarine (1997). Modelos Lineares Generalizados 43 1.8.6 Gráfico da variável adicionada Suponha novamente o modelo de regressão dado em (1.13), em que ω é agora uma variável adicional qualquer. Definindo Z = (X,ω), mostra-se facilmente que a estimativa de mı́nimos quadrados de θ = (βT ,γ)T é dada por θ̂ = (ZTZ)−1ZTy. Em particular mostra- se, após alguma álgebra, que γ̂ = ωT (I −H)y ωT (I− H)ω = ωT r ωT (I − H)ω . Isto é, γ̂ é o coeficiente da regressão linear passando pela origem do vetor de reśıduos r = (I − H)y sobre o novo reśıduo υ = (I − H)ω. Portanto, um gráfico de r contra υ pode fornecer informações sobre a evidência dessa regressão, indicando quais observações que estão contribuindo para a relação e quais observações que estão se desviando da mesma. Esse gráfico, conhecido como gráfico da variável adicionada, pode revelar quais pontos que estão influenciando (e de que maneira) a inclusão da nova variável no modelo. Para ilustrar a construção do gráfico da variável adicionada, vamos supor novamente o modelo com duas covariáveis e dois fatores. O gráfico da variável adicionada para avaliar a influência das observações no coeficiente de cov1, pode ser constrúıdo com os comandos fit < − lm( resp ∼ cov2 + A + B) r < − resid(fit) fit1 < − lm( cov1 ∼ cov2 + A + B) v < − resid(fit1) plot(v,r, xlab= ‘‘residuo v ’’, ylab= ‘‘residuo r ’’) 1.8.7 Seleção de modelos Existem vários procedimentos para a seleção de modelos de regressão, embora nenhum deles seja consistente, ou seja, mesmo para amostras grandes selecione com probabilidade um as variáveis explicativas com coeficiente de regressão não nulo. Os procedimentos mais conhecidos são maior R2p, menor s 2 p, Cp, forward, backward, stepwise e AIC (vide, por exemplo, Neter et al., 1996, Cap. 8), além de outros métodos que usam computação intensiva. Alguns desses métodos serão descritos brevemente a seguir. 44 Caṕıtulo 1 Método forward Inicia-se o método pelo modelo µ = α. Ajusta-se então para cada variável explicativa o modelo µ = α + βjxj , (j = 1, . . . , q). Testa-se H0 : βj = 0 contra H1 : βj 6= 0. Seja P o menor ńıvel descritivo dentre os q testes. Se P ≤ PE , a variável correspondente entra no modelo. Supor que X1 tenho sido escolhida. Então, no passo seguinte ajusta-se os modelos µ = α+ β1x1 + βjxj , (j = 2, . . . , q). Testa-se H0 : βj = 0 contra H1 : βj 6= 0. Seja P o menor ńıvel descritivo dentre os (q− 1) testes. Se P ≤ PE , a variável correspondente entra no modelo. Repetir o procedimento até que ocorra P > PE. Método backward Inicia-se o procedimento pelo modelo µ = α+ β1x1 + · · ·+ βqxq. Testa-se H0 : βj = 0 contra H1 : βj 6= 0 para j = 1, . . . , q. Seja P o maior ńıvel descritivo dentre os q testes. Se P > PS, a variável correspondente sai do modelo. Supor que X1 tenho sáıdo do modelo. Então, ajusta-se o modelo µ = α+ β2x2 + · · ·+ βqxq. Testa-se H0 : βj = 0 contra H1 : βj 6= 0 para j = 2, . . . , q. Seja P o maior ńıvel descritivo dentre os (q − 1) testes. Se P > PS, então a variável correspondente sai do modelo. Repetir o procedimento até que ocorra P ≤ PS. Método stepwise É uma mistura dos dois procedimentos acima. Inicia-se o processo com o modelo µ = α. Após duas variáveis terem sido inclúıdas no modelo, verifica-se se a primeira não sai do modelo. O processo continua até que nenhuma variável seja inclúıda ou seja retirada do modelo. Geralmente adota-se 0, 15 ≤ PE, PS ≤ 0, 25. Uma sugestão seria usar PE = PS = 0, 20. Modelos Lineares Generalizados 47 de regressão normal linear considerando m = 100. Para rodar o programa é preciso apenas colocar modelo ajustado em fit.model. Dáı, deve-se bater source(‘‘envel.norm ’’) em que envel.norm é o nome do arquivo externo em que deve estar o programa para gerar os envelopes (vide Apêndice). 1.8.9 Bandas de confiança Uma banda de confiança de coeficiente 1−α pode ser constrúıda para µ(z) = zTβ, ∀z ∈ IRp (vide, por exemplo, Casella e Straederman, 1980). Temos que β̂−β ∼ Np(0, σ2(XTX)−1). Logo, uma banda de confiança de coeficiente 1−α para a média µ(z), ∀z ∈ IRp, fica dada por zT β̂ ± σ√cα{zT (XTX)−1z}1/2, ∀z ∈ IRp, em que cα é tal que Pr{χ2p ≤ cα} = 1−α. É importante observar que z é um vetor p× 1 que varia livremente no IRp enquanto X é uma matriz fixa. 1.9 Extensão para os MLGs 1.9.1 Pontos de alavanca A idéia que está por trás do conceito de ponto de alavanca (vide, por exemplo, Hoaglin e Welsch, 1978; Cook e Weisberg, 1982; Emerson, Hoaglin e Kempthorne, 1984; St. Laurent e Cook, 1992 e Wei, Hu e Fung, 1998) é de avaliar a influência de yi sobre o próprio valor ajustado ŷi. Essa influência pode ser bem representada pela derivada ∂ŷi/∂yi que coincide, como foi visto na Seção 1.8.2, com hii no caso normal linear. Recentemente, Wei, Hu e Fung (1998) propuseram uma forma bastante geral para ∂ŷ/∂y quando a resposta é cont́ınua e que pode ser aplicada em diversas situações de estimação. No caso de MLGs a matriz (n× n) ∂ŷ/∂y pode ser obtida da forma geral ∂ŷ ∂y = {Dβ(−L̈ββ)−1L̈βy}|β̂, em que Dβ = ∂µ/∂β, L̈ββ = ∂ 2L(β)/∂β∂βT e L̈βy = ∂ 2L(β)/∂β∂yT . No caso de MLGs com ligação canônica mostra-se facilmente que ∂ŷ ∂y = V̂X(XT V̂X)−1XT . 48 Caṕıtulo 1 Outra definição de ponto de alavanca que tem sido muito utilizada na classe dos MLGs embora não coincida com a expressão acima, exceto no caso de resposta cont́ınua e ligação canônica, é constrúıda fazendo uma analogia entre a solução de máxima verossimilhança para β̂ num MLG e a solução de mı́nimos quadrados de um regressão normal ponderada. Para ver isso, note que na convergência do processo iterativo dado em (1.5), tem-se o seguinte: β̂ = (XTŴX)−1XTŴz, em que z = η̂+Ŵ−1/2V̂−1/2(y−µ̂). Portanto, β̂ pode ser interpretado como a solução de mı́nimos quadrados da regressão linear de Ŵ1/2z contra as colunas de Ŵ1/2X. A matriz de projeção da solução de mińınimos quadrados da regressão linear de z contra X com pesos W fica dada por H = W1/2X(XTWX)−1XTW1/2, que sugere a utilização dos elementos da diagonal principal de Ĥ para detectar-se a presença de pontos de alavanca nesse modelo de regressão normal ponderada. Essa ex- tensão para MLGs foi proposta por Pregibon (1981). Moolgavkar, Lustbader e Venzon (1984) estendem a proposta de Pregibon para modelos não-lineares e sugerem o uso dos elementos da diagonal principal da matriz de projeção no plano tangente à solução de máxima verossimilhança µ(β̂) para avaliar pontos de alavanca. Hosmer e Lemeshow (1989) mostram, contudo, que o uso da diagnonal principal da matriz de projeção H deve ser feito com algum cuidado em regressão loǵıstica e que as interpretações são diferentes daquelas do caso normal linear. 1.9.2 Reśıduos A definição de um reśıduo studentizado para os MLGs pode ser feita analogamente à regressão normal linear como veremos a seguir. Todavia, não necessariamente as pro- priedades continuam valendo. Assim, torna-se importante a definição de outros tipos de reśıduo cujas propriedades sejam conhecidas ou pelo menos estejam mais próximas das propriedades de t∗i . Uma primeira proposta seria considerar o reśıduo ordinário da solução de mı́nimos quadrados da regressão linear ponderada de z contra X, que é definido por r∗ = Ŵ1/2[z− η̂] = V̂−1/2(y − µ̂). Se assumirmos que Var(z) ∼= Ŵ−1φ−1, temos aproximadamente Modelos Lineares Generalizados 49 Var[r∗] ∼= φ−1(I − Ĥ). Logo, podemos definir o reśıduo padronizado tSi = φ1/2(yi − µ̂i) √ V̂i(1 − ĥii) , em que hii é o i-ésimo elemento da diagonal principal da matriz H. Fica fácil mostrar que r∗ = (I − Ĥ)Ŵ1/2z, isto é, Ĥ desempenha o papel de matriz de projeção ortogonal local, como na regressão normal linear em que W é identidade. No entanto, na prática, η̂ não é fixo nem conhecido, bem como z não segue distribuição normal. Uma implicação desse fato é que as propriedades de t∗i não são mais verificadas para tSi . Williams (1984) mostra através de estudos de Monte Carlo que a distribuição de tSi é em geral assimétrica, mesmo para grandes amostras. Outros reśıduos cujas distribuições poderiam estar mais próximas da normalidade têm sido sugeridos para os MLGs. Por exemplo, o reśıduo de Anscombe tAi = φ1/2{ψ(yi) − ψ(µ̂i)} V̂ 1/2(µ̂i)ψ′(µ̂i) , em que ψ(·) é uma transformação utilizada para normalizar a distribuição de Y . Para os MLGs essa transformação é definida por ψ(µ) = ∫ V −1/3(µ)dµ. Em particular para os MLGs, a função ψ(µ) vale µ, ∫ µ−1/3(1 − µ)−1/3dµ, 3 2 µ2/3, 3µ1/3 e logµ para a normal, binomial, Poisson, gamma e normal inversa, respectivamente. Contudo, os reśıduos mais utilizados em modelos lineares generalizados são definidos a partir dos componentes da função desvio. A versão padronizada (vide McCullagh, 1987; Davison e Gigli, 1989) é a seguinte: tDi = d∗(yi; µ̂i) √ (1 − ĥii) = φ1/2d(yi; µ̂i) √ (1 − ĥii) , em que d(yi; µ̂i) = ± √ 2{yi(θ̂0i − θ̂i) + (b(θ̂i) − b(θ̂0i ))}1/2. O sinal de d(yi; µ̂i) é o mesmo de yi − µ̂i. Williams (1984) verificou através de simulações que a distribuição de tDi tende a estar mais próxima da normalidade do que as distribuições dos demais reśıduos. McCullagh (1987, p. 214) mostra para os MLGs que a distribuição de probabilidades de d∗(Yi;µi) + ρ3i/6 √ 1 + (14ρ23i − 9ρ4i)/36 52 Caṕıtulo 1 Essa aproximação, introduzida por Pregibon (1981), é dada por β1(i) = β̂ − r̂Pi √ ω̂iφ−1 (1 − ĥii) (XTŴX)−1xi. (1.19) Logo, substituindo a expressão acima em (1.18), obtém-se LDi ∼= { ĥii (1 − ĥii) } t2Si. A distância de Cook aproximada fica facilmente obtida com o comando LD < − h*(ts^ 2)/(1 - h) A validade da aproximação de um passo tem sido investigada por alguns pesquisadores. A constatação é que a mesma em geral subestima o verdadeiro valor de LDi, no entanto é suficiente para chamar a atenção dos pontos aberrantes e influentes. 1.9.4 Influência local Cook (1986) mostra que a extensão do método de influência local para os MLGs segue diretamente quando a ligação é canônica. Nesse caso, o vetor dmax para avaliar a influência local das observações nas estimativas dos parâmetros é o autovetor correspondente ao maior autovalor da seguinte matriz n× n: A = diag(r̂P )Ĥdiag(r̂P ), em que r̂P = (r̂P1, . . . , r̂Pn) T e r̂Pi = φ 1/2(yi − µ̂i)/V̂ 1/2 é o i-ésimo reśıduo de Pearson avaliado em β̂. Para obter dmax, a maneira mais simples é construir a matriz A e extrair o seu au- tovetor correspondente ao maior autovalor. Os comandos são os seguintes: A < − diag(rp)%*% H %*% diag(rp) Lmax < − eigen(A)$val[1] dmax < − eigen(A)$vec[,1] dmax < − dmax/sqrt(Lmax) dmax < − abs(dmax) Por outro lado, se o interesse é detectar as observações influentes na estimativa de um coeficiente particular, associado por exemplo à variável explicativa X1, o vetor dmax fica dado por dTmax = ( v1r̂P1√ λmax , . . . , vnr̂Pn√ λmax ) , Modelos Lineares Generalizados 53 em que v1, . . . , vn são agora obtidos da regressão linear de X1 contra as colunas de X2 com matriz de pesos V̂, isto é v = V̂1/2X1 − V̂1/2X2(XT2 V̂X2)−1XT2 V̂X1. Para ligação não canônica os resultados continuam valendo desde que a matriz observada de Fisher seja substitúıda pela matriz de informação de Fisher. 1.9.5 Gráfico da variável adicionada Apresentamos a seguir a versão do gráfico da variável adicionada para os MLGs. Suponha um MLG com p parâmetros, β1, . . . , βp, e que um parâmetro adicional γ está sendo inclúıdo no modelo. O interesse é testar H0 : γ = 0 contra H1 : γ 6= 0. Seja η(β, γ) o preditor linear com p+ 1 parâmetros, isto é η(β, γ) = XTβ + γZ. A função escore para γ é dada por Uγ(β) = ∂L(β, γ) ∂γ = φ1/2ZTW1/2rP , em que Z = (z1, . . . , zn) T . De resultados anteriores temos que Var(γ̂) = φ−1[ZTW1/2MW1/2Z]−1, em que M = I − H. Logo, Var(γ̂) = φ−1(RTWR)−1 com R = Z − XC e C = (XTWX)−1XTWZ. Portanto, a estat́ıstica de escore para testar H0 : γ = 0 contra H1 : γ 6= 0 fica dada por ξSR = (r̂ T PŴ 1/2Z)2/(ZTŴ1/2M̂Ŵ1/2Z), em que Ŵ, r̂P e M̂ são avaliados em β̂ (sob H0). Sob H0, ξSR ∼ χ21 quando n→ ∞. Mostra-se (Wang, 1985), que a estat́ıstica de escore acima coincide com a estat́ıstica F de uma regressão linear ponderada para testar a inclusão da variável Z no modelo. Nessa regressão linear, o gráfico da variável adicionada é formado pelos reśıduos r̂P e υ = φ1/2(I − Ĥ)Ŵ1/2Z. O reśıduo υ pode ser obtido facilmente após a regressão linear ponderada (com pesos Ŵ) de Z contra X. Note que γ̂ = (υTυ)−1υT r. Logo, o gráfico de r̂P contra υ pode revelar quais observações estão contribuindo mais na significância de γ. A principal dificuldade para construir o gráfico da variável 54 Caṕıtulo 1 adicionada em MLGs é a obtenção do reśıduo υ, uma vez que o reśıduo r̂P é obtido facil- mente como já vimos anteriormente. Para ilustrar o cálculo de υ num modelo particular, suponha que temos duas covariáveis e dois fatores e que o interesse é construir o gráfico da variável adicionada correspondente à covariável cov1. Precisamos inicialmente ajustar o modelo com os dois fatores e a outra covariável e computar a matriz Ŵ cujos valores serão armazenados em W. Lembrando que Ŵ é a matriz estimada de pesos. Supondo, por exemplo, que temos um modelo de Poisson com ligação canônica, os passos para construir o gráfico são os seguintes: fit.poisson < − glm( resp ∼ cov2 + A + B, family=poisson) w < − fit.poisson$weights W < − diag(w) rp < − resid(fit.poisson, type =‘‘pearson ") X < − model.matrix(fit.poisson) H < − solve(t(X)%*%W%*%X) H < − sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W) v < − sqrt(W)%*%cov1 - H%*%sqrt(W)%*%cov1 plot(v, rp, xlab=‘‘Residuo v ’’, ylab=‘‘Residuo rp ’’) 1.9.6 Seleção de modelos Os métodos de seleção de modelos descritos na Seção 1.8.4 podem ser estendidas dire- tamente para os MLGs. Algumas observações, contudo, se fazem necessárias. Nos casos de regressão loǵıstica e de Poisson o teste da razão de verossimilhanças, pelo fato de ser obtido pela diferença de duas funções desvio, aparece como o mais indicado. Para os casos de regressão normal e gama o teste F, por não exigir a estimativa de máxima verossim- ilança do parâmetro de dispersão, é o mais indicado. Isso não impede que outros testes sejam utilizados. Já o método de Akaike pode ser expresso numa forma mais simples em função do desvio do modelo. Nesse caso, o critério consiste em encontrar o modelo tal que a quantidade abaixo seja minimizada AIC = Dp + 2p, em que Dp denota o desvio do modelo e p o número de parâmetros. Os métodos stepwise e de Akaike estão dispońıveis no S-Plus. O método stepwise está dispońıvel apenas para modelos normais lineares. O comando stepwise é definido por stepwise(Xvar, Modelos Lineares Generalizados 57 Propomos inicialmente um modelo normal linear simples em que Y denote a renda e X a escolaridade. O modelo fica portanto dado por yi = α+ βxi + i, i = 1, . . . , 27, com a suposição de que i ∼ N(0, σ2), sendo os erros mutuamente independentes. Escolaridade R en da 4 5 6 7 8 40 0 60 0 80 0 12 00 (a) DF Indice A la va nc a 0 5 10 15 20 25 0. 05 0. 15 0. 25 (b) DF Indice D is ta nc ia d e C oo k 0 5 10 15 20 25 0 1 2 3 4 5 (c) DF Valores Ajustados R es id uo S tu de nt iz ad o 400 600 800 1000 1200 -2 0 2 4 (d) DF Figura 1.5: Reta ajustada do modelo aditivo e gráficos de diagnóstico para o exemplo sobre escolaridade e renda. As estimativas dos parâmetros (desvio padrão) são dadas por α̂ = −381, 28 (69, 40) e β̂ = 199, 82 (13, 03), indicando que o coeficiente angular da reta é altamente significativo. Essa estimativa pode ser interpetada como o incremento esperado na renda média domi- ciliar de uma unidade da federação se o tempo de escolaridade médio domiciliar naquela unidade for acrescido de um ano. A estimativa de σ2 é dada por s2 = 77, 22, enquanto que o coeficiente de determinação foi de R2 = 0, 904. O ajuste do modelo e a exibição dos resultados podem ser obtidos com os comandos abaixo attach(censo.dat) fit1.censo < − lm(renda ∼ escolar) 58 Caṕıtulo 1 summary(fit1.censo) Ou, alternativamente, transformando o arquivo censo.dat num arquivo do tipo data frame, através dos comandos censo.dat < − data.frame(censo.dat) fit1.censo < − lm(renda ∼ escolar, data=censo.dat) summary(fit1.censo) Escolaridade Lo g( R en da ) 4 5 6 7 8 6. 0 6. 5 7. 0 (a) DF Indice A la va nc a 0 5 10 15 20 25 0. 05 0. 15 0. 25 (b) DF Indice D is ta nc ia d e C oo k 0 5 10 15 20 25 0. 0 0. 1 0. 2 0. 3 0. 4 0. 5 (c) MA Valores Ajustados R es id uo S tu de nt iz ad o 6.0 6.2 6.4 6.6 6.8 7.0 7.2 -2 -1 0 1 2 (d) RO MA MT Figura 1.6: Reta ajustada do modelo multiplicativo e gráficos de diagnóstico para o exemplo sobre escolaridade e renda. Pela Figura 1.5 onde são apresentados alguns gráficos de diagnóstico além da reta ajustada aos dados nota-se uma forte discrepância do Distrito Federal que aparece como ponto de alavanca, influente e aberrante. Além disso, nota-se pela Figura 1.5d ind́ıcios de heterocedasticidade, ou seja, um aumento da variabilidade com o aumento da escolar- idade. Isso pode também ser notado na Figura 1.5a. Assim, pode-se propor um modelo alternativo, por exemplo, com efeitos multiplicativos conforme dado abaixo logyi = α+ βxi + i, i = 1, . . . , 27, Modelos Lineares Generalizados 59 com a suposição de que i ∼ N(0, σ2), sendo os erros mutuamente independentes. Na Figura 1.6 tem-se o ajuste do modelo acima aos dados bem como alguns gráficos de diagnóstico que destacam DF como ponto de alavanca e MA como ponto influente além de aberrante. A Tabela 1.7 faz uma análise confirmatória em que verifica-se poucas variações nas estimativas dos parâmetros com a eliminação dessas unidades da federação. Finalmente, na Figura 1.7 tem-se os gráficos de diagnóstico para o modelo com efeitos aditivos (Figura 1.7a) e com efeitos multiplicativos (Figura 1.7b) e nota-se uma melhor acomodação e distribuição dos pontos dentro do envelope gerado no segundo caso. Percentis da N(0,1) R e s id u o S tu d e n ti z a d o -2 -1 0 1 2 -2 0 2 4 (a) Percentis da N(0,1) R e s id u o S tu d e n ti z a d o -2 -1 0 1 2 -3 -2 -1 0 1 2 3 (b) Figura 1.7: Gráficos normais de probabilidades para os modelos ajustados aditivo (a) e multiplicativo (b). Tabela 1.7 Estimativas de algumas quantidades com todos os pontos e quando as observações mais discrepantes são exclúıdas. Estimativa Com todos Exclúıdo Exclúıdo Exclúıdos os pontos DF MA DF e MA α̂ 5,065 (0,075) 4,982 (0,067) 5,028 (0,065) 5,006 (0,077) β̂ 0,264 (0,014) 0,279 (0,013) 0,271 (0,012) 0,274 (0,015) s2 0,069 0,075 0,069 0,076 R2 93,7% 95,1% 95,4% 93,4% 62 Caṕıtulo 1 Na Figura 1.8a apresentamos o gráfico do número de bactérias sobreviventes contra o tempo de exposição. Nota-se uma tendência decrescente e quadrática. Supondo que as amostras do produto enlatado submetidos à temperatura de 300oF têm o mesmo tamanho, pode-se pensar, em prinćıpio, que Yi ∼ P (µi), em que Yi denota o número de bactérias sobreviventes na i-ésima amostra i = 1, . . . , n. Como para µi grande é razoável assumir que Yi segue uma distribuição aproximadamente normal (vide Seção 2.2.1), propomos inicialmente os seguintes modelos: yi = α + βtempoi + i e yi = α + βtempoi + γtempo 2 i + i, em que i ∼ N(0, σ2). As estimativas dos parâmetros são apresentadas na Tabela 1.10. Pelos gráficos de envelopes (Figuras 1.8b e 1.8c) nota-se ind́ıcios de que a distribuição dos erros pode estar incorrretamente especificada. A maioria dos reśıduos assume valor negativo. Nota-se a presença de um ponto aberrante, observação # 1. Uma outra tentativa seria aplicar à resposta a transformação raiz quadrada que é conhecida no caso da Poisson como estabilizadora da variância além de manter a aproximação normal (vide Seção 2.2.1). Logo, podemos pensar em adotar os seguintes modelos alternativos: √ yi = α + βtempoi + i e √ yi = α + βtempoi + γtempo 2 i + i, em que i ∼ N(0, σ2). As estimativas dos parâmetros encontram-se na Tabela 1.10. Nota-se uma melhora na qualidade do ajuste, particularmente no segundo caso. Porém, ainda há ind́ıcios pelos gráficos de envelopes (Figuras 1.8d e 1.8e) de violação nas su- posições para os modelos, além da presença da observação # 1 como ponto aberrante. Decidimos, então, propor um modelo log-linear de Poisson em que assumimos Yi ∼ P (µi) e logµi = α + βtempoi. As estimativas dos parâmetros são também apresentadas na Tabela 1.10. Pelo gráfico de envelope (Figura 1.8f) não há evidências de que o modelo esteja mal ajustado. Nota-se também que a observação #1 foi acomodada dentro do envelope gerado. Parece, portanto, que esse último modelo é o que melhor se ajusta aos dados dentre os modelos propostos. O modelo ajustado fica então dado por µ̂(x) = e5,30−0,23x, Modelos Lineares Generalizados 63 em que x denota o tempo de exposição. Logo, se diminuirmos de uma unidade o tempo de exposição a variação no valor esperado fica dada por µ̂(x− 1) µ̂(x) = e0,23 = 1, 259. Ou seja, o número esperado de sobreviventes aumenta aproximadamente 25,9%. Tempo So br ev iv en te s 2 4 6 8 10 12 50 10 0 15 0 (a) 1 Percentis da N(0,1) R es id uo S tu de nt iz ad o -1 0 1 -2 0 2 4 6 8 (b) Percentis da N(0,1) R es id uo S tu de nt iz ad o -1 0 1 -2 0 2 4 6 8 (c) Percentis da N(0,1) R es id uo S tu de nt iz ad o -1 0 1 -2 0 2 4 6 (d) Percentis da N(0,1) R es id uo S tu de nt iz ad o -1 0 1 -2 0 2 4 (e) Percentis da N(0,1) C om po ne nt e do D es vi o -1 0 1 -3 -2 -1 0 1 2 (f) Figura 1.8: Diagrama de dispersão e gráficos normal de probabilidades para o exemplo sobre sobrevivência de bactérias. Tabela 1.10 Estimativas de algumas quantidades para os cinco modelos propostos. Estimativa Linear-Y Quadrático-Y Linear- √ Y Quadrático- √ Y Poisson α̂ 142,20(11,26) 181,20(11,64) 12,57(0,38) 13,64(0,51) 5,30(0,06) β̂ -12,48(1,53) -29,20(4,11) -0,82(0,05) -1,27(0,18) -0,23(0,01) γ̂ 1,29(0,31) 0,04(0,01) R2 86,9% 95,5% 96,1% 97,8% Desvio 8,42 (10 g.l.) 64 Caṕıtulo 1 1.10.4 Estudo seriado com ratos O exemplo a seguir provém de um estudo seriado com um tipo de tumor maligno para avaliar a influência da série (passagem do tumor) na morte (caquexia) de um certo tipo de rato (vide Paula, Barbosa e Ferreira, 1989; Paula et al., 1992). Os dados estão descritos no arquivo canc4.dat. Um total de 204 animais teve o tumor inoculado num determinado momento da série. Para cada animal, além do grupo de passagem, foram observadas as variáveis presença de massa tumoral, caquexia e o tempo de observação (em dias). Esses dados são resumidos na Tabela 1.11. Para inserir os dados diretamente no S-Plus e armazená-los no arquivo canc4a.dat, devemos fazer canc4a.dat < − scan(what=list(obs=0,rd=0)) 1: 6 2597 13 3105 8 2786 2: 12 1613 3 411 1 232 Agora, precisamos introduzir os fatores grupo de passagem e massa tumoral fnames < − list(gp=c(“ P0-P6 ", “ P7-P18", “ P19-P28"), mt=c(“ sim", “ nao")) Para informar o sistema a ordem em que os dados foram lidos, pode-se usar o comando fac.design. Em seguida, fazemos o emparelhamento rato.design < − fac.design(c(3,2), fnames, rep=1) attach(canc4a.dat) rato.df < − data.frame(obs,rd,rato.design) As informações completas sobre os dados estão armazenadas no arquivo rato.df. Para uma verificação basta bater rato.df Podemos agora (opcionalmente) criar uma matriz modelo no padrão dos MLGs attach(rato.df) gp < − C(gp,treatment) mt < − C(mt,treatment) Vamos supor que Oij , o número de ratos caquéticos no ńıvel i de massa tumoral e grupo de passagem j, segue uma distribuição de Poisson de média λijtij , i = 1, 2 e j = 1, 2, 3. Note que λij denota a taxa de caquexia (número médio de mortes por unidade de tempo) e tij o total de ratos-dias no ńıvel (i, j). Considere inicialmente o modelo log-linear logλij = α + βi + γj, Modelos Lineares Generalizados 67 1: 3.03 3.19 3.46 5.88 6.43 2: 5.53 4.26 5.22 6.74 9.97 3: 5.60 4.47 5.69 6.90 10.39 4: 9.30 4.53 6.54 6.98 13.55 5: . . . Denotaremos por Tij o tempo até a perda da velocidade para o j-ésimo motor de tipo i, i = 1, . . . , 5 e j = 1, . . . , 10. Na tabela abaixo são apresentadas as médias, desvios padrão e coeficientes de variação amostrais para os cinco tipos de turbina e como pode-se notar os coeficientes de variação variam menos que os desvios padrão. Isso sugere que uma distribuição gama com coeficiente de variação constante pode ser mais apropriada para explicar o tempo de duração do que uma distribuição normal com variância constante. Tipo I Tipo II Tipo III Tipo IV Tipo V Média 10,69 6,05 8,64 9,80 14,71 D.Padrão 4,82 2,91 3,29 5,81 4,86 C. Variação 45,09% 48,10% 38,08% 59,29% 33,04% Vamos assumir então que Tij segue uma distribuição gama de média µi e parâmetro de dispersão φ−1. Para comparar os cinco grupos utilizaremos inicialmente o modelo abaixo (modelo gama com ligação canônica) µ−1i = µ+ βi, em que β1 = 0. É importante observar que os resultados seriam os mesmos se fosse utilizada qualquer outra ligação. Para ajustar o modelo no S-Plus precisamos definir antes o fator tipo de turbina e fazer o emparelhamento dos dados com os ńıveis do mesmo. Os comandos são apresentados abaixo fnames < − list(tipo=c(“ I ", “ II ", “ III ", “ IV ", “ V ")) turbina.design < − fac.design(5,fnames,rep=10) attach(turbina.dat) turbina.df < − data.frame(tempo, turbina.design) turbina.df Os boxplots correspondentes aos tempos dos cinco grupos (vide Figura 1.10a) são obtidos com os comandos attach(turbina.df) 68 Caṕıtulo 1 plot.factor(turbina.df) Os passos para o ajuste do modelo são dados a seguir tipo < − C(tipo,treatment) fit.turbina < − glm(tempo ∼ tipo, family=Gamma) summary(fit.turbina) Indice D is ta nc ia d e C oo k 0 10 20 30 40 50 0. 0 0. 5 1. 0 1. 5 47 49 (a) Preditor Linear R es id uo C om po ne nt e do D es vi o 6 8 10 12 14 -2 0 2 4 (b) 1 47 49 Figura 1.9: Distância de Cook (a) e componente do desvio contra preditor linear (b) para o exemplo sobre desempenho de turbinas de avião. O desvio do modelo foi de D∗(y; µ̂) = 8, 861 × 5, 804 = 51, 43, com 45 graus de liber- dade, que leva a P = 0, 236 indicando um ajuste adequado. As estimativas dos parâmetros deram µ̂ = 0, 094 (0, 013), β̂2 = 0, 072 (0, 027), β̂3 = 0, 022 (0, 021), β̂4 = 0, 008 (0, 019) e β̂5 = −0, 025 (0, 017), indicando para o tipo II um tempo médio de sobrevivência signi- ficativamente menor do que os demais. Para o tipo V notamos um tempo médio maior do que os demais enquanto que os outros três tipos apresentam tempos médios significa- tivamente não diferentes. Esses resultados confirmam a análise descritiva apresentada na Figura 1.10a. A estimativa de máxima verossimilhança (desvio padrão aproximado) do parâmetro de dispersão foi de φ̂ = 5, 804(1, 129)), indicando que as distribuições dos tem- pos de sobrevivência não devem ser muito assimétricas. Na Figura 1.9 tem-se o gráfico da distância de Cook (Figura 1.9a) e o gráfico do componente do desvio padronizado contra o preditor linear (Figura 1.9b). Nota-se um forte destaque para a observação #49 que corresponde ao valor 25,46 para o tempo de duração de um dos motores de tipo IV. Esse valor, como mostra o boxplot correspondente na Figura 1.10 destoa dos demais tempos. Modelos Lineares Generalizados 69 A eliminação da observação #49 aumenta a significância marginal de β4, embora esse efeito continue não significativo a 10%. O gráfico normal de probabilidades com envelope para os componentes padronizados do desvio é apresentado na Figura 1.10b. Notamos, pelo gráfico, que não há ind́ıcios de afastamentos sérios da suposição de distribuição gama para os tempos de sobrevivência dos motores bem como para a suposição de homogeneidade de coeficiente de variação para os cinco grupos. A sequência de comandos para construir o gráfico normal de probabili- dades com envelopes é dada no Apêndice. É assumido que os resultados do ajuste estão guardados no objeto fit.model. 5 1 0 1 5 2 0 2 5 te m p o I II III IV V tipo (a) Percentis da N(0,1) C o m p o n e n te d o D e s v io -2 -1 0 1 2 -3 -2 -1 0 1 2 (b) Figura 1.10: Box-plot (a) e gráfico normal de probabilidades (b) para o exemplo sobre desempenho de turbinas de avião. A fim de facilitar as interpretações dos resultados de um modelo gama ou mesmo fazer comparações com o modelo normal linear, pode-se propor uma ligação identidade ao invés de ligação rećıproca. No exemplo das turbinas a parte sistemática do modelo ficaria dada por µi = µ+ βi, em que β1 = 0. Para ajustar o modelo no S-Plus deve-se fazer o seguinte: fit1.turbina < glm(tempo ∼ tipo, family=Gamma(link=identity)) 72 Caṕıtulo 1 Percentis da N(0,1) R e s id u o S tu d e n ti z a d o -2 -1 0 1 2 -2 0 2 4 (a) Percentis da N(0,1) R e s id u o S tu d e n ti z a d o -2 -1 0 1 2 -3 -2 -1 0 1 2 3 (b) Figura 1.12: Gráficos normais de probabilidades com todos os pontos (a) e sem o estado de WY (b), para o exemplo sobre consumo de combust́ıvel. Portanto, podemos dizer que para cada aumento de uma unidade na renda, o consumo médio de combust́ıvel diminui 0,07 unidades. Para cada aumento de 1% na porcentagem de motoristas licenciados o consumo médio de combust́ıvel aumenta 13,75 unidades, e para cada aumento de 1% no imposto do combust́ıvel o consumo médio diminui 29,48 unidades. Na Figura 1.11 temos alguns gráficos de diagnóstico e como podemos notar há um forte destaque para o estado de WY, que aparece como influente (Figura 1.11b) e aberrante (Figura 1.11c). Outros estados, tais como CT, NY, SD, TX e NV (Figura 1.11a) apare- cem como remotos no subespaço gerado pelas colunas da matrix X, embora não sejam confirmados como influentes. Não há ind́ıcios pela Figura 1.11d de heterocedasticidade. Pelo gráfico de envelope (Figura 1.12a) não há ind́ıcios fortes de afastamentos sérios da suposição de normalidade para os erros, apesar da influência no gráfico do estado de WY. O gráfico de envelope sem esse estado (Figura 1.12b) confirma esse suposição. Analisando os dados referentes ao estado de WY notamos que o mesmo tem uma taxa de 7% (abaixo da média de 7,67%), uma renda per-capita anual de US$ 4345 (ligeiramente acima da média de US$ 4241,83), uma proporção de motoristas licenciados de 0,672 (acima da média de 0,570), porém um consumo médio de combust́ıvel muito alto 968 (quando a média nacional era de 576,77). Talvez as longas distâncias do estado tenham obrigado os motoristas a um consumo alto de combust́ıvel. A eliminação desse estado muda substacialmente algumas estimativas, embora não mude as tendências. A estimativa Modelos Lineares Generalizados 73 da variável licença cai 13,2%, a estimativa do intercepto aumenta 27,8%, o s2 cai 17,1% e o R2 aumenta 4,1%. As demais estimativas não sofrem grandes variações. 1.11 Exerćıcios 1. Seja Y uma variável aleatória com distribuição binomial negativa, isto é, Y é o número de ensaios até a ocorrência do r-ésimo sucesso, em que π é a probabilidade de sucesso em cada ensaio. Mostre que a função de probabilidades de Y pode ser expressa na forma exponencial. Calcule µ e V (µ). Use a forma abaixo para a função de probabilidades de Y f(y; π, r) = ( y − 1 r − 1 ) πr(1 − π)(y−r), em que y = r, r + 1, . . .. 2. Considere a seguinte função densidade de probabilidade: f(y; θ, φ) = φa(y, φ) π(1 + y2)1/2 exp[φ{yθ + (1 − θ2)1/2}], em que 0 < θ < 1, −∞ < y < ∞, φ > 0 e a(·, ·) é uma função normalizadora. (i) Mostre que essa distribuição pertence à famı́lia exponencial; (ii) encontre E(Y ) = µ e V (µ); (iii) obtenha o reśıduo de Pearson e (iv) encontre a função desvio supondo uma amostra de n variáveis aleatórias independentes. 3. Mostre que a distribuição logaŕıtmica, com função de probabilidades f(y; ρ) = ρy/{−ylog(1 − ρ)}, em que y = 1, 2, . . . e 0 < ρ < 1, pertence à famı́lia exponencial. Calcule µ e V (µ). 4. Considere a distribuição estável cuja densidade é dada por f(y; θ, φ) = a(y, φ)exp[φ{θ(y + 1) − θlogθ}], em que θ > 0, −∞ < y < ∞, φ−1 > 0 é o parâmetro de escala e a(·, ·) é uma função normalizadora. Mostre que essa distribuição pertence à famı́lia exponencial. Encontre µ e V (µ). Obtenha a função desvio supondo uma amostra de n variáveis aleatórias independentes. 74 Caṕıtulo 1 5. Encontre a função desvio para as distribuições binomial negativa e logaŕıtmica. Mostre que o desvio da distribuição gama para o caso i.i.d é dado por D∗(y; µ̂) = 2nφlog(ȳ/ỹ), em que ỹ é a média geométrica das observações. 6. (Paula e Cordeiro, 1986). Suponha o modelo g(µ;λ) = η, em que η = Xβ com λ univariado. Mostre que o processo iterativo para estimar (βT , λ) é o mesmo de um MLG com parte sistemática g(µ, λ) = Xβ + Λλ, em que a matriz modelo é dada por X̃ = [X,Λ] e Λ = ∂η/∂λ. Particularize esse processo iterativo para as ligações Box-Cox e de Aranda-Ordaz. 7. Supor o modelo normal linear com parte sistemática dada por ηi = β1(x1i − x̄1) + β2(x2i− x̄2). Sabe-se que a correlação amostral entre x1 e x2 é dada por corr(x1, x2) = ∑n i=1(x1i − x̄1)(x2i − x̄2)/(n − 1)s1s2, em que s1 e s2 são os respectivos desvios- padrão amostrais de x1 e x2. Calcule a correlação corr(β̂1, β̂2). Discuta e tente explicar a relação entre as duas correlações. Use o fato de que det(XTX)−1 > 0. 8. Suponha o modelo de análise de variância com erros normais yij = α + βi + ij , em que ij ∼ N(0, σ2), i = 1, . . . , p e j = 1, . . . , ni. Supor β1 = 0. Mostre que Var(rij) = σ 2(1 − 1/ni). 9. Considere o modelo normal linear yi = x T i β + i, i = 1, . . . , n, em que i são mutuamente independentes tais que i ∼ N(0, σ2). Considere uma nova observação y(z) (que não está na amostra) e que satisfaz y(z) = zTβ + , em que  ∼ N(0, σ2). Mostre que um intervalo de confiança de coeficiente 1 − α para y(z) pode ser dado por [ŷ(z) ± tn−p(1 − α 2 )s{1 + zT (XTX)−1z}1/2], em que ŷ(z) = zT β̂, tn−1(1 − α2 ) é o percentil (1 − α2 ) da distribuição t de Student com n− p graus de liberdade e s2 é o erro quadrático médio do modelo ajustado. 10. Suponha agora o modelo de regressão normal linear simples yi = α + βxi + i, i = 1, . . . , n. Modelos Lineares Generalizados 77 20. O conjunto de dados descrito na tabela abaixo refere-se a um estudo cujo objetivo foi tentar prever o preço de venda de um imóvel (em US$ mil) dada a área total (em pés quadrados) numa região de Eugene, EUA (Gray, 1989). Esses dados estão armazenados no arquivo externo reg1.dat. Área 800 950 910 950 1200 1000 1180 1000 1380 1250 Preço 30,6 31,5 33,3 45,9 47,4 48,9 51,6 53,1 54,0 54,3 Área 1500 1200 1600 1650 1600 1680 1500 1780 1790 1900 Preço 55,2 55,2 56,7 57,9 58,5 59,7 60,9 60,9 62,4 63,0 Área 1760 1850 1800 1700 1370 2000 2000 2100 2050 1990 Preço 64,5 66,0 66,3 67,5 68,4 68,4 68,7 69,6 70,5 74,7 Área 2150 2050 2200 2200 2180 2250 2400 2350 2500 2500 Preço 75,0 75,3 79,8 80,7 80,7 83,4 84,0 86,1 87,0 90,3 Área 2500 2500 2680 2210 2750 2500 2400 3100 2100 4000 Preço 96,0 101,4 105,9 111,3 112,5 114,0 115,2 117,0 129,0 165,0 Tente inicialmente ajustar uma regressão normal linear para explicar o preço dada a renda. Faça uma análise de diagnóstico e proponha algum modelo alternativo (se for o caso) a fim de reduzir as eventuais influências de observações discrepantes bem como afastamentos de outras suposições feitas para o modelo. Interprete as estimativas obtidas para os coeficientes do modelo proposto. 21. (Pregibon, 1982). Mostre que o teste de escore para testar que o i-ésimo ponto é aberrante num MLG é dado por t2Si. Sugestão : chame η = x Tβ+γz, em que z é um vetor n× 1 de zeros com 1 na i-ésima posição. Qual a distribuição nula assintótica de t2Si? 22. Mostrar que a expressão para AIC no modelo normal linear com σ2 desconhecido pode ser expressa na forma equivalente AIC = nlog{D(y; µ̂)/n} + 2p, em que D(y; µ̂) = ∑n i=1(yi − µ̂i)2. 78 Caṕıtulo 1 23. Sejam Yi ∼ FE(µ1, φ1), i = 1, . . . , m, e Yi ∼ FE(µ2, φ2), i = m+ 1, . . . , n, variáveis aleatórias mutuamente independentes. Encontre a estimativa comum de máxima verossimilhança para φ1 e φ2 sob a hipótese H0 : φ1 = φ2. Particularize para os casos gama e normal. 24. No arquivo reg3.dat são descritas as seguintes variáveis referente a 50 estados norte-americanos: (i) nome (nome do estado), (ii) pop (população estimada em julho de 1975), (iii) renda (renda per-capita em 1974), (iv) tt analf (porporção de analfabetos em 1970), (v) expvida (expectativa de vida em anos 1969-70), (vi) crime (taxa de criminalidade por 100000 habitantes 1976), (vii) estud (porcentagem de estudantes que concluem o segundo grau 1970), (viii) temp (número de dias do ano com temperatura abaixo de zero grau Celsus na cidade mais importante do estado) e (ix) area (área do estado em milhas quadradas). Tente explicar e variável expvida usando um modelo de regressão normal linear dadas as variáveis explicativas renda, analf, crime, estud, temp e dens, em que dens=pop/area. Aplique o método stepwise de seleção de modelos. Faça uma análise completa de diagnóstico com o modelo selecionado. Interprete os resultados. 25. (Neter et el., 1996, p. 449) No arquivo vendas.dat são descritas informações a respeito das vendas no ano anterior de um tipo de telhado de madeira em 26 filiais de uma rede de lojas de construção. As variáveis estão colocadas na seguinte ordem: (i) telhados, total de telhados vendidos (em mil metros quadrados), (ii) gastos, gastos pela loja com promoções do produto (em mil US$), (iii) clientes, número de clientes cadastrados na loja (em milhares), (iv) marcas, número de marcas concor- rentes do produto e (v) potencial, potencial da loja (quanto maior o valor maior o potencial). Um dos objetivos do estudo com esse conjunto de dados é tentar prever o número esperado de telhados vendidos dadas as variáveis explicativas. Faça inicial- mente uma análise descritiva construindo, por exemplo, os diagramas de dispersão de cada variável explicativa contra a variável resposta telhados. Calcule também as correlações entre as variáveis. Use os métodos stepwise e AIC para selecionar um modelo de regressão normal linear. Se o modelo selecionado for diferente pelos dois métodos, adote algum critério para escolher um dos modelos. Interprete os coeficientes estimados do modelo selecionado. Faça uma análise de diagnóstico para verificar se existem afastamentos sérios das suposições feitas para o modelo e se Modelos Lineares Generalizados 79 existem observações discrepantes. 26. (Wood, 1973). No arquivo reg4.dat estão os dados referentes à produção de gasolina numa determinada refinaria segundo três variáveis observadas durante o processo e uma quarta variável que é uma combinação das três primeiras. A resposta é o número de octanas do produto produzido. A octanagem é a propriedade que de- termina o limite máximo que a gasolina, junto com o ar, pode ser comprimida na câmara de combustão do véıculo sem queimar antes de receber a centilha vinda das velas. As melhores gasolinas têm uma octanagem alta. Em grandes refinarias, o aumento de um octana na produção de gasolina pode representar um aumento de alguns milhões de dolares no custo final da produção. Assim, torna-se impor- tante o controle dessa variável durante o processo de produção. Use o método stepwise para selecionar as variáveis explicativas significativas. Faça uma análise de diagóstico com o modelo selecionado. Comente. 27. (Narula e Stangenhaus, 1988, p. 32) No arquivo imoveis.dat são apresentados dados relativos a uma amostra de 27 imóveis. Na ordem são apresentados os valores das seguintes variáveis: (i) imposto do imóvel (em 100 dolares), (ii) área do terreno (em 1000 pés quadrados), (iii) área constrúıda (em 1000 pés quadrados), (iv) idade da residência (em anos) e (v) preço de venda do imóvel (em 1000 dolares). Ajuste um modelo normal linear do preço de venda contra as demais variáveis. Use o método AIC para selecionar as variáveis explicativas. Faça uma análise de diagnóstico com o modelo selecionado. Interprete os coeficientes estimados. 28. (Paula e Oshiro, 2001). O espinhel de fundo é definido como um método de pesca passivo, sendo utilizado em todo o mundo em operações de pesca de diferentes magnitudes, da pesca artesanal a modernas pescarias mecanizadas. É adequado para capturar peixes com distribuição dispersa ou com baixa densidade, além de ser posśıvel utilizá-lo em áreas irregulares ou em grandes profundidades. É um dos métodos que mais satisfazem às premissas da pesca responsável, com alta seletivi- dade de espécies e comprimentos, alta qualidade do pescado, consumo de energia baixo e pouco impacto sobre o fundo oceânico. No arquivo pesca.dat estão parte dos dados de um estudo sobre a atividade das frotas pesqueiras de espinhel de fundo baseadas em Santos e Ubatuba no litoral paulista. A espécie de peixe considerada é o peixe-batata pela sua importância comercial e ampla distribuição espacial. As 82 Caṕıtulo 1 uma reparametrização tipo casela de referência em que µ11 = α, µ1j = α + βj, µ21 = α+ γ e µ2j = α + γ + βj j = 2, 3, 4. Voltagem(kV) Temperatura (oC) 200 250 300 350 170 439 572 315 258 904 690 315 258 1092 904 439 347 1105 1090 628 588 180 959 216 241 241 1065 315 315 241 1065 455 332 435 1087 473 380 455 Procure responder de que forma os ńıveis de voltagem e temperatura afetam o tempo médio de resistência dos vidros. Faça também uma análise de diagnóstico. 32. (Ryan e Joiner, 1994, p. 299). No arquivo trees.dat é apresentado um conjunto de dados que tem sido analisado sob diversos pontos de vista por vários pesquisadores (ver, por exemplo, Jørgensen, 1989). As variáveis observadas são o diâmetro (d), a altura (h) e o volume (v) de uma amostra de 31 cerejeiras numa floresta do estado da Pensilvânia, EUA. A relação entre diâmetro, altura e volume de uma árvore depende da forma da mesma e pode-se considerar duas possibilidades v = 1 4 πd2h para forma ciĺındrica e v = 1 12 πd2h para forma cônica. Em ambos os casos a relação entre logv, logd e logh é dada por logv = a + blogd+ clogh. Supor inicialmente o modelo linear v = α+ βd+ γh+ , em que  ∼ N(0, σ2). Faça uma análise de diagnóstico e verifique se é posśıvel melhorar o modelo, por exemplo incluindo algum termo quadrático. 33. (Neter et al., 1996, p. 613). Os dados do arquivo store.dat referem-se a uma amostragem feita por uma determinada loja com seus clientes, que foram divididos Modelos Lineares Generalizados 83 segundo 110 áreas da cidade onde a loja está instalada. Para cada área foram observadas as seguintes variáveis: (i) número de clientes da área que frequentaram a loja num determinado peŕıodo, (ii) número de domićılios, (iii) renda média anual por domićılio (em US$), (iv) idade média dos domićılios (em anos), (v) distância entre a área e o concorrente mais próximo (em milhas) e (vi) distância entre a área e a loja (em milhas). Proponha um modelo log-linear de Poisson para explicar a primeira variável, dadas as demais. Use o método AIC para selecionar as variáveis explicativas. Interprete o modelo ajustado através de razões de médias. Faça uma análise de diagnóstico com o modelo ajustado. Interprete os resultados e trace o perfil da loja. 34. (Agresti, 1990, pgs. 122-123). Cinquenta e quatro indiv́ıduos considerados idosos são submetidos a um exame psiquiátrico para avaliar a ocorrência ou não de sin- toma de caduquice. Acredita-se que o escore obtido num exame psicológico feito previamente esteja associado com a ocorrência ou não do sintoma. Os dados são apresentados abaixo (score: escala no exame psicológico e resp: ocorrência (=1) ou não ocorrência (=0) do sintoma). Score Resp Score Resp Score Resp Score Resp Score Resp 9 1 7 1 7 0 17 0 13 0 13 1 5 1 16 0 14 0 13 0 6 1 14 1 9 0 19 0 9 0 8 1 13 0 9 0 9 0 15 0 10 1 16 0 11 0 11 0 10 0 4 1 10 0 13 0 14 0 11 0 14 1 12 0 15 0 10 0 12 0 8 1 11 0 13 0 16 0 4 0 11 1 14 0 10 0 10 0 14 0 7 1 15 0 11 0 16 0 20 0 9 1 18 0 6 0 14 0 (i) Ajustar um modelo loǵıstico para explicar a probabilidade de ocorrência do sintoma em função do escore. Interpretar os resultados. (ii) Faça os gráficos de tDi , tGi, t 2 Si e LDi contra os valores ajustados. Construa envelopes com os reśıduos tDi e tGi. Interprete os gráficos e identifique os pontos discrepantes. 84 Caṕıtulo 2 Modelos para Dados Binários 87 Como em geral a porcentagem de indiv́ıduos doentes é muito menor do que a por- centagem de não-doentes, é bastante razoável num estudo cujo objetivo é avaliar a asso- ciação entre algum fator particular e uma certa doença, que a quantidade de doentes na amostra seja a maior posśıvel. Assim, a amostragem retrospectiva, em que os indiv́ıduos são escolhidos separadamente nos estratos D e D̄, pode ser mais conveniente do que os demais procedimentos amostrais. Um cuidado, entretanto, deve-se ter nesses estudos. É importante que os doentes (casos) sejam comparáveis aos não-doentes (controles) se- gundo outros fatores (fatores potenciais de confundimento), possivelmente associados com a doença. Nos estudos prospectivos, em que a amostragem é feita nos estratos A e B, esse tipo de problema pode ser controlado, embora em geral seja necessário um longo peŕıodo até a obtenção de um número suficiente de doentes para uma análise estat́ıstica mais representativa. Como as inferências para os estudos retrospectivos e prospectivos são idênticas, tratare- mos apenas o caso retrospectivo. Assim, assumimos que no estrato D são amostrados n1 indiv́ıduos e no estrado D̄ são amostrados n2 indiv́ıduos. O número observado de indiv́ıduos com presença de A nos estratos D e D̄ será denotado por y1 e y2, respectiva- mente. Os dados resultantes dessa amostragem podem ser resumidos conforme a tabela abaixo. Fator Doença A B Total D y1 n1 − y1 n1 D̄ y2 n2 − y2 n2 Discutimos nas seções seguintes a abordagem clássica para analisar a tabela acima. 2.2.2 Modelo probabiĺıstico não-condicional Denotaremos por Y1 e Y2 o número de indiv́ıduos com presença de A nos estratos D e D̄, respectivamente. Será também assumido que essas variáveis são binomiais independentes de parâmetros (n1, π1) e (n2, π2), respectivamente. Logo, a função de probabilidades conjunta de (Y1, Y2) fica dada por f(y;π) = ( n1 y1 )( n2 y2 ) πy11 π y2 2 (1 − π1)n1−y1(1 − π2)n2−y2 , (2.3) em que y = (y1, y2) T e π = (π1, π2) T . Seguindo a notação da seção anterior, temos que π1 = P1/(P1 + P3), 1 − π1 = P3/(P1 + P3), π2 = P2/(P2 + P4) e 1 − π2 = P4/(P2 + P4). 88 Caṕıtulo 2 Assim, mostra-se que ψ = P1P4 P3P2 = π1(1 − π2) π2(1 − π1) , e consequentemente que π1 = π2ψ/{π2ψ + 1 − π2}. A expressão (2.3) pode então ser expressa apenas em função de (ψ, π2), f(y;π) ∝ exp { y1logψ + (y1 + y2)log ( π2 1 − π2 )} (1 − π2)n {ψπ2 + 1 − π2}n1 , (2.4) em que n = n1 + n2. As estimativas de máxima verossimilhança de π1 e π2 são dadas por π̃1 = y1/n1 e π̃2 = y2/n2, respectivamente. Logo, a estimativa de m.v. não-condicional de ψ fica ψ̃ = y1(n2 − y2)/y2(n1 − y1). Note que E(ψ̃) = ∞, o que impossibilita qualquer tipo de inferência para pequenas amostras. Por outro lado, para n1 e n2 grandes, ψ̃ segue uma distribuição normal de média ψ e variância assintótica VarA(ψ̃) = ψ 2 { 1 n1π1(1 − π1) + 1 n2π2(1 − π2) } . Formalmente, podemos dizer que sob condições gerais de regularidade e assumindo que n1 n → a > 0, quando n→ ∞, vale o resultado assintótico √ n(ψ̃ − ψ) →d N(0,VI(ψ)), em que VI(ψ) = ψ 2{1/aπ1(1 − π1) + 1/(1− a)π2(1− π2)}. A variância assintótica VI(ψ) é consistentemente estimada por nVarA(ψ̃). Alguns autores preferem trabalhar com logψ em vez de ψ. Assim, podemos mostrar, sob condições gerais de regularidade, que a estimativa não-condicional logψ̃ segue para grandes amostras uma distribuição normal de média logψ e variância assintótica VarA(logψ̃) = {1/n1π1(1 − π1) + 1/n2π2(1 − π2)}. Isso é equivalente a dizer que √ n(logψ̃ − logψ) →d N(0, ψ−2VI(ψ)). Esse resultado será útil na construção de intervalos de confiança para ψ. 2.2.3 Modelo probabiĺıstico condicional Devido aos problemas inferenciais com o modelo não-condicional para pequenas amostras, a utilização de um modelo condicional, cuja construção será discutida a seguir, tem sido a solução encontrada sob o ponto de vista clássico para fazer inferências a respeito de ψ. Assim, aplicando o teorema da fatorização para a função de probabilidades (2.4), mostra-se que o conjunto de estat́ısticas (Y1, Y1 + Y2) é suficiente minimal para o vetor Modelos para Dados Binários 89 de parâmetros [logψ, log{π2/(1 − π2)}]. Logo, a distribuição de (Y1, Y2) condicionada a Y1 + Y2 = m, deverá resultar numa função de probabilidades que depende apenas do parâmetro de interese ψ. Essa distribuição resultante (vide Cornfield, 1956), tem sido largamente utilizada em pequenas amostras. Alguns autores questionam, entretanto, o procedimento adotado, uma vez que a estat́ıstica Y1 + Y2 não é ancilar para ψ; isto é, contém informações a respeito do parâmetro ψ. O condicionamento de (Y1, Y2) em Y1 + Y2 = m produz o modelo caracterizado pela famı́lia de distribuições hipergeométricas não-centrais, definida por f(y1|m;ψ) = ( n1 y1 )( n2 m−y1 ) ψy1 ∑ t ( n1 t )( n2 m−t ) ψt , (2.5) em que 0 < ψ < ∞ e t varia de max(0, m − n2) a min(n1, m). Em particular, quando ψ = 1, a expressão (2.5) reduz-se à conhecida distribuição hipergeométrica central, dada por f(y1|m;ψ = 1) = ( n1 y1 )( n2 m−y1 ) ( n1+n2 m ) , cuja média e variância são, respectivamente, E(1) = E(Y1|m;ψ = 1) = mn1 n e V(1) = Var(Y1|m;ψ = 1) = n1n2(n−m)m n2(n− 1) . Para o modelo condicional (2.5) o logaritmo da função de verossimilhança fica dado por L(ψ) ∝ y1logψ − log { ∑ t ( n1 t )( n2 m− t ) ψt } . Denotaremos por ψ̂ a estimativa de m.v. condicional. Essa estimativa pode ser expressa como a solução positiva da equação y1 = E(Y1|m; ψ̂). Note que o momento de ordem r da distribuição condicional, E(Y r1 |m;ψ), é dado por E(Y r1 |m;ψ) = Pr(ψ)/P0(ψ), em que Pr(ψ) = ∑ t t r ( n1 t )( n2 m−t ) ψt, r = 1, 2, . . . e P0(ψ) = ∑ t ( n1 t )( n2 m−t ) ψt. Assim, a equação de máxima verossimilhança para obter ψ̂, fica reescrita na forma y1 − P1(ψ̂) P0(ψ̂) = 0. (2.6) Com o aumento de n1, n2, m e n − m, torna-se impraticável obter ψ̂ através de (2.6), uma vez que essa equação contém polinômios em ψ̂ de grau bastante elevado. Uma sáıda,