Comment gérer la séparation parfaite dans la régression logistique?

163

Si vous avez une variable qui sépare parfaitement les zéros de la variable cible, R affichera le message d’alerte suivant: "séparation parfaite ou quasi parfaite":

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

Nous obtenons toujours le modèle mais les estimations des coefficients sont gonflées.

Comment gérez-vous cela dans la pratique?

utilisateur333
la source
4
question
user603
1
question connexe et démo sur la régularisation ici
Haitao Du

Réponses:

100

Une solution à ce problème consiste à utiliser une forme de régression pénalisée. En fait, c'est la raison pour laquelle certaines des formes de régression pénalisées ont été développées (bien qu'elles se soient avérées avoir d'autres propriétés intéressantes.

Installez et chargez le paquet glmnet dans R et vous êtes prêt à partir. L'un des aspects les moins conviviaux de glmnet est que vous ne pouvez alimenter que des matrices, pas les formules comme nous en avons l'habitude. Cependant, vous pouvez regarder model.matrix et autres pour construire cette matrice à partir d'un nom de données.fr et d'une formule ...

Maintenant, lorsque vous vous attendez à ce que cette séparation parfaite ne soit pas simplement un sous-produit de votre échantillon, mais qu'elle puisse être vraie dans la population, vous ne voulez surtout pas vous en occuper: utilisez cette variable de séparation simplement comme prédicteur de votre résultat, et non employant un modèle de toute sorte.

Nick Sabbe
la source
20
Vous pouvez également utiliser une interface de formule pour glmnet via le package caret.
Zach
"Maintenant, quand vous vous attendez à ..." Question à ce sujet. J'ai une étude de cas / contrôle sur la relation avec le microbiome. Nous avons également un traitement qui se trouve presque uniquement parmi les cas. Cependant, nous pensons que le traitement pourrait également affecter le microbiome. Est-ce un exemple de votre mise en garde? Hypothétiquement, nous pourrions trouver beaucoup plus de cas qui n'utilisent pas le traitement si nous essayons, mais nous avons ce que nous avons.
Abalter
142

