# Logistic regression performance with high number of predictors

I'm trying to understand the behavior of logistic regression in high dimensional problems (i.e. when you are fitting a logistic regression to data with a high number of predictor variables).

Every time I fit a logistic regression with a high number of predictors, I get the following warning in R:

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


I read somewhere that this is called the Hauck-Donner phenomenon and it is due to the fact that the fitted probabilities are too close to 0/1.

However, my hypothesis is that as you increase the number of predictors, the variances of your fitted probabilities has to increase in the logistic regression. This is because your log odds estimates is essentially a sum of random variables, and a sum of 100 similarly distributed random variables will (likely) have a larger variance than a sum of 10 rvs. Therefore, when you apply logistic regression to very high dimensional problems, your fitted probabilities will be closer to 0/1 (because higher variance), and therefore your coefficient estimates must be biased (incorrect) because of this problem? Is this hypothesis correct?

I've created a simulation with the following code to try to answer this:

genLogit <- function(n,dimens){
dimens <- floor(dimens/2)*2 #make sure dimens is even

xdata <- cbind(replicate(dimens/2,runif(n)),replicate(dimens/2,runif(n,-1,0)))
ydata <- apply(xdata,1,sum)

prob <- exp(ydata)/(1+exp(ydata))
runis <- runif(length(prob))

ydata <- ifelse(runis<prob,1,0)

model <- glm(ydata~.,data = data.frame(cbind(ydata,xdata)),family = binomial(link =
return(summary(model))
}


What the code is doing is basically simulating a logistic regression from the following:

$$\log\left(\frac{p}{1-p}\right) = U_1+U_1+\ldots+U_1 + U_{-1}+U_{-1}+\ldots+U_{-1}$$

where $U_1 = \text{Unif}(0,1)$ and $U_{-1} = \text{Unif}(-1,0)$. You can vary the number of predictors in the model as well as the number of data points generated. Then the function will fit a logistic regression to the simulated data, and you can examine the coefficients, residual deviance, fit, etc.

I understand that I all of my predictors have the same variance (which is not necessarily true when dealing with real data) but is this simulation still sufficient to prove my hypothesis?

• Your code is not complete. – user88 May 1 '12 at 9:17
• Your code looks like its prone to produce quasi-complete separation: Because of the first half being in [0,1] and the second half being in [-1,0] and the way you generate the zeros and ones by the deterministic rule ydata <- ifelse(runis<prob,1,0). – Momo May 3 '12 at 18:45

I think we should give the word to Venables and Ripley, page 198 in MASS:

There is one fairly common circumstance in which both convergence problems and the Hauck-Donner phenomenon can occur. This is when the fitted probabilities are extremely close to zero or one. Consider a medical diagnosis problem with thousands of cases and around fifty binary explanatory variables (which may arise from coding fewer categorical factors); one of these indicators is rarely true but always indicates that the disease is present. Then the fitted probabilities of cases with that indicator should be one, which can only be achieved by taking $\hat\beta_i = \infty$. The result from glm will be warnings and an estimated coefficient of around +/- 10.

Besides potential numerical difficulties there is no formal problem with probabilities being estimated numerically to 0 or 1. However, the $t$-test, which is based on a quadratic approximation, for testing the hypothesis $\beta_i = 0$ can become a poor approximation of the likelihood ratio test, and the $t$-test may appear insignificant though in reality the hypothesis is definitely wrong. As I understand it, this it what the warning is about.

With many predictors a situation like the one Venables and Ripley describes may easily occur; one predictor is mostly not informative, but in certain cases it is a strong predictor for a case.

• I think the MASS example is a little bit different because there's just one significant predictor. I'm wondering whether logistic regression is biased if there are many significant predictors? When I run my code with a high number of predictors, the fitted coefficients start to diverge from 1. – Michael May 1 '12 at 12:59
• Is it true that if you have a set of predictors X and then add another predictor Y to the set (so now the set of predictors is X+Y) then that will increase the variance of the log odds? Wouldn't that increase the variance of the fitted probabilities and be a cause of the Hauck-Donner phenomenon? – Michael May 1 '12 at 13:58
• @Michael, what you observe, I believe, is the increased variance of the fitted coefficients when you increase the number of predictors. It is not a bias phenomenon. Hauck-Donner refers specifically to the problem with the inaccurate t-test statistic. – NRH May 1 '12 at 16:44
• This is a fairly common problem but even a tiny bit of regularization can fix it. Standardized coefficients in the +/- 5 range will shift probability from .01 to .99. For sensible priors, see the bayesglm function in the ARM package and Gelman's paper on the subject: stat.columbia.edu/~gelman/research/published/priors11.pdf – Tristan May 2 '12 at 4:17

While the Hauck-Donner effect is closely related to it, I think the problem in your case is (quasi)-complete separation. This refers to the phenomenon that a certain combination (including interactions) of predictors will give rise to a subset of the observations where you observe only zeros or ones (basically a combination of predictor values will separate the two classes). Then the maximum likelihood estimate will not exist (it will be infinite, which makes its standard error rather large too). This is what V&R write in the quote by @NRH. If you have many predictors, especially categorical ones, this just becomes more likely to happen for a particular combination of predictors. The HD effect then occurs for the Wald test in such a situation. A canonical treatment of qcs is Albert and Anderson 1984 article in Biometrika.

You might want to look at the noverlap package (no longer on CRAN) in R which contains utilities to deal with quasi-complete separation, or the brglm package, that can deal with qcs.

I'm not quite sure how to explain your problem, but I can offer a potential solution -- try using an R package called glmnet instead. I have used glmnet for both linear and logistic regression. In one particular problem, I had approximately 1,200 cases (i.e. N = 1,200) and about 110 predictors. The package gave me great results.

Of course, it's worth pointing out that glmnet is primarily used for penalized logistic regression, but since the package lets you select the degree of penalty to apply, I'm sure you can set the penalty to zero to obtain the results of regular logistic regression (i.e. one with no penalty). In any case, the package was made specifically for problems in high dimensions (even N << P), and this seems to be the underlying source of your problem. I highly recommend glmnet.

• I think glmnet is a great tool and would give reasonable results in this circumstance. However, you can't set the penalty to Zero. Glmnet uses a coordinate descent algorithm, but it still will have diverging coefficient values just the same (and even has some exception code to handle this if you look at the source code). You can however just include a very small penalty. But why not cross validate and pick a good lambda? – Shea Parkes May 2 '12 at 3:49
• That is not the way I read the documentation. Setting the penalty to zero should give you ordinary logistic regression results. True, that won't solve the OP's problem, but I do not believe it is correct that "you cannot set the penalty to zero." – DWin Aug 20 '13 at 18:22