0
$\begingroup$

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?

Improve this question
$\endgroup$
2
  • $\begingroup$ The posts bearing the tag [hauck-donner effect] will be very helpful in answering this question. $\endgroup$
    – Sycorax ♦
    Jan 19, 2016 at 18:50
  • $\begingroup$ 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. $\endgroup$
    – dsaxton
    Jan 19, 2016 at 19:02

2 Answers 2

Reset to default
3
$\begingroup$

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))
Improve this answer
$\endgroup$
2
$\begingroup$

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.

Improve this answer
$\endgroup$

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