I'm trying to simulate a logistic regression. My goal is showing that if Y=1 is rare, than the intercept is biased. In my R script I define the logistic regression model through the latent variable's approach (see for example pp. 140 http://gking.harvard.edu/files/abs/0s-abs.shtml):

x   <- rnorm(10000)

b0h <- numeric(1000)
b1h <- numeric(1000)

for(i in 1:1000){
  eps <- rlogis(10000)
  eta <- 1+2*x+eps
  y   <-numeric(10000)
  y   <- ifelse (eta>0,1,0)

  m      <- glm(y~x,family=binomial)
  b0h[i] <- coef(m)[1]
  b1h[i] <- coef(m)[2]


The problem here is that I don't know how to force the observations y to be balanced before (50:50), then unbalanced (90:10). As we can see with the function table(), in my script the proportion of ones is random.


How to solve this problem?

  • $\begingroup$ If you are truly simulating from a logistic regression model then the proportion has to be random. $\endgroup$
    – Xi'an
    Dec 9 '14 at 10:57
  • $\begingroup$ Thanks for your answer @Xi'an. I know that the proportion has to be random, but I just want to show that if the class is skewed (one class is rare), than if I use logistic regression the intercept is biased (this is what I think to have understood from gking.harvard.edu/files/abs/0s-abs.shtml. How to force the class to be skewed in my Rscript? $\endgroup$
    – Luca Dibo
    Dec 9 '14 at 11:15
  • $\begingroup$ Then you have to sample from another model. $\endgroup$
    – Xi'an
    Dec 9 '14 at 11:46
  • $\begingroup$ Related: How to simulate artificial data for logistic regression? $\endgroup$ Dec 9 '14 at 19:23
  • 1
    $\begingroup$ You're never going to show this with a simulation! All you could possibly discover is that a certain amount of bias persists in the largest simulations your computer can handle. I am having a hard time understanding why the intercept should be biased asymptotically--you actually contradict yourself with this statement (since "asymptotically unbiased" implies the bias becomes unmeasurable with sufficiently large samples)--so I would guess it might have something to do with the nature of the asymptotics. $\endgroup$
    – whuber
    Dec 9 '14 at 23:02

Logistic regression doesn't really have an error term. Alternatively, you can think of the response distribution (the binomial) as having its random component intrinsically 'built-in' (for more, it may help to read my answer here: Difference between logit and probit models). As a result, I think it is conceptually clearer to generate data for simulations directly from a binomial parameterized as the logistic transformation of the structural component of the model, rather than use the logistic as a sort of error term.

From there, if you want the long run probability that $Y = 1$ to be $.5$ or $.1$, you just need your structural component to be balanced around $0$ (for $.5$), or $-2.197225$ (for $.1$). I got those values by converting the response probability to the log odds:
$$ \log(\text{odds}(Y=1)) = \frac{\exp(Pr(Y = 1))}{(1+\exp(Pr(Y = 1))} $$ The most convenient way to do this will be to use those values for your intercept ($\beta_0$) and have your slope be $0$. (Alternatively, you can use any two parameter values, $\beta_0$ and $\beta_1$, that you like such that, given your $X$ values, the log mean odds equals, e.g., $-2.197225$.) Here is an example with R code:

lo2p = function(lo){      # this function will perform the logistic transformation
  odds = exp(lo)          #   of the structural component of the data generating
  p    = odds / (1+odds)  #   process

N     = 1000              # these are the true values of the DGP
beta0 = -2.197225
beta1 = 0

set.seed(8361)            # this makes the simulation exactly reproducible
x     = rnorm(N)
lo    = beta0 + beta1*x
p     = lo2p(lo)          # these will be the parameters of the binomial response

b0h   = vector(length=N)  # these will store the results
b1h   = vector(length=N)
y1prp = vector(length=N)  # (for the proportion of y=1)

for(i in 1:1000){         # here is the simulation
  y        = rbinom(n=N, size=1, prob=p)
  m        = glm(y~x, family=binomial)
  b0h[i]   = coef(m)[1]
  b1h[i]   = coef(m)[2]
  y1prp[i] = mean(y)

mean(b0h)                 # these are the results
# [1] -2.205844
# [1] -0.0003422177
# [1] 0.100036

enter image description here

enter image description here

  • $\begingroup$ You don't want to control the mean log odds to achieve a given proportion of successes; you want to control the (weighted) mean probability. The two means are not quite equivalent. $\endgroup$
    – whuber
    Dec 9 '14 at 19:32
  • $\begingroup$ @whuber, hmmm I think I see your point: because of the intervening non-linear transformation, the mean of p(Y=1) for a large sample of simulated datasets could diverge from lo(Y=1) because of the asymmetry in the transformations of the datasets in the tails. Is that what you're getting at? (It does seem to work awfully well w/ lo(Y=1)=.1, but maybe that's not extreme enough.) $\endgroup$ Dec 9 '14 at 19:44
  • $\begingroup$ Yes, that's what I mean. Because the logistic transformation is nearly linear for values close to zero (around $-1$ to $1$), the difference really shows up only for more extreme values--precisely the case you are discussing in your answer! For instance, let $p_1=0.01$ and $p_2=0.19$. Although their mean is $0.1$, their mean log odds corresponds to a probability of only $0.046$. Rare events have more "leverage" in the averaging when it is done in terms of log odds. $\endgroup$
    – whuber
    Dec 9 '14 at 19:55
  • $\begingroup$ @whuber, do you think I should delete this answer? $\endgroup$ Dec 9 '14 at 19:59
  • 1
    $\begingroup$ @whuber, I changed that phrasing, let me know if you think it needs more. BTW, I am aware of replicate, but I don't tend to use such here; my goal is to make my R code as simple & transparent as possible, even if clunky, in the hopes that even those who don't know R could get the gist. $\endgroup$ Dec 9 '14 at 21:24

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