# 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/… Commented Dec 25, 2012 at 15:31
• @ocram: you are right; this is a duplicate question! Commented Dec 25, 2012 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. Commented Dec 26, 2012 at 7:42
• @zorbar: it was in my response to your question (now deleted). Commented Dec 26, 2012 at 9:23
• @user603: I probably missed your response; Thanks anyway Commented Dec 26, 2012 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! Commented Dec 25, 2012 at 16:09
• I had a question regarding the way you simulate data. When we simulate data for linear regression, we also simulate a noise (\epsilon). I understand that logistic function is a function of the expectation which by itself cancels the noise out. Is that the reason that you don't have any noise in your z ?
– Sam
Commented Oct 28, 2014 at 17:08
• @Sepehr Linear regression assumes a Gaussian distribution. The "noise" is just an interpretation of the variability around the mean, but this is not a noise added to the response : the response is written as $\text{mean response} + \text{noise}$, this is just an interpretation. The logistic regression assumes the response has a binomial distribution, and similarly there's no noise added to the response. Commented Oct 29, 2014 at 8:10
• @StéphaneLaurent, exactly. I completely get. Thanks very much for your answer.
– Sam
Commented Nov 1, 2014 at 17:53
• @StéphaneLaurent How should I generate the data in the case where the coefficients have l2 regularization on them? Commented Jan 19, 2021 at 17:40

LogisticRegression is suitable for fitting if probabilities or proportions are provided as the targets, not only 0/1 outcomes.

import numpy as np
import pandas as pd
def logistic(x, b, noise=None):
L = x.T.dot(b)
if noise is not None:
L = L+noise
return 1/(1+np.exp(-L))

x = np.arange(-10., 10, 0.05)
bias = np.ones(len(x))
X = np.vstack([x,bias]) # Add intercept
B =  [1., 1.] # Sigmoid params for X

# True mean
p = logistic(X, B)
# Noisy mean
pnoisy = logistic(X, B, noise=np.random.normal(loc=0., scale=1., size=len(x)))
# dichotomize pnoisy -- sample 0/1 with probability pnoisy
dichot = np.random.binomial(1., pnoisy)

pd.Series(p, index=x).plot(style='-')
pd.Series(pnoisy, index=x).plot(style='.')
pd.Series(dichot, index=x).plot(style='.')


Here we have three potential targets for logistic regression. p which is the true/target proportion/probability, pnoisy which is p with normal noise added in the log odds scale, and dichot, which is pnoisy treated as a parameter to the binomial PDF, and sampled from that. You should test all 3 -- I found some open source LR implementations can't fit p.

Depending on your application, you may prefer pnoisy.

In practice, you should also consider how the noise is likely to be shaped in you target application and try to emulate that.