Vous avez plusieurs options:

  1. Supprimer une partie du biais.

    (a) En pénalisant la probabilité selon la suggestion de @ Nick. La logistique des packages dans R ou l' FIRTHoption dans SAS PROC LOGISTICimplémente la méthode proposée par Firth (1993), "Réduction du biais du maximum de vraisemblance estimé", Biometrika , 80 , 1 .; ce qui élimine le biais de premier ordre des estimations du maximum de vraisemblance. ( Ici, @Gavin recommande le brglmpaquet, que je ne connais pas bien, mais je suppose qu'il implémente une approche similaire pour les fonctions de liaison non canoniques, par exemple probit.)

    (b) En utilisant des estimations sans biais médianes dans la régression logistique conditionnelle exacte. Paquet elrm ou logistiX dans R, ou la EXACTdéclaration dans SAS PROC LOGISTIC.

  2. Exclure les cas où la catégorie ou la valeur de prédicteur causant la séparation se produit. Celles-ci peuvent bien être hors de votre portée; ou digne d'une enquête plus poussée et ciblée. (Le paquetage safeBinaryRegression est pratique pour les trouver.)

  3. Re-cast le modèle. C'est généralement quelque chose que vous auriez fait auparavant si vous y aviez pensé, car c'est trop complexe pour la taille de votre échantillon.

    (a) Supprimez le prédicteur du modèle. Dicey, pour les raisons données par @Simon: "Vous supprimez le prédicteur qui explique le mieux la réponse".

    (b) En regroupant les catégories de prédicteurs / en regroupant les valeurs de prédicteurs. Seulement si cela a du sens.

    (c) Ré-exprimer le prédicteur sous forme de deux (ou plus) facteurs croisés sans interaction. Seulement si cela a du sens.

  4. 5212

  5. Ne fais rien. (Mais calculez les intervalles de confiance en fonction des probabilités du profil, car les estimations de Wald de l'erreur type seront très erronées.) Une option souvent négligée. Si le but du modèle est simplement de décrire ce que vous avez appris sur les relations entre les prédicteurs et la réponse, il n’est pas honteux de citer un intervalle de confiance pour un odds ratio de, par exemple, de 2,3 à la hausse. (En fait, il pourrait sembler risqué de citer des intervalles de confiance fondés sur des estimations non biaisées excluant les rapports de cotes les mieux corroborés par les données.) Des problèmes surviennent lorsque vous essayez de prédire à l'aide d'estimations ponctuelles et le prédicteur sur lequel la séparation se produit submerge les autres.

  6. Utilisez un modèle de régression logistique caché, comme décrit dans Rousseeuw & Christmann (2003), "Robustesse contre la séparation et les valeurs aberrantes dans la régression logistique", Computational Statistics & Data Analysis , 43 , 3 et mis en œuvre dans le package R hlr . (@ user603 suggère cela. ) Je n'ai pas lu le document, mais ils disent dans l'abstrait "un modèle légèrement plus général est proposé, selon lequel la réponse observée est fortement liée mais différente de la réponse vraie non observable", ce qui suggère de: moi, ce ne serait peut-être pas une bonne idée d’utiliser la méthode à moins que cela ne paraisse plausible.

  7. « Modifier quelques observations choisies au hasard de 1 à 0 ou 0 à 1 parmi les variables présentant une séparation complète »: @ robertf de commentaire . Cette suggestion semble découler du fait que la séparation est considérée comme un problème en soi plutôt que comme le symptôme d’un manque d’informations dans les données, ce qui pourrait vous amener à préférer d’autres méthodes à l’estimation avec le maximum de vraisemblance ou à limiter les déductions à celles que vous pouvez utiliser. précision raisonnable - approches qui ont leurs propres mérites et ne sont pas simplement des "solutions" pour la séparation. (Hormis le fait qu’elle est ad hoc sans vergogne , il est désagréable pour la plupart des analystes qui posent la même question de la même donnée, émettant les mêmes hypothèses, de donner des réponses différentes en raison du résultat d’un tirage au sort ou autre.)

Scortchi
la source
1
@Scortchi Il existe une autre option (hérétique). Qu'en est-il de changer quelques observations choisies au hasard de 1 à 0 ou de 0 à 1 parmi les variables présentant une séparation complète?
RobertF
@RobertF: Merci! Je n'avais pas pensé à celui-ci - si vous avez des références concernant ses performances, je vous en serais reconnaissant. Avez-vous rencontré des personnes l'utilisant dans la pratique?
Scortchi
@Scortchi - Non, il est fait référence à des chercheurs ajoutant des données artificielles pour éliminer la séparation complète, mais je n'ai trouvé aucun article sur la modification sélective des données. Je n'ai aucune idée de l'efficacité de cette méthode.
RobertF
1
@tatami: Tous les programmes (nombreux?) ne mettent pas en garde sur la séparation en soi, ce qui peut être délicat à détecter lorsqu'il s'agit d'une combinaison linéaire de plusieurs variables, mais à propos de l'échec de la convergence et / ou de valeurs ajustées proches de zéro ou d'un - je toujours vérifier ces.
Scortchi
2
@Scortchi: très beau résumé dans votre réponse. Personnellement, je privilégie l'approche bayésienne, mais il convient de mentionner la belle analyse du phénomène général d'un point de vue fréquentiste dans projecteuclid.org/euclid.ejs/1239716414 . L'auteur propose des intervalles de confiance unilatéraux qui peuvent être utilisés même en présence d'une séparation complète dans la régression logistique.
Cyan
55

Ceci est une extension des réponses de Scortchi et de Manoel, mais comme vous semblez utiliser RI, je pensais pouvoir vous fournir du code. :)

