Docsity
Docsity

Prepare-se para as provas
Prepare-se para as provas

Estude fácil! Tem muito documento disponível na Docsity


Ganhe pontos para baixar
Ganhe pontos para baixar

Ganhe pontos ajudando outros esrudantes ou compre um plano Premium


Guias e Dicas
Guias e Dicas

Modelos Lineares Generalizados, Notas de estudo de Estatística

Modelos lineares para dados não normais

Tipologia: Notas de estudo

Antes de 2010

Compartilhado em 16/09/2008

walmes-zeviani-8
walmes-zeviani-8 🇧🇷

3 documentos

1 / 153

Documentos relacionados


Pré-visualização parcial do texto

Baixe Modelos Lineares Generalizados e outras Notas de estudo em PDF para Estatística, somente na Docsity! Modelos Lineares Generalizados - da teoria à prática - M. Antónia Amaral Turkman DEIO/FC e CEAUL, Universidade de Lisboa e Giovani Loiola Silva DM/IST e CMA, Universidade Técnica de Lisboa 0Este trabalho foi parcialmente financiado por FCT - PRAXIS XXI - FEDER Prefácio Os Modelos Lineares Generalizados foram introduzidos no ińıcio dos anos 70, tendo um impacto muito grande no desenvolvimen- to da estat́ıstica aplicada. No prinćıpio o seu uso esteve confinado a um grupo restrito de investigadores, dada a falta de bibliografia acesśıvel e à complexidade inicial do GLIM, primeiro software to- talmente dirigido à aplicação da metodologia. Foram precisos cerca de 20 anos para que a teoria e prática dos Modelos Lineares Ge- neralizados passasse da “cadeira” do investigador para o domı́nio público. Este avanço só foi posśıvel a partir do momento em que o software existente passou a ser um pouco mais “amigável”. Ac- tualmente, a maioria dos pacotes estat́ısticos de maior expansão já contém módulos adequados ao estudo destes modelos. Pode pois dizer-se que, neste momento, o conhecimento adequado da meto- dologia dos Modelos Lineares Generalizados é imprescind́ıvel para qualquer indiv́ıduo que utilize métodos estat́ısticos. Dáı a escolha deste tema para um novo mini-curso, a ser ministrado por ocasião do VIII Congresso Anual da Sociedade Portuguesa de Estat́ıstica. A importância dos Modelos Lineares Generalizados não é ape- nas de ı́ndole prática. Do ponto de vista teórico a sua importância advém, essencialmente, do facto de a metodologia destes modelos constituir uma abordagem unificada de muitos procedimentos es- i iv Índice 2.1.2 Função de ligação canónica e estat́ısticas su- ficientes . . . . . . . . . . . . . . . . . . . . . 32 2.2 Estimação dos Parâmetros do Modelo . . . . . . . . . 36 2.2.1 Método iterativo de mı́nimos quadrados pon- derados . . . . . . . . . . . . . . . . . . . . . 36 2.2.2 Estimação do parâmetro de dispersão . . . . . 39 2.2.3 Propriedades assintóticas dos estimadores de máxima verosimilhança . . . . . . . . . . . . . 41 2.2.4 Existência e unicidade dos EMV . . . . . . . . 44 2.3 Testes de Hipóteses . . . . . . . . . . . . . . . . . . . 45 2.3.1 Teste de Wald . . . . . . . . . . . . . . . . . . 48 2.3.2 Teste de razão de verosimilhanças . . . . . . . 50 2.3.3 Estat́ıstica de Rao . . . . . . . . . . . . . . . 51 2.4 Quasi-verosimilhança . . . . . . . . . . . . . . . . . . 53 3 Selecção e Validação de Modelos 59 3.1 Qualidade de Ajustamento . . . . . . . . . . . . . . . 61 3.1.1 Função desvio . . . . . . . . . . . . . . . . . . 61 3.1.2 Estat́ıstica de Pearson generalizada . . . . . . 66 3.2 Selecção de Modelos . . . . . . . . . . . . . . . . . . 67 3.3 Análise de Reśıduos . . . . . . . . . . . . . . . . . . . 72 3.3.1 Matriz de projecção generalizada . . . . . . . 73 3.3.2 Definições de reśıduos . . . . . . . . . . . . . 74 3.3.3 Análise informal dos reśıduos . . . . . . . . . 79 3.4 Observações Discordantes . . . . . . . . . . . . . . . 84 Índice v 3.4.1 Medida de repercussão . . . . . . . . . . . . . 85 3.4.2 Medida de influência . . . . . . . . . . . . . . 86 3.4.3 Medida de consistência . . . . . . . . . . . . . 87 4 Aplicações I: Modelos Discretos 91 4.1 Modelos de Regressão Loǵıstica . . . . . . . . . . . . 92 4.1.1 Selecção do modelo loǵıstico . . . . . . . . . . 94 4.1.2 Avaliação e interpretação do modelo seleccio- nado . . . . . . . . . . . . . . . . . . . . . . . 98 4.2 Modelos de Dose-resposta . . . . . . . . . . . . . . . 102 4.3 Modelos Log-lineares . . . . . . . . . . . . . . . . . . 107 5 Aplicações II: Modelos Cont́ınuos 113 5.1 Modelos de Regressão Gama . . . . . . . . . . . . . . 114 5.2 Modelos de Sobrevivência . . . . . . . . . . . . . . . 123 A Programas do S-plus 127 A.1 Exemplo 4.1 . . . . . . . . . . . . . . . . . . . . . . . 127 A.2 Exemplo 5.1 . . . . . . . . . . . . . . . . . . . . . . . 129 B Programas do GLIM 133 B.1 Exemplo 4.2 . . . . . . . . . . . . . . . . . . . . . . . 133 B.2 Exemplo 4.3 . . . . . . . . . . . . . . . . . . . . . . . 136 B.3 Exemplo 5.2 . . . . . . . . . . . . . . . . . . . . . . . 139 Caṕıtulo 1 Introdução aos Modelos Lineares Generalizados Em muitos estudos estat́ısticos, quer sejam de natureza expe- rimental ou observacional, somos confrontados com problemas em que o objectivo principal é o de estudar a relação entre variáveis, ou mais particularmente, analisar a influência que uma ou mais va- riáveis (explicativas), medidas em indiv́ıduos ou objectos, têm sobre uma variável de interesse a que damos o nome de variável resposta. O modo como, em geral, o estat́ıstico aborda tal problema é através do estudo de um modelo de regressão que relacione essa variável de interesse com as variáveis ditas explicativas. O modelo linear normal, “criado” no ińıcio do século XIX por Legendre e Gauss, dominou a modelação estat́ıstica até meados do século XX, embora vários modelos não lineares ou não normais tenham entretanto sido desenvolvidos para fazer face a situações que não eram adequadamente explicadas pelo modelo linear nor- mal. São exemplo disso, tal como referem McCullagh and Nelder 1 4 1. Introdução aos Modelos Lineares Generalizados experimentais, sendo as componentes Yi do vector aleatório Y = (Y1, . . . , Yn) T independentes. Irá ser útil, no desenrolar da teoria, a representação dos dados em (1.1) na forma matricial y =   y1 y2 ... yn   X =   x11 x12 ... x1k x21 x22 ... x2k ... ... ... xn1 xn2 ... xnk   . (1.2) Em muitas situações práticas, principalmente quando as variáveis explicativas são de natureza qualitativa, há muitos indiv́ıduos na amostra que partilham do mesmo vector de covariáveis, significando isto que a matriz X tem vários grupos de linhas idênticas. Assim pode ter interesse em apresentar os dados, não desagrupados como em (1.2), mas de uma forma agrupada. Suponhamos então que podemos associar os indiv́ıduos em g grupos distintos, de tal modo que os nj indiv́ıduos do grupo j (j = 1, ..., g com ∑g j=1 nj = n) partilhem do mesmo vector de cova- riáveis, digamos xj = (xj1, ..., xjk) T . Os dados passarão a ser então representados por y =   y1 y2 ... yg   Xg =   x11 x12 ... x1k x21 x22 ... x2k ... ... ... xg1 xg2 ... xgk   , (1.3) onde yj, j = 1, ..., g, representa a média das variáveis respostas dos indiv́ıduos que pertencem ao j-ésimo grupo e não existem linhas idênticas em Xg. O agrupamento dos dados é particularmente importante, e tem significado especial, em situações em que as covariáveis são todas de natureza qualitativa. 1.2. A Famı́lia Exponencial 5 1.2 A Famı́lia Exponencial Como já se disse anteriormente, os modelos lineares generalizados pressupõem que a variável resposta tenha uma distribuição perten- cente a uma famı́lia particular, a famı́lia exponencial. A definição que vamos aqui apresentar é a adequada para os modelos para a va- riável resposta que interessa considerar no âmbito dos MLG. Veja- se, e.g., Cox and Hinkley (1974), para uma definição mais geral de famı́lia exponencial k-paramétrica e suas propriedades. Definição 1 (Famı́lia exponencial) Diz-se que uma variável aleatória Y tem distribuição pertencen- te à famı́lia exponencial de dispersão (ou simplesmente famı́lia ex- ponencial) se a sua função densidade de probabilidade (f.d.p.) ou função massa de probabilidade (f.m.p.) se puder escrever na forma f(y|θ, φ) = exp {yθ − b(θ) a(φ) + c(y, φ) } , (1.4) onde θ e φ são parâmetros escalares, a(·), b(·) e c(·, ·) são funções reais conhecidas. ♦ Na definição anterior, θ é a forma canónica do parâmetro de loca- lização e φ é um parâmetro de dispersão suposto, em geral, conheci- do. Neste caso a distribuição descrita em (1.4) faz parte da famı́lia exponencial uniparamétrica. Quando o parâmetro φ é desconheci- do a distribuição pode ou não fazer parte da famı́lia exponencial bi-paramétrica, tal como é geralmente definida (veja, e.g., Cox and Hinkley, 1974). Admite-se, ainda, que a função b(·) é diferenciável e que o suporte da distribuição não depende dos parâmetros. Neste caso prova-se que a famı́lia em consideração obedece às condições 6 1. Introdução aos Modelos Lineares Generalizados habituais de regularidade1. 1.2.1 Valor médio e variância Seja ℓ(θ;φ, y) = ln(f(y|θ, φ)). Defina-se a função score S(θ) = ∂ℓ(θ;φ, Y ) ∂θ . (1.5) Sabe-se que para famı́lias regulares se tem E(S(θ)) = 0 E(S2(θ)) = E [( ∂ℓ(θ;φ,Y ) ∂θ )2] = −E [ ∂2ℓ(θ;φ,Y ) ∂θ2 ] (1.6) e portanto como, no caso em que f(y|θ, φ) é dado por (1.4), ℓ(θ;φ, y) = yθ − b(θ) a(φ) + c(y, φ), obtém-se S(θ) = Y − b′(θ) a(φ) ∂S(θ) ∂θ = −b ′′(θ) a(φ) , (1.7) onde b′(θ) = ∂b(θ) ∂θ e b′′(θ) = ∂ 2b(θ) ∂θ2 . Assim de (1.6) e (1.7) E(Y ) = µ = a(φ)E(S(θ)) + b′(θ) = b′(θ) (1.8) var(Y ) = a2(φ)var(S(θ)) = a2(φ) b′′(θ) a(φ) = a(φ)b′′(θ). (1.9) Vê-se assim que a variância de Y é o produto de duas funções; uma, b′′(θ), que depende apenas do parâmetro canónico θ (e por- tanto do valor médio µ), a que se dá o nome de função de variância 1Para um estudo de condições de regularidade necessárias no desenvolvimen- to do estudo que se vai fazer, deve consultar-se um livro avançado de Estat́ıstica. Aconselha-se, por exemplo, Sen and Singer (1993). 1.2. A Famı́lia Exponencial 9 De (1.8) e (1.9) obtém-se directamente E(Y ) = b′(θ) = π, var(Y ) = b′′(θ)a(φ) = π(1 − π) m . O parâmetro canónico é a função logit, ln ( π 1−π ) . Exemplo 1.3 Gama Se Y tem distribuição gama com parâmetros de forma ν e de escala ν/µ (Y ∼ Ga(ν, ν/µ)), a sua f.d.p. é f(y|ν, µ) = 1 Γ(ν) (ν µ )ν yν−1 exp ( − ν µ y ) = exp { ν ( − y µ − lnµ ) + (ν − 1) ln y − ln Γ(ν) + ν ln ν } = exp{ν(θy + ln(−θ)) + (ν − 1) ln y − ln Γ(ν) + ν ln ν} com y > 0 e θ = − 1 µ . Temos novamente uma f.d.p. da forma (1.10) com θ = −1 µ , b(θ) = − ln(−θ), c(y, φ) = (ν − 1) ln y + ν ln ν − ln Γ(ν), b′(θ) = −1 θ , b′′(θ) = V (µ) = 1 θ2 = µ2, a(φ) = φ ω , φ = 1 ν , ω = 1. Novamente de (1.8) e (1.9) tem-se que E(Y ) = −1 θ = µ, var(Y ) = µ2 ν . A função de variância é, neste caso, V (µ) = µ2 e o parâmetro de dispersão é 1 ν . Na tabela 1.1 apresentamos uma lista das principais distribuições que pertencem à famı́lia exponencial com a respectiva caracteri- zação. 10 1. In tro d u ção aos M o d elos L in eares G en eralizad os Tabela 1.1: Algumas distribuições da famı́lia exponencial. distribuição normal binomial Poisson gama gaussiana inversa Notação N(µ, σ2) B(m, π)/m P (λ) Ga(ν, ν µ ) IG(µ, σ2) Suporte (−∞,+∞) {0, 1 m , ..., 1} {0, 1, ...} (0,+∞) (0,+∞) θ µ ln ( π 1−π ) lnλ − 1 µ − 1 2µ2 a(φ) σ2 1 m 1 1 ν σ2 φ σ2 1 1 1 ν σ2 ω 1 m 1 1 1 c(y, φ) −1 2 (y 2 φ + ln(2πφ)) ln ( m my ) − ln y! ν ln ν − ln Γ(ν) −1 2 {ln(2πφy3) +(ν − 1) ln y + 1 yφ } b(θ) θ 2 2 ln(1 + eθ) eθ − ln(−θ) −(−2θ)1/2 b′(θ) = E(Y ) θ π = e θ 1+eθ λ = eθ µ = −1 θ µ = (−2θ)−1/2 b′′(θ) = V (µ) 1 π(1 − π) λ µ2 µ3 var(Y ) σ2 π(1−π) m λ µ 2 ν µ3σ2 1.3. Descrição do Modelo Linear Generalizado 11 1.3 Descrição do Modelo Linear Gene- ralizado Os modelos lineares generalizados são uma extensão do modelo linear clássico Y = Zβ + ε, onde Z é uma matriz de dimensão n × p de especificação do mo- delo (em geral a matriz de covariáveis X com um primeiro vector unitário), associada a um vector β = (β1, . . . , βp) T de parâmetros, e ε é um vector de erros aleatórios com distribuição que se supõe Nn(0, σ 2I). Estas hipóteses implicam obviamente que E(Y|Z) = µ com µ = Zβ, ou seja, o valor esperado da variável resposta é uma função linear das covariáveis. A extensão mencionada é feita em duas direcções. Por um la- do, a distribuição considerada não tem de ser normal, podendo ser qualquer distribuição da famı́lia exponencial; por outro lado, embo- ra se mantenha a estrutura de linearidade, a função que relaciona o valor esperado e o vector de covariáveis pode ser qualquer função diferenciável. Assim os MLG são caracterizados pela seguinte estrutura: 1. Componente aleatória Dado o vector de covariáveis xi as variáveis Yi são (condi- cionalmente) independentes com distribuição pertencente à famı́lia exponencial da forma (1.4) ou (1.10), com E(Yi|xi) = µi = b ′(θi) para i = 1, .., n e, possivelmente, um parâmetro de dispersão φ não dependente de i. 14 1. Introdução aos Modelos Lineares Generalizados 1.4 Exemplos de Modelos Lineares Ge- neralizados Nesta secção remos apresentar alguns dos modelos lineares gene- ralizados mais comuns nas aplicações. Convém distinguir modelos para três tipos de respostas: (i) de natureza cont́ınua, (ii) de na- tureza dicotómica, ou na forma de proporções e (iii) na forma de contagens. Por essa razão apresentamos os exemplos agrupados de acordo com essa divisão. 1.4.1 Modelos para respostas cont́ınuas Modelo normal Já se referiu anteriormente que os MLG correspondem a uma generalização do modelo de regressão linear. Com efeito, se tivermos n respostas independentes Yi ∼ N(µi, σ2), i = 1..., n onde µi = z T i β = p∑ j=1 zijβj, o modelo considerado é um modelo linear generalizado, dado que: • as variáveis resposta são independentes, • a distribuição é da forma (1.10), com θi = zTi β, φ = σ2 e ωi = 1. • o valor esperado µi está relacionado com o preditor linear ηi = zTi β através da relação µi = ηi, • a função de ligação é a função identidade, que é, neste caso a função de ligação canónica. 1.4. Exemplos de Modelos Lineares Generalizados 15 Para este modelo podemos escrever Yi = z T i β + εi, i = 1, ..., n onde os εi são i.i.d. N(0, σ 2). Note-se ainda que a formulação apresentada inclui facilmente o caso especial em que Yi ∼ N(µi, σ2i ), com σ2i = σ 2 ωi , onde ωi é um peso conhecido associado à i-ésima observação. O modelo normal (modelo linear clássico), introduzido atrás, pressupõe, como se sabe, que a variância das resposta seja cons- tante. Contudo, na prática, surgem por vezes situações de variáveis resposta de natureza cont́ınua, em que a variância não é constante. Uma transformação que se usa com frequência para estabilizar a va- riância, é a transformação logaŕıtmica, a qual é posśıvel se as respos- tas forem positivas. Admitindo então que o logaritmo das respostas Yi segue uma distribuição normal, pode considerar-se um modelo de regressão linear clássico para o logaritmo das respostas. Neste caso, ter-se-á a relação ηi = E{ln(Yi)} = zTi β. Por várias razões e, em particular, se há necessidade das conclusões serem apresentadas na escala original das respostas, então é mais conveniente assumir que lnE(Yi) = z T i β, ou seja, que E(Yi) = exp(z T i β). Se, por outro lado, assumirmos que a variância aumenta com o valor médio de modo que o coeficiente de variação se mantém constante, o modelo gama passa a ser um modelo adequado para as respostas (McCullagh and Nelder, 1989). Modelo gama Admitindo então que as respostas são variáveis aleatórias Yi ∼ Ga(ν, ν µi ) independentes, com µi = exp(z T i β) obtém-se um modelo linear generalizado (o modelo de regressão gama), visto que: • as variáveis resposta são independentes, 16 1. Introdução aos Modelos Lineares Generalizados • a distribuição é da forma (1.10), com θi = − 1 exp(zTi β) , φ = 1 ν e ωi = 1, • o valor esperado µi está relacionado com o preditor linear ηi = zTi β através da relação µi = exp(ηi), • a função de ligação é a função logaŕıtmica. Neste caso podemos escrever o modelo na forma Yi = exp(z T i β) ǫi, i = 1, ..., n com ǫi i.i.d. Ga(ν, ν). Um exemplo de aplicação deste modelo é feito no caṕıtulo de aplicações práticas. Note-se que a função de ligação considerada não é a função de ligação canónica. A função de ligação canónica obtém-se quando ηi = θi, o que neste caso corresponde a ter − 1µi = z T i β. A função de ligação canónica é então a função rećıproca. Dado que µi > 0, a utilização do modelo gama com função de ligação canónica implica que se têm de impor restrições aos valores posśıveis para os parâmetros βj da estrutura linear. Nelder (1966) apresenta um exemplo de aplicação de regressão gama com função de ligação canónica. Modelo gaussiano inverso Suponhamos que Yi ∼ IG(µi, σ2), isto é f(yi|µi, σ) = 1√ 2πσ2y3i exp { − (yi − µi) 2 2µ2iσ 2yi } , yi > 0 e µi = (exp{zTi β}) 1 2 , para i = 1, ..., n. Neste caso obtém-se, como facilmente se verifica, um modelo li- near generalizado com função de ligação canónica. 1.4. Exemplos de Modelos Lineares Generalizados 19 O modelo linear generalizado, obtido pela associação do modelo binomial para as respostas, com a função de ligação complementar log-log conduz ao modelo de regressão complementar log-log. A utilização de uma ou outra função de ligação, e consequen- temente, a escolha do modelo de regressão a utilizar, depende da situação em causa. Em geral, a adaptabilidade dos modelos probit e loǵıstico é bastante semelhante, já que as funções correspondentes não se afastam muito uma da outra após um ajustamento adequa- do dos correspondentes preditores lineares. O modelo complementar log-log pode dar respostas diferentes já que a função complementar log-log, mesmo após o ajustamento do preditor linear η, se distancia das anteriores, tendo um crescimento mais abrupto (ver, Fahrmeir and Tutz, 1994, pg. 27). A função de ligação complementar log-log é mais utilizada para análise de dados sobre incidência de doenças. Nos caṕıtulos dedicados a aplicações práticas apresentaremos exemplos de modelos de regressão loǵıstico e probit. Relação com modelos lineares latentes Variáveis aleatórias binárias podem ser consideradas como re- sultantes da dicotomização de variáveis aleatórias cont́ınuas. Com efeito, se Z for uma variável aleatória cont́ınua com função de dis- tribuição FZ(·), e se, em vez de se observar Z, se observar apenas se Z está acima ou abaixo um determinado valor cŕıtico r, a variável aleatória Y = { 1, se Z ≤ r 0, se Z > r, é uma variável aleatória de Bernoulli com probabilidade de sucesso π = P (Y = 1) = FZ(r). Z será assim uma variável aleatória latente não observada. 20 1. Introdução aos Modelos Lineares Generalizados Os modelos para dados binários apresentados podem, deste mo- do, ser explicados como resultantes de modelos lineares latentes. Com efeito, se Z = zT α + σε, onde σ é um parâmetro de escala desconhecido, αT = (α1, ..., αp), z = (1, z2, ..., zp) T é um vector de especificação e ε tiver distribuição F (·) (e.g., loǵıstica, normal reduzida, ou de extremos) então π = P (Y = 1) = P (Z ≤ r) = P (zT α + σε ≤ r) = P (ε ≤ r − α1 σ − p∑ j=2 zj αj σ ) = F (zT β), com β1 = r−α1 σ e βj = −αjσ , j = 2, ..., p. A abordagem de um modelo de dados binários, na perspectiva de um modelo linear latente, permite a interpretação dos parâmetros β’s em função desse modelo; no entanto, se σ for desconhecido os efeitos das covariáveis no modelo linear (αj , j = 2, .., p) só são conhe- cidos a menos do factor 1 σ e portanto só se pode atribuir significado aos valores relativos dos parâmetros (e.g., β2/β3) e não aos seus valores absolutos. Dados agrupados Até aqui estivemos a supor que os dados se encontram numa for- ma não agrupada. Em muitas situações de interesse acontece, como já foi referido na secção 1.1, que vários indiv́ıduos, ou unidades expe- rimentais, partilham do mesmo vector de covariáveis, podendo então os indiv́ıduos ser agrupados de acordo com os diferentes padrões 1.4. Exemplos de Modelos Lineares Generalizados 21 de covariáveis. Neste caso consideramos como variável resposta a frequência relativa de sucessos do grupo, i.e., Y i = 1 ni ∑ni j=1 Yij, onde ni é o número de indiv́ıduos no i-ésimo grupo. Como as frequências absolutas têm distribuição B(ni, πi), as frequências relativas são dis- tribúıdas de acordo com B(ni, πi)/ni, i.e. P (Y i = yi|πi) = ( ni niyi ) π niyi i (1 − πi)ni−niyi yi = 0, 1 ni , ..., 1. Repare-se que, como ainda se tem E(Y i) = πi, podem ainda con- siderar-se as mesmas funções de ligação que se consideraram para o caso em que as respostas são binárias. A mesma metodologia pode também ser aplicada no caso em que as respostas individuais Yi (não agrupadas) são B(ni, πi) bastando para tal considerar como resposta a variável Yi ni . Sobredispersão ou Extra variação binomial Um fenómeno que ocorre com frequência nas aplicações é as res- postas apresentarem uma variância superior à variância explicada pelo modelo binomial. Este fenómeno, denominado de sobredis- persão ou extra variação binomial, pode ser devido ao facto de exis- tir heterogeneidade entre os indiv́ıduos não explicada pelas cova- riáveis, ou pelo facto de haver correlação entre as respostas. Esta última situação acontece quando, por exemplo, as respostas corres- pondem a indiv́ıduos da mesma famı́lia, ou a indiv́ıduos que co- mungam dos mesmos factores ambientais, formando assim grupos naturais, embora a heterogeneidade não explicada também produza correlação entre as respostas. Este problema pode ser resolvido se se introduzir um parâmetro φ > 1 de sobredispersão de tal modo 24 1. Introdução aos Modelos Lineares Generalizados • Ajustamento dos modelos, • Selecção e validação dos modelos. Na formulação do modelo há que entrar em consideração com (i) escolha da distribuição para a variável resposta. Para isso há necessidade de examinar cuidadosamente os dados; por exemplo, a distribuição gama e normal inversa são apropriadas para modelar dados de natureza cont́ınua e que mostram assimetrias; por vezes pode haver necessidade de transformar previamente os dados, etc. Assim, uma análise preliminar dos dados, é fundamental para que se possa fazer uma escolha adequada da famı́lia de distribuições a considerar. (ii) escolha das covariáveis e formulação apropriada da matriz de especificação. Aqui há que entrar em linha de conta com o proble- ma espećıfico em estudo e, muito particularmente, ter em atenção a codificação apropriada das variáveis de natureza qualitativa, crian- do nomeadamente, caso se revele necessário, variáveis mudas, para correctamente definir variáveis de natureza policotómica. (iii) escolha da função de ligação. A escolha de uma função de li- gação compat́ıvel com a distribuição do erro proposto para os dados deve resultar de uma combinação de considerações a priori sobre o problema em causa, exame intensivo dos dados, facilidade de inter- pretação do modelo, etc. Existem funções de ligação que produzem propriedades estat́ısticas desejadas para o modelo, como iremos ver na secção 2.1.2, mas a conveniência matemática por si só não deve determinar a escolha da função de ligação. A fase do ajustamento do modelo (ou modelos) passa pela estimação dos parâmetros do modelo, isto é, pela estimação dos co- eficientes β’s associados às covariáveis, e do parâmetro de dispersão 1.5. Metodologia dos Modelos Lineares Generalizados 25 φ caso ele esteja presente. É importante também nesta fase estimar parâmetros que representam medidas da adequabilidade dos valores estimados, obter intervalos de confiança e realizar testes de bonda- de de ajustamento. O problema da inferência em modelos lineares generalizados será tratado no caṕıtulo 2. Nos problemas em que a metodologia dos MLG tem cabimento, ou seja em problemas de regressão, há em geral um número con- siderável de posśıveis variáveis explicativas. A fase de selecção e validação dos modelos tem por objectivo encontrar submodelos com um número moderado de parâmetros que ainda sejam adequa- dos aos dados, detectar discrepâncias entre os dados e os valores preditos, averiguar da existência de outliers ou/e observações influ- entes, etc. Na selecção do melhor modelo para explicar o problema em estudo, devem ainda ser ponderados três factores: adequabili- dade, parcimónia e interpretação. Um bom modelo é aquele que consegue atingir um equiĺıbrio entre esses três factores. O problema da selecção e validação dos modelos será tratado no caṕıtulo 3. Caṕıtulo 2 Inferência De modo a poder aplicar a metodologia dos modelos lineares generalizados a um conjunto de dados há necessidade, após a for- mulação do modelo que se pensa adequado, de proceder à realização de inferências sobre esse modelo. A inferência com MLG é, essenci- almente, baseada na verosimilhança. Com efeito, não só o método da máxima verosimilhança é o método de eleição para estimar os parâmetros de regressão, como também testes de hipóteses sobre os parâmetros do modelo e de qualidade de ajustamento são, em geral, métodos baseados na verosimilhança. Os métodos inferenciais que vamos discutir neste caṕıtulo pres- supõem que o modelo está completamente e correctamente espe- cificado, de acordo com a formulação apresentada na secção 1.3. Contudo, por vezes, essa suposição não é realista. É o caso, por exemplo, em que se verifica que há sobredispersão num modelo de Poisson ou binomial e que portanto há necessidade de alterar a va- riância através da introdução de um parâmetro de sobredispersão, como se viu na secção 1.4. Nessa situação já não é posśıvel especi- 27 30 2. Inferência a chamar de log-verosimilhança) é dado por lnL(β) = ℓ(β) = n∑ i=1 ωi(yiθi − b(θi) φ + c(yi, φ, ωi) = n∑ i=1 ℓi(β), (2.4) onde ℓi(β) = ωi(yiθi − b(θi) φ + c(yi, φ, ωi) (2.5) é a contribuição de cada observação yi para a verosimilhança. Admitindo que se verificam certas condições de regularidade (ver, e.g., Sen and Singer, 1993) os estimadores de máxima verosimi- lhança para β são obtidos como solução do sistema de equações de verosimilhança ∂ℓ(β) ∂βj = n∑ i=1 ∂ℓi(β) ∂βj = 0, j = 1, ..., p. Para obter estas equações escrevemos: ∂ℓi(β) ∂βj = ∂ℓi(θi) ∂θi ∂θi(µi) ∂µi ∂µi(ηi) ∂ηi ∂ηi(β) ∂βj sendo ∂ℓi(θi) ∂θi = ωi(yi − b′(θi)) φ = ωi(yi − µi) φ , ∂µi ∂θi = b′′(θi) = ωivar(Yi) φ , ∂ηi(β) ∂βj = zij . Assim ∂ℓi(β) ∂βj = ωi(yi − µi) φ φ ωivar(Yi) ∂µi ∂ηi zij (2.6) 2.1. Estimação 31 e as equações de verosimilhança para β são n∑ i=1 (yi − µi)zij var(Yi) ∂µi ∂ηi = 0, j = 1, ..., p. (2.7) A função score, tal como já foi definida em (1.5), é o vector p- dimensional s(β) = ∂ℓ(β) ∂β = n∑ i=1 si(β), onde si(β) é o vector de componentes ∂ℓi(β) ∂βj obtidas em (2.6). O elemento genérico de ordem j da função score é então n∑ i=1 (yi − µi)zij var(Yi) ∂µi ∂ηi . (2.8) A matriz de covariância da função score, I(β) = E[−∂S(β) ∂β ] é co- nhecida por matriz de informação de Fisher. Para obter a matriz de informação de Fisher temos de considerar o valor esperado das segundas derivadas de ℓi(β). Tem-se, para famı́lias regulares, que −E ( ∂2ℓi ∂βj∂βk ) = E ( ∂ℓi ∂βj ∂ℓi ∂βk ) = E [((Yi − µi)zij var(Yi) ∂µi ∂ηi )((Yi − µi)zik var(Yi) ∂µi ∂ηi )] = E [(Yi − µi)2zijzik (var(Yi))2 (∂µi ∂ηi )2] = zijzik var(Yi) (∂µi ∂ηi )2 e, portanto, o elemento genérico de ordem (j, k) da matriz de infor- mação de Fisher é − n∑ i=1 E ( ∂2ℓi ∂βj∂βk ) = n∑ i=1 zijzik var(Yi) (∂µi ∂ηi )2 . (2.9) 32 2. Inferência Na forma matricial temos I(β) = ZTWZ , (2.10) onde W é a matriz diagonal de ordem n cujo i-ésimo elemento é ̟i = ( ∂µi ∂ηi )2 var(Yi) = ωi ( ∂µi ∂ηi )2 φV (µi) . (2.11) A última igualdade em (2.11) aparece devido à relação entre a função de variância V (µ) = b′′(θ) e a variância de Y , nomeada- mente var(Y ) = φ ω V (µ), referida em (1.9), na secção 1.2.1. 2.1.2 Função de ligação canónica e estat́ısticas suficientes Quando a função de ligação é a canónica, isto é, quando θi = ηi = z T i β, a verosimilhança pode escrever-se na forma L(β) = exp { n∑ i=1 ωi(z T i βyi − b(θi)) φ + n∑ i=1 c(yi, φ, ωi) } = exp { p∑ j=1 ( n∑ i=1 ωiyizij φ ) βj − n∑ i=1 ωib(θi) φ + n∑ i=1 c(yi, φ, ωi) } o que mostra, pelo teorema da factorização, que se φ for conhe- cido, a estat́ıstica suficiente mı́nima para o modelo (para o vector parâmetro β) tem dimensão p e é dada pelo vector ∑n i=1 ωiyizi = ( ∑n i=1 ωiyizi1, . . . , ∑n i=1 ωiyizip) T . Se φ for desconhecido e a famı́lia ainda se puder escrever na forma da famı́lia exponencial, então aque- le vector é uma componente da estat́ıstica suficiente mı́nima. Vejamos alguns exemplos. 2.1. Estimação 35 ou na forma matricial a ZTy = ZT µ̂. Se, como é habitual, a primeira coluna da matriz Z for um vector unitário, i.e., zi1 = 1, i = 1, ..., n, então a equação anterior implica, por exemplo, que a soma dos valores observados yi é igual à soma das estimativas dos valores esperados µi. Relativamente à matriz de informação de Fisher temos, pelo facto das segundas derivadas da log-verosimilhança ℓi ∂2ℓi ∂βi∂βk = −ωizij φ ∂µi ∂βk não dependerem das observações yi, que −E [ ∂2ℓi ∂βi∂βk ] = ωizij φ ∂µi ∂βk , ou seja a matriz de informação de Fisher I(β) coincide com −H(β), o simétrico da matriz Hessiana. Note-se contudo, que embora as funções de ligação canónica con- duzam a propriedades estat́ısticas desejáveis para o modelo, tais como, suficiência, facilidade de cálculo, unicidade das estimativas de máxima verosimilhança e, por vezes, interpretação simples, não há razão para, à partida, escolher a função de ligação canónica e nem sempre é com ela que se obtêm os melhores resultados. Por exemplo, no modelo gama, a utilização da função de ligação canónica obriga a impor restrições aos parâmetros. Como já se disse anteriormente, o aspecto importante a ter em consideração na escolha da ligação é a adaptabilidade e adequabilidade do modelo. 36 2. Inferência 2.2 Estimação dos Parâmetros do Mo- delo Os estimadores de máxima verosimilhança (EMV) de β são ob- tidos como solução das equações de verosimilhança (2.7). A solução não corresponde necessariamente a um máximo global da função ℓ(β). Contudo, em muitos modelos a função log-verosimilhança ℓ(β) é côncava, de modo que o máximo local e global coincidem. Para funções estritamente côncavas os estimadores de máxima ve- rosimilhança são mesmo únicos, quando existem. O problema da existência e unicidade dos estimadores de máxima verosimilhança será referido na secção 2.2.4. Partindo do prinćıpio que existe so- lução e que ela é única, subsiste ainda um problema com o cálculo das estimativas de máxima verosimilhança. É que as equações (2.7) não têm, em geral, solução anaĺıtica e, portanto, a sua resolução implica o recurso a métodos numéricos. Uma das razões que fez com que os MLG tivessem sucesso, foi o facto da possibilidade de usar um único algoritmo (sugerido por Nelder e Wedderburn, 1972), para resolver (2.7), havendo apenas a necessidade de proceder a pequenos ajustamentos de acordo com a distribuição de probabilidade e a função de ligação em consideração. Além disso o algoritmo proposto opera através de uma sequência de problemas de mı́nimos quadrados ponderados para os quais existem técnicas numéricas bem testadas. 2.2.1 Método iterativo de mı́nimos quadrados ponderados O esquema iterativo para a resolução numérica das equações de verosimilhança que se vai apresentar, é baseado no método de scores 2.2. Estimação dos Parâmetros do Modelo 37 de Fisher. Seja β̂ (0) uma estimativa inicial para β. O processo de scores de Fisher procede com o cálculo das sucessivas iteradas através da relação: β̂ (k+1) = β̂ (k) + [ I(β̂(k)) ]−1 s(β̂ (k) ), onde I(·)−1, a inversa (que se supõe existir) da matriz de informação de Fisher dada em (2.10) e s(·), o vector de scores, são calculados para β = β̂ (k) . A diferença existente entre este algoritmo e o algoritmo de New- ton-Raphson para resolver sistemas de equações não lineares, reside na utilização da matriz de informação de Fisher em vez da matriz Hessiana. A vantagem desta substituição deve-se ao facto de, em geral, ser mais fácil calcular a matriz de informação I, para além de ser sempre uma matriz semi-definida positiva. A expressão anterior pode ser ainda escrita na forma [ I(β̂(k)) ] β̂ (k+1) = [ I(β̂(k)) ] β̂ (k) + s(β̂ (k) ). (2.12) Atendendo a (2.8) e (2.9), o lado direito da equação (2.12) é um vector com elemento genérico de ordem l dado por: p∑ j=1 [ n∑ i=1 zijzil var(Yi) (∂µi ∂ηi )2] β (k) j + n∑ i=1 (yi − µi)zil var(Yi) ∂µi ∂ηi e, portanto, na forma matricial tem-se I(β̂(k))β̂(k+1) = ZTW (k)u(k), onde u(k) é um vector com elemento genérico u (k) i = p∑ j=1 zijβ (k) j + (yi − µ(k)i ) ∂η (k) i ∂µ (k) i = η (k) i + (yi − µ(k)i ) ∂η (k) i ∂µ (k) i (2.13) 40 2. Inferência verosimilhança, as estimativas de máxima verosimilhança para os parâmetros µi são dadas por µ̂i = h(z T i β̂), onde a função h(·) é a inversa da função de ligação. Por outro lado, como se tem que var(Yi) = b ′′(θi) φ ωi = V (µi)φ ωi , i = 1, ..., n, então E [ωi(Yi − µi)2 V (µi) ] = φ i = 1, ..., n. Pela lei fraca dos grandes números, se 1 n2 n∑ i=1 ω2iE(Yi − µi)4 [V (µi]2 −→ 0, quando n→ ∞, então 1 n n∑ i=1 ωi(Yi − µi)2 V (µi) P−→ φ. Segue-se então que se V (·) é uma função cont́ınua e µ̂i P−→ µi para todo o i, então 1 n− p n∑ i=1 ωi(Yi − µ̂i)2 V (µ̂i) P−→ φ. Assim podemos estimar φ por: φ̂ = 1 n− p n∑ i=1 ωi(yi − µ̂i)2 V (µ̂i) , (2.15) sendo φ̂ um estimador consistente de φ. Como, por outro lado se tem que, para grandes valores de n, 1 φ n∑ i=1 ωi(Yi − µ̂i)2 V (µ̂i) = n∑ i=1 ωi(Yi − µ̂i)2 var(Yi) 2.2. Estimação dos Parâmetros do Modelo 41 tem uma distribuição aproximada de um χ2 com n − p graus de liberdade, também se conclui que φ̂ é assintoticamente centrado. A Estat́ıstica do lado direito da equação (2.15) é conhecida como estat́ıstica de Pearson generalizada, a qual também é útil para julgar da qualidade de ajustamento de um modelo. Esta estimativa de φ é mais simples e produz, em geral, maior estabilidade numérica que a de máxima verosimilhança. Note-se que, no caso do modelo de regressão linear normal , ωi = 1 e V (µ̂i) = 1 e portanto φ̂ coincide com a estimativa habitual de σ2. 2.2.3 Propriedades assintóticas dos estimadores de máxima verosimilhança Na secção anterior vimos como, formulado um modelo linear ge- neralizado, podemos proceder à estimação por máxima verosimi- lhança do vector parâmetro β dos coeficientes de regressão. Para fazer inferências sobre estes parâmetros, nomeadamente obter in- tervalos de confiança e fazer testes de hipóteses, há necessidade de conhecer a distribuição de amostragem de β̂. Em geral, não é posśıvel nos MLG, obter as distribuições de amostragem exactas para os estimadores de máxima verosimilhança dos β’s. Iremos ape- lar então, para resultados conhecidos da teoria assintótica, que se verificam quando os modelos em estudo satisfazem certas condições de regularidade ; essas condições são, com efeito, satisfeitas para os MLG. Não iremos entrar em detalhes de natureza teórica, deixando ao cuidado do leitor interessado a leitura, por exemplo, de Fahrmeir and Kaufmann (1985), onde são estabelecidas, com o rigor adequa- do, condições gerais que garantem a consistência e a normalidade assintótica do estimador β̂. 42 2. Inferência Como vimos, o estimador de máxima verosimilhança, β̂, de β obtém-se como solução de s(β̂) = 0, onde s(β) é o vector score. Também é sabido que, em condições de regularidade , este vector aleatório é tal que E(S(β)) = 0, cov(S(β)) = E(S(β)S(β)T ) = I(β), onde I(β) é a matriz de informação de Fisher já definida. Por outro lado, pelo teorema do limite central, temos a garantia de que, pelo menos assintoticamente, S(β) tem uma distribuição normal multivariada, i.e., S(β) a∼ Np(0, I(β)) e que, portanto, para grandes amostras e supondo que o modelo com o vector parâmetro β especificado é verdadeiro, a estat́ıstica S(β)TI−1(β)S(β) tem uma distribuição assintótica de um qui-qua- drado, i.e., S(β)TI−1(β)S(β) a∼ χ2p. Se desenvolvermos S(β) em série de Taylor em torno de β̂ e retivermos apenas os termos de 1a ordem, obtemos a relação: S(β) ≈ S(β̂) + ∂S(β) ∂β | β=β̂ (β − β̂) (2.16) Atendendo a que S(β̂) = 0 e substituindo a matriz de informação observada −H(β̂) = ∂S(β) ∂β | β=β̂ pela matriz de informação de Fisher, isto é, fazendo −H(β̂) = I(β) em (2.16) (o que admitimos ser apro- ximadamente válido para grandes amostras), obtemos β̂ − β ≈ I−1(β)S(β) (2.17) 2.3. Testes de Hipóteses 45 cações, é saber se a verosimilhança tem um máximo na fronteira do espaço admisśıvel para o vector parâmetro, já que a existência de tal máximo pode levar a problemas de natureza computacional. Não há, no entanto, uma teoria geral sobre o problema da existên- cia e unicidade de estimadores de máxima verosimilhança para os modelos lineares generalizados, pois que nem todos os modelos têm propriedades comuns no que diz respeito a esta questão. Há, contu- do, resultados obtidos para casos particulares; Haberman (1974) dedicou o seu estudo a modelos log-lineares e binomiais; Silva- pulle (1981) apresentou condições necessárias e suficientes para a existência de estimadores de máxima verosimilhança para os mo- delos binomiais com funções de ligação gerais, dedicando particular atenção aos modelos loǵıstico e probit; Wederburn (1976) estudou condições de existência e unicidade dos estimadores de máxima vero- similhança nos modelos normal, binomial, Poisson e gama. As con- clusões desse estudo encontram-se resumidas na tabela 2.12. Mais referências podem ser encontradas em Fahrmeir and Tutz (1994). 2.3 Testes de Hipóteses A maior parte dos problemas de inferência relacionados com tes- tes de hipóteses sobre o vector parâmetro β, podem ser formulados em termos de hipóteses lineares da forma: H0 : Cβ = ξ versus H1 : Cβ 6= ξ, (2.18) onde C é uma matriz q × p, com q ≤ p, de caracteŕıstica completa q e ξ é um vector de dimensão q previamente especificado. Casos especiais importantes de (2.18) são: 2Esta tabela é reproduzida do trabalho de Wederburn (1976). 46 2. Inferência Tabela 2.1: Propriedades das estimativa de máxima verosimi- lhança de β para várias distribuições e funções de ligação; F significa existência de estimativa finita; I significa existência de estimativa no interior do espaço paramétrico; U significa unicidade. a) Modelos para os quais µ ≥ 0 ou µ > 0 se a função de ligação não estiver definida para µ = 0 função de ligação normal Poisson gama µλ(λ < −1) I F,I F,I µλ(−1 ≤ λ < 0) I F,I F,I,U lnλ I F,I,U F,I,U µλ(0 < λ < 1) F F,I,U F,I µ F,U F,I,U F,I µλ(λ > 1) F F,I F,I b) Modelos para a distribuição binomial função de ligação µ I,U sin−1 √ µ I,U Φ−1(µ) F,I,U ln µ 1−µ F,I,U ln{− ln(1 − µ)} F,I,U c) Modelos para a distribuição normal sem restrições a µ função de ligação µ F,I,U eµ F,I 2.3. Testes de Hipóteses 47 • Hipótese da nulidade de uma componente do vector parâmetro, nomeadamente H0 : βj = 0 versus H1 : βj 6= 0, para algum j, sendo neste caso q = 1, C = (0, ..., 0, 1, 0, ..., 0) e ocupando o 1 a j-ésima posição e ξ = 0. • Hipótese da nulidade de um subvector do vector parâmetro, H0 : βr = 0 versus H1 : βr 6= 0, para algum subvector de r componentes de β. Se tivermos, por exemplo H0 : (β1, ..., βr) T = (0, ..., 0)T , então q = r e C = ( Ir Or×(p−r) ) ξ = 0r, onde Ir é a matriz identidade de dimensão r, Or×(p−r) é uma matriz de zeros de dimensão r × (p − r) e 0r é o vector nulo de dimensão r. Estas hipóteses correspondem a testar submodelos do modelo original, importante na selecção de covariáveis, como se irá ver no caṕıtulo 3; a 1a corresponde a testar um submodelo com todas as covariáveis do modelo original à excepção da covariável zj relati- va ao parâmetro de regressão βj e a segunda corresponde a testar um modelo sem as r covariáveis relativas aos parâmetros supostos nulos sob a hipótese H0. Uma das situações em que isso acontece surge, por exemplo, quando uma variável é policotómica tomando, digamos, r + 1 valores distintos. Nesse caso é, como já se disse, aconselhável construir r variáveis dicotómicas para a representar, havendo nesse caso r parâmetros β’s associados a ela. Assim, para 50 2. Inferência A estat́ıstica de Wald é, em geral, a mais utilizada para testar hipóteses nulas sobre componentes individuais, embora também se use para testar hipóteses nulas do tipo βr = 0 quando o subvector βr representa o vector correspondente a uma recodificação de uma variável policotómica4. 2.3.2 Teste de razão de verosimilhanças A estat́ıstica de Wilks ou estat́ıstica de razão de verosimilhanças é definida por Λ = −2 ln maxH0 L(β) maxH0∪H1 L(β) = −2{ℓ(β̃) − ℓ(β̂)} (2.20) onde β̃, o estimador de máxima verosimilhança restrito, é o valor de β que maximiza a verosimilhança sujeito às restrições impostas pela hipótese Cβ = ξ. O teorema de Wilks (e.g., Cox and Hinkley, 1974) estabelece que, sob certas condições de regularidade , a estat́ıstica Λ tem, sob H0, uma distribuição assintótica de um χ2 sendo o número de graus de liberdade igual à diferença entre o número de parâmetros a estimar sob H0 ∪H1 (neste caso p) e o número de parâmetros a estimar sob H0 (neste caso p− q). Assim, sob H0, Λ = −2{ℓ(β̃) − ℓ(β̂)} a∼ χ2q. De acordo com o teste de razão de verosimilhanças a hipótese nula H0 : Cβ = ξ é rejeitada a favor de H1 : Cβ 6= ξ, a um ńıvel 4Por exemplo, o software SPSS usa a estat́ıstica de Wald para estas situações. 2.3. Testes de Hipóteses 51 de significância α, se o valor observado da estat́ıstica Λ for superior ao quantil de probabilidade 1 − α de um χ2q. Exemplo 2.4 Suponhamos que queremos testar H0 : βr = 0 versus H1 : βr 6= 0, onde βr é um subvector de β de dimensão r. O uso da estat́ıstica de razão de verosimilhanças não envolve grandes dificuldades computacionais, já que, para calcular a es- tat́ıstica Λ basta usar o método iterativo de mı́nimos quadrados ponderados para obter, (i) a estimativa de máxima verosimilhança, β̂0, do vector parâme- tro β0 que corresponde ao subvector de β sem as componentes que constituem βr, e a respectiva log-verosimilhança ℓ(β̂0), (ii) a estimativa de máxima verosimilhança β̂, do vector parâme- tro β e a respectiva log-verosimilhança ℓ(β̂), ou seja ajustar dois modelos aos dados (em que um é um submodelo do outro). Hipóteses mais gerais da forma (2.18) requerem, contudo, mais trabalho computacional. A estat́ıstica de razão de verosimilhanças é mais utilizada para comparar modelos que estão encaixados, isto é, modelos em que um é submodelo do outro, tal como iremos ver no caṕıtulo seguinte. 2.3.3 Estat́ıstica de Rao Seja novamente β̃ o estimador de máxima verosimilhança de β sujeito à restrição imposta pela hipótese nula Cβ = ξ. 52 2. Inferência A Estat́ıstica de Rao ou Estat́ıstica score para testar (2.18) é definida por U = [S(β̃)]TI−1(β̃)S(β̃) (2.21) A ideia por trás da sugestão desta estat́ıstica é a seguinte: Se β̂ é o estimador de β não restrito , isto é, calculado sem quaisquer restrições, então sabemos que s(β̂) = 0. Se substituirmos β̂ pelo estimador de máxima verosimilhança sob H0, isto é por β̃ s(β̃) será significativamente diferente de zero se H0 não for verdadeira. A estat́ıstica U mede assim a “distância” entre s(β̃) e 0. Tal como para os outros testes, usando a estat́ıstica score rejeita- se H0 a um ńıvel de significância α se o valor observado de U for superior ao quantil de probabilidade 1 − α de um χ2q . A estat́ıstica score é útil em situações em que já se calculou um estimador restrito para β. Tem a vantagem em relação à estat́ıstica de razão de verosimilhanças de não requerer o cálculo do estimador não restrito. Além disso, tal como a estat́ıstica de Wald pode ser utilizada em modelos com parâmetro de sobredispersão, já que para o seu cálculo só há necessidade de conhecer os momentos de 1a e 2a ordem. Como se viu, sob H0 as distribuições assintóticas das três es- tat́ısticas são idênticas. A qualidade da aproximação das distri- buições exactas das estat́ısticas Λ,W e U para a distribuição as- sintótica, depende da dimensão da amostra n e da forma do lo- garitmo da função de verosimilhança. Fahrmeir and Tutz (1994) exemplificam a situação com um modelo de Poisson com função de ligação linear. Se a log-verosimilhança for uma função quadrática em β então as três estat́ısticas coincidem. Para amostras pequenas pode ha- ver uma diferença considerável no valor destas estat́ısticas. Uma 2.4. Quasi-verosimilhança 55 Definição 2 Seja Ξ o espaço paramétrico, isto é µ ∈ Ξ. Diz-se que a função Q : Ξ → IR, definida por Q(µ, y) = ∫ µ y u(t, y)dt = ∫ µ y y − t φV (t) dt é uma função de quasi-verosimilhança (correctamente seria de quasi-log-verosimilhança). ♦ No caso em que temos n observações de variáveis aleatórias in- dependentes, definimos Q(µ, y) = n∑ i= Q(µi, yi) Esta função além de partilhar de muitas das propriedades formais que o logaritmo de uma função verosimilhança pode mesmo ser uma função de log-verosimilhança. Prova-se que se existir uma função de log-verosimilhança ℓ tal que ∂ℓ ∂µ = y − µ φV (µ) , com E(Y ) = µ e var(Y ) = φV (µ), então ℓ tem a estrutura corres- pondente a uma função de log-verosimilhança da famı́lia exponen- cial. Nos MLG sabemos que µi = h(ηi) = h(z T i β). Assim a função de quasi-verosimilhança Q(µ, y) é dada por Q(µ, y) = n∑ i=1 yi − h(zTi β) var(Yi) , e, se igualarmos a zero as derivadas de Q(µi, yi) em ordem a βj , para j = 1, ...p, obtemos o sistema de equações n∑ i=1 (yi − µi) V (µi) ∂µi ∂βj = 0 j = 1, ..., p 56 2. Inferência = n∑ i=1 (yi − µi)zij V (µi) ∂µi ∂ηi = 0 (2.23) as quais não dependem de φ e coincidem com as obtidas em (2.7). À função s∗(β) = ∂Q ∂β damos o nome de função quasi-score ou função de estimação generalizada e às equações em (2.23) damos o nome de equações de quasi-verosimilhança, sendo as estimativas resultantes, estimativas de quasi-máxima verosimilhança. Note-se que quando a função de variância é igual a 1 o método reduz-se ao método dos mı́nimos quadrados. As propriedades assintóticas dos estimadores de quasi-verosimi- lhança β̂ ∗ podem ser obtidas sob condições de regularidade seme- lhantes às necessárias para os estimadores de máxima verosimi- lhança. Em particular pode mostrar-se que β̂ ∗ a∼ Np ( β, (I∗)−1(β̂∗)V (β̂∗)(I∗)−1(β̂∗) ) , onde I∗(β) = E ( − ∂s ∗(β) ∂βT ) e V (β) = cov(s∗(β)). Comparando este resultado com o obtido quando os modelos estão completamente especificados, vemos que, essencialmente, ape- nas a matriz de covariância de β̂ ∗ tem de ser corrigida. Assim, o método de quasi-verosimilhança permite a obtenção de estimadores consistentes e assintoticamente normais para β, com apenas uma perda de eficiência. Para que esta perda de eficiência seja peque- na é necessário que a estrutura de variância proposta seja o mais próxima posśıvel da verdadeira estrutura de variância. Também é posśıvel proceder a testes de hipóteses semelhantes aos da secção anterior através do uso de estat́ısticas de Wald e de 2.4. Quasi-verosimilhança 57 Rao modificadas. Por exemplo, a Estat́ıstica de Wald modificada para testar H0 : Cβ = ξ versus H1 : Cβ 6= ξ, é dada por Wm = (Cβ̂ ∗ − ξ)T [CACT ]−1(Cβ̂∗ − ξ), onde A = (I∗)−1(β̂∗)V (β̂∗)(I∗)−1(β̂∗) é a matriz de covariância corrigida. Wm ainda tem uma distribuição assintótica χ2r, onde r é a caracteŕıstica da matriz C. Veja-se detalhes em White (1982). Tratamento semelhante pode ser feito para o caso em que o mo- delo se encontra completamente especificado, mas os parâmetros φ são distintos. Veja-se, para o efeito, Shao (1998, pg. 246-247). 3.1. Qualidade de Ajustamento 61 que ainda se ajusta adequadamente aos dados. Este modelo embora adaptando-se aos dados e podendo até ser adequado para réplicas do estudo, pode esconder caracteŕısticas ainda importantes dos dados. • Modelo corrente Em geral trabalha-se com modelos encaixados, isto é, passa- se do modelo maximal para o modelo minimal por exclusão de termos da desvio. O modelo corrente, é qualquer modelo com q parâmetros linearmente independentes situado entre o modelo maximal e o modelo minimal, e que está a ser sujeito a investigação. 3.1 Qualidade de Ajustamento 3.1.1 Função desvio O modelo saturado é útil para julgar da qualidade de ajusta- mento de um determinado modelo em investigação, que passamos a designar por M , através da introdução de uma medida da distância dos valores ajustados µ̂ com esse modelo e dos correspondentes va- lores observados y. Essa medida de discrepância entre o modelo saturado e o modelo corrente, é baseada na estat́ıstica de razão de verosimilhanças de Wilks referida na secção 2.3.25. Como vimos na secção 2.1.1, o logaritmo da função de verosimi- lhança (função log-verosimilhança) de um modelo linear generaliza- 5Seguindo a sugestão de Cordeiro (1986) traduzimos o termo “deviance” por desvio. 62 3. Selecção e Validação de Modelos do é dada por lnL(β) = ℓ(β) = n∑ i=1 ωi[yiq(µi) − b(q(µi))] φ + c(yi, φ, ωi) em que se substituiu θi por q(µi), para fazer salientar, na função log-verosimilhança , a relação funcional existente entre θi e µi. Como para o modelo saturado - que passamos a designar por S - se tem µ̂i = yi, o máximo da função log-verosimilhança para este modelo é ℓS(β̂S) = n∑ i=1 ωi[yiq(yi) − b(q(yi))] φ + c(yi, φ, ωi). Por outro lado, se designarmos por µ̂i a estimativa de máxima verosimilhança de µi, para i = 1, .., n, o máximo da função log- verosimilhança para o modelo em investigação com, digamos, m parâmetros na desvio é ℓM(β̂M) = n∑ i=1 ωi[yiq(µ̂i) − b(q(µ̂i))] φ + c(yi, φ, ωi). Os ı́ndices em β̂ e ℓ correspondem ao modelo em relação ao qual são calculados. Se compararmos o modelo em investigaçãoM com o modelo satu- rado S através da estat́ıstica de razão de verosimilhanças, obtemos D∗(y; µ̂) = −2(ℓM(β̂M) − ℓS(β̂S)) = −2 ∑ i ωi φ {[ yiq(µ̂i) − b(q(µ̂i)) ] − [ yiq(yi) − b(q(yi)) ]} = D(y; µ̂) φ . (3.1) A D∗(y; µ̂) definido em (3.1) damos o nome de desvio reduzido; ao numerador D(y; µ̂) damos o nome de desvio para o modelo corrente. Note-se que o desvio é só função dos dados. 3.1. Qualidade de Ajustamento 63 Como se pode observar de (3.1) o desvio pode ser decomposto D(y; µ̂) = ∑ i 2ωi { yi(q(yi) − q(µ̂i)) − b(q(yi)) + b(q(µ̂i)) } = ∑ i di na soma de parcelas di que medem a diferença dos logaritmos das verosimilhanças observada e ajustada para cada observação. A so- ma destas componentes é assim uma medida da discrepância total entre as duas log-verosimilhanças. É fácil de verificar que o desvio é sempre maior ou igual a zero, e decresce à medida que covariáveis vão sendo adicionadas ao modelo nulo, tomando obviamente o valor zero para o modelo saturado. Uma outra propriedade importante do desvio é a aditividade pa- ra modelos encaixados. Com efeito, suponhamos que temos dois modelos intermédios M1 e M2 estando M2 encaixado em M1, is- to é, são modelos do mesmo tipo, mas o modelo M2 contém me- nos parâmetros na desvio que o modelo M1. Se designarmos por D(y; µ̂j) o desvio do modelo Mj , j = 1, 2, então a estat́ıstica da razão de verosimilhanças para comparar estes dois modelos resume- se a −2(ℓM2(β̂2) − ℓM1(β̂1)) = D(y; µ̂2) −D(y; µ̂1) φ . Dos resultados do caṕıtulo anterior sabe-se que, sob a hipótese do modelo M1 ser verdadeiro, então D(y; µ̂2) −D(y; µ̂1) φ a∼ χ2p1−p2, onde pj , representa a dimensão do vector β para o modelo Mj , j = 1, 2. A comparação de modelos encaixados, pode então ser feitas à custa da diferença dos desvios de cada modelo. 66 3. Selecção e Validação de Modelos Outros exemplos da função desvio para modelos lineares genera- lizados encontram-se na tabela 3.1. Tabela 3.1: Expressões da função desvio para alguns modelos. normal ∑ i(yi − µ̂i)2 Poisson 2[ ∑ i yi ln yi µ̂i −∑i(yi − µ̂i)] binomial 2[ ∑ imiyi ln yi µ̂i + ∑ imi(1 − yi) ln 1−yi)1−µ̂i ] gama 2 ∑ i { − ln yi µ̂i + yi−µ̂i µ̂i } gaussiana inversa ∑ i (yi−µ̂i) 2 yiµ̂2i 3.1.2 Estat́ıstica de Pearson generalizada Outra medida da adequabilidade de modelos é a estat́ıstica de Pearson generalizada já definida na secção 2.2.2, X2 = ∑ i ωi(yi − µ̂i)2 V (µ̂i) . (3.2) Para a distribuição normal, a estat́ıstica X2 coincide, tal como o desvio, com a soma dos quadrados dos reśıduos. Para os modelos Poisson e Binomial coincide com a estat́ıstica original de Pearson. Novamente é costume usar X2 para testar a adequabilidade de um modelo comparando o valor observado com o quantil de probabili- dade 1−α de uma distribuição de χ2 com n−p graus de liberdade. Contudo, tal como acontece com o desvio, a aproximação pelo χ2 da distribuição de X2 pode ser, em certos modelos, má mesmo pa- ra grandes amostras, havendo a necessidade de agrupar os dados 3.2. Selecção de Modelos 67 o mais posśıvel, garantindo ao mesmo tempo que o número de ob- servações em cada grupo, digamos ni não seja pequeno. Assim, as estat́ısticas a usar para testar a adequabilidade do modelo em consideração devem ser do tipo D(y; µ̂) = g∑ i=1 2ωi { yi(q(yi) − q(µ̂i)) − b(q(yi)) + b(q(µ̂i)) } , e X2 = g∑ i=1 ωi(yi − µ̂i)2 V (µ̂i) , onde g é o número de grupos, e o número de observações em ca- da grupo é suficientemente grande em todos os grupos. Neste ca- so, a suposição de que ambas as estat́ısticas têm uma distribuição aproximada de um φχ2 com g − p graus de liberdade, já é menos problemática.7 A propriedade da aditividade da função desvio faz com que esta seja preferida, em relação à estat́ıstica de Pearson, como uma me- dida da discrepância, embora esta última tenha a vantagem de ter uma interpretação mais directa. 3.2 Selecção de Modelos Como já se disse, em problemas práticos que requerem uma análise estat́ıstica via modelos lineares generalizados, há geralmente um número elevado de covariáveis que podem ser potencialmen- te importantes para explicar a variabilidade inerente aos dados. Também tem, frequentemente, interesse investigar a influência de posśıveis interacções entre as covariáveis. Isto implica obviamen- te a existência de um número elevado de modelos a considerar de 7Consulte-se, e.g., Fahrmeir and Tutz (1994, pg. 48). 68 3. Selecção e Validação de Modelos modo a escolher um modelo posśıvel para explicar o fenómeno em estudo. Se pensarmos num modelo maximal como aquele que entra em linha de conta, na sua desvio, com todas as posśıveis covariáveis e interacções de interesse entre elas, um submodelo deste é qual- quer modelo que é obtido dele por exclusão de algum ou alguns dos termos da desvio. A existência de um número elevado de modelos a considerar, quer se parta do modelo maximal, quer se parta do modelo minimal, traz problemas de ordem combinatória - o número de combinações posśıveis torna-se rapidamente não manejável - e de ordem estat́ıstica - como decidir sobre o equiĺıbrio entre o efeito da inclusão ou exclusão de um termo na discrepância entre y e µ̂ e a complexidade de um modelo maior? Há pois necessidade de esta- belecer uma estratégia para a selecção do melhor, ou dos melhores modelos, já que raramente se pode falar na existência de um único “melhor modelo”. Ter um submodelo M1 de um modelo M , com vector parâmetro β de dimensão p, corresponde a ter um modelo com vector de parâmetros β1 que é um subvector de β. Sem perda de genera- lidade, podemos assumir a partição de β em (β1,β2) T . Assim, a adequabilidade de um submodelo pode ser testada formalmente co- mo H0 : β2 = 0, versus H1 : β2 6= 0. Esta hipótese pode ser testada usando a metodologia explicada na secção 2.3. Designemos, como anteriormente, a função score, a matriz de in- formação de Fisher e a sua inversa para o modelo especificado em H1, respectivamente, por s(β), I(β) e A(β). Então, em conformi- dade com a partição relativa ao vector β, podemos particioná-las 3.2. Selecção de Modelos 71 modelo em H0 é, AIC = −2ℓ(β̃1, 0, φ̃) + 2r, onde r = dim(β1). Um valor baixo para AIC é considerado como representativo de um melhor ajustamento e na selecção de modelos devemos ter como objectivo a minimização de AIC. Note-se a se- guinte relação existente entre AIC e o desvio reduzido relativo ao modelo especificado por H0 (estamos a supor que o parâmetro φ ou é conhecido, ou é substitúıdo por uma estimativa consistente) AICr = −2ℓ(β̃1, 0) + 2ℓ(β̂S) − 2ℓ(β̂S) + 2r = D∗r + 2r − 2ℓ(β̂S) onde o indice r serve para especificar o modelo em consideração e S, como habitualmente, refere-se ao modelo saturado. Cordeiro (1986) sugere ainda a seguinte modificação do critério de Akaike para seleccionar modelos, C∗r = D ∗ r + 2r − n = AICr + 2ℓ(β̂S) − n. Um gráfico de C∗r contra r fornece uma boa indicação para compa- ração de modelos. Se o modelo for verdadeiro é de esperar que C∗r seja próximo de r. Se tivermos dois modelos encaixados M1 e M2 com, digamos r1 e r2 parâmetros respectivamente, onde r1 > r2, vem AICr1 −AICr2 = C∗r1 − C∗r2 = D∗r1 −D∗r2 + 2(r1 − r2) e, supondo que o modelo M2 é verdadeiro, tem-se (Cordeiro, 1986) E(AICr1 −AICr2) = r1 − r2 + 0(n−1). Na comparação de modelos sucessivamente mais ricos, o declive esperado do segmento de recta 72 3. Selecção e Validação de Modelos que une AICr1 e AICr2 deve estar próximo de 1 supondo o modelo menor M2 verdadeiro. Pares de modelos que exibem declive maior do que 1 são indicação de que o modelo maior não é significativa- mente melhor que o modelo menor. 3.3 Análise de Reśıduos A análise de reśıduos é útil, não só para uma avaliação local da qualidade de ajustamento de um modelo no que diz respeito à escolha da distribuição, da função de ligação e de termos do predi- tor linear, como também para ajudar a identificar observações mal ajustadas, i.e., que não são bem explicadas pelo modelo. Um reśıduo Ri deve exprimir a discrepância entre o valor obser- vado yi e o valor µ̂i ajustado pelo modelo. No modelo linear normal em que o vector das respostas Y se pode escrever na forma Y = µ + ε = Zβ + ε, ε ∼ Nn(0, σ2I) tem-se β̂ = (ZTZ)−1ZTy e o vector dos reśıduos é naturalmente dado por R = y − ŷ onde ŷ = µ̂ = Zβ̂ é o vector dos valores ajustados. No caso dos modelos lineares generalizados não existe necessariamente uma componente εi para o qual o reśıduo Ri seja uma estimativa; faz portanto sentido, como iremos ver, considerar outras definições de reśıduos. Outra quantidade de interesse na análise dos reśıduos, no caso do modelo linear normal, é a matriz de projecção H = Z(ZTZ)−1ZT a qual é tal que ŷ = Hy (e que por essa razão se designa, em inglês, por hat matrix). Esta matriz é simétrica e idempotente. Os seus elementos hij são uma medida da influência exercida por yj em ŷi. A influência exercida por yi em ŷi é reflectida pelo elemento da 3.3. Análise de Reśıduos 73 diagonal principal hii. Como se tem que ∑ hii = p e 0 ≤ hii ≤ 1, um ponto é considerado influente se hii > 2p n (Hoaglin and Welsch, 1978). 3.3.1 Matriz de projecção generalizada Como vimos em (2.14), o processo iterativo conduz, no modelo linear generalizado a β̂ = ( ZTWZ )−1 ZTWu. Esta equação é idêntica à que se obteria para os estimadores de mı́nimos quadrados ponderados para o problema de regressão u = Zβ + ε, ou, alternativamente, a solução de minimos quadrados para o mo- delo linear u0 = Z0β + ε̃, onde u0 = W 1 2 u, Z0 = W 1 2Z e, portanto, a matriz de projecção correspondente é H = Z0(Z T 0 Z0) −1ZT0 = W 1 2Z(ZTW 1 2W 1 2Z)−1ZTW 1 2 = W 1 2ZI−1(β)ZTW 1 2 . (3.3) O lado direito de (3.3) advém do facto de I(β) = ZTWZ, tal como se viu em (2.10). Esta matriz H é também simétrica e idempotente e pode ser vista como uma matriz de projecção para a qual ainda se tem traço(H)=caracteŕıstica(H). Os elementos da diagonal princi- pal desta matriz são ainda tais que 0 ≤ hii ≤ 1 e valores elevados de hii correspondem a pontos extremos. Contudo, em contraste com o modelo linear normal, esta matriz não depende apenas da matriz Z, 76 3. Selecção e Validação de Modelos onde V (x) é a função de variância. Outro tipo de reśıduo é baseado na função desvio. Podemos usar a contribuição da i-ésima observação para a função desvio definida em (3.1), D(y; µ̂) = ∑ i 2ωi { yi(q(yi) − q(µ̂i)) − b(q(yi)) + b(q(µ̂i)) } = ∑ i di, para dar uma nova definição de reśıduo. Assim o desvio residual correspondente à i-ésima observação é definido por RDi = δi √ di, (3.7) onde δi = sinal(yi − µ̂i). O desvio residual padronizado é também obtido dividindo o desvio residual RDi por √ φ̂(1 − hii), ou seja R⋆Di = RDi√ φ̂(1 − hii) . (3.8) Exemplo 3.3 Reśıduos no modelo normal Para o modelo normal é fácil de verificar que os três tipos de reśıduos coincidem. Com efeito, atendendo a que para este modelo V (x) = 1, tem-se que A(x) = x e, portanto o reśıduo de Pearson (puro, i.e., não padronizado), RPi = yi − µ̂i, coincide com o reśıduo de Anscombe. Por outro lado, dado que di = (yi − µ̂i)2, tem-se que o desvio residual é também dado por yi − µ̂i Exemplo 3.4 Reśıduos no modelo Poisson No modelo Poisson tem-se, como se sabe, V (x) = x. Deste mo- do ∫ V −1/3(x)dx = 3 2 x2/3 e portanto os reśıduos de Pearson e de Anscombe puros são, respectivamente RPi = yi − µ̂i µ̂ 1/2 i RAi = 3(y 2/3 i − µ̂2/3i ) 2µ̂ 1/6 i . 3.3. Análise de Reśıduos 77 Consultando a tabela 3.1 facilmente se obtém para o desvio residual RDi = δi2 1/2(yi ln yi µ̂i − yi + µ̂i)1/2, onde δi = sinal(yi − µ̂i). Exemplo 3.5 Reśıduos no modelo binomial No modelo binomial tem-se que V (x) = x(1 − x) e ωi = mi. Deste modo A(x) = ∫ x−1/3(1 − x)−1/3dx e o reśıduo de Anscombe é dado por: RAi = m 1/2 i [A(yi) − A(µ̂i)] [µ̂i(1 − µ̂i)]1/6 . Cox and Snell (1968) calculam este reśıduo através da função beta incompleta. Facilmente se obtêm as seguintes expressões para os reśıduos de Pearson e desvio residual RPi = m 1/2 i (yi−µ̂i) [µ̂i(1−µ̂i)]1/2 RDi = δi [ 2mi ( ln yi µ̂i + (1 − yi) ln ( 1−yi 1−µ̂i ))]1/2 , onde δi = sinal(yi − µ̂i). Cordeiro (1986) faz um estudo comparativo entre os reśıduos de Anscombe e os desvios residuais para os modelos Poisson, gama e gaussiano inverso. Pierce and Schafer (1986) introduzem uma generalização dos reśıduos de Anscombe e sugerem outro tipo de reśıduos destinados a estabilizar a variância. Na tabela 3.2 apresentamos um quadro resumo dos três tipos de reśıduos para os modelos que temos vindo a considerar. 78 3. S elecção e V alid ação d e M o d elos Tabela 3.2: Expressões dos reśıduos para alguns modelos. RPi R A i R D i normal yi − µ̂i yi − µ̂i yi − µ̂i Poisson yi−µ̂i µ̂ 1/2 i 3(y 2/3 i −µ̂ 2/3 i ) 2µ̂ 1/6 i δi2 1/2(yi ln yi µ̂i − yi + µ̂i)1/2 binomial m 1/2 i (yi−µ̂i) [µ̂i(1−µ̂i)]1/2 m 1/2 i [A(yi)−A(µ̂i)] [µ̂i(1−µ̂i)]1/6 δi[2mi(ln yi µ̂i + (1 − yi) ln 1−yi1−µ̂i )] 1/2 gama yi−µ̂i µ̂i 3(y 1/3 i −µ̂ 1/3 i ) µ̂ 1/3 i δi[2(ln µ̂i yi + yi−µ̂i µ̂i )]1/2 gaussiana inversa yi−µ̂i µ̂ 3/2 i µ̂ −1/2 i ln yi µ̂i yi−µ̂i y 1/2 i µ̂i δi = sinal(yi − µ̂i) e A(x) = ∫ [x(1 − x)]−1/3dx. 3.3. Análise de Reśıduos 81 lado uma tendência negativa indica a situação contrária, is- to é, a variância cresce demasiadamente rápido em relação à média e, portanto, uma função de variância do tipo, e.g., V (µ) ∝ µ, deve ser substitúıda por uma função de variância do tipo V (µ) ∝ µλ, λ < 1. É posśıvel fazer uma avaliação formal da função de variância e estimar, e.g., a quantidade λ adequada. Para tal veja-se, e.g., a secção 4.3 de Fahrmeir and Tutz (1994), ou a secção 12.6.2 de McCullagh and Nelder (1989). • Avaliação da função de ligação Para uma avaliação informal da adequabilidade da função de ligação, pode fazer-se uma representação gráfica da variável dependente ajustada definida em (2.13) u = η̂ + D̂(y − µ̂), contra η̂, onde D̂ significa que a matriz diagonal D com ele- mento genérico ∂ηi ∂µi é calculada para os valores estimados. Se os pontos se distribúırem, aproximadamente, sobre uma linha recta, então a função de ligação é adequada; se, por outro la- do, se observar uma curvatura para cima, isso significa que é necessário usar uma função de ligação com potência superior; uma curvatura para baixo é sinónimo da situação contrária. Existem métodos formais para a avaliação da adequabilidade da função de ligação. Hinkley (1985) sugere considerar η̂2 como uma nova covariável a adicionar ao preditor linear e verificar se há um decĺınio da função desvio. Veja-se ainda a secção 12.6.3 de McCullagh and Nelder (1989). • Averiguação da adequabilidade da escala em que as covariáveis estão representadas 82 3. Selecção e Validação de Modelos Uma má escolha da escala em que uma ou mais covariáveis estão representadas pode afectar a adequabilidade do modelo de modo a parecer, erradamente, que outro tipo de anomalias estão presentes, tal como, por exemplo, uma má escolha da função de ligação; é pois importante saber distinguir em que situação nos encontramos. O objectivo, neste caso, é o de averiguar se um termo no preditor linear, do tipo, e.g., βx, deve ser substitúıdo por um termo do tipo βh(x, θ), onde h(·, θ) é uma transformação apropriada da covariável em questão (do tipo, por exemplo, Box-Cox); Uma representação gráfica adequada é a dos reśıduos parciais Rparciais,i contra os valores observados da covariável xi, sendo o vector de reśıduos parciais definidos por Rparcial = u− η̂ + γ̂x, onde u− η̂ é o vector dos reśıduos medidos na escala linear, u é a variável dependente ajustada já definida e γ̂ é a estimativa do parâmetro para a variável explicativa em consideração. Se a escala de x é satisfatória o gráfico deve ser aproximada- mente linear. • Avaliação da omissão de uma covariável Para averiguar se uma covariável que se omitiu, digamos z∗, deve ou não ser inclúıda no modelo, deve fazer-se uma repre- sentação gráfica dos reśıduos aumentados contra z∗. Esses são definidos do modo que passamos a descrever.8 8Baseado na secção 7.6 de Cordeiro (1986). 3.3. Análise de Reśıduos 83 Suponhamos que a componente sistemática correcta deve con- ter uma covariável adicional, isto é, g(µ) = Zβ + h(z∗; γ), onde h(·; γ) pode representar (i) um termo adicional em uma ou mais covariáveis originais, e.g., um termo quadrático ou uma interacção; (ii) uma contribuição linear ou não linear de alguma covariável omitida, e.g., h(z∗; γ) = z∗γ ou h(z∗; γ) = γ z∗ . O objectivo é definir reśıduos R̃ para o modelo ajustado η = Zβ tal que E(R̃) = h(z∗; γ). Se isto acontecer, um gráfico de R̃i contra z ∗ i exibirá, despre- zando a variação aleatória, a função h(z∗; γ). De um modo semelhante ao que se fez para definir os reśıduos parciais, os reśıduos aumentados são obtidos através do acrés- cimo [Z(ZT ŴZ)−1ZT Ŵ ]h(z∗; γ̂) aos reśıduos medidos na es- cala linear R = u − η̂ = [I − Z(ZTŴZ)−1ZT Ŵ ]u. Assim o vector dos reśıduos aumentados é dado por R̃ = [I − Z(ZTŴZ)−1ZT Ŵ ]u + Z(ZT ŴZ)−1ZT Ŵ ]h(z∗; γ̂). A análise gráfica dos reśıduos aumentados pode ser bastante útil na selecção de covariáveis, quando se tem muitas cova- riáveis a considerar. A formação da componente sistemática pode ser feita, passo a passo, com a introdução de uma única covariável de cada vez pelo método descrito. 86 3. Selecção e Validação de Modelos através de Z, como também do modelo adaptado através de W . Assim, uma observação extrema, isto é, com um valor elevado para uma ou mais das covariáveis, não tem necessariamente uma reper- cussão elevada se o seu peso (elemento correspondente de W ) for pequeno. Gráficos de hii contra µ̂i, ou contra i, são geralmente úteis na identificação de pontos com repercussão elevada. 3.4.2 Medida de influência Um indicador da influência da i-ésima observação (yi, zi) no vec- tor estimado β̂, pode ser calculado pela diferença β̂ − β̂(i), onde β̂(i) e β̂ representam, respectivamente, as estimativas de máxima verosimilhança do vector parâmetro β obtidas da amostra sem a observação (yi, zi) e da amostra com todas as observações. Se β̂(i) for substancialmente diferente β̂, então a observação (yi, zi) pode ser considerada influente. No caso dos modelos lineares normais a medida de influência utilizada é a sugerida por Cook (1977) Di = (β̂− β̂(i))T (ZTZ)(β̂ − β̂(i))/ps 2. Dado que, agora cov(β̂) = (ZTWZ)−1, tal como se viu em 2.2.3, é natural considerar como generalização da medida de influência de Cook Di = (β̂ − β̂(i))T (ZTWZ)(β̂ − β̂(i)) pφ̂ . (3.9) Contudo, no caso dos MLG, a estimação de β̂(i) necessita do recurso a métodos iterativos. O processo é, pois, computacionalmente caro para poder ser feito para todas as observações; alternativamente pode obter-se uma aproximação a β̂(i) fazendo apenas o 1 o passo do 3.4. Observações Discordantes 87 processo iterativo, usando β̂ como valor inicial; recorrendo a (2.14), a estimativa a um passo é dada por β̂(i),1 = I−1(i) (β̂)ZT(i)W(i)(β̂)u(β̂), (3.10) onde o ı́ndice (i) refere que os cálculos são feitos sem a i-ésima ob- servação. Para o cálculo aproximado de Di usa-se esta aproximação de β̂(i) em (3.9). Uma fórmula mais simples para β̂(i),1 é dada por Williams (1987), nomeadamente β̂(i),1 = β̂ −̟ 1/2 i (1 − hii)−1/2R⋆Pi (ZTWZ)−1zi, (3.11) onde ̟i, elemento de W está definido em (2.11). Aplicações do modelo linear generalizado mostram que a utili- zação de β̂(i),1 em (3.9) em vez do verdadeiro valor β̂(i) subestima o valor de Di. Contudo, segundo Williams (1987), o ponto importan- te é que a aproximação considerada geralmente identifica os casos mais influentes. Após serem identificados estes casos, pode então obter-se exactamente a influência da sua omissão. 3.4.3 Medida de consistência A possibilidade de uma determinada observação (yi, zi) ser in- consistente pode também ser averiguada adaptando o modelo sem essa observação e calculando os reśıduos da observação eliminada em relação ao correspondente valor predito µ̂(i) = h(z T i β̂(i)). Os reśıduos assim obtidos são chamados reśıduos de eliminação (“dele- tion residuals”). Por exemplo, a correspondente expressão para os reśıduos de eliminação de Pearson é R⋆P(i) = (yi − µ̂(i))wi [φ̂V (µ̂(i))(1 + h(ii))]1/2 , 88 3. Selecção e Validação de Modelos onde h(ii) = z T i (Z T (i)W(i)Z(i)) −1zi e o sinal + no denominador de R ⋆P (i) aparece devido ao facto de yi e µ̂(i) serem agora independentes. Novamente, o cálculo exacto destes reśıduos de eliminação pa- ra todas as observações é computacionalmente dispendioso e é ne- cessário encontrar outra alternativa. Em geral a solução encontrada é o cálculo desses reśıduos usando as estimativas obtidas após o 1o passo do processo iterativo, tal como se referiu na secção anterior. Williams (1987) sugere a utilização de outro tipo de reśıduos que ele denomina por reśıduo de verosimilhança. Para o efeito, seja Gi a redução operada no desvio reduzido quando se elimina do modelo a i-ésima observação, i.e. Gi = D ∗(y; µ̂) −D∗(y(i), µ̂(i)) = φ−1[di + ∑ j 6=i dj − ∑ j 6=i d(i),j], onde, y(i) designa o vector y sem o elemento yi e dj foi definido na secção 3.1.1. Assim, a contribuição de yi para Gi é exactamente φ −1di. Por ou- tro lado, Williams (1987) mostra que, usando a aproximação de β̂(i) dada em (3.11), e fazendo um desenvolvimento em série de Taylor da função desvio, o decréscimo de ∑ j 6=i dj quando µ̂ é substitúıdo por µ̂(i) é aproximadamente φ hii(R ⋆P i ) 2. Deste modo Gi pode ser aproximado por R2Gi definido por R2Gi = φ −1di + hii(R ⋆P i ) 2 = (1 − hii)(R⋆Di )2 + hii(R⋆Pi )2 (3.12) O reśıduo de verosimilhança é então dado por R∗Gi = δi √ (1 − hii)(R⋆Di )2 + hii(R⋆Pi )2, onde δi = sinal(yi − µ̂i). 92 4. Aplicações I: Modelos Discretos no ou benigno. Na secção 4.2 encontra-se um exemplo de dados na forma de proporção com a adopção de três MLG: loǵıstico, probit e complementar log-log. Na última secção analisa-se um conjunto de dados em forma de contagens através de modelos log-lineares que desempenham um papel importante na análise de dados categori- zados. 4.1 Modelos de Regressão Loǵıstica Na subsecção 1.4.2 apresentámos potenciais modelos lineares ge- neralizados para analisar dados binários, sendo o modelo loǵıstico o mais popular desses modelos, provavelmente, devido à simplicidade da sua implementação computacional. A adopção do modelo estru- tural (1.11) para a probabilidade de sucesso caracteriza o modelo de regressão loǵıstica. Exemplo 4.1 Processo Infeccioso Pulmonar No sector de Anatomia e Patologia do Hospital Heliópolis (São Paulo/Brasil) realizou-se um estudo retrospectivo com 175 pacien- tes entre 1970 e 1982, cujos dados se encontram em Paula et al. (1984). O objectivo principal desse estudo era avaliar a associação entre algumas variáveis histológicas e o tipo, maligno ou benigno, do Processo Infeccioso Pulmonar (PIP). Nesse estudo de caso-controle, os casos foram todos os pacien- tes diagnosticados, no peŕıodo e hospital há pouco mencionados, como portadores do PIP de origem maligna (71 pacientes). Os con- troles foram formados por uma amostra de 104 pacientes de uma população de 270, os quais foram também diagnosticados na mesma época e local e tiveram confirmado o PIP de origem benigna. 4.1. Modelos de Regressão Loǵıstica 93 A observação de cada um dos pacientes fez-se através de va- riáveis histológicas nos fragmentos de tecidos retirados da região pulmonar. Dessas variáveis somente as intensidades de histiócitos- linfócitos (HL) e de fibrose-frouxa (FF) foram consideradas impor- tantes na discriminação dos dois tipos de PIP. Porém, o conjunto de covariáveis será formado por dois factores potenciais de confun- dimento, sexo e idade. A descrição da codificação destas variáveis encontra-se na tabela 4.1. Tabela 4.1: Variáveis do processo infeccioso pulmonar. variáveis codificação tipo de PIP (Y ) 1=maligno 0=benigno idade (x1) em anos sexo (x2) 1=masculino 0=feminino intensidade de histiócitos-linfócitos 1=alta(3,4) HL (x3) 0=baixa(1,2) intensidade de fibrose frouxa 1=alta(3,4) FF (x4) 0=baixa(1,2) Fonte: Paula et al. (1984). Note-se que os dados sobre a variável resposta binária Y e o vec- tor de covariáveis x = (x1, . . . , x4) T não foram obtidos prospectiva- mente. Isto é, os dados sobre os casos e controles resultam de uma amostragem directa de um modelo para P (x | Y ), Y = 0, 1, con- trariamente aos dados prospectivos que estão associados ao modelo π(x) = P (Y | x). Entretanto, a sua análise pode ser processada de 94 4. Aplicações I: Modelos Discretos modo análogo àquele previsto para os dados de um estudo prospec- tivo, visto que o uso de um modelo prospectivo loǵıstico revela-se conveniente pelo facto da metodologia de máxima verosimilhança para os dados retrospectivos envolver ainda um modelo numa for- ma loǵıstica. Para maiores detalhes, veja Silva (1992, sec. 2.8). De acordo com a caracteŕıstica deste estudo retrospectivo, o modelo de regressão loǵıstico (1.11) será ajustado a estes dados, i.e., a probabilidade de tipo maligno do PIP para o i-ésimo pa- ciente (πi) está relacionada com o seu vector de covariáveis zi = (1, xi1, . . . , xi4) T através de πi = exp(zTi β) 1 + exp(zTi β) , (4.1) onde β = (β0, β1, . . . , β4) T e i = 1, . . . , 175. Os parâmetros de regressão βj , j 6= 0, são estimados directamente deste modelo, en- quanto β0 pode ser estimado posteriormente, se as probabilidades de selecção amostral dos tipos de PIP, φ1 e φ0, forem conhecidas. 4.1.1 Selecção do modelo loǵıstico Para a selecção de covariáveis que formem o “melhor” modelo loǵıstico usaremos um método de selecção stepwise baseados em p-values relativos aos testes de razão de verosimilhanças de Wilks entre modelos com inclusão ou exclusão de covariáveis, ou mesmo de suas interacções. Neste caso, o grau de importância de uma co- variável é medido pelo p-value do teste da razão de verosimilhanças entre os modelos que a incluem e a excluem. Quanto menor for este valor tanto mais importante será considerada a covariável. Como a covariável mais importante por este critério não é necessariamente significativa do ponto de vista estat́ıstico, há que impor um limite
Docsity logo



Copyright © 2024 Ladybird Srl - Via Leonardo da Vinci 16, 10126, Torino, Italy - VAT 10816460017 - All rights reserved