Regressão Logística Continuaremos com os modelos de regressão não linear para duas situações importantes onde as variáveis respostas são discretas e os erros não são normalmente distribuídos. O modelo de regressão não linear logístico é utilizado quando a variável resposta é qualitativa com dois resultados possíveis, por exemplo, sobrevivência de enxertos de ameixeiras (sobrevive ou não sobrevive). Este modelo pode ser estendido quando a variável resposta qualitativa tem mais do que duas categorias; por exemplo, a pressão sanguínea pode ser classificada como alta, normal e baixa. O modelo de regressão logístico, aqui estudado, pode ser utilizado para analisar dados observacionais ou experimentais no delineamento inteiramente casualizado. 1 Modelos de regressão com variáveis respostas binárias Em muitos estudos a variável resposta qualitativa tem duas possibilidades e, assim, pode ser representada pela variável indicadora, recebendo os valores 0 (zero) e 1 (um). Exemplos: 1) o objetivo da análise é verificar se uma firma tem ou não um departamento de relação industrial, de acordo com o seu tamanho. A variável resposta têm duas possibilidades: a firma tem ou não tem o departamento. Estes resultados pode m ser codificados como 1 e 0 (ou vice-versa). 2) Num estudo sobre a participação das esposas no mercado de trabalho, como função da idade da esposa, número de filhos e rendimento do marido, a variável resposta Y foi definida do seguinte modo: a mulher participa no mercado de trabalho ou não. Novamente, estas respostas podem ser codificadas como 1 e 0, respectivamente. 2 Interpretação da função de resposta quando a variável resposta é binária Vamos considerar o modelo de regressão linear simples: Yi 0 1 X i i 1 Yi 0 A resposta esperada é dada por: E(Yi ) 0 1 X i (1) Considere Yi uma v.a. Bernoulli com distribuição de probabilidade: Yi 1 P(Yi 1) i Yi 0 P(Yi 0) 1 i 3 Pela definição de valor esperado, obtemos: E(Yi ) i (2) Igualando-se (1) e (2), obtemos: E(Yi ) 0 1 X i i (3) Assim, a resposta média, quando a variável resposta é uma variável binária (1 ou 0), sempre representa a probabilidade de Y = 1, para o nível da variável preditora Xi. Na figura a seguir, a variável indicadora Y corresponde se a firma tem ou não tem um departamento de relação industrial e, a variável preditora X corresponde ao tamanho da firma. A função resposta mostra a probabilidade que uma firma de um dado tamanho tenha um departamento de relação industrial. 4 Probabilidade de uma firma ter departamento E(Y) 1 E(Y ) 0 1 X 0 Tamanho da firma X 5 Problemas quando a variável resposta é binária 1. Os erros não tem distribuição normal. Cada erro i Yi ( 0 1 X i ) pode assumir um dos dois valores: Yi 1 i 1 0 1 X i Yi 0 i 0 1 X i 2. Variâncias heterogêneas. A variância de Yi para o modelo de regressão linear simples é: 2 (Yi ) E[(Yi E (Yi ))2 ] (1 )2 i (0 i )2 (1 i ) 2 (Yi ) i (1 i ) E (Yi )(1 E (Yi )) 6 Como: i Yi i ( i constante) Temos: 2 (i ) i (1 i ) ( 0 1 X i ) (1 0 1 X i ) Depende de Xi 3. Restrição na função resposta. Como a função resposta representa probabilidades quando a variável resposta é binária, então: 0 E (Y ) 1 (4) A restrição na resposta média de apresentar valores no intervalo [0;1], freqüentemente é inapropriada, ou mesmo impossível, para uma função de resposta linear. Para o exemplo do departamento de relação industrial, o uso da função de resposta linear, sujeita a restrição na resposta média, requer probabilidade 0 (zero) na resposta média para todas as firmas pequenas e, uma probabilidade 1 (um) na resposta média, para todas as firmas grandes. 7 Este modelo freqüentemente não representa bem a situação em estudo. Ao invés, um modelo onde as probabilidades 0 e 1 são encontradas assintoticamente, como mostra a figura a seguir, é, de modo geral, mais apropriada. E(Y ) exp( 100.1 X ) 1exp( 100.1 X ) Figura: Função resposta logística 8 E(Y ) exp(100,1 X ) 1exp(100,1 X ) Figura: Função resposta logística 9 Função resposta logística com uma única variável preditora Considerações teóricas e práticas sugerem que quando a variável resposta é binária, a forma da função resposta será frequentemente curvilínea. Nas duas figuras anteriores, temos 2 funções respostas adequadas para uma variável resposta binária. Elas tem assíntotas em 0 e 1 e, assim, estão de acordo com a restrição (4). As funções respostas das figuras são denominadas funções logísticas, cuja expressão é: E (Y ) exp( 0 1 X ) 1 exp( 0 1 X ) (5) Forma equivalente: E(Y ) 1 exp 0 1 X 1 10 Propriedade da função logística Uma propriedade interessante é que a função logística pode ser linearizada. Denotando-se E(Y) por , pois a resposta média é a probabilidade quando a variável resposta é binária. Fazendo-se a transformação: 1 ' loge obtemos: ' 0 1 X (6) (7) Esta transformação é chamada de transformação logit da probabilidade . A razão /(1- ) na transformação logit é chamada de Odds (Chance). A função resposta transformada (7) é denominada como função resposta logit, e ’ é denominada de resposta média logit. Observe em (7) que: - ’ para -X. 11 Usos da função logística • Descritivo: descrever a natureza do relacionamento entre a resposta média (isto é, a probabilidade de comprar, por exemplo) e uma (ou mais) variáveis regressoras. • Preditivo: saber se uma pessoa irá comprar um automóvel no próximo ano, dado o seu rendimento. 12 Variável “threshold”. Exemplo: considere a força necessária para quebrar blocos de concreto, medida em libras por polegada ao quadrado. Assume-se que cada bloco tenha a sua variável threshold Ti, ou seja, ele irá quebrar-se se uma força igual ou maior do que Ti for aplicada e, não irá quebrar-se se uma força menor do que Ti for aplicada. Um bloco pode ser testado com apenas uma força; não é possível determinar a variável threshold para cada bloco, mas apenas se a variável threshold está acima ou abaixo da particular força aplicada ao bloco. Com estas considerações, temos: Yi 1 sempre que Ti X i Yi 0 sempre que Ti X i (Quebra) (Não Quebra) Segue-se para uma dada força Xi aplicada a um bloco selecionada ao acaso: i P(Yi 1 | X i ) P(Ti X i ) 13 A P(T X) é a distribuição de probabilidade acumulada da variável threshold de todos os blocos na população. Considerando esta distribuição como sendo a logística temos: exp( 0 1 X ) 1 exp( 0 1 X ) P(T X ) encontramos a função resposta logística (5). Outro exemplo de variável threshold: tolerância dos insetos a um inseticida. Uma função de resposta curvilínea com a mesma forma da função logística (5), é obtida transformando por meio da distribuição normal acumulada. Esta transformação é chamada de transformação probit. O modelo de regressão probit é menos flexível do que a regressão logística pois não pode ser diretamente aplicada com mais de uma variável preditora. A distribuição de probabilidade acumulada é dada por: PT X 0 1 X outra função de resposta curvilínea é a transformação complemento log log da probabilidade dada por: loge ( loge (1 )) Diferentemente das transformações logit e probit esta transformação não é simétrica em torno de =0,5. 14 Regressão logística com uma única variável preditora Modelo de regressão logística simples Quando a variável resposta é binária, tomando os valores 1 e 0, com probabilidades e 1-, respectivamente, Y é uma variável Bernoulli com parâmetro E(Y)= . O modelo na sua forma usual é dado por: Yi E (Yi ) i onde: 0 1 X i ) E (Yi ) i 1exp( exp( 0 1 X i ) (8) Função de verossimilhança Temos: P (Yi 1) i P (Yi 0) 1 i 15 A distribuição de probabilidade (Bernoulli) é dada por: f i (Yi ) iYi (1 i )1Yi Yi 0,1; i 1,2,.., n (9) Como as observações Yi são independentes, a conjunta fica: n n g (Y1 ,..., Yn ) f i (Yi ) iYi (1 i )1Yi i 1 (10) i 1 Aplicando o logaritmo, fica: n loge g Y1,...,Yn i 1 log (1 ) Yi loge 1i i n i 1 e i (11) Sabemos que E(Yi)=i para uma variável binária e, de (5) temos: 1 (1 exp( 0 1 X i )) 1 (12) Além disso, considerando (6) e (7), a função de verossimilhança é dada por: n n loge L( 0 , 1 ) Yi ( 0 1 X i ) loge (1 exp(0 1 X i )) i 1 (13) i 1 16 Estimadores de máxima verossimilhança Novamente, não existe uma solução analítica para os valores 0 e 1 que maximizam a função de verossimilhança (13). Métodos numéricos são necessários para encontrar as estimativas de máxima verossimilhança, b0 e b1. Encontradas as estimativas b0 e b1, substitui-se esses valores em (8) para encontrar os valores ajustados. O valor ajustado para o i-ésimo valor é dado por: exp( b0 b1 X i ) ˆ i 1exp( b0 b1X i ) (14) A função de resposta ajustada é dada por: b0 b1 X ) ˆ 1exp( exp( b0 b1 X ) (15) Se usarmos a transformação logit (6), a função resposta ajustada é dada por: ˆ ' b0 b1 X onde: ˆ ' loge 1ˆˆ (16) (17) 17 Exemplo. Um analista está estudando o efeito do tempo de experiência em programação computacional sobre a habilidade para completar, dentro de um determinado tempo, um tarefa difícil. Vinte e cinco (25) programadores foram selecionadas para o estudo. A variável preditora, X, corresponde ao meses de experiência. Os resultados foram (experiência em meses, sucesso na tarefa, valores ajustados): 1 4 2 9 6 2 5 1 8 4 1 8 1 2 2 2 6 3 0 1 1 3 0 0 0 0 1 1 0 0 0 1 0 1 0 1 0 . 3 1 0 2 6 2 0 . 8 3 5 2 6 3 0 . 1 0 9 9 9 6 0 . 7 2 6 6 0 2 0 . 4 6 1 8 3 7 0 . 0 8 2 1 3 0 0 . 4 6 1 8 3 7 0 . 2 4 5 6 6 6 0 . 6 2 0 8 1 2 0 . 1 0 9 9 9 6 0 . 8 5 6 2 9 9 0 . 2 1 6 9 8 0 0 . 8 5 6 2 9 9 500 . 0 9 5 1 5 4 2 010 . 5 4 2 4 0 4 1 300 . 2 7 6 8 0 2 900 . 1 6 7 1 0 0 3 210 . 8 9 1 6 6 4 2 400 . 6 9 3 3 7 9 1 310 . 2 7 6 8 0 2 1 900 . 5 0 2 1 3 4 400 . 0 8 2 1 3 0 2 810 . 8 1 1 8 2 5 2 210 . 6 2 0 8 1 2 810 . 1 4 5 8 1 5 18 Os tempos de experiência são bastante variados, como mostra a primeira coluna dos dados. Para todas as pessoas foi dada a mesma tarefa e os resultados do seu sucesso é mostrado na segunda coluna. Os resultados são codificados como: Y=1 se a tarefa foi completada com sucesso no tempo permitido, e Y=0 se a tarefa não foi completada com sucesso. O diagrama de dispersão é dado na figura a seguir. Somente indica que a habilidade para completar a tarefa com sucesso parece aumentar com a experiência. 19 Ajuste do modelo de regressão logístico (8). Os resultados da análise foram obtidos usando o SAS. Variable INTERCEPT EXPERIE DF 1 1 Parameter Estimates Estimate Std Error -3.0597 1.2594 0.1615 0.0650 Chi-Sq 5.9029 6.1760 Pr > Chi-Sq 0.0151 0.0129 A função resposta logística ajustada (15) é: 3, 0597 0 ,1615 X ) ˆ 1exp( exp( 3, 0597 0 ,1615 X ) Parameter Estimates Variable DF Estimate Std Error Chi-Sq TE RCE T 1 coluna -3 597dos dados. 1.2 594 Os valores ajustadosIN são dados naP terceira da .0 matriz Exercício: EXPERIE 1 .1615 .0650 obtenha a resposta média ajustada para i=1onde X10 =14. Interprete o 0 resultado. Resp. 0,310. 5.902 6.176 Interpretação: este valor ajustado é a estimativa da probabilidade de que uma pessoa com 14 meses de experiência tenha sucesso para completar a tarefa. 20 1,0 0,9 0,8 0,7 ˆ 0,6 0,5 0,4 0,3 0,2 0,1 0,0 0 5 10 15 20 25 30 35 Tempo de experiência 21 Interpretação de b1 Considere o valor da função resposta ajustada (16) em X=Xj ˆ ' ( X j ) b0 b1 X j (18) Por exemplo, para X1=14, temos: ' ˆ (14) 3,059 0,1615(14) 0,798 Considere, também, o valor da função resposta ajustada (16) para X=Xj+1 ˆ ' ( X j 1) b0 b1( X j 1) (19) Por exemplo, para X1=15, temos: ˆ ' (15) 3,059 0,1615(15) 0,6365 22 A diferença entre os dois valores fica: ˆ ' ( X j 1) ˆ ' ( X j ) b1 (20) De acordo com (17), (18) é o logaritmo da chance (odds) estimada quando X=Xj e denominamos por loge(chance1). Da mesma forma, (19) é o logaritmo da chance(odds) estimada quando X=Xj+1 e denominamos por loge(chance2). Observação: odds 1ˆˆ Assim, a diferença entre os dois valores ajustados pode ser dado por: loge (chance2 ) loge (chance1 ) log chance2 e chance1 b 1 No exemplo, chance1 e chance2, valem: ˆ ˆ chance1 expln 0,4502 1 ˆ 1 ˆ ˆ ˆ chance2 expln 0,5291 1 ˆ 1 ˆ 23 Aplicando o anti-logaritmo em cada lado, vemos que a razão das chances estimada, denominada de razão das chances (odds ratio), é dada por: ^ OR chance2 chance1 exp(b1 ) (21) Exemplo: para os dados da tarefa computacional, o valor da razão das chances é: ^ OR exp(0,1615) 1,175 Interpretação: a chance aumenta em 17,5% para cada mês adicional de experiência. Da mesma forma: chance2 0,5291 OR 1,1753 chance1 0,4502 24 Em geral, a razão das chances estimada quando existe uma diferença de c unidades em X é exp(cb1). Exemplo: desejamos comparar indivíduos com 10 meses e 25 meses de experiência, assim c=15 meses, a razão das chances é estimada por exp(15*0,1615)=11,3, portanto, isto indica que a chance de uma pessoa com experiência terminar a tarefa aumenta mais de 11 vezes quando comparado com uma pessoa com pouca experiência. Regressão logística com várias variáveis preditoras Modelo No modelo (8) substituímos: 0 1 X por: 0 1 X1 ... p1 X p1 (22) 25 Em termos matriciais, temos: 0 1 . β ( px1) . . p 1 1 X 1 X2 X (px1) . . X ( p 1) 1 X i1 X i2 Xi . (px1) . X i ,( p 1) (23) Temos que: β X 0 1 X1 ... p1 X p1 (24) ' β 'Xi 0 1 X i1 ... p1 X i, p1 (25) 26 A função (5) pode ser generalizada como: E (Y ) exp(β ' X ) 1 exp(β ' X ) (26) Uma forma equivalente é dada por: 1 E(Y ) (1 exp(β X)) ' (27) A transformação logit dada em (6) agora resulta em: ' β'X (28) Formulação do modelo: Sejam Yi variáveis aleatórias independentemente distribuídas segundo uma Bernoulli com valores esperados E(Yi)=i, onde: E (Yi ) i 1exp(β 'X ) exp(β ' Xi ) i (29) 27 As variáveis X podem ser variáveis preditoras diferentes, ou algumas podem representar efeitos de curvatura e/ou interação. Também, as variáveis preditoras podem ser quantitativas, ou elas podem ser qualitativas e representadas por variáveis indicadoras. Esta flexibilidade torna o modelo de regressão logístico múltiplo muito atrativo. Ajustando o modelo A função log-verossimlhança (13) estende-se diretamente para o modelo de regressão logística múltipla, dada por: n n loge L(β ) Yi (β ' Xi ) loge (1 exp(β ' Xi )) i 1 (30) i 1 Métodos numéricos devem ser utilizados para encontrar os valores de 0, 1,..., p-1 que maximizam (30). As estimativas de máxima verossimilhança serão denotadas por b0, b1,...,bp-1. 28 A função resposta logística ajustada e os valores ajustados são dados por: ˆ exp(b' X ) 1 exp(b' X ) (1 exp(b X) ) ' -1 ' -1 ˆ i 1exp(b'X ) (1 exp(b Xi ) ) exp(b' Xi ) (31) (32) i Exemplo: um estudo na área da saúde está investigando um surto epidêmico de uma doença transmitida por um mosquito, indivíduos foram aleatoriamente selecionados em dois setores de uma cidade para determinar se a pessoa tinha recentemente contraído a doença em estudo. Isto foi verificado por um entrevistador, que fez certas questões específicas para saber se o entrevistado apresentou sintomas da doença durante um período específico. A variável resposta Y foi codificada como 1 se a doença estava presente, e 0 em caso contrário. Três variáveis preditoras foram incluídas no estudo: idade, status sócio-econômico da família e o setor da cidade. A idade (X1) é uma variável quantitativa. O status sócioeconômico é uma variável com 3 categorias. Esta variável é representada por duas variáveis indicadoras (X2 e X3) do seguinte modo: 29 Classe X2 X3 Alta 0 0 Média 1 0 Baixa 0 1 A variável setor da cidade também é uma variável categorizada. Como existiam apenas dois setores na cidade, uma variável indicadora (X4) foi usada, definida como X4=0 para o setor 1 e X4=1 para o setor 2. A razão para a escolha da classe social alta ser tomada como referência é que é esperado que esta classe teria a menor taxa de doença entre as classes sociais. Fazendo-se esta classe como referência, a razão das chances associados com os coeficientes de regressão 2 e 3 espera-se serem maiores do que 1, facilitando a interpretação. Pela mesma razão, o setor 1, onde a epidemia foi menos severa, foi escolhida como referência para a variável indicadora X4. 30 Parte dos dados: Observação (1) Idade i 1 2 3 4 5 6 ... 98 Xi1 33 35 6 60 18 26 ... 35 (2) (3) Status sócio-econômico Xi2 0 0 0 0 0 0 ... 0 Xi3 0 0 0 0 1 1 ... 1 (4) Setor Xi4 0 0 0 0 0 0 ... 0 (5) Status doença Yi 0 0 0 0 1 0 ... 0 (60 Valores ajustados ˆ i .209 .219 .106 .371 .111 .136 ... .171 O primeiro propósito da análise foi verificar a força de associação entre as variáveis preditoras e a probabilidade de uma pessoa ter contraído a doença. Modelo ajustado Onde: E(Y ) (1 exp(β 'X))1 β 'X 0 1 X1 2 X 2 3 X 3 4 X 4 31 Estimativas de máxima verossimilhança da função de regressão logística - dados de doenças Estimativas dos coeficientes, desvios padrões e razão das chances Coeficientes de Estimativas dos Estimativas dos Estimativas das regressão coeficientes de desvios padrões razões das chances regressão -2,3129 0,6426 0 0,02975 0,0135 1,030 1 0,4088 0,5990 1,505 2 -0,30525 0,6041 0,737 3 1,5747 0,5016 4,829 4 Odds ratio da classe alta: e-2- 3=e-0,4+0,31=0,899. Interpretar. Estimativa da matriz de variâncias-covariâncias aproximadas 0,4129 0,0057 0,1836 0,2010 0,1632 0,00018 0,00115 0.00073 0,00034 2 s (b) 0,3588 0,1482 0,0129 0 , 3650 0 , 0623 0,2516 32 A função resposta logística estimada é dada por: ˆ (1 exp(2,3129 0,02975X1 0,4088X 2 0,30525X 3 1,5747X 4 ))1 Interpretação das razões das chances (odds ratios): ^ OR 1,03 A chance de uma pessoa estar doente aumenta cerca de 3% com cada ano adicional de idade (X1), para dado status sócioeconômico e setor da cidade (constantes). ^ OR 4,829 A chance de uma pessoa no setor 2 (X4) que tenha contraído a doença é quase 5 vezes maior para uma pessoa do setor 1, dado a idade e o status sócio-econômico. Exercício: encontre o valor ajustado para o caso i=1, onde X11=33, X12=0, X13=0, X14=0. Resposta: ˆ1 0,209 Interpretação: é a estimativa da probabilidade de uma pessoa com 33 anos de idade, da classe alta, do setor 1, contrair a doença. 33 Construção de Modelos: Seleção de Variáveis Preditoras Nesta seção vamos considerar o processo de seleção de variáveis explanatórias via o método passo a passo (stepwise), e a validação do modelo de regressão logístico. Método passo a passo (stepwise) para construção do modelo Usa-se o procedimento de regressão stepwise para adicionar ou remover variáveis explanatórias do modelo, assim como efeitos de curvatura e de interação. O método é idêntico ao modelo de regressão linear. 34 Teste se vários k=0 Aqui, o nosso interesse é verificar se um subconjunto das variáveis X podem ser retiradas do modelo de regressão logística múltiplo, isto é, vamos testar se os coeficientes de regressão k são iguais a zero. Para este fim nós vamos usar o Teste da Razão de Verossimilhança, que é baseado na estatística chamada de Deviance do modelo. Deviance do modelo: Definição: a deviance (desvio) de um modelo de pesquisa compara o logaritmo da verossimilhança deste modelo com o logaritmo da verossimilhança do modelo completo. Um modelo completo é um modelo que se ajusta completamente aos dados, isto é, para cada observação tem-se um parâmetro. A deviance, para o modelo de regressão logístico (29), é dada por: n DEV ( X 0 , X 1 ,..., X p 1 ) 2 Yi loge (ˆ i ) 1 Yi loge 1 ˆ i i 1 onde ˆi é o i-ésimo valor ajustado do modelo de regressão logístico (32). 35 Deviance pequena Deviance grande A explicação do modelo ajustado (de pesquisa) é praticamente igual ao do modelo completo, ou seja, podemos usar o modelo ajustado (pesquisa), pois, geralmente tem menos parâmetros, ele é mais simples. A explicação do modelo ajustado (de pesquisa) é pobre, ou seja, não podemos usar o modelo ajustado (pesquisa). Deviance Parcial Para cada modelo ajustado (ou de pesquisa) podemos calcular a sua deviance (desvio). A diferença entre as deviances de dois modelos de pesquisa é denominada de deviance parcial e, através dela, é possível testarmos se determinada(s) variável(eis) explanatória(s) pode(m) ser retirada(s) do modelo. 36 A seguir mostraremos o processo de teste usando a deviance parcial. Vamos considerar o modelo logístico completo com função resposta dada por: 1 exp β X ' C 1 (33) β X 0 1 X1 ... p1 X p1 ' C Calcula-se as estimativas de máxima verossimilhança (bC) e a deviance deste modelo, a qual é representada por: DEV(X0, X1, ...,Xp-1). Vamos, agora, considerar que desejamos testar as seguintes hipóteses: H 0 : q q 1 ... p 1 0 H a : pelo menosum βk é diferentede zero Os p-q coeficientes são testados. 37 O modelo de regressão logístico reduzido tem a seguinte função resposta: 1 exp β X ' R 1 (33) β X 0 1 X1 ... q1 X q1 ' R Calcula-se as estimativas de máxima verossimilhança (bR) e a deviance deste modelo, a qual é representada por: DEV(X0, X1, ...,Xq-1). Interpretação: Se a deviance (residual) do modelo reduzido não é muito maior do que a deviance (residual) do modelo completo, a nossa conclusão é que as variáveis Xq, Xq+1,...,Xp-1, podem ser retiradas do modelo de regressão logístico múltiplo. Uma grande diferença entre as duas deviances (residuais) significa que as variáveis preditoras Xq, Xq+1,...,Xp-1, devem ser mantidas no modelo, pois elas melhoram muito o ajuste do modelo (a explicação do modelo). 38 A diferença entre as duas deviances é a deviance parcial e é dada por: DEV X q , X q 1 ,..., X p 1 | X 0 , X 1 ,..., X q 1 DEV X 0 , X 1 ,..., X q 1 - DEV X 0 , X 1 ,..., X p 1 (34) A deviance parcial dada em (34) segue, aproximadamente, para um n razoavelmente grande, uma distribuição de qui-quadrado com p-q graus de liberdade. Os graus de liberdade correspondem a diferença nos graus de liberdade do erro para os dois modelos ajustados: (n-q)-(n-p)=p-q. Regra de decisão usando a aproximação pelo Qui-Quadrado DEV X q , X q 1 ,..., X p 1 | X 0 , X 1 ,..., X q 1 2 1 ; p q não rejeitarH 0 DEV X q , X q 1 ,..., X p 1 | X 0 , X 1 ,..., X q 1 2 1 ; p q rejeitarH 0 39 Ilustração do uso da deviance parcial: Vamos considerar um modelo de regressão logístico com: β' X 0 1 X1 2 X 2 3X3 e a sua deviance residual é calculada. Hipóteses em teste: H 0 : 2 3 0 H a : pelo menosum βk é diferentede zero Então, vamos ajustar um modelo de regressão logístico com: β' X 0 1 X1 e vamos obter a deviance residual deste modelo. A deviance parcial necessária para verificar as hipóteses, é dada por: DEV X 2 , X 3 | X 0 , X1 DEV X 0 , X1 DEV X 0 , X1, X 2 , X 3 40 Exemplo: continuação do exemplo de surto de uma doença. O modelo foi ajustado com três variáveis explanatórias: idade, classe sócioeconômica e setor da cidade. A deviance para este modelo é dada por: DEV X 0 , X1 , X 2 , X 3 , X 4 101,054 H 0 : 1 0 Hipótese: H a : β1 0 Assim, vamos ajustar um modelo com: β' X 0 2 X 2 3X3 4 X 4 e a sua deviance vale: DEV X 0 , X 2 , X 3 , X 4 106,204 41 A deviance parcial vale: DEV X 1 | X 0 , X 2 , X 3 , X 4 DEV X 0 , X 2 , X 3 , X 4 DEV X 0 , X 1 , X 2 , X 3 , X 4 DEV X 1 | X 0 , X 2 , X 3 , X 4 106,204 101,054 5,15 Temos que 2(0,05,1)=3,81. Como 5,15>3,81, rejeitamos a hipótese nula e, portanto, a variável X1 deve permanecer no modelo. 42 Vamos, usando o pacote estatístico SAS testar a seguinte hipótese: H 0 : 2 3 0 H a : pelo menosum k 0 Model Fit Statistics Intercept and Criterion -2 Log L Covariates 101.054 (Deviance residual do modelo completo) Model Fit Statistics Intercept and Criterion -2 Log L Covariates 102.259 (Deviance residual do modelo de pesquisas) 43 Residual Chi-Square Test Chi-Square DF Pr > ChiSq 1.2213 2 0.5430 De acordo com o teste podemos retirar do modelo as variáveis X2 e X3 (sócioeconômicas). Ainda podemos considerar um modelo com as interações: β F X 0 1 X 1 2 X 2 3X3 4 X 4 5 X 1 X 2 6 X1X3 ' 7 X 1 X 4 8 X 2 X 4 9 X 3 X 4 Vamos, usando o pacote estatístico SAS testar a seguinte hipótese: H 0 : 5 6 7 8 9 0 H a : pelo m enosum k 0 βR X 0 1 X1 2 X 2 3X3 4 X 4 ' 44 The LOGISTIC Procedure Residual Chi-Square Test Chi-Square 6.6450 DF 5 Pr > ChiSq 0.2484 A conclusão é de que não é necessário a inclusão de interações de primeira ordem no modelo de regressão logístico. Validação do modelo Novos dados (uma nova amostra) ou, então, uma amostra reservada dos dados, deveria ser usada para verificar se o mesmo modelo pode ser usado com estes dados novos, se os coeficientes de regressão e os erros padrões são similares, e se as mesmas conclusões inferenciais seriam obtidas. 45 Diagnóstico do Modelo Verificação do ajuste do modelo Verificar o ajuste da parte linear do modelo de regressão logístico e identificar deviance residual que são valores extremos (outlying) Observação: observações outlying são observações bem separadas dos resto dos dados. Geralmente são identificadas com resíduos grandes. Elas tem um efeito muito grande sobre a função de regressão de mínimos quadrados ajustada. Pontos cruciais: Verificar se a função resposta estimada é monotônica e em forma sigmoidal (de S) Verificar a presença de outliers, pontos influentes e se o modelo de regressão logístico ajustado é adequado. 46 Verificação do Ajuste do Modelo Procedimento: ˆ ' b' X i. Criar classes com valores similares de escala logito) ; (valores ajustados na ii. Estas classes devem ter aproximadamente o mesmo número de casos (sugestão: usar classes igualmente espaçadas); iii. Para cada classe computar a proporção de 1’s, representada por pj; iv. Fazer um gráfico das proporções versus os pontos médios das classes dos valores de ˆ ' b' X Exemplo: continuação do exemplo de uma tarefa de programação. Os valores estimados são dados por: ˆ ´ 3,0597 0,1615* experiência Estes valores, ordenados crescentemente, são dados na tabela. 47 Obs experiencia sucesso valor ajustado 1 4 0 -2.4137 2 4 0 -2.4137 3 5 0 -2.2522 4 6 0 -2.0907 5 6 0 -2.0907 6 8 1 -1.7677 7 9 0 -1.6062 8 11 0 -1.2832 9 12 0 -1.1217 10 13 0 -0.9602 11 13 1 -0.9602 12 14 0 -0.7987 13 18 1 -0.1527 14 18 0 -0.1527 15 19 0 0.0088 16 20 1 0.1703 17 22 1 0.4933 18 22 1 0.4933 19 24 0 0.8163 20 25 1 0.9778 21 28 1 1.4623 22 29 0 1.6238 23 30 1 1.7853 24 30 1 1.7853 25 32 1 2.1083 48 Classes (j) ˆ ' ˆ ' Ponto médio nj pj Freqüências de uns 1 -2,42 ├ -1,514 -1,967 7 0,143 1 2 -1,514 ├ -0,608 -1,061 5 0,200 1 3 -0,608 ├ 0,298 -0,155 4 0,500 2 4 0,298 ├ 1,204 0,751 4 0,750 3 5 1,204 ├ 2,11 1,657 5 0,800 4 Observação: as classes são igualmente espaçadas. 49 Conclusão: a função resposta é monotônica e em forma de uma sigmóide. 50 Gráfico half-normal de probabilidade com envelope simulado É útil para verificar se algum(ns) valor(es) da deviance residual é discrepante(outlying), e para verificar se a parte linear do modelo de regressão logístico é adequada. Num gráfico half-normal de probabilidades o k-ésimo resíduo, em valor absoluto, ordenado é colocado num gráfico com o seguinte percentil: k n 1/ 8 z 2n 1 / 2 (35) (é o valor de z que dá uma área acumulada de (k+n-1/8)/(2n+1/2)) Outliers aparecerão no alto, à direita do gráfico, como pontos separados dos outros. O envelope simulado é uma faixa, cujos resíduos devem cair dentro desta faixa se o modelo é adequado (ajustado, correto). 51 • • • • • • Passos para construir o gráfico half-normal de probabilidades com envelope simulado Para cada uma das n observações, gerar um experimento Bernoulli (0 ou 1), onde o parâmetro da Bernoulli para a i-ésima observação é ˆ , a i probabilidade estimada da resposta Yi=1de acordo com o modelo ajustado originalmente; Ajustar um modelo de regressão logístico para as n novas respostas onde a variável preditora mantém seus valores originais, e obtenha as deviances residuais. Ordenar as deviances residuais tomadas em valor absoluto em ordem crescente. Repetir os dois primeiros passos 18 vezes; Agrupe as menores deviances residuais absolutas a partir dos 19 grupos e determine os valores mínimo, máximo e médio desses 19 resíduos; Repita o passo anterior agrupando os segundos menores resíduos absolutos, depois os terceiros menores resíduos absolutos, e assim por diante. Represente os valores mínimo, médio e máximo de cada um dos n grupos de resíduos ordenados versus o correspondente valor esperado em (35) em um gráfico half-normal de probabilidades para os valores de deviances residuais absolutas ordenadas da amostra original e ligue os pontos por linhas retas. 52 Exemplo: continuação do exemplo de uma tarefa de programação. 53 Variáveis: obs experiência sucesso probabilidade_estimada simulação1...simulação19 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 14 29 6 25 18 4 18 12 22 6 30 11 30 5 20 13 9 32 24 13 19 4 28 22 8 0 0 0 1 1 0 0 0 1 0 1 0 1 0 1 0 0 1 0 1 0 0 1 1 1 0.31026 0.83526 0.11000 0.72660 0.46184 0.08213 0.46184 0.24567 0.62081 0.11000 0.85630 0.21698 0.85630 0.09515 0.54240 0.27680 0.16710 0.89166 0.69338 0.27680 0.50213 0.08213 0.81183 0.62081 0.14582 0 1 0 0 0 0 1 0 0 0 1 0 1 0 1 0 0 0 0 1 1 0 1 1 0 0 0 0 0 1 0 1 0 1 0 1 1 1 0 1 1 0 1 1 0 0 0 1 1 0 0 1 0 1 0 0 0 1 0 1 1 0 1 0 1 0 0 1 1 0 0 0 1 1 1 0 1 0 0 0 0 0 0 1 0 1 0 1 0 1 0 0 1 1 1 1 1 1 1 0 1 1 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 1 0 0 0 0 1 1 1 1 1 0 1 0 1 1 0 1 0 1 0 1 0 1 0 0 1 1 0 1 1 1 1 0 0 1 0 0 0 0 1 0 1 0 1 1 1 0 1 1 0 1 0 0 0 0 1 1 0 0 1 0 1 1 0 1 0 1 0 1 0 1 0 0 1 0 1 1 0 0 0 1 1 0 0 0 0 0 1 0 1 0 1 1 1 1 1 0 0 0 0 1 0 0 1 0 1 1 0 0 0 0 1 1 0 1 0 1 0 1 0 1 0 1 1 0 1 1 0 0 0 1 1 0 0 1 0 1 1 1 0 0 1 0 1 0 1 1 0 1 0 1 0 0 0 0 1 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 0 1 0 0 1 1 0 1 0 1 1 0 0 1 0 1 0 0 1 1 1 0 1 1 0 0 1 0 0 1 1 0 0 0 1 0 0 1 1 0 1 1 0 0 1 0 0 1 1 1 0 1 0 0 1 1 0 1 0 1 1 0 0 1 1 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 0 1 1 0 0 1 1 1 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 1 1 0 1 0 0 1 0 1 1 0 0 0 1 0 1 0 0 0 0 1 1 1 1 0 1 0 0 1 1 0 1 1 0 0 1 0 0 1 1 1 1 1 0 0 1 1 0 1 0 1 0 0 54 0 1 0 1 1 0 0 0 1 0 1 0 1 1 1 0 0 1 1 1 1 0 1 1 0 Deviannce residual ordenadas crescentemente para as amostras simuladas desvio1 desvio2 desvio3 desvio4 desvio5 desvio6 desvio7 desvio8 desvio9 desvio10 desvio11 desvio12 desvio13 desvio14 desvio15 desvio16 desvio17 desvio18 desvio19 0,34 0,39 0,42 0,33 0,47 0,20 0,33 0,08 0,55 0,21 0,58 0,02 0,38 0,14 0,28 0,16 0,54 0,35 0,11 0,34 0,43 0,46 0,37 0,50 0,25 0,33 0,11 0,55 0,25 0,63 0,04 0,38 0,19 0,28 0,16 0,57 0,41 0,15 0,37 0,43 0,46 0,40 0,50 0,25 0,36 0,11 0,58 0,25 0,63 0,04 0,41 0,19 0,31 0,18 0,60 0,41 0,15 0,40 0,47 0,49 0,40 0,54 0,28 0,38 0,12 0,61 0,27 0,67 0,05 0,44 0,21 0,34 0,20 0,60 0,44 0,17 0,40 0,47 0,49 0,41 0,54 0,31 0,40 0,12 0,64 0,27 0,70 0,07 0,45 0,24 0,42 0,25 0,62 0,47 0,20 0,46 0,47 0,49 0,44 0,54 0,41 0,40 0,13 0,68 0,28 0,71 0,07 0,45 0,34 0,43 0,28 0,67 0,59 0,29 0,50 0,51 0,53 0,45 0,58 0,46 0,46 0,14 0,71 0,32 0,77 0,07 0,52 0,34 0,46 0,35 0,70 0,59 0,29 0,58 0,51 0,53 0,45 0,58 0,55 0,46 0,16 0,71 0,32 0,77 0,09 0,53 0,36 0,52 0,39 0,70 0,59 0,30 0,62 0,55 0,57 0,48 0,58 0,55 0,48 0,16 0,72 0,34 0,80 0,12 0,56 0,38 0,52 0,44 0,71 0,63 0,35 0,66 0,60 0,67 0,54 0,62 0,57 0,51 0,17 0,78 0,41 0,89 0,12 0,57 0,41 0,55 0,44 0,77 0,68 0,38 0,71 0,65 0,71 0,59 0,80 0,63 0,53 0,23 0,84 0,47 0,92 0,14 0,61 0,44 0,57 0,47 0,79 0,68 0,38 0,72 0,75 0,77 0,69 0,86 0,63 0,55 0,28 0,88 0,50 0,92 0,18 0,77 0,44 0,60 0,49 0,83 0,78 0,47 0,72 0,81 0,77 0,71 0,91 0,67 0,69 0,30 0,88 0,56 0,97 0,19 0,78 0,52 0,62 0,58 0,87 0,83 0,47 0,77 0,87 0,88 0,77 0,91 0,73 0,75 0,37 0,92 0,59 1,01 0,25 0,78 0,56 0,66 0,58 0,90 0,84 0,51 0,82 0,87 0,89 0,82 0,91 0,76 0,82 0,41 1,05 0,67 1,05 0,30 0,83 0,63 0,72 0,65 0,95 0,89 0,58 0,92 0,87 0,89 0,82 0,91 0,80 0,92 0,49 1,05 0,72 1,08 0,30 0,84 0,66 0,81 0,72 1,11 0,95 0,61 1,17 0,94 0,95 0,84 1,03 0,83 0,92 0,53 1,20 0,72 1,09 0,40 0,96 0,75 1,00 0,75 1,19 1,08 0,70 1,18 1,00 1,01 0,91 1,15 0,98 1,08 0,53 1,21 0,75 1,27 0,49 1,09 0,84 1,00 0,75 1,21 1,08 0,76 1,31 1,14 1,21 0,97 1,21 1,06 1,11 0,59 1,25 0,84 1,31 0,62 1,11 0,99 1,03 0,83 1,25 1,21 0,80 1,32 1,14 1,21 1,04 1,27 1,12 1,20 0,71 1,25 0,89 1,32 0,65 1,19 0,99 1,03 0,92 1,32 1,35 0,86 1,37 1,29 1,29 1,23 1,40 1,15 1,24 1,05 1,42 1,10 1,37 0,65 1,24 1,26 1,08 0,96 1,32 1,42 0,97 1,39 1,52 1,50 1,23 1,60 1,15 1,63 1,05 1,47 1,10 1,55 0,77 1,42 1,49 1,16 1,10 1,46 1,42 1,09 1,45 1,67 1,57 1,56 1,67 1,61 1,68 1,48 1,62 1,37 1,56 0,81 1,71 1,56 1,51 1,25 1,52 1,70 1,40 1,80 1,76 1,87 1,85 1,74 2,03 1,72 1,66 1,68 1,68 1,69 1,65 1,79 1,61 1,81 1,30 1,70 1,70 1,61 1,86 2,06 2,01 2,33 1,80 2,03 1,85 1,91 1,88 2,49 1,74 2,05 2,04 2,03 2,39 2,80 2,00 1,85 2,42 55 percentil mínimo médio máximo valor original 0,031027 0,02 0,31 0,58 0,408 0,080746 0,04 0,34 0,63 0,414 0,130666 0,04 0,35 0,63 0,414 0,180913 0,05 0,37 0,67 0,447 0,231622 0,07 0,39 0,70 0,48 0,282934 0,07 0,43 0,71 0,483 0,335002 0,07 0,46 0,77 0,557 0,387995 0,09 0,48 0,77 0,557 0,442101 0,12 0,50 0,80 0,605 0,497535 0,12 0,54 0,89 0,646 0,554542 0,14 0,59 0,92 0,699 0,613411 0,18 0,64 0,92 0,751 0,67449 0,19 0,68 0,97 0,8 0,738194 0,25 0,72 1,01 0,805 0,805048 0,30 0,77 1,05 0,86 0,875709 0,30 0,82 1,11 0,976 0,95104 0,40 0,90 1,20 0,976 1,032197 0,49 0,96 1,27 1,106 1,120793 0,59 1,04 1,31 1,113 1,219191 0,65 1,09 1,35 1,181 1,331064 0,65 1,20 1,42 1,24 1,462645 0,77 1,32 1,63 1,538 1,625949 0,81 1,51 1,71 1,603 1,849703 1,30 1,72 2,03 1,9 Observação 2 1,74 2,08 2,80 1,962 Observação 25 2,245251 56 Envelope simulado deviance residual 3,00 2,50 2,00 1,50 1,00 0,50 0,00 0 0,5 1 1,5 2 2,5 percentil mínimo médio máximo Deviance 57 INTERPRETAÇÃO: Todos os pontos caem dentro do envelope simulado, assim, não há necessidade de se aplicar medidas remediadoras, contra outliers, por exemplo. Além disso, a maioria das deviances residuais estão próximas da linha média, indicando que o modelo de regressão logística é adequado aos dados. As observações 2 e 25 estão no lado direito alto, mas estão dentro da faixa. 58