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

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?

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.

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. $\frac{5}{2}$$\frac{1}{2}$

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

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:

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

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

warning$\infty$

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.

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)

$\chi^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

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.

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.

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.