class: center, middle, inverse, title-slide .title[ # MODELOS LINEARES GENERALIZADOS STA13829 ] .subtitle[ ## Modelos para Dados de Contagem ] .author[ ### Nátaly A. Jiménez Monroy ] .institute[ ### LECON/DEST - UFES ] --- [//]: <> (https://pkg.garrickadenbuie.com/extra-awesome-xaringan/intro/index.html#1) [//]: <> (https://pkg.garrickadenbuie.com/xaringanthemer/articles/xaringanthemer.html) [//]: <> (https://www.biostatistics.dk/talks/CopenhagenRuseRs-2019/index.html#1) [//]: <> (https://rstudio-education.github.io/sharing-short-notice/#1) [//]: <> (https://www.kirenz.com/slides/xaringan-demo-slides.html#1) [//]: <> (https://github.com/yihui/xaringan/issues/26) [//]: <> (https://github.com/emitanaka/anicon) [//]: <> (https://github.com/mitchelloharawild/icons) [//]: <> (https://slides.yihui.org/2020-genentech-rmarkdown.html#1) [//]: <> (https://github.com/gadenbuie/xaringanExtra) [//]: <> (class: center, middle, animated, slideInRight) class: animated, slideInRight <style> body {text-align: justify} </style> <!-- Justify text. --> # Modelo de Poisson * Adequado para os casos em que o interesse está na contagem do número de vezes que ocorre um evento em determinado espaço ou intervalo de tempo. -- * A distribuição de Poisson é usada com frequência para modelar esses dados. Se `\(Y\)` representa o número de ocorrências, sua função de probabilidade pode ser escrita como $$ f(y)=\frac{\mu^ye^{-\mu}}{y!},\quad y=0,1,\ldots, $$ onde `\(\mu\)` é o número médio de ocorrências. Pode mostrar-se que `\(\text{E}(Y)=\mu\)` e `\(\text{Var}(Y)=\mu\)`. -- * O parâmetro `\(\mu\)` precisa de uma definição cuidadosa. Frequentemente é definido como uma **<span style="color:orange">taxa</span>**, por exemplo, número de clientes que compram um produto para cada 100 clientes que entram em uma loja. Para acidentes de carro, a taxa pode ser definida de diversas formas: acidentes por 1000 habitantes, acidentes por 1000 habilitados, acidentes por 1000 automóveis, acidentes por 1000 km percorridos, etc. A escala de tempo deve ser incluída na definição, por exemplo, acidentes por 1000 km percorridos por ano. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## O Modelo Sejam `\(Y_1,\ldots,Y_n\)` variáveis aleatórias independentes distribuídas tais que `\(Y_i\sim \text{Po}(\mu_i)\)`, com componente sistemática dada por `\(g(\mu_i)=\eta_i\)`, em que `\(\eta_i=\mathbb{x}_i^t\mathbb{\beta}\)`, `\(\mathbb{x}_i=(x_{i1},\ldots,x_{ip})^t\)` e `\(\mathbb{\beta}=(\beta_1,\ldots,\beta_p)^t\)`. As funções de ligação mais usadas são a logarítmica `\((g(\mu_i)=\log(\mu_i))\)`, raiz quadrada `\((g(\mu_i)=\sqrt\mu_i)\)` e identidade `\((g(\mu_i)=\mu_i)\)`. -- No caso das unidades experimentais serem observadas em tempos distintos `\(t_i, i=1,\ldots,n\)` e for assumido que `\(Y_i\sim \text{Po}(\lambda_it_i)\)`, a componente sistemática do modelo de ligação logarítmica fica dada por $$ \log(\mu_i)=\log(t_i) + \mathbb{x}_i^t\mathbb{\beta}, $$ em que `\(\log(t_i)\)` desempenha o papel de **<span style="color:orange">offset</span>** e isto deve ser considerado no ajuste. -- Os tempos `\(t_i, i=1,\ldots,n\)` também podem ser incluídos como valores da variável explicativa `\(\log(T_i)\)`. Neste caso, a componente sistemática é dada por $$ \log(\mu_i)=\theta \log(t_i) + \mathbb{x}_i^t\mathbb{\beta}, $$ O teste `\(\text{H}_0: \theta=0\)` vs. `\(\text{H}_1: \theta\neq 0\)` verifica se `\(\log(t_i)\)` deve ser considerado no modelo como *offset*. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - I Consideramos os dados apresentados em Neter *et al*. (1996, p. 613) sobre o perfil dos clientes de uma determinada loja, oriundos de 110 áreas de uma cidade. O objetivo do estudo é relacionar o número esperado de clientes em cada área com as variáveis explicativas: número de domicílios (em mil), renda média anual (em mil USD), idade média dos domicílios (em anos), distância ao concorrente mais próximo (em milhas) e distância à loja (em milhas). -- **Unidade experimental**: área. -- **Modelo**: Seja `\(Y_i, i=1,\ldots,110\)` o número de clientes da `\(i-\)`ésima área que foram à loja no período determinado. Assumimos que `\(Y\sim \text{Po}(\mu_i)\)` com componente sistemática dada por -- $$ \log(\mu_i)=\alpha+ \beta_1 \color{red}{\textbf{n.domicilios}_i} + \beta_2 \text{renda}_i + \beta_3 \text{idade}_i + \beta_4 \text{dist.concorrente}_i + \beta_5 \text{dist.loja}_i. $$ > A variável número de domicílios (n.domicilios) deve ser incluída no modelo dado que as áreas têm diferentes números de domicílios. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - II <div class="datatables html-widget html-fill-item" id="htmlwidget-d9216f678575484bc60f" style="width:100%;height:auto;"></div> <script type="application/json" data-for="htmlwidget-d9216f678575484bc60f">{"x":{"filter":"none","vertical":false,"extensions":["Buttons"],"data":[[9,6,28,11,4,4,0,14,16,13,9,14,5,9,9,7,4,26,32,26,11,12,3,15,12,9,14,10,22,8,3,11,2,6,11,10,0,15,9,16,29,6,26,13,0,8,8,21,12,26,3,7,5,9,13,28,10,12,6,11,12,11,9,8,9,8,6,6,15,15,12,10,8,3,10,10,15,6,9,21,13,8,6,8,21,7,19,13,24,7,3,8,9,18,12,8,16,5,8,11,12,14,17,17,6,15,10,6,4,6],[606,641,505,866,599,520,354,483,1034,456,19,530,337,586,1113,525,377,1127,877,1007,657,302,603,556,635,386,1011,925,898,731,584,439,153,1069,443,392,828,159,830,234,1004,643,741,306,180,644,109,809,722,1006,786,1041,524,725,482,666,450,667,921,412,526,523,1066,1001,669,582,872,758,782,551,201,730,738,469,898,780,622,391,531,566,410,719,684,865,1031,862,758,1141,1289,674,683,650,406,966,1103,312,787,416,528,919,482,781,120,693,348,780,752,817,268,519],[41393,23635,55475,64646,31972,41755,46014,34626,85207,33021,39198,38794,30855,28852,120065,32229,36828,90302,51707,89860,60513,42191,28736,49129,29308,26734,57862,70030,46027,32202,32871,29564,46806,59805,42555,36998,85664,21238,47972,33246,45927,58315,69177,40886,44588,47347,31791,42740,59175,48862,54678,59835,51756,34817,29942,68684,64790,58535,42919,40722,42120,28647,61464,70136,34595,30878,39366,61563,38412,41045,23864,38647,58387,37242,38337,68201,41066,40873,54655,49826,29013,78082,57506,47118,72373,67787,40305,50026,98701,58195,47991,63123,39051,114633,55773,43393,61765,33348,44541,40795,55972,33140,19673,36190,25768,53974,71814,54429,34022,52850],[3,18,27,31,7,23,26,1,13,32,22,5,1,7,9,3,15,26,27,55,32,54,41,33,42,14,54,36,44,43,13,18,21,22,53,7,4,4,40,26,24,8,9,27,14,35,9,17,35,48,20,40,39,18,14,25,3,25,13,32,30,43,40,29,38,58,52,31,26,2,43,9,13,40,32,5,46,19,40,1,50,31,51,46,48,1,15,45,8,54,57,15,29,38,44,41,53,48,31,8,9,30,21,6,42,47,1,47,54,43],[3.04,1.95,6.54,1.67,0.72,2.24,0.77,3.51,4.23,3.07,2.96,2.77,1.33,2.98,3.58,1.27,1.92,5.83,5.19,5.03,4.38,3.41,0.34,4.78,2.53,4.99,4.6,4.58,3.03,5.15,1.47,3.67,0.84,2.5,2.62,1.03,1.3,2.98,2.28,3.95,4.9,0.78,6.61,4.53,0.88,2.94,4.37,4.1,2.38,5.04,3.59,1.68,0.57,1.88,3.17,5.78,4.35,2.78,2.48,2.47,4.29,2.69,1.15,2.58,4.06,1.91,0.73,3.08,2.72,3.62,4.8,0.67,2.01,1.42,2.63,4.12,4.48,1.67,2.32,3.06,2.68,2.7,2.13,2.17,6.27,2.1,3.95,2.79,5.87,4.3,1.54,3.17,3.11,6.33,4.58,2.25,5.39,1.48,4.91,2.97,2.91,1.42,2.65,4.7,1.43,4.21,3.13,1.9,1.2,2.92],[6.32,8.890000000000001,2.05,5.81,8.109999999999999,6.81,9.27,7.92,4.4,6.03,6.09,6.08,9.859999999999999,8.640000000000001,5.26,7.56,8.91,1.74,3.66,2.03,8.300000000000001,5.21,8.289999999999999,3.89,6.17,9.699999999999999,3.94,8.66,5.6,9.67,8.02,5.1,9.18,9.43,5.75,7.74,9.66,8.66,9.26,4.61,2.69,6.26,0.87,2.68,9.380000000000001,7.69,9.31,4.75,5.09,2.21,8.52,7.59,9.1,7.96,6.91,2.55,6.03,5.59,7.69,9.43,6.15,7.54,8.25,9.67,8.779999999999999,6.86,8.67,8.33,6.71,7.45,8.74,7.92,6.6,8.369999999999999,9.56,6.69,4.1,6.9,5.69,4.03,7.58,4.89,8.31,9.06,1.75,8.630000000000001,5.58,6.18,2.73,6.4,9.52,9.460000000000001,9.619999999999999,2.22,8.68,6.43,3.37,7.66,9.67,7.79,5.85,5.71,6.25,9.539999999999999,7.11,6.41,5.47,9.9,9.51,8.619999999999999]],"container":"<table class=\"compact\">\n <thead>\n <tr>\n <th>n.clientes<\/th>\n <th>n.domic<\/th>\n <th>renda<\/th>\n <th>idade<\/th>\n <th>dist.concorr<\/th>\n <th>dist.loja<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"dom":"tBp","buttons":[{"extend":"csv","filename":"Stores"},{"extend":"excel","filename":"Stores"},{"extend":"pdf","filename":"Stores"}],"pageLength":12,"columnDefs":[{"className":"dt-right","targets":[0,1,2,3,4,5]},{"name":"n.clientes","targets":0},{"name":"n.domic","targets":1},{"name":"renda","targets":2},{"name":"idade","targets":3},{"name":"dist.concorr","targets":4},{"name":"dist.loja","targets":5}],"order":[],"autoWidth":false,"orderClasses":false,"lengthMenu":[10,12,25,50,100]}},"evals":[],"jsHooks":[]}</script> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - III <img src="Modelos_dados_de_contagem_files/figure-html/unnamed-chunk-2-1.png" width="47%" style="display: block; margin: auto;" /> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - IV **Modelo** ``` ## ## Call: ## glm(formula = n.clientes ~ n.domic + renda + idade + dist.concorr + ## dist.loja, family = poisson(log), data = store) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 2.942e+00 2.072e-01 14.198 < 2e-16 *** ## n.domic 6.058e-04 1.421e-04 4.262 2.02e-05 *** ## renda -1.169e-05 2.112e-06 -5.534 3.13e-08 *** ## idade -3.726e-03 1.782e-03 -2.091 0.0365 * ## dist.concorr 1.684e-01 2.577e-02 6.534 6.39e-11 *** ## dist.loja -1.288e-01 1.620e-02 -7.948 1.89e-15 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 422.22 on 109 degrees of freedom ## Residual deviance: 114.99 on 104 degrees of freedom ## AIC: 571.02 ## ## Number of Fisher Scoring iterations: 4 ``` ``` ## Analysis of Deviance Table ## ## Model: poisson, link: log ## ## Response: n.clientes ## ## Terms added sequentially (first to last) ## ## ## Df Deviance Resid. Df Resid. Dev ## NULL 109 422.22 ## n.domic 1 42.662 108 379.56 ## renda 1 0.807 107 378.75 ## idade 1 0.316 106 378.43 ## dist.concorr 1 195.949 105 182.49 ## dist.loja 1 67.500 104 114.99 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - V **Observações:** * O desvio do modelo foi 114,99 com 104 graus de liberdade, o que resulta em um nível descritivo de 0,22 indicando um bom ajuste do modelo. -- * O número esperado de clientes na loja cresce com o aumento do número de domicílios na área e da distância ao concorrente mais próximo. -- * O número esperado de clientes na loja diminui com o aumento da renda média e da idade média dos domicílios, bem como da distância da área à loja. -- * Se a distância ao concorrente mais próximo aumentar uma milha, esperamos um aumento relativo no número de clientes de `\(\exp(0,168) = 1,2\)`. Ou seja, uma aumento em média de 20%. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - VI **Diagnóstico do modelo: Análise gráfica de resíduos.** Código diponível em https://www.ime.usp.br/~giapaula/diag_pois. <img src="Modelos_dados_de_contagem_files/figure-html/unnamed-chunk-4-1.png" width="40%" style="display: block; margin: auto;" /> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - VII **Diagnóstico do modelo: envelope simulado.** Código disponível em https://www.ime.usp.br/~giapaula/envel_pois. <img src="Modelos_dados_de_contagem_files/figure-html/unnamed-chunk-5-1.png" width="45%" style="display: block; margin: auto;" /> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - VIII >**Observações:** >* Dentre as observações destacadas pelos gráficos de diagnóstico, duas apresentam algumas variações desproporcionais nas estimativas dos parâmetros, porém sem mudança inferencial. -- >* Não há indícios de que a ligação usada seja inapropriada (Figura 4 da análise de resíduos). -- >* O gráfico de envelope não mostra indícios de afastamentos relevantes da suposição de distribuição de Poisson para o número de clientes. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - I Consideramos dados apresentados em Venables & Ripley (1999, Caps 6 e 7) sobre um estudo sociológico desenvolvido na Austrália com 146 estudantes de `\(8^a\)` série e ensino médio com o objetivo de comparar a ausência na escola segundo os fatores: ano que o estudante está cursando (F0: `\(8^a\)` série, F1: `\(1^o\)` ano do ensino médio, F2: `\(2^o\)` ano do ensino médio, F3: `\(3^o\)` ano do ensino médio), etnia (A: aborígine, N: não aborígine), desempenho escolar (SL: insuficiente, AL: suficiente) e sexo (M: masculino, F: feminino). -- **Modelo**: Seja `\(Y_{ijklm}\)`, `\(i=1,\ldots,4\)`; `\(j,k,l=1,2\)` e `\(m=1,\ldots,146\)` o número de faltas em um determinado período referentes ao `\(i-\)`ésimo ano, da etnia `\(j\)`, com desempenho escolar `\(k\)` e pertencente ao `\(l-\)`ésimo sexo. Assumimos que `\(Y\sim \text{Po}(\mu_{ijkl})\)` com componente sistemática dada por -- $$ \log(\mu_{ijklm})=\alpha+ \beta_1 \text{ano}_i + \beta_2 \text{etnia}_j + \beta_3 \text{desempenho}_k + \beta_4 \text{sexo}_l $$ --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - II <div class="datatables html-widget html-fill-item" id="htmlwidget-38b1ecd4faabbb4699cf" style="width:100%;height:auto;"></div> <script type="application/json" data-for="htmlwidget-38b1ecd4faabbb4699cf">{"x":{"filter":"none","vertical":false,"extensions":["Buttons"],"data":[["A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N","N"],["M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","M","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F","F"],["F0","F0","F0","F0","F0","F0","F0","F0","F1","F1","F1","F1","F1","F2","F2","F2","F2","F2","F2","F2","F2","F2","F2","F2","F3","F3","F3","F3","F3","F3","F3","F0","F0","F0","F0","F0","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F2","F2","F2","F2","F2","F2","F2","F2","F2","F3","F3","F3","F3","F3","F3","F3","F3","F3","F0","F0","F0","F0","F0","F0","F0","F0","F0","F1","F1","F1","F1","F1","F1","F1","F1","F1","F2","F2","F2","F2","F2","F2","F2","F2","F2","F2","F3","F3","F3","F3","F3","F3","F3","F0","F0","F0","F0","F0","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F1","F2","F2","F2","F2","F2","F2","F2","F2","F2","F2","F3","F3","F3","F3","F3","F3","F3","F3","F3","F3"],["SL","SL","SL","AL","AL","AL","AL","AL","SL","SL","SL","AL","AL","SL","SL","SL","SL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","SL","AL","AL","AL","AL","SL","SL","SL","SL","SL","SL","SL","SL","SL","SL","AL","AL","AL","AL","AL","SL","SL","SL","SL","SL","SL","SL","SL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","SL","SL","SL","AL","AL","AL","AL","AL","AL","SL","SL","SL","SL","SL","SL","SL","AL","AL","SL","SL","SL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","SL","AL","AL","AL","AL","SL","SL","SL","SL","SL","SL","SL","SL","SL","SL","SL","AL","AL","AL","AL","AL","AL","SL","SL","SL","SL","SL","SL","SL","SL","SL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL","AL"],[2,11,14,5,5,13,20,22,6,6,15,7,14,6,32,53,57,14,16,16,17,40,43,46,8,23,23,28,34,36,38,3,5,11,24,45,5,6,6,9,13,23,25,32,53,54,5,5,11,17,19,8,13,14,20,47,48,60,81,2,0,2,3,5,10,14,21,36,40,6,17,67,0,0,2,7,11,12,0,0,5,5,5,11,17,3,4,22,30,36,8,0,1,5,7,16,27,0,30,10,14,27,41,69,25,10,11,20,33,5,7,0,1,5,5,5,5,7,11,15,5,14,6,6,7,28,0,5,14,2,2,3,8,10,12,1,1,9,22,3,3,5,15,18,22,37]],"container":"<table class=\"compact\">\n <thead>\n <tr>\n <th>Etnia<\/th>\n <th>Sexo<\/th>\n <th>Ano<\/th>\n <th>Desempenho<\/th>\n <th>Faltas<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"dom":"tBp","buttons":[{"extend":"csv","filename":"Ausencias"},{"extend":"excel","filename":"Ausencias"},{"extend":"pdf","filename":"Ausencias"}],"pageLength":12,"columnDefs":[{"className":"dt-right","targets":4},{"name":"Etnia","targets":0},{"name":"Sexo","targets":1},{"name":"Ano","targets":2},{"name":"Desempenho","targets":3},{"name":"Faltas","targets":4}],"order":[],"autoWidth":false,"orderClasses":false,"lengthMenu":[10,12,25,50,100]}},"evals":[],"jsHooks":[]}</script> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - III <img src="Modelos_dados_de_contagem_files/figure-html/unnamed-chunk-7-1.png" width="47%" style="display: block; margin: auto;" /> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - IV **Modelo** ``` ## ## Call: ## glm(formula = Faltas ~ Etnia + Sexo + Ano + Desempenho, family = poisson(log), ## data = ausencia) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 2.71538 0.06468 41.980 < 2e-16 *** ## EtniaN -0.53360 0.04188 -12.740 < 2e-16 *** ## SexoM 0.16160 0.04253 3.799 0.000145 *** ## AnoF1 -0.33390 0.07009 -4.764 1.90e-06 *** ## AnoF2 0.25783 0.06242 4.131 3.62e-05 *** ## AnoF3 0.42769 0.06769 6.319 2.64e-10 *** ## DesempenhoSL 0.34894 0.05204 6.705 2.02e-11 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 2073.5 on 145 degrees of freedom ## Residual deviance: 1696.7 on 139 degrees of freedom ## AIC: 2299.2 ## ## Number of Fisher Scoring iterations: 5 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - V **Diagnóstico do modelo: Análise gráfica de resíduos.** Código diponível em https://www.ime.usp.br/~giapaula/diag_pois. <img src="Modelos_dados_de_contagem_files/figure-html/unnamed-chunk-9-1.png" width="40%" style="display: block; margin: auto;" /> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - VI **Diagnóstico do modelo: envelope simulado.** Código disponível em https://www.ime.usp.br/~giapaula/envel_pois. <img src="Modelos_dados_de_contagem_files/figure-html/unnamed-chunk-10-1.png" width="45%" style="display: block; margin: auto;" /> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - VII **Observações:** * Embora haja ganho no poder explicativo do modelo com as variáveis explicativas, o ajuste do modelo não foi satisfatório (Desvio muito próximo do desvio nulo). ``` ## Analysis of Deviance Table ## ## Model 1: Faltas ~ 1 ## Model 2: Faltas ~ Etnia + Sexo + Ano + Desempenho ## Resid. Df Resid. Dev Df Deviance Pr(>Chi) ## 1 145 2073.5 ## 2 139 1696.7 6 376.83 < 2.2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` -- * Há fortes indícios de que a ligação usada é inapropriada. -- * O gráfico de envelope mostra que a distribuição assumida não é adequada. Isto pode ser devido à suposição restritiva de que a média e a variância são iguais. Na prática, muitas vezes a variância é maior que a média. Este fenômeno é conhecido como **<span style="color:orange">sobredispersão</span>**. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Modelo Binomal Negativo Assumimos que `\(Y_1,\ldots,Y_n\)` são variáveis aleatórias independentes tais que `\(Y_i\sim \text{BN}(\mu_i,\phi)\)`. A função de probabilidade é dada por $$ f(y_i;\mu_i,\phi)=\frac{\Gamma(\phi+y_i)}{\Gamma(y_i+1)\Gamma(\phi)}\left\(\frac{\mu_i}{\mu_i+\phi}\right\)^{y_i}\left\(\frac{\phi}{\mu_i+\phi}\right\)^{\phi}, \quad y_i=0,1,2,\ldots . $$ Temos que `\(\text{E}(Y_i)=\mu_i\)` e `\(\text{Var}(\mu_i)+\mu_i^2/\phi\)`. Assumimos que a parte sistemática é dada por `\(g(\mu_i)=\eta_i=\mathbb{x}_i^t\mathbb{\beta}\)`, em que `\(\mathbb{x}_i=(x_{i1},\ldots,x_{ip})^t\)` e `\(\mathbb{\beta}=(\beta_1,\ldots,\beta_p)^t\)`. Assim como nos modelos de Poisson, as funções de ligação mais usadas são a logarítmica `\((g(\mu_i)=\log(\mu_i))\)`, raiz quadrada `\((g(\mu_i)=\sqrt\mu_i)\)` e identidade `\((g(\mu_i)=\mu_i)\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - I Voltando aos dados dos estudantes australianos do *Exemplo 2*, assumamos agora que `\(Y\sim \text{BN}(\mu_{ijkl},\phi)\)` com componente sistemática dada por -- $$ \log(\mu_{ijklm})=\alpha+ \beta_1 \text{ano}_i + \beta_2 \text{etnia}_j + \beta_3 \text{desempenho}_k + \beta_4 \text{sexo}_l. $$ --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - I ``` ## ## Call: ## glm.nb(formula = Faltas ~ Etnia + Sexo + Ano + Desempenho, data = ausencia, ## init.theta = 1.274892646, link = log) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 2.89458 0.22842 12.672 < 2e-16 *** ## EtniaN -0.56937 0.15333 -3.713 0.000205 *** ## SexoM 0.08232 0.15992 0.515 0.606710 ## AnoF1 -0.44843 0.23975 -1.870 0.061425 . ## AnoF2 0.08808 0.23619 0.373 0.709211 ## AnoF3 0.35690 0.24832 1.437 0.150651 ## DesempenhoSL 0.29211 0.18647 1.566 0.117236 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for Negative Binomial(1.2749) family taken to be 1) ## ## Null deviance: 195.29 on 145 degrees of freedom ## Residual deviance: 167.95 on 139 degrees of freedom ## AIC: 1109.2 ## ## Number of Fisher Scoring iterations: 1 ## ## ## Theta: 1.275 ## Std. Err.: 0.161 ## ## 2 x log-likelihood: -1093.151 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - II ``` ## ## Call: ## glm.nb(formula = Faltas ~ Etnia + Ano, data = ausencia, init.theta = 1.249142793, ## link = log) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 3.0382 0.1957 15.529 < 2e-16 *** ## EtniaN -0.5611 0.1547 -3.628 0.000286 *** ## AnoF1 -0.3855 0.2274 -1.695 0.090019 . ## AnoF2 0.1846 0.2313 0.798 0.424744 ## AnoF3 0.2550 0.2407 1.060 0.289332 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for Negative Binomial(1.2491) family taken to be 1) ## ## Null deviance: 192.04 on 145 degrees of freedom ## Residual deviance: 167.84 on 141 degrees of freedom ## AIC: 1107.8 ## ## Number of Fisher Scoring iterations: 1 ## ## ## Theta: 1.249 ## Std. Err.: 0.157 ## ## 2 x log-likelihood: -1095.801 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - III A retirada das variáveis não significativas *Sexo* e *Desempenho* não parece ter produzido uma melhora significativa no modelo. Para confirmar isto, realizamos a análise dos desvios: ``` r anova(fit.bn1,fit.bn2,test="Chisq") ``` ``` ## Likelihood ratio tests of Negative Binomial Models ## ## Response: Faltas ## Model theta Resid. df 2 x log-lik. Test ## 1 Etnia + Ano 1.249143 141 -1095.801 ## 2 Etnia + Sexo + Ano + Desempenho 1.274893 139 -1093.151 1 vs 2 ## df LR stat. Pr(Chi) ## 1 ## 2 2 2.650185 0.2657784 ``` >Podemos testar se há interação entre os fatores *Ano* e *Etnia* na tentativa de melhorar o modelo. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - IV ``` ## ## Call: ## glm.nb(formula = Faltas ~ Etnia * Ano, data = ausencia, link = log, ## init.theta = 1.357377778) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 2.6280 0.2495 10.535 < 2e-16 *** ## EtniaN 0.1311 0.3455 0.379 0.70436 ## AnoF1 0.1784 0.3195 0.558 0.57664 ## AnoF2 0.8267 0.3172 2.606 0.00916 ** ## AnoF3 0.3708 0.3337 1.111 0.26650 ## EtniaN:AnoF1 -0.9916 0.4394 -2.257 0.02403 * ## EtniaN:AnoF2 -1.2392 0.4466 -2.775 0.00552 ** ## EtniaN:AnoF3 -0.1763 0.4636 -0.380 0.70380 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for Negative Binomial(1.3574) family taken to be 1) ## ## Null deviance: 205.56 on 145 degrees of freedom ## Residual deviance: 167.84 on 138 degrees of freedom ## AIC: 1102.6 ## ## Number of Fisher Scoring iterations: 1 ## ## ## Theta: 1.357 ## Std. Err.: 0.174 ## ## 2 x log-likelihood: -1084.638 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - V ``` ## Likelihood ratio tests of Negative Binomial Models ## ## Response: Faltas ## Model theta Resid. df 2 x log-lik. Test df LR stat. ## 1 Etnia + Ano 1.249143 141 -1095.801 ## 2 Etnia * Ano 1.357378 138 -1084.638 1 vs 2 3 11.16345 ## Pr(Chi) ## 1 ## 2 0.01087407 ``` >Pelo teste da razão de verossimilhanças a interação foi estatisticamente significativa para o modelo, a um nível de significância de 5%. >O desvio do modelo final foi 167,84 em 138 graus de liberdade. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - VI **Diagnóstico do modelo: Análise gráfica de resíduos.** Código diponível em https://www.ime.usp.br/~giapaula/diag_nbin. <img src="Modelos_dados_de_contagem_files/figure-html/unnamed-chunk-17-1.png" width="40%" style="display: block; margin: auto;" /> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - VII **Diagnóstico do modelo: envelope simulado.** Código disponível em https://www.ime.usp.br/~giapaula/envel_nbin. <img src="Modelos_dados_de_contagem_files/figure-html/unnamed-chunk-18-1.png" width="40%" style="display: block; margin: auto;" /> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Tabelas de contingência: Delineamentos Amostrais Aqui consideramos experimentos com informações coletadas apenas para variáveis binárias. Os valores obtidos desses experimentos são representados na forma de uma tabela de contingência `\(2\times 2\)` da forma | Covariável || Resposta | || Total | |:-------------:||:-------------:|:------------:||:--------:| | || `\(C_1\)` | `\(C_2\)` || | | `\(A_1\)` || `\(n_{11}\)` | `\(n_{12}\)` || `\(n_{1+}\)` | | `\(A_2\)` || `\(n_{21}\)` | `\(n_{22}\)` || `\(n_{2+}\)` | | Total || `\(n_{+1}\)` | `\(n_{+2}\)` || `\(n\)` | * Totais marginais: `\(n_{1+}\)`, `\(n_{2+}\)`, `\(n_{+1}\)`, `\(n_{+2}\)`. * Total geral ou tamanho amostral: `\(n\)`. > A forma como esses totais são considerados (fixos ou aleatórios) define o modelo probabilístico que será considerado. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Produto de Binomiais Independentes - I * Experimento planejado de forma que duas amostras aleatórias `\((A_1 \text{ e } A_2)\)` de tamanhos `\(n_{1+}\)` e `\(n_{2+}\)`, respectivamente, sejam obtidas a partir de duas populações independentes. * Observa-se, para os elementos em `\(A_1\)` e `\(A_2\)`, quantos classificam na categoria `\(C_1\)` e quantos em `\(C_2\)`. * Definem-se as variáveis aleatórias * `\(N_{11}\)` com distribuição `\(Bin(n_{1+},p_{11})\)`, associada à primeira linha da tabela, em que `\(p_{11}\)` representa a probabilidade de qualquer elemento da população `\(A_1\)` ser classificado na categoria `\(C_1\)`. * `\(N_{21}\)` com distribuição `\(Bin(n_{2+},p_{21})\)`, associada à segunda linha da tabela, em que `\(p_{21}\)` representa a probabilidade de qualquer elemento da população `\(A_2\)` ser classificado na categoria `\(C_1\)`. -- * **Observação:** `\(\sum_{j=1}^2 p_{1j}=1\)` e `\(\sum_{j=1}^2 p_{2j}=1\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Produto de Binomiais Independentes - II * A distribuição conjunta associada à tabela é o produto de duas binomiais independentes: `$$P(N_{11}=n_{11},N_{21}=n_{21})={n_{1+}\choose n_{11}} p_{11}^{n_{11}}(1-p_{11})^{n_{1+}-n_{11}}{n_{2+}\choose n_{21}} p_{21}^{n_{21}}(1-p_{21})^{n_{2+}-n_{21}}.$$` -- * Caso a variável resposta apresente mais de duas categorias `\((r>2)\)`, a distribuição associada será o produto de multinomias independentes, sendo uma multinomial associada a cada linha da tabela. -- * Esse esquema amostral é equivalente a um delineamento amostral estratificado. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Produto de Binomiais Independentes - III Na literatura há controvérsias sobre qual verossimilhança usar: a expressão do produto de binomiais independentes, ou a verossimilhança condicional proveniente da distribuição obtida ao se considerar que ambas as marginais da tabela são fixas. Fisher, Cox, Yates, Mantel e outros advogam pela **<span style="color:orange">abordagem condicional</span>**: Para obter a verossimilhança condicional temos de calcular: `$$\textrm{P}(X=n_{11}|n_{1+},n_{2+},n_{+1},n_{+2},n)=\textrm{P}(X=n_{11}|X+Y=n_{+1}),$$` onde `\(X\sim \textrm{Bin}(n_{1+},p_{11})\)` e `\(Y\sim \textrm{Bin}(n_{2+},p_{21})\)`. Então <div class="math"> \[\begin{align*} \textrm{P}&(X=n_{11}|X+Y=n_{+1})=\dfrac{\textrm{P}(X=n_{11},X+Y=n_{+1})}{\textrm{P}(X+Y=n_{+1})}\\ &=\dfrac{\textrm{P}(X=n_{11},Y=n_{+1}-n_{11})}{\textrm{P}(X+Y=n_{+1})}\\ &=\dfrac{\textrm{P}(X=n_{11})\textrm{P}(Y=n_{+1}-n_{11})}{\sum_{i=1}^{n_{+1}}\textrm{P}(X=i,Y=n_{+1}-i)}\\ &=\dfrac{\binom{n_{1+}}{n_{11}}p_{11}^{n_{11}}(1-p_{11})^{n_{1+}-n_{11}}\binom{n_{2+}}{n_{+1}-n_{11}}p_{21}^{n_{+1}-n_{11}}(1-p_{21})^{n_{2+}-n_{+1}+n_{11}}}{\sum_i\binom{n_{1+}}{i}p_{11}^i(1-p_{11})^{n_{1+}-i}\binom{n_{2+}}{n_{+1}-i}p_{21}^{n_{+1}-i}(1-p_{21})^{n_{2+}-n_{+1}+i}}, \end{align*}\] </div> onde o somatório varia de `\(\textrm{Máx}(0,n_{+1}-n_{2+})\)` a `\(\textrm{Min}(n_{+1},n_{1+})\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Produto de Binomiais Independentes - IV Fazendo `\(q_{11}=1-p_{11}\)` e `\(q_{21}=1-p_{21}\)`, temos <div class="math"> \[\begin{align*} \textrm{P}(X=n_{11}|X+Y=n_{+1})&=\dfrac{\binom{n_{1+}}{n_{11}}p_{11}^{n_{11}}q_{11}^{-n_{11}}\binom{n_{2+}}{n_{+1}-n_{11}}p_{21}^{-n_{11}}q_{21}^{n_{11}}}{\sum_i\binom{n_{1+}}{i}p_{11}^i q_{11}^{-i}\binom{n_{2+}}{n_{+1}-i}p_{21}^{-i}q_{21}^i}\\ &=\dfrac{\binom{n_{1+}}{n_{11}}\binom{n_{2+}}{n_{+1}-n_{11}}\psi^{n_{11}}}{\sum_i\binom{n_{1+}}{i}p_{11}^i \binom{n_{2+}}{n_{+1}-i}\psi^i}, \end{align*}\] </div> onde `\(\psi=\dfrac{p_{11}q_{21}}{p_{21}q_{11}}\)`. Sob `\(H_0\)`, `\(\psi=1\)` e `$$\textrm{P}(X=n_{11}|X+Y=n_{+1})=\dfrac{\binom{n_{1+}}{n_{11}}\binom{n_{2+}}{n_{+1}-n_{11}}}{\binom{n}{n_{+1}}},$$` que é uma **<span style="color:orange">Distribuição Hipergeométrica</span>**. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Multinomial * Experimento delineado de forma que se obtenha uma amostra aleatória de tamanho `\(n\)` fixo de uma população. * Observa-se quantos elementos dessa amostra classificam-se em cada uma das possíveis combinações da tabela: `\((A_1,C_1)\)`, `\((A_1,C_2)\)`, `\((A_2,C_1)\)` ou `\((A_2,C_2)\)`. -- * O modelo associado é dado por: `$$P(N_{11}=n_{11}, N_{12}=n_{12}, N_{21}=n_{21}, N_{22}=n_{22})=n!\prod_{i,j=1}^2 \frac{p_{ij}^{n_{ij}}}{n_{ij}!},$$` onde `\(n_{ij}\geq 0\)`, `\(\sum_{i=1}^2 n_{ij}= n\)` e `\(p_{ij}=1\)`. -- * Este delineamento amostral é equivalente ao processo de amostragem aleatória simples, em que se seleciona uma amostra aleatória de tamanho `\(n\)` a partir de uma população suficientemente grande. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Produto de Poisson * Experimento delineado de forma que apenas a duração do experimento é estabelecida. * Os tamanhos amostrais são aleatórios, dado que serão conhecidos apenas ao final do experimento. * Assume-se `\(N_{ij}\)`, `\((i,j=1,2)\)`, independentes com distribuição de Poisson de média `\(\mu_{ij}=t\lambda_{ij}\)`. A taxa média por unidade de tempo é representada por `\(\lambda_{ij}\)` e a duração do tempo representada por `\(t\)`. -- * O modelo probabilístico associado é `$$P(N_{11}=n_{11}, N_{12}=n_{12}, N_{21}=n_{21}, N_{22}=n_{22})=\prod_{i,j=1}^2 \frac{e^{-\mu_{ij}}\mu_{ij}^{n_{ij}}}{n_{ij}!}\quad \mu_{ij}>0.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Hipergeométrico * Considerado em situações em que os indivíduos são inicialmente alocados aleatoriamente em dois grupos considerados equivalentes, observando-se para cada um deles uma resposta binária. * Assumindo que os dois grupos não diferem significativamente, `\(n_{+1}\)` é considerado fixo. * Dessa forma, `\(n_{1+}\)`, `\(n_{2+}\)` e `\(n_{+2}\)` são fixos. Dado o condicionamento nas marginais fixas, o conhecimento de qualquer uma das frequências `\(n_{ij}\)` determina as restantes. -- * Sob essa suposição, o modelo probabilístico associado é dado por `$$P(N_{11}=n_{11})=\frac{{n_{1+}\choose n_{11}}{n_{2+}\choose n_{+1}-n_{11}}}{{n\choose n_{+1}}}$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Observações * Os esquemas amostrais anteriores são os mais comuns. * Nem todas as tabelas de contingência são geradas sob esses esquemas amostrais. * Por exemplo, há casos em que as amostras não são independentes. Ou casos em que sequer é considerada uma amostra aleatória. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Modelos Log-lineares * Modelam os padrões de associação e de interação entre variáveis categóricas. * São adequados para esquemas amostrais de Poisson, Multinomiais e Produto de Multinomiais. * Apropriados para situações em que não há distinção clara entre resposta e variáveis explicativas ou quando há mais de uma resposta. ><span style="color:orange">**Essa é a principal diferença entre os modelos logísticos e os loglineares**</span>. * Assumem que as variáveis discretas são nominais, mas é possível fazer ajustes para considerar ordinalidade e com dados pareados. * São mais gerais que os modelos logit, mas alguns modelos loglineares têm correspondência direta com os modelos logit. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelos Loglineares para tabelas de duas vias - I Sejam `\(\mu_{ij}\)` as contagens esperadas, `\(E(n_{ij})\)`, em uma tabela `\(I\times J\)` formada pelo cruzamento das variáveis `\(A\)` e `\(B\)`. O objetivo é modelar as contagens `\(\mu_{ij}=n\pi_{ij}\)`. Sejam `\(N=I\times J\)`, as contagens das caselas, observações independentes de uma variável aleatória Poisson, `\(n_{ij}=\text{Po}(\mu_{ij})\)`. Um modelo loglinear análogo à ANOVA de duas vias com interação é dado por `$$\log(\mu_{ij})=\lambda+\lambda_i^A+\lambda_j^B+\lambda_{ij}^{AB},$$` onde `\(i=1,\ldots,I; j=1,\ldots,J\)` são os níveis das variáveis categóricas `\(A\)` e `\(B\)`. Com o objetivo de garantir a unicidade das estimativas, impõem-se as restrições `\(\sum_i \lambda_i=\sum_j \lambda_j=\sum_i\sum_j \lambda_{ij}=0\)`. >Esse modelo é superparametrizado porque o termo `\(\lambda_{ij}\)` já tem `\(I\times J\)` parâmetros correspondentes às médias das caselas `\(\mu_{ij}\)`. A constante `\(\lambda\)` e os efeitos principais `\(\lambda_i\)` e `\(\lambda_j\)` correspondem a `\(1+I+J\)` parâmetros adicionais. -- Nesse contexto, há dois tipos de modelos a serem considerados * Modelo de independência (A, B). * Modelo saturado (AB). --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Modelo de Independência - I A independência pode ser escrita em termos de probabilidades das caselas como o produto das probabilidades marginais `$$\pi_{ij}=\pi_{i+}\pi_{+j}\quad i=1,\ldots,I, j=1,\ldots,J$$` e em termos das frequências das caselas `$$\mu_{ij}=n\pi_{ij}=n\pi_{i+}\pi_{+j}\quad i=1,\ldots,I, j=1,\ldots,J.$$` Tomando logaritmos obtemos o modelo loglinear de independência `$$\log(\mu_{ij})=\lambda + \lambda_i^A + \lambda_j^B,$$` onde `\(A\)` e `\(B\)` denotam as duas variáveis categóricas. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo de Independência - II **Observações:** * `\(\lambda\)` representa o efeito global, ou média geral dos logaritmos das contagens esperadas, garante que `\(\sum_i\sum_j \mu_{ij}=n\)`, isto é, as contagens esperadas sob o modelo ajustado totalizam o tamanho amostral `\(n\)`. * `\(\lambda_i^A\)` representa o efeito principal da variável A. Garante que `\(\sum_j \mu_{ij}=n_{i+}\)`, isto é, os totais marginais sob o modelo ajustado totalizam as contagens marginais observadas. Representa o efeito da classificação na linha `\(i\)`. * `\(\lambda_j^B\)` representa o efeito principal da variável B. Garante que `\(\sum_i \mu_{ij}=n_{+j}\)`. Representa o efeito da classificação na coluna `\(j\)`. * Restrições: `\(\lambda_I^A=\lambda_J^B=0\)`, ou alternativamente, `\(\sum_i \lambda_i^A=\sum_j \lambda_j^B=0\)`. -- Os valores estimados por Máxima Verossimilhança para as contagens são os mesmos valores esperados sob o teste de independência em tabelas de duas vias, i.e., `\(E(\mu_{ij})=n_{i+}n_{+j}/n\)`. Portanto, as estatísticas `\(\chi^2\)` e `\(G^2\)` (deviance) para independência são testes de bondade do ajuste para o modelo loglinear de independência, onde as hipóteses são: o modelo de independência é o correto vs. o modelo saturado é o correto. Esse modelo também implica que todos os odds ratios são iguais a 1. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo de Independência - III Dependendo do tipo de restrição imposta, podemos obter diferentes estimativas dos parâmetros (i.e, diferentes valores dos `\(\lambda s\)`). Entretanto, as diferenças (log-odds) são únicos: `\(\lambda_i^A - \lambda_{i'}^A\)` e `\(\lambda_j^B - \lambda_{j'}^B\)`, onde o subíndice `\(i\)` denota um nível da variável categórica `\(A\)` e `\(i'\)` denota outro nível da mesma varíável; similarmente para `\(B\)`. Dessa forma, os odds ratios também são únicos. `$$\begin{align*} \log(odds)&=\log\left(\frac{\mu_{i1}}{\mu_{i2}}\right)=\log(\mu_{i1})-\log(\mu_{i2})\\ &=(\lambda+\lambda_i^A+\lambda_1^B)-(\lambda+\lambda_i^A+\lambda_2^B)\\ &=(\lambda_1^B-\lambda_2^B)\\ \end{align*}$$` `$$\Longrightarrow odds=\exp(\lambda_1^B-\lambda_2^B)$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo de Independência - IV Sob o modelo de independência, espera-se que o log-odds seja 0, isto é, que o OR=1. Assim, `$$\begin{align*} \log(odds)&=\log\left(\frac{\mu_{11}\mu_{22}}{\mu_{12}\mu_{21}}\right)\\ &=\log(\mu_{11})+\log(\mu_{22})-\log(\mu_{12})-\log(\mu_{21})\\ &=\lambda+\lambda_1^A+\lambda_1^B+\lambda+\lambda_2^A+\lambda_2^B-\lambda-\lambda_1^A-\lambda_2^B-\lambda-\lambda_2^A-\lambda_1^B\\ &=0 \end{align*}$$` -- >**Observação:** O odds ratio mede a força da associação e depende apenas do termo de interação `\({\lambda_{ij}^{AB}}\)`, que claramente não aparece neste modelo. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Saturado No modelo saturado as `\(N=I\times J\)` contagens nas caselas ainda são assumidas como observações independentes de uma variável aleatória Poisson, mas nãos e assume independência entre `\(A\)` e `\(B\)`. O modelo é representado como `$$\log(\mu_{ij})=\lambda+\lambda_i^A+\lambda_j^B+\lambda_{ij}^{AB}.$$` Nesse modelo também se aplicam as restrições semelhantes às da ANOVA para contornar a superparametrização. **Estimativas dos parâmetros:** O odds ratio está diretamente relacionado com o termo de interação. Por exemplo, para uma tabela `\(2\times 2\)`: `$$\begin{align*} \log(OR)&=\log\left(\frac{\mu_{11}\mu_{22}}{\mu_{12}\mu_{21}}\right)\\ &=\lambda+\lambda_1^A+\lambda_1^B+\lambda_{11}^{AB}+\lambda+\lambda_2^A+\lambda_2^B+\lambda_{22}^{AB}-\lambda-\lambda_1^A-\lambda_2^B-\lambda_{12}^{AB}-\lambda-\lambda_2^A-\lambda_1^B-\lambda_{21}^{AB}\\ &=\lambda_{11}^{AB}+\lambda_{22}^{AB}-\lambda_{12}^{AB}-\lambda_{21}^{AB}. \end{align*}$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo - I Um estudo duplo-cego investigou o efeito terapeutico da vitamina C para tratar resfriados comuns. O estudo foi conduzido durante um período de duas semanas em uma amostra de 280 esquiadores franceses. Uma observação foi perdida por inconsistências. | Tratamento || Resfriado | Sem resfriado || Total | |:----------:||----------:|:-------------:||:-----:| | || | || | | Placebo || 31 | 109 || 140 | | Vit. C || 17 | 122 || 139 | | Total || 48 | 231 || 279 | ``` r dados<-matrix(c(31, 17, 109, 122), ncol=2, dimnames=list(Tratamento=c("Placebo", "VitC"), Resfriado=c("Com", "Sem"))) dados ``` ``` ## Resfriado ## Tratamento Com Sem ## Placebo 31 109 ## VitC 17 122 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - II ``` r library(vcd) assocstats(dados) ``` ``` ## X^2 df P(> X^2) ## Likelihood Ratio 4.8717 1 0.027301 ## Pearson 4.8114 1 0.028272 ## ## Phi-Coefficient : 0.131 ## Contingency Coeff.: 0.13 ## Cramer's V : 0.131 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - III ``` r OR<-oddsratio(dados, log=FALSE) OR ``` ``` ## odds ratios for Tratamento and Resfriado ## ## [1] 2.041015 ``` ``` r confint(OR) ``` ``` ## 2.5 % 97.5 % ## Placebo:VitC/Com:Sem 1.070353 3.89193 ``` >A chance de pegar resfriado é aproximadamente o dobro para aqueles que tomaram o placebo, quando comparados com os que tomaram vitamina C. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - IV ### Modelo de Independência ``` r dados<-as.table(dados) mod.indep<-loglin(dados, list(1,2), fit=TRUE, param=TRUE) ``` ``` ## 2 iterations: deviation 0 ``` ``` r mod.indep[-3]#eliminando a impressão do terceiro item da lista (graus de liberdade) ``` ``` ## $lrt ## [1] 4.871697 ## ## $pearson ## [1] 4.811413 ## ## $margin ## $margin[[1]] ## [1] "Tratamento" ## ## $margin[[2]] ## [1] "Resfriado" ## ## ## $fit ## Resfriado ## Tratamento Com Sem ## Placebo 24.08602 115.91398 ## VitC 23.91398 115.08602 ## ## $param ## $param$`(Intercept)` ## [1] 3.963656 ## ## $param$Tratamento ## Placebo VitC ## 0.003584245 -0.003584245 ## ## $param$Resfriado ## Com Sem ## -0.7856083 0.7856083 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - V **log-odds estimados:** `$$\begin{align*} \log(\mu_{11})&=3.963656+0.003584245-0.7856083\\ \cdots\\ \log(\mu_{22})&=3.963656-0.003584245+0.7856083\\ \end{align*}$$` **odds de resfriado:** `$$\begin{align*} \exp(\lambda_1^R-\lambda_2^R)&=\exp(-0.7856083-0.7856083)\\ &=\exp(-1.5712)\\ &=0.208 \end{align*}$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - VI ### Modelo Saturado ``` r mod.sat<-loglin(dados, list(c(1,2)), fit=TRUE, param=TRUE) ``` ``` ## 2 iterations: deviation 0 ``` ``` r mod.sat ``` ``` ## $lrt ## [1] 0 ## ## $pearson ## [1] 0 ## ## $df ## [1] 0 ## ## $margin ## $margin[[1]] ## [1] "Tratamento" "Resfriado" ## ## ## $fit ## Resfriado ## Tratamento Com Sem ## Placebo 31 109 ## VitC 17 122 ## ## $param ## $param$`(Intercept)` ## [1] 3.940642 ## ## $param$Tratamento ## Placebo VitC ## 0.1220252 -0.1220252 ## ## $param$Resfriado ## Com Sem ## -0.8070421 0.8070421 ## ## $param$Tratamento.Resfriado ## Resfriado ## Tratamento Com Sem ## Placebo 0.1783618 -0.1783618 ## VitC -0.1783618 0.1783618 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - VII **Contagens esperadas sob o modelo ajustado** ``` ## Resfriado ## Tratamento Com Sem ## Placebo 31 109 ## VitC 17 122 ``` **Parâmetros estimados (efeitos principais)** ``` ## $`(Intercept)` ## [1] 3.940642 ## ## $Tratamento ## Placebo VitC ## 0.1220252 -0.1220252 ## ## $Resfriado ## Com Sem ## -0.8070421 0.8070421 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - VIII **Parâmetros estimados (Interação)** ``` ## $Tratamento.Resfriado ## Resfriado ## Tratamento Com Sem ## Placebo 0.1783618 -0.1783618 ## VitC -0.1783618 0.1783618 ``` >Temos então: `$$\begin{align*} \log(\mu_{11})&=3.940642+0.1220252-0.8070421+0.1783618\\ \cdots\\ \log(\mu_{22})&=3.940642-0.1220252+0.8070421+0.1783618\\ \end{align*}$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - IX ### Alternativa: ``` r dados1<-as.data.frame(dados) dados1 ``` ``` ## Tratamento Resfriado Freq ## 1 Placebo Com 31 ## 2 VitC Com 17 ## 3 Placebo Sem 109 ## 4 VitC Sem 122 ``` ``` r modelo.ind<-glm(dados1$Freq~dados1$Tratamento+dados1$Resfriado, family=poisson()) modelo.ind[1] ``` ``` ## $coefficients ## (Intercept) dados1$TratamentoVitC dados1$ResfriadoSem ## 3.181631652 -0.007168489 1.571216700 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - X ``` r summary(modelo.ind) ``` ``` ## ## Call: ## glm(formula = dados1$Freq ~ dados1$Tratamento + dados1$Resfriado, ## family = poisson()) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 3.181632 0.156179 20.372 <2e-16 *** ## dados1$TratamentoVitC -0.007168 0.119738 -0.060 0.952 ## dados1$ResfriadoSem 1.571217 0.158626 9.905 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 135.4675 on 3 degrees of freedom ## Residual deviance: 4.8717 on 1 degrees of freedom ## AIC: 34.004 ## ## Number of Fisher Scoring iterations: 4 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - XI ### Modelo Saturado ``` r mod.sat<-loglin(dados, list(c(1,2)), fit=TRUE, param=TRUE) ``` ``` ## 2 iterations: deviation 0 ``` ``` r mod.sat ``` ``` ## $lrt ## [1] 0 ## ## $pearson ## [1] 0 ## ## $df ## [1] 0 ## ## $margin ## $margin[[1]] ## [1] "Tratamento" "Resfriado" ## ## ## $fit ## Resfriado ## Tratamento Com Sem ## Placebo 31 109 ## VitC 17 122 ## ## $param ## $param$`(Intercept)` ## [1] 3.940642 ## ## $param$Tratamento ## Placebo VitC ## 0.1220252 -0.1220252 ## ## $param$Resfriado ## Com Sem ## -0.8070421 0.8070421 ## ## $param$Tratamento.Resfriado ## Resfriado ## Tratamento Com Sem ## Placebo 0.1783618 -0.1783618 ## VitC -0.1783618 0.1783618 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - XII **Contagens esperadas sob o modelo ajustado** ``` ## Resfriado ## Tratamento Com Sem ## Placebo 31 109 ## VitC 17 122 ``` **Parâmetros estimados (efeitos principais)** ``` ## $`(Intercept)` ## [1] 3.940642 ## ## $Tratamento ## Placebo VitC ## 0.1220252 -0.1220252 ## ## $Resfriado ## Com Sem ## -0.8070421 0.8070421 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - XIII ``` r anova(modelo.ind) ``` ``` ## Analysis of Deviance Table ## ## Model: poisson, link: log ## ## Response: dados1$Freq ## ## Terms added sequentially (first to last) ## ## ## Df Deviance Resid. Df Resid. Dev ## NULL 3 135.468 ## dados1$Tratamento 1 0.004 2 135.464 ## dados1$Resfriado 1 130.592 1 4.872 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo - XIV ``` r modelo.sat<-glm(dados1$Freq~dados1$Tratamento*dados1$Resfriado, family=poisson()) modelo.sat ``` ``` ## ## Call: glm(formula = dados1$Freq ~ dados1$Tratamento * dados1$Resfriado, ## family = poisson()) ## ## Coefficients: ## (Intercept) ## 3.4340 ## dados1$TratamentoVitC ## -0.6008 ## dados1$ResfriadoSem ## 1.2574 ## dados1$TratamentoVitC:dados1$ResfriadoSem ## 0.7134 ## ## Degrees of Freedom: 3 Total (i.e. Null); 0 Residual ## Null Deviance: 135.5 ## Residual Deviance: -1.021e-14 AIC: 31.13 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Modelos Loglineares para tabelas de três vias Seja `\(\mu_{ijk}\)` que representa a média dos níveis `\(i\)`, `\(j\)` e `\(k\)` das variáveis `\(A\)`, `\(B\)` e `\(C\)` respectivamente. O modelo é dado por `$$\log(\mu_{ijk})=\lambda+\lambda_i^A+\lambda_j^B+\lambda_k^C+\lambda_{ij}^{AB}+\lambda_{ik}^{AC}+\lambda_{jk}^{BC}+\lambda_{ijk}^{ABC}.$$` -- Há diversos modelos que podem ser testados e ajustados, são eles: * Saturado * Independência completa * Independência conjunta * Independência condicional * Associação Homogênea --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Saturado - I É o modelo default e serve para testar a bondade do ajuste dos outros modelos. Sejam `\(N=I\times J\times K\)` as contagens nas caselas, assumidas como observações independentes de uma variável aleatória Poisson. O modelo é dado por `$$\log(\mu_{ijk})=\lambda+\lambda_i^A+\lambda_j^B+\lambda_k^C+\lambda_{ij}^{AB}+\lambda_{ik}^{AC}+\lambda_{jk}^{BC}+\lambda_{ijk}^{ABC}.$$` Sujeito a `\(\lambda_I^A=\lambda_J^B=\lambda_K^C=\lambda_{Ij}^{AB}=\cdots=\lambda_{ijK}^{ABC}=0.\)` * `\(\lambda\)` representa o efeito global o grande média (na escala logarítmica) das contagens esperadas, garante que `\(\sum_i\sum_j\sum_k \mu_{ijk}=n.\)` * `\(\lambda_i^{A}\)`, `\(\lambda_j^{B}\)` e `\(\lambda_k^{C}\)` representam os efeitos principais das variáveis `\(A\)`, `\(B\)` e `\(C\)`, ou desvios da média global. * `\(\lambda_{ij}^{AB}\)`, `\(\lambda_{ik}^{AC}\)` e `\(\lambda_{jk}^{BC}\)` representam a interação/associação entre duas variáveis quando controlada a terceira (i.e., odds ratios condicionais, testes de associação parcial) e reflete o afastamento da independência. * `\(\lambda_{ijk}^{ABC}\)` representa a associação/interação entre três variáveis e reflete o afastamento da independência. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Saturado - II **Observações:** * Se há termo de interação significante, normalmente não olhamos para os termos de ordem inferior e apenas interpretamos os termos de ordem superior. * O modelo saturado tem ajuste perfeito, com `\(G^2=0\)` e 0 graus de liberdade, pois o número de caselas é igual ao número de parâmetros únicos no modelo. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - I Lembremos das informações sobre admissões nas pós-graduações nos seis maiores departamentos da U.C. Berkeley em 1973. | Departamento || Homens admitidos || Homens rejeitados || Mulheres admitidas || Mulheres rejeitadas | |:------------:||:---------------:||:------------------:||:--------------------:||:---------------------:| | A || 512 || 313 || 89 || 19 | | B || 353 || 207 || 17 || 8 | | C || 120 || 205 || 202 || 391 | | D || 139 || 279 || 131 || 244 | | E || 53 || 138 || 94 || 299 | | F || 22 || 351 || 24 || 317 | Denotemos por `\(D=\)` Departamento, `\(S=\)` Sexo e `\(A=\)` Admissão (admitido ou rejeitado). Seja `\(Y\)` a frequência observada em uma casela particular da tabela de três vias. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - II ``` r UCBAdmissions ``` ``` ## , , Dept = A ## ## Gender ## Admit Male Female ## Admitted 512 89 ## Rejected 313 19 ## ## , , Dept = B ## ## Gender ## Admit Male Female ## Admitted 353 17 ## Rejected 207 8 ## ## , , Dept = C ## ## Gender ## Admit Male Female ## Admitted 120 202 ## Rejected 205 391 ## ## , , Dept = D ## ## Gender ## Admit Male Female ## Admitted 138 131 ## Rejected 279 244 ## ## , , Dept = E ## ## Gender ## Admit Male Female ## Admitted 53 94 ## Rejected 138 299 ## ## , , Dept = F ## ## Gender ## Admit Male Female ## Admitted 22 24 ## Rejected 351 317 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - III ``` r berkeley<-as.data.frame(UCBAdmissions) berkeley$Gender = relevel(berkeley$Gender, ref='Female') berkeley$Dept = relevel(berkeley$Dept, ref='F') berkeley.sat = glm(Freq~Admit*Gender*Dept, family=poisson(), data=berkeley) summary(berkeley.sat) ``` ``` ## ## Call: ## glm(formula = Freq ~ Admit * Gender * Dept, family = poisson(), ## data = berkeley) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 3.17805 0.20412 15.569 < 2e-16 *** ## AdmitRejected 2.58085 0.21171 12.190 < 2e-16 *** ## GenderMale -0.08701 0.29516 -0.295 0.7682 ## DeptA 1.31058 0.23001 5.698 1.21e-08 *** ## DeptB -0.34484 0.31700 -1.088 0.2767 ## DeptC 2.13021 0.21591 9.866 < 2e-16 *** ## DeptD 1.69714 0.22204 7.644 2.11e-14 *** ## DeptE 1.36524 0.22870 5.969 2.38e-09 *** ## AdmitRejected:GenderMale 0.18890 0.30516 0.619 0.5359 ## AdmitRejected:DeptA -4.12505 0.32968 -12.512 < 2e-16 *** ## AdmitRejected:DeptB -3.33462 0.47817 -6.974 3.09e-12 *** ## AdmitRejected:DeptC -1.92041 0.22876 -8.395 < 2e-16 *** ## AdmitRejected:DeptD -1.95888 0.23781 -8.237 < 2e-16 *** ## AdmitRejected:DeptE -1.42370 0.24250 -5.871 4.33e-09 *** ## GenderMale:DeptA 1.83670 0.31672 5.799 6.66e-09 *** ## GenderMale:DeptB 3.12027 0.38572 8.090 5.99e-16 *** ## GenderMale:DeptC -0.43376 0.31687 -1.369 0.1710 ## GenderMale:DeptD 0.13907 0.31938 0.435 0.6632 ## GenderMale:DeptE -0.48599 0.34151 -1.423 0.1547 ## AdmitRejected:GenderMale:DeptA 0.86318 0.40267 2.144 0.0321 * ## AdmitRejected:GenderMale:DeptB 0.03113 0.53349 0.058 0.9535 ## AdmitRejected:GenderMale:DeptC -0.31382 0.33741 -0.930 0.3523 ## AdmitRejected:GenderMale:DeptD -0.10691 0.34013 -0.314 0.7533 ## AdmitRejected:GenderMale:DeptE -0.38908 0.36500 -1.066 0.2864 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 2.6501e+03 on 23 degrees of freedom ## Residual deviance: 1.2612e-13 on 0 degrees of freedom ## AIC: 207.06 ## ## Number of Fisher Scoring iterations: 3 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - IV **Observações:** * Os efeitos principais de `\(D\)`, `\(S\)` e `\(A\)` são difíceis de interpretar e não muito úteis, devido a que há associações de duas e três vias significantes * Por exemplo, o coeficiente estimado da associação SA é `\(0.1889\)`. Dessa forma, o odds ratio estimado para o Departamento `\(F\)` (departamento de referência) é `\(\exp(0.1889)=1.208\)`. A categoria de referência para S é "Mulher" e para A é "admitido". A tabela S `\(\times\)` A no Departamento F, i.e., a tabela parcial é dada por | Dep. F || Rejeitado | Admitido | |:--------:||:----------:|:----------:| | Homem || 351 | 22 | | Mulher || 317 | 24 | De onde o odds ratio estimado é `$$OR=\frac{351\times 24}{317\times 22}=1.208.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - V * A estatística de Wald para o coeficiente SA é dada por `$$z=0.1889/0.3052,$$` que corresponde à estatística `\(\chi^2\)` `\((0.62^2=0.38)\)` com valor `\(p=0.5359\)` e indica que o odds ratio SA para o departamento F não é significativamente diferente de 1 (log-odds ratio não é significativamente diferente de 0). * Para obter o odds ratio SA para qualquer outro departamento, devemos combinar os coeficientes SA com um dos coeficientes DSA. Por exemplo, os odds ratio SA para o departamento A é `$$\exp(0.1889+0.8632)=2.864.$$` * A estatística `\(z\)` indica que o odds ratio SA para o departamento A é significativamente diferente do odds ratio SA no departamento F. * O modelo saturado se torna mais complicado na medida em que o número de variáves aumenta. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo de Independência Completa * É o modelo mais restritivo. Todas as variáveis são assumidas conjuntamente independentes, independentemente de qualquer condicionamento. * Assume-se que as `\(N=I\times J\times K\)` contagens nas caselas são observações independentes de uma variável aleatória Poisson. * Assume-se também que não há interações parciais `\(\lambda_{ij}^{AB}=\lambda_{ik}^{AC}=\lambda_{jk}^{BC}\)`, para todo `\(i,j,k\)`; nem interações triplas `\(\lambda_{ijk}^{ABC}=0\)`, para todo `\(i,j,k.\)` * Restrições devem ser impostas para evitar a superparametrização. * A estrutura do modelo é dada por `$$\log(\mu_{ijk})=\lambda + \lambda_i^A + \lambda_j^B + \lambda_k^C.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - VI ``` r berkeley.ind = glm(Freq~Admit + Gender + Dept, family=poisson(), data=berkeley) summary(berkeley.ind) ``` ``` ## ## Call: ## glm(formula = Freq ~ Admit + Gender + Dept, family = poisson(), ## data = berkeley) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 4.72072 0.04553 103.673 < 2e-16 *** ## AdmitRejected 0.45674 0.03051 14.972 < 2e-16 *** ## GenderMale 0.38287 0.03027 12.647 < 2e-16 *** ## DeptA 0.26752 0.04972 5.380 7.44e-08 *** ## DeptB -0.19927 0.05577 -3.573 0.000352 *** ## DeptC 0.25131 0.04990 5.036 4.74e-07 *** ## DeptD 0.10368 0.05161 2.009 0.044533 * ## DeptE -0.20098 0.05579 -3.602 0.000315 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 2650.1 on 23 degrees of freedom ## Residual deviance: 2097.7 on 16 degrees of freedom ## AIC: 2272.7 ## ## Number of Fisher Scoring iterations: 5 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - VII ``` r berkeley.ind$dev ``` ``` ## [1] 2097.671 ``` O Deviance do modelo, com 16 graus de liberdade indica que o modelo não ajusta bem. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Independência Conjunta * Com três variáveis, há três formas de obter independência conjunta. O pressuposto é que uma variável é independente das outras duas, mas essas duas podem ter qualquer associação arbitrária. * Se A é cojuntamente independente de B e C, denotado por (A, BC), então A é independente de B e de C, e podemos fatorar a distribuição conjunta como o produto da distribuiçao marginal de A e da conjunta de (B, C). * Sem perda de generalidade, para o caso (A, BC), assume-se que as `\(N=I\times J\times K\)` contagens são observações de uma variável aleatória Poisson. * Assume-se também que não há interações envolvendo A: `\(\lambda_{ij}^{AB}=\lambda_{ik}^{AC}=0\)`, para todo `\(i, j, k\)`; nem ha interação tripla `\(\lambda_{ijk}^{ABC}\)`, para todo `\(i,j,k.\)` * A estrutura do modelo é dada por `$$\log(\mu_{ijk})=\lambda + \lambda_i^A + \lambda_j^B + \lambda_k^C + \lambda_{jk}^{BC}.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - VIII Consideremos o modelo onde a admissão é conjuntamente independente de departamento e sexo, denotado por (A, DS): ``` r berkeley.conj = glm(Freq~Admit + Gender + Dept + Gender*Dept, family=poisson(), data=berkeley) berkeley.conj$coef ``` ``` ## (Intercept) AdmitRejected GenderMale DeptA ## 4.88451279 0.45673941 0.08969594 -1.14975125 ## DeptB DeptC DeptD DeptE ## -2.61300665 0.55331192 0.09504355 0.14192713 ## GenderMale:DeptA GenderMale:DeptB GenderMale:DeptC GenderMale:DeptD ## 1.94355622 3.01936502 -0.69106516 0.01646425 ## GenderMale:DeptE ## -0.81123213 ``` ``` r berkeley.conj$dev ``` ``` ## [1] 877.0564 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - IX * O modelo implica que a associação entre D e S não depende do nível da variável A. Dessa forma, a associação entre o departamento e sexo é independente da decisão de aceitar ou rejeitar. * O primeiro coeficiente estimado para a associação DS indica que o odds ratio estimado entre sexo e departamento (especificamente A vs F) é `\(\exp(1.9436)=6,98\)` com IC de 95% dado por `$$\left[\exp(1.695);\exp(2.192)\right]=(5,45; 8,95).$$` * A estatística de bondade do ajuste (Deviance=877.06) indica que o modelo não ajusta bem, pois o valor "Deviance/gl" é muito maior que 1. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo de independência Condicional * Para três variáveis, há três formas de obter independência condicional. Assume-se que duas variáveis são independentes, dada a terceira variável. Por exemplo, se A e B são condicionalmente independentes, dado C, denotado por (AC, BC), então a distribuição de (AB), dado C, pode ser fatorada como o produto das duas marginais, dado C. * Sem perda de generalidade, para o caso (AC, BC), assume-se que as `\(N=I\times J\times K\)` contagens são observações de uma variável aleatória Poisson. * Assume-se também que não há interação parcial entre A e B: `\(\lambda_{ij}^{AB}=0\)`, para todo `\(i, j\)`; nem ha interação de três vias `\(\lambda_{ijk}^{ABC}\)`, para todo `\(i,j,k.\)` * A estrutura do modelo é dada por `$$\log(\mu_{ijk})=\lambda + \lambda_i^A + \lambda_j^B + \lambda_k^C + \lambda_{ik}^{AC} + \lambda_{jk}^{BC}.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - X Consideremos o modelo onde a admissão e departamento são conjuntamente independentes, dado o sexo, denotado por (AS, DS): ``` r berkeley.cond = glm(Freq~Admit + Gender + Dept + Admit*Gender + Dept*Gender, family=poisson(), data=berkeley) berkeley.cond$coef ``` ``` ## (Intercept) AdmitRejected GenderMale ## 4.63964796 0.83048640 0.47267109 ## DeptA DeptB DeptC ## -1.14975125 -2.61300665 0.55331192 ## DeptD DeptE AdmitRejected:GenderMale ## 0.09504355 0.14192713 -0.61035238 ## GenderMale:DeptA GenderMale:DeptB GenderMale:DeptC ## 1.94355622 3.01936502 -0.69106516 ## GenderMale:DeptD GenderMale:DeptE ## 0.01646425 -0.81123213 ``` ``` r berkeley.cond$dev ``` ``` ## [1] 783.607 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Associação Homogênea * O modelo de associações homogêneas também é conhecido como **Modelo sem interações triplas**. Denotado por (AB, AC, BC), a única restrição que esse modelo impõe é que nenhuma associação condicional por pares depende do valor da terceira variável. Por exemplo, o odds ratio condicional entre A e B, dado que C está no primeiro nível, deve ser o mesmo odds ratio condicional entre A e B, dado C no segundo nível, e assim por diante. * Assume-se que as `\(N=I\times J\times K\)` contagens são observações de uma variável aleatória Poisson. * Não há interação tripla `\(\lambda_{ijk}^{ABC}\)`, para todo `\(i,j,k.\)` * O modelo é dado por `$$\log(\mu_{ijk})=\lambda + \lambda_i^A + \lambda_j^B + \lambda_k^C + \lambda_{ij}^{AB} + \lambda_{jk}^{BC} + \lambda_{ik}^{AC}.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - XI ``` r berkeley.homog = glm(Freq~(Admit + Gender + Dept)^2, family=poisson(), data=berkeley) berkeley.homog$coef ``` ``` ## (Intercept) AdmitRejected GenderMale ## 3.137357882 2.624558572 -0.003730812 ## DeptA DeptB DeptC ## 1.135552319 -0.342488552 2.222782088 ## DeptD DeptE AdmitRejected:GenderMale ## 1.743871669 1.480917951 0.099870088 ## AdmitRejected:DeptA AdmitRejected:DeptB AdmitRejected:DeptC ## -3.306480056 -3.263082125 -2.043882034 ## AdmitRejected:DeptD AdmitRejected:DeptE GenderMale:DeptA ## -2.011873587 -1.567174318 2.002319157 ## GenderMale:DeptB GenderMale:DeptC GenderMale:DeptD ## 3.077139541 -0.662813561 0.043994838 ## GenderMale:DeptE ## -0.792866728 ``` ``` r berkeley.homog$dev ``` ``` ## [1] 20.20428 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - XII * O modelo implica que a associação condicional entre sexo e departamento não depende o valor fixado de admissão. A associação condicional entre sexo e admissão não dependem do departamento fixado, e a associação entre departamento e admissão não dependem do sexo. * O interesse está nas interações de ordem superior. Por exemplo, o primeiro coeficiente do cunjunto "AdmitRejected:GenderMale", é o log odds ratio condicional entre sexo e admissão. Nesse caso, para um departamento fixado, o odds de um homem ser rejeitado é `\(\exp(0.0999)=1.1506\)` vezes o odds de uma mulher ser rejeitada. * Embora nessa interpretação o departamento esteja fixado (dessa forma a comparação é feita entre indivíduos no mesmo departamento), não interessa em qual departamento estamos focando, pois todos conduzem ao mesmo resultado sob esse modelo. * O modelo ainda não parece ajustar adequadamente (Deviance = 20.204, com 5 gl). Mas comparado com os anteriores, parece ser melhor. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Seleção do Modelo * Considerando que os modelos aqui estudados são hierárquicos, onde cada um é caso especial de outro, usamos a estatística de razão de verossimilhanças para medir a redução no ajuste do modelo menor (hipótese nula), com respeito ao modelo maior (hipótese alternativa). Os graus de liberdade desses testes são a diferença entre os números de parâmetros envolvidos nos dois modelos. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - XIII Resumo de todos os possíveis modelos | Modelo | gl | `\(G^2\)` | valor p | |:-----------:|:----------:|:----------:|:----------:| | (D, S, A) | 16 | 2097.671 | <0.01 | | (DS, A) | 11 | 877.056 | <0.01 | | (D, SA) | 15 | 2004.222 | <0.01 | | (DA, S) | 11 | 1242.350 | <0.01 | | (DS, SA) | 10 | 783.607 | <0.01 | | (DS, DA) | 6 | 21.736 | <0.01 | | (DA, SA) | 10 | 1148.901 | <0.01 | | (DS, DA, SA) | 5 | 20.204 | <0.01 | | (DSA) | 0 | 0.00 | | --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - XIV **Observações:** * O modelo saturado seria o melhor, considerando que o qualquer modelo reduzido tem ajuste significativamente pior. * Por exemplo, para testar (DS, DA), que assume que sexo e admissão são condicionalmente independentes, dado o departamento, as hipóteses seriam `$$\begin{align*} H_0&: (DS, DA)\\ H_1&: (DS, DA, SA) \end{align*}$$` A estatística de teste pode ser calculada diretamente das estatísticas deviance dos modelos ajustados. Neste caso, `$$G^2=783.607-20.204=763.4$$` Com distribuição `\(\chi^2\)` com `\(10-5=5\)` graus de liberdade, que é altamente significante. Caso o modelo de independência condicional não tivesse sido rejeitado, esse poderia ser colocado na hipótese alternativa para testar o modelo reducido de independência conjunta (DS, A), e assim por diante. --- class: animated, hide-logo, bounceInDown ## Política de proteção aos direitos autorais > <span style="color:grey">O conteúdo disponível consiste em material protegido pela legislação brasileira, sendo certo que, por ser o detentor dos direitos sobre o conteúdo disponível na plataforma, o **LECON** e o **NEAEST** detém direito exclusivo de usar, fruir e dispor de sua obra, conforme Artigo 5<sup>o</sup>, inciso XXVII, da Constituição Federal e os Artigos 7<sup>o</sup> e 28<sup>o</sup>, da Lei 9.610/98. A divulgação e/ou veiculação do conteúdo em sites diferentes à plataforma e sem a devida autorização do **LECON** e o **NEAEST**, pode configurar violação de direito autoral, nos termos da Lei 9.610/98, inclusive podendo caracterizar conduta criminosa, conforme Artigo 184<sup>o</sup>, §1<sup>o</sup> a 3<sup>o</sup>, do Código Penal. É considerada como contrafação a reprodução não autorizada, integral ou parcial, de todo e qualquer conteúdo disponível na plataforma.</span> .pull-left[ <img src="images/logo_lecon.png" width="50%" style="display: block; margin: auto;" /> ] .pull-right[ <img src="images/logo_neaest.png" width="50%" style="display: block; margin: auto;" /> ] <br></br> .center[ [https://lecon.ufes.br](https://lecon.ufes.br/)         [https://analytics.ufes.br](https://analytics.ufes.br) ] <font size="2"><span style="color:grey">Material elaborado pela equipe LECON/NEAEST: Alessandro J. Q. Sarnaglia, Bartolomeu Zamprogno, Fabio A. Fajardo, Luciana G. de Godoi e Nátaly A. Jiménez.</span></font>