Como lidar com a separação perfeita na regressão logística?

Se você tiver uma variável que separa perfeitamente os zeros e os da variável de destino, R produzirá a seguinte mensagem de aviso "separação perfeita ou quase perfeita":

Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 

Ainda temos o modelo, mas as estimativas do coeficiente são infladas.

Como você lida com isso na prática?

Uma solução para isso é utilizar uma forma de regressão penalizada. De fato, esta é a razão original pela qual algumas das formas de regressão penalizadas foram desenvolvidas (embora tenham revelado outras propriedades interessantes).

Instale e carregue o pacote glmnet no R e você estará pronto para começar. Um dos aspectos menos amigáveis ​​do glmnet é que você só pode alimentar matrizes, não fórmulas como estamos acostumados. No entanto, você pode olhar model.matrix e similares para construir essa matriz a partir de um data.frame e uma fórmula ...

Agora, quando você espera que essa separação perfeita não seja apenas um subproduto da sua amostra, mas possa ser verdadeira na população, você especificamente não deseja lidar com isso: use essa variável de separação simplesmente como o único preditor do seu resultado, não empregando um modelo de qualquer tipo.

Você tem várias opções:

1. Remova alguns dos preconceitos.

   (a) Penalizando a probabilidade de acordo com a sugestão de @ Nick. O pacote logistf em R ou a FIRTHopção no SAS PROC LOGISTICimplementam o método proposto em Firth (1993), "Redução de viés das estimativas de máxima verossimilhança", Biometrika , 80 , 1 .; que remove o viés de primeira ordem das estimativas de probabilidade máxima. ( Aqui, o @Gavin recomenda o brglmpacote, com o qual não estou familiarizado, mas deduzi que implementa uma abordagem semelhante para funções de link não canônicas, como probit.)

   (b) Usando estimativas isentas de mediana na regressão logística condicional exata. Pacote elrm ou Logistix em R, ou a EXACTdeclaração em SAS de PROC LOGISTIC.

2. Exclua os casos em que a categoria ou valor do preditor que está causando a separação ocorre. Isso pode estar fora do seu escopo; ou digno de uma investigação mais aprofundada. (O pacote safeBinaryRegression do R é útil para encontrá-los.)

3. Relance o modelo. Normalmente, isso é algo que você teria feito antes, se tivesse pensado nisso, porque é muito complexo para o tamanho da sua amostra.

   (a) Remova o preditor do modelo. Arriscado, pelas razões apresentadas pelo @Simon: "Você está removendo o preditor que melhor explica a resposta".

   (b) Ao recolher categorias de preditores / agrupar os valores dos preditores. Somente se isso fizer sentido.

   (c) Re-expressar o preditor como dois (ou mais) fatores cruzados sem interação. Somente se isso fizer sentido.

4. Use uma análise bayesiana de acordo com a sugestão de Manoel . Embora pareça improvável que você queira apenas por causa da separação, vale a pena considerar por outros méritos. O artigo que ele recomenda é Gelman et al (2008), "Uma distribuição anterior padrão pouco informativa para modelos logísticos e outros modelos de regressão", Ann. Appl. Estado. , 2 , 4 : o padrão em questão é um Cauchy independente anterior para cada coeficiente, com média de zero e escala de ; a ser usado após padronizar todos os preditores contínuos para ter uma média de zero e um desvio padrão de . Se você pode elucidar priors fortemente informativos, tanto melhor.$\frac{5}{2}$$\frac{1}{2}$

5. Fazer nada. (Mas calcule os intervalos de confiança com base nas probabilidades do perfil, pois as estimativas de erro padrão de Wald estarão muito erradas.) Uma opção frequentemente negligenciada. Se o objetivo do modelo é apenas descrever o que você aprendeu sobre as relações entre preditores e resposta, não há vergonha em citar um intervalo de confiança para uma razão de chances de, digamos, 2,3 ou mais. (De fato, pode parecer suspeito citar intervalos de confiança com base em estimativas imparciais que excluem as razões de chances mais suportadas pelos dados.) Os problemas surgem quando você está tentando prever o uso de estimativas pontuais, e o preditor em que a separação ocorre envolve os outros.

6. Use um modelo de regressão logística oculto, conforme descrito em Rousseeuw e Christmann (2003), "Robustez contra separação e outliers na regressão logística", Computational Statistics & Data Analysis , 43 , 3 e implementado no pacote R hlr . (@ user603 sugere isso. ) Eu não li o artigo, mas eles dizem no resumo "um modelo um pouco mais geral é proposto sob o qual a resposta observada está fortemente relacionada, mas não é igual à resposta verdadeira inobservável", o que sugere Para mim, talvez não seja uma boa ideia usar o método, a menos que isso pareça plausível.

