logistic regression simulation doesn't converge

Question

I want to do a logistic regression simulation using R

I use this code

set.seed(666)
age = rnorm(60)         
blood_pressure = rnorm(60)
race = sample(c(rep(1,30),rep(0,30)))
inactivity = sample(c(rep(1,30),rep(0,30)))
weight = rnorm(60)

z=1+1*age+blood_pressure*2+3*weight+3*inactivity+0*race
pr=exp(z)/(1+exp(z))
y=rbinom(60,1,pr)

df = data.frame(y=y,age,blood_pressure,inactivity,weight,race)
glm(y~age+blood_pressure+inactivity+weight+race,data=df,family=binomial(link='logit'),control = list(maxit = 50))

I got very strange result from it.

Coefficients:
   (Intercept)             age  blood_pressure      inactivity          weight            race  
        -39.75           46.64          106.65          143.52          229.75          100.87

And it says the model doesn't converge.

Does someone know what's wrong and how to fix it?

The posts bearing the tag [hauck-donner effect] will be very helpful in answering this question. — Sycorax, Jan 19 '16 at 18:50
Another thing you might try is fitting a model without an intercept. It's possible that `race` and `inactivity` together may be behaving too much like an intercept, but this is just a guess. — dsaxton, Jan 19 '16 at 19:02

score 3 · Answer 1 · edited Jan 19 '16 at 18:50

What else would you expect from setting the seed to 666?

But more seriously, the issue is that you have a very small sample size (60). As such, the chances that you can find a linear combination of your covariates such that observations with the linear combination below some threshold are all 0's, while observations with the linear combination above some threshold are all 1's is reasonably high. If you can find this, then at least one of your estimated coefficients will be infinity, causing all kinds of problems.

To illustrate, try increasing your sample size:

n <- 600
set.seed(555)
age = rnorm(n)         
blood_pressure = rnorm(n)
race = sample(c(rep(1,n/2),rep(0,n/2)))
inactivity = sample(c(rep(1,n/2),rep(0,n/2)))
weight = rnorm(n)

z=1+1*age+blood_pressure*2+3*weight+3*inactivity+0*race
pr=exp(z)/(1+exp(z))
y=rbinom(n,1,pr)

df = data.frame(y=y,age,blood_pressure,inactivity,weight,race)
glm(y~age+blood_pressure+inactivity+weight+race,data=df,
    family=binomial(link='logit'),control = list(maxit = 50))

score 2 · Answer 2 · edited Jan 19 '16 at 18:51

Look at the summary of your estimates and play with individual variables. Currently the summary returns:

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)
(Intercept)        -39.75  982476.12       0        1
age                 46.64  258055.84       0        1
blood_pressure     106.65  356355.42       0        1
inactivity         143.52  579221.32       0        1
weight             229.75  500237.07       0        1
race               100.87 1032211.12       0        1

This is clearly wrong. When you remove race you get

               Estimate Std. Error z value Pr(>|z|)
(Intercept)       4.368      4.399   0.993    0.321
age               2.399      2.039   1.177    0.239
blood_pressure    9.120      6.596   1.383    0.167
inactivity        8.654      5.854   1.478    0.139
weight           15.164     10.705   1.417    0.157

Hence including race made the model misbehave. You may have either multicollinearity issues, or the problem that all your observations are 0 on one side of a hyperplane, and 1 on the other side, as Cliff AB suggested.

logistic regression simulation doesn't converge

2 Answers2