Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

If you have a variable which perfectly separates zeroes and ones in target variable, R will yield the following "perfect or quasi perfect separation" warning message:

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

We still get the model but the coefficient estimates are inflated.

How do you deal with this in practice?

share|improve this question
3  
related question –  user603 Dec 13 '12 at 14:56

4 Answers 4

up vote 25 down vote accepted

A solution to this is to utilize a form of penalized regression. In fact, this is the original reason some of the penalized regression forms were developed (although they turned out to have other interesting properties.

Install and load package glmnet in R and you're mostly ready to go. One of the less user-friendly aspects of glmnet is that you can only feed it matrices, not formulas as we're used to. However, you can look at model.matrix and the like to construct this matrix from a data.frame and a formula...

Now, when you expect that this perfect separation is not just a byproduct of your sample, but could be true in the population, you specifically don't want to handle this: use this separating variable simply as the sole predictor for your outcome, not employing a model of any kind.

share|improve this answer
1  
You can also use a formula interface for glmnet through the caret package. –  Zach Nov 4 '13 at 17:48

You've several options:

  1. Remove some of the bias.

    (a) By penalizing the likelihood as per @Nick's suggestion. Package logistf in R or the FIRTH option in SAS's PROC LOGISTIC implement the method proposed in Firth (1993), "Bias reduction of maximum likelihood estimates", Biometrika, 80,1.; which removes the first-order bias from maximum likelihood estimates. (Here @Gavin recommends the brglm package, which I'm not familiar with, but I gather it implements a similar approach for non-canonical link functions e.g. probit.)

    (b) By using median-unbiased estimates in exact conditional logistic regression. Package logistiX in R or the EXACT statement in SAS's PROC LOGISTIC.

  2. Exclude cases where the predictor category or value causing separation occurs. These may well be outside your scope; or worthy of further, focused investigation. (The R package safeBinaryRegression is handy for finding them.)

  3. Re-cast the model. Typically this is something you'd have done beforehand if you'd thought about it, because it's too complex for your sample size.

    (a) Remove the predictor from the model. Dicey, for the reasons given by @Simon: "You're removing the predictor that best explains the response".

    (b) By collapsing predictor categories / binning the predictor values. Only if this makes sense.

    (c) Re-expressing the predictor as two (or more) crossed factors without interaction. Only if this makes sense.

  4. Use a Bayesian analysis as per @Manoel's suggestion. Though it seems unlikely you'd want to just because of separation, worth considering on its other merits.The paper he recommends is Gelman et al (2008), "A weakly informative default prior distribution for logistic & other regression models", Ann. Appl. Stat., 2, 4: the default in question is an independent Cauchy prior for each coefficient, with a mean of zero & a scale of $\frac{5}{2}$; to be used after standardizing all continuous predictors to have a mean of zero & a standard deviation of $\frac{1}{2}$. If you can elucidate strongly informative priors, so much the better.

  5. Do nothing. (But calculate confidence intervals based on profile likelihoods, as the Wald estimates of standard error will be badly wrong.) An often over-looked option. If the purpose of the model is just to describe what you've learnt about the relationships between predictors & response, there's no shame in quoting a confidence interval for an odds ratio of, say, 2.3 upwards. (Indeed it could seem fishy to quote confidence intervals based on unbiased estimates that exclude the odds ratios best supported by the data.) Problems come when you're trying to predict using point estimates, & the predictor on which separation occurs swamps the others.

  6. Use a hidden logistic regression model, as described in Rousseeuw & Christmann (2003),"Robustness against separation and outliers in logistic regression", Computational Statistics & Data Analysis, 43, 3, and implemented in the R package hlr. (@user603 suggests this.) I haven't read the paper, but they say in the abstract "a slightly more general model is proposed under which the observed response is strongly related but not equal to the unobservable true response", which suggests to me it mightn't be a good idea to use the method unless that sounds plausible.

share|improve this answer
    
Outstanding answer, @Scortchi! –  abaumann Feb 28 at 11:00

This is an expansion of Scortchi and Manoel's answers, but since you seem to use R I thought I'd supply some code. :)

I believe the easiest and most straightforward solution to your problem is to use a Bayesian analysis with non-informative prior assumptions as proposed by Gelman et al (2008). As Scortchi mentions, Gelman recommends to put a Cauchy prior with mean 0.0 and scale 2.5 on each coefficient (normalized to have mean 0.0 and a SD of 0.5). This will regularize the coefficients and pull them just slightly towards zero. In this case it is exactly what you want. Due to having very wide tails the Cauchy still allows for large coefficients (as opposed to the short tailed Normal), from Gelman:

enter image description here

How to run this analysis? Use the bayesglm function in arm package that implements this analysis!

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

Does not work that well... Now the Bayesian version:

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, no?

References

Gelman et al (2008), "A weakly informative default prior distribution for logistic & other regression models", Ann. Appl. Stat., 2, 4 http://projecteuclid.org/euclid.aoas/1231424214

share|improve this answer
1  
No. Too simple. Can you explain what you have just done? What is the prior that bayesglm uses? If ML estimation is equivalent to Bayesian with a flat prior, how do non-informative priors help here? –  StasK Feb 28 at 2:14
1  
Added some more info! The prior is vague but not flat. It has some influence as it regularizes the estimates and pull them slightly towards 0.0 which is what I believe you want in this case. –  Rasmus Bååth Feb 28 at 9:57

Be careful with this warning message from R. Take a look at this blog post by Andrew Gelman, and you will see that it is not always a problem of perfect separation, but sometimes a bug with glm. It seems that if the starting values are too far from the maximum-likelihood estimate, it blows up. So, check first with other software, like Stata.

If you really have this problem, you may try to use Bayesian modeling, with informative priors.

But in practice I just get rid of the predictors causing the trouble, because I don't know how to pick an informative prior. But I guess there is a paper by Gelman about using informative prior when you have this problem of perfect separation problem. Just google it. Maybe you should give it a try.

share|improve this answer
4  
The problem with removing predictors is that you're removing the predictor that best explains the response, which is usually what you aiming to do! I would argue that this only makes sense if you've overfit your model, for example by fitting too many complicated interactions. –  Simon Byrne Jun 29 '11 at 12:36
3  
Not a bug, but a problem with the initial estimates being too far from the MLE, which won't arise if you don't try to choose them yourself. –  Scortchi Sep 1 '13 at 15:44
    
I understand this, but I do think this is a Bug in the algorithm. –  Manoel Galdino Sep 19 '13 at 14:38
2  
Well I don't want to quibble about the definition of 'bug'. But the behaviour's neither unfathomable nor unfixable in base R - you don't need to "check with other software". If you want to deal automatically with many non-convergence problems, the glm2 package implements a check that the likelihood's actually increasing at each scoring step, & halves the step size if it isn't. –  Scortchi Oct 7 '13 at 11:36

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.