# How to simulate artificial data for logistic regression?

I know I'm missing something in my understanding of logistic regression, and would really appreciate any help.

As far as I understand it, the logistic regression assumes that the probability of a '1' outcome given the inputs, is a linear combination of the inputs, passed through an inverse-logistic function. This is exemplified in the following R code:

#create data:
x1 = rnorm(1000)           # some continuous variables
x2 = rnorm(1000)
z = 1 + 2*x1 + 3*x2        # linear combination with a bias
pr = 1/(1+exp(-z))         # pass through an inv-logit function
y = pr > 0.5               # take as '1' if probability > 0.5

#now feed it to glm:
df = data.frame(y=y,x1=x1,x2=x2)
glm =glm( y~x1+x2,data=df,family="binomial")


and I get the following error message:

Warning messages: 1: glm.fit: algorithm did not converge 2: glm.fit: fitted probabilities numerically 0 or 1 occurred

I've worked with R for some time now; enough to know that probably I'm the one to blame.. what is happening here?

-
The way you simulate your data looks weird to me. If you want, for an alternative more standard way, you can have a look here: stats.stackexchange.com/questions/12857/… –  ocram Dec 25 '12 at 15:31
@ocram: you are right; this is a duplicate question! –  user603 Dec 25 '12 at 15:34
I did run an erroneous simulation, as @Stéphane Laurent explained. However, the problem was perfect separation in logistic regression, a problem I was not familiar with, and that I was rather surprised to learn about. –  zorbar Dec 26 '12 at 7:42
@zorbar: it was in my response to your question (now deleted). –  user603 Dec 26 '12 at 9:23
@user603: I probably missed your response; Thanks anyway –  zorbar Dec 26 '12 at 9:53

No. The response variable $y_i$ is a Bernoulli random variable taking value $1$ with probability $pr(i)$.

> set.seed(666)
> x1 = rnorm(1000)           # some continuous variables
> x2 = rnorm(1000)
> z = 1 + 2*x1 + 3*x2        # linear combination with a bias
> pr = 1/(1+exp(-z))         # pass through an inv-logit function
> y = rbinom(1000,1,pr)      # bernoulli response variable
>
> #now feed it to glm:
> df = data.frame(y=y,x1=x1,x2=x2)
> glm( y~x1+x2,data=df,family="binomial")

Call:  glm(formula = y ~ x1 + x2, family = "binomial", data = df)

Coefficients:
(Intercept)           x1           x2
0.9915       2.2731       3.1853

Degrees of Freedom: 999 Total (i.e. Null);  997 Residual
Null Deviance:      1355
Residual Deviance: 582.9        AIC: 588.9

-
You're right - I've missed this step. thanks a lot for your help! –  zorbar Dec 25 '12 at 16:09