class: center, middle, inverse, title-slide .title[ # MODELOS LINEARES GENERALIZADOS STA13829 ] .subtitle[ ## Regressão Logística Binária ] .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. --> # Características * Melhor opção para os casos onde a variável dependente é binária. * Identifica-se o **<span style="color:orange">sucesso</span>** como o resultado mais importante da resposta ou aquele que se pretende relacionar com as demais variáveis de interesse. * O modelo pode ser estendido quando a variável resposta qualitativa tem mais do que duas categorias. Por exemplo, pressão sanguínea classificada como alta, normal e baixa. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo 1 - I Em um estudo prospectivo para se quantificar a influência de fatores de risco na ocorrência de doenças cardíacas, a variável resposta será dicotômica (presença ou ausência da doença) e variáveis como peso relativo, idade, sexo, nível de colesterol, hemoglobina, consumo de cigarro serão controladas. * A análise estatística deve responder às perguntas - Quais fatores de risco são os mais importantes? - É possível prever com base neles se um paciente está mais ou menos suscetível a contrair uma doença cardíaca? --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo 1 - II ``` r library(aplore3); head(chdage, n=15) ``` ``` ## id age agegrp chd ## 1 1 20 20-39 No ## 2 2 23 20-39 No ## 3 3 24 20-39 No ## 4 4 25 20-39 No ## 5 5 25 20-39 Yes ## 6 6 26 20-39 No ## 7 7 26 20-39 No ## 8 8 28 20-39 No ## 9 9 28 20-39 No ## 10 10 29 20-39 No ## 11 11 30 30-34 No ## 12 12 30 30-34 No ## 13 13 30 30-34 No ## 14 14 30 30-34 No ## 15 15 30 30-34 No ``` -- A tabela mostra os primeiros 15 casos de um conjunto de 100 indivíduos analisados em um estudo que relacionou a idade em anos (Age) com presença (1)/ausência (0) de evidências significantes de doença coronária (CDH). --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo 1 - III ``` r plot(as.integer(chd)-1 ~ age, pch = 20, main = "Diagrama de dispersão CHD x Idade", ylab = "CHD", xlab = "Idade (anos)", data = chdage) ``` <img src="Modelo_logistico_files/figure-html/unnamed-chunk-2-1.png" width="35%" style="display: block; margin: auto;" /> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo 1 - IV ####Tabela de frequências relativas por grupo ``` ## chd ## agegrp No Yes Sum ## 20-39 9 1 10 ## 30-34 13 2 15 ## 35-39 9 3 12 ## 40-44 10 5 15 ## 45-49 7 6 13 ## 50-54 3 5 8 ## 55-59 4 13 17 ## 60-69 2 8 10 ## Sum 57 43 100 ``` -- >**Observação:** À medida que a idade aumenta, a proporção de indivíduos com evidências de doenças coronárias aumenta. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Exemplo 1 - V ``` ## 20-39 30-34 35-39 40-44 45-49 50-54 55-59 60-69 ## 0.1000000 0.1333333 0.2500000 0.3333333 0.4615385 0.6250000 0.7647059 0.8000000 ``` <img src="Modelo_logistico_files/figure-html/unnamed-chunk-4-1.png" width="38%" style="display: block; margin: auto;" /> -- >**Observação**: O gráfico assemelha-se ao da função de distribuição acumulada de uma variável aleatória. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ### Observações * Em regressão, o principal interesse é relacionado ao valor médio da variável resposta, dado o valor da variável independente. Esta quantidade é a média condicional `\(E(Y|x)\)`, onde `\(Y\)` denota a variável resposta e `\(x\)` denota o valor da variável independente. -- * Em regressão linear, assume-se que a média condicional pode ser expressa como uma equação linear em `\(x\)` (ou alguma transformação de `\(x\)` ou `\(Y\)`) tais que `$$\textrm{E}(Y|x)=\beta_0+\beta_1 x.$$` Isto implica que é possível assumir qualquer valor para `\(E(Y|x)\)` quando `\(x\)` varia entre `\(-\infty\)` e `\(\infty\)`. -- * Com dados dicotômicos a média condicional deve estar entre 0 e 1 `\((0\leq \textrm{E}(Y|x)\leq 1)\)`. -- * Considere-se a quantidade `\(\pi(x)=\textrm{E}(Y|x)\)` para representar a média condicional de `\(Y\)` dado `\(x\)`, quando a distribuição logística é usada. A forma específica do modelo de regressão logística que aqui se usa é: `$$\pi(x)=\dfrac{e^{\beta_0+\beta_1 x}}{1+e^{\beta_0+\beta_1 x}}.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Transformação Logit A transformação de `\(\pi(x)\)` é denominada **<span style="color:orange">logit</span>**. Esta transformação é definida como `$$\begin{align*} g(x)&=\ln \left(\dfrac{\pi(x)}{1-\pi(x)}\right)\\ &=\beta_0+\beta_1 x \end{align*}$$` -- > **Observação**: A função `\(g(x)\)` tem várias propriedades desejáveis do modelo de regressão linear. O logit é linear em seus parâmetros, pode ser contínuo e pode variar de `\(-\infty\)` e `\(\infty\)`, dependendo do domínio de valores de `\(x\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Distribuição condicional da variável resposta `$$y=E(Y|x)+\epsilon,$$` * **Modelo de regressão linear**: assume-se que `\(\epsilon\sim \textrm{N}(0,\sigma^2)\)`. Então a distribuição condicional da variável resposta dado `\(x\)` será `\(N(\textrm{E}(Y|x),\sigma^2)\)`. -- * **Modelo de regressão logístico**: o erro só pode assumir dois valores: `$$\begin{align*} \textrm{Se } y&=1, \textrm{ então } \epsilon=1-\pi(x), \textrm{ com probabilidade } \pi(x)\\ \textrm{Se } y&=0, \textrm{ então } \epsilon=-\pi(x), \textrm{ com probabilidade } 1-\pi(x). \end{align*}$$` Então `\(\epsilon\)` tem distribuição com média zero e variância `\(\pi(x)(1-\pi(x))\)`. Ou seja, a distribuição condicional da variável resposta segue a distribuição binomial, com probabilidade dada pela média condicional `\(\pi(x)\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Regressão logística simples - I Consideramos o modelo logístico linear simples em que `\(\pi(x)\)`, a probabilidade de sucesso dado o valor `\(x\)` de uma variável explicativa qualquer, é definida tal que: `$$\ln \left(\dfrac{\pi(x)}{1-\pi(x)}\right)=\beta_0+\beta_1 x,$$` em que `\(\beta_0\)` e `\(\beta_1\)` são parâmetros desconhecidos. -- Então, `$$\pi(x)=\dfrac{e^{\beta_0+\beta_1 x}}{1+e^{\beta_0+\beta_1 x}}.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Regressão logística simples - II **Observações**: * Quando `\(x\rightarrow \infty\)`, `\(\pi(x)\uparrow 1\)`, se `\(\beta_1 > 0\)`. * Quando `\(x\rightarrow \infty\)`, `\(\pi(x)\downarrow 0\)`, se `\(\beta_1 < 0\)`. -- * Denomina-se **<span style="color:orange">chance</span>** o termo `$$\frac{\pi(x)}{1-\pi(x)}.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Relação entre a variável resposta e um preditor - I <img src="images/logisticasimples.png" width="80%" style="display: block; margin: auto;" /> --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Relação entre a variável resposta e um preditor - II * O modelo pode ser aplicado para analisar a associação entre uma determinada doença e a ocorrência ou não de um fato particular. -- * São amostrados, independentemente, `\(n_1\)` indivíduos com presença do fator `\((x = 1)\)` e `\(n_2\)` indivíduos com ausência do fator `\((x = 0)\)`. `\(\pi(x)\)` é a probabilidade do desenvolvimento da doença após um certo período fixo. -- * A chance de desenvolvimento da doença para um indivíduo com presença do fator fica dada por: `$$\dfrac{\pi(1)}{1-\pi(1)}=e^{\beta_0+\beta_1}$$` -- * A chance de desenvolvimento da doença para um indivíduo com ausência de fator é: `$$\dfrac{\pi(0)}{1-\pi(0)}=e^{\beta_0}$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Relação entre a variável resposta e um preditor - III * A razão de chances é dada por `$$\Psi=\dfrac{\frac{\pi(1)}{1-\pi(1)}}{\frac{\pi(0)}{1-\pi(0)}}=\dfrac{\pi(1)(1-\pi(0))}{\pi(0)(1-\pi(1))}=e^{\beta_1}.$$` -- * Mesmo que a amostragem seja retrospectiva, isto é, mesmo que sejam amostrados `\(n_1\)` indivíduos doentes e `\(n_2\)` indivíduos não doentes, o resultado acima continua valendo. -- >**Observação**: Essa é uma das grandes vantagens da regressão logística, a possibilidade de interpretação direta dos coeficientes como medidas de associação. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Verossimilhança - I A distribuição de probabilidade (Bernoulli) é dada por `$$f(Y_i|x_i)=\pi(x_i)^{y_i}(1-\pi(x_i))^{1-y_i},\quad y_i=0,1;\quad i=1,2,\ldots,n.$$` Dado que o pressuposto é de observações independentes, a função de verossimilhança fica dada por `$$\textrm{L}(\beta)=\prod_{i=1}^n \pi(x_i)^{y_i}(1-\pi(x_i))^{1-y_i}$$` -- Aplicando logaritmo tem-se `$$\log \textrm{L}(\beta)=\sum_{i=1}^n \left[y_i \log \left(\dfrac{\pi(x_i)}{1-\pi(x_i)}\right)\right]+\sum_{i=1}^n \log (1-\pi(x_i)).$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Verossimilhança - II Mas, `$$1-\pi(x_i)=(1+e^{\beta_0+\beta_1 x})^{-1}.$$` -- Dessa forma, `$$\log \textrm{L}(\beta_0,\beta_1)=\sum_{i=1}^n y_i (\beta_0+\beta_1 x_i)-\sum_{i=1}^n \log (1+e^{\beta_0+\beta_1 x_i}).$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Estimação dos parâmetros Não existe uma solução analítica para os valores `\(\beta_0\)` e `\(\beta_1\)` que maximizam a função de verossimilhança. Métodos numéricos são necessários para encontrar as estimativas de máxima verossimilhança, as quais denotamos, respectivamente, por `\(\hat{\beta_0}\)` e `\(\hat{\beta_1}\)`. O valor ajustado para o i-ésimo valor é dado por: `$$\hat{\pi}(x_i)=\dfrac{e^{\hat{\beta}_0 + \hat{\beta}_1 x_i}}{1+e^{\hat{\beta}_0 + \hat{\beta}_1 x_i}}$$` -- O logit estimado é dado por `$$\hat{g}(x_i)=\hat{\beta}_0+\hat{\beta}_1 x_i.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - VI ``` ## ## Call: ## glm(formula = chd ~ age, family = binomial, data = chdage) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -5.30945 1.13365 -4.683 2.82e-06 *** ## age 0.11092 0.02406 4.610 4.02e-06 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 136.66 on 99 degrees of freedom ## Residual deviance: 107.35 on 98 degrees of freedom ## AIC: 111.35 ## ## Number of Fisher Scoring iterations: 4 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - VII * Estimativas de máxima verossimilhança: `$$\hat{\beta}_0=-5,30945\qquad \hat{\beta}_1=0,11092$$` -- * Valores ajustados: `$$\hat{\pi}(x)=\dfrac{e^{-5,30945+0,11092 \times \text{Age}}}{1+e^{-5,30945+0,11092 \times \text{Age}}}$$` -- * O logit estimado é dado por `$$\hat{g}(x)=-5,30945+0,11092\times \text{Age}$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Teste da razão de verossimilhanças A comparação de valores observados e preditos usando a função de verossimilhança é baseada na expressão: `$$D=-2\ln\left(\dfrac{\textrm{verossimilhança do modelo ajustado}}{\textrm{verossimilhança do modelo saturado}}\right).$$` -- Particularmente, em regressão logística: `$$D= -2\sum_{i=1}^n \left[y_i \ln \left(\dfrac{\hat{\pi}_i}{y_i}\right)+(1-y_i)\ln \left(\dfrac{1-\hat{\pi}_i}{1-y_i}\right)\right],$$` onde `\(\hat{\pi}_i=\hat{\pi}(x_i)\)`. -- >**Observação**: A estatística `\(D\)` é chamada desvio **<span style="color:orange">deviance</span>** e faz o mesmo papel que a soma dos quadrados dos resíduos na regressão linear, ou seja, é uma estatística que auxilia na comparação de valores observados e preditos. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Saturado - I É o modelo mais complexo possível. Este modelo tem um parâmetro para cada observação e, portanto, tem um ajuste perfeito. Nesse caso, `\(\hat{\pi}_i=y_i\)` e a log-verossimilhança é dada por `$$\begin{align*} \log (\textrm{modelo saturado})&=\log\left(\prod_{i=1}^n y_i^{y_i}(1-y_i)^{1-y_i}\right)\\ &=\sum_{i=1}^n y_i\log(1-y_i)+(1-y_i)\log(1-y_i)\\ &=0 \end{align*}$$` -- Daí, o desvio é `$$D=-2 \ln (\textrm{verossimilhança do modelo ajustado}).$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Modelo Saturado - II Para estudar a significância de uma variável independente no modelo, comparamos o valor `\(D\)` obtido com e sem a variável independente na equação, `$$G=D(\textrm{modelo sem a variável})-D(\textrm{modelo com a variável})$$` `\(G\)` possui o mesmo papel do numerador no teste F parcial na regressão linear. Como a verossimilhança do modelo saturado é a mesma nos dois modelos, obtemos que: `$$G=-2\ln\left(\dfrac{\textrm{verossimilhança sem a variável}}{\textrm{verossimilhança com a variável}}\right).$$` -- Assintoticamente e sob a hipótese nula, tem-se que `\(G\)` tem uma distribuição qui-quadrado com 1 grau de liberdade. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - VIII Verificamos se a variável "Age" é significativa para explicar "CHD" através do teste de razão de verossimilhanças. `$$\begin{align*} \textrm{H}_0: \beta_1&=0\\ \textrm{H}_1:\beta_1&\neq 0 \end{align*}$$` -- Logo, `$$\begin{align*} G&=D(\textrm{modelo sem a variável})-D(\textrm{modelo com a variável})\\ &=136,66-107,35\\ &=29,31 \end{align*}$$` Dado que `\(\textrm{P}(\chi^2_1>29,31)=6,167658\times 10^{-8}\)`, rejeitamos `\(H_0\)`. Ou seja, a variável "Age" é significante para o modelo. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - IX **Análise do devio no R** ``` r anova(mod1.3) ``` ``` ## Analysis of Deviance Table ## ## Model: binomial, link: logit ## ## Response: chd ## ## Terms added sequentially (first to last) ## ## ## Df Deviance Resid. Df Resid. Dev ## NULL 99 136.66 ## age 1 29.31 98 107.35 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Teste de Wald Para `\(\beta_1=0\)`, `$$W=\dfrac{\hat{\beta}_1}{\widehat{ep}(\hat{\beta}_1)}$$` O valor-p é definido como `\(\textrm{P}(|Z|>|W|)\)`, onde `\(Z\)` denota uma variável aleatória com distribuição normal padrão. -- >**Observação**: Hauck & Donner (1977) estudaram o desempenho do teste e Wald e observaram que com frequência ele não rejeita a hipótese nula quando o coeficiente é significativo. Eles recomendam o uso do teste da razão de verossimilhança para verificar se realmente o coeficiente não é significativo quando o teste de Wald não rejeita a hipótese nula. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - X Verificamos se a variável "Age" é significativa para explicar "CHD" através do teste de razão Wald. `$$\begin{align*} \textrm{H}_0: \beta_1&=0\\ \textrm{H}_1:\beta_1&\neq 0 \end{align*}$$` -- Temos, `$$\begin{align*} \hat{\beta}_1&=0,11092\\ \widehat{ep}(\hat{\beta}_1)&=0,02406 \end{align*}$$` Daí, `\(W=4,610\)`. O valor-p é `\(\textrm{P}(|Z|>4,610)<0,001\)`. Logo, confirmamos que a variável "Age" é importante para o modelo. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Teste de Score Para `\(\beta_1=0\)`, `$$ST=\dfrac{\sum_{i=1}^n x_i(y_i-\bar{y})}{\sqrt{\bar{y}(1-\bar{y})\sum_{i=1}^n (x_i-\bar{x})^2}}$$` No caso univariado, este teste é baseado na distribuição condicional das equações de verossimilhança. O p-valor é definido como sendo `\(\textrm{P}(|Z|>ST)\)`, onde Z denota uma variável aleatória com distribuição normal padrão. -- **Exemplo**: No caso do exemplo CHD de Hosmer & Lemeshow, `$$ST=\dfrac{296,66}{\sqrt{3333,742}}=5,14.$$` `\(\textrm{P}(|Z|>5,14)<0,0001\)`. Logo, a variável "Age" é significante para o modelo. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - XI **Verificação da função de ligação** ``` r eta2<-(mod1.3$linear.predictors)^2 chdage.1<-cbind(chdage,eta2) mod1.4 <- glm(chd ~ age + eta2, family = binomial, data = chdage.1) anova(mod1.3,mod1.4, test = "LRT") ``` ``` ## Analysis of Deviance Table ## ## Model 1: chd ~ age ## Model 2: chd ~ age + eta2 ## Resid. Df Resid. Dev Df Deviance Pr(>Chi) ## 1 98 107.35 ## 2 97 107.29 1 0.065015 0.7987 ``` -- >Não temos evidências que a função de ligação canônica está inadequada. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - XII **Diagnóstico do modelo: Análise gráfica de resíduos.** Código diponível em https://www.ime.usp.br/~giapaula/diag_bino. <img src="Modelo_logistico_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 1 - XIII **Diagnóstico do modelo: envelope simulado.** Código disponível em https://www.ime.usp.br/~giapaula/envel_bino. <img src="Modelo_logistico_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 ## Intervalos de confiança Os intervalos de confiança para os parâmetros são baseados em seus respectivos testes de Wald. Assim, * O intervalo de `\(100(1-\alpha)\%\)` de confiança para `\(\beta_0\)` é dado por `$$\left[\hat{\beta}_0-z_{\alpha/2}\widehat{ep}(\hat{\beta}_0),\hat{\beta}_0+z_{\alpha/2}\widehat{ep}(\hat{\beta}_0)\right].$$` -- * O intervalo de `\(100(1-\alpha)\%\)` de confiança para `\(\beta_1\)` é dado por `$$\left[\hat{\beta}_1-z_{\alpha/2}\widehat{ep}(\hat{\beta}_1),\hat{\beta}_1+z_{\alpha/2}\widehat{ep}(\hat{\beta}_1)\right].$$` Onde `\(z_{\alpha/2}\)` é o quantil de uma normal padrão dado por `\(\textrm{P}(z>z_{\alpha/2})=\alpha/2\)`, `\(\widehat{ep}(\hat{\beta}_0)\)` e `\(\widehat{ep}(\hat{\beta}_1)\)` denotam o estimador baseado no desvio padrão de `\(\hat{\beta}_0\)` e `\(\hat{\beta}_1\)`, respectivamente. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Intervalo de confiança para o logit - I O logit é a parte linear do modelo de regressão logístico, dado por `$$\hat{g}(x)=\hat{\beta}_0+\hat{\beta}_1 x.$$` -- O estimador da variância do logit é dado por `$$\widehat{\textrm{Var}}\left[\hat{g}(x)\right]=\widehat{\textrm{Var}}(\hat{\beta_0})+x^2 \widehat{\textrm{Var}}(\hat{\beta_1})+2x\widehat{\textrm{Cov}}(\hat{\beta_0},\hat{\beta_1}),$$` onde a matriz de variância-covariância de `\(\hat{\beta}_0\)` e `\(\hat{\beta}_1\)` corresponde ao inverso da matriz de informação de Fisher, calculada em `\(\hat{\beta}_0\)` e `\(\hat{\beta}_1\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Intervalo de confiança para o logit - II **Intervalo de confiança de `\(100(1-\alpha)\%\)` para `\(g(x)\)`**: `$$\left[\hat{g}(x)-z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(x)]},\hat{g}(x)+z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(x)]}\right].$$` -- O intervalo de `\(100(1-\alpha)\%\)` de confiança para `\(\pi(x)\)` é dado por `$$\left[\frac{\exp \left[\hat{g}(x)-z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(x)]}\right]}{1+\exp \left[\hat{g}(x)-z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(x)]}\right]},\frac{\exp \left[\hat{g}(x)+z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(x)]}\right]}{1+\exp \left[\hat{g}(x)+z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(x)]}\right]}\right].$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - XIV ``` r confint(mod1.3, level=.95, type="LR") ``` ``` ## Waiting for profiling to be done... ``` ``` ## 2.5 % 97.5 % ## (Intercept) -7.72587162 -3.2461547 ## age 0.06693158 0.1620067 ``` ``` r exp(coef(mod1.3))[-1] ``` ``` ## age ## 1.117307 ``` ``` r exp(confint(mod1.3))[-1, ] ``` ``` ## Waiting for profiling to be done... ``` ``` ## 2.5 % 97.5 % ## 1.069222 1.175868 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - XV * Logit para uma pessoa de 50 anos: `$$\hat{g}(50)=-5,31 + 0,1111\times 50=0,240.$$` -- * Variância estimada: `$$\widehat{\textrm{Var}}\left[\hat{g}(50)\right]=1,2852+(50)^2 \times 0,00058+2(50)\times (-0,0267)=0,0650$$` -- * Intervalo de `\(95\%\)` de confiança para o logit: `$$0,240\pm 1,96\sqrt{0,0650}=[-0,260,0,740]$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 1 - XVI Temos `$$\hat{\pi}(50)=\dfrac{e^{\hat{g}(50)}}{1+e^{\hat{g}(50)}}=\dfrac{e^{-5,31 + 0,1111\times 50}}{1+e^{-5,31 + 0,1111\times 50}}=0,56$$` -- O intervalo de confiança para o valor ajustado é dado por `$$\left[\frac{\exp \left[0,240-1,96\times 0,2550\right]}{1+\exp \left[0,240-1,96\times 0,2550\right]};\frac{\exp \left[0,240+1,96\times 0,2550\right]}{1+\exp \left[0,240+1,96\times 0,2550\right]}\right]$$` `$$\left[\frac{\exp \left[-0,260\right]}{1+\exp \left[-0,260\right]};\frac{\exp \left[0,740\right]}{1+\exp \left[0,740\right]}\right]$$` `$$\left[0,435;0,677\right]$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn # Regressão Logística com várias variáveis preditoras O logit do modelo de regressão logístico múltiplo é dado por: `$$g(\mathbf{x})=\beta_0+\beta_1 x_1+\cdots +\beta_p x_p,$$` em que `\(\mathbf{x}=(x_1,x_2,\ldots,x_p)^t\)` e o modelo de regressão logística é `$$\pi(\mathbf{x})=\dfrac{e^{g(\mathbf{x})}}{1+e^{g(\mathbf{x})}}$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Verossimilhança A obtenção da função de verossimilhança se dá de maneira similar à encontrada para o modelo de regressão logístico simples. Assim como naquele caso, em modelo de regressão logístico múltiplo, são necessários métodos numéricos para obter as estimativas de máxima verossimilhança do vetor de parâmetros regressores. A função resposta logística ajustada e os valores ajustados são dados por: `$$\begin{align*} \hat{\pi}(\mathbf{x})&=\dfrac{\exp(\hat{\beta}^t\mathbf{x})}{1+\exp(\hat{\beta}^t\mathbf{x})}=[1+\exp(-\hat{\beta}^t\mathbf{x})]^{-1}\\ \hat{\pi}(x_i)&=\dfrac{\exp(\hat{\beta}^t x_i)}{1+\exp(\hat{\beta}^t x_i)}=[1+\exp(-\hat{\beta}^t x_i)]^{-1} \end{align*}$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - I Um estudo 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 setor da cidade. A idade `\((X_1)\)` é uma variável quantitativa. O status sócio-econômico é uma variável com 3 categorias. Esta variável é representada por duas variáveis indicadoras, `\(X_2\)` e `\(X_3\)`, assim: | Classe || `\(X_2\)` | `\(X_3\)` || |:-------:||--------:|:-----------:|| | Alta || 0 | 0 || | Média || 1 | 0 || | Baixa || 0 | 1 || --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - II A variável setor da cidade também é uma variável categorizada. Como existiam apenas dois setores na cidade, uma variável indicadora `\((X_4)\)` foi usada, definida como `\(X=0\)` para o setor 1 e `\(X=1\)` para o setor 2. A razão para a escolha da classe social alta ser tomada como referência é que se espera que esta classe tenha as menores taxas de doença entre as classes sociais. Fazendo-se esta classe como referência, espera-se que a razão das chances associada aos coeficientes de regressão 2 e 3 sejam maiores do que 1, facilitando a interpretação. Pela mesma razão, o setor 1, onde a epidemia foi menos severa, foi escolhido como referência para a variável indicadora `\(X_4\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - III
--- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - IV ####Modelo Ajustado ``` ## ## Call: ## glm(formula = Dengue ~ Idade + Classe + Setor, family = binomial, ## data = dengue) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -2.293933 0.436769 -5.252 1.5e-07 *** ## Idade 0.026991 0.008675 3.111 0.001862 ** ## Classe2 0.044609 0.432490 0.103 0.917849 ## Classe3 0.253433 0.405532 0.625 0.532011 ## Setor2 1.243630 0.352271 3.530 0.000415 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 236.33 on 195 degrees of freedom ## Residual deviance: 211.22 on 191 degrees of freedom ## AIC: 221.22 ## ## Number of Fisher Scoring iterations: 3 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - V #### Razões de chances estimadas ``` r exp(cbind(OR=coef(fit.model1), confint(fit.model1))) ``` ``` ## Waiting for profiling to be done... ``` ``` ## OR 2.5 % 97.5 % ## (Intercept) 0.1008689 0.04094164 0.2286073 ## Idade 1.0273586 1.01029414 1.0454528 ## Classe2 1.0456186 0.44264364 2.4354849 ## Classe3 1.2884413 0.58303368 2.8813353 ## Setor2 3.4681814 1.75627062 7.0281934 ``` -- **Alguma interpretações**: * `\(\hat{\psi}=1,03\)` (A chance de uma pessoa estar doente aumenta cerca de 3% com cada ano adicional de idade `\((X_1)\)`, para dado status sócio-econômico e setor da cidade (constantes).) -- * `\(\hat{\psi}=3,47\)` (A chance de uma pessoa no setor 2 `\((X_4)\)` contrair a doença é quase 4 vezes maior que para uma pessoa do setor 1, dado a idade e o status sócio-econômico constantes.) --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - VI #### **Probabilidades estimadas** `$$\hat{\pi}=(1+\exp (-2,293+0,027X_1+0,044X_2+0,0,253X_3+1,243X_4))^{-1}.$$` -- >Por exemplo, o valor ajustado para o caso `\(i=1\)`, onde `\(X_{11}=33\)`, `\(X_{12}=0\)`, `\(X_{13}=0\)`, `\(X_{14}=0\)` é `\(\hat{\pi}=0,197\)`. -- >**Interpretação**: A probabilidade estimada de uma pessoa com 33 anos de idade, da classe alta, do setor 1, contrair a doença é de aproximadamente 0,20. -- **No R:** ``` r perfil<-data.frame(Idade=33,Classe="1",Setor="1") pred.prob<-predict(fit.model1,perfil,type = "response") pred.prob ``` ``` ## 1 ## 0.197304 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - VI ####Teste de razão de verossimilhanças As hipóteses de interesse são: `$$\begin{align*} \textrm{H}_0&:\mathbf{\beta}=\mathbf{\beta}^0\\ \textrm{H}_1&:\mathbf{\beta}\neq\mathbf{\beta}^0, \end{align*}$$` em que `\(\mathbf{\beta}^0\)` é um vetor `\((p+1)-\)`dimensional. -- De forma similar ao modelo de regressão logístico simples, considera-se a estatística `\(G\)`. Para verificar a significância dos coeficientes, considera-se, sob a hipótese nula, o vetor `\(\mathbf{\beta}^0\)` igual `\(\mathbf{0}\)`. -- Assintoticamente, `\(G\)` tem distribuição qui-quadrado com `\(p+1\)` graus de liberdade. Daí, rejeita-se `\(H_0\)` se `\(P(\chi^2_{p+1}>G)<\alpha\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - VII ####Teste de Wald O análogo do Teste de Wald é obtido pela expressão `$$W=\mathbf{\beta}'\left[\widehat{\textrm{Var}}(\hat{\mathbf{\beta}})\right]^{-1}\mathbf{\beta}=\mathbf{\beta}'(\mathbf{X}\Sigma\mathbf{X})\mathbf{\beta},$$` onde `\(\Sigma\)` é dada por `$$\left[\begin{array}{cccc} \pi_1(1-\pi_1)&0&\cdots&0\\ 0&\pi_2(1-\pi_2)&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\pi_n(1-\pi_n) \end{array}\right].$$` -- Sob `\(H_0\)` a estatística `\(W\)` possui distribuição qui-quadrado com `\(p+1\)` graus de liberdade. Como este teste exige a execução de operações matriciais e a obtenção de `\(\mathbf{\beta}\)`, então não há vantagens computacionais sobre o teste de razão de verossimilhança para testar a significância do modelo. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - VIII ####Intervalo de confiança para o logit `$$\hat{g}(\mathbf{X})=\hat{\mathbf{\beta}}\mathbf{X}.$$` O estimador da variância do logit é dado por `$$\widehat{\textrm{Var}}\left[\hat{g}(\mathbf{x})\right]=\sum_{j=0}^p x_j^2\widehat{\textrm{Var}}(\hat{\beta_j})+2\sum_{j=0}^p\sum_{k=j+1}^p x_jx_k\widehat{\textrm{Cov}}(\hat{\beta_j},\hat{\beta_k}).$$` -- Equivalentemente, `$$\widehat{\textrm{Var}}\left[\hat{g}(\mathbf{x})\right]=\mathbf{x}'\widehat{\textrm{Var}}(\mathbf{\beta})\mathbf{x}=\mathbf{x}'(\mathbf{X}'\Sigma\mathbf{X})^{-1}\mathbf{x}.$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - IX **Intervalo de confiança de `\(100(1-\alpha)\%\)` para `\(g(x)\)`**: `$$\left[\hat{g}(\mathbf{x})-z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(\mathbf{x})]},\hat{g}(\mathbf{x})+z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(\mathbf{x})]}\right].$$` -- O intervalo de `\(100(1-\alpha)\%\)` de confiança para `\(\pi(\mathbf{x})\)` é dado por `$$\left[\frac{\exp \left[\hat{g}(\mathbf{x})-z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(\mathbf{x})]}\right]}{1+\exp \left[\hat{g}(\mathbf{x})-z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(\mathbf{x})]}\right]},\frac{\exp \left[\hat{g}(\mathbf{x})+z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(\mathbf{x})]}\right]}{1+\exp \left[\hat{g}(\mathbf{x})+z_{\alpha/2}\sqrt{\widehat{\textrm{Var}}[\hat{g}(\mathbf{x})]}\right]}\right].$$` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Teste de qualidade de ajuste de Hosmer e Lemeshow - I * Definida comparando-se o número observado com o número esperado de sucessos de `\(g\)` grupos formados. -- * O primeiro grupo deve conter `\(n_1'\)` elementos correspondentes às `\(n_1'\)` menores probabilidades ajustadas, as quais são denotadas por `\(\hat{\pi}_{(1)}\leq \hat{\pi}_{(2)}\leq \ldots \leq \hat{\pi}_{(n_1')}\)`. O segundo grupo deve conter os `\(n_2'\)` elementos correspondentes às seguintes probabilidades ajustadas `\(\hat{\pi}_{(n_1'+1)}\leq \hat{\pi}_{(n_1'+2)}\leq \ldots \leq \hat{\pi}_{(n_1'+n_2')}\)`. E assim, até o último grupo, que deve conter as `\(n_g'\)` maiores probabilidades ajustadas `\(\hat{\pi}_{(n_1'+\ldots+n_{g-1}'+1)}\leq \hat{\pi}_{(n_1'+\ldots+n_{g-1}'+2)}\leq \ldots \leq \hat{\pi}_{(n)}\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Teste de qualidade de ajuste de Hosmer e Lemeshow - II * O número observado de sucessos no primeiro grupo formado é dado por `\(O_1=\sum_{j=1}^{n_1'} y_{(j)}\)`, em que `\(y_{(j)}=0\)` se o elemento correspondente é fracasso e `\(y_{(j)}=1\)` se é sucesso. -- * Em geral, `$$O_i=\sum_{j=n_1'+\ldots+n_{i-1}+1}^{n_1'+\ldots+n_{i}} y_{(j)},\quad 2\leq i\leq g.$$` -- * A estatística de Hosmer-Lemeshow é definida como `$$\hat{C}=\sum_{i=1}^g \dfrac{(O_i-n_i'\hat{\pi}_i)^2}{n_i'\hat{\pi}_i(1-\hat{\pi}_i)},$$` em que `\(\hat{\pi}_1=(1/n_1')\sum_{j=1}^{n_1'}\hat{\pi}_{(j)}\)` e `\(\hat{\pi}_i=(1/n_i')\sum_{n_1'+\cdots+n_{i-1}'+1}^{n_1'+\cdots+n_{i}'}\hat{\pi}_{(j)}\)`, `\(2\leq i\leq g\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Teste de qualidade de ajuste de Hosmer e Lemeshow - III * Hosmer e Lemeshow sugerem a formação de `\(g=10\)` grupos de mesmo tamanho (aproximadamente), de modo que o primeiro grupo contenha `\(n_1'\)` elementos correspondentes às `\((n/10)\)` menores probabilidades, e assim por diante até o último grupo com `\(n_{10}'\)` elementos correspondentes às `\((n/10)\)` maiores probabiidades. -- * Quando não há empates é relativamente simples montar os 10 grupos com tamanhos aproximadamente iguais. -- * Quando há empates, pode ser necessário que dois indivíduos com a mesma configuração de covariáveis sejam alocados em grupos adjacentes a fim de que os grupos formados não tenham tamanhos muito desiguais. -- * A distribuição assintótica de `\(\hat{C}\)` é `\(\chi^2_{(g-2)}\)`. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - X ``` r require(ResourceSelection) ``` ``` ## Carregando pacotes exigidos: ResourceSelection ``` ``` ## ResourceSelection 0.3-6 2023-06-27 ``` ``` r hosmer=hoslem.test(dengue$Dengue,fitted(fit.model1),g=10) hosmer ``` ``` ## ## Hosmer and Lemeshow goodness of fit (GOF) test ## ## data: dengue$Dengue, fitted(fit.model1) ## X-squared = 196, df = 8, p-value < 2.2e-16 ``` --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - XI **Verificação da função de ligação** ``` r eta2<-(fit.model1$linear.predictors)^2 dengue.1<-cbind(dengue,eta2) fit.model2 <- glm(Dengue ~ Idade + factor(Classe) + factor(Setor) + eta2, family = binomial, data = dengue.1) anova(fit.model1,fit.model2, test = "LRT") ``` ``` ## Analysis of Deviance Table ## ## Model 1: Dengue ~ Idade + Classe + Setor ## Model 2: Dengue ~ Idade + factor(Classe) + factor(Setor) + eta2 ## Resid. Df Resid. Dev Df Deviance Pr(>Chi) ## 1 191 211.22 ## 2 190 210.42 1 0.8024 0.3704 ``` -- >Não temos evidências que a função de ligação canônica está inadequada. --- [//]: <> (class: center, middle, animated, slideInRight/ class: animated slideInRight fadeOutLeft) class: animated, fadeIn ## Exemplo 2 - XII **Diagnóstico do modelo: Análise gráfica de resíduos.** <img src="Modelo_logistico_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 ## Exemplo 1 - XIII **Diagnóstico do modelo: envelope simulado.** <img src="Modelo_logistico_files/figure-html/unnamed-chunk-19-1.png" width="45%" style="display: block; margin: auto;" /> --- 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>