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

163

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?

user333
fonte
4
pergunta
user603
1
pergunta relacionada e demonstração sobre a regularização aqui
Haitao Du

Respostas:

100

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.

Nick Sabbe
fonte
20
Você também pode usar uma interface de fórmula para glmnet através do pacote de interpolação.
Zach
"Agora, quando você espera ..." Pergunta sobre isso. Eu tenho um estudo de caso / controle olhando para o relacionamento com o microbioma. Também temos um tratamento que quase só é encontrado entre os casos. No entanto, acreditamos que o tratamento também possa afetar o microbioma. Este é um exemplo de sua ressalva? Hipoteticamente, poderíamos encontrar muitos outros casos que não usariam o tratamento se tentássemos, mas temos o que temos.
abalter 27/09
142

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. 15212

  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.)

Scortchi
fonte
1
@ Scortchi Existe outra opção (herética). Que tal mudar algumas observações selecionadas aleatoriamente de 1 para 0 ou 0 para 1 entre variáveis ​​que exibem separação completa?
robertf
@RobertF: Obrigado! Eu não tinha pensado nisso - se você tiver alguma referência a respeito de seu desempenho, eu ficaria agradecido. Você já se deparou com pessoas que o usam na prática?
Scortchi
@ Scortchi - Não, existem referências a pesquisadores que adicionam dados artificiais para eliminar a separação completa, mas não encontrei nenhum artigo sobre modificação seletiva dos dados. Não tenho idéia de quão eficaz esse método seria.
robertf
1
@tatami: Nem todos os programas (muitos?) alertam sobre a separação em si, o que pode ser difícil de detectar quando se trata de uma combinação linear de várias variáveis, mas sobre falhas de convergência e / ou valores ajustados próximos de nada ou nenhum - eu sempre verifique estes.
Scortchi
2
@ Scortchi: resumo muito bom em sua resposta. Pessoalmente, sou a favor da abordagem bayesiana, mas vale a pena mencionar a bela análise do fenômeno geral de um ponto de vista freqüentista em projecteuclid.org/euclid.ejs/1239716414 . O autor oferece alguns intervalos de confiança unilaterais que podem ser usados ​​mesmo na presença de separação completa na regressão logística.
Cyan
55

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:

insira a descrição da imagem aqui

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

Rasmus Bååth
fonte
6
Não . Simples demais . Você pode explicar o que acabou de fazer? Qual é o prior que bayesglmusa? Se a estimativa de ML é equivalente a Bayesiana com um plano simples, como os priores não informativos ajudam aqui?
StasK 28/02
5
Adicionado mais algumas informações! O prior é vago, mas não plano. Ele tem alguma influência, pois regulariza as estimativas e as puxa levemente para 0,0, o que acredito que você deseja neste caso.
Rasmus Bååth
> m = bayesglm (correspondência ~., família = binomial (link = 'logit'), dados = df) Mensagem de aviso: probabilidades ajustadas numericamente ocorreram 0 ou 1 Não é bom!
Chris
Como iniciante, tente uma regularização um pouco mais forte aumentando os prior.dfpadrões para 1.0e / ou diminuindo os prior.scalepadrões 2.5, talvez comece a tentar:m=bayesglm(match ~. , family = binomial(link = 'logit'), data = df, prior.df=5)
Rasmus Bååth
1
O que exatamente estamos fazendo quando aumentamos o prior.df no modelo. Existe um limite para o quão alto queremos ir? Meu entendimento é que ele restringe o modelo para permitir a convergência com estimativas precisas de erro?
hamilthj
7

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

Mike Hunter
fonte
3

warning

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.

insira a descrição da imagem aqui

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)

χ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

AdamO
fonte
2
Isso é ótimo, mas eu tenho algumas queixas, é claro: (1) O teste da razão de verossimilhança não usa a matriz de informações; é apenas o teste de Wald que falha e falha catastroficamente na presença de separação. (2) Não estou familiarizado com os estimadores de uma etapa, mas a estimativa da inclinação aqui parece absurdamente baixa. (3) Um intervalo de confiança não é isento de mediana. O que você vincula nessa seção é o intervalo de confiança de meio p. (4) Você pode obter intervalos de confiança invertendo os testes LR ou de pontuação. ...
Scortchi
... (5) Você pode executar o teste de pontuação em R, apresentando o argumento test="Rao"para a anovafunção. (Bem, as duas últimas são anotações, não
brincadeiras
@ scortchi bom saber que anova tem testes de pontuação padrão! Talvez uma implementação manual seja útil. Os ICs não são medianos e não-imparciais, mas os ICs para o estimador não-mediano fornecem inferência consistente para parâmetros de contorno. O meio p é um estimador. O p pode ser transformado para uma razão de chances b / c, é invariável para transformações um para um. O teste LR é consistente para os parâmetros de contorno?
Adamo
Somente a hipótese nula não deve conter parâmetros em um limite para o teorema de Wilks, embora os testes de pontuação e LR sejam aproximados em amostras finitas.
Scortchi
2

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.

Manoel Galdino
fonte
8
O problema com a remoção de preditores é que você está removendo o preditor que melhor explica a resposta, que geralmente é o que você deseja fazer! Eu argumentaria que isso só faz sentido se você superajustar seu modelo, por exemplo, ajustando muitas interações complicadas.
Simon Byrne
4
Não é um bug, mas um problema com as estimativas iniciais muito distantes do MLE, que não surgirão se você não tentar escolhê-las.
Scortchi
Eu entendo isso, mas acho que isso é um bug no algoritmo.
Manoel Galdino
5
Bem, eu não quero discutir sobre a definição de 'bug'. Mas o comportamento não é insondável nem impossível de fixar na base R - você não precisa "verificar com outro software". Se você deseja lidar automaticamente com muitos problemas de não convergência, o glm2pacote implementa uma verificação de que a probabilidade está realmente aumentando a cada etapa da pontuação e reduz pela metade o tamanho da etapa, se não estiver.
Scortchi
3
Existe (no CRAN) o pacote R, safeBinaryRegression projetado para diagnosticar e corrigir esses problemas, usando métodos de otimização para verificar com certeza se há separação ou quase-separação. Tente!
Kjetil b halvorsen
2

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.

StayLearning
fonte
-1

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.

yash
fonte
2
Não vejo que esta resposta acrescente muito à discussão. A abordagem bayesiana é amplamente abordada nas respostas anteriores, a remoção de preditores "problemáticos" também já é mencionada (e desencorajada). A seleção de variáveis ​​por etapas raramente é uma ótima idéia, até onde eu sei.
einar