Je pense que la solution la plus simple et la plus simple à votre problème consiste à utiliser une analyse bayésienne avec des hypothèses antérieures non informatives, comme proposé par Gelman et al. (2008). Comme Scortchi le mentionne, Gelman recommande de placer un coefficient de Cauchy antérieur avec une médiane de 0,0 et une échelle de 2,5 sur chaque coefficient (normalisé avec une moyenne de 0,0 et un écart-type de 0,5). Cela régularisera les coefficients et les tirera légèrement vers zéro. Dans ce cas, c'est exactement ce que vous voulez. En raison de ses queues très larges, le Cauchy autorise toujours des coefficients élevés (par opposition à la Normal à queue courte), de Gelman:

entrez la description de l'image ici

Comment faire cette analyse? Utilisez la bayesglmfonction dans le paquet arm qui implémente cette analyse!

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

Ne fonctionne pas très bien ... Maintenant la version bayésienne:

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 simple, non?

Références

Gelman et al (2008), "Une distribution a priori faiblement informative pour les modèles de régression logistique et autres", Ann. Appl. Stat., 2, 4 http://projecteuclid.org/euclid.aoas/1231424214

Rasmus Bååth
la source
6
Non, trop simple. Pouvez-vous expliquer ce que vous venez de faire? Quel est le prieur qui bayesglmutilise? Si l'estimation de ML est équivalente à Bayesian avec un a priori plat, comment des priors non informatifs peuvent-ils aider ici?
StasK
5
Ajout de quelques informations supplémentaires! Le prieur est vague mais pas plat. Cela a une certaine influence en régularisant les estimations et en les tirant légèrement vers 0,0, ce que je crois que vous voulez dans ce cas.
Rasmus Bååth Le
> m = bayesglm (match ~., famille = binôme (link = 'logit'), data = df) Message d'avertissement: probabilités ajustées numériquement 0 ou 1 survenues Pas terrible!
Chris
En entrée, essayez une régularisation un peu plus forte en augmentant les prior.dfvaleurs par défaut 1.0et / ou en diminuant les prior.scalevaleurs par défaut 2.5, commencez peut-être à essayer:m=bayesglm(match ~. , family = binomial(link = 'logit'), data = df, prior.df=5)
Rasmus Bååth
1
Que faisons-nous exactement lorsque nous augmentons prior.df dans le modèle. Y a-t-il une limite à la hauteur à laquelle nous voulons aller? Si j'ai bien compris, cela contraint le modèle à permettre la convergence avec des estimations d'erreur précises?
Hamilthj
7

L'article de Paul Allison est l'une des explications les plus complètes des problèmes de "séparation quasi complète". Il écrit sur les logiciels SAS, mais les problèmes qu’il aborde peuvent être générés par tous les logiciels:

  • La séparation complète se produit chaque fois qu'une fonction linéaire de x peut générer des prédictions parfaites de y

  • La séparation quasi-complète se produit lorsque (a) il existe un vecteur de coefficient b tel que bxi ≥ 0 lorsque yi = 1 et bxi ≤ 0 * lorsque ** yi = 0 et que cette égalité est valable pour au moins un cas dans chaque catégorie de la variable dépendante. En d’autres termes, dans le cas le plus simple, pour toute variable indépendante dichotomique d’une régression logistique, s’il existe un zéro dans la table 2 × 2 constituée de cette variable et de la variable dépendante, l’estimation ML du coefficient de régression n’existe pas.

Allison discute de nombreuses solutions déjà mentionnées, notamment la suppression de variables à problème, la réduction des catégories, l'inaction, la régression logistique exacte , l'estimation bayésienne et l'estimation du maximum de vraisemblance pénalisé.

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

Mike Hunter
la source
3

warning

Avec des données générées le long des lignes de

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

L'avertissement est fait:

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

ce qui reflète bien évidemment la dépendance inhérente à ces données.