7. "Altere algumas observações selecionadas aleatoriamente de 1 para 0 ou 0 para 1 entre variáveis ​​que exibem separação completa": comentário de RobertF . Essa sugestão parece surgir de considerar a separação como um problema em si, e não como um sintoma de escassez de informações nos dados, o que pode levar você a preferir outros métodos à estimativa de probabilidade máxima ou limitar as inferências àquelas que você pode fazer. precisão razoável - abordagens que têm seus próprios méritos e não são apenas "correções" para a separação. (Além de ser descaradamente ad hoc , é desagradável para a maioria os analistas que fazem a mesma pergunta dos mesmos dados, fazendo as mesmas suposições, devem dar respostas diferentes devido ao resultado de um sorteio ou algo assim.)

Essa é uma expansão das respostas de Scortchi e Manoel, mas como você parece usar o RI, pensei em fornecer algum código. :)

Acredito que a solução mais fácil e direta para o seu problema é usar uma análise bayesiana com premissas prévias não informativas, conforme proposto por Gelman et al (2008). Como Scortchi menciona, Gelman recomenda colocar um Cauchy anterior com mediana de 0,0 e escala de 2,5 em cada coeficiente (normalizado para ter média de 0,0 e DP de 0,5). Isso regularizará os coeficientes e os puxará levemente para zero. Nesse caso, é exatamente o que você deseja. Por possuir caudas muito largas, o Cauchy ainda permite grandes coeficientes (em oposição ao Normal de cauda curta), de Gelman:

Como executar esta análise? Use a bayesglmfunção no pacote arm que implementa essa análise!

library(arm)

set.seed(123456)
# Faking some data where x1 is unrelated to y
# while x2 perfectly separates y.
d <- data.frame(y  =  c(0,0,0,0, 0, 1,1,1,1,1),
                x1 = rnorm(10),
                x2 = sort(rnorm(10)))

fit <- glm(y ~ x1 + x2, data=d, family="binomial")

## Warning message:
## glm.fit: fitted probabilities numerically 0 or 1 occurred 

summary(fit)
## Call:
## glm(formula = y ~ x1 + x2, family = "binomial", data = d)
##
## Deviance Residuals: 
##       Min          1Q      Median          3Q         Max  
## -1.114e-05  -2.110e-08   0.000e+00   2.110e-08   1.325e-05  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -18.528  75938.934       0        1
## x1              -4.837  76469.100       0        1
## x2              81.689 165617.221       0        1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1.3863e+01  on 9  degrees of freedom
## Residual deviance: 3.3646e-10  on 7  degrees of freedom
## AIC: 6
## 
## Number of Fisher Scoring iterations: 25

Não funciona tão bem ... Agora a versão bayesiana:

fit <- bayesglm(y ~ x1 + x2, data=d, family="binomial")
display(fit)
## bayesglm(formula = y ~ x1 + x2, family = "binomial", data = d)
##             coef.est coef.se
## (Intercept) -1.10     1.37  
## x1          -0.05     0.79  
## x2           3.75     1.85  
## ---
## n = 10, k = 3
## residual deviance = 2.2, null deviance = 3.3 (difference = 1.1)

Super simples, não?

Referências

Gelman et al (2008), "Uma distribuição prévia padrão pouco informativa para modelos de regressão logística e outros", Ann. Appl. Stat., 2, 4 http://projecteuclid.org/euclid.aoas/1231424214

Uma das explicações mais completas sobre questões de "separação quase completa" com a máxima probabilidade é o artigo de Paul Allison. Ele está escrevendo sobre o software SAS, mas os problemas que ele aborda são generalizáveis ​​para qualquer software:

• A separação completa ocorre sempre que uma função linear de x pode gerar previsões perfeitas de y

• A separação quase completa ocorre quando (a) existe algum vetor de coeficiente b tal que bxi ≥ 0 sempre que yi = 1 e bxi ≤ 0 * sempre que ** yi = 0 e essa igualdade se aplica a pelo menos um caso em cada categoria do variável dependente. Em outras palavras, no caso mais simples, para qualquer variável independente dicotômica em uma regressão logística, se houver um zero na tabela 2 × 2 formada por essa variável e a variável dependente, a estimativa de ML para o coeficiente de regressão não existe.

Allison discute muitas das soluções já mencionadas, incluindo exclusão de variáveis ​​problemáticas, categorias em colapso, não fazendo nada, alavancando a regressão logística exata , estimativa bayesiana e estimativa de probabilidade máxima penalizada.

http://www2.sas.com/proceedings/forum2008/360-2008.pdf

warning$\infty$

Com os dados gerados ao longo das linhas de

x <- seq(-3, 3, by=0.1)
y <- x > 0
summary(glm(y ~ x, family=binomial))

O aviso é feito:

Warning messages:
1: glm.fit: algorithm did not converge 
2: glm.fit: fitted probabilities numerically 0 or 1 occurred 

o que reflete muito obviamente a dependência que é construída nesses dados.

