4.4. Regressão Logística Simples
4.1.1 Modelo Estatístico
Um modelo de regressão logística simples é usado para o caso de regressão com uma variável explicativa.
Suponha uma amostra de $n$ observações independentes da terna $(x_i,m_i,y_i),$ $i=1,2,\ldots,n$, sendo que:
-
$x_i$ é o valor da variável explicativa;
-
$m_i$ é a quantidade de itens verificados na amostra (número de ensaios);
-
$y_i$ número de ocorrência de um evento (exemplo: quantidade de peças não conforme) em $m_i$ ensaios; e
-
$n$ é o tamanho da amostra.
Com isso, assumimos que a variável resposta tem distribuição de probabilidade binomial $( Y_i \sim B(m_i,\pi_i))$, tal que
$$P[Y_i=y_i]=\binom{m_i}{y_i}\pi_i^{y_i}(1-\pi_i)^{m_i-y_i}.$$
Para adequarmos a resposta média ao modelo linear usamos a função de ligação
$$\pi(x_i)=\frac{e^{\beta_0+\beta_1x_i}}{1+e^{\beta_0+\beta_1x_i}}, i=1, \ldots, n,$$
que pode ser escrita como
$$\ln\left(\frac {\pi_i} {1 - \pi_i}\right)= \beta_0+\beta_1x_i.$$
As figuras a seguir ilustram a forma do modelo logístico para $\beta_1$ positivo e negativo.
$$\pi(x_i)=\frac{e^{10-0.1x_i}}{1+e^{10-0.1x_i}}$$
Figura 4.4.1: Modelo logístico com $\beta_1$ positivo.
Figura 4.4.2: Modelo logístico com $\beta_1$ negativo.
Neste caso, utilizamos o método da máxima verossimilhança para estimarmos os parâmetros $(\beta_0 ,\beta_1)$. De forma genérica, o método de máxima verossimilhança nos fornece valores para os parâmetros desconhecidos que maximizam a probabilidade de se obter determinado conjunto de dados.
Assumindo que $(x_0, m_0 , y_0), \ldots, (x_n, m_n, y_n)$ são independentes, a função de verossimilhança é da seguinte forma
$$P [Y=y_1,\ldots,y_n|\beta_0,\beta_1]~=~\prod_{i=1}^n \binom{m_i} {y_i} \pi_i^{y_i}(1-\pi_i)^{m_i-y_i}$$
$$=\prod_{i=1}^n \binom{m_i} {y_i} \pi_i^{y_i}(1-\pi_i)^{m_i}(1-\pi_i)^{-y_i}$$
$$=\prod_{i=1}^n \binom{m_i} {y_i} \frac{\pi_i^{y_i}(1-\pi_i)^{m_i}}{(1-\pi_i)^{y_i}}$$
$$=\prod_{i=1}^n \binom{m_i} {y_i} \left(\frac{\pi_i}{1-\pi_i}\right)^{y_i}{(1-\pi_i)^{m_i}}$$
Ignorando o termo constante $\displaystyle\binom{m_i} {y_i},$ que não depende de $x_i,$ e tomando o logaritmo $(\ln)$ em ambos os lados da expressão anterior, temos
$$L~(\beta_0,\beta_1|(x_i;m_i;y_i))~=~\ln \left(\frac{\pi_i}{1-\pi_i}\right)^{y_i}+\ln {(1-\pi_i)^{m_i}}$$
$$L~(\beta_0,\beta_1|(x_i;m_i;y_i))~=~{y_i} \ln \left(\frac{\pi_i}{1-\pi_i}\right)+ {m_i} \ln(1-\pi_i) \tag{4.1.1.1}$$
Detalhando $\ln \left(\displaystyle\frac{\pi_i}{1-\pi_i}\right),$ e considerando que,
$$\pi_i~=~\frac{e^{\beta_0+\beta_1 x_i}}{1+e^{\beta_0+\beta_1 x_i}}, {temos}$$
$$\ln \left(\displaystyle\frac{\pi_i}{1-\pi_i}\right)~=~\ln \left(\displaystyle\frac{\displaystyle\frac{e^{\beta_0+\beta_1x_i}}{1+e^{\beta_0+\beta_1x_i}}}{1-\displaystyle\frac{e^{\beta_0+\beta_1 x_i}}{1+e^{\beta_0+\beta_1 x_i}}}\right)$$
$$=\ln \left(\frac{\displaystyle\frac{e^{\beta_0+\beta_1 x_i}}{1+e^{\beta_0+\beta_1 x_i}}}{\displaystyle\frac{1}{1+e^{\beta_0+\beta_1 x_i}}}\right)$$
Assim a expressão (4.1.1.1), pode ser reescrita como:
$$L~(\beta_0,\beta_1|(x_i;m_i;y_i))~=~\sum^n_{i=1}\left[ y_i ~\ln \left(\frac{\pi_i}{1-\pi_i} \right)+m_i \ln (1-\pi_i)\right]$$
$$=\sum^n_{i=1}\left[ y_i ~\left(\beta_0 + \beta_1 x_i \right)+ m_i \ln \left(1 - \frac{e^{\beta_0+\beta_1 x_i}}{1+e^{\beta_0+\beta_1x_i}}\right)\right]$$
$$=\sum^n_{i=1}\left[y_i~(\beta_0+\beta_1 x_i)+m_i \ln \left(\frac{1}{1+e^{\beta_0+\beta_1x_i}}\right)\right]$$
$$=\sum^n_{i=1}\left[y_i~(\beta_0+\beta_1 x_i) + m_i (\ln 1 - \ln (1+e^{\beta_0+\beta_1 x_i}))\right]$$
$$=\sum^n_{i=1}y_i~(\beta_0+\beta_1 x_i) - \sum^n_{i=1} m_i \ln (1+e^{\beta_0+\beta_1 x_i})$$
Portanto,
$$L~(\beta_0,\beta_1|(x_i;m_i;y_i))=\sum^n_{i=1}y_i~(\beta_0+\beta_1x_i) - \sum^n_{i=1}m_i \ln (1+e^{\beta_0+\beta_1 x_i}) \tag{4.1.1.2}$$
Para simplificar a notação faremos $L~(\beta_0,\beta_1|(x_i;m_i;y_i))= L~(\beta_0,\beta_1).$
4.1.2 Estimação dos Parâmetros do modelo
Para ajustar um modelo de regressão devemos estimar os parâmetros $\beta_0$ e $\beta_1$ do modelo.
Os estimadores de máxima verossimilhança para os parâmetros $\beta_0$ e $\beta_1$ são os valores de $\widehat{\beta_0}$ e $\widehat{\beta_1}$ que maximizam o logaritmo da função de verossimilhança. A função de verossimilhança tem máximo, pois $0 < P[Y_i = y_i \mid x_i] < 1,$ pois a função logaritmo é estritamente crescente.
Para maximizar a função de verossimilhança basta derivarmos em relação aos parâmetros do modelo, da seguinte forma
$$\frac{\partial}{\partial\beta_0} L(\beta_0,\beta_1)~=~\sum^n_{i=1}y_i-\sum^n_{i=1}{m_i} \frac{e^{\beta_0+\beta_1x_i}}{1+e^{\beta_0+\beta_1x_i}}$$
$$\frac{\partial}{\partial\beta_1} L(\beta_0,\beta_1)~=~\sum^n_{i=1}y_i~x_i-\sum^n_{i=1}{m_i x_i} \frac{e^{\beta_0+\beta_1x_i}}{1 +e^{\beta_0+\beta_1x_i}}$$
Igualando estas derivadas a zero e substituindo os parâmetros ($\beta_0,\beta_1$) pelos estimadores $(\widehat{\beta_0},\widehat{\beta_1}),$
$$\sum^n_{i=1}y_i-\sum^n_{i=1}{m_i}\dfrac{e^{\widehat{\beta_0}+\widehat{\beta_1}x_i}}{1+e^{\widehat{\beta_0}+\widehat{\beta_1}x_i}}=0$$
$$\sum^n_{i=1}y_i~x_i-\sum^n_{i=1}{m_i x_i}\dfrac{e^{\widehat{\beta_0}+\widehat{\beta_1}x_i}}{1 +e^{\widehat{\beta_0}+\widehat{\beta_1}x_i}} = 0$$
Porém estas equações são não-lineares nos parâmetros e para resolvê-las é preciso recorrer a métodos numéricos interativos, como Newton-Raphson (Gourieroux e Monfort, 1995). Este método é definido expandindo-se a função U($\pmb{\beta}$) em torno do ponto inicial $\pmb{\beta}^{(0)}$, tal que
$$U(\pmb{\beta})\approx U(\pmb{\beta}^{(0)})+U^\prime (\pmb{\beta}^{(0)})(\pmb{\beta}-\pmb{\beta}^{(0)}), \tag{4.1.2.1}$$
sendo que U($\pmb{\beta}$) são as derivadas de primeira ordem do logaritmo da função de verossimilhança em relação aos parâmetros do modelo e $U^\prime (\pmb{\beta}) $ são as derivadas de ordem 2 do logaritmo da função de verossimilhança.
Se repetirmos o processo (4.1.2.1) chegaremos ao processo iterativo
$$\pmb{\beta}^{(m+1)}=\pmb{\beta}^{(m)}+[-U^\prime (\pmb{\beta}^{(m)})]^{-1}U^\prime (\pmb{\beta}^{(m)}),$$
sendo que $m=0,1,\ldots$
Como a matriz $-U^\prime (\pmb{\beta})$ pode não ser positiva definida, e portanto não invertível, ela é substituída pela matriz de informação de Fisher. Assim
$$\pmb{\beta}^{(m+1)}=\pmb{\beta}^{(m)}+[{-I}(\pmb{\beta}^{(m)})^{-1}]U^\prime (\pmb{\beta}^{(m)}), \quad m=0,1,\ldots \tag{4.1.2.2}$$
A matriz de informação de Fisher, para o modelo logístico com uma variável, tem a seguinte forma:
$$ I(\widehat{\beta})=- \begin{bmatrix} \dfrac{\partial^2}{\partial\beta_0^2} \ln L(\beta_0,\beta_1) \quad \dfrac{\partial^2}{\partial\beta_0 \beta_1} \ln L(\beta_0,\beta_1) \cr \cr \dfrac{\partial^2}{\partial\beta_0\beta_1} \ln L(\beta_0,\beta_1) \quad \dfrac{\partial^2}{\partial \beta_1^2} \ln L(\beta_0,\beta_1) \end{bmatrix} $$
$$=\begin{bmatrix}\sum\limits_{i=1}^n m_i\dfrac{e^{\beta_0+\beta_1x_i}}{(1+e^{\beta_0+\beta_1x_i})^2}~~\sum\limits_{i=1}^n m_i x_i\dfrac{e^{\beta_0+\beta_1x_i}}{(1+e^{\beta_0+\beta_1x_i})^2} \cr \cr \sum\limits_{i=1}^n m_i x_i \dfrac{e^{\beta_0+\beta_1x_i}}{(1+e^{\beta_0+\beta_1x_i})^2} ~~\sum\limits_{i=1}^n m_i x_i^2\dfrac{e^{\beta_0+\beta_1x_i}}{(1+e^{\beta_0+\beta_1x_i})^2} \end{bmatrix} \tag{4.1.2.3}$$
Após obter as estimativas dos parâmetros do modelo é possível calcular as probabilidades estimadas
$$\widehat{\pi}_i=\dfrac{e^{\widehat{\beta_0}+\widehat{\beta_1}x_i}}{1+e^{\widehat{\beta_0}+\widehat{\beta_1}x_i}} \tag{4.1.2.4}$$
Exemplo 4.1.2.1
Realizou-se um treinamento com 30 funcionários de um determinado setor de uma empresa. O objetivo do treinamento foi determinar o menor número de horas de treinamento necessários para ocorrência do menor número de erros de montagem. A Tabela 4.4.1 mostra os dados coletados no treinamento.
| Horas de Treinamento (X) | Erros de Montagem (Y) | Peças Analisadas (M) | Horas de Treinamento (X) | Erros de Montagem (Y) | Peças Analisadas (M) |
|---|---|---|---|---|---|
| 30 | 2 | 200 | 17 | 9 | 200 |
| 30 | 2 | 200 | 17 | 10 | 200 |
| 30 | 2 | 200 | 16 | 10 | 200 |
| 29 | 2 | 200 | 15 | 11 | 200 |
| 28 | 3 | 200 | 13 | 11 | 200 |
| 27 | 4 | 200 | 12 | 12 | 200 |
| 26 | 5 | 200 | 11 | 12 | 200 |
| 26 | 5 | 200 | 11 | 12 | 200 |
| 25 | 6 | 200 | 11 | 13 | 200 |
| 24 | 6 | 200 | 10 | 13 | 200 |
| 23 | 8 | 200 | 10 | 13 | 200 |
| 20 | 8 | 200 | 9 | 13 | 200 |
| 20 | 8 | 200 | 8 | 13 | 200 |
| 20 | 8 | 200 | 8 | 14 | 200 |
| 17 | 9 | 200 | 5 | 14 | 200 |
Tabela 4.4.1: Dados do Treinamento.
Neste exemplo, a variável explicativa é denominada ($X$), a variável resposta ($Y$) e ($M$) é o número de ensaios (repetições). Como estes dados têm distribuição binomial é possível ajustar um modelo logístico.
Para estimar os parâmetros deste modelo utilizamos o método iterativo de Newton-Raphson. Como ponto inicial consideramos $ \beta = (\beta_0, \beta_1) = (0,0) $ e utilizando a expressão 4.1.2.2, tem-se os seguintes passos do processo iterativo:
[1ª] Iteração:
$$\widehat{\beta}= \begin{bmatrix}0 \cr 0 \cr \end{bmatrix}+ \begin{bmatrix}1500 \quad 27400 \cr 27400 \quad 589700 \cr \end{bmatrix}^{-1}\times\begin{bmatrix}-2742 \cr -50994 \cr \end{bmatrix}=\begin{bmatrix}-1,64228866 \cr -0,01016668 \cr \end{bmatrix}$$
[2ª] Iteração:
$$\widehat{\beta}=\begin{bmatrix}-1,64228866 \cr -0,01016668 \cr \end{bmatrix}+\begin{bmatrix}716,4305 \quad 12774,56 \cr 12774,5605 \quad 270061,59 \cr \end{bmatrix}^{-1}\times\begin{bmatrix}-574,4487 \cr -10968,1864 \cr \end{bmatrix}=\begin{bmatrix} -2,13823758 \cr -0,02732075 \cr \end{bmatrix}$$
[3ª] Iteração:
$$\widehat{\beta}=\begin{bmatrix}-2,13823758 \cr -0,02732075 \cr \end{bmatrix}+\begin{bmatrix}379,0670 \quad 6399,839 \cr 6399,8389 \quad 129715,680 \cr \end{bmatrix}^{-1} \times\begin{bmatrix}-149,8268 \cr -3036,3513 \cr \end{bmatrix}=\begin{bmatrix}-2,13856968 \cr -0,05071211 \cr \end{bmatrix}$$
[4ª] Iteração:
$$\widehat{\beta}=\begin{bmatrix}-2,13856968 \cr -0,05071211 \cr \end{bmatrix}+\begin{bmatrix}269,9642 \quad 4225,006 \cr 4225,0061 \quad 80422,963 \cr \end{bmatrix}^{-1}\times\begin{bmatrix}-27,2826 \cr -622,3350 \cr \end{bmatrix}=\begin{bmatrix}-2,02583693 \cr -0,06437278 \cr \end{bmatrix}$$
[5ª] Iteração:
$$\widehat{\beta}=\begin{bmatrix}-2,02583693 \cr -0,06437278 \cr \end{bmatrix}+\begin{bmatrix}246,7687 \quad 3697,78 \cr 3697,7797 \quad 67728,94 \cr \end{bmatrix}^{-1}\times\begin{bmatrix}-2,370660 \cr -59,028702 \cr \end{bmatrix}=\begin{bmatrix}-2,00685123 \cr -0,06628088 \cr \end{bmatrix}$$
[6ª] Iteração:
$$\widehat{\beta}=\begin{bmatrix}-2,00685123 \cr -0,06628088 \cr \end{bmatrix}+\begin{bmatrix}244,5458 \quad 3642,619 \cr 3642,6187 \quad 66355,737 \cr \end{bmatrix}^{-1}\times\begin{bmatrix} -0,03166171 \cr -0,79013813 \cr \end{bmatrix}=\begin{bmatrix}-2,00658851 \cr -0,06630721 \cr \end{bmatrix}$$
O processo foi interrompido na 6ª iteração, pois, a partir dela, o resultado se estabiliza. Assim, para os dados da Tabela 4.1.2.1 os estimadores dos parâmetros $\beta=(\beta_0,\beta_1),$ são respectivamente:
$$\widehat{\beta_0} = -2,0066$$
$$\widehat{\beta_1} = -0,0663$$
O estimador $\widehat{\beta_0}$ é o intercepto do modelo e $\widehat{\beta_1}$ é o coeficiente da variável explicativa (Horas de Treinamento).
A probabilidade estimada de ocorrer Erros de Montagem em relação as Horas de Treinamento foi obtida a partir da expressão 4.1.2.4. Estas estimativas são apresentados na Tabela 4.4.2. Verificamos pela Figura 4.4.3 que quanto maior o número de Horas em Treinamento menor será a probabilidade de erro.
| Horas de Treinamento | Probabilidade do Erro | Horas de Treinamento | Probabilidade do Erro |
|---|---|---|---|
| 30 | 0,01806045 | 17 | 0,04173392 |
| 30 | 0,01806045 | 17 | 0,04173392 |
| 30 | 0,01806045 | 16 | 0,04446776 |
| 29 | 0,01927472 | 15 | 0,04737184 |
| 28 | 0,02056892 | 13 | 0,05372869 |
| 27 | 0,02194807 | 12 | 0,05720136 |
| 26 | 0,02341749 | 11 | 0,06088404 |
| 26 | 0,02341749 | 11 | 0,06088404 |
| 25 | 0,02498277 | 11 | 0,06088404 |
| 24 | 0,02664982 | 10 | 0,06478753 |
| 23 | 0,02842486 | 10 | 0,06478753 |
| 20 | 0,03446517 | 9 | 0,06892291 |
| 20 | 0,03446517 | 8 | 0,07330157 |
| 20 | 0,03446517 | 8 | 0,07330157 |
| 17 | 0,04173392 | 5 | 0,08801434 |
Tabela 4.4.2: Probabilidade de ocorrência do erro.
Figura 4.4.3: Gráfico da probabilidade ajustada.
4.1.2.1 Interpretação dos parâmetros do modelo
A interpretação dos parâmetros de um modelo de regressão logística é obtida comparando a probabilidade de sucesso com a probabilidade de fracasso, usando a função odds ratio - OR (razão de chances). Essa função é obtida a partir da função odds.
$$g(x)=\displaystyle\frac{\pi(x)}{[1-\pi(x)]}=\displaystyle\frac{\displaystyle\frac{e^{\beta_0+\beta_1x_i}}{1+e^{\beta_0+\beta_1x_i}}}{1-\displaystyle\frac{e^{\beta_0+\beta_1x_i}}{1+e^{\beta_0+\beta_1x_i}}}= \displaystyle\frac{\displaystyle\frac{e^{\beta_0+\beta_1x_i}}{1+e^{\beta_0+\beta_1x_i}}}{\displaystyle\frac{1}{1+e^{\beta_0+\beta_1x_i}}}= e^{\beta_0+\beta_1x_i}.$$
Assim, ao tomarmos dois valores distintos da variável explicativa, $x_j$ e $x_{j+1}$, obtemos
$$OR=\frac{g(x_{j+1})}{g(x_{j})}=\frac{e^{\beta_0+\beta_1~x_{j+1}}}{e^{\beta_0+\beta_1~x_{j}}}. \tag{4.1.2.1.1}$$
Temos ainda que:
$$\ln(OR)=\ln\left[\frac{g(x_{j+1})}{g(x_{j+1})} \right]=\ln\left[g(x_{j+1}) \right]-\ln\left[g(x_{j}) \right]$$
$$=\beta_0 + \beta_1 x_{j+1} -\beta_0 -\beta_1 x_{j} = \beta_1(x_{j+1}-x_j).$$
Fazendo $x_{j+1}-x_j=1$ unidade, então $$\ln(OR) = \ln(e^{{\beta_1}})={\beta_1}. $$
Assim, temos o quão provável o resultado ocorrerá entre os indivíduos $x_{j+1}$ em relação aos indivíduos $x_j$, fazendo, portanto, algumas análises:
$$\beta_1> 0~~\Rightarrow~~ OR> 1 \Rightarrow ~~ \pi(x_{j+1})~~>~~\pi(x_{j})$$
$$\beta_1< 0~~\Rightarrow~~ OR< 1 \Rightarrow ~~ \pi(x_{j+1})~~ <~~\pi(x_{j})$$
Veja “variáveis independentes categóricas” quando a variável explicativa é categórica.
Exemplo 4.1.2.1.1
Considerando os dados e a estimativa de $\beta_1$ do Exemplo 4.2.1.1, vamos calcular o Odds Ratio. Assim, temos:
$$\text{OR}~(\hat{\beta_1}) = e^{-0,0663} = 0,936.$$
Como o Odds Ratio é menor que 1, a probabilidade de cometer Erros de Montagem tende a diminuir quando aumentam as Horas de Treinamento.
4.1.2.2 Estimativa dos desvios padrão
No modelo de regressão logístico o desvio padrão dos estimadores é obtido a partir da matriz de informação de Fisher. Podemos ainda obter a matriz de informação de Fisher $I(\widehat{\beta})$ para o modelo logístico a partir dos dados, da seguinte forma, $I(\widehat{\beta}) = X^{^\prime } V X,$ sendo que
$$X=\begin{bmatrix}1~~x_{1} \cr 1~~x_{2} \cr \vdots \quad \vdots \cr 1~~x_{n} \cr \end{bmatrix}, {V} = \text{diag}~[m_1 \hat{\pi}_1 (1 - \hat{\pi}_1),\ldots, m_n\hat{\pi}_n (1 - \hat{\pi}_n)],$$
$m_i$ é o número de repetições para cada elemento da amostra, $i=1,\ldots,n.$
As variâncias e covariâncias dos estimadores $\widehat{\beta}=(\widehat{\beta_0}, \widehat{\beta_1})$ são obtidos, invertendo a matriz de informação de Fisher, isto é, calculando $\widehat{\Sigma} = I^{-1}(\widehat{\beta}).$
O $j$-ésimo elemento da diagonal principal da matriz $\widehat{\Sigma}$ é a variância do estimador $\widehat{\beta}_j,$ denominada $\widehat{\sigma}^2(\widehat{\beta}_j).$ Os demais elementos da matriz $\widehat{\Sigma}$ são as covariâncias entre $(\widehat{\beta}_j;\widehat{\beta}_u),$ $j\neq u.$
Desta forma o desvio padrão é definido como: $$\widehat{DP}(\widehat{\beta}_j)= \sqrt{\widehat{\sigma}^2(\widehat{\beta}_j)}.$$
Exemplo 4.1.2.2.1
Vamos considerar os dados e as estimativas do Exemplo 4.2.1.1.
Para o cálculo dos desvios padrão, primeiramente é preciso obter a matriz $\widehat{\Sigma}=I(\widehat{\beta})^{-1}=(X^\prime V X)^{-1}$
$$X^{^\prime } = \begin{bmatrix}1 \quad 1 \quad 1 \quad 1 \quad 1 \ldots 1 \cr 30 \quad 30 \quad 30 \quad 30 \quad 29\ldots 5 \cr \end{bmatrix}_{(2 \times 30)} \quad \text{e}$$
$${V} = \hbox{diag}~ [200 * 0,01806045 * 0,9819395 ; \ldots ; 200 * 0,08801434 * 0,9119857]$$
$$=\begin{bmatrix}3,546855 \quad 0 \quad \ldots \quad 0 \cr 0 \quad 3,546855 \quad \ldots \quad 0 \cr \vdots \qquad \vdots \qquad \ddots \quad \qquad \vdots \cr 0 \quad 0 \quad \ldots \quad 16,053563 \cr \end{bmatrix}_{(30 \times 30)}$$
portanto a matriz de informação de Fisher será:
$$I(\widehat{\beta})=\begin{bmatrix}244,5157 \quad 3641,87 \cr 3641,8697 \quad 66337,09 \cr \end{bmatrix}$$
e a matriz de variâncias e covariâncias:
$$\Sigma=\begin{bmatrix}0,022432 \quad -0,001232 \cr -0,001232 \quad 0,000083 \cr \end{bmatrix}$$
Da matriz $\Sigma$ concluímos que as estimativas dos desvios padrão dos estimadores $\widehat{\beta}=(\widehat{\beta_0},\widehat{\beta_1}),$ são:
$$\widehat{{DP}}(\widehat{\beta_0}) = \sqrt{\widehat{\sigma}^2(\widehat{\beta_0})}= \sqrt{0,022432}=0,1498$$
$$\widehat{{DP}}(\widehat{\beta_1})=\sqrt{\widehat{\sigma}^2(\widehat{\beta_1})} = \sqrt{0,000083} =0,0091.$$
4.1.3 Inferência em um modelo logístico simples
Após estimar os coeficientes, temos interesse em assegurar a significância das variáveis no modelo. Isto geralmente envolve formulação e teste de uma hipótese estatística para determinar se a variável independente no modelo é significativamente relacionada com a variável resposta. Para isso, temos os testes de hipóteses. Os testes de hipóteses mais utilizados são os testes da Razão da Verossimilhança, Wald e Escore. A seguir, temos a abordagem de cada um deles.
4.1.3.1 Teste de Wald
O teste de Wald é obtido por comparação entre a estimativa de máxima verossimilhança do parâmetro ($\widehat{\beta_1}$) e a estimativa de seu erro padrão. A razão resultante, sob a hipótese $\text{H}_0:~\beta_1=0,$ tem distribuição normal padrão.
A estatística do teste Wald para a regressão logística é $W_j=\displaystyle\frac{\widehat{\beta_1}}{\widehat{DP}(\widehat{\beta_1})}.$
O p-valor é definido como
$$\text{p-valor}=\mathbb{P}[|Z|>|W_j||H_0]=2\mathbb{P}[Z>|W_j||H_0],$$
sendo que Z denota a variável aleatória da distribuição normal padrão.
Hauck e Donner (1977) examinaram o desempenho do teste de Wald e descobriram que ele se comporta de maneira estranha, em determinadas situações; frequentemente não rejeitando a hipótese nula quando o coeficiente é significativo. Eles recomendam a utilização do teste da razão de verossimilhança para testar se realmente o coeficiente não é significativo quando o teste de Wald não rejeita a hipótese nula.
Exemplo 4.1.3.1.1
Vamos agora testar se os parâmetros do Exemplo 4.2.1.1 é significativo para o modelo. Para isso, precisamos dos valores dos desvios padrão calculados no Exemplo 4.1.2.2.1.
Os valores da estatística do teste de Wald, para as hipóteses ${H}_0:\beta_j=0$ e ${H}_1:\beta_j \neq 0,$ $j=0,1,$ são:
$$W_0=\frac{\widehat{\beta_0}}{\widehat{DP}(\widehat{\beta_0})}=\frac{-2,0065}{0,149772}= -13,397$$
$$W_1=\frac{\widehat{\beta_1}}{\widehat{DP}(\widehat{\beta_1})}=\frac{-0,0663}{0,009092}= -7,291$$
Para estas hipóteses, os valores de p são:
Para $\widehat{\beta_0} = P(|Z|> 13,396) = 0,000.$
Para $\widehat{\beta_1} = P(|Z|> 7,29) = 0,000.$
Como o p-valor é menor que o nível de significância $\alpha = 0,05$ em ambos os casos, concluímos que os parâmetros $\beta_0$ e $\beta_1$ são significativos no modelo.
4.1.3.2 Teste da Razão de Verossimilhança
Na regressão linear o interesse está no valor da SQR. Um valor alto da SQR sugere que a variável independente é importante, caso contrário, a variável independente não é útil na predição da variável resposta.
Na regressão logística a ideia é a mesma: comparar os valores observados da variável resposta com os valores preditos obtidos dos modelos com e sem a variável em questão. A comparação dos observados com os valores preditos é baseado no log da verossimilhança. Para entender melhor essa comparação, é útil pensar em um valor observado da variável resposta também como sendo um valor predito resultante de um modelo saturado. Um modelo saturado é aquele que contém tantos parâmetros quanto observações.
A comparação dos observados com os valores preditos usando a função de verossimilhança é baseada na seguinte expressão:
$$D=-2ln\frac{\hbox{(verossimilhança do modelo ajustado)}}{\hbox{(verossimilhança do modelo saturado)}}.$$
Com o propósito de assegurar a significância de uma variável independente, comparamos o valor da D com e sem a variável na equação. A mudança em D devido a inclusão da variável no modelo é obtida da seguinte maneira: $$G=D \hbox{(modelo~sem~a~variável)} - D \hbox{(modelo~com~a~variável)}.$$
Podemos então escrever a estatística G como:
$$G=-2ln\frac{\hbox{(verossimilhança sem a variável)}}{\hbox{(verossimilhança com a variável)}}.$$
ou ainda: $$G=-2ln(L_s)+2ln(L_c),$$
em que $L_s$ é a verossimilhança do modelo sem a covariável e $L_c$ é a verossimilhança do modelo com a covariável.
Queremos testar:
$$\begin{cases} H_{0}:\beta_1=0 \cr H_1:\beta_1\neq 0 \end{cases}$$
Sob a hipótese nula, a estatística $G$ tem distribuição chi-quadrado com 1 grau de liberdade.
Exemplo 4.1.3.2.1
Vamos considerar o Exemplo 4.1.2.1 para verificar se a variável “horas de treinamento” é significativa para explicar o erro na montagem, através do teste da razão de verossimilhança (TRV).
O valor do log da verossimilhança do modelo apenas com o intercepto ($L_s$) é -1064,183 e do modelo com a variável ($L_c$) é -1035,089.
Assim, o valor da estatística teste é:
$$G=-2(-1064,183)-(-2(-1035,089))= 58,188.$$
O p-valor $P(\chi^{2}_{1}>\ G=58,188)< 0,0001$.
Rejeitamos a hipótese nula. Assim, pelo TRV, temos que a variável horas de treinamento é significativa para o modelo
4.1.3.3 Teste Score
A estatística teste para o Teste Score é:
$$ST=\frac{\displaystyle\sum_{i=1}^{n}x_i(y_i-\bar{y})}{(\overline{y}(1-\overline{y})\displaystyle\sum_{i=1}^{n}(x_i-\bar{x})^2)^{1/2}},$$
em que $\bar{y}=\hat{\pi}$ (proporção de sucessos na amostra).
No Teste Score também temos o interesse em testar:
$$\begin{cases} H_{0}:\beta_1=0 \cr H_1:\beta_1\neq 0 \end{cases} $$
O p-valor é definido como $P(|Z|> |ST|)$, sendo que Z denota a variável aleatória da distribuição normal padrão.
4.1.4 Intervalos de Confiança
4.1.4.1 Intervalo de Confiança para os parâmetros
A base da construção das estimativas do intervalo de confiança para os parâmetros é a mesma teoria estatística que usamos para os testes de significância do modelo. Em particular, um intervalo de confiança para a inclinação e intercepto são baseados em seus respectivos testes de Wald. O intervalo de confiança de $100(1-\alpha)$% para o parâmetro $\beta_1$ é:
$$IC(\beta_1,1-\alpha)= [\hat{\beta_1}-z_{1-\alpha/2}DP(\hat{\beta_1}); \quad \hat{\beta_1}+z_{1-\alpha/2}DP(\hat{\beta_1})].$$
E para o intercepto:
$$IC(\beta_0,1-\alpha)= [\hat{\beta_0}-z_{1-\alpha/2}DP(\hat{\beta_0}); \quad \hat{\beta_0}+z_{1-\alpha/2}DP(\hat{\beta_0})],$$
em que $z_{1-\alpha/2}$ é o ponto da normal padrão correspondente a $100(1-\alpha/2)$%.
4.1.4.2 Intervalo de Confiança para Logito
A logito é a parte linear do modelo de regressão logística. O estimador para logito é: $$\hat{g}(x)=\hat{\beta_0}+\hat{\beta_1}x.$$
O estimador da variância do estimador da logito requer a obtenção da variância da soma. No caso é:
$$\hat{Var}[\hat{g}(x)]=\hat{Var}(\hat{\beta_0})+x^2\hat{Var}(\hat{\beta_1})+2x\hat{Cov}(\hat{\beta_0},\hat{\beta_1}). \tag{4.1.4.2.1}$$
O intervalo de confiança para a logito é: $$IC(g(x),1-\alpha)= [\hat{g}(x)-z_{1-\alpha/2}DP(\hat{g}(x)); \quad \hat{g}(x)+z_{1-\alpha/2}DP(\hat{g}(x))],$$
em que $DP(\hat{g}(x))$ é a raiz quadrada de 4.1.4.2.1 e $z_{1-\alpha/2}$ é o ponto da normal padrão.
4.1.4.3 Intervalo de Confiança para os valores ajustados
O estimador do logito e seu intervalo de confiança fornece o estimador dos valores ajustados. O intervalo de confiança dos valores ajustados é dado por:
$$IC(\pi,1-\alpha)= \left[\dfrac{e^{\hat{g}(x)-z_{1-\alpha/2}DP(\hat{g}(x))}}{1+e^{\hat{g}(x)-z_{1-\alpha/2}DP(\hat{g}(x))}}; \quad \dfrac{e^{\hat{g}(x)+z_{1-\alpha/2}DP(\hat{g}(x))}}{1+e^{\hat{g}(x)+z_{1-\alpha/2}DP(\hat{g}(x))}}\right]. \tag{4.1.4.2.2}$$
4.1.4.4 Intervalo de Confiança para a Odds Ratio
Sejam os limites do intervalo de confiança para $\beta_1$:
$\beta_{I}=\hat{\beta_1}-z_{1-\alpha/2}DP(\hat{\beta_1})$ e $\beta_{S}=\hat{\beta_1}+z_{1-\alpha/2}DP(\hat{\beta_1}).$
O intervalo de confiança para a Odds Ratio é:
$$IC(\hbox{Odds Ratio},1-\alpha)= [e^{\beta_{I}}; \quad ~e^{\beta_{S}}]. \tag{4.1.4.2.3}$$
Exemplo 4.1.4.1
Vamos construir intervalo de confiança de Wald para $\beta_0$ e $\beta_1$ considerando as estimativas dos parâmetros e dos seus desvios padrão já calculados nos exemplos 4.1.2.1 e 4.1.2.2.1.
IC para $\beta_0$
$$IC(\beta_0,0,95) = [-2,0066-1,96 * 0,1498; \quad -2,0066+1,96 * 0,1498],$$
$$IC(\beta_0, 0,95)= [-2,0066-0,2936; \quad -2,0066+0,2936],$$
$$IC(\beta_0,0,95)= [-2,3002; \quad -1,713].$$
IC para $\beta_1$
$$IC(\beta_1,0,95)= [-0,0663-1,96 * 0,0091; \quad -0,0663+1,96 * 0,0091].$$
$$IC(\beta_1,0,95)= [-0,0663-0,01783; \quad -0,0663+0,01783].$$
$$IC(\beta_1,0,95)= [-0,08413; \quad -0,04847].$$
Exemplo 4.1.4.2
Vamos construir intervalo de 95% de confiança para a logito.
Consideramos a primeira observação, com $X$ (horas de treinamento) $=30$. O estimador para a logito é:
$$\hat{g}(x=30)=-2,0066-0,0663 * 30=-3,9956.$$
O estimador da variância é:
$$\hat{Var}[\hat{g}(x=30)]=0,022432+900 * 0,000083+60 * (-0,001232)$$
$$\hat{Var}[\hat{g}(x=30)]=0,022432+0,0747-0,07392$$
$$\hat{Var}[\hat{g}(x=30)]=0,022432+0,0747-0,07392$$
$$\hat{Var}[\hat{g}(x=30)]=0,0232$$
Assim, o intervalo de 95% de confiança para logito para $x=30$ é:
$$IC(g(x=30),0,95)= [-3,9956-1,96 * 0,1523; \quad -3,9956+1,96 * 0,1523],$$
$$IC(g(x=30),0,95)= [-3,9956-0,2985; \quad -3,9956+0,2985],$$
$$IC(g(x=30),0,95)= [-4,2941; \quad -3,6971],$$
Fazendo analogamente para cada valor da variável explicativa, temos na Tabela 4.4.3 os intervalos de confiança para a logito.
| X | Limite Inferior | Limite Superior |
|---|---|---|
| 30 | -4,294 | -3,697 |
| 30 | -4,294 | -3,697 |
| 30 | -4,294 | -3,697 |
| 29 | -4,212 | -3,647 |
| 28 | -4,130 | -3,596 |
| 27 | -4,048 | -3,546 |
| 26 | -3,966 | -3,495 |
| 26 | -3,966 | -3,495 |
| 25 | -3,885 | -3,443 |
| 24 | -3,804 | -3,391 |
| 23 | -3,724 | -3,339 |
| 20 | -3,489 | -3,176 |
| 20 | -3,489 | -3,176 |
| 20 | -3,489 | -3,176 |
| 17 | -3,266 | -3,002 |
| 17 | -3,266 | -3,002 |
| 17 | -3,266 | -3,002 |
| 16 | -3,195 | -2,940 |
| 15 | -3,127 | -2,875 |
| 13 | -2,999 | -2,738 |
| 12 | -2,938 | -2,666 |
| 11 | -2,880 | -2,592 |
| 11 | -2,880 | -2,592 |
| 11 | -2,880 | -2,592 |
| 10 | -2,823 | -2,517 |
| 10 | -2,823 | -2,517 |
| 9 | -2,767 | -2,440 |
| 8 | -2,713 | -2,361 |
| 8 | -2,713 | -2,361 |
| 5 | -2,554 | -2,122 |
Tabela 4.4.3: Intervalo de confiança para a logito
Exemplo 4.1.4.3
Sabemos pelo exemplo 4.1.4.2 que o intervalo de 95% de confiança para $g(30)$ é $[-4,294; -3,697]$. Assim, através da equação 4.1.4.2.2, temos o seguinte intervalo de confiança para o valor ajustado:
$$IC(\pi(30),0,95)= \left[\frac{e^{-4,294}}{1+e^{-4,294}}; \quad \frac{e^{-3,697}}{1+e^{-3,697}}\right]$$
$$IC(\pi(30),0,95)= [0,01346; \quad 0,02420]$$
Fazendo analogamente para cada valor da variável explicativa, temos na Tabela 4.4.4 os intervalos de confiança para os valores ajustados.
| X | Limite Inferior | Limite Superior |
|---|---|---|
| 30 | 0,01346 | 0,02420 |
| 30 | 0,01346 | 0,02420 |
| 30 | 0,01346 | 0,02420 |
| 29 | 0,01460 | 0,02541 |
| 28 | 0,01583 | 0,02669 |
| 27 | 0,01716 | 0,02804 |
| 26 | 0,01859 | 0,02947 |
| 26 | 0,01859 | 0,02947 |
| 25 | 0,02013 | 0,03097 |
| 24 | 0,02179 | 0,03257 |
| 23 | 0,02356 | 0,03426 |
| 20 | 0,02963 | 0,04006 |
| 20 | 0,02963 | 0,04006 |
| 20 | 0,02963 | 0,04006 |
| 17 | 0,03677 | 0,04735 |
| 17 | 0,03677 | 0,04735 |
| 17 | 0,03677 | 0,04735 |
| 16 | 0,03934 | 0,05023 |
| 15 | 0,04199 | 0,05341 |
| 13 | 0,04748 | 0,06076 |
| 12 | 0,05030 | 0,06500 |
| 11 | 0,05317 | 0,06964 |
| 11 | 0,05317 | 0,06964 |
| 11 | 0,05317 | 0,06964 |
| 10 | 0,05612 | 0,07470 |
| 10 | 0,05612 | 0,07470 |
| 9 | 0,05913 | 0,08021 |
| 8 | 0,06223 | 0,08617 |
| 8 | 0,06223 | 0,08617 |
| 5 | 0,07213 | 0,10700 |
Tabela 4.4.4: Intervalo de Confiança para os valores ajustados
Exemplo 4.1.4.4
O intervalo de 95% de confiança para o Odds Ratio da variável explicativa (Horas de Treinamento), dado pela equação 4.1.4.2.3, utilizando os limites do intervalo de confiança para $\beta_1$ do exemplo 4.1.4.1 é:
$$IC(\hbox{Odds Ratio};0,95)= [\exp(-0,08413); \exp(-0,04847)]$$
$$IC(\hbox{Odds Ratio};0,95)= [0,919; 0,953].$$
4.2 Regressão Logística Múltipla
Apresentamos na Seção 4.1 o modelo de regressão logística considerando apenas uma variável explicativa. Assim como no modelo de regressão linear, podemos ajustar um modelo para a variável resposta levando em conta mais de uma variável explicativa (covariável), o que chamamos de Modelo de Regressão Logística Múltipla.
Assim, um modelo de regressão logística múltipla é usado para o caso de regressão com mais de uma variável explicativa.
Considere um conjunto de $p$ variáveis independentes denotadas como um vetor ${X}=(X_1,X_2,X_3,\ldots,X_p)$. Neste caso a função de ligação é
$$g(X)=\ln\left(\frac{\pi(X)}{1-\pi(X)}\right) = \beta_0+\beta_1 X_1+\beta_2 X_2 + \ldots + \beta_p X_p \tag{4.2.1}$$
e
$$E[Y]=\pi(X)=\frac{e^{g(X)}}{1+e^{g(X)}}. \tag{4.2.2}$$
4.2.1 Estimativas dos parâmetros do modelo
Para obter as estimativas para vetor $\beta=(\beta_0,\beta_1,\ldots,\beta_p)$ dos parâmetros do modelo e a matriz de covariâncias será utilizado o método de máxima verossimilhança, da seguinte forma:
$$L~(\beta_0,\beta_1,\ldots,\beta_p|(x_i;m_i;y_i))=\sum^n_{i=1} \left[ y_i~ g(X) -m_i \ln (1+e^{g(X)})\right], \tag{4.2.1.1}$$
em que g(X) é dado por 4.2.1.
Derivando (4.2.1.1) em relação aos parâmetros, temos:
$$\frac{\partial}{\partial\beta_0} L(\beta_0,\beta_1,\ldots,\beta_p) = \sum^n_{i=1} y_i - \sum^n_{i=1} m_i \frac{e^{g(X)}}{1 + e^{g(X)}}.$$
$$\frac{\partial}{\partial\beta_j} ~L(\beta_0,\beta_1,\ldots,\beta_p)=\sum^n_{i=1} y_i~x_j - \sum^n_{i=1} m_i x_j \frac{e^{g(X)}}{1 + e^{g(X)}}.$$
Igualando a zero e substituindo os parâmetros pelos estimadores tem-se as seguintes equações:
$$\sum^n_{i=1}~y_i~(1+ e^{g(X)}) - \sum^n_{i=1}m_i~e^{g(X)} = 0$$
$$\sum^n_{i=1}~y_i~x_i~(1+ e^{g(X)})-\sum^n_{i=1}m_i~x_i~e^{g(X)} = 0$$
A solução dessas equações fornece as estimativas dos parâmetros do modelo, utilizando processos iterativos, análoga à estimação dos parâmetros do modelo de regressão logística simples.
Após obter as estimativas dos parâmetros do modelo podemos calcular as probabilidades ajustadas:
$$\widehat{\pi}_i = \frac{e^{\hat{g}(X_i)}}{1 + e^{\hat{g}(X_i)}},$$
em que
$$\hat{g}(X_i)=\hat{\beta_0}+\hat{\beta_1} X_{1i} + \ldots + \hat{\beta_p} X_{pi}.$$
Exemplo 4.2.1.1:
Em uma empresa são fabricadas peças de ferro que são moldadas em moldes de areia. Entre as peças produzidas as que apresentam uma grande quantidade de areia incrustada são consideradas Refugo. A Volatilidade da Areia e o coeficiente RFV (Resistência ao Fluido Verde) influenciam na quantidade de areia incrustada. A partir dos dados da Tabela (4.4.5) avaliar a significância das variáveis e interação entre elas.
| Observação | Quantidade Produzida | Refugo | Volatilidade | RFV |
|---|---|---|---|---|
| 1 | 832 | 270 | 1,906 | 5,642 |
| 2 | 996 | 152 | 1,766 | 7,63 |
| 3 | 1224 | 289 | 1,673 | 5,253 |
| 4 | 712 | 2 | 1,982 | 5,223 |
| 5 | 2072 | 11 | 2 | 5,064 |
| 6 | 544 | 14 | 2,12 | 5,395 |
| 7 | 700 | 5 | 2,085 | 6,138 |
| 8 | 3840 | 47 | 1,97 | 5,82 |
| 9 | 1940 | 101 | 2,15 | 4,498 |
| 10 | 1005 | 17 | 2,37 | 6,478 |
| 11 | 1260 | 26 | 2,37 | 5,826 |
| 12 | 1815 | 308 | 2,597 | 6,052 |
| 13 | 1340 | 79 | 2,44 | 5,839 |
| 14 | 1485 | 134 | 2,473 | 5,08 |
| 15 | 1585 | 127 | 2,493 | 5,313 |
| 16 | 1095 | 83 | 2,43 | 5,21 |
| 17 | 1370 | 81 | 3,42 | 5,04 |
| 18 | 1405 | 58 | 3,607 | 5,2222 |
Tabela 4.4.5: Dados do molde de areia.
O seguinte modelo é proposto para o exemplo
$$\hbox{Probabilidade de refugo}=\pi_i=\dfrac{e^{\beta_0 + \beta_1 x_1 +\beta_2 x_2+\beta_3 x^2_1 + \beta_4 x^2_2 + \beta_5 x_1x_2}}{1+e^{\beta_0 + \beta_1 x_1 +\beta_2 x_2+\beta_3 x^2_1 + \beta_4 x^2_2 + \beta_5 x_1x_2}}$$
As estimativas dos parâmetros do modelo são:
$\widehat{\beta_0}$ = 30,3389 (Intercepto)
$\widehat{\beta_1}$ = -14,1496 (Volatilidade)
$\widehat{\beta}_2$ = -6,26927 (RFV)
$\widehat{\beta}_3$ = 1,101 (Efeito quadrático da Volatilidade)
$\widehat{\beta}_4$ = 0,2819 (Efeito quadrático do RFV)
$\widehat{\beta}_5$ = 1,5376 (Interação entre a Volatilidade e o RFV)
E as estimativas do odds ratio, como mostradas no caso do modelo de regressão logística simples, são:
$OR(\widehat{\beta_1})$ = $\exp(-14,1496) = 0,00$
$OR(\widehat{\beta}_2)$ = $\exp(-6,26927) = 0,00$
$OR(\widehat{\beta}_3)$ = $\exp(1,101) = 3,01$
$OR(\widehat{\beta}_4)$ = $\exp(0,2819) = 1,33$
$OR(\widehat{\beta}_5)$ = $\exp(1,5376) = 4,65$
O vetor de probabilidade ajustada é
$$\widehat{\pi} = (0,076; 0,128; 0,131; 0,075; 0,076; 0,061; 0,072; 0,071; 0,068; 0,101; 0,061; 0,081;$$
$ 0,062; 0,040; 0,044; 0,044; 0,043; 0,079).$
4.2.1.1 Estimativa do desvio padrão
O método de estimar as variâncias e covariâncias dos coeficientes estimados segue a teoria da estimação de máxima verossimilhança. Essa teoria assegura que os estimadores são obtidos da matrix de derivadas segunda parciais da função log de verossimilhança. Essas derivadas parciais tem a seguinte forma geral:
$$\frac{\partial^2L(\beta)}{\partial\beta_j^2}=-\sum_{i=1}^nm_ix_{ij}^{2}\pi_i(1-\pi_i). \quad \quad \tag{4.2.1.1.1}$$
$$\frac{\partial^2L(\beta)}{\partial\beta_j\partial\beta_l}=-\sum_{i=1}^nm_ix_{ij}x_{il}\pi_i(1-\pi_i), \quad \tag{ 4.2.1.1.2}$$
para j,l=0,1,…,p em que $\pi_i$ simplifica $\pi(x_i)$. Seja a matriz $(p+1)\times(p+1)$ que contém os termos negativos de (4.2.1.1.1) e (4.2.1.1.2) denotada por $I(\beta)$ que é matriz de informação de Fisher. As variâncias e covariâncias dos coeficientes estimados são obtidos da inversa da matriz: $Var(\beta)=I^{-1}(\beta)$. A variância de $\beta_j$, $Var(\beta_j)$, é o $j^{th}$ elemento da diagonal da matriz e a $Cov(\beta_j,\beta_l)$ é obtida através do elemento da matriz referente a linha de $\beta_j$ e coluna de $\beta_l$ ou vice-versa, já que $Cov(\beta_j,\beta_l)=Cov(\beta_l,\beta_j).$ Os estimadores das variâncias e covariâncias, $\hat{Var}(\hat{\beta})$, são obtidos de $Var(\beta)$ em $\hat{\beta}$.
Ainda, a matriz de informação de Fisher estimada pode ser obtida por: $\hat{I}(\hat{\beta})=X^\prime VX$.
$$X=\begin{bmatrix}1 \quad x_{11} \quad \ldots \quad x_{1p} \cr 1 \quad x_{21} \quad \ldots \quad x_{2p} \cr \vdots \cr 1 \quad x_{n1} \quad \ldots \quad x_{np} \cr \end{bmatrix}_{n \times (p+1)} $$
e
$$V= \begin{bmatrix} m_1\widehat{\pi}_1(1-\widehat{\pi}_1) \quad 0 \quad ~~\ldots \quad 0 \cr 0 \quad ~~ m_2\widehat{\pi}_2(1-\widehat{\pi}_2) \quad ~~\ldots \quad 0 \cr \vdots \quad ~~\vdots \quad ~~\ddots \quad ~~\vdots \cr 0 \quad 0 \quad ~~\ldots \quad ~~m_n\widehat{\pi}_n(1-\widehat{\pi}_n) \cr \end{bmatrix}_{n \times n}$$
Dessa forma, o desvio padrão do coeficiente $\beta_j$ é: $$\hat{DP}(\hat{\beta}_j)=\sqrt{\hat{Var}(\hat{\beta_j})}. \tag{4.2.1.1.3}$$
Exemplo 4.2.1.1.1
Para os dados do Exemplo 4.2.1.1, vamos calcular a estimativa do desvio padrão de cada parâmetro considerado no exemplo em questão.
O primeiro passo é obter a matriz de informação de Fisher, conforme 4.1.2.3,
$$X=\begin{bmatrix}1 \quad 1,91 \quad ~5,64 \quad ~3,63 \quad ~31,83 \quad 10,75 \cr 1 \quad 1,77 \quad ~7,63 \quad ~3,12 \quad ~58,22 \quad 13,47 \cr 1 \quad 1,67 \quad ~5,25 \quad ~2,80 \quad ~27,59 \quad ~8,79 \cr 1 \quad 1,98 \quad ~5,22 \quad ~3,93 \quad ~27,28 \quad 10,35 \cr 1 \quad ~2,00 \quad ~5,06 \quad ~4,00 \quad ~25,64 \quad 10,13 \cr 1 \quad ~2,12 \quad ~5,39 \quad ~4,49 \quad ~29,11 \quad 11,44 \cr 1 \quad ~2,08 \quad ~6,14 \quad ~4,35 \quad ~37,68 \quad 12,80 \cr 1 \quad 1,97 \quad ~5,82 \quad ~3,88 \quad ~33,87 \quad 11,47 \cr 1 \quad ~2,15 \quad ~4,50 \quad ~4,62 \quad ~20,23 \quad ~9,67 \cr 1 \quad ~2,37 \quad ~6,48 \quad ~5,62 \quad ~41,96 \quad 15,35 \cr 1 \quad ~2,37 \quad ~5,83 \quad ~5,62 \quad ~33,94 \quad 13,81 \cr 1 \quad ~2,60 \quad ~6,05 \quad ~6,74 \quad ~36,63 \quad 15,72 \cr 1 \quad ~2,44 \quad ~5,84 \quad ~5,95 \quad ~34,09 \quad 14,25 \cr 1 \quad ~2,47 \quad ~5,08 \quad ~6,12 \quad ~25,81 \quad 12,56 \cr 1 \quad ~2,49 \quad ~5,31 \quad ~6,22 \quad ~28,23 \quad 13,25 \cr 1 \quad ~2,43 \quad ~5,21 \quad ~5,90 \quad ~27,14 \quad 12,66 \cr 1 \quad ~3,42 \quad ~5,04 \quad 11,70 \quad ~25,40 \quad 17,24 \cr 1 \quad ~3,61 \quad ~5,22 \quad 13,01 \quad ~27,27 \quad 18,84 \cr \end{bmatrix}_{18 \times 6} $$
$$V=\begin{bmatrix}832\times 0,079518 \times (1 - 0,079518) \quad \quad 0 \quad ~\ldots \quad 0 \cr 0 \quad ~996\times 0,118598 \times (1-0,118598) \quad ~\ldots \quad 0 \cr \vdots \quad \quad \vdots \quad ~\ddots \cr 0 \quad 0 \quad ~\ldots \quad 1405\times 0,047904 \times (1-0,047904) \cr \end{bmatrix}_{18 \times 18}$$
ou
$$V= \begin{bmatrix} 58,71 \quad 0 \quad ~\ldots \quad 0 \cr 0 \quad 111,16 \quad ~\ldots \quad 0 \cr \vdots \quad ~\vdots \quad ~\ddots \quad ~\vdots \cr 0 \quad 0 \quad ~\ldots \quad 101,81 \cr \end{bmatrix}_{18 \times 18}$$
$$\Sigma(\beta) = (X' V X)^{-1} =\left[ \begin{array}{ccc}13,47 \quad ~~-4,57 \quad ~~-3,05 \quad ~~ 0,23 \quad ~~ 0,15 \quad ~~ 0,63 \cr -4,57 \quad ~1,84 \quad ~~ 0,93 \quad ~~-0,11 \quad ~~-0,04 \quad ~~-0,23 \cr -3,05 \quad ~~ 0,93 \quad 0,73 \quad ~~-0,04 \quad ~~-0,04 \quad ~~-0,14 \cr 0,23 \quad ~~-0,11 \quad ~~-0,04 \quad 0,01 \quad ~~ 0,00 \quad ~~ 0,01 \cr 0,15 \quad ~~-0,04 \quad ~~-0,04 \quad ~~ 0,00 \quad 0,00 \quad ~~ 0,01 \cr 0,63 \quad ~~-0,23 \quad ~~-0,14 \quad ~~ 0,01 \quad ~~ 0,01 \quad 0,03 \cr \end{array} \right]$$
Para $\beta_0=\sqrt{\sigma^2(\widehat{\beta_0})}=\sqrt{13,46607} =3,67$
Para $\beta_1=\sqrt{\sigma^2(\widehat{\beta_1})}=\sqrt{1,836394}=1,36$
Para $\beta_2=\sqrt{\sigma^2(\widehat{\beta}_2)}=\sqrt{0,731351}= 0,86$
Para $\beta_3=\sqrt{\sigma^2(\widehat{\beta}_3)} = \sqrt{0,0104028}=0,10$
Para $\beta_4=\sqrt{ \sigma^2(\widehat{\beta}_4)}=\sqrt{0,0021488}=0,05$
Para $\beta_5=\sqrt{\sigma^2(\widehat{\beta}_5)}=\sqrt{0,0331131}=0,18$
4.2.2 Inferência no modelo de regressão logística múltipla
Assim como na Inferência em um modelo de regressão logística simples, podemos testar a significância dos parâmetros pelo teste de Wald, score e Razão de Verossimilhança (TRV).
4.2.2.1 Teste da Razão de Verossimilhança
O teste da razão de verossimilhança para a significância dos p coeficientes das variáveis independentes do modelo é realizado da mesma maneira que no modelo de regressão logística simples. A estatística teste G é dada por:
$$G=-2ln\left[\frac{\hbox{(verossimilhança sem a variável)}}{\hbox{(verossimilhança com a variável)}}\right].$$
ou ainda:
$$G=-2ln(L_s)+2ln(L_c),$$
em que $L_s$ é a verossimilhança do modelo sem a covariável e $L_c$ é a verossimilhança do modelo com a covariável.
No caso da regressão múltipla, temos o interesse em saber se pelo menos uma variável é significativa para o modelo. Sob a hipótese nula, os p coeficientes são iguais a zero , assim, a estatística G tem distribuição Qui-Quadrado com p graus de liberdade. Nesse caso $L_c$ é a verossimilhança do modelo com as p variáveis explicativas e $L_s$ é a verossimilhança do modelo apenas com o intercepto.
4.2.2.2 Teste de Wald
Vamos considerar a seguinte hipótese:
$$\begin{cases} H_{0}:\beta_j=0 \cr H_1:\beta_j\neq 0 \end{cases} $$
Para testar a hipótese acima, a estatística de Wald é obtida da seguinte forma:
$$W_j=\frac{\widehat{\beta}_j}{\widehat{DP}(\widehat{\beta}_j)}. \tag{4.2.2.2.1}$$
onde o desvio padrão é obtido de (4.2.1.1.3). Se não rejeitarmos $H_0$ temos que a variável $X_j$ não explica a variável resposta.
De forma equivalente, teste Wald também pode ser obtido pela multiplicação dos seguintes vetores:
$$W=\widehat{\beta}^\prime [\widehat{I}(\widehat{\beta})]^{-1}\widehat{\beta}=\widehat{\beta}^\prime (X^\prime V X)\widehat{\beta}.$$
Em que $\hat{I}(\beta)$ é a matriz da informação de Fisher estimada, apresentada na estimativa do desvio padrão do modelo de regressão logística simples.
Exemplo 4.2.2.1:
Para os dados do Exemplo 4.2.1.1 vamos testar a significância dos parâmetros pelo teste da Razão de Verossimilhança e Wald.
Lembrando que o modelo proposto é:
$$\hbox{Probabilidade de refugo}=\pi_i=\dfrac{e^{\beta_0 + \beta_1 x_1 +\beta_2 x_2+\beta_3 x^2_1 + \beta_4 x^2_2 + \beta_5 x_1x_2}}{1+e^{\beta_0 + \beta_1 x_1 +\beta_2 x_2+\beta_3 x^2_1 + \beta_4 x^2_2 + \beta_5 x_1x_2}},$$
em que $x_1$ é a variável Volatilidade e $x_2$ é a variável RFV.
- TRV (teste da Razão de Verossimilhança)
Pelo TRV, vamos testar primeiramente se pelo menos uma variável é significativa para o modelo. Para isso, precisamos do log da verossimilhança dos modelos com e sem as variáveis em consideração.
O log da verossimilhança sem as variáveis é -1025,81 e com as variáveis é -928,15. Assim, o valor da estatística teste é:
$$G=-2(-1025,81)-(-2(-928,15))= 195,33.$$
O p-valor $P(\chi^{2}_{5}>\ G=195,33)= 0$.
Rejeitamos a hipótese nula. Assim, pelo TRV, temos que pelo menos uma variável testada é significativa para o modelo.
Assim, vamos testar individualmente cada uma das variáveis.
- Para a variável Volatilidade:
O log da verossimilhança sem as variáveis é -982,12 e com as variáveis é -928,15. Assim, o valor da estatística teste é:
$$G=-2(-982,12 )-(-2(-928,15))= 107,94.$$
O p-valor $P(\chi^{2}_{1}>\ G=107,94)=0$.
Rejeitamos a hipótese nula. Assim, pelo TRV, temos que a variável Volatilidade é significativa para o modelo.
- Para a variável RFV:
O log da verossimilhança sem as variáveis é -954,66 e com as variáveis é -928,15. Assim, o valor da estatística teste é:
$$G=-2(-954,66)-(-2(-928,15))= 58,02.$$
O p-valor $P(\chi^{2}_{1}>\ G=58,02)< 0,0001$.
Rejeitamos a hipótese nula. Assim, pelo TRV, temos que a variável RFV é significativa para o modelo.
- Para a variável Volatilidade ao quadrado:
O log da verossimilhança sem as variáveis é -983,77 e com as variáveis é -928,15. Assim, o valor da estatística teste é: $$G=-2(-983,77)-(-2(-928,15))= 111,24.$$
O p-valor $P(\chi^{2}_{1}>\ G=111,24)= 0$.
Pelo TRV rejeitamos a hipótese nula e por isso temos que a variável Volatilidade ao quadrado é significativa para o modelo.
- Para a variável RFV ao quadrado:
O log da verossimilhança sem as variáveis é -946,26 e com as variáveis é -928,15. Assim, o valor da estatística teste é:
$$G=-2(-946,26)-(-2(-928,15.))= 36,22.$$
O p-valor $P(\chi^{2}_{1}>\ G=36,22)< 0,0001$.
Pelo TRV rejeitamos a hipótese nula e por isso temos que a variável RFV ao quadrado é significativa para o modelo.
- Variável Volatilidade * RFV
O log da verossimilhança sem as variáveis é -936,31 e com as variáveis é -928,15. Assim, o valor da estatística teste é:
$$G=-2(-963,31)-(-2(-928,15))= 70,32.$$
O p-valor $P(\chi^{2}_{1}>\ G=70,32)= 0$.
Pelo TRV rejeitamos a hipótese nula e por isso temos que a interação de Volatilidade com RFV é significativa para o modelo.
- Teste de Wald
A estatística do teste Wald é dada pela expressão (4.2.2.2.1), sendo
$W_0=\cfrac{\widehat{\beta_0}}{\widehat{DP}(\widehat{\beta_0})}=\cfrac{30,3389}{3,66962}= 8,27$
$W_1=\cfrac{\widehat{\beta_1}}{\widehat{DP}(\widehat{\beta_1})}=\cfrac{-14,1496}{1,35514}= -10,44$
$W_2=\cfrac{\widehat{\beta}_2}{\widehat{DP}(\widehat{\beta}_2)}=\cfrac{-6,26927}{0,855191}= -7,73$
$W_3=\cfrac{\widehat{\beta}_3}{\widehat{DP}(\widehat{\beta}_3)}=\cfrac{1,10101}{0,1019943}= 10,79$
$W_4=\cfrac{\widehat{\beta}_4}{\widehat{DP}(\widehat{\beta}_4)}=\cfrac{0,2819}{0,046355}=6,08$
$W_5=\cfrac{\widehat{\beta}_5}{\widehat{DP}(\widehat{\beta}_5)}=\cfrac{1,53769}{0,18197}= 8,45$
A estimativa do desvio padrão utilizada para o cálculo da estatística de Wald está na “Estimação do Desvio Padrão”.
Para avaliar a significância dos coeficientes as seguintes hipóteses são testadas
$$\begin{cases} H_0:\beta_j=0 \cr H_1:\beta_j\neq 0\end{cases} \quad j=0,1,2,3,4,5$$
O p-valor definido como $P(|Z|> w_j)$, $j=0, 1, 2, 3, 4, 5,$ quando Z denota a variável aleatória da distribuição normal padrão.
Assim, temos:
Para $\beta_0 = P(|Z| > 8,27) = 0,000$
Para $\beta_1 = P(|Z| > 10,44) = 0,000$
Para $\beta_2 = P(|Z| > 7,73) = 0,000$
Para $\beta_3 = P(|Z| > 10,79) = 0,000$
Para $\beta_4 = P(|Z| > 6,08) = 0,000$
Para $\beta_5 = P(|Z| > 8,45) = 0,000$
Avaliando o p-valor para os $\beta_j,$ $j=0,1,2,3,4,5,$ rejeita-se $H_0$ ao nível de significância $\alpha = 0,05$ e conclui-se que os parâmetros são significativos no modelo.
A interação entre as variáveis e seus efeitos quadráticos podem ser visualizadas na Figura 4.4.4.
Figura 4.4.4: Superfície de Resposta para a probabilidade de peças refugadas.
Verificamos na Figura 4.4.4 que valores altos de Volatilidade e valores baixos de RFV resultam em menores proporções de refugo.
4.2.3 Intervalos de Confiança na Regressão Logística Múltipla
Discutimos intervalos de confiança para os parâmetros, logito, valores ajustados e odds ratio no modelo de regressão logística simples. Os métodos utilizados no caso de apenas uma covariável são os mesmos no caso da regressão múltipla.
4.2.2.1 Intervalo de Confiança para os parâmetros
O intervalo de confiança para um parâmetro $\beta_j$ é baseado em seu respectivo teste de Wald. O intervalo de confiança de $100(1-\alpha)$% para o parâmetro $\beta_j$ é:
$$IC(\beta_j,1-\alpha)= [\hat{\beta_j}-z_{1-\alpha/2}DP[(\hat{\beta_j})]; \quad \hat{\beta_j}+z_{1-\alpha/2}DP[(\hat{\beta_j})]].$$
Em que a obtenção da estimativa do parâmetro $\beta_j$ e de seu desvio padrão é mostrado em 4.2.1.
Exemplo 4.2.2.1
O modelo proposto é:
$$\hbox{Probabilidade de refugo}=\pi_i=\dfrac{e^{\beta_0 + \beta_1 x_1 +\beta_2 x_2+\beta_3 x^2_1 + \beta_4 x^2_2 + \beta_5 x_1x_2}}{1+e^{\beta_0 + \beta_1 x_1 +\beta_2 x_2+\beta_3 x^2_1 + \beta_4 x^2_2 + \beta_5 x_1x_2}},$$
em que $x_1$ é a variável Volatilidade e $x_2$ é a variável RFV.
Assim, vamos construir intervalo de confiança de Wald para os parâmetros considerando as estimativas dos parâmetros e dos seus desvios padrão já calculados nos exemplos 4.2.1.1 e 4.2.1.1.1.
- Para $\beta_0$
$$IC(\beta_0,0,95)= [30,339-1,96 * 3,67; \quad 30,339+1,96 * 3,67],$$
$$IC(\beta_0,0,95)= [30,339-7,1932; \quad 30,339+7,1932],$$
$$IC(\beta_0,0,95)= [23,14; \quad 37,53].$$
- Para $\beta_1$
$$IC(\beta_1,0,95)= [-14,1496-1,96 * 1,35514; \quad -14,1496 + 1,96 * 1,35514],$$
$$IC(\beta_1,0,95)= [-14,1496-2,65607; \quad -14,1496+2,65607],$$
$$IC(\beta_1,0,95)= [-16,8; \quad -11,50].$$
- Para $\beta_2$
$$IC(\beta_2,0,95)= [-6,26927-1,96 * 0,855191; \quad -6,26927+1,96 * 0,855191],$$
$$IC(\beta_2,0,95)= [-6,26927-1,67617; \quad -6,26927+1,67617],$$
$$IC(\beta_2,0,95)= [-7,94; \quad -4,59].$$
- Para $\beta_3$
$$IC(\beta_3,0,95)= [1,10101-1,96 * 0,101994; \quad 1,10101 + 1,96 * 0,101994],$$
$$IC(\beta_3,0,95)= [1,10101- 0,199908; \quad 1,10101+0,199908],$$
$$IC(\beta_3,0,95)= [0,90; \quad 1,30].$$
- Para $\beta_4$
$$IC(\beta_4,0,95)= [0,2819-1,96 * 0,046355; \quad 0,2819+1,96 * 0,046355],$$
$$IC(\beta_4,0,95)= [0,2819-0,090855; \quad 0,2819+0,090855],$$
$$IC(\beta_4,0,95)= [0,19; \quad 0,37].$$
- Para $\beta_5$
$$IC(\beta_5,0,95)= [1,53769 - 1,96 * 0,18197; \quad 1,53769 + 1,96 * 0,18197],$$
$$IC(\beta_5,0,95)= [1,53769 - 0,35666; \quad 1,53769 +0,35666],$$
$$IC(\beta_5,0,95)= [1,18; \quad 1,89].$$
4.2.2.2 Intervalo de Confiança para Logito
A logito é a parte linear do modelo de regressão logística. O estimador para logito é:
$$\hat{g}(X)=\hat{\beta_0}+\hat{\beta_1} X_1+\hat{\beta}_2 X_2 + \ldots +\hat{\beta}_p X_p \tag{4.2.2.2.1}$$
O estimador da variância do estimador da logito requer a obtenção da variância da soma. O caso geral é:
$$\hat{Var}[\hat{g}(x)]=\sum_{j=0}^{p}x_j^{2}\hat{Var}(\hat{\beta}_j)+\sum_{j=0}^{p}\sum_{k=j+1}^{p}2x_{j}x_{k}\hat{Cov}(\hat{\beta}_j,\hat{\beta}_k). \tag{4.2.2.2.2}$$
O intervalo de confiança para a logito é: $$IC(g(x),1-\alpha)= [\hat{g}(x)-z_{1-\alpha/2}DP[\hat{g}(x)]; \quad \hat{g}(x)+z_{1-\alpha/2}DP[\hat{g}(x)]],$$
em que $DP[\hat{g}(x)]$ é a raiz quadrada de 4.2.2.2.2 e $z_{1-\alpha/2}$ é o ponto da normal padrão.
4.2.2.3 Intervalo de Confiança para os valores ajustados
O estimador do logito e seu intervalo de confiança fornece o estimador dos valores ajustados. O intervalo de confiança dos valores ajustados é dado por:
$$IC(\pi,1-\alpha)= \left[\dfrac{e^{\hat{g}(x)-z_{1-\alpha/2}DP[\hat{g}(x)]}}{1+e^{\hat{g}(x)-z_{1-\alpha/2}DP[\hat{g}(x)]}}; \quad \dfrac{e^{\hat{g}(x)+z_{1-\alpha/2}DP[\hat{g}(x)]}}{1+e^{\hat{g}(x)+z_{1-\alpha/2}DP[\hat{g}(x)]}}\right]. \tag{4.2.2.3.1}$$
Exemplo 4.2.2.3
Vamos calcular o intervalo de confiança para o valor predito para a primeira observação que tem os seguintes valores das variáveis explicativas: Volatilidade=1,906, RFV=5,64, Volatilidade ao quadrado=3,6328, RFV ao quadrado=31,832, Interação Volatilidade e RFV=10,7536.
$$\hat{g}(x)=30,3389-14,1496 * 1,906-6,2693 * 5,64+1,1010 * 3,6328+$$
$$+0,2819 * 31,832+1,5377 * 10,7536=-2,49.$$
$$DP[\hat{g}(x)]=0,003.$$
Assim, para a primeira observação, o intervalo de 95% de confiança dos valores preditos é:
$$IC(\pi,0,95)= \left[\dfrac{e^{-2,49-(1,96 * 0,003)}}{1+e^{-2,49-(1,96 * 0,003)}}; \quad \dfrac{e^{-2,49+(1,96 * 0,003)}}{1+e^{-2,49+(1,96 * 0,003)}}\right].$$
$$IC(\pi,0,95)=[0,07; 0,08].$$
De maneira análoga para as outras observações, o Intervalo de 95% de Confiança para os valores ajustados estão na Tabela 4.4.6.
| Probabilidade Ajustada | LI | LS |
|---|---|---|
| 0,076399885 | 0,070450452 | 0,082349318 |
| 0,127990657 | 0,107519976 | 0,148461338 |
| 0,131046731 | 0,115743286 | 0,146350177 |
| 0,074710904 | 0,070150594 | 0,079271214 |
| 0,07579158 | 0,070721074 | 0,080862087 |
| 0,060636055 | 0,056779495 | 0,064492615 |
| 0,071903929 | 0,066629812 | 0,077178046 |
| 0,071022609 | 0,065386177 | 0,076659042 |
| 0,067959741 | 0,05838268 | 0,077536802 |
| 0,101027497 | 0,089405118 | 0,112649877 |
| 0,060889841 | 0,05658044 | 0,065199242 |
| 0,080881681 | 0,072272112 | 0,089491249 |
| 0,061890438 | 0,057329439 | 0,066451437 |
| 0,040125901 | 0,035798741 | 0,044453061 |
| 0,044062237 | 0,040031392 | 0,048093083 |
| 0,04366689 | 0,039566594 | 0,047767186 |
| 0,04284874 | 0,036425365 | 0,049272115 |
| 0,078649613 | 0,06629152 | 0,091007705 |
Tabela 4.4.6: Intervalo de Confiança dos Valores Ajustados
4.2.2.4 Intervalo de Confiança para a Odds Ratio
Sejam os limites do intervalo de confiança para $\beta_j$:
$\beta_{I}=\hat{\beta_j}-z_{1-\alpha/2}DP[(\hat{\beta_j})]$ e $\beta_{S}=\hat{\beta_j}+z_{1-\alpha/2}DP[(\hat{\beta_j})],$
em que j=1,2,…p.
O intervalo de confiança para a Odds Ratio é:
$$IC(\hbox{Odds Ratio},1-\alpha)= [e^{\beta_{I}}; \quad ~e^{\beta_{S}}]. \tag{4.2.2.4.1}$$
Exemplo 4.2.2.4
Vamos calcular o intervalo de confiança para os dados do exemplo 4.2.1.1.
- Para Volatilidade
$$\left[\exp(-14,1496-1,96 * 1,35514) \right.~;~~\left.\exp(-14,1496 + 1,96 * 1,35514)\right]$$
$$\left[\exp(-14,1496- 2,65607)\right.~;~~\left.\exp(-14,1496+ 2,65607)\right]$$
$$\left[5,027942e-08; \quad 1,019585e-05\right].$$
- Para RFV
$$\left[\exp(-6,26927-1,96 * 0,855191)\right.~;~~\left.\exp(-6,26927+1,96 * 0,855191)\right]$$
$$\left[\exp(-6,26927-1,67617)\right.~;~~\left.\exp(-6,26927+1,67617)\right]$$
$$\left[0,00035427; \quad 0,010121\right].$$
- Para o Efeito quadrático da Volatilidade
$$\left[\exp(1,10101-1,96 * 0,101994); \quad \exp(1,10101 + 1,96 * 0,101994) \right]$$
$$\left[\exp(1,10101- 0,199908) ; \quad \exp(1,10101+0,199908)\right]$$
$$\left[2,4623 ; \quad 3,6726 \right].$$
- Para o Efeito quadrático do RFV
$$\left[\exp(0,2819-1,96 * 0,046355) ; \quad \exp(0,2819+1,96 * 0,046355)\right]$$
$$\left[\exp(0,2819-0,090855) ; \quad \exp(0,2819+0,090855)\right]$$
$$\left[1,2105 ; \quad 1,4517\right].$$
- Para a Interação entre a Volatilidade e o RFV
$$\left[\exp(1,53769 - 1,96 * 0,18197) ; \quad \exp(1,53769 + 1,96 * 0,18197)\right]$$
$$\left[\exp(1,53769 - 0,35666) ; \quad \exp(1,53769 +0,35666)\right]$$
$$\left[3,2577 ; \quad 6,64822\right].$$
4.2.4 Variáveis independentes categóricas
4.2.4.1 Variável Independente Dicotômica
Vamos considerar variáveis explicativas assumindo valores discretos. Por exemplo, considere uma única variável explicativa $X$ com níveis A ou B. Para estimar os parâmetros do modelo, estas variáveis são substituídas por valores numéricos, por exemplo 0 para o nível A e 1 para o nível B. Esses valores são chamados de codificação das categorias da variável qualitativa. Assim, a mesma teoria desenvolvida para o modelo logístico considerando variáveis explicativas contínuas pode ser aplicado para esses casos. Através dos valores da Tabela 4.4.7, a odds Ratio será interpretado por:
$$OR=\cfrac{\cfrac{\pi(1)}{1 - \pi(1)}}{\cfrac{\pi(0)}{1 - \pi(0)}} = \cfrac{\cfrac{\cfrac{e^{\beta_0+\beta_1}}{1+e^{\beta_0+\beta_1}}}{\cfrac{1}{1+e^{\beta_0+\beta_1}}}}{\cfrac{\cfrac{e^{\beta_0}}{1+e^{\beta_0}}}{\cfrac{1}{1+e^{\beta_0}}}}=exp(\beta_1) \tag{4.2.4.1.1}$$
| x=1 | x=0 | |
|---|---|---|
| y=1 | $\pi(1)=\cfrac{exp(\beta_0+\beta_1)}{1+exp(\beta_0+\beta_1)}$ | $\pi(0)=\cfrac{exp(\beta_0)}{1+exp(\beta_0)}$ |
| y=0 | $1-\pi(1)=\cfrac{1}{1+exp(\beta_0+\beta_1)}$ | $1-\pi(0)=\cfrac{1}{1+exp(\beta_0)}$ |
| Total | 1 | 1 |
Tabela 4.4.7: Valores do Modelo de Regressão Logística quando a variável independente é dicotômica.
Portanto podemos fazer algumas análises,
$$OR> 1 \quad \Rightarrow~~\cfrac{\pi_{1}}{1-\pi_{1}}~~>~~\cfrac{\pi_{0}}{1-\pi_{0}}~~\Rightarrow~~\pi_B~~> ~~\pi_A$$
$$OR< 1 \quad \Rightarrow ~~ \cfrac{\pi_1}{1 - \pi_{1}}~~<~~\cfrac{\pi_{0}}{1 - \pi_{0}}~~\Rightarrow~~ \pi_B~~<~~\pi_A$$
Antes de fazer conclusões a respeito da odds ratio encontrada, é importante considerar o efeito que a codificação da variável tem no valor da estimativa da odds. Vimos que $\hat{OR}=exp(\hat{\beta_1})$. Mas esse resultado é válido apenas quando a variável independente é codificada em 0 ou 1.
A estimativa do log da odds ratio para uma variável independente com dois níveis a e b é:
$$ln[\hat{OR}(a,b)]=\hat{g}(x=a)-\hat{g}(x=b)$$
$$=(\hat{\beta_0}+\hat{\beta_1}.a)-(\hat{\beta_0}+\hat{\beta_1}.b)$$
$$=\hat{\beta_1}(a-b).$$
Assim, a estimativa da odds ratio é dada por:
$$\hat{OR}(a,b)=exp[\hat{\beta_1}(a-b)].$$
O método da codificação também influencia nos limites do intervalo de confiança. Assim, o intervalo de confiança é:
$$IC(OR,1-\alpha)= [exp(\hat{\beta_1}(a-b)-z_{1-\alpha/2}|a-b|DP[(\hat{\beta_1})]); \quad exp(\hat{\beta_1}(a-b)+z_{1-\alpha/2}|a-b|DP[(\hat{\beta_1})])].$$
Em que $|a-b|$ é o valor absoluto de $(a-b)$. Desde que possamos controlar a codificação das variáveis dicotômicas (duas categorias), Hosmer e Lemeshow (1976) recomendam a codificação em 0 ou 1.
4.2.4.2 Variáveis com mais de duas categorias
Suponha que ao invés de duas categorias, a variável independente tem mais de duas categorias ($k> 2$). Essa variável poderia ser, por exemplo, a raça: branco, preto, amarelo ou outros. Vamos supor que a Tabela 4.4.8 é a classificação cruzada da variável raça (x) em relação á variável resposta (y). Vale frisar que os dados são fictícios.
| branco | preto | amarelo | outros | Total | |
|---|---|---|---|---|---|
| y=1 | 5 | 20 | 15 | 10 | 50 |
| y=0 | 20 | 10 | 10 | 10 | 50 |
| 25 | 30 | 25 | 20 | 100 |
Tabela 4.4.8: Classificação cruzada dos dados hipotéticos da raça em relação a variável resposta binária.
Assim como no caso da variável independente dicotômica, precisamos codificar a variável explicativa. Mas quando temos mais de duas categorias, temos as variáveis codificadas que são chamadas de variáveis dummies. Além disso, no caso de variável com mais de duas categorias temos que especificar a de referência em que a odds ratio das demais categorias são comparadas com ela. Na Tabela 4.4.9 temos a codificação da variável raça.
D: Variáveis Dummies
| Raça | D1 | D2 | D3 |
|---|---|---|---|
| Amarelo | 0 | 0 | 0 |
| Branco | 1 | 0 | 0 |
| Outros | 0 | 1 | 0 |
| Preto | 0 | 0 | 1 |
Tabela 4.4.9: Codificação da variável raça
Vale ressaltar que temos k categorias da variável explicativa, temos então k-1 variáveis dummies. Assim, o preditor linear é dado por:
$$\beta_0+\beta_1D1+\beta_2D2+\beta_3D3.$$
A variável de referência é aquela em que os valores de todas as dummies é 0, que no caso é a raça amarelo.
A comparação da raça branco com amarelo é:
$$ln[\hat{OR}\hbox{(Branco,Amarelo)}]=\hat{g}\hbox{(Branco)}-\hat{g}\hbox{(Amarelo)}=$$
$$=[\hat{\beta_0}+\hat{\beta_1}(D1=1)+\hat{\beta}_2(D2=0)+\hat{\beta}_3(D3=0)]-[\hat{\beta_0}+\hat{\beta_1}(D1=0)+\hat{\beta}_2(D2=0)+\hat{\beta}_3(D3=0)]$$
$$=[\hat{\beta_0}+\hat{\beta_1}]-[\hat{\beta_0}]$$
$$=\hat{\beta_1}.$$
Assim, a estimativa da Odds Ratio(Branco, Amarelo) é:
$$\hat{OR}\hbox{(Branco,Amarelo)}=exp(\hat{\beta_1}).$$
Agora comparando a raça outros em relação à amarelo, temos:
$$ln[\hat{OR}\hbox{(Outros,Amarelo)}]=\hat{g}\hbox{(Outros)}-\hat{g}\hbox{(Amarelo)}=$$
$$=[\hat{\beta_0}+\hat{\beta_1}(D1=0)+\hat{\beta}_2(D2=1)+\hat{\beta}_3(D3=0)]-[\hat{\beta_0}+\hat{\beta_1}(D1=0)+\hat{\beta}_2(D2=0)+\hat{\beta}_3(D3=0)]$$
$$=[\hat{\beta_0}+\hat{\beta}_2]-[\hat{\beta_0}]$$
$$=\hat{\beta}_2.$$
Assim, a estimativa da Odds Ratio(Outros, Amarelo) é:
$$\hat{OR}\hbox{(Outros,Amarelo)}=exp(\hat{\beta}_2).$$
Analogamente, temos que:
$$\hat{OR}\hbox{(Preto,Amarelo)}=exp(\hat{\beta}_3).$$
Exemplo 4.2.3.1
Um experimento DOE foi proposto para avaliar a influência das váriáveis em um sistema de medição. Após a realização de um brainstorming verificou-se que as variáveis mais influentes no modelo foram o preparador e a sonda . Para a análise deste sistema de medição realizou-se um experimento completo com dois fatores e dois níveis. Os dados coletados estão organizados na Tabela 4.4.10.
| Concordância | Não-Concordância | Sonda | Preparador | Total |
|---|---|---|---|---|
| 928 | 32 | A | A | 960 |
| 862 | 98 | B | B | 960 |
| 830 | 130 | A | B | 960 |
| 932 | 28 | B | A | 960 |
Tabela 4.4.10: Dados do sistema de medição.
Para este exemplo, o seguinte modelo foi proposto:
$$\hbox{Probabilidade de falha}=\pi_i = \frac{\exp[\beta_0 + \beta_1~x_1 +\beta_2~x_2]}{1+\exp[\beta_0 + \beta_1~x_1 +\beta_2~x_2]}$$
Considerando a variável $x_1$ = Preparador e $x_2$ = Sonda, temos:
| Níveis | Sonda(x1) | Preparador(x2) |
|---|---|---|
| A | 0 | 0 |
| B | 1 | 1 |
Tabela 4.4.11: Codificação das variáveis qualitativas.
A partir da função de máxima verossimilhança, os parâmetros do modelos foram obtidos.
$\widehat{\beta_0}$ =$-3,30347$ (Intercepto)
$\widehat{\beta_1}$ =$-0,27926$ (Sonda - Nível B)
$\widehat{\beta}_2$ =$1,43136$ (Preparador - Nível B)
Da tabela (4.4.11) calculamos a probabilidade de falha do sistema de medição.
Para $x_1 = 0$ e $x_2 = 0$ ( Nível A - Sonda e Nível A - Preparador)
$$\hbox{Probabilidade de falha}=\pi_i = \frac{\exp[-3,30347 -0,27926~(0) +1,43136~(0)]}{1+\exp[-3,30347 -0,27926~(0) +1,43136~(0)]} = 0,03544$$
Para $x_1 = 1$ e $x_2 = 0$ ( Nível B - Sonda e Nível A - Preparador)
$$\hbox{Probabilidade de falha}=\pi_i = \frac{\exp[-3,30347 -0,27926~(1) +1,43136~(0)]}{1+\exp[-3,30347 -0,27926~(1) +1,43136~(0)]} = 0,027041$$
Para $x_1 = 0$ e $x_2 = 1$ ( Nível A - Sonda e Nível B - Preparador)
$$\hbox{Probabilidade de falha}=\pi_i = \frac{\exp[-3,30347 -0,27926~(0) +1,43136~(1)]}{1+\exp[-3,30347 -0,27926~(0) +1,43136~(1)]} = 0,13329$$
Para $x_1 = 1$ e $x_2 = 1$ ( Nível B - Sonda e Nível B - Preparador)
$$\hbox{Probabilidade de falha}=\pi_i = \frac{\exp[-3,30347 -0,27926~(1) +1,43136~(1)]}{1+\exp[-3,30347 -0,27926~(1) +1,43136~(1)]} = 0,104202$$
E as estimativas da probabilidade de ocorrência do evento de interesse são:
| Sonda | Preparador | Probabilidade |
|---|---|---|
| A | A | 0,035452 |
| B | B | 0,104202 |
| A | B | 0,133298 |
| B | A | 0,027048 |
Tabela 4.4.12: Probabilidade do evento.
As estimativas da odds ratio são:
Sonda: $\hat{OR}(B,A) = exp(\hat{\beta_1})=0,76$
Preparador: $\hat{OR}(B,A)=exp(\hat{\beta}_2)= 4,18$.
Antes de calcular o intervalo de confiança para a odds ratio, é preciso calcular as estimativas para o desvio padrão, por meio da matriz $\Sigma(\beta)=[X’VX]^{-1}$:
$$\Sigma(\beta) = \left[ \begin{bmatrix} 1 \quad 1 \quad 1 \quad 1 \cr 0 \quad 1 \quad 0 \quad 1 \cr 0 \quad 1 \quad 1 \quad 0 \cr \end{bmatrix}{4x3} \times \begin{bmatrix} 960 * 0,035(1-0,035) \quad 0 \quad \ldots \quad 0 \cr 0 \quad 960 * 0,104(1-0,104) \quad \ldots \quad 0 \cr 0 \quad 0 \quad 960 * 0,133(1-0,133) \quad 0 \cr \quad 0 \quad \ldots \quad 960 * 0,027(1-0,027) \end{bmatrix}{4x4}\times \right.$$
$$\left. \times \begin{bmatrix} 1 \quad 0 \quad 0 \cr 1 \quad 1 \quad 1 \cr 1 \quad 0 \quad 1 \cr 1 \quad 1 \quad 0 \cr \end{bmatrix}{4x3}^{-1}= \begin{bmatrix} \mathbf{0,020177} \quad -0,00681 \quad -0,01713 \cr -0,00681 \quad \mathbf{0,015664} \quad -0,00019 \cr -0,01713 \quad -0,00019 \quad \mathbf{0,022204} \end{bmatrix}{2x2} \right]$$
As estimativas dos desvios padrão são:
Para $\beta_0 = \sqrt{ \sigma^2(\widehat{\beta_0}) } = \sqrt{0,0201 } =0,142045$
Para $\beta_1 = \sqrt{ \sigma^2(\widehat{\beta_1}) } = \sqrt{0,015662 } = 0,125156$
Para $\beta_2 = \sqrt{ \sigma^2(\widehat{\beta}_2) } = \sqrt{0,02220 } = 0,149009$
Os resultados da estatística do teste Wald são:
$W_0=\cfrac{\widehat{\beta_0}}{\widehat{DP}(\widehat{\beta_0})}=\cfrac{-3,30347}{0,142045}=-23,26$
$W_1=\cfrac{\widehat{\beta_1}}{\widehat{DP}(\widehat{\beta_1})}=\cfrac{-0,279268}{0,125156}= -2,23$
$W_2=\cfrac{\widehat{\beta}_2}{\widehat{DP}(\widehat{\beta}_2)}=\cfrac{1,43136}{0,149009}= 9,61$
A significância dos parâmetros são testados considerando as seguintes hipóteses:
$$\begin{cases} H_0:\beta_j=0 \cr H_1:\beta_j \neq 0 \cr \end{cases} \quad ~j=0, 1, 2.$$
Os valores de p-valor para o teste de hipótese são:
Para $\beta_0 = P(|Z| > 23,26) = 0,000$
Para $\beta_1 = P(|Z| > 2,23) = 0,026$
Para $\beta_2 = P(|Z| > 7,73) = 0,000$
Dessa forma, rejeitamos $H_0$ a um nível de significância $\alpha = 0,05$ e concluímos assim que os parâmetros $\beta_0$, $\beta_1$ e $\beta_2$ são significativos no modelo.
Os intervalos de confiança para a Odds Ratio das variáveis são:
Para $\beta_1$ ( Sonda - tipo B )
$$\left[\exp(-0,2792 - 1,96 * (0,1251)) ; \quad \exp(-0,2792 + 1,96 * (0,1251)) \right]$$
$$\left[\exp(-0,2792 -0,2451) ; \quad \exp(-0,2792 +0,2451)\right]$$
$$\left[\exp(-0,5243) ; \quad \exp(-0,0341) \right]$$
$$\left[0,591 ; \quad 0,97\right]$$
Para $\beta_2$ ( Preparador - tipo B )
$$\left[\exp(1,4313 - 1,96 * (0,149)) ; \quad \exp(1,4313 + 1,96 * (0,149)) \right]$$
$$\left[\exp(1,4313 -0,2920) ; \quad \exp(1,4313 +0,2920) \right]$$
$$\left[\exp(1,1393) ; \quad \exp(1,723) \right]$$
$$\left[3,124 ; \quad 5,60 \right]$$
Como o Intervalo de Confiança da Odds Ratio do nível B da sonda em relação ao nível A é menor que 1, concluímos que a sonda B apresenta menos chance de falha no sistema de medição quando comparada a sonda A.
Como o Intervalo de Confiança da Odds Ratio do nível B do preparador em relação ao nível A é maior que 1, concluímos que o preparador B apresenta maior chance de falha no sistema de medição quando comparado ao preparador A.
4.2.5 Seleção de Variáveis
A abordagem tradicional na construção de modelos estatísticos é encontrar o modelo mais parcimonioso que explica os dados. Quanto mais variáveis no modelo, maior se torna a estimativa do erro e mais dependente o modelo fica dos dados observados.
Existem algumas técnicas para auxiliar na seleção de variáveis para um modelo de Regressão Logística. Vamos considerar duas delas: Método de Seleção Stepwise e Todos os Modelos Possíveis.
4.2.5.1 Seleção Stepwise
O método Stepwise para a seleção de variáveis é muito usado em regressão linear.
Qualquer procedimento para seleção or exclusão de variáveis de um modelo é baseado em um algoritmo que checa a importância das variáveis, incluindo ou excluindo-as do modelo se baseando em uma regra de decisão. A importância da variável é definida em termos de uma medida de significância estatística do coeficiente associado à variável para o modelo. Essa estatística depende das suposições do modelo. No Stepwise da regressão linear um teste F é usado desde que os erros tenham distribuição normal. Na regressão logística os erros seguem distribuição binomial e a significância é assegurada via Teste da Razão de Verossimilhança. Assim, em cada passo do procedimento a variável mais importante, em termos estatísticos, é aquela que produz a maior mudança no logaritmo da verossimilhança em relação ao modelo que não contém a variável.
Temos a seguir o algoritmo do método de Stepwise passo a passo;
- Passo 0: Suponha que temos p variáveis explicativas candidatas ao modelo. Esse passo começa com o ajuste apenas do intercepto e seja $L_0$ o log da verossimilhança desse ajuste. Após isso, ajustamos os p modelos com apenas uma variável explicativa. $L_{j}^{(0)}$ é o log da verossimilhança do modelo contendo a variável $x_j$. O valor do teste da Razão de Verossimilhança do modelo contendo $x_j$ versus o modelo com apenas o intercepto é $G_j^{(0)}=-2(L_0-L_j^{(0)})$ e o p-valor é $p_j^{(0)}=P[\chi_{v}^2> G_j^{(0)}]$, em que v=1 se $x_j$ é contínuo e v=k-1 se $x_j$ é categórico com k categorias.
Seja $p_{e_1}^{(0)}$ o p-valor associado ao teste da variável $x_1$ e além disso, $p_{e_1}^{(0)}=min(p_j^{(0)})$, ou seja, é o menor p-valor de todos os testes da Razão de Verossimilhança. Se $p_{e_1}^{(0)}<\alpha_e$ então $x_1$ entra no modelo.
$\alpha_e$ é o nível de significância para verificar se a variável entra ou não no modelo. Lee e Koval (1997) examinaram a questão da significância para a regressão logística. Os resultados mostram que a escolha $\alpha_e=0,05$ é muito rigorosa e por isso pode excluir variáveis importantes do modelo. É indicado escolher um valor de $\alpha_e$ entre 0,15 e 0,20.
-
Passo 1: Ajustamos agora o modelo contendo $x_1$. Seja $L_{e_1}^{(1)}$ o log da verossimilhança desse modelo. Para verificar se p-1 variáveis são importantes para o modelo, uma vez que $x_1$ está nele, ajustamos p-1 modelos de regressão contendo $x_1$ e $x_j$. O log da verossimilhança é denotada por $L_{e_1;j}^{(1)}$ e a estatística teste da Razão de Verossimilhança é: $G_j^{(1)}=-2(L_{e_1}^{(1)}-L_{e_1;j}^{(1)})$. Suponha que $p_{e_2}^{(1)}=min(p_j^{(1)})$, ou seja, o p-valor do teste associado à variável $x_2$ é o menor dentre todos os outros. Se esse p-valor for menor que $\alpha_e$ $x_2$ entra no modelo e prosseguimos para o passo 2, caso contrário o algoritmo termina e apenas a variável $x_1$ foi incluída no modelo.
-
Passo 2: O passo 2 começa com o ajuste do modelo contendo $x_1$ e $x_2$. É possível que, com a inclusão de $x_2$, $x_1$ passa a ser não significativa para o modelo. Por isso, nesse passo testamos a significância de uma das variáveis dado que a outra está no modelo. Seja $L_{-e_j}^{(2)}$ o log da verossimilhança do modelo com $x_j$ removido. A estatística da razão de verossimilhança é $G_{-e_j}^{(2)}=-2(L_{-e_j}^{(2)}-L_{e_1;e_2}^{(2)})$ e $p_{-e_j}^{(2)}$ é o seu p-valor. Para decidir se removemos ou não alguma das variáveis, selecionamos o maior p-valor. Vamos supor que o p-valor da variável $x_2$ seja o maior, ou seja, $p_{r_2}^{(2)}=max(p_{-e_1}^{(2)},p_{-e_2}^{(2)})$. Para ver se a variável $x_2$ sai do modelo, comparamos o seu p-valor com um nível de significância de saída, ou seja, $\alpha_r$. Se $p_{r_2}^{(2)}>\alpha_r$ então $x_2$ sai do modelo. Contudo, determinamos o valor de $\alpha_r$ que seja maior que $\alpha_e$ para não acontecer do programa remover e incluir a mesma variável no modelo em sucessivos passos. No método Stepwise em modelos de regressão logística, é recomendado $\alpha_e=0,15$ e $\alpha_r=0,2.$
Após a verificação da eliminação da variável $x_2$ e supondo que ela não saiu do modelo, ou seja, seu p-valor não é maior que $\alpha_r$, ajustamos p-2 modelos de regressão logística contendo $x_{e_1}$,$x_{e_2}$ e $x_j$ com j=1,2,..,p. Seja $x_3$ a variável associada ao menor p-valor no teste da razão de verossimilhança, isto é, $p_{e_3}^{(2)}=min(p_j^{(2)}).$ Se o p-valor é menor que $\alpha_e$, $x_3$ entra no modelo e prosseguimos para o passo 3, caso contrário paramos por aqui e o modelo é explicado por $x_1$ e $x_2$.
-
Passo 3: O procedimento do passo 3 é idêntico ao do passo 2. O programa ajusta o modelo incluindo a variável selecionada durante os passos anteriores, checando se a inclusão dessa variável fez com que as outras variáveis do modelo perdessem a significância. O processo continua dessa mesma maneira até o último passo, o passo S.
-
Passo S: Esse passo ocorre quando todas as p variáveis entraram no modelo ou se todas variáveis no modelo tem p-valor para sair menor que $\alpha_r$ e as variáveis não incluídas no modelo tem p-valor para entrar maior que $\alpha_e$.
4.2.5.2 Todos os modelos possíveis
No contexto de regressão linear múltipla, discutimos os critérios $R^2_p$, $R^2_a$, $C_p$, $AIC$, $BIC$ e $PRESS$ para a seleção de variáveis levando em conta todos os modelos possíveis. Na regressão logística, os critérios $AIC$ e $BIC$ são facilmente adaptados. Por essa razão, focaremos nesses dois critérios para a avaliação de todos os modelos possíveis. Essas medidas são calculadas da seguinte maneira:
$$AIC=-2log_eL(\beta)+2p.$$
$$BIC=-2log_eL(\beta)+plog_e(n).$$
Em que $log_eL(\beta)$ é o log da verossimilhança do modelo (veja as estimativas dos parâmetros no modelo de regressão logística múltipla).
Um outro critério é $-2log_eL(\beta)$ que não é penalizado quando adicionamos preditores, já que ele não depende do número de covariáveis no modelo. Já os critérios AIC e BIC se deixam levar pelo número de variáveis explicativas, já que no cálculo do AIC temos a soma de $2p$ no valor do critério e no BIC $plog_e(n)$.
4.3 Medidas da qualidade do ajuste do modelo
Podemos avaliar a qualidade do ajuste do modelo fazendo um gráfico do modelo ajustado, mas existem medidas, como a deviance, a estatística qui-quadrada de Pearson e o teste de Hosmer e Lemeshow que ajudam a verificar a qualidade do modelo.
Para os testes de Deviance e qui-quadrado de Pearson, precisamos identificar, dentre as observações, aquelas com os mesmos valores entre si para todos as covariáveis. Assim, seja $\mathbf {X}=(X_1,X_2,…,X_p)$ o conjunto de variáveis explicativas do modelo. Seja também J o número distinto de valores que x assume na amostra (combinações de níveis diferentes das covariáveis, chamado de Covariate Patterns). O número de indivíduos na amostra com valores iguais $\mathbf{x}=\mathbf{x}_j$ é denotado por: $m_j$. Já $y_j$ é o número de indivíduos positivos, isto é, y=1, em $m_j$. Consequentemente, $\sum y_j=n_1$ que representa o total de indivíduos com y=1.
Já para o teste de Hosmer e Lemeshow não é necessário a existência de combinações de valores das covariáveis.
A seguir apresentamos cada um desses testes para avaliar a qualidade do ajuste.
4.3.1 Resíduos de Pearson
Na regressão linear, medidas da qualidade do ajuste são funções dos resíduos definido como a diferença entre o observado e valores ajustados $(y-\hat{y})$. Na regressão logística existem diferentes formas de calcular essa diferença. Para entendermos que os valores ajustados na regressão logística são calculados para cada Covariate Pattern (veja a definição em “Medidas de Qualidade do Ajuste do Modelo”) e depende das estimativas de probabilidade para cada Covariate Pattern, vamos denotar os valores ajustados para o j-ésimo Covariate Pattern como $\hat{y}_j$ em que:
$$\hat{y}_j=m_j\hat{\pi}_j,$$
em que $m_j$ é o número de observações na Covariate Pattern j e $\hat{\pi}_j$ é a probabilidade ajustada dos indivíduos em j.
A medida de Pearson para a diferença entre o observado e predito é:
$$r(y_j,\hat{\pi}_j)=\cfrac{(y_j-m_j\hat{\pi}_j)}{\sqrt{m_j\hat{\pi}_j(1-\hat{\pi}_j)}},$$
em que $y_j$ o número de indivíduos em j com y=1.
Assim, a estatística qui-quadrado de Pearson é dada por:
$$\chi^2=\sum_{j=1}^{J}{r(y_j,\hat{\pi}_j)}^2.$$
Podemos dizer que $\chi^2$ se aproxima, assintoticamente, $\chi^2_{J-(p+1)},$ $p$ o número de covariáveis do modelo ajustado e $J$ é o número de Covariate Pattern. Porém, em alguns casos essa aproximação é ruim (veja Restrições do teste).
Exemplo 4.3.1.1
- Regressão Logística Simples: Exemplo 4.1.2.1
Para calcular a estatística qui-quadrado de Pearson para testar se o modelo ajustado é adequado, temos que calcular, para cada grupo, o resíduo de Pearson.
Para o primeiro grupo (como são dados resumidos, a primeira “observação”), temos que:
$$\hat{y}_1=m_1\hat{\pi}_1=200 * 0,01806=3,612, \quad y_1=2.$$
$$r(y_1,\hat{\pi}_1)=\cfrac{(y_1-m_1\hat{\pi}_1)}{\sqrt{m_1\hat{\pi}_1(1-\hat{\pi}_1)}}=\cfrac{(2-3,612)}{\sqrt{200 * 0,01806(0,98194)}}=\frac{-1,612}{1,8833}=-0,8559.$$
De maneira análoga calculamos os resíduos de Pearson para os outros grupos e obtemos os seguintes valores:
| Grupos (observações resumidas) | Resíduos de Pearson |
|---|---|
| 1 | -0,85599 |
| 2 | -0,85599 |
| 3 | -0,85599 |
| 4 | -0,954 |
| 5 | -0,55487 |
| 6 | -0,18804 |
| 7 | 0,147992 |
| 8 | 0,147992 |
| 9 | 0,454625 |
| 10 | 0,294173 |
| 11 | 0,98504 |
| 12 | 0,429086 |
| 13 | 0,429086 |
| 14 | 0,429086 |
| 15 | 0,230969 |
| 16 | 0,230969 |
| 17 | 0,584557 |
| 18 | 0,379551 |
| 19 | 0,507824 |
| 20 | 0,079737 |
| 21 | 0,170431 |
| 22 | -0,05228 |
| 23 | -0,05228 |
| 24 | 0,24343 |
| 25 | 0,012207 |
| 26 | 0,012207 |
| 27 | -0,219 |
| 28 | -0,45045 |
| 29 | -0,17915 |
| 30 | -0,89921 |
Tabela 4.4.13: Resíduos de Pearson
Com os valores dos resíduos da Tabela 4.4.13, a estatística qui-quadrado de Pearson é:
$$\chi^2=\sum_{j=1}^{30}{r(y_j,\hat{\pi}_j)}^2=7,3535.$$
O p-valor $P(\chi^2_{28} > 7,3535)=0,99$ e por isso não rejeitamos a hipótese de que o modelo ajustado é adequado.
- Regressão Logística Múltipla: Exemplo 4.2.1.
Para obtermos a estatística qui-quadrado de Pearson para testar se o modelo ajustado é adequado, temos que calcular, para cada grupo, o resíduo de Pearson.
Assim, para o primeiro grupo (como são dados resumidos, a primeira “observação”), temos que:
$$\hat{y}_1=m_1\hat{\pi}_1=832 * 0,0764=63,56, \quad y_1=270.$$
$$r(y_1,\hat{\pi}_1)=\cfrac{(y_1-m_1\hat{\pi}_1)}{\sqrt{m_1\hat{\pi}_1(1-\hat{\pi}_1)}}=\cfrac{(270-63,56)}{\sqrt{832 * 0,0764(0,9236)}}=\frac{206,44}{7,66214}=26,94.$$
De maneira análoga calculamos os resíduos de Pearson para os outros grupos e obter os seguintes valores:
| Grupos (observações resumidas) | Resíduos de Pearson |
|---|---|
| 1 | 26,94226 |
| 2 | 2,325756 |
| 3 | 10,89268 |
| 4 | -7,29709 |
| 5 | -12,1222 |
| 6 | -3,41077 |
| 7 | -6,6327 |
| 8 | -14,1813 |
| 9 | -2,78226 |
| 10 | -8,84807 |
| 11 | -5,97549 |
| 12 | 13,87763 |
| 13 | -0,44592 |
| 14 | 9,839347 |
| 15 | 6,995839 |
| 16 | 5,203155 |
| 17 | 2,974618 |
| 18 | -5,2033 |
Tabela 4.4.14: Resíduos de Pearson
Com os valores dos resíduos da Tabela 4.4.14, temos que a estatística qui-quadrado de Pearson é:
$$\chi^2=\sum_{j=1}^{18}{r(y_j,\hat{\pi}_j)}^2=1830,15.$$
O p-valor $(P(\chi^2_{12} > 1830,15))$ é 0 e por isso rejeitamos a hipótese de que o modelo ajustado é adequado. Assim, temos indícios pelo teste qui-quadrado de Pearson que o modelo ajustado não é adequado.
4.3.2 Deviance
A soma dos quadrados dos resíduos (SQE) da regressão linear, no modelo logístico é denominada deviance (D), definida como:
$$d(y_j,\widehat{\pi}_j)=\pm\left[2\left[y_j~\ln \left(\frac{y_j}{m_j\widehat{\pi}_j} \right)+(m_j - y_j) ~\ln \left(\frac{(m_j -y_j)}{m_j\widehat{\pi}_j(1-\widehat{\pi}_j)} \right) \right]\right]^{1/2}. \tag{4.3.2.1}$$
Em que o sinal + ou - é o mesmo sinal que $(y_j-m_j\hat{\pi}_j)$. Além disso, $m_j$ é o número de observações na Covariate Pattern j, $\hat{\pi}_j$ é a probabilidade ajustada dos indivíduos no grupo j e
$y_j$ o número de indivíduos em j com y=1.
A estatística teste é:
$$D=\sum_{j=1}^J{d(y_j,\hat{\pi}_j)}^2.$$
A estatística D, sob a suposição que o modelo ajustado é correto, tem distribuição assintótica $\chi^2$ com $J-(p+1)$ graus de liberdade.
Restrições do teste Qui-Quadrado de Pearson e Deviance
Quando os resultados assintóticos são baseados no fato em que n cresce mantendo os $m_j$ pequenos (ou tornando os $m_j$ pequenos), chamamos de n-assintótico.
Quando $J< n$ é fixado e n cresce fazendo com que os $m_j$ também cresçam, os resultados assintóticos baseados no $m_j$ são chamados de m-assintótico.
Assim, quando $J\approx n$, sob n-assintótico, os p-valores calculados dos testes Deviance e qui-quadrado de Pearson, usando a distribuição $\chi^2$ com $(J-p-1)$ graus de liberdade estão incorretos.
Uma saída é criar agrupamentos que permite o uso de resultados m-assintóticos.
Vale ressaltar que para a realização do teste de qualidade de ajuste de Hosmer e Lemeshow não é necessário existir Covariate Pattern.
Exemplo 4.3.2.1
- Regressão Logística Simples: Exemplo 4.1.2.1
Para obtermos a estatística Deviance para testar se o modelo ajustado é adequado, temos que calcular pela equação (4.4.15), para cada grupo (como são dados resumidos, a primeira “observação”), o resíduo Deviance.
Assim, temos os valores dos resíduos Deviance:
| Grupos (observações) | Resíduos Deviance |
|---|---|
| 1 | -0,93425 |
| 2 | -0,93425 |
| 3 | -0,93425 |
| 4 | -1,05002 |
| 5 | -0,58267 |
| 6 | -0,19086 |
| 7 | 0,146409 |
| 8 | 0,146409 |
| 9 | 0,440923 |
| 10 | 0,288473 |
| 11 | 0,929236 |
| 12 | 0,418699 |
| 13 | 0,418699 |
| 14 | 0,418699 |
| 15 | 0,228179 |
| 16 | 0,228179 |
| 17 | 0,56749 |
| 18 | 0,372416 |
| 19 | 0,495682 |
| 20 | 0,079443 |
| 21 | 0,169152 |
| 22 | -0,0524 |
| 23 | -0,0524 |
| 24 | 0,240936 |
| 25 | 0,012201 |
| 26 | 0,012201 |
| 27 | -0,22097 |
| 28 | -0,45869 |
| 29 | -0,18041 |
| 30 | -0,92968 |
Tabela 4.4.15: Resíduos Deviance
Com os valores dos resíduos da Tabela 4.4.15, temos que a estatística Deviance é:
$$D=\sum_{j=1}^{30}{d(y_j,\hat{\pi}_j)}^2=7,8722.$$
O p-valor $P(\chi^2_{28} > 7,8722)=0,99$. Assim, não rejeitamos a hipótese de que o modelo é adequado.
- Regressão Logística Múltipla: Exemplo 4.2.1.1
Para obtermos a estatística Deviance para testar se o modelo ajustado é adequado, temos que calcular pela equação (4.3.2.1) os resíduos Deviance.
Assim, os valores dos resíduos são:
| Grupos (observações) | Resíduos Deviance |
|---|---|
| 1 | 20,72168 |
| 2 | 2,26703 |
| 3 | 9,965709 |
| 4 | -9,65112 |
| 5 | -15,6351 |
| 6 | -3,83042 |
| 7 | -8,40627 |
| 8 | -17,3252 |
| 9 | -2,89436 |
| 10 | -10,7692 |
| 11 | -6,8788 |
| 12 | 12,25469 |
| 13 | -0,44928 |
| 14 | 8,503884 |
| 15 | 6,305022 |
| 16 | 4,730169 |
| 17 | 2,818639 |
| 18 | -5,68624 |
Tabela 4.4.16: Resíduos Deviance
Com os valores dos resíduos da Tabela 4.4.16, temos que a estatística Deviance é:
$$D=\sum_{j=1}^{18}{d(y_j,\hat{\pi}_j)}^2=1753,714.$$
O p-valor $P(\chi^2_{12} > 1753,714)=0$. Assim, rejeitamos a hipótese de que o modelo é adequado.
4.3.3 Teste de Hosmer e Lemeshow
O teste de Hosmer-Lemeshow é muito utilizado em regressão logística com a finalidade de testar a bondade do ajuste, em outras palavras, o teste comprova se o modelo proposto pode explicar bem o que se observa. O teste avalia o modelo ajustado através das distâncias entre as probabilidades ajustadas e as probabilidades observadas.
A bondade do teste é baseada na divisão da amostra segundo suas probabilidades ajustadas com base nos valores dos parâmetros estimados pela regressão logística. Os valores ajustados são dispostos do menor para o maior, e em seguida, separados em g grupos de tamanho aproximadamente igual. Hosmer e Lemeshow (1980) propõe que seja utilizado $g=10$.
Suponha que o agrupamento foi feito em 10 grupos. Desta forma o primeiro grupo será formado pelas observações com probabilidades ajustadas de até $10%$. O segundo grupo será formado pelas observações com probabilidades ajustadas de até $10%$ seguintes, ou seja, probabilidades entre $10$ e $20%$. Assim a divisão é feita até que se obtenha os 10 grupos.
Na literatura há pouca orientação sobre como escolher o número de grupos. As simulações mostradas em Hosmer e Lemeshow (1980) foram baseadas no uso de $g> p+1$, em que p é o número de covariáveis do modelo ajustado. Se as frequências esperadas em alguns dos grupos forem muito pequenas a estatística do teste de Hosmer-Lemeshow é calculada, entretanto pode não ser confiável. Neste caso, devemos especificar um número menor de grupos, contudo não se pode utilizar menos de 3 grupos, pois com $(g<3)$ a estatística do teste é impossibilitada de ser calculada.
Antes do cálculo da estatística teste, é necessário estimar a frequência esperada dentro de cada grupo, para isso dividimos a variável resposta, que é dicotômica. Para Y=1, a frequência esperada estimada é dada pela soma das probabilidades estimadas de todos os indivíduos dentro daquele grupo. Para Y=0, a frequência esperada estimada é dada pela soma de 1-probabilidade estimada de todos os indivíduos dentro daquele grupo. A seguir temos o exemplo de uma tabela que deve ser preenchida com as frequências esperadas para realizar o cálculo da estatística de Hosmer e Lemeshow.
| Decil | Y=1 Obs | Y=1 Esp | Y=0 Obs | Y=0 Esp | Total |
|---|---|---|---|---|---|
| 1 | |||||
| 2 | |||||
| 3 | |||||
| 4 | |||||
| 5 | |||||
| 6 | |||||
| 7 | |||||
| 8 | |||||
| 9 | |||||
| 10 |
Tabela 4.4.17: exemplo a ser preenchido para realizar o cálculo da estatística de Hosmer e Lemeshow
Tendo as frequências esperadas calculamos a estatística de Hosmer e Lemeshow, $\hat{C}$, que é obtida da seguinte forma:
$$\hat{C}=\sum_{k=1}^{g}\frac{({o_{k}-n^{'}_{k}\bar{\pi}_{k}})^2}{n^{'}_{k}\pi_{k}(1-\pi_{k})},$$
em que:
-
$n^{'}_{k}$ é o número de indivíduos no k-ésimo grupo.
-
$\bar{\pi_{k}}=\sum_{j=1}^{C_{k}}\frac{m_{j}\bar{\pi}_{j}}{n^{'}_{k}}$
-
$C_k$: o número total de combinações de níveis dentro do k-ésimo decil.
-
$O_k=\sum_{j=1}^{C_{k}}y_j$: número total de respostas dentro do grupo k.
Hosmer e Lemeshow (1980) mostrou por simulação que a estatística do teste segue, aproximadamente, uma distribuição qui-quadrado com g−2 graus de liberdade, quando o modelo está especificado corretamente.
4.4 Diagnóstico do Modelo
4.4.1 Leverage
Assim como no modelo linear, uma métrica para diagnosticar outliers é a leverage (diagonal da matriz chapéu). No modelo linear, a matriz chapéu é definida por:
$$H=X(X^\prime X)^{-1}X^\prime.$$
Os elementos da diagonal principal da matriz H são denominados $h_{ii}.$
No modelo de regressão logística, J é o número de diferentes combinações de níveis das covariáveis.
$X_{(j)(p+1)}$: matriz de todos os valores de todas combinações de níveis para todas as covariáveis.
$$H=V^{1/2}X(X^\prime VX)^{-1}X^\prime V^{1/2}.$$
$v_{j}=m_{j}\hat{\pi}(x_j)(1-\hat{\pi}(x_j)), \quad \quad \quad ~j=1,2,…,J.$
$h_i = m_j \hat{\pi} _j (1-\hat{\pi} _j) x^{^\prime }_{j}(X^\prime VX)^{-1}x^{^\prime}_j,$
em que $b_j=x^{^\prime }_{j}(X^\prime VX)^{-1}x^{^\prime }_j.$
Assim, $h_j=v_jb_j$.
No caso da regressão linear, se $h_j> 2(p+1)/n$ então a observação j é considerada um outlier. A diferença é que neste caso a matriz chapéu depende apenas de valores das covariáveis enquanto que na regressão logística, essa matriz também depende das probabilidades estimadas (através de $v_j$).
A consequência disso é que uma observação pode estar fora do comum nas covariáveis, mas não possui um grande valor de $h_j$, se a probabilidade ajustada estiver próxima de 0 ou 1.
4.4.2 Resíduo de Pearson
O resíduo de Pearson é dado por:
$$r_{Pi}=\frac{Y_i-\hat{\pi}_i}{\sqrt{\hat{\pi}_i(1-\hat{\pi}_i)}}$$
4.4.3 Resíduo de Pearson Studentizado
O resíduo de Pearson Studentizado tem a seguinte forma:
$$r_{SPi}=\frac{r_{Pi}}{\sqrt{1-h_i}},$$
em que $h_i$ é a diagonal da matriz chapéu.
4.4.4 Resíduo Deviance
O resíduo deviance é:
$$d_i=\pm\sqrt{-2[Y_ilog_e(\hat{\pi}_i)+(1-Y_i)log_e(1-\hat{\pi}_i)}.$$
4.4.5 Diagnóstico de Influência
Além de detectar outliers, é importante também detectar pontos influentes, ou seja, pontos que afetam de forma significativa o ajuste do modelo.
Assim como no modelo linear, podemos utilizar a distância de Cook para avaliar a influência geral da observação i nas estimativas dos coeficientes da regressão. A métrica da distância de Cook para a observação i é:
$$D_i=\frac{r^{2}_{SPi}h_i}{(p+1)(1-h_i)}.$$
Notamos que $ D_i $ é grande quando o resíduo $ e_i $ e/ou a leverage $ h_{ii} $ são grandes. O critério utilizado para destacar as observações influentes é $ D_i~>1 $, desta forma, todas as distâncias maiores que 1 são consideradas influentes.
4.5 Predição
Quando temos a variável resposta binária (1 se o indivíduo for evento e 0 caso contrário) é necessário escolher uma regra de predição ($\hat{Y}=0 \hbox{ ou } 1$), já que $\hat{\pi}$ está entre 0 e 1. É intuitivo pensar que se o valor de $\hat{\pi}_i$ for grande, $\hat{Y}_i=1$ e se $\hat{\pi}_i$ for pequeno, $\hat{Y}_i=0$. Mas como determinar o ponto que para os valores acima dele o indivíduo é classificado como evento ($\hat{Y}_i=1$) e valores abaixo dele o indivíduo é classificado como não evento ($\hat{Y}_i=0$) ? Esse ponto é conhecido como ponto de corte.
Uma forma bastante utilizada para determinar o ponto de corte é através da Curva ROC (Receiver Operating Characteristic Curve). A curva ROC plota $P(\hat{Y}=1|Y=1)$ (chamado de sensibilidade) versus $1-P(\hat{Y}=0|Y=0)$ (chamado de 1-especificidade) para todos os possíveis pontos de corte entre 0 e 1.
A seguir temos um exemplo da curva ROC:
Figura 4.4.5: Curva ROC
A escolha do ponto de corte deve ser baseada em uma combinação ótima tanto da sensibilidade quanto da especificidade, pois partimos do suposto que classificar o indivíduo como evento dado que ele é não evento (falso positivo) e classificar o indivíduo como não evento dado que ele é evento (falso negativo) traz prejuízos equivalentes para o pesquisador. Pela análise da curva ROC, escolhemos o ponto de corte referente a combinação da sensibilidade e 1-especificidade que mais se aproxima do canto superior esquerdo do gráfico.
4.5.1 Métricas de Desempenho da Predição
Após o ajuste de um modelo e a determinação do ponto de corte, é importante avaliar o poder de discriminação do modelo, isto é, discriminar os eventos dos não eventos.
Para essa avaliação, métricas foram criadas. São elas: Acurácia, Sensibilidade, Especificidade, Verdadeiro Preditivo Positivo e Verdadeiro Preditivo Negativo. Mas antes de entrar especificamente nessas medidas, precisamos apresentar a matriz de confusão, dada na Tabela 4.5.1. Seu funcionamento é simples: é uma tabela de contingência em que na linha está o valor previsto e na coluna o valor observado (valor verdadeiro).
| Valor Observado (valor verdadeiro) | ||
|---|---|---|
| Y=1 | Y=0 | |
| Valor Predito Y=1 | VP (verdadeiro positivo) | FP (falso positivo) |
| Valor Predito Y=0 | FN (falso negativo) | VN (verdadeiro negativo) |
Tabela 4.4.18: Matriz de Confusão
4.5.1.1 Acurácia
É a proporção de predições corretas, sem considerar o que é positivo e o que negativo e sim o acerto total. É dada por:
$$ACC=(VP+VN)/(P+N),$$
em que p é o número de eventos (Y=1, chamado aqui de positivo) e n é o número de não eventos (Y=0, chamado aqui de negativo).
4.5.1.2 Sensibilidade
É a proporção de verdadeiros positivos, ou seja, avalia a capacidade do modelo classificar um indivíduo como evento $(\hat{Y}=1)$ dado que realmente ele é evento (Y=1):
$$SENS=VP/(VP+FN).$$
4.5.1.3 Especificidade
É a proporção de verdadeiros negativos, isto é, avalia a capacidade do modelo predizer um indivíduo como não evento $(\hat{Y}=0)$ dado que ele realmente é não evento (Y=0).
$$ESPEC=VN/(VN+FP).$$
4.5.1.4 Verdadeiro Preditivo Positivo
É a proporção de verdadeiros positivos em relação a todas as predições positivas, isto é, o indivíduo ser evento (Y=1) dado que o modelo classificou o indivíduo como evento $(\hat{Y}=1)$
$$VPP=VP/(VP+FP).$$
4.5.1.5 Verdadeiro Preditivo Negativo
É a proporção de verdadeiros negativos em relação a todas predições negativas, ou seja, o indivíduo ser não evento (Y=0) dado que o modelo o classificou como não evento $(\hat{Y}=0)$.
$$VPN=VN/(VN+FN).$$