En R, le test de Wald se trouve avec summary.glmou avec waldtestdans l’ lmtestemballage. Le test du rapport de vraisemblance est effectué avec anovaou avec lrtestdans l' lmtestemballage. Dans les deux cas, la matrice d'information a une valeur infinie et aucune inférence n'est disponible. Au contraire, R ne produit sortie, mais vous ne pouvez pas faire confiance. L'inférence que R produit typiquement dans ces cas a des valeurs de p très proches de un. En effet, la perte de précision dans le OU est de plusieurs ordres de grandeur inférieure à la perte de précision dans la matrice de variance-covariance.

Quelques solutions décrites ici:

Utilisez un estimateur en une étape,

Beaucoup de théories soutiennent le faible biais, l’efficacité et la généralisabilité des estimateurs à une étape. Il est facile de spécifier un estimateur en une étape dans R et les résultats sont généralement très favorables pour la prédiction et l'inférence. Et ce modèle ne divergent jamais, car l’itérateur (Newton-Raphson) n’a tout simplement pas la chance de le faire!

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

Donne:

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

Vous pouvez donc voir que les prévisions reflètent le sens de la tendance. Et l'inférence suggère fortement les tendances que nous croyons être vraies.

entrez la description de l'image ici

effectuer un test de score,

La statistique Score (ou Rao) diffère du rapport de vraisemblance et des statistiques de Wald. Il ne nécessite pas d'évaluation de la variance sous l'hypothèse alternative. Nous ajustons le modèle sous le null:

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

Dans les deux cas, vous avez une inférence pour un OU de l'infini.

et utilisez des estimations non biaisées médianes pour un intervalle de confiance.

Vous pouvez produire un IC à 95% non biaisé et non médian pour le rapport de cotes infini en utilisant une estimation médiane sans biais. Le paquet epitoolsen R peut le faire. Et je donne un exemple d'application de cet estimateur ici: Intervalle de confiance pour l'échantillonnage de Bernoulli

AdamO
la source
2
C’est formidable, mais j’ai bien sûr quelques doutes: (1) Le test du rapport de vraisemblance n’utilise pas la matrice d’information; seul le test de Wald y parvient et échoue de manière catastrophique en cas de séparation. (2) Je ne connais pas du tout les estimateurs en une étape, mais l'estimation de la pente ici semble absurdement basse. (3) Un intervalle de confiance n'est pas sans biais médian. Ce que vous associez à cette section est l’intervalle de confiance moyen. (4) Vous pouvez obtenir des intervalles de confiance en inversant les tests de LR ou de score. ...
Scortchi
... (5) Vous pouvez effectuer le test de score dans R en donnant l'argument test="Rao"à la anovafonction. (
Et
@scortchi bon à savoir anova a des tests de score par défaut! Peut-être qu'une implémentation manuelle est utile. Les IC ne sont pas des valeurs sans biais médianes, mais les IC de l'estimateur sans valeur médiane fournissent une inférence cohérente pour les paramètres limites. Le milieu p est un tel estimateur. Le p peut être transformé en un rapport de cotes b / c il est invariant aux transformations un à un. Le test LR est-il cohérent pour les paramètres de limite?
AdamO
Seule l'hypothèse nulle ne doit pas contenir de paramètres à une limite pour que le théorème de Wilks puisse s'appliquer, bien que les tests score et LR soient approximatifs dans des échantillons finis.
Scortchi
2

Soyez prudent avec ce message d'avertissement de R. Jetez un coup d'œil à ce billet d'Andrew Gelman sur le blog et vous verrez que ce n'est pas toujours un problème de séparation parfaite, mais parfois un bug glm. Il semble que si les valeurs de départ sont trop éloignées de l'estimation du maximum de vraisemblance, elle explose. Alors, vérifiez d’abord avec d’autres logiciels, comme Stata.

Si vous avez vraiment ce problème, vous pouvez essayer d’utiliser la modélisation bayésienne, avec des priors informatifs.