Em R, o teste de Wald é encontrado com summary.glmou com waldtesta lmtestembalagem. O teste da razão de verossimilhança é realizado com anovaou com lrtesto lmtestpacote. Nos dois casos, a matriz de informações é infinitamente avaliada e nenhuma inferência está disponível. Em vez disso, R não produzir saída, mas você não pode confiar nele. A inferência que R normalmente produz nesses casos tem valores de p muito próximos de um. Isso ocorre porque a perda de precisão no OR é de ordens de magnitude menor que a perda de precisão na matriz de variância-covariância.

Algumas soluções descritas aqui:

Use um estimador de uma etapa,

Há muita teoria apoiando o baixo viés, eficiência e generalização dos estimadores de uma etapa. É fácil especificar um estimador de uma etapa em R e os resultados geralmente são muito favoráveis ​​para previsão e inferência. E esse modelo nunca diverge, porque o iterador (Newton-Raphson) simplesmente não tem a chance de fazê-lo!

fit.1s <- glm(y ~ x, family=binomial, control=glm.control(maxit=1))
summary(fit.1s)

Dá:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.03987    0.29569  -0.135    0.893    
x            1.19604    0.16794   7.122 1.07e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Assim, você pode ver as previsões refletindo a direção da tendência. E a inferência é altamente sugestiva das tendências que acreditamos serem verdadeiras.

faça um teste de pontuação,

A estatística Score (ou Rao) difere da razão de verossimilhança e estatística wald. Não requer uma avaliação da variância sob a hipótese alternativa. Nós ajustamos o modelo sob o nulo:

mm <- model.matrix( ~ x)
fit0 <- glm(y ~ 1, family=binomial)
pred0 <- predict(fit0, type='response')
inf.null <- t(mm) %*% diag(binomial()$variance(mu=pred0)) %*% mm
sc.null <- t(mm) %*% c(y - pred0)
score.stat <- t(sc.null) %*% solve(inf.null) %*% sc.null ## compare to chisq
pchisq(score.stat, 1, lower.tail=F)

$\chi^2$

> pchisq(scstat, df=1, lower.tail=F)
             [,1]
[1,] 1.343494e-11

Nos dois casos, você tem inferência para um OR do infinito.

e use estimativas imparciais imparciais para um intervalo de confiança.

Você pode produzir um IC 95% não-singular mediano, não-singular para o índice de chances infinitas, usando estimativa não-verbal mediana. O pacote epitoolsem R pode fazer isso. E eu dou um exemplo de implementação deste estimador aqui: Intervalo de confiança para amostragem de Bernoulli

Tenha cuidado com esta mensagem de aviso da R. Dê uma olhada neste post de Andrew Gelman, e você verá que nem sempre é um problema de separação perfeita, mas às vezes é um erro glm. Parece que se os valores iniciais estiverem muito longe da estimativa de probabilidade máxima, ele explodirá. Portanto, verifique primeiro com outro software, como o Stata.

Se você realmente tiver esse problema, tente usar a modelagem bayesiana, com prévios informativos.

Mas, na prática, acabo com os preditores que causam o problema, porque não sei como escolher um informativo prévio. Mas acho que há um artigo de Gelman sobre o uso de informações informativas antes, quando você tem esse problema do problema de separação perfeita. Basta pesquisar no Google. Talvez você deva tentar.

Não tenho certeza se concordo com as declarações da sua pergunta.

Acho que essa mensagem de aviso significa que, para alguns dos níveis observados de X em seus dados, a probabilidade ajustada é numericamente 0 ou 1. Em outras palavras, na resolução, ela aparece como 0 ou 1.

Você pode executar predict(yourmodel,yourdata,type='response')e encontrará zeros ou zeros como probabilidades previstas.

Como resultado, acho que não há problema em usar apenas os resultados.

Entendo que este é um post antigo, no entanto, continuarei respondendo a isso, pois lutei dias com ele e ele pode ajudar outras pessoas.

A separação completa acontece quando as variáveis ​​selecionadas para ajustar-se ao modelo podem diferenciar com muita precisão entre 0 e 1 ou sim e não. Toda a nossa abordagem da ciência de dados é baseada na estimativa de probabilidade, mas falha neste caso.

Etapas de retificação: -

1. Use bayesglm () em vez de glm (), caso a variação entre as variáveis ​​seja baixa

2. Às vezes, usar (maxit = "algum valor numérico") junto com o bayesglm () pode ajudar

3. Terceira e mais importante verificação das variáveis ​​selecionadas para o ajuste do modelo, deve haver uma variável cuja colinearidade com a variável Y (saída) seja muito alta, descarte essa variável do seu modelo.

Como no meu caso, eu tinha dados de rotatividade de telecomunicações para prever a rotatividade dos dados de validação. Eu tinha uma variável nos meus dados de treinamento que poderia diferenciar muito entre sim e não. Depois de soltá-lo, consegui o modelo correto. Além disso, você pode usar stepwise (fit) para tornar seu modelo mais preciso.