---

O problema da (quase) separação na análise de dados binários

---

\[\mathbf{Prof.\ Cesar\ Augusto\ Taconeli}\]

\[\mathbf{Modelos\ Lineares\ Generalizados\ (CE225)}\]

\[\mathbf{Outubro/2017}\]

---

Definição

O problema da separação, na regressão para dados binários, se caracteriza pela existência deum hiperplano no espaço das covariáveis que separe perfeitamente as observações com \(y=0\) das observações com \(y=1\).

---

A toy example

x <- 1:6
y <- c(1, 1, 1, 0, 0, 0)
plot(x, y, pch = 20)
abline(v = 3.5, col = 'red', lty = 2, lwd = 1.5)

Observe que a reta \(x=3.5\) separa completamente observações com cada valor da resposta. Vamos ajustar um MLG com distribuição binomial.

ajuste <- glm(y ~ x, family = 'binomial')
summary(ajuste)

Call:
glm(formula = y ~ x, family = "binomial")

Deviance Residuals: 
         1           2           3           4           5           6  
 2.110e-08   2.110e-08   1.052e-05  -1.052e-05  -2.110e-08  -2.110e-08  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)    165.32  407521.43       0        1
x              -47.23  115264.41       0        1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8.3178e+00  on 5  degrees of freedom
Residual deviance: 2.2152e-10  on 4  degrees of freedom
AIC: 4

Number of Fisher Scoring iterations: 25

confint.default(ajuste)

                2.5 %   97.5 %
(Intercept) -798562.0 798892.6
x           -225961.3 225866.9

predict(ajuste, type = 'response')

           1            2            3            4            5 
1.000000e+00 1.000000e+00 1.000000e+00 5.537833e-11 2.220446e-16 
           6 
2.220446e-16 

Observe que as estimativas pontuais e os erros padrões (principalmente) são excessivamente grandes. Como resultado, o resultado do teste de Wald é contraditório e os intervalos de confiança baseados nesse teste são absurdamente largos. Isso é resultado de estimativas para \(\pi\) iguais a zero e iguais a 1.

Nesse caso, se fixarmos \(\hat{\beta}_0 = -3.5\hat{\beta}_1\) (tal que = 0.5 em \(x=3.5\)) e fizermos \(\hat{\beta} \rightarrow -\infty\), temos uma sequência de ajustes

### Sequência de ajustes com maior log-verossimilhança.
inv_logit <- function(beta0, beta1, x) exp(beta0 + beta1 * x)/(exp(beta0 + beta1 * x) + 1)
curve(inv_logit(3.5, -1, x), 0, 6, las = 1, ylim = c(-0.1, 1.1), ylab = 'P(y=1)')
curve(inv_logit(7, -2, x), 0, 6, add = TRUE, col = 'green')
curve(inv_logit(14, -4, x), 0, 6, add = TRUE, col = 'red')
curve(inv_logit(56, -16, x), 0, 6, add = TRUE, col = 'orange')
curve(inv_logit(280, -80, x), 0, 6, add = TRUE, col = 'blue')
points(x, y, cex = 1.2, pch = 20)