Mais en pratique, je me débarrasse simplement des prédicteurs à l’origine du problème, car je ne sais pas comment choisir un préalable informatif. Mais je suppose qu’il existe un article de Gelman sur l’utilisation informative avant lorsque vous avez ce problème de séparation parfaite. Il suffit de google. Peut-être devriez-vous essayer.

Manoel Galdino
la source
8
Le problème avec la suppression des prédicteurs est que vous supprimez le prédicteur qui explique le mieux la réponse, ce qui correspond généralement à ce que vous souhaitez faire! Je dirais que cela n’a de sens que si vous surévalez votre modèle, par exemple en adaptant trop d’interactions complexes.
Simon Byrne
4
Pas un bug, mais un problème avec les estimations initiales trop éloignées du MLE, qui ne se poseront pas si vous n'essayez pas de les choisir vous-même.
Scortchi
Je comprends cela, mais je pense que c'est un bogue dans l'algorithme.
Manoel Galdino
5
Eh bien, je ne veux pas contester la définition de «bug». Mais le comportement n'est ni insondable ni corrigible en base R - vous n'avez pas besoin de "vérifier avec un autre logiciel". Si vous souhaitez traiter automatiquement de nombreux problèmes de non-convergence, le glm2package implémente un contrôle indiquant que la probabilité augmente réellement à chaque étape du score, et divise par deux la taille de l'étape, si ce n'est pas le cas.
Scortchi
3
Il existe (sur CRAN) le package R safeBinaryRegression conçu pour diagnostiquer et résoudre de tels problèmes, en utilisant des méthodes d'optimisation pour vérifier de manière certaine s'il existe une séparation ou une quasi-séparation. Essayez le!
kjetil b halvorsen
2

Je ne suis pas sûr d’être d’accord avec les affirmations de votre question.

Je pense que ce message d’alerte signifie, pour certains des X observés niveau dans vos données, la probabilité ajustée est numériquement égale à 0 ou 1. En d’autres termes, à la résolution, elle indique 0 ou 1.

Tu peux courir predict(yourmodel,yourdata,type='response') et vous y trouverez 0 ou / et 1 comme probabilités prédites.

En conséquence, je pense qu'il est correct d'utiliser simplement les résultats.

StayLearning
la source
-1

Je comprends que c’est un vieux billet, mais je vais quand même répondre, car j’ai eu du mal à le faire pendant des jours et que cela peut aider les autres.

La séparation complète se produit lorsque les variables sélectionnées pour s’adapter au modèle peuvent différencier très précisément les valeurs 0 et 1 ou les valeurs oui et non. Toute notre approche de la science des données est basée sur une estimation de probabilité mais elle échoue dans ce cas.

Étapes de rectification: -

  1. Utilisez bayesglm () au lieu de glm (), lorsque la variance entre les variables est faible

  2. Parfois, utiliser (maxit = "une valeur numérique") avec bayesglm () peut aider

3. Troisième et plus important contrôle pour les variables sélectionnées pour l'ajustement du modèle, il doit exister une variable pour laquelle la multi-colinéarité avec la variable Y (sortie) est très élevée, ignorez cette variable de votre modèle.

Comme dans mon cas, j'avais des données de désabonnement des télécommunications pour prédire le désabonnement des données de validation. J'avais une variable dans mes données d'entraînement qui pouvait très bien faire la différence entre le oui et le non. Après l'avoir laissé tomber, je pourrais obtenir le bon modèle. De plus, vous pouvez utiliser pas à pas (ajustement) pour rendre votre modèle plus précis.

yash
la source
2
Je ne vois pas que cette réponse ajoute beaucoup à la discussion. L’approche bayésienne est couverte en détail dans les réponses précédentes, la suppression des prédicteurs "problématiques" est également déjà mentionnée (et découragée). La sélection variable par étapes est rarement une bonne idée à ma connaissance.
einar